# The excess of cool supergiants from contemporary stellar evolution models
defies the metallicity-independent Humphreys-Davidson limit

Avishai Gilkis1, Tomer Shenar2, Varsha Ramachandran3, Adam S. Jermyn4
agilkis<EMAIL_ADDRESS>Laurent Mahy2, Lidia M.
Oskinova3, Iair Arcavi1,5 and Hugues Sana2
1 The School of Physics and Astronomy, Tel Aviv University, Tel Aviv 6997801,
Israel
2 Institute of Astrophysics, KU Leuven, Celestijnenlaan 200 D, 3001 Leuven,
Belgium
3 Institut für Physik und Astronomie, Universität Potsdam, Karl-Liebknecht-
Str. 24 / 25, D-14476 Potsdam, Germany
4 Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, New
York, NY 10010, USA
5 CIFAR Azrieli Global Scholars program, CIFAR, Toronto, Canada

###### Abstract

The Humphreys-Davidson (HD) limit empirically defines a region of high
luminosities ($\log_{10}\left(L/\mathrm{L}_{\odot}\right)\gtrsim 5.5$) and low
effective temperatures ($T_{\rm eff}\lesssim 20\,{\rm kK}$) on the
Hertzsprung-Russell Diagram in which hardly any supergiant stars are observed.
Attempts to explain this limit through instabilities arising in near- or
super-Eddington winds have been largely unsuccessful. Using modern stellar
evolution we aim to re-examine the HD limit, investigating the impact of
enhanced mixing on massive stars. We construct grids of stellar evolution
models appropriate for the Small and Large Magellanic Clouds (SMC, LMC), as
well as for the Galaxy, spanning various initial rotation rates and convective
overshooting parameters. Significantly enhanced mixing apparently steers
stellar evolution tracks away from the region of the HD limit. To quantify the
excess of over-luminous stars in stellar evolution simulations we generate
synthetic populations of massive stars, and make detailed comparisons with
catalogues of cool ($T_{\mathrm{eff}}\leq 12.5\,\mathrm{kK}$) and luminous
($\log_{10}\left(L/\mathrm{L}_{\odot}\right)\geq 4.7$) stars in the SMC and
LMC. We find that adjustments to the mixing parameters can lead to agreement
between the observed and simulated red supergiant populations, but for hotter
supergiants the simulations always over-predict the number of very luminous
($\log_{10}\left(L/\mathrm{L}_{\odot}\right)\geq 5.4$) stars compared to
observations. The excess of luminous supergiants decreases for enhanced
mixing, possibly hinting at an important role mixing has in explaining the HD
limit. Still, the HD limit remains unexplained for hotter supergiants.

###### keywords:

stars: evolution – stars: massive

††pubyear: 2021††pagerange: The excess of cool supergiants from contemporary
stellar evolution models defies the metallicity-independent Humphreys-Davidson
limit–References

##  1 Introduction

The upper-right part of the Hertzsprung-Russell diagram (HRD) features a stark
absence of observed stars (Fig. 1), a phenomenon termed the Humphreys-Davidson
(HD) limit (Humphreys & Davidson, 1979). While the luminous blue variables
(LBVs) venture into this region during outbursts, with a few exceptions, cool
supergiants (CSGs), which comprise red, yellow, and blue supergiants (RSGs,
YSGs, BSGs) with effective temperatures $T_{\rm eff}\lesssim
12.5\,\mathrm{kK}$ and luminosities $\log_{10}\left(L_{\rm
max}/\mathrm{L}_{\odot}\right)\gtrsim 5.5$ (i.e., initial masses $M_{\rm
i}\gtrsim 30\,\mathrm{M}_{\odot}$) are not observed. In contrast, hundreds of
main sequence progenitors with $M_{\rm i}\gtrsim 30\,\mathrm{M}_{\odot}$,
ranging all the way to $\approx 150\,\mathrm{M}_{\odot}$ and perhaps more,
were directly observed (Crowther et al., 2010; Bestenlehner et al., 2011;
Almeida et al., 2017; Shenar et al., 2017; Tehrani et al., 2019; Mahy et al.,
2020). This implies two possibilities: either stars with $M_{\rm i}\gtrsim
30\,\mathrm{M}_{\odot}$ skip the CSG phase altogether, or they experience this
phase very briefly, making it observationally rare.

Despite various attempts, no evolutionary models are currently capable of
reproducing the observed absence of stars beyond the HD limit (Davies,
Crowther & Beasor, 2018; Schootemeijer & Langer, 2018). The HD limit has
consequences not only for our understanding of the evolution of the
progenitors of Wolf-Rayet (WR) stars and black holes (BHs), but also for
estimates of the likelihood of binary interaction in the upper-mass end of
stars. The reliability of our predictions of gravitational-wave (GW) events
are thus severely limited as long as the HD limit has not been sufficiently
understood.

|
---|---
Figure 1: HRD positions of populations of massive stars in the SMC (left
panel) and LMC (right panel), based on analyses of apparently single and
binary WR stars (Hainich et al., 2014; Hainich et al., 2015; Shenar et al.,
2016, 2019; Crowther et al., 2002; Tramper et al., 2015), YSGs (Neugent et
al., 2010, 2012), RSGs (Davies et al., 2018; Neugent et al., 2012), LBVs
(Humphreys et al., 2016; Kalari et al., 2018), BSGs (Trundle et al., 2004;
Trundle & Lennon, 2005; Urbaneja et al., 2017), and populations of OB-type
stars (Ramachandran et al., 2018a, b, 2019; Schneider et al., 2018; Hunter et
al., 2008; Castro et al., 2018; Dufton et al., 2019; Bouret et al., 2013). The
typical error bar is shown in the top left corner. Evolutionary tracks (black
solid lines), accounting for rotation with
$\varv_{\mathrm{rot,init}}\,\approx$ 100 km s-1 for the LMC and $\approx$ 180
km s-1 for the SMC (Brott et al., 2011; Köhler et al., 2015), are labeled with
their initial mass. The red dashed lines mark the observational limit reported
by Humphreys & Davidson (1979). The diagrams are largely complete for
supergiants with $T_{\rm eff}\lesssim 12.5\,\mathrm{kK}$, but not necessarily
complete for hotter stars.

It is commonly assumed that the HD limit is a result of radiative instability,
tracing a region in which stars reach their Eddington luminosity and become
LBVs (e.g., Lamers & Fitzpatrick, 1988; Glatzel & Kiriakidis, 1993). However,
attempts to prove this by means of direct radiative transfer calculations in
stellar models have had only limited success. Ulmer & Fitzpatrick (1998)
showed that the so-called modified Eddington limit mimics the shape of the HD
limit, but that it lies roughly one magnitude above the observed HD limit.
More recently, Sanyal et al. (2017) showed that stars with $M_{\rm i}\gtrsim
30\,\mathrm{M}_{\odot}$ and a solar metallicity content ($Z=Z_{\odot}$) reach
the Eddington limit in their interiors and undergo envelope inflation.
However, both these studies and others clearly predict that $L_{\rm max}$
should grow with decreasing $Z$. This is because the opacities in the interior
of stars are strongly correlated with the content of metals within them. The
lower $Z$ is, the weaker the outward radiative pressure becomes. Hence, from
this perspective, environments with a lower ambient $Z$ are expected to allow
the stability of more massive stars, and higher $L_{\rm max}$.

This prediction does not seem to be confirmed by observations. Humphreys &
Davidson (1979) originally reported a maximum RSG luminosity at around
$\log_{10}\left(L/\mathrm{L}_{\odot}\right)\approx 5.7$, later revised
slightly downwards to $5.6$ by Levesque et al. (2005). Similar values are
reported by Massey & Evans (2016) and Drout, Massey & Meynet (2012) for the
super-solar metallicity environment of the Andromeda galaxy (M31) and the
approximately solar $Z$-environment of the Triangulum galaxy (M33).

This trend continues with the Small and Large Magellanic Clouds (SMC, LMC),
which have metallicities of $Z\approx 0.15,0.4\,Z_{\odot}$, respectively (Korn
et al., 2000; Hunter et al., 2007; Trundle et al., 2007). Figure 1 shows the
HRD positions of analysed massive stars in the SMC and LMC, adapted from
Ramachandran et al. (2018b, 2019). These populations are thought to be
complete for the WRs as well as the luminous
($\log_{10}\left(L/\mathrm{L}_{\odot}\right)\gtrsim 4.7$) YSGs and RSGs, but
far from complete for the OB-type stars. Figure 1 illustrates strongly the
absence of RSGs with $\log_{10}\left(L/\mathrm{L}_{\odot}\right)\gtrsim 5.5$
in both galaxies, with the exception of two peculiar objects in the LMC: HD
33579 (Wolf, 1972) and WOH G 64 (Ohnaka et al., 2008).

Davies et al. (2018) performed a statistical comparison of analyzed RSG
populations in the SMC and LMC, and showed compelling evidence that
$\log_{10}\left(L_{\rm max}/\mathrm{L}_{\odot}\right)\lesssim 5.5$ in both
galaxies. Moreover, no evidence is found for higher $L_{\rm max}$ in the SMC
compared to the LMC. Davies et al. (2018) illustrated why this is very likely
a genuine physical fact rather than a statistical one, and how this stands in
tension with recent rotating and non-rotating Geneva evolution models (Ekström
et al., 2012; Georgy et al., 2013).

The challenge of explaining the HD limit solely by radiative instabilities
might imply that other processes are involved. It was already recognized in
the past that rotationally induced mixing can strongly hinder the redward
evolution of massive stars towards the RSG phase (e.g., Maeder & Meynet,
2000). However, initial rotation rates in the excess of $\approx
200\,\mathrm{km}\,\mathrm{s}^{-1}$ (at which strong rotational mixing starts)
seem to be inconsistent with observed rotation rates of massive stars
(Ramírez-Agudelo et al., 2015, 2017; Ramachandran et al., 2019), and models
using realistic rotation rates do not reproduce the apparently $Z$-independent
HD limit (Davies et al., 2018).

Higgins & Vink (2020) recently explored this problem by investigating the
impact of mixing parameters on RSG stellar models. They found that enhanced
semi-convection can reproduce the observed HD limit. However, a solution to
the HD limit problem needs to consider YSGs and cool BSGs as well, as we argue
in our study.

In this paper we tackle the puzzle of the HD limit from a different angle,
quantifying the duration spent beyond the HD limit and the expected number of
stars for various stellar evolution models. For this purpose we study models
with enhanced mixing and construct synthetic populations based on these
models. Moreover, unlike previous studies that considered RSGs alone, we
consider simultaneously RSGs, YSGs, and cool BSGs ($T_{\mathrm{eff}}\lesssim
12.5\,\mathrm{kK}$). As our study shows, considering CSGs as a whole is vital
for correctly assessing the discrepancy between models and observations.

The paper is organized as follows. In Section 2 we describe our numerical
approach. In Section 3 we show the effect of enhanced mixing on stellar
evolution tracks. In Section 4 we make comparisons between observations and
synthetic populations that we generate from our stellar evolution tracks. In
Section 5 we examine the impact of modelling assumptions. We discuss our
results in relation to previous studies, and summarise, in Section 6.

##  2 Numerical method

We use the Modules for Experiments in Stellar Astrophysics code (MESA, version
10398, Paxton et al. 2011; Paxton et al. 2013, 2015, 2018) to evolve stellar
models with $39$ different zero-age main sequence (ZAMS) masses between
$M_{\mathrm{ZAMS}}=4\,\mathrm{M}_{\odot}$ and
$M_{\mathrm{ZAMS}}=107\,\mathrm{M}_{\odot}$. Three different initial
compositions are used, appropriate for SMC, LMC and Milky Way (MW) stars, as
listed in Table 1 (following Hainich et al. 2019).

Table 1: Initial compositions in terms of mass fractions for stellar evolution calculations. | SMC | LMC | MW
---|---|---|---
H | $0.75328$ | $0.74392$ | $0.72031$
He | $0.24448$ | $0.25072$ | $0.26646$
C | $0.00039621$ | $0.00094807$ | $0.0023401$
N | $0.00011606$ | $0.00027772$ | $0.00068549$
O | $0.00096042$ | $0.0022981$ | $0.0056725$
Ne | $0.00021246$ | $0.00050838$ | $0.0012548$
Mg | $0.00011456$ | $0.00027412$ | $0.00067661$
Si | $0.00012179$ | $0.00029144$ | $0.00071934$
S | $0.000056647$ | $0.00013555$ | $0.00033457$
Ar | $0.000013447$ | $0.000032178$ | $0.000079424$
Ca | $0.00001175$ | $0.000028117$ | $0.0000694$
Fe | $0.00023666$ | $0.0005663$ | $0.0013978$
$Z$ | $0.00224$ | $0.00536$ | $0.01323$

