# Hydrodynamic simulations of the KT Eridani nova super-remnant

M. W. Healy-Kalesh,1 M. J. Darnley,1 M. M. Shara,2 K. M. Lanzetta,3 J. T.
Garland,2 and S. Gromoll4
1Astrophysics Research Institute, Liverpool John Moores University, IC2
Liverpool Science Park, Liverpool, L3 5RF, UK
2Department of Astrophysics, American Museum of Natural History, Central Park
West at 79th Street, New York, NY 10024, USA
3Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY
11794-3800, USA
4Amazon Web Services, 410 Terry Ave. N, Seattle, WA 98109, USA
E-mail<EMAIL_ADDRESS>(MWH-K)

(Accepted 2023 October 10. Received 2023 October 09; in original form 2023
August 19)

###### Abstract

A nova super-remnant (NSR) is an immense structure associated with a nova that
forms when frequent and recurrent nova eruptions sweep up surrounding
interstellar material (ISM) into a high density and distant shell. The
prototypical NSR, measuring over 100 pc across, was discovered in 2014 around
the annually erupting nova M 31N 2008-12a. Hydrodynamical simulations
demonstrated that the creation of a dynamic NSR by repeated eruptions
transporting large quantities of ISM is not only feasible but that these
structures should exist around all novae, whether the white dwarf (WD) is
increasing or decreasing in mass. But it is only the recurrent nova (RNe) with
the highest WD masses and accretion rates that should host observable NSRs. KT
Eridani is, potentially, the eleventh RNe recorded in the Galaxy and is also
surrounded by a recently unveiled H$\alpha$ shell tens of parsecs across,
consistent with a NSR. Through modelling the nova ejecta from KT Eri, we
demonstrate that such an observable NSR could form in approximately 50,000
years, which fits with the proper motion history of the nova. We compute the
expected H$\alpha$ emission from the KT Eri NSR and predict that the structure
might be accessible to wide-field X-ray facilities.

###### keywords:

stars: individual (KT Eri) – novae, cataclysmic variables – hydrodynamics –
ISM: general

††pubyear: 2023††pagerange: Hydrodynamic simulations of the KT Eridani nova
super-remnant–References

##  1 Introduction

Classical novae (CNe) are a type of cataclysmic variable that can suddenly
increase in luminosity by over five orders of magnitude before fading to
quiescence (Shara et al., 2012a; Shara et al., 2012b). CNe arise from
accreting binaries (Walker, 1954) in which hydrogen-rich material streams from
a less-evolved companion star (typically subdwarf, subgiant, or red giant;
Darnley et al., 2012) onto a white dwarf (WD) via an accretion disk (see,
e.g., Warner, 1995). As more matter accumulates, the layer at the base of the
accreted envelope is subjected to immense temperatures and compression as a
result of the WD’s intense surface gravity and consequently becomes electron
degenerate. This leads to an irreversible cascade of thermonuclear reactions
in the form of a thermonuclear runaway (TNR) that rips through the accreted
envelope (Starrfield et al., 1972). The TNR pushes the envelope towards the
Fermi temperature whereby electron degeneracy is lifted and the material
expands in response to the high temperatures; this is the nova eruption
(Starrfield et al., 1976). Around 10${}^{-4}\ \text{M}_{\odot}$ of material,
travelling between a few hundred and a few thousand kilometers per second, is
ejected during the eruption (O’Brien et al., 2001) forming an expanding nova
shell (see, for example, Wade, 1990; Slavin et al., 1995; Woudt & Ribeiro,
2014; Takeda et al., 2022).

The configuration of the nova binary remains unaltered following an eruption.
As such, accretion will resume after the re-establishment of the disk (Worters
et al., 2007) and the whole process will repeat. Recurrent novae (RNe) are the
resultant phenomenon. Recent theoretical models (Yaron et al., 2005; Kato et
al., 2014) show that all nova systems can accommodate these repeated eruptions
and therefore all CNe are inherently recurrent; they simply have inter-
eruptions periods that are far longer than the ${\sim}$100 years of modern
astronomical data. However, observationally, RNe are defined as exhibiting
more than one eruption from the same system; ten such RNe are located within
the Galaxy (Schaefer, 2010; Darnley, 2021), with the other twenty-three
residing in M 31 and the LMC (Darnley & Henze, 2020). As with CNe, nova shells
exist around RNe; a key difference being collisions between consecutive
ejecta, which can produce shock-heating and clumping of material as evidenced
around the Galactic RN, T Pyxidis (Shara et al., 1997; Toraskar et al., 2013).

M 31N 2008-12a (or 12a) is the most extreme RN known (see, e.g., Darnley et
al., 2016; Henze et al., 2018; Darnley & Henze, 2020; Darnley, 2021, and
references therein). Situated in the Andromeda galaxy, 12a has experienced a
nova eruption annually since at least 2008 (Healy-Kalesh et al. in prep)
equating to a mean recurrence period of 0.99 years (Darnley & Henze, 2020;
Darnley, 2021). These rapid-fire eruptions are produced as a result of the
most massive WD in any known nova system ($1.38\ \text{M}_{\odot}$; Kato et
al., 2015) accreting material from its companion at an extreme rate
(${\sim}10^{-6}\ \text{M}_{\odot}\ \text{yr}^{-1}$; Darnley et al., 2017). The
exceptional nature of 12a was further showcased through the discovery of a
vast surrounding elliptical structure (Darnley et al., 2015). With a projected
size of $90\times 134$ parsecs (Darnley et al., 2019) far larger than any
single-eruption nova shell (e.g., Bode et al., 2004; Shara et al., 2012a;
Shara et al., 2012b) as well as the majority of supernova remnants (Stil &
Irwin, 2001), the first nova super-remnant (NSR) had been uncovered (Darnley
et al., 2019).

Darnley et al. (2019) demonstrated, through hydrodynamical modelling
(utilising the MORPHEUS code; Vaytet, O’Brien & Bode, 2007), that an extensive
NSR could be created through the continual sweeping up of local interstellar
medium by the rapid-fire eruptions of 12a over the system’s lifetime.
Extensive simulations presented in Healy-Kalesh et al. (2023) exploring the
influence of several system parameters (including ISM density, accretion rate,
WD temperature, and initial WD mass) on NSR growth found that all RNe should
be surrounded by vast dynamic NSRs. However, 12a had remained the only nova
known to host a NSR.

KT Eridani (Nova Eridani 2009) was a bright ($V\sim 5.42$), very fast-fading
($t_{3}\sim 13.6$ d) nova (Schaefer et al., 2022) discovered on November 25
2009 (Yamaoka et al., 2009; Itagaki, 2009). Comprehensive photometric coverage
of the pre- and post-maximum light curve with the Solar Mass Ejection Imager
(SMEI; Jackson et al., 1997) instrument on the Coriolos satellite captured the
elusive pre-maximum halt phase of novae evolution (Hounsell et al., 2010a, b).
Though observationally classified as a CN, KT Eri is often regarded as the
eleventh RN in the Galaxy on account of the system’s characteristics, rather
than the detection of a second eruption (Pagnotta & Schaefer, 2014; Schaefer
et al., 2022). Specifically it is very fast-fading; has a small outburst
amplitude; exhibits spectra containing triple-peaked H$\alpha$ emission with
FWHM of 3200–3400 km s-1 and high excitation lines (He ii); has two light
curve plateaus and has a WD mass of $1.25\pm 0.03\text{M}_{\odot}$
(>1.25$\text{M}_{\odot}$) orbiting a subgiant companion
(T${}_{\text{eff}}=6200$ K) with an orbital period of $P=2.61595$ d (Pagnotta
& Schaefer, 2014; Schaefer et al., 2022). Furthermore, Schaefer et al. (2022)
use their determinations of the WD mass and accretion rate ($3.5\times
10^{-7}\ \text{M}_{\odot}\ \text{yr}^{-1}$) to derive a recurrence period of
40–50 years.

Narrowband imaging of the surroundings of KT Eri, obtained with the Condor
Array Telescope (Lanzetta et al., 2023), reveal a vast structure, ${\sim}50$
parsecs in diameter, coincident with KT Eri (Shara et al. 2023), akin to 12a’s
NSR (Darnley et al., 2019).

In this Paper, we model nova ejecta from a fixed mass WD in an attempt to
replicate the NSR observed surrounding KT Eri. In Section 2 we present our
hydrodynamic simulations before presenting the results in Section 3. We then
discuss consistencies between our modelling and the observed NSR in Section 4
before concluding our paper in Section 5.

Table 1: Parameters for each simulation. Columns record the simulation number, recurrence period, ejecta velocity, number of eruptions to grow the nova super-remnant to 25 parsecs (the angular radius of the observed NSR around KT Eri - see Shara et al. 2023), the elapsed (cumulative) time of the simulation, and the total kinetic energy released. Runs 1⋆ – 5⋆ have the same characteristics as Runs 1 – 5 but the single component ejecta has been replaced by multicomponent ejecta as described in Section 2.4. Run 4† has the same characteristics as Run 4 however the recurrence period differs from that given as set by a Gaussian distribution (see Section 3.4 for details). Run # | $\text{P}_{\text{rec}}$ | $v_{\text{ej}}$ | Number of | Elapsed time | Total KE
---|---|---|---|---|---
| (yrs) | (km/s) | eruptions | (years) | (erg)
01 |  005 | 6,000 | 10,113 | $5.1\times 10^{4}$ | $6.5\times 10^{48}$
02 |  010 | 6,000 |  05,076 | $5.1\times 10^{4}$ | $6.5\times 10^{48}$
03 |  020 | 6,000 |  02,543 | $5.1\times 10^{4}$ | $6.5\times 10^{48}$
04 |  050 | 6,000 |  01,019 | $5.1\times 10^{4}$ | $6.5\times 10^{48}$
05 | 100 | 6,000 |  0,0510 | $5.1\times 10^{4}$ | $6.5\times 10^{48}$
06 |  005 | 5,000 | 11,435 | $5.8\times 10^{4}$ | $5.1\times 10^{48}$
07 |  010 | 5,000 |  05,739 | $5.8\times 10^{4}$ | $5.1\times 10^{48}$
08 |  020 | 5,000 |  02,875 | $5.8\times 10^{4}$ | $5.1\times 10^{48}$
09 |  050 | 5,000 |  01,152 | $5.8\times 10^{4}$ | $5.1\times 10^{48}$
10 | 100 | 5,000 |  0,0576 | $5.8\times 10^{4}$ | $5.1\times 10^{48}$
11 |  005 | 4,000 | 13,262 | $6.7\times 10^{4}$ | $3.8\times 10^{48}$
12 |  010 | 4,000 |  06,656 | $6.7\times 10^{4}$ | $3.8\times 10^{48}$
13 |  020 | 4,000 |  03,334 | $6.7\times 10^{4}$ | $3.8\times 10^{48}$
14 |  050 | 4,000 |  01,336 | $6.7\times 10^{4}$ | $3.8\times 10^{48}$
15 | 100 | 4,000 |  0,0668 | $6.7\times 10^{4}$ | $3.8\times 10^{48}$
01⋆ |  005 | variable | 10,113 | $5.1\times 10^{4}$ | $6.5\times 10^{48}$
02⋆ |  010 | variable |  05,076 | $5.1\times 10^{4}$ | $6.5\times 10^{48}$
03⋆ |  020 | variable |  02,543 | $5.1\times 10^{4}$ | $6.5\times 10^{48}$
04⋆ |  050 | variable |  01,019 | $5.1\times 10^{4}$ | $6.5\times 10^{48}$
05⋆ | 100 | variable |  0,0510 | $5.1\times 10^{4}$ | $6.5\times 10^{48}$
04† |  050 | 6,000 |  01,033 | $5.2\times 10^{4}$ | $6.6\times 10^{48}$