###  2.1 Microphysics

The equation of state employed by MESA is a blend of the following equations
of state: OPAL (Rogers & Nayfonov, 2002), SCVH (Saumon, Chabrier & van Horn,
1995), HELM (Timmes & Swesty, 2000), and PC (Potekhin & Chabrier, 2010).
Radiative opacities are primarily from OPAL (Iglesias & Rogers, 1993, 1996),
with low-temperature data from Ferguson et al. (2005) and the high-
temperature, Compton-scattering dominated regime according to Buchler & Yueh
(1976). Electron conduction opacities follow Cassisi et al. (2007).

We use the built-in MESA nuclear reaction network approx21. The Joint
Institute for Nuclear Astrophysics REACLIB reaction rates (Cyburt et al.,
2010) are used, with additional tabulated weak reaction rates (Fuller, Fowler
& Newman 1985; Oda et al. 1994; Langanke & Martínez-Pinedo 2000) and screening
via the prescriptions of Salpeter (1954), Dewitt, Graboske & Cooper (1973),
Alastuey & Jancovici (1978) and Itoh et al. (1979). The formulae of Itoh et
al. (1996) are used for thermal neutrino loss rates.

###  2.2 Wind mass loss

For cool phases ($T_{\mathrm{eff}}\leq 10\,000\,\mathrm{K}$), the mass-loss
prescription of de Jager, Nieuwenhuijzen & van der Hucht (1988) is employed,
regardless of the surface hydrogen mass fraction $X_{\mathrm{s}}$. Wind mass
loss is according to Vink, de Koter & Lamers (2001) for hot
($T_{\mathrm{eff}}\geq 11\,000\,\mathrm{K}$) hydrogen-rich phases
($X_{\mathrm{s}}\geq 0.4$) of the evolution, most notably the main sequence.
If the surface hydrogen mass fraction is low but non-negligible ($0.1\leq
X_{\mathrm{s}}<0.4$) then the empirical mass-loss rate relation of Nugis &
Lamers (2000) is used. For hot hydrogen-deficient models
($X_{\mathrm{s}}<0.1$) we follow Yoon (2017) and Woosley (2019), whose
prescriptions are based on the mass-loss rates of Hainich et al. (2014) and
Tramper, Sana & de Koter (2016). For
$10\,000\,\mathrm{K}<T_{\mathrm{eff}}<11\,000\,\mathrm{K}$ we interpolate
between the two regimes described above. The mass-loss rate is multiplied by a
factor $\eta_{\mathrm{w}}$. We mostly use $\eta_{\mathrm{w}}=1$, but in one
set of models we explore the impact of boosted mass loss on our results by
taking $\eta_{\mathrm{w}}=2$ (Section 5.1) throughout the entire evolution. We
note that the metallicity dependence of the mass-loss rate during the main
sequence, which follows Vink et al. (2001), is $\dot{M}\propto Z^{0.85}$. The
fit by de Jager et al. (1988), which prevails during most of the CSG phases,
has no implicit dependence on metallicity.

###  2.3 Mixing and rotation

The models have initial rotation velocities at the equator of
$V_{\mathrm{i}}=100$, $200$, and $300\,\mathrm{km}\,\mathrm{s}^{-1}$. The
shellular approximation where the angular velocity is constant over isobars
(Meynet & Maeder, 1997) is used for rotating models in MESA (Paxton et al.,
2013). Mixing processes induced by rotation are implemented in a diffusion
approximation (Paxton et al., 2013) and principally follow Heger, Langer &
Woosley (2000). Transport of angular momentum and chemical mixing caused by
internal magnetic fields follows the new prescription of Fuller, Piro & Jermyn
(2019).

The efficiency of rotational mixing is set by two parameters, the ratio of the
turbulent viscosity to the diffusion coefficient $f_{c}$ (input parameter
am_D_mix_factor in MESA), and the ratio between the actual molecular weight
gradient and the value used for computing the mixing coefficients $f_{\mu}$
(input parameter am_gradmu_factor in MESA). These parameters are calibrated to
give the observed nitrogen enhancement for evolved stars. We use the
calibration of Heger et al. (2000), $f_{c}=1/30$ and $f_{\mu}=0.05$ for most
models. However, as discussed by Chieffi & Limongi (2013), this calibration is
not unique (also Potter, Tout & Eldridge 2012). In one set of models we use
alternative values of $f_{c}=0.2$ and $f_{\mu}=1$ (Chieffi & Limongi, 2013).

Convective regions are defined by the Ledoux criterion (except for one set of
models which employs the Schwarzschild criterion) and treated according to
Henyey, Vardya & Bodenheimer (1965) with a mixing-length parameter of
$\alpha_{\mathrm{MLT}}=1.5$. Semiconvective mixing in regions which are Ledoux
stable but Schwarzschild unstable follows Langer, Fricke & Sugimoto (1983)
with an efficiency parameter $\alpha_{\mathrm{sc}}=1$ or
$\alpha_{\mathrm{sc}}=100$. Thermohaline mixing is according to Kippenhahn,
Ruschenplatt & Thomas (1980) with an efficiency parameter of
$\alpha_{\mathrm{th}}=1$. We use the so-called MLT++ implementation for
efficient energy transport in convective regions (Paxton et al., 2013).

Mixing above the convective core boundary is extended in two approaches.
First, for most models, a step overshoot approach is taken (e.g. Shaviv &
Salpeter 1973; Maeder & Meynet 1987), where the convective region is extended
by a fraction $\alpha_{\mathrm{ov}}$ of the pressure scale height $H_{P}$. The
lowest value for $\alpha_{\mathrm{ov}}$ that we use is
$\alpha_{\mathrm{ov}}=0.1$ (e.g. Ekström et al. 2012; Higgins & Vink 2020),
followed by $\alpha_{\mathrm{ov}}=0.335$ as calibrated by Brott et al. (2011),
with higher values starting at $\alpha_{\mathrm{ov}}=0.5$ (e.g. Vink et al.
2010; Higgins & Vink 2019), and additional higher values of
$\alpha_{\mathrm{ov}}=0.8$, $\alpha_{\mathrm{ov}}=1$ and
$\alpha_{\mathrm{ov}}=1.2$111These parameter values might not be physical, but
rather used as proxy to investigate the issue at hand.. In a second approach
we use an exponential core overshooting (e.g. Herwig 2000), where the mixing
efficiency decays smoothly outside the core, rather than dropping abruptly as
in the step overshoot approach. We follow the prescription of Jermyn, Tout &
Chitre (2018) to compute the fraction $f_{\mathrm{ov}}$ used for the decay
scale,

$\displaystyle f_{\mathrm{ov}}=\frac{2}{\ln
150-2\ln\left(v_{\mathrm{c}}/c_{\mathrm{s}}\right)+\left(\frac{5}{2}\left(H_{P}/r_{\mathrm{c}}\right)-1\right)},$
(1)

where $v_{\mathrm{c}}/c_{\mathrm{s}}$ is the mass-averaged convective core
Mach number, $H_{P}$ is taken at the top of the core, and $r_{\mathrm{c}}$ is
the convective core radius. The value of $f_{\mathrm{ov}}$ is updated
according to equation (1) at the end of every evolution step. This
prescription corresponds to a meridional circulation driven by anisotropy in
the heat flux emerging from the convective core. The anisotropy is
rotationally-induced, but because the convective turnover time is so long in
the core this effect saturates at slower angular velocities than any we
consider here.

##  3 Stellar evolution tracks

A total of $5967$ models were evolved, for $39$ initial masses
($M_{\mathrm{ZAMS}}=4$, $5$, $6$, $7$, $8$, $9$, $10$, $11$, $12$, $13$, $14$,
$15$, $16$, $17$, $19$, $20$, $22$, $23$, $25$, $27$, $29$, $31$, $33$, $36$,
$39$, $42$, $45$, $48$, $52$, $56$, $60$, $64$, $69$, $74$, $80$, $86$, $92$,
$99$, and $107\,\mathrm{M}_{\odot}$), $3$ initial rotation velocities
($V_{\mathrm{i}}=100$, $200$, and $300\,\mathrm{km}\,\mathrm{s}^{-1}$), $3$
compositions (SMC, LMC, and MW), and $17$ sets of modelling assumptions as
detailed in Table 2. The evolution reached the end of core helium burning for
the lower masses ($M_{\mathrm{ZAMS}}<10\,\mathrm{M}_{\odot}$) or core carbon
burning for the higher masses ($M_{\mathrm{ZAMS}}\geq
10\,\mathrm{M}_{\odot}$). We present several sub-sets of the results to
highlight our main findings.

Table 2: Modelling assumptions for all sets of models computed. # | stability | overshooting | $\alpha_{\mathrm{sc}}$ | $f_{c}$ | $f_{\mu}$ | $\eta_{\mathrm{w}}$
---|---|---|---|---|---|---
| criterion |  |  |  |  |
$1$ | Ledoux | $\alpha_{\mathrm{ov}}=0.1$ | $1$ | $1/30$ | $0.05$ | $1$
$2$ | Ledoux | $\alpha_{\mathrm{ov}}=0.335$ | $1$ | $1/30$ | $0.05$ | $1$
$3$ | Ledoux | $\alpha_{\mathrm{ov}}=0.5$ | $1$ | $1/30$ | $0.05$ | $1$
$4$ | Ledoux | $\alpha_{\mathrm{ov}}=0.8$ | $1$ | $1/30$ | $0.05$ | $1$
$5$ | Ledoux | $\alpha_{\mathrm{ov}}=1.0$ | $1$ | $1/30$ | $0.05$ | $1$
$6$ | Ledoux | $\alpha_{\mathrm{ov}}=1.2$ | $1$ | $1/30$ | $0.05$ | $1$
$7$ | Ledoux | $f_{\mathrm{ov,JTC}}$ | $1$ | $1/30$ | $0.05$ | $1$
$8$ | Ledoux | $\alpha_{\mathrm{ov}}=0.1$ | $100$ | $1/30$ | $0.05$ | $1$
$9$ | Ledoux | $\alpha_{\mathrm{ov}}=0.335$ | $100$ | $1/30$ | $0.05$ | $1$
$10$ | Ledoux | $\alpha_{\mathrm{ov}}=0.5$ | $100$ | $1/30$ | $0.05$ | $1$
$11$ | Ledoux | $\alpha_{\mathrm{ov}}=0.8$ | $100$ | $1/30$ | $0.05$ | $1$
$12$ | Ledoux | $\alpha_{\mathrm{ov}}=1.0$ | $100$ | $1/30$ | $0.05$ | $1$
$13$ | Ledoux | $\alpha_{\mathrm{ov}}=1.2$ | $100$ | $1/30$ | $0.05$ | $1$
$14$ | Ledoux | $f_{\mathrm{ov,JTC}}$ | $100$ | $1/30$ | $0.05$ | $1$
$15$ | Schwarzschild | $\alpha_{\mathrm{ov}}=0.335$ | $-$ | $1/30$ | $0.05$ | $1$
$16$ | Ledoux | $\alpha_{\mathrm{ov}}=0.335$ | $100$ | $0.2$ | $1.0$ | $1$
$17$ | Ledoux | $\alpha_{\mathrm{ov}}=0.335$ | $100$ | $1/30$ | $0.05$ | $2$
|
---|---
|
Figure 2: Hertzsprung-Russell diagrams for models with a semiconvective mixing efficiency of $\alpha_{\mathrm{sc}}=100$ and an overshoot parameter of $\alpha_{\mathrm{ov}}=0.335$ (left) and $\alpha_{\mathrm{ov}}=1.2$ (right) for LMC (top) and MW (bottom) initial compositions, with an initial rotation velocity of $200\,\mathrm{km}\,\mathrm{s}^{-1}$. Initial masses are in the range $4\,\mathrm{M}_{\odot}\leq M_{\mathrm{ZAMS}}\leq 107\,\mathrm{M}_{\odot}$. Models every $50\,000\,\mathrm{yr}$ are marked. The thick magenta line maps the HD limit. |
---|---
|
Figure 3:  Hertzsprung-Russell diagrams for models with a semiconvective
mixing efficiency of $\alpha_{\mathrm{sc}}=100$ and an overshoot parameter of
$\alpha_{\mathrm{ov}}=0.335$ (left) and $\alpha_{\mathrm{ov}}=1.2$ (right) for
LMC (top) and MW (bottom) initial compositions, with an initial rotation
velocity of $100\,\mathrm{km}\,\mathrm{s}^{-1}$. Initial masses are in the
range $4\,\mathrm{M}_{\odot}\leq M_{\mathrm{ZAMS}}\leq
107\,\mathrm{M}_{\odot}$. The line colour follows the surface Eddington
factor. Black crosses show the points in the evolution where core hydrogen
burning ends, and red asterisks mark the end of the evolution, when core
carbon burning ends. The thick magenta line maps the HD limit.

Figure 2 shows evolutionary tracks for LMC and MW compositions and an initial
rotation velocity of $V_{\mathrm{i}}=200\,\mathrm{km}\,\mathrm{s}^{-1}$, for
$\alpha_{\mathrm{sc}}=100$ and step overshooting with
$\alpha_{\mathrm{ov}}=0.335$ and $\alpha_{\mathrm{ov}}=1.2$. Points at
intervals of $50\,000\,\mathrm{yr}$ along the evolution are marked with dots,
to illustrate the duration of different phases, and the relative expected
number of stars at each part of the evolution for each initial mass. It can be
seen that for $\alpha_{\mathrm{ov}}=0.335$ numerous models spend a long time
beyond the HD limit, while for $\alpha_{\mathrm{ov}}=1.2$ the evolutionary
tracks almost do not cross the limit, and then only for a brief time if they
do. The value of $\alpha_{\mathrm{ov}}=1.2$ is probably excessive, as a gap is
present between the tracks and the diagonal part of the HD limit. We do not
suggest that $\alpha_{\mathrm{ov}}=1.2$ is a reasonable choice, but rather use
it to demonstrate the strong effect of overshooting. To test the models a
statistical analysis is discussed in Section 4.

Figure 3 shows evolutionary tracks for LMC and MW compositions and an initial
rotation velocity of $V_{\mathrm{i}}=100\,\mathrm{km}\,\mathrm{s}^{-1}$ (two
times slower compared to the models in Figure 2), for the same mixing
parameters as in Figure 2 ($\alpha_{\mathrm{sc}}=100$, and
$\alpha_{\mathrm{ov}}=0.335$ or $\alpha_{\mathrm{ov}}=1.2$), with the colour
along the tracks showing the ratio between the stellar luminosity and the
Eddington luminosity $\Gamma_{\mathrm{Edd}}$ (computed by MESA, taking into
account the gas opacity). Compared to the tracks in Figure 2, the models with
$\alpha_{\mathrm{ov}}=1.2$ in Figure 3 evolve somewhat further beyond the HD
limit, because of the slightly reduced rotational mixing. We point out in
Figure 3 the possibility of stellar models to reside beyond the HD limit while
not exceeding their Eddington luminosity. Also marked are the terminal age
main sequence (TAMS222The TAMS is defined as the point where the central
hydrogen mass fraction drops below $0.01$.) and final positions on the
HRD333The final position is marked only for models which ignited carbon in
their centre, and are therefore considered to be CCSN progenitors.. For
$\alpha_{\mathrm{ov}}=0.335$, the TAMS location is beyond the HD limit for the
highest initial masses. For $\alpha_{\mathrm{ov}}=1.2$, the TAMS moves redward
for the lower masses, but for the higher masses it can move bluewards. There
are less pre-SN red supergiants for $\alpha_{\mathrm{ov}}=1.2$, though for
such a large extent of overshooting lower initial masses might give rise to
CCSNe whose progenitors are red supergiants.

The duration of the evolutionary stage when a star is located on the HRD above
the HD limit is shown in Figure 4 for various parameters. For each
composition, initial mass, and overshooting scheme, there are $3$ different
initial rotation velocities, and we take the longest time from the three
tracks. For $\alpha_{\mathrm{ov}}=0.335$, massive stars spend $\ga
100\,000\,\mathrm{yr}$ beyond the HD limit. The highest masses with MW
composition do not cross the HD limit thanks to mass removal by winds (besides
the models with $\alpha_{\mathrm{ov}}=0.1$), but as the rate depends on the
metallicity, this effect is not present in the SMC and LMC models, except for
those with significantly enhanced overshooting. With increasing
$\alpha_{\mathrm{ov}}$, the LMC and MW models spend less time beyond the HD
limit.

---
Figure 4: Time spent beyond the HD limit as function of ZAMS mass for SMC
(top), LMC (middle) and MW (bottom) composition, for different core
overshooting prescriptions. For each initial mass, the maximal value out of
all initial rotation rates is taken. Points labeled as $f_{\mathrm{ov,JTC}}$
use the exponential overshooting coefficient $f_{\mathrm{ov}}$ as function of
core properties described by Jermyn et al. 2018.

For the SMC the models with $30\la M_{\mathrm{ZAMS}}/\mathrm{M}_{\odot}\la 80$
always cross the HD limit and spend a long time beyond it. This is because
crossing the HD limit is a result of the combination of mixing and mass loss.
With increased mixing, envelope material is used as fuel in the core, both
increasing the core mass and luminosity and therefore the mass loss, while at
the same time decreasing the envelope mass which needs to be removed by winds.

The models which use the overshooting prescription of Jermyn et al. (2018)
give quantitatively similar results to those with $\alpha_{\mathrm{ov}}=1$,
with $f_{\mathrm{ov}}$ varying slightly with time and initial mass but
generally close to $f_{\mathrm{ov}}\approx 0.1$ as shown in Figure 5. Jermyn
et al. (2018) motivate their prescription as describing enhanced rotational
mixing caused by core anisotropy. In that sense there is a physical motivation
for enhanced mixing. However, Jermyn et al. (2018) caution that their
prescription appears to overestimate $f_{\mathrm{ov}}$ relative to what Claret
& Torres (2017) infer from observations. Furthermore, while the usage of
$f_{\mathrm{ov}}$ is convenient, their mixing mechanism does not exhibit an
exponentially decaying geometry like the $f_{\mathrm{ov}}$ prescription
assumes and so this implementation does not accurately describe the physics.
Our results imply that there is a strong motivation to further investigate the
mixing mechanism proposed by Jermyn et al. (2018).

Figure 5: Exponential decay scale as fraction of the pressure scale height as
function of time for several models with MW composition and an initial
rotation velocity of $V_{\mathrm{i}}=100\,\mathrm{km}\,\mathrm{s}^{-1}$ as
calculated by using equation (1).

##  4 Statistical comparisons between synthetic populations and observations

###  4.1 Our sample compilation

We compile the luminosities and temperatures of all known CSGs in the
Magellanic Clouds with a luminosity of
$\log_{10}\left(L/\mathrm{L}_{\odot}\right)=4.7$ or higher. Unlike Davies et
al. (2018), who only considered RSGs, we consider RSGs, YSGs, and cool BSGs.
The reason is that the HD limit stretches over a large region on the HRD
populated by these various spectral types. For simplicity, we only consider
here the horizontal part of the HD limit, which extends up to effective
temperatures of $\approx 12.5\,\mathrm{kK}$. Following standard convention
(e.g., Drout et al., 2012), we define RSGs, YSGs, and cool BSGs in the
temperature regimes $T_{\rm eff}\leq 4\,800\,$K, $4\,800\,{\rm K}<T_{\rm
eff}<7\,500\,{\rm K}$, and $7\,500\leq T_{\rm eff}\leq 12\,500\,{\rm K}$,
respectively.

For the LMC, we cross-match the RSG list of Davies et al. (2018) with the RSG-
YSG list of Neugent et al. (2012). For targets that appear in both
compilations, we adopt temperatures and luminosities from Davies et al.
(2018), who derived these parameters from a complete spectral energy
distribution (SED) fitting. When temperatures are not specified, we use
calibrations between spectral types and temperatures by Tabernero et al.
(2018) to derive the temperature. Since both studies claim to be complete for
$\log_{10}\left(L/\mathrm{L}_{\odot}\right)\geq 4.7$, we only include objects
exceeding this threshold.

The catalogue of Neugent et al. (2012) also includes stars hotter than
$7\,500\,\mathrm{K}$, which are considered BSGs. Since their study is not
necessarily complete for BSGs, we extended our catalogue by retrieving the
Gaia DR2 catalogue (Gaia Collaboration, 2018) centered on the LMC with a
search radius of $5.5$ degrees. To identify cool BSG candidates, we filtered
all LMC stars in the Gaia DR2 catalogue with Gaia fluxes fulfilling the
criteria $-0.1<Bp-Rp<0.6\,$mag and $G<15.5\,$mag, accounting for a typical
reddening value of $E_{B-V}=0.09$ (e.g. Fitzpatrick & Garmany, 1990). We then
cross-matched our list with the SIMBAD catalogue to retrieve spectral types
using the Vizier X-match service. Main references are Sanduleak (1970),
Ardeberg et al. (1972), Stock et al. (1976), Evans et al. (2006), and Urbaneja
et al. (2017). All identified targets were classified before and have spectral
types consistent with BSGs, and the majority of those were included in
previous spectroscopic analyses of CSGs in the LMC.

All stars with spectral types earlier than B7 in our final list were removed,
including a few WR stars. Finally We included four LBVs from the compilation
given by Smith (2019).

For all remaining objects, we extracted radial velocities (RVs) and proper
motions (PMs) from the SIMBAD database (Wenger et al., 2000). The sources of
the RVs were predominantly the Gaia DR2 catalogue (Gaia Collaboration, 2018),
Massey & Olsen (2003), Neugent et al. (2012), Fehrenbach (1972), and
Fehrenbach & Duflot (1982). PMs originate in the Gaia DR2 catalogue for all
sources but ten, for which they are retrieved from Gaia DR1 (Lindegren et al.,
2016). The mean PM is $1.78\,{\rm mas}\,{\rm yr}^{-1}$ with a standard
deviation of $0.3\,{\rm mas}\,{\rm yr}^{-1}$, which reflects the measurement
limit of Gaia. There are $15$ outliers with PMs larger than $4\,{\rm
mas}\,{\rm yr}^{-1}$ within their respective errors, and they are omitted from
our sample to ensure that we do not include foreground Galactic objects.

When available, $T_{\rm eff}$ and $\log L$ values for the cool BSGs were
adopted from Urbaneja et al. (2017) and Smith (2019). Otherwise, we used
spectral-type calibrations by Fitzpatrick & Garmany (1990) to derive the
effective temperatures and estimated the extinction parameters based on the
expected intrinsic colours. We then used bolometric corrections following
Flower (1996) and Torres (2010), assuming a distance of $49.97\,\mathrm{kpc}$
(Pietrzyński et al., 2013). The final list for the LMC comprises $375$ stars:
$265$ RSGs, $39$ YSGs, and $71$ cool BSGs (four of which are LBVs).

For the SMC, we repeat this procedure using the RSGs listed by Davies et al.
(2018) and the YSGs listed by Neugent et al. (2010). The RVs, PMs, and
spectral types are again extracted using SIMBAD. The RVs and PMs originate
predominantly from the Gaia DR2 catalogue, but also from Massey & Olsen
(2003), Neugent et al. (2010), and González-Fernández et al. (2015). The
spectral types are retrieved from Feast et al. (1960), Dubois et al. (1977),
Humphreys (1983), Lennon (1997), and Dufton et al. (2000). Again, all objects
earlier than B7 are removed. We identify about ten outliers with
$\mathrm{PM}>4\,{\rm mas}\,{\rm yr}^{-1}$, which are removed from our sample.
When available, $T_{\rm eff}$ and $\log L$ values for the cool BSGs were
adopted from Dufton et al. (2000). Otherwise, we use calibrations by Evans &
Howarth (2003) to derive the effective temperatures, make the same assumptions
as Neugent et al. (2010) regarding the reddening and the distance towards the
SMC, and use the same relations as above to derive the bolometric corrections
and luminosities. The final list comprises $179$ stars: $140$ RSGs, $7$ YSGs,
and $32$ cool BSGs.

We note that accurate derivations of $T_{\rm eff}$ and $\log L$ should rely on
the fitting of SEDs. However, given the statistical nature of our study, the
calibrations used above should be sufficient for our purpose (e.g. Neugent et
al., 2010, 2012).

###  4.2 Population synthesis

We construct synthetic populations by generating random initial masses
according to a Salpeter IMF. The initial rotation velocity is chosen according
to the observed distributions for the LMC (Ramachandran et al., 2018b) and the
SMC (Ramachandran et al., 2019). For both the initial mass and velocity, the
nearest values available in our models are used to chose a stellar evolution
track to follow, rather than interpolating between tracks. This results in
$58\%$ of SMC models being assigned an initial rotation velocity of
$100\,\mathrm{km}\,\mathrm{s}^{-1}$, $31\%$ getting
$200\,\mathrm{km}\,\mathrm{s}^{-1}$ and $11\%$ with
$300\,\mathrm{km}\,\mathrm{s}^{-1}$. For the LMC the corresponding percentages
are $80\%$, $19\%$ and $1\%$.