##  2 Simulations

As with previous work simulating NSRs (Darnley et al., 2019; Healy-Kalesh et
al., 2023), we utilised MORPHEUS (Vaytet et al., 2007) to model the evolution
of the NSR shell associated with KT Eri. MORPHEUS integrates three codes,
namely the one-dimensional Asphere (Vaytet et al., 2007), two-dimensional
Novarot (Lloyd et al., 1997) and the three-dimensional CubeMPI (Wareing et
al., 2006), developed by the Manchester University–LJMU Nova groups to
construct an MPI-OpenMP Eulerian second-order Godunov simulation code.

In Healy-Kalesh et al. (2023), the authors followed the growth of NSRs created
by nova ejecta from evolving (mass growth or reduction) WDs. Those simulations
were concerned with the extreme evolution of NSRs as the WD approached the
Chandrasekhar mass. For KT Eri, the much longer recurrence periods are close
to those during early ‘steady state’ NSR evolution. Therefore, we evolve a NSR
from a fixed mass WD (mass accumulation efficiency = 0) in a similar manner to
the simulations of the 12a NSR by Darnley et al. (2019).

We adopted the same mass loss rate of $2\times 10^{-8}\ \text{M}_{\odot}\
\text{yr}^{-1}$ and wind velocity of 20 km s-1 from the companion as used in
the previous studies (see Darnley et al., 2019; Healy-Kalesh et al., 2023, for
details). While KT Eri has a companion with an effective temperature that is
consistent with a subgiant (Schaefer et al., 2022), we use this mass loss rate
from the companion between nova eruptions predominantly for its help with
computational time. The eruptions from the central nova are then represented
by a instantaneous increase in ejecta velocity and mass-loss rate tuned to
match the eruption of KT Eri (as described in Section 1), punctuating the wind
mass-loss at intervals given by a set recurrence period. As in previous
studies, we assume one-dimensional spherical symmetry for these simulations
(largely for computational reasons). For our simulations, we have chosen a
resolution of 1000 au/cell as this is able to showcase the gross structure of
the NSR during its full evolution (see Section 3.1.1 and Figure 2 for a
comparison between different spatial resolutions), while allowing for feasible
computational times. Furthermore, the resolution we have chosen is so much
greater than the orbital separation of the WD and the companion in the KT Eri
system (10.4 R⊙; Schaefer et al., 2022) that we can ignore any ejecta-
donor/accretion disk interaction.

###  2.1 Radiative cooling

The role of radiative cooling in the formation of the 12a NSR was first
explored in Darnley et al. (2019). They found that cooling had little
influence on NSR evolution as the annual eruptions were highly energetic
throughout the full evolution and, did not allow sufficient time for the NSR
to cool. In contrast, Healy-Kalesh et al. (2023) found that radiative cooling
has a large influence on formation of NSRs (including its radial size and
shell structure) created by eruptions from an evolving WD; this is a result of
the long inter-eruption periods and low energy eruptions during the early
‘steady state’ growth permitting sufficient time for effective cooling. As
such, we incorporate the same radiative cooling in this work as employed in
previous work (Vaytet et al., 2007; Darnley et al., 2019; Healy-Kalesh et al.,
2023) whereby the Raymond et al. (1976) cooling curve within Morpheus is
adopted.

###  2.2 Local ISM density

The density of the surrounding ISM plays a pivotal role in shaping the growing
NSR (Healy-Kalesh et al., 2023). As such, we need to estimate the ISM density
of the immediate surroundings of KT Eri. KT Eri is situated at a Gaia distance
of $5110^{+920}_{-430}$ pc (Schaefer et al., 2022) and lies ${\sim}3$ kpc
below the Galactic plane (Shara et al. 2023). The scale height of gas in the
peripheries of Milky Way-like galaxies is ${\sim}800$ pc (Gensior et al.,
2023) therefore KT Eri lies approximately $4\times$ this scale height from the
Galactic plane. This indicates that the gas density in the region where KT Eri
is located is ${\sim}2$ per cent of the gas in the plane of the Galaxy in the
neighbourhood of KT Eri. Direct measurements by the Voyager 1 spacecraft
(Kurth et al., 2023) find the local ISM electron density (and therefore likely
hydrogen) in the Galactic plane to be ${\sim}0.1\ \text{cm}^{-3}$. Therefore,
we predict the ISM density around KT Eri to be $2\times 10^{-3}$ H atom per
cubic centimetre. This ISM density is utilised for all models in this study.
For consistency with models in Healy-Kalesh et al. (2023), we refer to this
ISM density by the number density $n=2\times 10^{-3}\ \text{cm}^{-3}$
(dropping the units) throughout.

###  2.3 Bulk ejecta simulations

The first set of simulations (Runs 1 – 15) follow the growth of a NSR from the
central nova ejecting a single bulk of material with every eruption. While a
simplistic interpretation of a nova eruption, Healy-Kalesh et al. (2023)
showed that the gross structure of a NSR is unaffected by the structure of
each ejecta (see Section 2.4 for the possible impact of a complex-multi-
component ejecta on NSR substructure).

We assume a static WD for our study as (i) suitable eruption models for KT Eri
do not yet exist, (ii) it is relatively poorly studied (compared to other RNe)
and (iii) its long recurrence period suggests KT Eri may be a relatively young
system compared to other RNe (e.g., RS Ophiuchi, U Scorpii and M 31N 2008-12a;
see e.g., Darnley, 2021).

For the static WD scenario, the total mass ejected during each nova eruption
can be found by taking the mass accretion rate of KT Eri ($3.57\times 10^{-7}\
\text{M}_{\odot}\ \text{yr}^{-1}$; Schaefer et al., 2022) over the course of
an inter-eruption period. The recurrence period of KT Eri is predicted to be
40 – 50 years by Schaefer et al. (2022), however, with only one observed
eruption of KT Eri, other recurrence periods are viable. As such, we chose to
sample a range of recurrence periods: 5, 10, 20, 50 and 100 years.

For the velocity of the nova ejecta, Arai et al. (2013) report that during its
nebular phase, KT Eri exhibited [O iii] $\lambda 4959$ and $\lambda 5007$
lines composed of multiple components with velocities of $-2000\ \text{km}\
\text{s}^{-1}$, $-1000\ \text{km}\ \text{s}^{-1}$, $+700\ \text{km}\
\text{s}^{-1}$ and $+1800\ \text{km}\ \text{s}^{-1}$. Additionally, Yamaoka et
al. (2009) report broad Balmer lines with the H$\alpha$ emission line FWHM
being approximately $3200-3400\ \text{km}\ \text{s}^{-1}$ alongside broad
emission lines in the $0.9\,\mu$m to $2.5\,\mu$m regime with velocities of
$4000\ \text{km}\ \text{s}^{-1}$. KT Eri is classified as a fast He/N nova,
therefore as well as an ejecta velocity of $4000\ \text{km}\ \text{s}^{-1}$
(in line with the velocities derived from spectra), we also chose to model
eruptions with feasible ejecta velocities of $5000\ \text{km}\ \text{s}^{-1}$
and $6000\ \text{km}\ \text{s}^{-1}$.

The ejecta mass, recurrence period and ejecta velocities outlined were adopted
to construct the characteristics of the nova ejecta in each of our models. The
combinations of parameters for Runs 1 – 15 are detailed in Table 1.

Figure 1: Left: The dynamics of the reference simulation (Run 4) NSR with
$3.57\times 10^{-5}\ \text{M}_{\odot}\ \text{yr}^{-1}$ and $n=2\times 10^{-3}$
after ${\sim}51$ kyr (1019 eruptions) with 1000 au resolution. The regions of
interest we refer to in the text are labelled. Right: The evolution of the
inner edge of the ejecta pile-up region and the inner and outer edges of the
NSR shell with respect to cumulative time. The estimated radius of the
observed shell around KT Eri is indicated.

###  2.4 Multicomponent ejecta

Nova eruptions are comprised of multiple components of ejecta: material being
ejected initially at slower velocities followed by faster ejecta at later
times (O’Brien et al., 1994). This leads to intra-ejecta shocks that heavily
contribute to the total optical luminosity of the nova event (Aydi et al.,
2020a, b; Murphy-Glaysher et al., 2022).

Accordingly, the second set of simulations (Runs 1⋆ – 5⋆) mimic the broad
characteristics of Runs 1 – 5, however the single bulk ejecta detailed in
Section 2.3 are replaced by a multicomponent ejecta. Specifically, each
eruption is broken down into twelve distinct components, with linearly
decreasing ejecta velocities from 6000 km s-1 to 0 km s-1 over the same
timescale as the 3 mag decline time ($t_{3}\sim 13.6$d).

The total kinetic energy contained within the nova eruptions influences the
growth of the NSR and its ultimate size (Healy, 2021; Healy-Kalesh et al.,
2023). Therefore, we scaled the total ejected mass from the combined
components of an eruption such that its total KE matched the total KE of the
eruption with bulk ejecta (as in Section 2.3). The details of these runs are
also presented in Table 1.

##  3 Results

###  3.1 Reference simulation: Run 4

####  3.1.1 Dynamics

We have selected Run 4 as our reference simulation given that the recurrence
period of 50 years is closest to the predicted value as given in Schaefer et
al. (2022). Each eruption from this system ejects $1.785\times 10^{-5}\
\text{M}_{\odot}$ at a velocity of $6000\ \text{km}\ \text{s}^{-1}$ into a
prepopulated ISM with a density of $n=2\times 10^{-3}$ ($3.34\times 10^{-27}\
\text{g}\ \text{cm}^{-3}$). We find that it would take this nova system
$5.1\times 10^{4}$ years over 1019 eruptions to grow a NSR to 25 pc – the
radial size of the shell surrounding KT Eri (Shara et al. 2023).

In the left-hand plot of Figure 1, we show the dynamics of the NSR after the
full $5.1\times 10^{4}$ years including its radial density, pressure,
temperature and velocity distribution. We also indicate the locations of the
cavity, ejecta pile-up region and the NSR shell (confined by its inner and
outer edge) that we refer to throughout the paper. The radial growth curve of
the outer and inner of the shell and the ejecta pile-up boundary in provided
in the right-hand plot of Figure 1.