The stellar age is chosen according to a uniform distribution, corresponding
to a constant star-formation rate (SFR)444This is a reasonable assumption as
we are interested in young stars, though in general the SFR is not constant in
the SMC and the LMC (Section 5.4).. The stellar properties are interpolated
from the evolutionary track according to the generated stellar age. If the
generated age is longer than the lifetime of the computed stellar evolution
track, the star is discarded.

Random stars are generated until the combined number of RSGs and YSGs
($T_{\mathrm{eff}}<7500\,\mathrm{K}$) in the luminosity interval
$4.7\leq\log_{10}(L/\mathrm{L}_{\odot})<5.2$ matches the observed number –
this is our normalisation. For the SMC this gives $134$ red and yellow
supergiants, and for the LMC $274$. An example of a synthetic population is
presented in Figure 6.

Figure 6: Example synthetic populations of stars with SMC (top) and LMC
(bottom) initial compositions.

###  4.3 Observations vs. simulations

For each of the $17$ sets of modelling assumptions (Table 2) we generate $25$
random realisations of such populations to get an error estimate for the
computed numbers. We count the number of stars generated which are over-
luminous, i.e. those with $\log_{10}\left(L/\mathrm{L}_{\odot}\right)\geq
5.4$. The results for the $850$ synthetic populations are summarised in Table
3 and Table 4, and in Figure 7 we show the number of over-luminous supergiants
in our synthetic populations for all models which employ step overshooting.
Tables 3 and 4 show that looking only at red supergiants leads to
underestimating the excess of over-luminous stellar models. The excess
generally decreases with increasing $\alpha_{\mathrm{ov}}$, although for
efficient semiconvective mixing ($\alpha_{\mathrm{sc}}=100$) the behaviour is
non-monotonic for the red and yellow regimes. The models employing the
overshooting prescription given by equation (1) produce very similar results
to those with $\alpha_{\mathrm{ov}}=1$.

|
---|---
|
|
Figure 7: Number of supergiant stars with $\log_{10}\left(L/\mathrm{L}_{\odot}\right)\geq 5.4$ in synthetic populations as function of the overshooting extent, for two values of the semiconvective mixing efficiency, as well as models employing the Schwarzschild stability criterion instead of Ledoux, for the SMC (left) and the LMC (right). For $\alpha_{\mathrm{sc}}=100$, models with doubled mass loss or alternative parameters for rotational mixing are also shown. The top panels show the number of over-luminous red ($T_{\mathrm{eff}}\leq 4\,800\,\mathrm{K}$) supergiant stars, the middle panels show red and yellow ($T_{\mathrm{eff}}<7\,500\,\mathrm{K}$) supergiant stars, and the bottom panels show cool ($T_{\mathrm{eff}}\leq 12\,500\,\mathrm{K}$) supergiant stars. Table 3: Excess of over-luminous stars in synthetic populations. | overshooting | $\alpha_{\mathrm{sc}}$ | $N_{\mathrm{excess,RSG+YSG}}$ | $N_{\mathrm{excess,RSG}}$
---|---|---|---|---
SMC | $\alpha_{\mathrm{ov}}=0.1$ | $1$ | $152\pm 15$ | $79\pm 10$
SMC | $\alpha_{\mathrm{ov}}=0.335$ | $1$ | $20\pm 5$ | $13\pm 4$
SMC | $\alpha_{\mathrm{ov}}=0.5$ | $1$ | $17\pm 5$ | $10\pm 4$
SMC | $\alpha_{\mathrm{ov}}=0.8$ | $1$ | $11\pm 3$ | $6\pm 2$
SMC | $\alpha_{\mathrm{ov}}=1.0$ | $1$ | $7\pm 2$ | $4\pm 2$
SMC | $\alpha_{\mathrm{ov}}=1.2$ | $1$ | $5\pm 2$ | $2\pm 1$
SMC | $f_{\mathrm{ov,JTC}}$ | $1$ | $8\pm 3$ | $4\pm 3$
SMC | $\alpha_{\mathrm{ov}}=0.1$ | $100$ | $9\pm 3$ | $3\pm 2$
SMC | $\alpha_{\mathrm{ov}}=0.335$ | $100$ | $10\pm 3$ | $3\pm 1$
SMC | $\alpha_{\mathrm{ov}}=0.5$ | $100$ | $18\pm 5$ | $7\pm 3$
SMC | $\alpha_{\mathrm{ov}}=0.8$ | $100$ | $10\pm 3$ | $5\pm 2$
SMC | $\alpha_{\mathrm{ov}}=1.0$ | $100$ | $8\pm 4$ | $4\pm 3$
SMC | $\alpha_{\mathrm{ov}}=1.2$ | $100$ | $4\pm 2$ | $1\pm 1$
SMC | $f_{\mathrm{ov,JTC}}$ | $100$ | $7\pm 2$ | $4\pm 2$
LMC | $\alpha_{\mathrm{ov}}=0.1$ | $1$ | $145\pm 15$ | $103\pm 12$
LMC | $\alpha_{\mathrm{ov}}=0.335$ | $1$ | $41\pm 6$ | $27\pm 5$
LMC | $\alpha_{\mathrm{ov}}=0.5$ | $1$ | $34\pm 7$ | $22\pm 4$
LMC | $\alpha_{\mathrm{ov}}=0.8$ | $1$ | $21\pm 5$ | $12\pm 4$
LMC | $\alpha_{\mathrm{ov}}=1.0$ | $1$ | $14\pm 4$ | $5\pm 2$
LMC | $\alpha_{\mathrm{ov}}=1.2$ | $1$ | $9\pm 4$ | $3\pm 2$
LMC | $f_{\mathrm{ov,JTC}}$ | $1$ | $14\pm 3$ | $6\pm 3$
LMC | $\alpha_{\mathrm{ov}}=0.1$ | $100$ | $33\pm 5$ | $7\pm 3$
LMC | $\alpha_{\mathrm{ov}}=0.335$ | $100$ | $44\pm 7$ | $10\pm 3$
LMC | $\alpha_{\mathrm{ov}}=0.5$ | $100$ | $32\pm 5$ | $16\pm 3$
LMC | $\alpha_{\mathrm{ov}}=0.8$ | $100$ | $20\pm 4$ | $11\pm 3$
LMC | $\alpha_{\mathrm{ov}}=1.0$ | $100$ | $12\pm 4$ | $5\pm 2$
LMC | $\alpha_{\mathrm{ov}}=1.2$ | $100$ | $7\pm 2$ | $2\pm 1$
LMC | $f_{\mathrm{ov,JTC}}$ | $100$ | $14\pm 3$ | $7\pm 3$
Table 4: Excess of over-luminous stars in additional synthetic populations. | model assumptions | $N_{\mathrm{excess,RSG+YSG}}$ | $N_{\mathrm{excess,RSG}}$
---|---|---|---
SMC | Schwarzschild | $4\pm 2$ | $1\pm 1$
SMC |  $f_{c}=0.2$, $f_{\mu}=1$ | $12\pm 4$ | $4\pm 2$
SMC | $\eta_{\mathrm{w}}=2$ | $10\pm 4$ | $2\pm 2$
LMC | Schwarzschild | $21\pm 5$ | $4\pm 2$
LMC |  $f_{c}=0.2$, $f_{\mu}=1$ | $51\pm 7$ | $12\pm 3$
LMC | $\eta_{\mathrm{w}}=2$ | $15\pm 3$ | $3\pm 2$

In Figure 8 we show the simulated luminosity distribution for a few sets of
mixing parameters compared to the observed distribution. For the RSGs, the
simulated distribution is similar to the observed distribution for most mixing
assumptions, with a small excess in the simulations for some cases, such as
the combination of $\alpha_{\mathrm{ov}}=0.335$ and $\alpha_{\mathrm{sc}}=1$.
The excess in the simulations increases when including also the YSGs,
especially as there is a non-negligible number of such stars in the simulated
populations, while the number of observed luminous YSGs is small. When
including also cool BSGs, i.e. all stars with $T_{\mathrm{eff}}\leq
12\,500\,\mathrm{K}$, the excess of over-luminous stars in the simulated
populations becomes even larger. An excess in the lower luminosities for BSGs
is acceptable as we do not expect the sample to be complete for these
temperatures, but higher luminosity stars
($\log_{10}\left(L/\mathrm{L}_{\odot}\right)\geq 5.4$) would definitely be
observed, and therefore the simulations cannot be taken to properly account
for the stellar population.

|
---|---
|
|
Figure 8: Number of supergiant stars as function of luminosity, for a sample
of five different sets of mixing parameters compared to the observed
distributions for the SMC (left) and the LMC (right). The top panels show the
number of RSGs ($T_{\mathrm{eff}}\leq 4\,800\,\mathrm{K}$), the middle panels
show RSGs and YSGs ($T_{\mathrm{eff}}<7\,500\,\mathrm{K}$), and the bottom
panels show cool ($T_{\mathrm{eff}}\leq 12\,500\,\mathrm{K}$) supergiant
stars.

In Figure 9 we show the ratio between the number of RSGs with
$4.7\leq\log_{10}\left(L/\mathrm{L}_{\odot}\right)<5.4$ to the number of YSGs
in the same luminosity range. This ratio is a good test for mixing in evolved
massive stars. Schootemeijer et al. (2019) discuss a similar ratio, but
between the numbers of RSGs and BSGs. We note also that the expansion and
colour of supergiants depends on the treatment of convection, and employing
MLT++ results in very luminous stars being hotter (Klencki et al., 2020). In
the present study we focus on the most luminous stars, whose numbers are small
compared to those with lower luminosities which mostly affect the YSG/RSG (or
BSG/RSG) ratio. Supergiants of all colours need to be taken into account, to
separate the issue of the number ratio of different colours from the issue of
the HD limit.

---
Figure 9: Number of YSG stars relative to the number of RSG stars in synthetic
populations as function of the overshooting extent, for two values of the
semiconvective mixing efficiency, for the SMC (top) and the LMC (bottom).

##  5 Sources of uncertainty and their impact

###  5.1 Mass loss

Alongside mixing, continuous or eruptive mass loss is a key parameter that is
responsible for stripping a star of its H-rich envelope. The larger the mass
loss, the less likely the star is to cross the HD limit. Motivated by the
apparent $Z$-independence of the HD limit, we chose to focus on mixing here.
However, to shed more light on the interplay between mass loss and mixing in
this context, we also provide a set of models with mass-loss rates that are
boosted by a factor of two. Recent theoretical and empirical determinations of
mass-loss rates during the OB and the RSG phases (e.g., Ramachandran et al.,
2019; Björklund et al., 2020; Beasor et al., 2020) rather suggest that
standard prescriptions such as those used here already lead to an
overestimation of the mass loss. Nevertheless, considering the poor
understanding of eruptive mass loss, we explore the impact of an increased
mass loss. We do this for the case of $\alpha_{\mathrm{ov}}=0.335$ and
$\alpha_{\mathrm{sc}}=100$.

The impact of boosting the mass loss throughout the stellar evolution by a
factor of two is shown in Figure 7. As could be anticipated, boosting the mass
loss helps reduce the excess of CSGs above the HD limit. However, even with
increased mass loss, the discrepancy remains significant. Moreover, while it
leads to an appreciable improvement in the case of the LMC, the improvement is
negligible in the SMC, which again opposes the apparent $Z$-independence of
the HD limit. To conclude, while continuous mass loss will play a role in
shaping the HD limit, there is little support that it alone can explain the
observations. It is well possible that a $Z$-independent rapid phase of mass
loss (e.g., in a CSG or LBV phase) needs to be invoked to avoid the
discrepancy between observations and theory, but a consistent physical
framework for implementing it is still missing.

###  5.2 Mixing prescriptions

The main uncertainty in stellar evolution modelling on which we focus in this
work is mixing in stellar interiors. Therefore, in addition to the efficiency
of semiconvection and the extent of overshooting, the effects of a couple of
other mixing assumptions were tested. While we mostly used the Ledoux
criterion for defining the convective boundary, in one set of models (with
$\alpha_{\mathrm{ov}}=0.335$) we used the Schwarzschild criterion instead. The
stability criterion employed affects not only the boundary of the convective
core, but also intermediate regions between the core boundary and the stellar
photosphere and the position of models in the HRD (Georgy, Saio & Meynet,
2014). The analysis of our synthetic populations generated from stellar
evolution tracks which employed the Schwarzschild criterion yields a somewhat
reduced excess of over-luminous CSGs compared to the analysis with the Ledoux
criterion and the same overshooting extent. Similarly to our results with
boosted mass loss (Section 5.1), the excess remains non-negligible, especially
when considering the entire relevant temperature range.