We can see in the top left-hand panel of the left plot of Figure 1 that,
alongside the outer edge of the NSR shell extending to a radius of 25 pc, the
inner edge of the shell reaches out to ${\sim}21.6$ pc; this corresponds to a
shell thickness of ${\sim}13.7$ per cent. As evident in the right-hand plot of
Figure 1, the thickness of the NSR shell remains approximately constant
throughout its full evolution: e.g., $14.2$ per cent and $13.6$ per cent after
100 eruptions (${\sim}4,950$ yr) and 500 eruptions (${\sim}24,970$ yr),
respectively. A constant NSR shell thickness is reminiscent of the NSR grown
around the 12a system in Darnley et al. (2019) and during the early ‘steady
state’ evolution from (Healy-Kalesh et al., 2023) .

Near to the central nova, the density is high as individual eruptions eject
mass into the established structure (see top left-hand panel of the left plot
in Figure 1). Yet beyond this, the density within the cavity drops to $n\simeq
3.4\times 10^{-6}$ (${\sim}600\times$ less dense than the surrounding ISM
density) at the ejecta pile-up boundary due to the high velocity eruptions
traversing the region with very little resistance. With such low densities in
the cavity, the newly ejected material experiences close to free-expansion and
therefore maintain velocities of ${\sim}5800\ \text{km}\ \text{s}^{-1}$ out to
${\sim}7.8$ parsec, as illustrated in the bottom left-hand panel of the left
plot in Figure 1. The pressure (see top right-hand panel) and temperatures
(see bottom right-hand panel) also drop precipitously across the cavity in
response to the sparse levels of material here.

After crossing the low-density cavity, the ejecta collide with the ejecta
pile-up region at ${\sim}7.8$ parsec, resulting in velocities dropping
considerably to ${\sim}1800\ \text{km}\ \text{s}^{-1}$. As shown in the left-
hand plot of Figure 1, while the density increases approximately four-fold at
the border, the pressure and temperatures increase substantially to $1.2\times
10^{-12}\ \text{g}\ \text{cm}^{-3}\ \text{s}^{-2}$ and $7.1\times 10^{8}$ K
(both over four orders of magnitude), respectively, as collisions between
incoming ejecta with the established pile-up region violently shock-heats the
material.

Away from the extreme interface of the cavity and pile-up region, the density
of material within the ejecta pile-up region remains relatively constant
(${\sim}2-5\times 10^{-29}\ \text{g}\ \text{cm}^{-3}$) out to the inner edge
of the NSR shell, as does the pressure ($1.2-1.6\times 10^{-12}\ \text{g}\
\text{cm}^{-3}\ \text{s}^{-2}$) and temperature ($4-7\times 10^{8}$ K), though
the temperature of the material drops appreciably ${\sim}5$ pc prior to the
shell. On the other hand, due to the ejecta interacting with both material
from preceding eruptions and reverse shock fronts, the velocity of the
material within the pile-up region drops continuously from ${\sim}1800\
\text{km}\ \text{s}^{-1}$ down to ${\sim}250\ \text{km}\ \text{s}^{-1}$ within
the NSR shell.

The NSR shell, as indicated in the top left-hand panel of the left plot in
Figure 1, is a high density (${\sim}1.3\times 10^{-26}\ \text{g}\
\text{cm}^{-3}$ at the outer edge) region comprised almost exclusively of ISM
material swept up by the frequent highly energetic nova eruptions. While only
${\sim}4\times$ more dense than the surrounding ISM, the density of the shell
is up to ${\sim}650\times$ and ${\sim}2300\times$ greater than the pile-up
region and cavity it encompasses, respectively. In the vicinity closest to the
inner edge of the shell, the temperature drops abruptly from ${\sim}2.3\times
10^{8}$ K to ${\sim}5.9\times 10^{6}$ K (at the inner edge of the shell) in
the space of approximately 0.5 pc. The temperature gradient between the inner
and outer edge (${\sim}1.9\times 10^{6}$ K) stems from the material within the
shell being more and more shielded from the high energy ejecta impacting the
inner edge, and therefore is able to cool. There is a small increase in
pressure near to the outer edge (${\sim}2\times 10^{-12}\ \text{g}\
\text{cm}^{-3}\ \text{s}^{-2}$) with the velocity remaining approximately
${\sim}250\ \text{km}\ \text{s}^{-1}$ across the whole shell.

We explore the role of spatial resolution in this study by re-simulating Run 4
with a resolution of 40 au/cell; in Figure 2 we compare this higher resolution
simulation with Run 4 (1000 au/cell). The spatial resolution chosen for this
study does not have an impact on the large-scale structure of the NSR. The
most significant difference seen in Figure 2 is the expected greater number of
resolved individual ejecta across the cavity and pile-up boundary.

Figure 2: Comparing the dynamics of Run 4 (resolution of 1000 au/cell) and Run
4 with a resolution of 40 au/cell.

####  3.1.2 X-ray luminosity

Figure 3: Synthetic X-ray luminosity (without absorption) of the NSR grown in
Run 4 (black), along with the density distribution of the NSR (grey). The
dashed line indicates an interpolation as described in the text.

To predict the X-ray luminosity of the NSR grown in Run 4, we adopted the same
method as used previous work (Vaytet, 2009; Vaytet et al., 2011; Darnley et
al., 2019; Healy-Kalesh et al., 2023). This involves computing the emission
measure contribution from each spherical shell of the NSR in Run 4 before
binning these contributions into logarithmically divided temperature bins. We
then pass these temperature-binned emission measures into XSPEC while
employing the APEC (Smith et al., 2001) model (without absorption) to obtain
the X-ray luminosity. The radial profile of the NSR X-ray luminosity is shown
in Figure 3. We have interpolated the likely X-ray luminosity for the region
between ${\sim}7.8$ pc and ${\sim}13.4$ pc (represented as a dashed line in
Figure 3) due to the temperatures in this region being too high ($>64$ keV)
for the APEC model to predict the X-ray emission.

We can see in Figure 3 that the X-ray luminosity of the NSR from Run 4 is
relatively high (between $\sim 10^{20-25}\ \text{erg}\ \text{s}^{-1}$) at the
centre as a result of the highly energetic nova eruptions. Beyond ${\sim}2$
pc, the X-ray luminosity drops away dramatically as the ejecta travel
unimpeded through the low density cavity. Once the ejecta crashes into the
high density pile-up region border (at ${\sim}7.8$ pc), significant shock-
heating occurs leading to a huge increase in X-ray emission ($10^{24}\
\text{erg}\ \text{s}^{-1}$). Shock-heating across the whole ejecta pile-up
region maintains the high X-ray luminosity up to the inner edge of the NSR
shell. Here, the density increases by over two orders of magnitude – this
triggers another substantial jump in X-ray emission, reaching a peak of
$\text{L}_{\text{X-ray}}\simeq 10^{30}\ \text{erg}\ \text{s}^{-1}$ at the edge
of the NSR shell.

The total X-ray luminosity from the whole NSR is
$\text{L}_{\text{X-ray}}=6.5\times 10^{32}\ \text{erg}\ \text{s}^{-1}$, and so
is much brighter than the predicted X-ray luminosities of the NSR surrounding
12a: $3\times 10^{29}\ \text{erg}\ \text{s}^{-1}$ (Darnley et al., 2019) and
$1\times 10^{31}\ \text{erg}\ \text{s}^{-1}$ (Healy-Kalesh et al., 2023). Even
still, $\text{L}_{\text{X-ray}}$ for Run 4 is over four orders of magnitude
fainter than the typical X-ray luminosities of novae during their super-soft
source phase (see, e.g., Henze et al., 2010, 2011), however, does exceed the
X-ray luminosity of some quiescent novae (see, e.g., RS Oph; Page et al.,
2022) therefore allowing for potential detection.

####  3.1.3 H$\alpha$ surface brightness

Figure 4: As in Figure 3 but with the simulated H$\alpha$ surface brightness.

As with the NSR around 12a (Darnley et al., 2019), the shell surrounding KT
Eri is visible through its prominent H$\alpha$ emission (Shara et al. 2023).
We compute the H$\alpha$ surface brightness of Run 4 by firstly following the
method described in Healy-Kalesh et al. (2023). Here, we place the simulated
NSR from Run 4 at the distance of KT Eri (${\sim}5110$ pc; Schaefer et al.,
2022) and applied the H$\alpha$ extinction ($A_{\text{H}\alpha}=0.208$) toward
the nova. The total predicted H$\alpha$ luminosity of the NSR from this radial
distribution for Run 4 is $\text{L}_{\text{H}\alpha}\simeq 1.1\times 10^{33}\
\text{erg}\ \text{s}^{-1}$ and the vast majority of this emission emanates
from the NSR shell.

The simulated H$\alpha$ surface brightness distribution, assuming a spherical
geometry, is shown in Figure 4. At all radii, the dominant contribution to the
surface brightness is the projected NSR shell. The surface brightness close to
the nova is ${\sim}6\times 10^{-16}\ \text{erg}\ \text{s}^{-1}\
\text{cm}^{-2}\ \text{arcsec}^{-2}$ and gradually increases with increasing
radius, reaching ${\sim}1.3\times 10^{-15}\ \text{erg}\ \text{s}^{-1}\
\text{cm}^{-2}\ \text{arcsec}^{-2}$ at the NSR shell inner edge and then peaks
at ${\sim}2.1\times 10^{-15}\ \text{erg}\ \text{s}^{-1}\ \text{cm}^{-2}\
\text{arcsec}^{-2}$ within the NSR shell at ${\sim}23.7$ pc. Any H$\alpha$
contribution from ISM beyond, or in front or behind, the shell is not
included. We note that the predicted H$\alpha$ surface brightness of the shell
is ${\sim}20-50$ times larger than that reported by Shara et al. (2023). This
discrepancy may result from the cooling package employed within our
simulations (see Section 2.1) whereby the cooling rates are potentially
overestimated for densities $\sim 1\ \text{cm}^{-3}$ and not scaled for the
lower ISM densities. A number of other factors could also have affected this
estimate, including the ISM density, the assumption that the ISM is purely
hydrogen, the assumed eruption history, but, most importantly, the assumption
of spherical symmetry of both the ejecta and the ISM.

Figure 5: Comparing the dynamics of the NSRs grown in Runs 1 – 15. Top left:
Dynamics of remnants grown in Runs 1 – 5 (including the reference simulation,
Run 4) with a constant $v_{\text{ej}}=6000\ \text{km}\ \text{s}^{-1}$ and
varying $\text{P}_{\text{rec}}$. Top right: Dynamics of remnants grown in Runs
6 – 10 with a constant $v_{\text{ej}}=5000\ \text{km}\ \text{s}^{-1}$ and
varying $\text{P}_{\text{rec}}$. Bottom left: Dynamics of remnants grown in
Runs 11 – 15 with a constant $v_{\text{ej}}=4000\ \text{km}\ \text{s}^{-1}$
and varying $\text{P}_{\text{rec}}$. Bottom right: Dynamics of remnants grown
in Runs 4, 9 and 14 with a constant $\text{P}_{\text{rec}}=50$ years and
varying ejecta velocity.

###  3.2 Bulk ejecta: Runs 1 – 3, 5 – 15

After establishing the dynamics of a reference simulation (Run 4), we explored
NSR evolution from a KT Eri-like system with varying recurrence periods and
ejecta velocities: Runs 1–3 and Runs 5–15. The parameters used in these runs
are provided in Table 1 and we show the results of Runs 1 – 15 (with the
inclusion of the reference simulation) in Figure 5.