Up until now we have discussed mixing only for convective regions and near
their boundaries. In models of rotating stars, considerable mixing occurs also
in radiative regions, owing to various instabilities (e.g., Heger et al.,
2000). In one set of models we changed the rotational mixing parameters in
radiative regions to $f_{c}=0.2$ and $f_{\mu}=1$, as described in Section 2.3.
Our analysis of the synthetic populations generated from tracks with these
alternative mixing parameters yields almost no change in the CSG excess (Fig.
7). This might be a result of the possible degeneracy of these parameters
(Chieffi & Limongi, 2013).

As an additional test of the sensitivity to the implementation of rotational
mixing, we also compare our results to the Geneva tracks with MW and SMC
compositions (Ekström et al., 2012; Georgy et al., 2013). In Figure 10 we show
the time spent beyond the HD limit for several of our modelling assumptions
and for the Geneva models, for which rotational mixing is treated as an
advective-diffusive process, compared to the diffusion approximation adopted
by MESA. The Geneva models spend time periods beyond the HD limit rather
similar to our tracks computed with MESA. While a quantitative comparison
using synthetic populations will give a more definitive answer, we surmise
that the computed excess of CSGs will be similar with the Geneva tracks.

---
Figure 10: Time spent beyond the HD limit as function of ZAMS mass for SMC
(top) and MW (bottom) composition, for stellar evolution tracks as described
in the inset.

###  5.3 Multiplicity

In principle, the presence of a close companion (orbital period $P\lesssim
3\,\mathrm{yr}$) would inhibit the expansion of its companion. Because a large
fraction of massive stars interact with a companion during their evolution
(Sana et al., 2012; Moe & Di Stefano, 2017), the effects of binaries on the
simulated populations of CSGs and the HD limit should be taken into account.
Given the many uncertain parameters that describe the properties and evolution
of binaries, we decided to not account for multiplicity here. One may naively
expect that binary interactions may be highly relevant for explaining the HD
limit, since they can prevent stars from becoming CSGs. However, because the
total number of RSGs and YSGs is used for the normalisation of our synthetic
populations, and because it itself will be affected by binary interactions,
the effect of stellar multiplicity is not trivial. If the binary separation
shrinks as the binary mass increases, then one may expect binary interactions
to help resolve the discrepancy. However, solid evidence for this at the
upper-mass end is currently lacking. Addressing this question in future
studies should be helpful to quantitatively constrain the impact of
multiplicity on the HD limit.

###  5.4 Star-formation history

There are strong indications that the recent star-formation history in the
Magellanic Clouds, especially in the SMC, was not continuous but is rather
characterised by various peaks over the past $100\,\mathrm{Myr}$ (Indu &
Subramaniam, 2011; Ramachandran et al., 2019; Fulmer et al., 2020). As such,
this can have an important impact on the synthetic populations obtained, and
hence on the final conclusions. However, it is unlikely that star-formation
history alone can explain the discrepancy in the SMC and the LMC
simultaneously. Moreover, the presence of WR stars in both galaxies that span
a substantial luminosity range (Hainich et al., 2014; Hainich et al., 2015;
Shenar et al., 2016, 2019) suggests that a similar distribution of CSGs across
the luminosity range could be anticipated. Given the many uncertainties,
repeating this investigation with detailed star-formation histories is beyond
the scope of our work, but should be explored in future studies.

##  6 Summary and discussion

We evolve numerous grids of stellar evolution tracks, for MW, LMC and SMC
compositions, with a variety of mixing prescriptions. We find that enhanced
mixing diverts stellar evolution tracks from the “forbidden region” defined by
the HD limit (Figs. 2-3). Based upon these grids of stellar models we
construct synthetic populations with initial compositions appropriate for the
LMC and SMC, and make quantitative comparisons with the observed populations
of cool ($T_{\mathrm{eff}}\leq 12.5\,\mathrm{kK}$) supergiant stars in these
galaxies (Section 4). We find that enhanced mixing reduces the excess of over-
luminous CSGs in our simulated stellar populations. While for RSGs there does
not seem to be a severe problem, the tension between observations and
simulations increases with increasing the upper temperature cutoff. We can
therefore consider the existence of a “Cool Supergiant Problem” as the
apparent mismatch between observed supergiants and the results of stellar
evolution models.

When considering only RSGs, we find that the dependence of the excess on the
overshooting extent is non-monotonic for efficient semiconvective mixing
($\alpha_{\mathrm{sc}}=100$). Recent studies by Schootemeijer et al. (2019)
and Higgins & Vink (2020) advocate for highly efficient mixing in regions of
semiconvection to account for the properties of supergiants in the LMC and
SMC. Schootemeijer et al. (2019) also claim that convective overshooting can
be constrained by the properties of the populations, such as the ratio between
RSGs and BSGs. We differ from Schootemeijer et al. (2019) in two regards:
($i$) We consider much higher values of $\alpha_{\mathrm{ov}}$, beyond an
apparent extremum around $\alpha_{\mathrm{ov}}\approx 0.5$; ($ii$)
Schootemeijer et al. (2019) focus on the bulk of the supergiant population,
while we investigate the stars with the highest luminosities, which are a
small fraction of the overall population.

Higgins & Vink (2020) suggest that the HD limit can be explained with
decreased overshoot mixing, $\alpha_{\mathrm{ov}}=0.1$, though they consider
only the red part of the evolution. We show that for red supergiants taking
$\alpha_{\mathrm{ov}}=0.1$ together with $\alpha_{\mathrm{sc}}=100$ indeed
gives a reasonable account of the supergiant populations in the LMC and SMC
(Figs. 7-8). This is explained by the significant effect of efficient
semiconvection in increasing the effective surface temperatures of the stellar
models. Thus, for efficient semiconvection, core He-burning stars tend to
appear as YSGs or BSGs instead of RSGs. However, this does not solve the
supergiant excess beyond the HD limit, but merely “sweeps it under the
carpet”. When we consider the entire temperature range to which the horizontal
HD limit applies ($T_{\rm eff}\lesssim 12.5\,\mathrm{kK}$) we find that there
is a large excess of stellar models in our simulated populations. As Higgins &
Vink (2020) did not consider the same temperature range as us, there is no
discrepancy between the results. It is important to realise that the HD limit
is not limited only to RSGs, and hence all relevant temperatures need to be
included.

It is also evident from our results that the extent of overshooting affects
the formation of WR stars. Enhanced mixing reduces the mass threshold for a
star to remove its envelope through winds, and therefore affects the necessity
of binary interactions in forming WR stars (Shenar et al., 2020), as well as
decreasing the single-star evolution WR luminosity threshold. For example,
Figs. 2-3 show that increasing $\alpha_{\mathrm{ov}}$ from $0.335$ to $1.2$
reduces the lowest luminosity of a WR star formed through single-star
evolution from $\log_{10}\left(L/\mathrm{L}_{\odot}\right)\approx 6$ to
$\approx 5.5$ for the LMC, and from $\approx 5.5$ to $\approx 5.25$ for the
MW.

Our results also have implications for energetic transient astrophysical
phenomena. Firstly, the initial mass threshold for core-collapse supernovae is
lower for enhanced mixing. Secondly, the position in the HRD where stellar
models end their lives is substantially hotter for enhanced mixing. This has
implications for the so-called “Red Supergiant Problem”, which is the apparent
discrepancy between the most luminous progenitor of a Type IIP supernova and
the most luminous known RSGs (Smartt et al., 2009; Horiuchi et al., 2014;
Davies & Beasor, 2018, 2020; Kochanek, 2020). Our models with lower
overshooting values exhibit such a problem, though for enhanced mixing the
more massive stars do not die as RSGs, in accordance with observations.

The formation of black hole binaries (which are progenitors of gravitational
wave events) depends on the expansion of massive stars and their interaction
with a companion. Klencki et al. (2020) discuss the role of metallicity in
stellar expansion and interaction, as with lower metallicity stars are
generally more compact because of the lower gas opacity. Smaller stellar radii
imply that a smaller fraction of stars will experience significant binary
interactions, such as common envelope evolution, that lead to short-period
binary black holes which will merge quickly enough to produce observable
gravitational wave events. Our models with enhanced mixing expand to smaller
radii in general, therefore affecting the evolution towards merging black
holes. Since a large fraction of massive stars interact with a companion
during their evolution

Klencki et al. (2020) also discuss the role of MLT++ in relation to the HD
limit. According to Klencki et al. (2020), the more limited expansion of
massive stellar models employing MLT++ is favorable in terms of explaining the
HD limit, and it might be more accurate for rather massive ($M\ga
50\,\mathrm{M}_{\odot}$) stars. We note that even though we use the favorable
MLT++ prescription we still find an excess of over-luminous stars in the
temperature range relevant for the HD limit.

It is important to stress that our study focuses on the upper-mass end. For
example, applying high $\alpha_{\rm ov}$ values as invoked here can suppress
blue loops (Stothers & Chin, 1991; Walmswell, Tout & Eldridge, 2015) on the
HRD (Fig. 2) and prevent lower-mass stars from becoming Cepheids (Anderson et
al., 2014, 2016), in contrast to observations. However, studies of
intermediate-mass stars suggest a mass-dependent convective overshoot extent
(Claret & Torres, 2016, 2017, 2018, 2019). These studies suggest an increase
up to $\alpha_{\mathrm{ov}}\approx 0.2$ at $M\approx 2\,\mathrm{M}_{\odot}$,
plateauing afterwards. A similar trend is explained by Jermyn et al. (2018).
Overshooting is less well-constrained for the higher masses which are relevant
for our study and for the existence of the HD limit. Models of massive stars
make various calibrations for the overshooting extent, with results ranging
from $\alpha_{\mathrm{ov}}=0.1$ calibrated for the mass range of
$1.3$-$9\,\mathrm{M}_{\odot}$ (Ekström et al., 2012) to
$\alpha_{\mathrm{ov}}=0.5$ for higher masses around $M\ga
30\,\mathrm{M}_{\odot}$ (Higgins & Vink, 2019). So while very high overshoot
values cannot be applied to the lowest masses in our grids
($M=4\,\mathrm{M}_{\odot}$), higher mass stars might exhibit behaviour
appropriate for enhanced mixing, with few real constraints on the extent.

In conclusion, we propose that internal mixing in massive stars might play an
important part in explaining the empiric HD limit. Enhanced mixing prevents
the redward evolution of stellar models towards or beyond the HD limit. We do
not suggest to adopt a higher overshooting parameter, but rather that our
results hint at a deficiency in the modelling of mixing in stellar interiors
(e.g., Aerts et al., 2019; Schootemeijer et al., 2019). The extent of core
overshooting might be highly mass-dependent, or rotational mixing is more
efficient for slowly rotating stars. Moreover, it is well possible that the
final explanation relies on a multitude of mechanisms (e.g., mixing, mass
loss, multiplicity, star formation history; see Sect. 5). We therefore
encourage future investigations of this problem that address these various
mechanisms. For now, internal mixing in massive stars remains an unresolved
issue in stellar modelling with broad implications, and the origin of the
Humphreys-Davidson limit remains uncertain.

## Acknowledgments

We acknowledge the constructive report by our referee Dr. Cyril Georgy. This
research has made use of the SIMBAD database and cross-match service operated
at CDS, Strasbourg, France. The research has received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme with grant number 772225: MULTIPLES. AG and
TS would like to thank Danny Lennon and Ben Davies for helpful input and
exchanges. The Flatiron Institute is supported by the Simons Foundation. IA is
a CIFAR Azrieli Global Scholar in the Gravity and the Extreme Universe Program
and acknowledges support from that program, from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation program
(grant agreement number 852097), from the Israel Science Foundation (grant
numbers 2108/18 and 2752/19), from the United States - Israel Binational
Science Foundation (BSF), and from the Israeli Council for Higher Education
Alon Fellowship.

## Data Availability Statement

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Aerts et al. (2019) Aerts C., Mathis S., Rogers T. M., 2019, ARA&A, 57, 35
  * Alastuey & Jancovici (1978) Alastuey A., Jancovici B., 1978, ApJ, 226, 1034
  * Almeida et al. (2017) Almeida L. A., et al., 2017, A&A, 598, A84
  * Anderson et al. (2014) Anderson R. I., Ekström S., Georgy C., Meynet G., Mowlavi N., Eyer L., 2014, A&A, 564, A100
  * Anderson et al. (2016) Anderson R. I., Saio H., Ekström S., Georgy C., Meynet G., 2016, A&A, 591, A8
  * Ardeberg et al. (1972) Ardeberg A., Brunet J. P., Maurice E., Prevot L., 1972, A&AS, 6, 249
  * Beasor et al. (2020) Beasor E. R., Davies B., Smith N., van Loon J. T., Gehrz R. D., Figer D. F., 2020, MNRAS, 492, 5994
  * Bestenlehner et al. (2011) Bestenlehner J. M., et al., 2011, A&A, 530, L14
  * Björklund et al. (2020) Björklund R., Sundqvist J. O., Puls J., Najarro F., 2020, arXiv e-prints, p. arXiv:2008.06066
  * Bouret et al. (2013) Bouret J. C., Lanz T., Martins F., Marcolino W. L. F., Hillier D. J., Depagne E., Hubeny I., 2013, A&A, 555, A1
  * Brott et al. (2011) Brott I., et al., 2011, A&A, 530, A115
  * Buchler & Yueh (1976) Buchler J. R., Yueh W. R., 1976, ApJ, 210, 440
  * Cassisi et al. (2007) Cassisi S., Potekhin A. Y., Pietrinferni A., Catelan M., Salaris M., 2007, ApJ, 661, 1094
  * Castro et al. (2018) Castro N., Oey M. S., Fossati L., Langer N., 2018, ApJ, 868, 57
  * Chieffi & Limongi (2013) Chieffi A., Limongi M., 2013, ApJ, 764, 21
  * Claret & Torres (2016) Claret A., Torres G., 2016, A&A, 592, A15
  * Claret & Torres (2017) Claret A., Torres G., 2017, ApJ, 849, 18
  * Claret & Torres (2018) Claret A., Torres G., 2018, ApJ, 859, 100
  * Claret & Torres (2019) Claret A., Torres G., 2019, ApJ, 876, 134
  * Crowther et al. (2002) Crowther P. A., Dessart L., Hillier D. J., Abbott J. B., Fullerton A. W., 2002, A&A, 392, 653
  * Crowther et al. (2010) Crowther P. A., Schnurr O., Hirschi R., Yusof N., Parker R. J., Goodwin S. P., Kassim H. A., 2010, MNRAS, 408, 731
  * Cyburt et al. (2010) Cyburt R. H., et al., 2010, ApJS, 189, 240
  * Davies & Beasor (2018) Davies B., Beasor E. R., 2018, MNRAS, 474, 2116
  * Davies & Beasor (2020) Davies B., Beasor E. R., 2020, MNRAS, 493, 468
  * Davies et al. (2018) Davies B., Crowther P. A., Beasor E. R., 2018, MNRAS, 478, 3138
  * Dewitt et al. (1973) Dewitt H. E., Graboske H. C., Cooper M. S., 1973, ApJ, 181, 439
  * Drout et al. (2012) Drout M. R., Massey P., Meynet G., 2012, ApJ, 750, 97
  * Dubois et al. (1977) Dubois P., Jaschek M., Jaschek C., 1977, A&A, 60, 205
  * Dufton et al. (2000) Dufton P. L., McErlean N. D., Lennon D. J., Ryans R. S. I., 2000, A&A, 353, 311
  * Dufton et al. (2019) Dufton P. L., Evans C. J., Hunter I., Lennon D. J., Schneider F. R. N., 2019, A&A, 626, A50
  * Ekström et al. (2012) Ekström S., et al., 2012, A&A, 537, A146
  * Evans & Howarth (2003) Evans C. J., Howarth I. D., 2003, MNRAS, 345, 1223
  * Evans et al. (2006) Evans C. J., Lennon D. J., Smartt S. J., Trundle C., 2006, A&A, 456, 623
  * Feast et al. (1960) Feast M. W., Thackeray A. D., Wesselink A. J., 1960, MNRAS, 121, 337
  * Fehrenbach (1972) Fehrenbach C., 1972, A&A, 19, 427
  * Fehrenbach & Duflot (1982) Fehrenbach C., Duflot M., 1982, A&AS, 48, 409
  * Ferguson et al. (2005) Ferguson J. W., Alexander D. R., Allard F., Barman T., Bodnarik J. G., Hauschildt P. H., Heffner-Wong A., Tamanai A., 2005, ApJ, 623, 585
  * Fitzpatrick & Garmany (1990) Fitzpatrick E. L., Garmany C. D., 1990, ApJ, 363, 119
  * Flower (1996) Flower P. J., 1996, ApJ, 469, 355
  * Fuller et al. (1985) Fuller G. M., Fowler W. A., Newman M. J., 1985, ApJ, 293, 1
  * Fuller et al. (2019) Fuller J., Piro A. L., Jermyn A. S., 2019, MNRAS, 485, 3661
  * Fulmer et al. (2020) Fulmer L. M., Gallagher J. S., Hamann W.-R., Oskinova L. M., Ramachandran V., 2020, A&A, 633, A164
  * Gaia Collaboration (2018) Gaia Collaboration 2018, VizieR Online Data Catalog, p. I/345
  * Georgy et al. (2013) Georgy C., et al., 2013, A&A, 558, A103
  * Georgy et al. (2014) Georgy C., Saio H., Meynet G., 2014, MNRAS, 439, L6
  * Glatzel & Kiriakidis (1993) Glatzel W., Kiriakidis M., 1993, MNRAS, 263, 375
  * González-Fernández et al. (2015) González-Fernández C., Dorda R., Negueruela I., Marco A., 2015, A&A, 578, A3
  * Hainich et al. (2014) Hainich R., et al., 2014, A&A, 565, A27
  * Hainich et al. (2015) Hainich R., Pasemann D., Todt H., Shenar T., Sand er A., Hamann W. R., 2015, A&A, 581, A21
  * Hainich et al. (2019) Hainich R., Ramachandran V., Shenar T., Sand er A. A. C., Todt H., Gruner D., Oskinova L. M., Hamann W. R., 2019, A&A, 621, A85
  * Heger et al. (2000) Heger A., Langer N., Woosley S. E., 2000, ApJ, 528, 368
  * Henyey et al. (1965) Henyey L., Vardya M. S., Bodenheimer P., 1965, ApJ, 142, 841
  * Herwig (2000) Herwig F., 2000, A&A, 360, 952
  * Higgins & Vink (2019) Higgins E. R., Vink J. S., 2019, A&A, 622, A50
  * Higgins & Vink (2020) Higgins E. R., Vink J. S., 2020, A&A, 635, A175
  * Horiuchi et al. (2014) Horiuchi S., Nakamura K., Takiwaki T., Kotake K., Tanaka M., 2014, MNRAS, 445, L99
  * Humphreys (1983) Humphreys R. M., 1983, ApJ, 265, 176
  * Humphreys & Davidson (1979) Humphreys R. M., Davidson K., 1979, ApJ, 232, 409
  * Humphreys et al. (2016) Humphreys R. M., Weis K., Davidson K., Gordon M. S., 2016, ApJ, 825, 64
  * Hunter et al. (2007) Hunter I., et al., 2007, A&A, 466, 277
  * Hunter et al. (2008) Hunter I., Lennon D. J., Dufton P. L., Trundle C., Simón-Díaz S., Smartt S. J., Ryans R. S. I., Evans C. J., 2008, A&A, 479, 541
  * Iglesias & Rogers (1993) Iglesias C. A., Rogers F. J., 1993, ApJ, 412, 752
  * Iglesias & Rogers (1996) Iglesias C. A., Rogers F. J., 1996, ApJ, 464, 943
  * Indu & Subramaniam (2011) Indu G., Subramaniam A., 2011, A&A, 535, A115
  * Itoh et al. (1979) Itoh N., Totsuji H., Ichimaru S., Dewitt H. E., 1979, ApJ, 234, 1079
  * Itoh et al. (1996) Itoh N., Hayashi H., Nishikawa A., Kohyama Y., 1996, ApJS, 102, 411
  * Jermyn et al. (2018) Jermyn A. S., Tout C. A., Chitre S. M., 2018, MNRAS, 480, 5427
  * Kalari et al. (2018) Kalari V. M., Vink J. S., Dufton P. L., Fraser M., 2018, A&A, 618, A17
  * Kippenhahn et al. (1980) Kippenhahn R., Ruschenplatt G., Thomas H. C., 1980, A&A, 91, 175
  * Klencki et al. (2020) Klencki J., Nelemans G., Istrate A. G., Pols O., 2020, A&A, 638, A55
  * Kochanek (2020) Kochanek C. S., 2020, MNRAS, 493, 4945
  * Köhler et al. (2015) Köhler K., et al., 2015, A&A, 573, A71
  * Korn et al. (2000) Korn A. J., Becker S. R., Gummersbach C. A., Wolf B., 2000, A&A, 353, 655
  * Lamers & Fitzpatrick (1988) Lamers H. J. G. L. M., Fitzpatrick E. L., 1988, ApJ, 324, 279
  * Langanke & Martínez-Pinedo (2000) Langanke K., Martínez-Pinedo G., 2000, Nuclear Physics A, 673, 481
  * Langer et al. (1983) Langer N., Fricke K. J., Sugimoto D., 1983, A&A, 126, 207
  * Lennon (1997) Lennon D. J., 1997, A&A, 317, 871
  * Levesque et al. (2005) Levesque E. M., Massey P., Olsen K. A. G., Plez B., Josselin E., Maeder A., Meynet G., 2005, ApJ, 628, 973
  * Lindegren et al. (2016) Lindegren L., et al., 2016, A&A, 595, A4
  * Maeder & Meynet (1987) Maeder A., Meynet G., 1987, A&A, 182, 243
  * Maeder & Meynet (2000) Maeder A., Meynet G., 2000, A&A, 361, 159
  * Mahy et al. (2020) Mahy L., et al., 2020, A&A, 634, A118
  * Massey & Evans (2016) Massey P., Evans K. A., 2016, ApJ, 826, 224
  * Massey & Olsen (2003) Massey P., Olsen K. A. G., 2003, AJ, 126, 2867
  * Meynet & Maeder (1997) Meynet G., Maeder A., 1997, A&A, 321, 465
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Neugent et al. (2010) Neugent K. F., Massey P., Skiff B., Drout M. R., Meynet G., Olsen K. A. G., 2010, ApJ, 719, 1784
  * Neugent et al. (2012) Neugent K. F., Massey P., Skiff B., Meynet G., 2012, ApJ, 749, 177
  * Nugis & Lamers (2000) Nugis T., Lamers H. J. G. L. M., 2000, A&A, 360, 227
  * Oda et al. (1994) Oda T., Hino M., Muto K., Takahara M., Sato K., 1994, Atomic Data and Nuclear Data Tables, 56, 231
  * Ohnaka et al. (2008) Ohnaka K., Driebe T., Hofmann K. H., Weigelt G., Wittkowski M., 2008, A&A, 484, 371
  * 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
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton B., et al., 2018, ApJS, 234, 34
  * Pietrzyński et al. (2013) Pietrzyński G., et al., 2013, Nature, 495, 76
  * Potekhin & Chabrier (2010) Potekhin A. Y., Chabrier G., 2010, Contributions to Plasma Physics, 50, 82
  * Potter et al. (2012) Potter A. T., Tout C. A., Eldridge J. J., 2012, MNRAS, 419, 748
  * Ramachandran et al. (2018a) Ramachandran V., Hainich R., Hamann W. R., Oskinova L. M., Shenar T., Sander A. A. C., Todt H., Gallagher J. S., 2018a, A&A, 609, A7
  * Ramachandran et al. (2018b) Ramachandran V., Hamann W. R., Hainich R., Oskinova L. M., Shenar T., Sander A. A. C., Todt H., Gallagher J. S., 2018b, A&A, 615, A40
  * Ramachandran et al. (2019) Ramachandran V., et al., 2019, A&A, 625, A104
  * Ramírez-Agudelo et al. (2015) Ramírez-Agudelo O. H., et al., 2015, A&A, 580, A92
  * Ramírez-Agudelo et al. (2017) Ramírez-Agudelo O. H., et al., 2017, A&A, 600, A81
  * Rogers & Nayfonov (2002) Rogers F. J., Nayfonov A., 2002, ApJ, 576, 1064
  * Salpeter (1954) Salpeter E. E., 1954, Australian Journal of Physics, 7, 373
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Sanduleak (1970) Sanduleak N., 1970, Contributions from the Cerro Tololo Inter-American Observatory, 89
  * Sanyal et al. (2017) Sanyal D., Langer N., Szécsi D., -C Yoon S., Grassitelli L., 2017, A&A, 597, A71
  * Saumon et al. (1995) Saumon D., Chabrier G., van Horn H. M., 1995, ApJS, 99, 713
  * Schneider et al. (2018) Schneider F. R. N., et al., 2018, A&A, 618, A73
  * Schootemeijer & Langer (2018) Schootemeijer A., Langer N., 2018, A&A, 611, A75
  * Schootemeijer et al. (2019) Schootemeijer A., Langer N., Grin N. J., Wang C., 2019, A&A, 625, A132
  * Shaviv & Salpeter (1973) Shaviv G., Salpeter E. E., 1973, ApJ, 184, 191
  * Shenar et al. (2016) Shenar T., et al., 2016, A&A, 591, A22
  * Shenar et al. (2017) Shenar T., et al., 2017, A&A, 598, A85
  * Shenar et al. (2019) Shenar T., et al., 2019, A&A, 627, A151
  * Shenar et al. (2020) Shenar T., Gilkis A., Vink J. S., Sana H., Sander A. A. C., 2020, A&A, 634, A79
  * Smartt et al. (2009) Smartt S. J., Eldridge J. J., Crockett R. M., Maund J. R., 2009, MNRAS, 395, 1409
  * Smith (2019) Smith N., 2019, MNRAS, 489, 4378
  * Stock et al. (1976) Stock J., Osborn W., Ibanez M., 1976, A&AS, 24, 35
  * Stothers & Chin (1991) Stothers R. B., Chin C.-W., 1991, ApJ, 374, 288
  * Tabernero et al. (2018) Tabernero H. M., Dorda R., Negueruela I., González-Fernández C., 2018, MNRAS, 476, 3106
  * Tehrani et al. (2019) Tehrani K. A., Crowther P. A., Bestenlehner J. M., Littlefair S. P., Pollock A. M. T., Parker R. J., Schnurr O., 2019, MNRAS, 484, 2692
  * Timmes & Swesty (2000) Timmes F. X., Swesty F. D., 2000, ApJS, 126, 501
  * Torres (2010) Torres G., 2010, AJ, 140, 1158
  * Tramper et al. (2015) Tramper F., et al., 2015, A&A, 581, A110
  * Tramper et al. (2016) Tramper F., Sana H., de Koter A., 2016, ApJ, 833, 133
  * Trundle & Lennon (2005) Trundle C., Lennon D. J., 2005, A&A, 434, 677
  * Trundle et al. (2004) Trundle C., Lennon D. J., Puls J., Dufton P. L., 2004, A&A, 417, 217
  * Trundle et al. (2007) Trundle C., Dufton P. L., Hunter I., Evans C. J., Lennon D. J., Smartt S. J., Ryans R. S. I., 2007, A&A, 471, 625
  * Ulmer & Fitzpatrick (1998) Ulmer A., Fitzpatrick E. L., 1998, ApJ, 504, 200
  * Urbaneja et al. (2017) Urbaneja M. A., Kudritzki R. P., Gieren W., Pietrzyński G., Bresolin F., Przybilla N., 2017, AJ, 154, 102
  * Vink et al. (2001) Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, A&A, 369, 574
  * Vink et al. (2010) Vink J. S., Brott I., Gräfener G., Langer N., de Koter A., Lennon D. J., 2010, A&A, 512, L7
  * Walmswell et al. (2015) Walmswell J. J., Tout C. A., Eldridge J. J., 2015, MNRAS, 447, 2951
  * Wenger et al. (2000) Wenger M., et al., 2000, A&AS, 143, 9
  * Wolf (1972) Wolf B., 1972, A&A, 20, 275
  * Woosley (2019) Woosley S. E., 2019, ApJ, 878, 49
  * Yoon (2017) Yoon S.-C., 2017, MNRAS, 470, 3970
  * de Jager et al. (1988) de Jager C., Nieuwenhuijzen H., van der Hucht K. A., 1988, A&AS, 72, 259