It is immediately apparent from the top left ($6000\ \text{km}\ \text{s}^{-1}$
in Runs 1 – 5), top right ($5000\ \text{km}\ \text{s}^{-1}$ in Runs 6 – 10)
and bottom left ($4000\ \text{km}\ \text{s}^{-1}$ in Runs 11 – 15) plot of
Figure 5 that varying the recurrence period does not affect the large-scale
structure of the NSR. Specifically, the cavity/pile-up region boundary and the
inner and outer edge of the shell extend out to the same radius in all runs
with the same ejecta velocity. However, the longer the period between
eruptions, the more structure we see in the cavity (along with a marginal
amount of fluctuation in the inner ejecta pile-up region) as a result of
resolving the individual eruptions. Within the cavity, the average density and
velocity remain the same as the recurrence period is altered, yet the pressure
and therefore temperature within the region increase moderately as
$\text{P}_{\text{rec}}$ becomes longer. Furthermore, the higher velocity
ejecta drive larger fluctuations in density, pressure and temperature, likely
because of the higher levels of KE being ejected from each nova eruption.

The bottom right plot of Figure 5 illustrates the fully grown NSR from Runs 4,
9 and 14: systems with the same recurrence period of 50 years but with ejecta
velocities of $6000\ \text{km}\ \text{s}^{-1}$, $5000\ \text{km}\
\text{s}^{-1}$ and $4000\ \text{km}\ \text{s}^{-1}$, respectively. As
previously found, the inner and outer edge of the NSR shells all extend to the
same distance of ${\sim}21.6$ pc and 25 pc, respectively, corresponding to a
thickness ${\sim}13.7$ per cent (as in the reference simulation in Section
3.1.1). Though, in contrast to the similarity in the dynamics of the NSRs when
sampling the recurrence period (with only the cavity showing minor
differences), varying the ejecta velocity produces more noticeable changes (as
shown in the bottom right-hand plot of Figure 5).

Figure 6: Comparing the dynamics of the NSRs grown in Runs 1⋆ – 5⋆. Left:
Dynamics of the NSR grown in the Run 4 (reference simulation) compared to Run
4⋆. Run 4 is denoted as the model with bulk ejecta and Run 4⋆ is denoted as
the model with multicomponent ejecta. Right: Dynamics of remnants grown with
multicomponent ejecta with varying recurrence periods (1⋆ – 5⋆).

Specifically, the radial extent of the cavity/pile-up region boundary
increases from ${\sim}7.8$ pc to ${\sim}8.1$ pc to ${\sim}8.3$ pc as the
ejecta velocity decreases ($6000\ \text{km}\ \text{s}^{-1}$ to $5000\
\text{km}\ \text{s}^{-1}$ to $4000\ \text{km}\ \text{s}^{-1}$) – this likely
results from the increasing number of eruptions (and therefore the longer
evolutionary times) needed to reach the radius of 25 pc also allowing more
time to move the cavity border further. As we would expect, the velocity of
the ejecta traversing the cavity is lower for Run 14 (${\sim}3860\ \text{km}\
\text{s}^{-1}$) than Run 4 (${\sim}5800\ \text{km}\ \text{s}^{-1}$) on account
of the lower initial velocity for the former – the small drop in ejecta
velocity seen in all runs is directly related to the small (but non-
negligible) resistance of the lower density material in the forming cavity.
The difference in the velocity profile for Runs 4, 9 and 14 is maintained
throughout the pile-up region and across the NSR shell.

In the bottom right-hand panel of the bottom right-hand plot in Figure 5 we
see that the temperature across the whole NSR is cooler for lower velocity
ejecta. For example, the temperature of the inner ejecta pile-up (on the
border with the cavity) for Run 14 (the cyan line) is ${\sim}3.2\times 10^{8}$
K whereas the same region in Run 4 (the red line) reaches ${\sim}7.1\times
10^{8}$ K. This large difference in temperature is also found in the NSR
shell, with the lower velocity ejecta leading to cooler shells.

###  3.3 Multi-component ejecta

Here, we consider the NSRs grown with ejecta composed of multiple components.
As detailed in Table 1, and by design (as the KE of each multi-component
ejecta matches the KE of the bulk ejecta), Run 4⋆ takes the same time
(${\sim}5.1\times 10^{4}$ years) as Run 4 (both with
$\text{P}_{\text{rec}}=50$ yr) to grow a NSR to match the size of observed KT
Eri shell. The nova in this model also exhibits an identical number of
eruptions (1,019) to match the NSR in Run 4 and as a result ejects the amount
of kinetic energy as the reference simulation.

We grew the NSR in Run 4⋆ until it reached the size of the observed shell
presented in Shara et al. (2023) as we have with all other models in this
study. Yet, despite reaching this same size (and having the same shell
thickness) as Run 4, the NSR in Run 4⋆ displays appreciably different
dynamics. As illustrated in the top left panel of the left-hand plot of Figure
6, the boundary between the cavity and the ejecta pile-up region is situated
${\sim}1.2$ parsecs further from the central nova compared to the equivalent
border in Run 4. Furthermore, the density within the cavity and pile-up region
is approximately one order of magnitude higher in the NSR grown with
multicomponent ejecta compared to the NSR grown with bulk ejecta. While the
pressure and velocity are identical across the ejecta pile-up and NSR shell,
the pressure in the cavity is higher ($\times 2.5$) in Run 4⋆ compared to Run
4 and the velocity is lower ($\times 0.75$). As shown in the bottom right
panel of the left-hand plot in Figure 6, the temperature in the cavity and
shell are very similar whereas there is a large drop in the temperature of the
whole pile-up region. The differences in the dynamics outlined above arise
from the lower density components of the ejecta in Run 4⋆ being more
susceptible to radiative cooling and also the slower gas, which takes longer
to traverse the cavity.

In the right-hand plot of Figure 6, we show the dynamics of the remnants grown
with multicomponent ejecta with varying recurrence periods in (Runs 1⋆ – 5⋆;
see Table 1 for details). In a similar manner to the NSRs grown with bulk
ejecta eruptions of the same ejecta velocity, we see that the remnants in Runs
1⋆ – 5⋆ all share the same (almost identical) structure beyond the ejecta
pile-up boundary. While the velocity remains the same in the cavity for all
multicomponent ejecta runs, the pressure and temperature in this region do
vary between runs. Additionally, the density within the cavity regions are
largely different: the shorter the recurrence period, the smoother the density
distribution. Large fluctuations apparent in the NSR cavity grown from a nova
with a long recurrence period reveal individual nova ejecta traversing this
region before their collision with the ejecta pile-up boundary.

Figure 7: Comparing the dynamics of the NSR from Run 4 (reference simulation
with smoothly distributed eruptions) with the dynamics of the NSR from Run 4†
with Gaussian distributed eruptions.

###  3.4 Randomly occurring eruptions

For Runs 1 – 15 and Runs 1⋆ – 5⋆, we have a recurrent nova evolving with
constant recurrence period. However, known RNe deviate from their average
recurrence period by up to 10 per cent (Henze et al., 2018; Darnley, 2021). To
implement this feature, we allowed each inter-eruption period to differ by
randomly selecting the inter-eruption times from a Gaussian distribution with
a mean recurrence period of 50 years ($\mu=50$) and a standard deviation of 10
per cent ($\sigma=0.1$) and constructed the corresponding nova ejecta.

As shown in Figure 7, random inter-eruption periods selected from a Gaussian
distribution ($\mu=50$, $\sigma=0.1$) between eruptions has no significant
impact on the structure of the NSR (Run 4†) compared to the reference
simulation (Run 4; with identical inter-eruptions times). The number of
eruptions (1,033) and total time of evolution ($5.2\times 10^{4}$ years)
closely match the values for the reference simulation (see Table 1), as does
the total kinetic energy. As the same amount of energy is being injected from
the central nova into the surrounding ISM, it is inevitable that a NSR with a
radius of 25 pc is created in a similar time frame.

The one difference between the two NSRs is the development of some structure
in the cavity and inner ejecta-pile up region. There are small fluctuations in
density, pressure and temperature in the cavity region whereas the more
significant fluctuations in velocity are located within the ejecta-pile up
region. As the material is being ejected from the nova at non-identical time
intervals in this scenario, a proportion of the ejecta will catch the
preceding ejecta at random points of the NSR - this gives rise to the non-
smooth distribution of material within the structure. However, as with all the
other simulations, the structure of the NSR shell is unaffected.

##  4 Discussion

For the recurrence periods sampled in our study, we find that the NSR shell
grown from ejecta with the same velocity reach a radius of 25 parsecs in the
same amount of time, with only the number of eruptions changing. As such, the
time taken for the NSR to grow to the observed size of the shell found
surrounding KT Eri (Shara et al. 2023) in our reference simulation (${\sim}50$
kyr) is fixed regardless of the period between eruptions. A 50 kyr timescale
also supports the assumption that the WD in the KT Eri system is relatively
young compared to other RNe.

Likewise, the shell thickness of ${\sim}14$ per cent for the reference
simulation NSR is consistent across all runs employing eruptions with bulk
ejecta. Though the ISM surrounding KT Eri is likely to be non-uniform and so
shapes the observed NSR into a complex structure, the predicted thickness of
${\sim}14$ per cent from our model is close to the thickness of northern edge
of the observed shell (Shara et al. 2023). Though, unlike in Healy-Kalesh et
al. (2023), we have not considered mass-loss occurring during prior phases of
KT Eri’s evolution. The compression of the NSR shell in the direction of
proper motion (see Fig 3 of Shara et al. 2023) may then be evidence of
interaction with a bow shock from this earlier phase. Also, magnetic fields,
not considered in this work, affect stellar winds which then have a large
impact on ISM structure (Wareing et al., 2016) and consequently the growth of
a NSR shell.

In Section 3.3, we investigated the impact of velocity variations within the
nova ejecta on the structure of the growing NSR. Even though we tuned the
model such that the final NSR shell extended to 25 parsecs, we did find that
the NSR substructure is considerably different compared to the NSR grown with
constant ejecta velocity eruptions. This is a consequence of the components
with slower, low-density gas being initially cooler and therefore losing
energy through radiative cooling more efficiently than the bulk ejecta in the
reference simulation (Run 4). In Healy-Kalesh et al. (2023), we demonstrated
that a ‘two-component’ (rather than multiple component as tested here) ejecta
did not have an discernible effect on the NSR structure as, in this scenario,
the ejecta still met the criteria whereby they were not able to cool
efficiently. The twelve-component ejecta (in Runs 1⋆ – 5⋆) are more complex
and therefore we see a more dramatic difference in the NSR sub-structure due
to more efficient cooling.

##  5 Conclusions

We have presented a set of simulations used to model a nova super-remnant
found in narrow-band imaging surrounding the Galactic classical nova, KT
Eridani (Shara et al. 2023). For the models with $v_{\text{ej}}=6000\
\text{km}\ \text{s}^{-1}$, a nova super-remnant was grown to match the radial
size of the observed structure (${\sim}25$ pc) in ${\sim}$50,000 years –
consistent with the proper motion history of KT Eri (Shara et al. 2023).
Additionally, the NSR shell thickness of ${\sim}14$ per cent found in the
models is in reasonable agreement with the observed shell.

An estimate of the integrated X-ray luminosity of the KT Eri nova super-
remnant suggests that the structure may be accessible to wide-field X-ray
facilities. Our prediction of the H$\alpha$ surface brightness of the NSR
shell is ${\sim}20-50$ times larger than that measured by Shara et al. (2023),
however, this will, at least, in part be due to several of our assumptions,
most importantly that of spherical symmetry. Exploring cooling conditions and
implementing complex ISM structure through differing mass-loss phases and
magnetic field shaping in future work will help to reconcile our model
predictions with observed parameters.

We strongly encourage further observations of the KT Eri nova super-remnant,
especially exploration of the predicted X-ray emission.

## Acknowledgements

The authors would like to thank our reviewer, Dr Christopher Wareing, for his
helpful suggestions that we will implement in future modelling. MWH-K
acknowledges a PDRA position funded by the UK Science and Technology
Facilities Council (STFC). MWH-K and MJD receive funding from STFC grant
number ST/S505559/1. MMS, KML and JTG acknowledge the support of NSF award
2108234. This work made use of the high performance computing facilities at
Liverpool John Moores University, partly funded by LJMU’s Faculty of
Engineering and Technology and by the Royal Society.

## Data Availability

The data in this study can be shared on reasonable request to the
corresponding author. This work was conducted with the Morpheus (Vaytet et
al., 2007) program and analysed using the Python libraries: Numpy (Harris et
al., 2020) and Matplotlib (Hunter, 2007).

## References

  * Arai et al. (2013) Arai A., Isoggai M., Imamura K., Ikeda Y., Arasaki T., Kitao E., Taguchi G., 2013, in Di Stefano R., Orio M., Moe M., eds, Proceedings of the International Astronomical Union Vol. 281, Binary Paths to Type Ia Supernovae Explosions. pp 119–120, doi:10.1017/S1743921312014809
  * Aydi et al. (2020a) Aydi E., et al., 2020a, Nature Astronomy, 4, 776
  * Aydi et al. (2020b) Aydi E., et al., 2020b, ApJ, 905, 62
  * Bode et al. (2004) Bode M. F., O’Brien T. J., Simpson M., 2004, ApJ, 600, L63
  * Darnley (2021) Darnley M. J., 2021, in The Golden Age of Cataclysmic Variables and Related Objects V. p. 44 (arXiv:1912.13209)
  * Darnley & Henze (2020) Darnley M. J., Henze M., 2020, Advances in Space Research, 66, 1147
  * Darnley et al. (2012) Darnley M. J., Ribeiro V. A. R. M., Bode M. F., Hounsell R. A., Williams R. P., 2012, ApJ, 746, 61
  * Darnley et al. (2015) Darnley M. J., et al., 2015, A&A, 580, A45
  * Darnley et al. (2016) Darnley M. J., et al., 2016, ApJ, 833, 149
  * Darnley et al. (2017) Darnley M. J., et al., 2017, ApJ, 849, 96
  * Darnley et al. (2019) Darnley M. J., et al., 2019, Nature, 565, 460
  * Gensior et al. (2023) Gensior J., Feldmann R., Mayer L., Wetzel A., Hopkins P. F., Faucher-Giguère C.-A., 2023, MNRAS, 518, L63
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Healy (2021) Healy M. W., 2021, PhD thesis, https://researchonline.ljmu.ac.uk/id/eprint/15515/1/2021healyphd.pdf
  * Healy-Kalesh et al. (2023) Healy-Kalesh M. W., Darnley M. J., Harvey É. J., Copperwheat C. M., James P. A., Andersson T., Henze M., O’Brien T. J., 2023, MNRAS, 521, 3004
  * Henze et al. (2010) Henze M., et al., 2010, A&A, 523, A89
  * Henze et al. (2011) Henze M., et al., 2011, A&A, 533, A52
  * Henze et al. (2018) Henze M., et al., 2018, ApJ, 857, 68
  * Hounsell et al. (2010a) Hounsell R., et al., 2010a, ApJ, 724, 480
  * Hounsell et al. (2010b) Hounsell R., et al., 2010b, The Astronomer’s Telegram, 2558, 1
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Itagaki (2009) Itagaki K., 2009, Central Bureau Electronic Telegrams, 2050, 1
  * Jackson et al. (1997) Jackson B., Buffington A., Hick P., Kahler S., Keil S., Altrock R., Simnett G., Webb D., 1997, Physics and Chemistry of the Earth, 22, 441
  * Kato et al. (2014) Kato M., Saio H., Hachisu I., Nomoto K., 2014, ApJ, 793, 136
  * Kato et al. (2015) Kato M., Saio H., Hachisu I., 2015, ApJ, 808, 52
  * Kurth et al. (2023) Kurth W. S., Burlaga L. F., Kim T., Pogorelov N. V., Granroth L. J., 2023, ApJ, 951, 71
  * Lanzetta et al. (2023) Lanzetta K. M., Gromoll S., Shara M. M., Berg S., Valls-Gabaud D., Walter F. M., Webb J. K., 2023, PASP, 135, 015002
  * Lloyd et al. (1997) Lloyd H. M., O’Brien T. J., Bode M. F., 1997, MNRAS, 284, 137
  * Murphy-Glaysher et al. (2022) Murphy-Glaysher F. J., et al., 2022, MNRAS, 514, 6183
  * O’Brien et al. (1994) O’Brien T. J., Lloyd H. M., Bode M. F., 1994, MNRAS, 271, 155
  * O’Brien et al. (2001) O’Brien T. J., Davis R. J., Bode M. F., Eyres S. P. S., Porter J. M., 2001, in Schilizzi R. T., ed., International Astronomical Union Vol. 205, Galaxies and their Constituents at the Highest Angular Resolutions. p. 260
  * Page et al. (2022) Page K. L., et al., 2022, MNRAS, 514, 1557
  * Pagnotta & Schaefer (2014) Pagnotta A., Schaefer B. E., 2014, ApJ, 788, 164
  * Raymond et al. (1976) Raymond J. C., Cox D. P., Smith B. W., 1976, ApJ, 204, 290
  * Schaefer (2010) Schaefer B. E., 2010, ApJS, 187, 275
  * Schaefer et al. (2022) Schaefer B. E., Walter F. M., Hounsell R., Hillman Y., 2022, MNRAS, 517, 3864
  * Shara et al. (1997) Shara M. M., Zurek D. R., Williams R. E., Prialnik D., Gilmozzi R., Moffat A. F. J., 1997, AJ, 114, 258
  * Shara et al. (2012a) Shara M. M., Mizusawa T., Zurek D., Martin C. D., Neill J. D., Seibert M., 2012a, ApJ, 756, 107
  * Shara et al. (2012b) Shara M. M., Mizusawa T., Wehinger P., Zurek D., Martin C. D., Neill J. D., Forster K., Seibert M., 2012b, ApJ, 758, 121
  * Slavin et al. (1995) Slavin A. J., O’Brien T. J., Dunlop J. S., 1995, MNRAS, 276, 353
  * Smith et al. (2001) Smith R. K., Brickhouse N. S., Liedahl D. A., Raymond J. C., 2001, ApJ, 556, L91
  * Starrfield et al. (1972) Starrfield S., Truran J. W., Sparks W. M., Kutter G. S., 1972, ApJ, 176, 169
  * Starrfield et al. (1976) Starrfield S., Sparks W. M., Truran J. W., 1976, in Eggleton P., Mitton S., Whelan J., eds, IAU Symposium Vol. 73, Structure and Evolution of Close Binary Systems. pp 155–172
  * Stil & Irwin (2001) Stil J. M., Irwin J. A., 2001, ApJ, 563, 816
  * Takeda et al. (2022) Takeda L., Diaz M., Campbell R. D., Lyke J. E., Lawrence S. S., Linford J. D., Sokolovsky K. V., 2022, MNRAS, 511, 1591
  * Toraskar et al. (2013) Toraskar J., Mac Low M.-M., Shara M. M., Zurek D. R., 2013, ApJ, 768, 48
  * Vaytet (2009) Vaytet N. M. H., 2009, PhD thesis, https://www.nbi.dk/~nvaytet/documents/thesis.pdf
  * Vaytet et al. (2007) Vaytet N. M. H., O’Brien T. J., Bode M. F., 2007, ApJ, 665, 654
  * Vaytet et al. (2011) Vaytet N. M. H., O’Brien T. J., Page K. L., Bode M. F., Lloyd M., Beardmore A. P., 2011, ApJ, 740, 5
  * Wade (1990) Wade R. A., 1990, in Cassatella A., Viotti R., eds, Lecture Notes in Physics, Berlin Springer Verlag Vol. 369, IAU Colloq. 122: Physics of Classical Novae. pp 179–187, doi:10.1007/3-540-53500-4_118
  * Walker (1954) Walker M. F., 1954, PASP, 66, 230
  * Wareing et al. (2006) Wareing C. J., O’Brien T. J., Zijlstra A. A., Kwitter K. B., Irwin J., Wright N., Greimel R., Drew J. E., 2006, MNRAS, 366, 387
  * Wareing et al. (2016) Wareing C. J., Pittard J. M., Falle S. A. E. G., Van Loo S., 2016, MNRAS, 459, 1803
  * Warner (1995) Warner B., 1995, Cataclysmic variable stars. Cambridge Astrophysics Series, Cambridge, New York: Cambridge University Press, 1995
  * Worters et al. (2007) Worters H. L., Eyres S. P. S., Bromage G. E., Osborne J. P., 2007, MNRAS, 379, 1557
  * Woudt & Ribeiro (2014) Woudt P. A., Ribeiro V. A. R. M., eds, 2014, Stella Novae: Past and Future Decades Astronomical Society of the Pacific Conference Series Vol. 490\. Astronomical Society of the Pacific, San Francisco
  * Yamaoka et al. (2009) Yamaoka H., et al., 2009, IAU Circ., 9098, 1
  * Yaron et al. (2005) Yaron O., Prialnik D., Shara M. M., Kovetz A., 2005, ApJ, 623, 398


# Physics-Informed Tensor Basis Neural Network for Turbulence Closure Modeling