direct numerical simulations shows generally a fairly good agreement,
suggesting that their influence is not significant in the cases that we have
analyzed. Nevertheless, some crucial aspects are worth to be discussed.
The second order response discussed in §5 includes several contributions,
among which the time- and azimuthal-averaged flows, i.e. steady and steady-
axisymmetric mean flows, associated with the two perfectly balanced counter-
rotating waves. Within the present WNL model, these steady forcing terms have
been considered non-resonant and the corresponding responses have been
therefore computed straightforwardly. However, from a mathematical
perspective, one must notice that existence of viscous modes, depending on the
associated damping coefficient, may induce large second order mean flow
responses when system (27) is inverted for the couple
$\left(m^{ij},\omega^{ij}\right)=\left(m,0\right)$ or $\left(0,0\right)$, i.e.
the system could be nearly singular if the damping coefficients of some of the
first viscous modes are sufficiently small and therefore close to neutrality.
It is fundamental to note that for
$\left(m^{ij},\omega^{ij}\right)=\left(0,0\right)$ the linearized governing
equations decouple into two sets of equations, one governing the vertical-
plane dynamics, $r$-$z$, and one governing the azimuthal velocity component
only. Therefore, the steady axisymmetric mean flow responses can be
interpreted as the sum of two contributions, i.e. an in-plane or toroidal mean
flow in the $r$-$z$ plane (analogous to that recently studied experimentally
by Périnet et al. (2017) with regard to a three-dimensional rectangular
container) and a pure azimuthal or drift mean flow component. When a standing
wave solution is considered, the second order toroidal mean flows produced by
the two counter rotating waves sum up, whereas the two opposite azimuthal mean
flows cancel out, leading to a net zero drift. However, through an amplitude
equation model capturing the infinite number of viscous modes at leading order
(together with the classic standing wave form), it has been shown (Fauve et
al., 1991; Knobloch et al., 2002; Martel et al., 2000; Vega et al., 2001) that
in some conditions the amplitude of one of the two azimuthal mean flow
component may exceed that of the opposite mean flow, hence inducing,
particularly at larger standing wave amplitudes, a secondary instability via
space-reflection symmetry breaking.
An exception is made by standing waves produced by the combination of two
axisymmetric, $m=0$, capillary–gravity waves, as that analyzed in §6 for $m=0$
and $n=2$. In this specific case, the resulting second order azimuthal mean
flow is zero by construction precluding the existence of such a secondary
instability and the overall mean flow acts toroidally (in the $r$-$z$ plane)
exclusively. In other words, although the influence of viscous modes should
not be neglected a priori, these arguments explain how standing wave dynamics
produced by axisymmetric capillary–gravity waves, as the one investigated
through DNS in §6 of the present manuscript, do not suffer from secondary
drift instability.

#### Supplementary Material

Wolfram Mathematica codes developed in this work for the automatized
linearization process and specifically for extraction of the WNL second order
forcing terms and third order resonating terms, discussed in §5, are available
to the readers as a supplementary material.

#### Acknowledgements

The authors wish to thank S. Gomé for his experimental observations that have
motivated this work as well as Dr. L. Siconolfi for his support.

#### Funding

We acknowledge the Swiss National Science Foundation under grant
200021_178971.

#### Declaration of Interests

The authors report no conflict of interest.