Leon Riccius, Atul Agrawal & Phaedon-Stelios Koutsourelakis
Professorship of Data-driven Materials Modeling
School of Engineering and Design
Technical University of Munich
Boltzmannstr. 15, 85748 Garching, Germany
{leon.riccius, atul.agrawal<EMAIL_ADDRESS>

###### Abstract

Despite the increasing availability of high-performance computational
resources, Reynolds-Averaged Navier-Stokes (RANS) simulations remain the
workhorse for the analysis of turbulent flows in real-world applications.
Linear eddy viscosity models (LEVM), the most commonly employed model type,
cannot accurately predict complex states of turbulence. This work combines a
deep-neural-network-based, nonlinear eddy viscosity model with turbulence
realizability constraints as an inductive bias in order to yield improved
predictions of the anisotropy tensor. Using visualizations based on the
barycentric map, we show that the proposed machine learning method’s
anisotropy tensor predictions offer a significant improvement over all LEVMs
in traditionally challenging cases with surface curvature and flow separation.
However, this improved anisotropy tensor does not, in general, yield improved
mean-velocity and pressure field predictions in comparison with the best-
performing LEVM.

##  1 Introduction

The incompressible Navier-Stokes equations are vital for describing fluid
motion at low Reynolds numbers, impacting fields like aircraft design and
ocean current modeling. Turbulent flows are prohibitively expensive to resolve
fully, and engineers often resort to reduced models such as RANS for
efficiency. These models employ closures such as the Launder-Sharma
$k-\epsilon$ Launder and Sharma [1974] or Wilcox’s $k-\omega$ Wilcox [2008],
which rely on linear assumptions, limiting their accuracy in complex flow
scenarios. Nonlinear models have been explored but face challenges. This
contribution builds upon the resurgence of turbulence modeling research,
instigated by data-driven approaches Duraisamy et al. [2019], in response to
the stagnation seen in the 2000s after earlier advancements.

This work combines additional flow features derived by Wang and colleagues
Wang and Xiao [2016], Wang et al. [2017b] with the neural network architecture
proposed by Ling et al. Ling et al. [2016] to give point-based estimates of
the anisotropy tensor appearing in the RANS closure. The training objectives
are supplemented by a loss term penalizing predictions that violate the
physical realizability constraints of turbulent states. The trained network
proposed was tested on unseen flow scenarios and used as a source term in the
RANS equations to produce estimates of the mean-flow quantities.

##  2 Physics-Informed Tensor Basis Neural Network

###  2.1 Reynolds stresses and realizability constraints

The Reynolds stress tensor $\tau_{ij}=\left\langle u_{i}u_{j}\right\rangle$,
where $u_{i}=U_{i}-\left\langle U_{i}\right\rangle$ arises by time-averaging
(mean indicated by $\langle\cdot\rangle$) of the Navier-Stokes equations and
depends on the fluctuations $u_{i}$ of the velocity field $U_{i}$. It can be
decomposed into an isotropic $\delta_{ij}$ and anisotropic $a_{ij}$ part which
is given by $a_{ij}=\tau_{ij}-\nicefrac{{2}}{{3}}k\delta_{ij}$.The turbulent
kinetic energy $k$ is the trace of the Reynolds stresses.

The anisotropic $a_{ij}$ has zero trace, and its normalized version $b_{ij}$
is referred to as anisotropy tensor, i.e.:

$\displaystyle b_{ij}=\frac{a_{ij}}{2k}=\frac{\tau_{ij}}{\left\langle
u_{k}u_{k}\right\rangle}-\frac{1}{3}\delta_{ij},$ (1)

The Reynolds stress tensor is a symmetric, positive semi-definite second-order
tensor with a non-negative determinant and trace. Following Schumann [1977],
the physical constraints, known also as realizability constraints, on the
anisotropy tensor are:

$\displaystyle-\frac{1}{3}\leq
b_{\alpha\alpha}\leq\frac{2}{3}\quad\forall\alpha\in\\{1,2,3\\},\quad\quad-\frac{1}{2}\leq
b_{\alpha\beta}\leq\frac{1}{2}\quad\forall\alpha\neq\beta.$ (2)

The barycentric map introduced in Banerjee et al. [2007] uses an eigenvalue
decomposition of $\boldsymbol{b}$ to define three fundamental states of
turbulence. All other states of turbulence can be expressed as a linear
combination of these three limiting states. The limiting states form a
triangle of all admissible states; its vertices are defined by the
realizability constraints (2). As demonstrated in Figure 1, each point in the
barycentric triangle corresponds to a unique color. The mapping is given in
the supplementary material 6.1.

(a) Characteristic regions

(b) RANS $k-\omega$

(c) DNS

(d) RGB mapping

Figure 1: Barycentric map representing nature of turbulence of RANS $k-\omega$
(a) and DNS (b)

###  2.2 Closure Model

While LEVMs assume $\boldsymbol{b}$ to be a linear function of the mean
velocity gradient, a more general class of turbulence models can be formulated
when dropping this assumption. The class of algebraic stress models is formed
by nonlinear eddy viscosity models (NLEVM), which determine the Reynolds
stresses from the local turbulent kinetic energy $k$, the eddy viscosity
$\epsilon$, and the mean velocity gradient. Pope Pope [1975] has shown that
every second-order tensor that can be formed from the normalized mean rate of
strain
$\hat{\boldsymbol{S}}=\nicefrac{{\epsilon}}{{2k}}(\nabla\left\langle\boldsymbol{U}\right\rangle+\nabla\left\langle\boldsymbol{U}\right\rangle^{T})$
and the normalized rate of rotation
$\hat{\boldsymbol{\Omega}}=\nicefrac{{\epsilon}}{{2k}}(\nabla\left\langle\boldsymbol{U}\right\rangle-\nabla\left\langle\boldsymbol{U}\right\rangle^{T})$
and fulfills these requirements is a linear combination of ten basis tensors
$\mathcal{T}_{ij}^{(n)}$. The most general form of a NLEVM is given by

$\displaystyle
b_{ij}=\sum_{n=1}^{10}G^{(n)}(\lambda_{1},...,\lambda_{5})\,{\mathcal{T}}_{ij}^{(n)}(\hat{S}_{ij},\hat{\Omega}_{ij}),$
(3)

where $G^{(n)}$ are the coefficients of the basis tensors and $\lambda_{k}$
are the tensor invariants dependent on $\hat{\boldsymbol{S}}$ and
$\hat{\boldsymbol{\Omega}}$. Ling et al. Ling et al. [2016] introduced the
Tensor Basis Neural Network (TBNN), which makes use of modern machine learning
methods to learn these functions $G^{(n)}$ from high-fidelity fluid simulation
data. Even though improved results compared to the $k-\epsilon$ model were
reported, extracting enough information from these five invariants has proven
difficult. This is especially true for flow cases with at least one direction
of homogeneity, where invariants $\lambda_{3}$ and $\lambda_{4}$ vanish for
the entire flow domain. Proof of this statement is given in the supplementary
material 6.2. It is, however, possible to include more features from local
flow quantities and derive more invariants while still employing the integrity
basis formed by $\mathcal{T}^{(n)}$. This work, in part, follows the research
of Wang et al. Wang et al. [2017a], who derived an extended feature set also
considering the gradients of the turbulent kinetic energy and the pressure.
The model was implemented in PyTorch Paszke et al. [2019] and can be accessed
via github.com/pkmtum/PI_TBNN.

####  2.2.1 Enforcing Realizability Constraints in Training

Ling et al. Ling et al. [2016] enforced the realizability constraints in post-
processing by simply projecting the points onto the closest boundary of the
barycentric triangle. We propose incorporating these constraints in the TBNN’s
training, which we then colloquially refer to as the Physics-Informed Tensor
Basis Neural Network (PI-TBNN). The inequalities given in Eq. (2) can be
transformed into contributions to the loss function via the penalty method.
The additional term reflects inductive bias about the problem structure and
essentially acts as a regularizer. In plain words, if training samples are
outside the domain of realizable turbulence states, the penalty term will
force them back in. The constraints are

$\begin{split}c_{1}(\boldsymbol{b})&=\min_{\alpha}(b_{\alpha\alpha})-\nicefrac{{1}}{{3}}<0\quad\forall\alpha\in\\{1,2,3\\},\\\
c_{2}(\boldsymbol{b})&=({3|\phi_{2}|-\phi_{2}})/2-\phi_{1}<0,\\\
c_{3}(\boldsymbol{b})&=\nicefrac{{1}}{{3}}-\phi_{2}<0,\end{split}\hskip
11.38092pt\begin{split}c_{4}(\boldsymbol{b})&=2|b_{12}|-\left(b_{11}+b_{22}+\nicefrac{{2}}{{3}}\right)<0,\\\
c_{5}(\boldsymbol{b})&=2|b_{13}|-\left(b_{11}+b_{33}+\nicefrac{{2}}{{3}}\right)<0,\\\
c_{6}(\boldsymbol{b})&=2|b_{23}|-\left(b_{22}+b_{33}+\nicefrac{{2}}{{3}}\right)<0,\\\
\end{split}$ (4)

where $\phi_{i}$ denotes the eigenvalues of $\boldsymbol{b}$ in order to
distinguish them from the invariants. $\phi_{1}$ and $\phi_{2}$ are the
largest and second-largest eigenvalues. The penalty term is then given by

$\displaystyle\mathcal{L}(\hat{\boldsymbol{b}}(\boldsymbol{\lambda},\boldsymbol{\theta}))=\beta\sum_{k=1}^{6}\max(0,c_{k}(\hat{\boldsymbol{b}}(\boldsymbol{\lambda},\boldsymbol{\theta}))),$
(5)

where $\boldsymbol{\lambda}$ is the collection of invariants,
$\boldsymbol{\theta}$ are the NN parameters, and
$\boldsymbol{\hat{b}}(\boldsymbol{\lambda},\boldsymbol{\theta})$ is the
predicted anisotropy tensor. The penalty coefficient $\beta$ determines its
impact on the loss function. The complete loss function, considering the MSE
loss, the regularization, and penalty terms, is given by

$\displaystyle E(\boldsymbol{\lambda}_{i},\boldsymbol{\theta})$
$\displaystyle=\frac{1}{D}\sum_{i=1}^{D}\|\boldsymbol{\hat{b}}(\boldsymbol{\lambda}_{i},\boldsymbol{\theta})-\boldsymbol{b}_{i}\|+\frac{\alpha}{2}\|\boldsymbol{\theta}\|^{2}_{2}+\frac{\beta}{D}\sum_{i=1}^{D}\sum_{k=1}^{6}\max(0,c_{k}(\boldsymbol{\hat{b}}(\boldsymbol{\lambda}_{i},\boldsymbol{\theta}))),$
(6)

where $\boldsymbol{b}_{i}$ are the high-fidelity responses. The number of data
points is denoted $D$. The coefficient $\alpha$ controls the degree of $L_{2}$
regularization.

##  3 Numerical Results

A total of four flow geometries were used as benchmarks in this work. The flow
over periodic hills (PH) Mellen et al. [2000], the converging-diverging
channel flow (CDC) Laval and Marquillie [2011], and the curved backward-facing
step (CBFS) Bentaleb et al. [2012]. These exhibit adverse pressure gradients
over curved surfaces, leading to flow separation and subsequent reattachment.
The data set was also extended to include the square duct (SD) flow case
Pinelli et al. [2010]. This scenario clearly illustrates the limitations of
LEVMs and is well-suited for investigating the forward propagation of the
predicted anisotropy tensor. All flow cases were replicated in OpenFOAM Weller
et al. [1998] to obtain the baseline RANS data. The flow case setup is
analogous to Leon Riccius [2021]. A detailed description of the data sets is
given in the supplementary material 6.3.

###  3.1 Anisotropy Tensor Prediction

The PI-TBNN was tested on flow cases it had not seen during training. They
differ either in geometry or Reynolds number from the training data. Figure 2
compares the anisotropy tensors for square duct (SD) using the barycentric
colormap. Only the TBNN with the extended feature set can accurately reproduce
the state of turbulence for this flow case. A similar picture arises on the PH
geometry in Figure 3, where the $k-\omega$ model cannot capture 1C turbulence
and axisymmetric expansion at the top and the bulk of the flow domain,
respectively. The PI-TBNN, however, does exhibit such characteristics. For all
test cases considered, the PI-TBNN achieves about 70% reduction of the RMSE
compared to the baseline $k-\omega$ model and 50% reduction of the RMSE
compared to Ling et al. [2016] (see Table 1).

(a) LEVM $k-\omega$

(b) PI-TBNN1

(c) PI-TBNN2

(d) DNS

Figure 2: Stress types of $k-\omega$ (a), PI-TBNN with FS1 (b), FS[1-3] (c), and DNS (d) for SD. Table 1: RMSE of $\boldsymbol{b}$ from RANS, PI-TBNN, and TBNN predictions for the three test cases. Flow case |  LEVM $k-\omega$ |  PI-TBNN, FS1 |  PI-TBNN, FS[1-3] |  TBNN Ling Ling et al. [2016]
---|---|---|---|---
Square duct |  $0.2175$ |  $0.0992$ |  $0.0663$ |  $0.14$
Periodic hills |  $0.1016$ |  $0.0628$ |  $0.0419$ |
Curved backward-facing step |  $0.1173$ |  $0.0619$ |  $0.0414$ |

(a) LEVM $k-\omega$

(b) PI-TBNN

(c) DNS

Figure 3: Visualization of stress types of LEVM $k-\omega$ (a), PI-TBNN (b),
and DNS (c)

###  3.2 Anisotropy Tensor Propagation

While some applications involving wall shear stress computations may directly
benefit from an improved prediction of $\boldsymbol{b}$, the quantities of
interest are usually the mean velocity and pressure fields. Hence, with the
PI-TBNN model, the Reynolds equations were solved for the mean velocity and
pressure fields. Table 2 shows that the PI-TBNN outperforms the $k-\epsilon$
model for the PH and CBFS geometries by a small margin and yields better in-
plane prediction for the square duct flow case. The most accurate model,
however, remains the $k-\omega$ model for all three test geometries.
Surprisingly, on the CBFS geometry, the PI-TBNN shows the largest
discrepancies in the region of the flow separation, as can be seen in Figure
4. The $k-\omega$ model even beats the mean field resulting from using the
anisotropy tensor from the DNS on the periodic hills test case, indicating
that an improved anisotropy tensor does not necessarily lead to improved mean
velocity and pressure fields (behavior also reported in Taghizadeh et al.
[2020], Duraisamy [2021]). Both Ling et al. [2016] and Wu et al. [2018]
reported improvements in the mean-field prediction over a LEVM but used the
$k-\epsilon$ as the baseline LEVM, which shows a larger discrepancy from the
ground truth than the $k-\omega$ model, i.e. the best performing model of its
class.

Table 2: RMSE of $\boldsymbol{U}$ for RANS with anisotropy tensor from $k-\omega$, $k-\epsilon$, PI-TBNN, and DNS. Reference velocity fields come from DNS. Flow case |  $k-\omega$ |  $k-\epsilon$ |  $\boldsymbol{b}_{\mathrm{PI-TBNN}}$ |  $\boldsymbol{b}_{\mathrm{DNS}}$
---|---|---|---|---
Square duct RMSE$(\boldsymbol{U})$ |  $0.0496$ |  $0.0667$ |  $0.0716$ |  $0.0322$
Square duct RMSE$([U_{2},U_{3}])$ |  $0.0066$ |  $0.0066$ |  $0.0045$ |  $0.0025$
Periodic hills RMSE$(\boldsymbol{U})$ |  $0.0375$ |  $0.0545$ |  $0.0541$ |  $0.0465$
Curved backward-facing step RMSE$(\boldsymbol{U})$ |  $0.0609$ |  $0.0883$ |  $0.0868$ |

(a) Periodic hills

(b) Curved backward-facing step

Figure 4: Streamwise mean velocity profiles at specific $x$-locations of PH
(a), and CBFS (b).

##  4 Conclusion

We introduced the PI-TBNN, which extended the TBNN framework with an extensive
feature set and an inductive bias in the form of a physics-informed addition
to the loss function. The addition of features was motivated both
analytically—showing that the number of distinct invariants for 2D flow
scenarios is three, not five—and empirically through improved predictions. It
has been shown that the new approach yields more accurate predictions of the
anisotropy tensor than the original TBNN of Ling et al. Ling et al. [2016].
The improvements were illustrated with the barycentric colormap and quantified
by comparing RMSEs. It is, however, limited in its predictive capabilities of
the mean velocity and pressure fields. While it still outperformed the widely
popular $k-\epsilon$ model on geometries with flow separation, it consistently
fails to compete with the $k-\omega$ model. The rather large discrepancy of
the RANS using the DNS anisotropy tensor indicates that it is more beneficial
to train models that not only aim at improving predictions of the Reynolds
stresses but instead target the mean field quantities directly, e.g., as in
Agrawal and Koutsourelakis [2023], Hayek et al. [2018], Parish and Duraisamy
[2016].

##  5 Broader Impact

Turbulence is a key physical characteristic of a broad range of fluid flows.
Understanding this phenomenon is crucial for complex designs, environmental
modeling, and many more engineering applications. Computational power has
increased massively in the past decades, enabling scale-resolving simulations
like direct numerical simulations of a number of canonical turbulent flows.
However, fast approximations like the RANS continue to remain essential for
industrial applications, whose accuracy hinges heavily on turbulence closure
models.

The presented research serves as an extension to the state-of-the-art data-
driven turbulence closure model proposed by Ling et al. [2016]. By combining a
deep neural network with an inductive bias informed by turbulence
realizability constraints, plus an extensive feature set, the PI-TBNN showed
considerable improvements. These improvements showcased through barycentric
colormap visualizations and quantified reductions in RMSE signify a step
forward in our ability to capture complex turbulence states, particularly in
scenarios involving surface curvature and flow separation.

However, the study also sheds light on the nuanced relationship between
improved anisotropy tensor predictions and the ultimate goal of predicting
mean velocity and pressure fields. While the PI-TBNN excels in enhancing
anisotropy tensor predictions, it does not consistently outperform the
established $k-\omega$ model in mean-field predictions, underscoring the need
for further exploration in this area.

We do not see any direct ethical concerns associated with this research. The
impact on society is primarily through the over-arching context of research
using machine learning to improve our general understanding of turbulence in
fluids.

## References

  * Agrawal and Koutsourelakis [2023] A. Agrawal and P.-S. Koutsourelakis.  A probabilistic, data-driven closure model for RANS simulations with aleatoric, model uncertainty.  _arXiv preprint_ , 7 2023.  doi: https://doi.org/10.48550/arXiv.2307.02432.  URL http://arxiv.org/abs/2307.02432.
  * Banerjee et al. [2007] S. Banerjee, R. Krahl, F. Durst, and C. Zenger.  Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches.  _Journal of Turbulence_ , 8:1–27, 2007.  ISSN 14685248.  doi: 10.1080/14685240701506896.  URL https://www-tandfonline-com.eaccess.ub.tum.de/doi/abs/10.1080/14685240701506896.
  * Bentaleb et al. [2012] Y. Bentaleb, S. Lardeau, and M. A. Leschziner.  Large-eddy simulation of turbulent boundary layer separation from a rounded step.  _Journal of Turbulence_ , 13:1–28, 1 2012.  ISSN 14685248.  doi: 10.1080/14685248.2011.637923.  URL https://www.tandfonline.com/doi/full/10.1080/14685248.2011.637923.
  * Duraisamy [2021] K. Duraisamy.  Perspectives on machine learning-augmented Reynolds-averaged and large eddy simulation models of turbulence.  _Physical Review Fluids_ , 6(5):050504, 5 2021\.  ISSN 2469990X.  doi: 10.1103/PHYSREVFLUIDS.6.050504/FIGURES/5/MEDIUM.  URL https://journals.aps.org/prfluids/abstract/10.1103/PhysRevFluids.6.050504.
  * Duraisamy et al. [2019] K. Duraisamy, G. Iaccarino, and H. Xiao.  Turbulence Modeling in the Age of Data.  _Annual Review of Fluid Mechanics_ , 51(1):357–377, 2019.  ISSN 0066-4189.  doi: 10.1146/annurev-fluid-010518-040547.
  * Hayek et al. [2018] M. E. Hayek, Q. Wang, and G. M. Laskowski.  Adjoint-based optimization of RANS eddy viscosity model for U-bend channel flow.  In _AIAA Aerospace Sciences Meeting, 2018_ , number 210059, Reston, Virginia, 1 2018. American Institute of Aeronautics and Astronautics Inc, AIAA.  ISBN 9781624105241.  doi: 10.2514/6.2018-2091.  URL https://arc.aiaa.org/doi/10.2514/6.2018-2091.
  * Launder and Sharma [1974] B. E. Launder and B. I. Sharma.  Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc.  _Letters in Heat and Mass Transfer_ , 1(2):131–137, 11 1974.  ISSN 00944548.  doi: 10.1016/0094-4548(74)90150-7.
  * Laval and Marquillie [2011] J. P. Laval and M. Marquillie.  Direct numerical simulations of converging–diverging channel flow.  In _ERCOFTAC Series_ , volume 14, pages 203–209. Springer Netherland, 2011.  ISBN 9789048196029.  doi: 10.1007/978-90-481-9603-6{\\_}21.
  * Leon Riccius [2021] Leon Riccius.  Machine Learning Augmented Turbulence Modelling for the Reynolds Stress Closure Problem, 3 2021.
  * Ling et al. [2016] J. Ling, A. Kurzawski, and J. Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  _Journal of Fluid Mechanics_ , 807:155–166, 2016.  ISSN 14697645.  doi: 10.1017/jfm.2016.615.
  * Mellen et al. [2000] C. P. Mellen, J. Frohlic, and W. Rodi.  Large Eddy Simulation of the flow over periodic hills.  _16th IMACS World Congress_ , pages 3–8, 2000.
  * Parish and Duraisamy [2016] E. J. Parish and K. Duraisamy.  A paradigm for data-driven predictive modeling using field inversion and machine learning.  _Journal of Computational Physics_ , 305:758–774, 2016\.  ISSN 10902716.  doi: 10.1016/j.jcp.2015.11.012.  URL http://www.elsevier.com/open-access/userlicense/1.0/.
  * Paszke et al. [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  PyTorch: An imperative style, high-performance deep learning library, 12 2019.  ISSN 23318422.  URL http://arxiv.org/abs/1912.01703.
  * Pinelli et al. [2010] A. Pinelli, M. Uhlmann, A. Sekimoto, and G. Kawahara.  Reynolds number dependence of mean flow structure in square duct turbulence.  _Journal of Fluid Mechanics_ , 644:107–122, 2010.  ISSN 00221120.  doi: 10.1017/S0022112009992242.  URL http://openaccess.city.ac.uk/.
  * Pope [1975] S. B. Pope.  A more general effective-viscosity hypothesis.  _Journal of Fluid Mechanics_ , 72(2):331–340, 11 1975.  ISSN 14697645.  doi: 10.1017/S0022112075003382.  URL http://www.journals.cambridge.org/abstract_S0022112075003382.
  * Schumann [1977] U. Schumann.  Realizability of Reynolds-stress turbulence models.  _Physics of Fluids_ , 20(5):721–725, 1977.  ISSN 10706631.  doi: 10.1063/1.861942.  URL https://doi.org/10.1063/1.861942.
  * Taghizadeh et al. [2020] S. Taghizadeh, F. D. Witherden, and S. S. Girimaji.  Turbulence closure modeling with data-driven techniques: physical compatibility and consistency considerations.  _New Journal of Physics_ , 22(9):093023, 9 2020\.  ISSN 1367-2630.  doi: 10.1088/1367-2630/ABADB3.  URL https://iopscience.iop.org/article/10.1088/1367-2630/abadb3https://iopscience.iop.org/article/10.1088/1367-2630/abadb3/meta.
  * Wang and Xiao [2016] J. X. Wang and H. Xiao.  Data-driven CFD modeling of turbulent flows through complex structures.  _International Journal of Heat and Fluid Flow_ , 62:138–149, 2016.  ISSN 0142727X.  doi: 10.1016/j.ijheatfluidflow.2016.11.007.
  * Wang et al. [2017a] J. X. Wang, J. Wu, J. Ling, G. Iaccarino, and H. Xiao.  A comprehensive physics-informed machine learning framework for predictive turbulence modeling, 2017a.
  * Wang et al. [2017b] J. X. Wang, J. L. Wu, and H. Xiao.  Physics-informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data.  _Physical Review Fluids_ , 2(3):34603, 2017b.  ISSN 2469990X.  doi: 10.1103/PhysRevFluids.2.034603.
  * Weller et al. [1998] H. G. Weller, G. Tabor, H. Jasak, and C. Fureby.  A tensorial approach to computational continuum mechanics using object-oriented techniques.  _Computers in Physics_ , 12(6):620, 12 1998.  ISSN 08941866.  doi: 10.1063/1.168744.  URL http://scitation.aip.org/content/aip/journal/cip/12/6/10.1063/1.168744.
  * Wilcox [2008] D. C. Wilcox.  Formulation of the k-$\omega$ turbulence model revisited.  In _AIAA Journal_ , volume 46, pages 2823–2838, 11 2008.  doi: 10.2514/1.36541.  URL https://arc.aiaa.org/doi/abs/10.2514/1.36541.
  * Wu et al. [2018] J. L. Wu, H. Xiao, and E. Paterson.  Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework.  _Physical Review Fluids_ , 7(3):74602, 2018.  ISSN 2469990X.  doi: 10.1103/PhysRevFluids.3.074602.

##  6 Supplementary Material

###  6.1 RGB colormap

As demonstrated in Figure 1, each point in the barycentric triangle
corresponds to a unique color. The mapping from the barycentric coordinates to
the RGB values follows

$\displaystyle\begin{bmatrix}\mathrm{R}\\\ \mathrm{G}\\\
\mathrm{B}\end{bmatrix}=\frac{1}{\max C_{ic}}\left(C_{1c}\begin{bmatrix}1\\\
0\\\ 0\end{bmatrix}+C_{2c}\begin{bmatrix}0\\\ 1\\\
0\end{bmatrix}+C_{3c}\begin{bmatrix}0\\\ 0\\\
1\end{bmatrix}\right)\quad\mathrm{for}\;i\in\\{1,2,3\\}.$ (7)

###  6.2 Scalar invariants

Two of the five scalar invariants ($\lambda_{3},\lambda_{4}$) are zero for
flows with one direction of homogeneity. The two invariants read

$\displaystyle\lambda_{3}=\mathrm{tr}(\boldsymbol{\hat{S}}^{3}),\hskip
85.35826pt\lambda_{4}=\mathrm{tr}(\boldsymbol{\hat{\Omega}}^{2}\boldsymbol{\hat{S}}).$
(8)

Assuming the flow is homogeneous in $z$-direction, the partial derivatives of
the mean velocity with respect to $z$ vanish, and the mean rate of strain and
rotation read

$\displaystyle\hat{S}_{ij}=\frac{1}{2}\frac{k}{\epsilon}\begin{bmatrix}2\frac{\partial\left\langle
U_{x}\right\rangle}{\partial x}&\frac{\partial\left\langle
U_{x}\right\rangle}{\partial y}+\frac{\partial\left\langle
U_{y}\right\rangle}{\partial x}\\\ \frac{\partial\left\langle
U_{y}\right\rangle}{\partial x}+\frac{\partial\left\langle
U_{x}\right\rangle}{\partial y}&2\frac{\partial\left\langle
U_{y}\right\rangle}{\partial y}\\\ \end{bmatrix},\hskip
7.11317pt\hat{\Omega}_{ij}=\frac{1}{2}\frac{k}{\epsilon}\begin{bmatrix}0&\frac{\partial\left\langle
U_{x}\right\rangle}{\partial y}-\frac{\partial\left\langle
U_{y}\right\rangle}{\partial x}\\\ \frac{\partial\left\langle
U_{y}\right\rangle}{\partial x}-\frac{\partial\left\langle
U_{x}\right\rangle}{\partial y}&0\\\ \end{bmatrix}.$ (9)

The incompressibility constraint of the Reynolds equations reduces to

$\displaystyle\frac{\partial\left\langle U_{i}\right\rangle}{\partial
x_{i}}=\mathrm{tr}\left(\frac{\partial\left\langle
U_{i}\right\rangle}{\partial x_{j}}\right)=\frac{\partial\left\langle
U_{x}\right\rangle}{\partial x}+\frac{\partial\left\langle
U_{y}\right\rangle}{\partial y}=0.$ (10)

When using the simplified expressions of $\boldsymbol{\hat{S}}$ and
$\boldsymbol{\hat{\Omega}}$ in combination with the incompressibility
constraint, invariant $\lambda_{3}$ is given by

$\displaystyle\mathrm{tr}(\hat{S}_{ij}^{3})$ $\displaystyle=$ (11)
$\displaystyle\footnotesize{\frac{k^{3}}{8\epsilon^{3}}}\mathrm{tr}$
$\displaystyle\left({\footnotesize\begin{bmatrix}\hat{S}_{11}(\hat{S}_{11}^{2}+\hat{S}_{12}^{2})+\hat{S}_{12}(\hat{S}_{11}\hat{S}_{12}+\hat{S}_{12}\hat{S}_{22})&\hat{S}_{12}(\hat{S}_{11}^{2}+\hat{S}_{12}^{2})+\hat{S}_{22}(\hat{S}_{11}\hat{S}_{12}+\hat{S}_{12}\hat{S}_{22})\\\
\hat{S}_{11}(\hat{S}_{11}\hat{S}_{12}+\hat{S}_{12}\hat{S}_{22})+\hat{S}_{12}(\hat{S}_{12}^{2}+\hat{S}_{22}^{2})&\hat{S}_{22}(\hat{S}_{12}^{2}+\hat{S}_{22}^{2})+\hat{S}_{12}(\hat{S}_{11}\hat{S}_{12}+\hat{S}_{12}\hat{S}_{22})\\\
\end{bmatrix}}\right)$
$\displaystyle=\frac{k^{3}}{8\epsilon^{3}}\left(\hat{S}_{11}(\hat{S}_{11}^{2}+\hat{S}_{12}^{2})+2\hat{S}_{12}^{2}\underbrace{(\hat{S}_{11}+\hat{S}_{22})}_{=0}+\hat{S}_{22}(\hat{S}_{22}^{2}+\hat{S}_{12}^{2})\right),$
(12)
$\displaystyle=\frac{k^{3}}{8\epsilon^{3}}\left(\hat{S}_{11}(\hat{S}_{11}^{2}+\hat{S}_{12}^{2})+\hat{S}_{22}(\hat{S}_{11}^{2}+\hat{S}_{12}^{2})\right),$
(13)
$\displaystyle=\frac{k^{3}}{8\epsilon^{3}}\left(\underbrace{(\hat{S}_{11}+\hat{S}_{22})}_{=0}(\hat{S}_{11}^{2}+\hat{S}_{12}^{2})\right)=0.$
(14)

The derivation of invariant $\lambda_{4}$ is more straightforward and thus
written in terms of the mean velocity gradient. It is given by

$\displaystyle\mathrm{tr}(\hat{\Omega}_{ij}^{2}\hat{S}_{jk})$
$\displaystyle=\mathrm{tr}\left(\frac{k^{3}}{8\epsilon^{3}}\left(\frac{\partial\left\langle
U_{x}\right\rangle}{\partial y}-\frac{\partial\left\langle
U_{y}\right\rangle}{\partial
x}\right)^{2}\begin{bmatrix}2\frac{\partial\left\langle
U_{x}\right\rangle}{\partial x}&\frac{\partial\left\langle
U_{x}\right\rangle}{\partial y}+\frac{\partial\left\langle
U_{y}\right\rangle}{\partial x}&0\\\ \frac{\partial\left\langle
U_{y}\right\rangle}{\partial x}+\frac{\partial\left\langle
U_{x}\right\rangle}{\partial y}&2\frac{\partial\left\langle
U_{y}\right\rangle}{\partial y}&0\\\ 0&0&0\\\ \end{bmatrix}\right)$ (15)
$\displaystyle=\frac{k^{3}}{4\epsilon^{3}}\left(\frac{\partial\left\langle
U_{x}\right\rangle}{\partial y}-\frac{\partial\left\langle
U_{y}\right\rangle}{\partial
x}\right)\underbrace{\left(\frac{\partial\left\langle
U_{x}\right\rangle}{\partial x}+\frac{\partial\left\langle
U_{y}\right\rangle}{\partial y}\right)}_{=0}=0.$ (16)

###  6.3 Data set

The high-fidelity direct numerical simulation (DNS)/large eddy simulation
(LES) data are used for the training. Out of the two available DNS for the
CDC, the one at ${\mathrm{Re}=7900}$ was used for estimating the
regularization parameters. The higher Reynolds number simulation at
${\mathrm{Re}=12600}$ was used for training and validation, along with the DNS
at ${\mathrm{Re}=2800}$ and LES at ${\mathrm{Re}=10595}$ for PH and the DNS
for SD at $\mathrm{Re}=\\{2000,2400,2900,3200\\}$. The data set was split into
training and validation sets at a ratio of $70/30$. Therefore, the total
number of data points available for training was $26600$.

The testing set consists of three flow geometries (SD, PH, and CBFS). Two of
these geometries, SD and PH, are also part of the training and validation set,
however, at different Reynolds numbers. The periodic hills case at
${\mathrm{Re}=5600}$ was selected to investigate the interpolation properties
— the ML model has previously seen this flow geometry and is expected to yield
good predictions. The square duct is a canonical flow case that clearly
illustrates the deficiencies of the LEVMs and is suitable to present
difficulties of propagating the Reynolds stresses to the flow field. Finally,
the curved backward-facing step tests the ML model’s extrapolation
capabilities.