## References

  * Basak et al. (2021) Basak, S., Farsoiya, P. K. & Dasgupta, R. 2021 Jetting in finite-amplitude, free, capillary-gravity waves. J. Fluid Mech. 909.
  * Batson et al. (2013) Batson, W., Zoueshtiagh, F. & Narayanan, R. 2013 The faraday threshold in small cylinders and the sidewall non-ideality. J. Fluid Mech. 729 (496-523), 9.
  * Bechhoefer et al. (1995) Bechhoefer, J., Ego, V., Manneville, S. & Johnson, B. 1995 An experimental study of the onset of parametrically pumped surface waves in viscous fluids. J. Fluid Mech. 288, 325–350.
  * Benjamin & Scott (1979) Benjamin, T. B. & Scott, J. C. 1979 Gravity-capillary waves with edge constraints. J. Fluid Mech. 92, 241–267.
  * Benjamin & Ursell (1954) Benjamin, T. B. & Ursell, F. J. 1954 The stability of the plane free surface of a liquid in vertical periodic motion. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 225 (1163), 505–515.
  * Bongarzone et al. (2021) Bongarzone, A., Bertsch, A., Renaud, P. & Gallaire, F. 2021 Impinging planar jets: hysteretic behaviour and origin of the self-sustained oscillations. J. Fluid Mech. 913.
  * Bostwick & Steen (2009) Bostwick, J. B. & Steen, P. H. 2009 Capillary oscillations of a constrained liquid drop. Physics of Fluids 21 (3), 032108.
  * Canuto et al. (2007) Canuto, C., Hussaini, M. Y., Quarteroni, A. & Zang, T. A. 2007 Spectral methods: evolution to complex geometries and applications to fluid dynamics. Springer Science & Business Media.
  * Case & Parkinson (1957) Case, K. M. & Parkinson, W. C. 1957 Damping of surface waves in an incompressible liquid. Journal of Fluid Mechanics 2 (2), 172–184.
  * Chen & Vinals (1999) Chen, P. & Vinals, J. 1999 Amplitude equation and pattern selection in faraday waves. Physical Review E 60 (1), 559\.
  * Ciliberto & Gollub (1985) Ciliberto, S. & Gollub, J. P. 1985 Chaotic mode competition in parametrically forced surface waves. J. Fluid Mech. 158, 381–398.
  * Cocciaro et al. (1993) Cocciaro, B., Faetti, S. & Festa, C. 1993 Experimental investigation of capillarity effects on surface gravity waves: non-wetting boundary conditions. Journal of Fluid Mechanics 246, 43–66.
  * Cocciaro et al. (1991) Cocciaro, B., Faetti, S. & Nobili, M. 1991 Capillarity effects on surface gravity waves in a cylindrical container: wetting boundary conditions. Journal of Fluid Mechanics 231, 325–343.
  * Das & Hopfinger (2008) Das, S. P. & Hopfinger, E. J. 2008 Parametrically forced gravity waves in a circular cylinder and finite-time singularity. J. Fluid Mech. 599, 205.
  * Davis (1974) Davis, S. H. 1974 On the motion of a fluid-fluid interface along a solid surface. J. Fluid Mech. 65 (1), 71–95.
  * Dodge et al. (1965) Dodge, F. T., Kana, D. D. & Abramson, H. N. 1965 Liquid surface oscillations in longitudinally excited rigid cylindrical containers. Aiaa Journal 3 (4), 685–695.
  * Douady (1990) Douady, S. 1990 Experimental study of the Faraday instability. Journal of fluid mechanics 221, 383–409.
  * Edwards & Fauve (1994) Edwards, W. S. & Fauve, S. 1994 Patterns and quasi-patterns in the faraday experiment. J. Fluid Mech 278, 123–148.
  * Faraday (1831) Faraday, M. 1831 Phil. Trans. .
  * Fauve et al. (1991) Fauve, S., Douady, S. & Thual, O. 1991 Drift instabilities of cellular patterns. Journal de Physique II 1 (3), 311–322.
  * Friedrichs (2012) Friedrichs, K. O. 2012 Spectral theory of operators in Hilbert space. Springer Science & Business Media.
  * Gluckman et al. (1993) Gluckman, B. J., Marcq, P., Bridger, J. & Gollub, J. P. 1993 Time averaging of chaotic spatiotemporal wave patterns. Physical review letters 71 (13), 2034.
  * Graham-Eagle (1983) Graham-Eagle, J. 1983 A new method for calculating eigenvalues with applications to gravity-capillary waves with edge constraints. In Mathematical Proceedings of the Cambridge Philosophical Society, , vol. 94, pp. 553–564. Cambridge University Press.
  * Gu & Sethna (1987) Gu, X. M. & Sethna, P. R. 1987 Resonant surface waves and chaotic phenomena. J. Fluid Mech. 183, 543–565.
  * Heinrichs (2004) Heinrichs, W. 2004 Spectral collocation schemes on the unit disc. J. Comp. Physics 199 (1), 66–86.
  * Henderson et al. (1992) Henderson, D., Hammack, J., Kumar, P. & Shah, D. 1992 The effects of static contact angles on standing waves. Physics of Fluids A: Fluid Dynamics 4 (10), 2320–2322.
  * Henderson & Miles (1990) Henderson, D. M. & Miles, J. W. 1990 Single-mode faraday waves in small cylinders. Journal of Fluid Mechanics 213, 95–109.
  * Henderson & Miles (1994) Henderson, D. M. & Miles, J. W. 1994 Surface-wave damping in a circular cylinder with a fixed contact line. J. Fluid Mech. 275, 285–299.
  * Hocking (1987) Hocking, L. M. 1987 The damping of capillary–gravity waves at a rigid boundary. Journal of fluid mechanics 179, 253–266.
  * Howell et al. (2000) Howell, D. R., Buhrow, B., Heath, T., McKenna, C., Hwang, W. & Schatz, M. F. 2000 Measurements of surface-wave damping in a container. Physics of Fluids 12 (2), 322–326.
  * Hsu (1977) Hsu, C. S. 1977 On nonlinear parametric excitation problems. In Advances in Applied Mechanics, , vol. 17, pp. 245–301. Elsevier.
  * Huh & Scriven (1971) Huh, C. & Scriven, L. E. 1971 Hydrodynamic model of steady movement of a solid/liquid/fluid contact line. Journal of colloid and interface science 35 (1), 85–101.
  * Ito et al. (1999) Ito, T., Tsuji, Y. & Kukita, Y. 1999 Interface waves excited by vertical vibration of stratified fluids in a circular cylinder. Journal of nuclear science and technology 36 (6), 508–521.
  * Jian & Xuequan (2005) Jian, Y. & Xuequan, E. 2005 Instability analysis of nonlinear surface waves in a circular cylindrical container subjected to a vertical excitation. European Journal of Mechanics-B/Fluids 24 (6), 683–702.
  * Jiang et al. (2004) Jiang, L., Perlin, M. & Schultz, W. W. 2004 Contact-line dynamics and damping for oscillating free surface flows. Physics of Fluids 16 (3), 748–758.
  * Keulegan (1959) Keulegan, G. H. 1959 Energy dissipation in standing waves in rectangular basins. Journal of Fluid Mechanics 6 (1), 33–50.
  * Kidambi (2009) Kidambi, R. 2009 Meniscus effects on the frequency and damping of capillary-gravity waves in a brimful circular cylinder. Wave Motion 46 (2), 144–154.
  * Kidambi (2013) Kidambi, R. 2013 Inviscid Faraday waves in a brimful circular cylinder. J. Fluid Mech. 724, 671–694.
  * Knobloch et al. (2002) Knobloch, E., Martel, C. & Vega, J. M. 2002 Coupled mean flow-amplitude equations for nearly inviscid parametrically driven surface waves. Annals of the New York Academy of Sciences 974 (1), 201–219.
  * Kovacic et al. (2018) Kovacic, I., Rand, R. & Sah, S. M. 2018 Mathieu’s Equation and Its Generalizations: Overview of Stability Charts and Their Features. Applied Mechanics Reviews 70 (2), 020802.
  * Kumar & Tuckerman (1994) Kumar, K. & Tuckerman, L. S. 1994 Parametric instability of the interface between two fluids. Journal of Fluid Mechanics 279, 49–68.
  * Lam & Caps (2011) Lam, K. D. N. T. & Caps, H. 2011 Effect of a capillary meniscus on the faraday instability threshold. The European Physical Journal E 34 (10), 1–5.
  * Lamb (1932) Lamb, H. 1932 Hydrodynamics, 6th edn. Cambridge University Press.
  * Landau & Lifshitz (1959) Landau, L. D. & Lifshitz, M. 1959 Fluid Mechanics, 1st edn. Pergamon Press.
  * Lauga et al. (2007) Lauga, E., Brenner, M. & Stone, H. 2007 Microfluidics: the no–slip boundary condition. Springer handbook of experimental fluid mechanics pp. 1219–1240.
  * Liu & Liu (2012) Liu, R. & Liu, Q. S. 2012 Nonmodal stability in hagen-poiseuille flow of a shear thinning fluid. Physical Review E 85 (6), 066318.
  * Martel & Knobloch (1997) Martel, C. & Knobloch, E. 1997 Damping of nearly inviscid water waves. Physical Review E 56 (5), 5544.
  * Martel et al. (2000) Martel, C., Knobloch, E. & Vega, J. M. 2000 Dynamics of counterpropagating waves in parametrically forced systems. Physica D: Nonlinear Phenomena 137 (1-2), 94–123.
  * Martel et al. (1998) Martel, C., Nicolas, J. A. & Vega, J. M. 1998 Surface-wave damping in a brimful circular cylinder. J. Fluid Mech. 360, 213–228.
  * Matthiessen (1868) Matthiessen, L. 1868 Ann. Phys., Lpz. 134, 107\.
  * Matthiessen (1870) Matthiessen, L. 1870 Ann. Phys., Lpz. 141, 375\.
  * Meliga et al. (2009) Meliga, P., Chomaz, J.-M. & Sipp, D. 2009 Global mode interaction and pattern selection in the wake of a disk: a weakly nonlinear expansion. Journal of Fluid Mechanics 633, 159–189.
  * Meliga et al. (2012) Meliga, P., Gallaire, F. & Chomaz, J. M. 2012 A weakly nonlinear mechanism for mode selection in swirling jets. J. Fluid Mech. 699, 216–262.
  * Meron (1987) Meron, E. 1987 Parametric excitation of multimode dissipative systems. Physical Review A 35, 11.
  * Meron & Procaccia (1986) Meron, E. & Procaccia, I. 1986 Low-dimensional chaos in surface waves: Theoretical analysis of an experiment. Physical Review A 34 (4), 3221.
  * Miles (1990) Miles, J. 1990 Capillary-viscous forcing of surface waves. J. Fluid Mech. 219, 635–646.
  * Miles (1991) Miles, J. 1991 The capillary boundary layer for standing waves. J. Fluid Mech. 222, 197–205.
  * Miles (1967) Miles, J. W. 1967 Surface-wave damping in closed basins. Proc. R. Soc. Lond 297, 459–475.
  * Miles (1984) Miles, J. W. 1984 Nonlinear Faraday resonance. Journal of Fluid Mechanics 146, 285–302.
  * Miles & Henderson (1990) Miles, J. W. & Henderson, D. 1990 Parametrically forced surface waves. Annual Review of Fluid Mechanics 22 (1), 143–165.
  * Miles & Henderson (1998) Miles, J. W. & Henderson, D. M. 1998 A note on interior vs. boundary-layer damping of surface waves in a circular cylinder. J. Fluid Mech. 364, 319–323.
  * Milner (1991) Milner, S. T. 1991 Square patterns and secondary instabilities in driven capillary waves. Journal of fluid mechanics 225, 81–100.
  * Müller et al. (1997) Müller, H. W., Wittmer, H., Wagner, C., Albers, J. & Knorr, K. 1997 Analytic stability theory for faraday waves and the observation of the harmonic surface response. Physic. Rev. Lett. 78 (12), 2357.
  * Nagata (1989) Nagata, M. 1989 Nonlinear faraday resonance in a box with a square base. J. Fluid Mech. 209, 265–284.
  * Navier (1823) Navier, C. L. M. H. 1823 Mémoire sur les lois du mouvement des fluides. Mém. Acad. R. des Sci. Inst. France 6 (1823), 389–440.
  * Nayfeh (1987) Nayfeh, A. H. 1987 Surface waves in closed basins under parametric and internal resonances. The Physics of Fluids 30 (10), 2976–2983.
  * Nayfeh (2008) Nayfeh, A. H. 2008 Perturbation Methods. Wiley.
  * Nayfeh & Mook (1995) Nayfeh, A. H. & Mook, D. T. 1995 Nonlinear oscillations. John Wiley & Sons.
  * Nicolás (2002) Nicolás, J. A. 2002 The viscous damping of capillary-gravity waves in a brimful circular cylinder. Physics of Fluids 14 (6), 1910–1919.
  * Nicolás (2005) Nicolás, J. A. 2005 Effects of static contact angles on inviscid gravity-capillary waves. Physics of Fluids 17 (2), 022101\.
  * Périnet et al. (2017) Périnet, N., Gutiérrez, P., Urra, H., Mujica, N. & Gordillo, L. 2017 Streaming patterns in faraday waves. J. Fluid Mech. 819, 285–310.
  * Perlin & Schultz (2000) Perlin, M. & Schultz, W. W. 2000 Capillary effects on surface waves. Ann. Rev. Fluid Mech. 32 (1), 241–274.
  * Picard & Davoust (2007) Picard, C. & Davoust, L. 2007 Resonance frequencies of meniscus waves as a physical mechanism for a dna biosensor. Langmuir 23 (3), 1394–1402.
  * Rajchenbach & Clamond (2015) Rajchenbach, J. & Clamond, D. 2015 Faraday waves: Their dispersion relation, nature of bifurcation and wavenumber selection revisited. Journal of Fluid Mechanics 777.
  * Rayleigh (1883a) Rayleigh, L. 1883a On maintained vibrations. Phil. Mag. 15, 229-235 (Scientific Papers) 2, 188–193.
  * Rayleigh (1883b) Rayleigh, L. 1883b On the crispations of fluid resting on a vibrating support. Phil. Mag. 16, 50-58 (Scientific Papers) 2, 212–219.
  * Shao et al. (2021a) Shao, X., Gabbard, C. T., Bostwick, J. B. & Saylor, J. R. 2021a On the role of meniscus geometry in capillary wave generation. Experiments in Fluids 62 (3), 1–4.
  * Shao et al. (2021b) Shao, X., Wilson, P., Saylor, J. R. & Bostwick, J. B. 2021b Surface wave pattern formation in a cylindrical container. J. Fluid Mech. 915.
  * Skeldon & Guidoboni (2007) Skeldon, A. C. & Guidoboni, G. 2007 Pattern selection for faraday waves in an incompressible viscous fluid. SIAM Journal on Applied Mathematics 67 (4), 1064–1100.
  * Sommariva (2013) Sommariva, A. 2013 Fast construction of fejér and clenshaw–curtis rules for general weight functions. Computers & Mathematics with Applications 65 (4), 682–693.
  * Stuart & Fauve (1993) Stuart, W. S. Edwards & Fauve, S. 1993 Parametrically excited quasicrystalline surface waves. Physical Review E 47 (2), R788.
  * Ting & Perlin (1995) Ting, C. L. & Perlin, M. 1995 Boundary conditions in the vicinity of the contact line at a vertically oscillating upright plate: an experimental investigation. J. Fluid Mech. 295, 263–300.
  * Tipton (2003) Tipton, C. R. 2003 Interfacial Faraday waves in a small cylindrical cell. University of Manchester.
  * Tipton & Mullin (2004) Tipton, C. R. & Mullin, T. 2004 An experimental study of faraday waves formed on the interface between two immiscible liquids. Physics of Fluids 16 (7), 2336–2341.
  * Torres et al. (1995) Torres, M., Pastor, G., Jiménez, I. & Espinosa, F. M. De 1995 Five-fold quasicrystal-like germinal pattern in the faraday wave experiment. Chaos, Solitons & Fractals 5 (11), 2089–2093.
  * Vega et al. (2001) Vega, J. M., Knobloch, E. & Martel, C. 2001 Nearly inviscid faraday waves in annular containers of moderately large aspect ratio. Physica D: Nonlinear Phenomena 154 (3-4), 313–336.
  * Viola et al. (2016) Viola, F., Arratia, C. & Gallaire, F. 2016 Mode selection in trailing vortices: harmonic response of the non-parallel Batchelor vortex. Journal of Fluid Mechanics 790, 523–552.
  * Viola et al. (2018) Viola, F., Brun, P.-T. & Gallaire, F. 2018 Capillary hysteresis in sloshing dynamics: a weakly nonlinear analysis. Journal of Fluid Mechanics 837, 788–818.
  * Viola & Gallaire (2018) Viola, F. & Gallaire, F. 2018 Theoretical framework to analyze the combined effect of surface tension and viscosity on the damping rate of sloshing waves. Physical Review Fluids 3 (9), 094801\.
  * Virnig et al. (1988) Virnig, J. C., Berman, A. S. & Sethna, P. R. 1988 On three-dimensional nonlinear subharmonic resonant surface waves in a fluid. part 2. experiment. Trans. ASME E: J. Appl. Mech. 55, 220–224.
  * Ward et al. (2019) Ward, K., Zoueshtiagh, F. & Narayanan, R. 2019 Faraday instability in double-interface fluid layers. Physical Review Fluids 4 (4), 043903.
  * Zhang & Vinals (1997) Zhang, W. & Vinals, J. 1997 Pattern formation in weakly damped parametric surface waves driven by two frequency components. Journal of Fluid Mechanics 341, 225–244.