11institutetext:  Geneva Observatory, Geneva University, CH-1290 Versoix,
Switzerland

# Impact of Various Accretion Rates on Population III Stars

Devesh Nandal et al

###### Abstract

Context.

Aims.  We compute models of accreting population III protostars to explore the
impact of different accretion rates on their structure and evolution. We
explore accretion rates from 10-6 to 10 and study the impact on the maximum
mass reached, reionization and finale fates of the resulting populations. We
use the stellar evolution code genec to compute models of the evolution of
accreting population III protostars.

study the formation of massive and supermassive stars in dark matter mini-
halos exploring the impact of different accretion rates.

We use the stellar evolution code genec to compute the evolution of accreting
protostars made of primordial compositions (i.e. resulting only from the Big-
Bang nucleosynthesis). We explore a range of accretion rates spanning seven
orders of magnitudes and study the impact on the maximum mass reached,
reionization and finale fates of the resulting populations.

We estimate in which conditions the infalling material might be optically
thick and thus affect the apparent effective temperature of these objects. We
also study some physical properties of the supermassive stars that can form in
atomic cooling halos.

We show that in case of accretion rates below 2.5x10-5 M⊙ per year, the
maximum mass reached would be below the minimum mass for producing pair
instability supernovae. In order to produce black holes with masses above 250
M⊙, the mininum accretion rate is 5x10-5 M⊙ per year. For accretion rates
below 10-2 M⊙ per year, the maximum mass that can be reached is mainly
determined by the mass at which the formation timescale becomes equal to the
core- Hydrogen burning timescale. For accretion rates above 10-2 M⊙ per year,
the maximum mass reached is due to reaching the general relativity instability
during the core H-burning phase. Supermassive stars, that reaches the GR
instability when they have masses around 100 000 M⊙, would appear as objects
with a luminosity around NN and an effective temeprature of NN. They may
present pulsational instability of… The impact on the reionization of normal
Pop III stars of the accretion rates is rather modest and would only affect
populations coming from accretion rates below about NN M⊙ per year. The
precursor of the supermassive stars that are the progenitors of the seed for
forming massive black holes at high redshift have maximum luminoisties around
1010 L⊙. They are pulsating objects with period of about 100 years and would
appear as long wavelength emitters.

Methods.

Results. Conclusions. binaries:close-stars; stars: abundances;
rotation;evolution

##  1 Introduction

The scarcity of massive stars in the current universe limits the study of
their evolutionary stages. In particular, the formation phase is quite
perplexing since the central protostellar object is directly concealed
throughout the pre-main sequence evolution (REF). Observations and theoretical
modeling provide valuable insight towards the stellar formation process (REF).
The process commences with an initial dynamical collapse of a molecular cloud
leading to its eventual fragmentation (REF Rees’ Paper). In the regions of
over-densities within the fragmented clouds, a hydrostatic stellar core is
formed and grows until the surrounding reservoir of matter is exhausted. The
growth of core may also be halted by various feedback mechanisms that repel
the infalling matter (REF - Lyman Werner UV Feedback, B field etc). This marks
the end of accretion phase and the central object emerges from the optically
thick media attaining its final mass. Depending on the timescale of accretion,
the protostar may undergo further contraction or might already have a central
temperature sufficient enough for hydrogen ignition and to disembark on its
main-sequence evolution (REF). Stellar mergers is another possible formation
scenario for massive stars (REF) but this study explores the accretion
scenario as it is generally agreed to be the most common formation channel
(REF).

Accretion rate is the underlying quantity that impacts the pre-main sequence
evolution and it scales with the free fall timescale of the cloud or
equivalently the sound velocity of the medium. Analytical and numerical
experiments hint towards the accretion rate being non-constant and varying as
a function of current mass or time and described by accretion laws ranging
from solar (REF) down to primordial metallicities (REF).

Accreting models are often computed using a pre-existing core with mass around
NN M⊙ at hydrostatic equilibrium. Although deuterium burning occurs within the
core, it is energetically insignificant and the main source of energy
production comes from the gravitational contraction of the object. Accretion
onto this object may occur via spherical infall (REF) or through a cold
accretion disc where the entropy on accreting mass is equivalent to the
entropy of the surface layers of the protostar (REF). The choice of accretion
geometry is important as the initial structure of the core has an impact on
the further stages of evolution (REF).

The pre-Main-Sequence track of a star is composed of two parts. The first part
is the accretion phase that follows the birthline and its timescale is given
by the formation timescale, $\tau_{\rm form}=M/\dot{M}_{\rm acc}$, with
$\dot{M}_{\rm acc}$ the accretion rate (here we assume a constant star
formation rate or a time-averaged value). The second phase is the contraction
phase towards the zero-age-main-sequence and the timescale is given by
$\tau_{\rm KH}$. For cases with $\tau_{\rm form}<\tau_{\rm KH}$, the pre-MS
phase is composed from a short accretion phase followed by a long contraction
phase; here the birtline will join the ZAMS in the HR diagram where $\tau_{\rm
form}\sim\tau_{\rm KH}$. For case with $\tau_{\rm form}>\tau_{\rm KH}$, the
contraction phase actually occurs during the accretion phase and therefore the
star can undergo core H-burning while it is still accreting mass (REF) in the
natal cloud. This has an impact on the star’s ionising power and
interestingly, it also opens up a way to constraint an upper mass limit of
stars. A star cannot be more massive than the limit given by equating its
formation timescale to its lifetime. Of course, the formation can be
interrupted by some feedback mechanisms (radiative feedback or centrifugal
feedback), however it has been argued that such limits can be routed (REF). To
obtain such limits, it is necessary to compute the birthlines of protostellar
objects, i.e. accreting tracks until the end of the evolution.

In the present work we investigate a few questions along the accretion
scenario in primordial environements. These questions are

  * •

What are the differences between the accretion scenario for the formation of
Pop III stars and that of solar composition stars? In the literature these two
questions are not often addressed together with the same stellar evolution
code and we take the opportunity to fill this gap here.

  * •

What is the impact of the formation process on the ionization rates of Pop III
stars?

  * •

What is the upper mass limit of stars of Pop III stars obtained for different
accretion rates?

  * •

How the final fate of Pop III may depend on the accretion rates during their
formation process?

In section 2,…

##  2 Stellar Models

###  2.1 Ingredients of the stellar models

We use the genec stellar evolution code (Eggenberger2008) to compute our
models. The population III models are constructed with a composition of X =
0.xx and Y = 0.xx, the same as in in Murphy2021a. The pre-Main-Sequence phase
and the core hydrogen burning phases are computed with accretion. Other
physical inputs including the equation of state, nuclear reaction rates and
opacities are the same as in Ekstrom2012. When a convective core is present,
we have used a step overshoot equal to xx. We do not consider the impact of
stellar rotation in this work.

The accretion is accounted for as in Haemmerle2014; Haemmerle2016a. Accretion
mode assumes the infall of matter via a geometrically thin cold disc which
implies that the specific entropy of the accreted matter is the same as that
of the stellar surface. Furthermore, any excessive buildup of entropy during
accretion is assumed to be radiated away adequately via poles
(Palla1992);(Hosokawa2013). A point of note here is regarding the treatment of
outer regions of star in GENEC. The outer region is defined as the region
where ionization is partial and convection is non-adiabatic, while interior is
assumed to be completely ionized and convection adiabatic. The border between
these two regions occurs at a fraction of the total mass that is around
0.999111This limits may change along the evolution. Typically along the
Hayashi track, because the envelope is cooler, complete ionization is reached
in deeper layer and a lower value of FITM may be chosen. When the star is in
the blue region of the HRD, smaller values of FITM can be chosen. The effects
of varying fitm are discussed in (Haemmerle2018).

Computation of the models begin by accreting primordial matter on hydrostatic
seeds at high rates. Composition of matter infalling is identical to the
hydrostatic seed and consists of a homogeneous distribution of hydrogen with a
mass fraction $X=0.7516$, of helium with $Y=0.2484$ and deuterium at
$X(D)=5x10^{-5}$ (Bernasconi1996; Behrend2001; Haemmerle2016a). For models
with accretion rate of $<0.01M_{\odot}/yr$ the mass of hydrostatic seed is
0.7M⊙, for rates 0.001, 0.01 and 0.1 M⊙/yr, the mass of the hydrostatic seed
is 2 M⊙ but for the higher accretion rates of 1 and 10 M⊙/yr, a 10 M⊙ seed is
needed to facilitate numerical convergence. The models initially begin their
lifetimes as fully convective objects and posses a flat entropy profile with
central temperatures at $6.6x10^{5}$ for 2 M⊙ seeds and $4.1x10^{5}$ for 10 M⊙
seeds.

The computation for all models is followed until the end of Hydrogen burning
as it accounts for nearly 90 percent of the evolution. Further evolution
becomes computationally challenging, especially for higher accretion rates
($>0.01M_{\odot}/yr$) and therefore will be pursued in future work.

Table 1: A description of stellar parameters for all models. The accretion rate is given by $\dot{M}[{M_{\odot}}/yr]$, Z is the metallicity of the model, $t_{\rm H}$ is the lifetime of core hydrogen burning in million years, $t_{red}/t_{total}$ is the ratio of time spent in red versus the total lifetime of star during main sequence and M${}_{end}[{M_{\odot}}]$ is the mass reached either at at the end of hydrogen burning; $*$ indicates the mass reached at the end of computation which is before the completion of H burning.  _Model_ | _$\dot{M}[{M_{\odot}}/yr]$_ | _Z_ | _$t_{\rm H}[Myr]$_ | _$t_{red}/t_{total}$_ | _$M_{end}[{M_{\odot}}]$_
---|---|---|---|---|---
yellow!40A | $10^{-6}$ | 0.000 | 1.040 | 2.261 | 2.796
brown!30B | $10^{-5}$ | 0.000 | 0.630 | 0.998 | 0.003
red!30C | $10^{-4}$ | 0.000 | 104046 | - | -
blue!25 D | $10^{-3}$ | 0.000 | 2.796 | 0.003 | 2797
green!30 E | $10^{-2}$ | 0.000 | 2.261 | 0.998 | 22616
cyan!30 F | $10^{-1}$ | 0.000 | 1.040 | 0.630 | 104046
orange!30 G | $1$ | 0.000 | 0.214 | 0.230 | 214092∗
purple!30 H | $10$ | 0.000 | 0.045 | 0.438 | 449661∗
yellow!60 I | $CH$ | 0.000 | 67327 | 9912 | 1471
brown!70 J | $Yoshida$ | 0.000 | 67327 | 9912 | 1471
magenta!30 K | $10^{-3}$ | 0.014 | 67327 | 9912 | 1471
green!60 L | $10^{-5}$ | 0.014 | 67327 | 9912 | 1471

###  2.2 The accretion rates

The accretion phase is driven by complex processes involving in-falls and
outflows. This is is likely true at all metallicities ranging from Pop III to
solar.

####  2.2.1 Primordial versus solar metallicities

A significant difference however between the case of primordial star formation
and those of slightly enriched ones is the averaged temperature in the forming
cloud. While at solar metallicity the typical temperature is 10 K, in
primordial regions it is around 1000 K. This changes the free fall timescale
and thus the accretion rate. We expect higher accretion rates at low $Z$ than
in solar metallicity regions.

In the present work, we want to investigate different accretion rates on the
pre-MS and MS evolution of massive Pop III stars. We shall start by
investigating accretion rates that are predicted by hydrodynamical simulations
of Pop III star formation and discuss, from a stellar physics point of view
the similarities as well as the differences with the pre-Main-Sequence (PMS)
phase at solar metallicity. In a next step we shall explore the impacts of
smaller and larger accretion rates discussing impact on the most massive stars
that can be formed, on the path of the accreting object in the HR diagram, its
ionizing power and on the nature of its final fate.

In case of the lower accretion rates, the results obtained might be a bit
academic in the sense that they are not supported by hydrodynamical
simulations. On the other hand, hydrodynamical simulations cannot resolve the
stars themselves. They provide hints only on mass fluxes through surfaces that
are much larger than the stars themselves. If the mass fluxes over those
surface is likely very reliable, it is actually not guarantee that all that
mass will finish locked into stars, so that the net actual accretion rates may
be smaller than the one given by these simulations.

Multi-D hydrodynamical simulations in general stops when only a protostellar
cloud is formed with typically a density expressed in term of $n_{\rm H}$, the
number of hydrogen atoms per cubic centimeter of the order of $10^{21}$. These
simulations cannot for instance model the disk that likely is formed, neither
jets that at least accompany star formation at solar metallicity, thus we can
legitimately wonder whether the whole mass flux predicted at much larger
distance will actually be the accretion rate on the young star? We know also
that the reality will likely host a diversity of situations. Thus it appears
legitimate to broaden the domain of accretion rates and study the potential
impacts on the ouputs of the models of star formation.

Thus, in the present paper we shall follow a kind of agnostic view allowing us
to explore much smaller and much larger accretion rates than classical ones.
This will allow to put in evidence a few interesting consequences coming
either from much lower or much larger accretion rates than classically
considered.

In the domain of the low accretion rate, we shall see that it actually may
impact the ionizing power of these first generations of stars. In the domain
of the very large accretion rates, supermassive stars are formed whose maximum
values can either be limited by their formation timescale that becomes longer
than. their lifetimes or by reaching conditions favoring the star to become
unstable (GR instablity). We shall discuss a few properties of these
supermassive stars. Other aspects of the physics of these supermassive stars
will be discussed in more details in an accompanyyng paper.

  * •

Range of accretion rates that you consider and why? Explain Yoshida and CH
rates

  * •

evolutionary state

  * •

assumptions / caveats

##  3 Formation process at different metallicities

Figure 1: Birthlines for different accretion rates and metalicities. The
dashed red lines correspond to the variable accretion rate obtained from
hydrodynamical simulations by Yoshida et al. 2006 for Pop III stars (see
text). The solid black track represents the model at solar metallicity with a
constant accretion rate of 10$-5$ M⊙ per year. The short dashed green track
corresponds to the zero metallicity model with constant accretion rate of
10$-5$ M⊙ per year. The black circle on each track marks the beginning of core
hydrogen burning and the black triangle marks the end of H burning. Masses at
the end of hydrogen burning are indicated in solar masses at the end of each
track. The solid blue numbered dots form the ZAMS of Murphy et al. (2020) at
zero metallicity and the solid yellow dots are taken from solar metallicity
models of Ekström et al. (2012). The positions of a few initial masses are
indicated along the ZAMS. The grey dotted lines are the iso-radii of the
indicated radii.  Figure 2:  I would swap the x and y axes The variation of
accretion rate with respect to luminosity for the Yoshida model at zero
metallicity throughout its entire evolution. The red and green markers
represent the mass of the model at the corresponding accretion rate and shows
that most of the mass accreted was during the time when accretion rate was
lower.  Figure 3: Plot depicting the ratio of Kelvin Helmholtz timescale to
accretion timescale versus the mass of models at zero metallicity. The red
line corresponds to the hydrodynamical accretion law obtained from Yoshida et
al. 2006 which is compared with the blue track representing a
$10^{-3}$$M_{\odot}$/yr and green track set at $10^{-5}$$M_{\odot}$/yr
accretion rate

In Figure 1, we compare typical birthlines at solar and primordial
composition. The one for primordial matter has been obtained with accretion
rates deduced from hydrodynamical simulations from Yoshida et al REF. It is
given by the expression below

$\dot{M_{\rm acc}=}$

The one for solar metallicity has been obtained for a constant accretion rate
of 10-5 M⊙ per year. Actually in REF where we studied more specifically the
pre-MS phase at solar metallicity, we used a more sophisticated accretion
rates, the so called Churchwll-Henning accretion rate). However the case of
constant accretion rate around 10-5 M⊙ per year is not very far from the case
obtained with the Churchwell-Henning accretion rate.

Comparing the red and black birthlines provide some comparisons between what
is expected to happen in primordial mini-haloes and what is expected to occur
in the solar neighborhood about the formation of stars. We can note the
following differences. The stronger accretion rate in primordial conditions
with respect to solar conditions produces two striking differences. First the
track for the accreting Pop III protostar is shifted to lower effective
temperatures with respect to the solar metallicity one during the first phase
of the evolution before the star contracts to the ZAMS. This is an effect
mainly due to the different accretion rate and not so much to the difference
metallicity. This can indeed be seen on the right panel of 1, where birthlines
obtained with a constant accretion rate of 10-3 M⊙ per year are plotted for a
Pop III (black lines) and a solar metallicity track (red line, IT WOULD BE
BETTER USE SAME COLORS AS IN LEFT PANEL). We see that the two tracks are
superposed, indicating that the change in metallicity has actually little
effect). When the accretion rate increases, the accreting star has longer
thermal adjustment timescale. The energy released by the internal contraction
cannot get out rapidly and part remains locked in the star keeping it more
inflated.

Looking at the left panel of Fig. 1, we see that the birthline at solar
metallicity and thus with the smaller accretion rate joins the ZAMS around a
mass of 7 M⊙, while this happens for a much higher mass in the Z=0 higher
accretion rate case (around 120 M⊙). This is again an effect due to the change
in the accretion rate that has already been discussed in previous works (REF).
Higher the accretion rate is, larger the mass reached during a Kelvin-
Helmholtz timescale. Of course we can wonder at which epoch the Kelvin-
Helmholtz timescale has to be estimated. In Fig. LABEL:timescales, we compare
accretion timescales with MS timescales as given by the models of Murphy et
al. and with the Kelvin-Helmholtz timescale estimated at the mass when the
accretion stops. The orange line (KH time) intersects the red line
($10^{-3}$$M_{\odot}$/yr) at the mass of 39$M_{\odot}$. This mass should
correspond to the mass at which the birthline in black in the right panel of
Fig. 1 joins the ZAMS. We see that the true mass is significantly larger, more
than 60$M_{\odot}$. Thus we see that using the KH timescale of the final mass
underestimates the mass at which the birthline joins the ZAMS. This is indeed
expected since the KH at the final mass corresponds to the contraction time of
a constant mass, while here we have a growing mass. Thus in previous phases,
the KH is larger tha the one one obtains for a constant 60$M_{\odot}$mass
model. Numerical computations are needed to infer this value. We see also
looking at the right panel of Fig. 1, that the mass at which the birthline
joins the ZAMS depends also on the metallicity. Indeed, the right panel of
Fig. 1 indicates that at solar metallicity this mass is around 32, while it is
60 at Z=0 metallicity. This comes from the fact that star at a metallicity Z=0
are more compact. They are more compact due to the fact that stars have lower
opacities but also because they lack CNO elements as catalyst of the H.burning
reaction in the core (see more on this point below). Smaller radii, keeping
the mass and luminosity equal, make larger KH timescales and thus more time
for the star to accumulate mass.

We can note also that along the high accretion rate birthline, a small hook
occurs at a luminosity around log L/L$\odot$ equal to 5, when the mass reaches
a value around 13.3. This occurs when the accretion rate decreases (see the
left panel of Fig. 2) and the ratio between the KH and the accretion timescale
decreases (see the right panel of Fig. 2). Interestingly, would the high
accretion rate corresponding to nearly 0.01 M⊙ per year remains active, likely
the track would evolve to the red part of the HR diagram. The decrease of this
accretion rate allows the accreting protostar to contract onto the ZAMS.

In Fig. LABEL:kip1 and LABEL:kip2, the evolution of the structure can be seen
as a function of time for different metallicities, accretion rates. We see
that changing just the metallicity, keeping the accretion rates the same
produces very similar results in term of the structure evolution. We note
however a few small differences (compare the panels on the left with the one
at the middle in Figs. LABEL:kip1 and LABEL:kip2. The beginning (CAN WE ZOOM
ON THE BEGINNING TO SEE BETTER WHAT DOES HAPPEN THERE?), Are Pop III more
compact at the stage? convective core appears at later time in Pop III (pp-
CNO). Energy generation in H-burning core is more concentrated in radius.

If now we compare the middle panel with the right panel, keeping the
metallicity the same but decreasing the accretion rate, we see that the
protostar gets a smaller radius. At a given mass how do the convective core
radius compares?

DISCUSS TRANSITION pp-CNO, Deuterium burning? SURFACE ABUNDANCES?

##  4 Impact of the accretion law on the ionizing power of Pop III stars

At solar metallicity, when the accretion rate is of order of 10-5 M⊙ per year,
we expect that this would produce a lack of young O-type star near the ZAMS.
Recent observations by REF show that indeed there seems to be a lack of such
young O-type stars. In that work, we have shown that this region devoid of
O-type star might be explained in the frame of the accreting scenario for the
formation of massive stars with a rate of around 10-5 M⊙ per year. This is
interestng since a single value may reproduce both the observed upper envelope
of pre-MS stars in the low luminosity range and the lower envelope in the
upper HR diagram for core H-burning stars. When accretion rates as those
obtained by Yoshida et al. REF, this effect would not actually be present. The
accretion rates are so high that when the star stops accreting (for reasons
that we shall not discuss here), it has consumed very small amounts of
hydrogen in the core and the star would still be near the ZAMS. This of course
would be very different in case the accretion rates would be much smaller. We
shall discuss this point in the next section.

In Fig. xx, we compare birthlines at Z=0.000 amd Z=0.014 for constant and
equal accretion rates. We note when tzhe accretion rates are similar, the
birthlines are very similar except that the ZAMS is shifted to the blue in
case of the Z=0.000 tracks. This shift to bluer position is due to two
effects. First the opacity in Z=0.000 models are lower (COMPARE OPACITIES IN
OUTER LAYERS WHEN THE BIRTHLINE JOIN THE ZAMS IN Z=0.014 AND Z=0.000 MODELS).
The other effect is that in solar metallicity model CNO burning is activated
while in the Pop III models pp chains are the main energy source. The pp
chains are less efficient in producing energy at a given temperature (GIVE MAY
BE AN ORDER OF MAGNITUDE). Thus the Pop III model has to reach higher central
temperature to release equivalent amount of energy by nuclear fusion, this
implies that the star has to contract more.

Figure 4: Different birthlines (see text) at zero metallicity in the HR
diagram obtained with various accretion laws. The tracks in various colours
are represented as follows: solid black line - $10^{0}$$M_{\odot}$/yr, dashed
red line - $10^{-1}$$M_{\odot}$/yr,mini dashed green line -
$10^{-2}$$M_{\odot}$/yr, dotted blue line - $10^{-3}$$M_{\odot}$/yr, mini
dotted and dashed cyan line - $10^{-4}$$M_{\odot}$/yr, dotted and dashed
magenta line - $10^{-5}$$M_{\odot}$/yr, orange line - $10^{-6}$$M_{\odot}$/yr,
olive line - variable accretion rate according to CH law,pink dashed line -
Accretion law from hydrodynamic simulations based on Yoshida et al. . The
solid black dotted line is the ZAMS of Murphy et al. (2021). The coloured
contours begin at the bottom left with 0.01 $R_{\odot}$and progress by a
factor of 10 with every changing shade, ending at 10000 $R_{\odot}$

In Fig. 4, birthlines for accretion rates varying between 10-6 up to 1 M⊙ per
year are presented per step of 1 dex. We also show in Fig. xx, how the mass
reached at the end of the MS phase varies as a function of the accretion rate.

As already mentioned in the introduction, we are aware that these accretions
rates might not all be very relevant for Pop III star formation. On the other
hand, the computations performed allow to have a broad view on how things
would changes when the acretion rate is varied over 6 orders of magnitude and
provide some interesting insights on the physics of star formation by
accretion. Moreover even if for the tracks shown in Fig. 4 we used an initial
composition for zero metallicity stars, we have shown above that at least some
characteristics as the mass at which the birtline joins the ZAMS or the mass
reached at the end of the MS phase does not much depend on metallicity.

Figure 5: Left panel: The abundance of carbon(blue), Nitrogen(green) and
Oxygen(red) versus the mass for a Yoshida et al. 2006 accretion rate model.
Middle panel:The abundance of carbon(blue), Nitrogen(green) and Oxygen(red)
versus the mass for a $10^{-5}$ $M_{\odot}$/yr model.Right panel: The
abundance of carbon(blue), Nitrogen(green) and Oxygen(red) versus the mass for
a 0.1 $M_{\odot}$/yr model.

When the luminoisty increases with time and the effective temperature remain
at much cooler values than the Zero Age Main-Sequence position. This is due to
the fact that the star has not yet reached the radiative equilibrium, namely a
structure such that the energy produced by nuclear burning allows to
compensate for the losses at the surface. A typical timescale equal to the
Kelvin-Helmholtz timescale, $\tau_{\rm KH}=GM^{2}/(RL)$,is needed to achieve
the radiative equilbrium, where $G$ is the gravitational constant, $M$, $R$
and $L$, the mass, radius and luminosity of the star. This timescale decreases
when the mass increases.

The tracks can be devided in two families: the first one comprise the tracks
that after a phase of accretion at large radius contract to the ZAMS
(contracting scheme). The second contain tghe tracks that at a given
luminoisty evolve to the red part of the HR diagram along the Hayashi track
(expanding scheme). What is the physics behind such a dichotomy?

The contracting scheme likely correspond to the case of formation of normal
stars. The masses at which the birthlines join the ZAMS for the different
accretion rates is indiacted in Table 1. As is well known, the mass increases
with the accretion rate (IS IT POSSIBLE O OBTAIN A FIT, AN ANALYTICAL
RELATION?). This point corresponds to the stage at which H-burning occurs at
the centre. Since we are dealling here with Pop III stars, the H-ignition
occurs through pp chains.

When evolution goes on from this stage on, the birthline incurve to the right
due to the fact that H-burning occurs in the core.(WILL CNO CYCLE SWITCH ON
ALREADY DURING THE ACCRETION PHASE?). The accreting core H-burning phase
occurs on still more extended domain of effective temperatures when the
accretion rate increases.

The properies of the high accretion rate tracks that follow an expanding
scheme show higher luminosities than the contracting family tracks in the low
luminous domain. Then they evolve to the right of the HR diagram, so in the
opposite direction of the contracting family tracks. Why and when this happens
depends on the accretion rate.

When the track reaches the Hayashi track, the star evolves vertically until
H-burning reactions begins to occur at the centre. Fro this stage the track
turn towards bluer regions of the HR diagram.

It is interesting to mention that in case accretion stops before the core
H-burning switches on in the core, then the further evolution will bring the
star to the blue (MAKE A TRIAL). In case H-burning is well active in the
central region and that accretion stops we have that star remains in the red
part of trhe HR diagram (to checks).

The effective temperature shown in Fig. 4 corresponds to the effective
temperature at the surface of the hydrostatic core. We can estimate the
optical depeth of the infalling matter in the following way…

Depending on the accretion rates, stars may begin to be unshielded by the
infalling material at different evolutionary stage during their MS phase. In
table 2, we indicate for each accretion rate, the mass at which the star
stopping accreting would appear on the ZAMS and then we indicate the value of
the mass fraction of hydrogen in the core when more massive stars would appear
when they are no longer accreting mass. We indicate also the fraction of the
MS lifetime during which those star would actually emit ionizing photons while
not being surrounded by the accreting material.

Lower is the accretion rate, smaller will be the time duration during which
the new star will be able to emit freely ionizing photons. For bona fide
accretion rate of Pop III star, this effect should however remain modest, but
can become important may be at higher metallicities.

(Murphy2021a)

###  4.1 The light element burning

Figure 6: The plot contains the abundance profile of deuterium in the interior
of the model at different phases in its evolution starting from A and I versus
the internal mass for Yoshida track. The plot shows that for the early models
A,B,C, deuterium is present in the star and can be observed at the surface but
for model D, the amount is almost extinct and becomes 0 for the superseding
models. This corresponds to the temperature of the atmosphere being higher
than 1 million K which immediately destroys deuterium as it is added to the
model by accretion.

Figure 7:  Left panel: The relative abundance of elements versus mass for a
model with accretion rate taken from Yoshida et al. 2006 .The blue curve
depicts Deuterium whereas the red curve is for Lithium. Middle panel: same as
left panel for an accretion rate of $10^{-5}$ $M_{\odot}$/yr Right panel:same
as left panel for an accretion rate of 0.1 $M_{\odot}$/yr

##  5 The upper mass limit of stars

Figure 8: A plot of mass versus accretion rates ranging from
$10^{-6}$$M_{\odot}$/yr to $1000$$M_{\odot}$/yr for zero metallicity models.
An additional solar metallicity model at $10^{-35}$$M_{\odot}$/yr is depicted
in green and red dot. The green shaded region relates corresponds to a
accretion rate that leads to the formation of a black hole with mass up to 150
$M_{\odot}$. The orange region marks the limits between 150-250
$M_{\odot}$where pair instability supernovae could exists and the region in
cyan is the mass domain of intermediate and super massive black holes.

In Fig. xx, the mass reached at the end of the MS phase is indicated as a
function of the accretion rate. The plot indicates also how such star will end
their lifetime. Interestingly, in order to obtain PISNe, the accretion rate
should increase a value of … which is indeed expected tio be reached for Pop
III stars, at leats in our present understanding. We see also that in the
upper mass range, when the accretion rates is above a limit of NN, the star
will collapse due to GR instability before any core H-burning occurs. The
minimum accretion rates allowing the complete core H-burning phase to occur
during the accretion process is…

##  6 Impact of the accretion rate on the final fates of primordial stars

Figure 9: Left panel: Eulerian kippenhahn diagram for a 0.1 $M_{\odot}$/yr
model. Right panel:The variation of Beta (ratio of gas pressure versus total
pressure) versus the mass for a 0.1 $M_{\odot}$/yr model.

Figure 10: Left panel: Log T versus log R plot for the 0.1 $M_{\odot}$/yr
model at the end of Hydrogen burning. The red and green points indicate the
enclosed mass of the interior corresponding to the enclosed radius. A point to
note here is that most of the mass(126000 $M_{\odot}$) is within a radius of
log 13 and only a fraction of mass (still 23000$M_{\odot}$) in in the outer
envelope implying the star is extremely bloated.The cyan region represents a
convective zone whereas yellow yone is the radiative region. Right panel: The
plot depicts the change in abundance of various elements inside the 0.1
$M_{\odot}$/yr model at the end of its Hydrogen burning versus the enclosed
mass.

###  6.1 Optical thickness of the infalling mass

The observational signature of supermassive stars is another aspect that
requires a closer evaluation. To begin this discussion we look at a similar
analogy; if a massive star going through the Wolf-Rayet phase is observed
during helium burning, the central object is not visible to the observer due
to the strong outflow of matter (rate of flow, I guess its $10^{-4}$
$M_{\odot}$/yr )which creates an optically thick medium around the star. If
the situation was reversed and the matter is accreting on a supermassive star,
will the object at the center be visible ? To answer this question, we begin
by comparing two accretion rates of supermassive stars who are placed in an
atomically cooled halo with a radius of 1 parsec. The density of such a medium
is obtained using the relation rho = (), where the velocity is the escape
velocity of mass from the central object. The value for opacity is calculated
at the surface of the star and is assumed to be constant, therefore it only
represents an upper estimate. The optical depth for a $10^{-3}$ $M_{\odot}$/yr
model is calculated at five different mass values and the results indicate
that the medium is optically thick. Furthermore, for an accretion rate of $1$
$M_{\odot}$/yr, the medium becomes even more optically thick which leads to
the conclusion that if accretion time was much larger than the lifetime of H
burning, the supermassive stars will never be visible thoughout their
lifetime.

Figure 11: Left panel: The plot containing density (blue) and Temperature
(red) versus radius of a 1 $M_{\odot}$/yr model at the time of GR instability
(close to but at the end of H burning). The mass of this object is 112680
$M_{\odot}$. The temperature here is expressed in Kelvin and the density is in
$gcm^{-3}$. The minimum value of radius is taken from the surface of the star
and the maximum value is 2.15 light years. Right panel:Variation of optical
depth from the surface to 100 parsecs for 1 $M_{\odot}$/yr model at different
masses during its evolution. The value for kappa is taken constant and
therefore this plot corresponds to the upper limit of such an optical depth.
Figure 12: Variation of optical depth from the surface to 100 parsecs for
$10^{-3}$ $M_{\odot}$/yr model at different masses during its evolution. The
value for kappa is taken constant and therefore this plot corresponds to the
upper limit of such an optical depth. Figure 13: A grid for the model computed
with accretion rate taken from Yoshida e al 2006 shows various parameters at 4
different times in the early evolutionary history (for now marked by the moel
number) versus the radius in centimeter. The first row shows the nuclear
energy production, the second row shows the relative abundance of Hydrogen,
Deuterium, Helium 3 and Lithium and the final row shows the central abundance
of C12, N15 and O16.

##  7 Discussion and future perspectives

Devesh: Expand on this section as much as you can.

  * •

You can mention that you only consider a constant accretion rate and state
that you leave this for future work.

  * •

Discuss possible observational signatures of PopIII SMS. Can we actually
observe them direcly in the future? Under what conditions (e.g. what accretion
rate?)

  * •

Add some content about ionising feedback. We have to do the calculations
first.

May want to touch on nucleosynthesis impact, speculate on the possible impact
of rotation

##  8 Conclusions

For example:

  1. 1.

We computed a grid of models with a wide variety of accretion rates and
masses….

  2. 2.

Our models suggest that the maximum possible mass of a population III star is
$\sim 2\times 10^{5}$ $\mathrm{M}_{\odot}$ (replace with correct value). This
is limited by the timescale of nuclear burning vs. the accretion timescale and
by the GR instability.

  3. 3.

We discuss the impact of the accretion rate on the ionising radiation produced
by population III stars.

  4. 4.

The final fate of these stars …

  5. 5.

The internal structure ..


# Tropical neural networks and its applications to classifying phylogenetic
trees

Ruriko Yoshida  Georgios Aliatimis  Keiji Miura

###### Abstract

Deep neural networks show great success when input vectors are in an Euclidean
space. However, those classical neural networks show a poor performance when
inputs are phylogenetic trees, which can be written as vectors in the tropical
projective torus. Here we propose tropical embedding to transform a vector in
the tropical projective torus to a vector in the Euclidean space via the
tropical metric. We introduce a tropical neural network where the first layer
is a tropical embedding layer and the following layers are the same as the
classical ones. We prove that this neural network with the tropical metric is
a universal approximator and we derive a backpropagation rule for deep neural
networks. Then we provide TensorFlow 2 codes for implementing a tropical
neural network in the same fashion as the classical one, where the weights
initialization problem is considered according to the extreme value
statistics. We apply our method to empirical data including sequences of
hemagglutinin for influenza virus from New York. Finally we show that a
tropical neural network can be interpreted as a generalization of a tropical
logistic regression.

##  1 Introduction

A neural network is a learning method, called deep learning, to learn data in
a way to mimic a brain system, i.e., which interconnects nodes, called
neurons, in a layered structure like a human brain system [13, 8, 15]. Recent
years deep neural networks show great success to process input data which lay
in the Euclidean space [11]. However, when input data are phylogenetic trees
or the time series with trends, represented as vectors in the tropical
projective torus [26, 27, 23, 38, 39, 40, 29, 42, 35], classical neural
networks show a poor performance. Therefore in this paper, we propose neural
networks which process an input data as vectors over the tropical projective
torus. The tropical projective torus denoted by
$\mathbb{R}^{d}/\mathbb{R}\mathbf{1}$ is the $d$-dimensional real numbers,
$\mathbb{R}^{d}$, mod by the vector with all ones, i.e.,
$\mathbf{1}:=(1,1,\ldots,1)\in\mathbb{R}^{d}$. Over the tropical projective
torus denoted by $\mathbb{R}^{d}/\mathbb{R}\mathbf{1}$, we define
$x:=(x_{1},x_{2},\ldots,x_{d})=(x_{1}+c,x_{2}+c,\ldots,x_{d}+c)\in\mathbb{R}^{d}/\mathbb{R}\mathbf{1}$
where $c\in\mathbb{R}$ [20]. Here we consider the tropical metric, also known
as the generalized Hilbert projective metric, over the tropical projective
torus as activation functions in a hidden layer of a neural network. It is
important to keep the invariance of the input vector under the one-vector
which is innate in the tropical projective torus [20, 17, 33, 18, 41, 30]. Our
strategy is to embed an input vector in the tropical projective torus into a
vector in the classical Euclidean space in the first layer. This is analogous
to the word embedding in the field of natural language processing [37, 28].
Then the following layers can be the same as the classical ones.

Although some previous works analyzed ReLU neural networks by using the
tropical geometry, the neural networks themselves are defined on a classical
Euclidean space [43, 1, 24]. In this paper, on the other hand, we consider a
tropical projective torus as an input space and keep the invariance under the
one-vector. That is, our work is truly tropical.

In this paper, we first introduce a tropical embedding layer. We use the
tropical embedding layer as the first layer of the classical neural networks
to keep the invariance in the tropical projective space. To check if this
tropical neural network has enough flexibility, we next prove that this
tropical neural network is a universal approximator. Then we derive a
backpropagation rule for the tropical neural networks. We provide TensorFlow 2
codes for implementing a tropical neural network in the same fashion as the
classical one, where the weights initialization problem is considered
according to the extreme value statistics. We show its applications to
phylogenomics, a new field in evolutionary biology which applies tools from
phylogenetic trees to genome data. Applications includes simulated data under
the multi-species coalescent model which is the most popular model to analyze
gene tree analysis on a genome data [21], and empirical data of influenza
virus data set collected from the state of New York [27]. Finally we briefly
show that a tropical neural network can be interpreted as a generalization of
a tropical logistic regression.

##  2 Tropical Embedding for Tropical Neural Networks

The classical neural networks only accept a input vector in an Euclidean space
in its original form. Thus they cannot accept a phylogenetic tree as an input
since a space of phylogenetic trees is not Euclidean [33, 3, 6], for example.
Therefore, we first consider a tropical embedding layer, which is analogous to
the word embedding in natural language processing [37, 28]. Once a
phylogenetic tree is embedded in an Euclidean space, a classical neural
network can is applied to analyzing it.

######  Definition 1 (tropical neural networks).

A tropical neural network is a network where a tropical embedding layer as the
first hidden layer is followed by a classical neural network (classical
layers).

######  Definition 2 (tropical embedding layer).

Let $x$ in $\mathbb{R}^{d}/\mathbb{R}\mathbf{1}$ be an input vector to the
tropical embedding layer. the activity of $j$-th neuron as an output of the
tropical embedding layer is given by

$z_{j}=\max_{i}(x_{i}+w^{(1)}_{ji})-\min_{i}(x_{i}+w^{(1)}_{ji}).$ (1)

######  Remark 3.

Note that no activation function is executed for $z$ as “max - min” operation
is somehow regarded as the activation function of the neurons in the first
hidden layer.

######  Remark 4.

There is a geometric interpretation: “max - min” operation measures the
distance between the points $x$ and $w^{(1)}_{j}$. Therefore $z(x)$ is
invariant along one vectors $\mathbf{1}$.

######  Remark 5.

There are alternative ways to attain the invariance such as

$z_{j}=\max_{i}(x_{i}+w^{(1)}_{ji})-\textrm{2nd}\max_{i}(x_{i}+w^{(1)}_{ji}).$
(2)

There is a geometric interpretation: “max - 2nd max” operation measures the
distance between a point $x$ and the tropical hyperplane whose normal vector
is $w^{(1)}$ [17]. Therefore $z(x)$ is invariant along one vectors
$\mathbf{1}$. You could even use $j$-th max in general. However, the
repertoire of functions never increase by using alternative ones. That is,
from the view point of universal approximation theorem, using Eq. (1)
suffices. In addition, Eq. (1) seems to perform better than the alternative
ones according to our numerical experiments (not shown). Therefore we solely
use Eq. (1) as a tropical embedding layer in what follows.

######  Remark 6.

Suppose $A\in\mathbb{Z}_{+}^{N\times d}$. Then we consider the ReLU such that

$\max\\{Ax+b,0\\}.$

Assume that $A{\bf 1}\not=0$. Suppose $x\in\mathbb{R}^{d}/\mathbb{R}{\bf 1}$.
Then we have $x:=x+c\cdot(1,\ldots,1)=x+c\cdot{\bf
1}\in\mathbb{R}^{d}/\mathbb{R}{\bf 1}$. Then for $c\ll 0$ and fixed $x$, we
have:

$\begin{array}[]{ccl}\max\\{Ax+cA{\bf 1}+b,0\\}&=&0.\\\ \end{array}$

As $c\to-\infty$, we have

$\frac{1}{1+\exp(-\max\\{Ax+cA{\bf 1}+b,0\\})}\to\frac{1}{1+1}=1/2$

for any $x\in\mathbb{R}^{d}/\mathbb{R}{\bf 1}$. Also for $c\gg 0$ and fixed
$x$, we have:

$\begin{array}[]{ccl}\max\\{Ax+cA{\bf 1}+b,0\\}&=&Ax+cA{\bf 1}+b.\\\
\end{array}$

As $c\to\infty$, we have

$\frac{1}{1+\exp(-(Ax+cA{\bf 1}+b))}\to 1$

for any $x\in\mathbb{R}^{d}/\mathbb{R}{\bf 1}$. Therefore, neural networks
with the ReLU cannot learn from observation in these cases. However, with the
activator function defined in Eq. (1), we have

$\max_{i}(x_{i}+c\cdot{\bf 1}+w^{(1)}_{ji})-\min_{i}(x_{i}+c\cdot{\bf
1}+w^{(1)}_{ji})=\max_{i}(x_{i}+w^{(1)}_{ji})-\min_{i}(x_{i}+w^{(1)}_{ji}).$

######  Remark 7.

Classical neural networks are not well-defined in the tropical projective
torus, since the neuron values are not invariant under transformations of the
form $x\to x+(c,\dots,c)$. Meanwhile, the tropical embedding layer of Eq. (1)
is invariant under such transformations.

##  3 Universal Approximation Theorems for Tropical Neural Networks

It is very important to check if the tropical embedding layer as in Eq. (1)
followed by classical layers has enough varieties to represent considerable
input-output relations [8]. In this section, we show that the tropical neural
network can approximate enough variety of functions so that we can safely use
it.

######  Definition 8.

The norm $\|\cdot\|_{q}$ for $q\geq 1$ is defined by

$\|f(x)\|_{q}=\int_{\mathbb{R}^{n}}|f(x)|^{q}dx$ (3)

The space $L^{q}(\mathbb{R}^{d}),(1<q<\infty),$ is the set of Lebesgue
integrable functions $f$ from $\mathbb{R}^{d}$ to $\mathbb{R}$ for which
$\|f(x)\|_{q}<\infty$.

######  Definition 9.

The space $C^{0}(\mathbb{R}^{d})$ is the set of continuous, compactly
suppported functions from $\mathbb{R}^{d}$ to $\mathbb{R}$.

######  Remark 10.

Note that $C^{0}(\mathbb{R}^{d})\subset L^{q}(\mathbb{R}^{d})$.

For the classical case, a universal approximation theorem for ReLU feedforward
neural nerworks has been proved in [4].

######  Theorem 11 (classical universal approximation theorem [4]).

Any function of $x_{j}$ for $j=1,\ldots,d$ in
$L^{q}(\mathbb{R}^{d}),(1<q<\infty),$ can be arbitrarily well approximated in
the $\|\cdot\|_{q}$ by a ReLU feedforward neural network with at most
$L=2(\lfloor\log_{2}d\rfloor+2)$ layers.

As the $d-1$ neurons in the tropical embedding layer can easily represent
$(x_{j}-x_{d})$ for $j=1,\ldots,d-1$ and Theorem 11 can be applied to the
second and later layers of a tropical neural network (that is equivalent to a
classical neural network), we can prove the following theorem.

######  Theorem 12 (tropical universal approximation theorem).

Any function of $(x_{j}-x_{d})$ for $j=1,\ldots,d-1$ in
$L^{q}(\mathbb{R}^{d}/\mathbb{R}\mathbf{1})\simeq
L^{q}(\mathbb{R}^{d-1}),(1<q<\infty),$ can be arbitrarily well approximated in
the $\|\cdot\|_{q}$ by a tropical neural network with at most
$L=2(\lfloor\log_{2}d\rfloor+2)+1$ layers (which include an tropical embedding
layer as the first layer).

###### Proof.

For any $f\in L^{q}(\mathbb{R}^{d-1})$, $\exists g\in C_{0}(\mathbb{R}^{d-1})$
such that $\|f-g\|_{q}<\epsilon/2$ [8]. Let $K$ be the support of $g$ and let
$M$ be $\max_{x\in K}\|x\|$. For $x\in K$, we can set
$w^{(1)}_{jj}=-w^{(1)}_{jd}=2M$ and $w^{(1)}_{ji}=0$ for $i\neq j,d$ to obtain
$z_{j}=x_{j}-x_{d}+4M$ for $j=1,\ldots,d-1$. This means that a neuron in the
first tropical embedding layer can represent $x_{j}-x_{d}$. Then $d-1$ neurons
can represent $z_{1}$, $z_{2}$, …, $z_{d-1}$. Finally, simply apply Theorem 11
to the classical neural network $F(z_{1},\ldots,z_{d-1})$ consisting of the
second and later layers of a tropical neural network to obtain
$\|g-F\|_{q}<\epsilon/2$. Taken together,
$\|f-F\|_{q}<\|f-g\|_{q}+\|g-F\|_{q}<\epsilon$. ∎

There is another type of classical universal approximation theorems.

######  Definition 13.

The width $d_{m}$ of a neural network is defined to be the maximal number of
neurons in a layer.

######  Theorem 14 (classical universal approximation theorem for width-
bounded ReLU networks [19]).

For any $f\in L^{1}(\mathbb{R}^{d})$ and any $\epsilon>0$, there exists a
classical neural networks $F(x)$ with ReLU activation functions with width
$d_{m}\leq d+4$ that satisfies

$\int_{\mathbb{R}^{d}}|f(x)-F(x)|dx<\epsilon.$ (4)

Again, as the $d-1$ neurons in the tropical embedding layer can easily
represent $(x_{j}-x_{d})$ for $j=1,\ldots,d-1$ and Theorem 14 can be applied
to the second and later layers of a tropical neural network (that is
equivalent to a classical neural network), we can prove the following theorem.

######  Theorem 15 (tropical universal approximation theorem with bounded
width).

For any function $f$ of $(x_{j}-x_{d})$ for $j=1,\ldots,n-1$ in
$L^{1}(\mathbb{R}^{d}/\mathbb{R}\mathbf{1})\simeq L^{1}(\mathbb{R}^{d-1})$ and
any $\epsilon>0$, there exists a tropical neural networks $F(x)$ with width
$d_{m}\leq d+4$ that satisfies

$\int_{\mathbb{R}^{d-1}}|f(x)-F(x)|dx<\epsilon.$ (5)

###### Proof.

For any $f\in L^{1}(\mathbb{R}^{d-1})$, $\exists g\in C_{0}(\mathbb{R}^{d-1})$
such that $\|f-g\|_{1}<\epsilon/2$ [8]. Let $K$ be the support of $g$ and let
$M$ be $\max_{x\in K}\|x\|$. For $x\in K$, we can set
$w^{(1)}_{jj}=-w^{(1)}_{jd}=2M$ and $w^{(1)}_{ji}=0$ for $i\neq j,d$ to obtain
$z_{j}=x_{j}-x_{d}+4M$ for $j=1,\ldots,d-1$. This means that a neuron in the
first tropical embedding layer can represent $x_{j}-x_{d}$. Then $d-1$ neurons
can represent $z_{1}$, $z_{2}$, …, $z_{d-1}$. Finally, simply apply Theorem 14
to the classical neural network $F(z_{1},\ldots,z_{d-1})$ consisting of the
second and later layers of a tropical neural network to obtain
$\|g-F\|_{q}<\epsilon/2$. Taken together,
$\|f-F\|_{q}<\|f-g\|_{q}+\|g-F\|_{q}<\epsilon$. ∎

##  4 Backpropagation Rule for Simplest Tropical Neural Networks

Figure 1: architecture of a simplest neural networks that accept a vector in
$\mathbb{R}^{d}/\mathbb{R}\mathbf{1}$

Here we demonstrate that the gradients of the loss function with respect to
weights exist for tropical neural networks. The gradient is computable through
the chain rule for differentials called backpropagation rule in the similar
fashion to the classical case. The gradients obtained in this way can
guarantee the successful update of the weights at each iteration of learning.

We consider a simplest three layer network whose weights in the first and the
second layers are denoted as $w^{(1)}\in\mathbb{R}^{d\times N}$ and
$w^{(2)}\in\mathbb{R}^{N\times 1}$. Suppose the activity in the first hidden
layer is given by Eq. 1 and the output of the network is given as

$y=\sum_{j}^{N}w^{(2)}_{j}z_{j}.$ (6)

Note that although here we derive a backpropagation rule for this regression
setting just for simplicity, the backpropagation rule can be derived in a
similar manner for the classification setting with a sigmoid function, too.

Below is the summary of parameters of the neural network.

  * •

$w^{(1)}$, $w^{(2)}$: weights in the first and second layers;

  * •

$z_{j}$: activation of $j$-th neuron in the hidden layer; and

  * •

$y$: activation of the output neuron.

######  Theorem 16.

The partial derivatives of the cost function
$Q:=\frac{1}{2}(y-y^{\textrm{true}})^{2}$ with respect to weights for the
above tropical neural network $y=f(x)$ are given by

$\frac{\partial Q}{\partial
w^{(1)}_{ji}}=(y-y^{\textrm{true}})w^{(2)}_{j}(\delta(i=i^{\textrm{max}}_{j})-\delta(i=i^{\textrm{min}}_{j})),$
(7)

where $i^{\textrm{max}}_{j}$ (or $i^{\textrm{min}}_{j}$) is the index $i$ for
which $(x_{i}+w^{(1)}_{ji})$ takes the maximum (or minimum) and

$\frac{\partial Q}{\partial w^{(2)}_{j}}=(y-y^{\textrm{true}})z_{j}.$ (8)

###### Proof.

Direct calculations. ∎

######  Example 17.

As a simplest example of Eq. (7), let us consider the three dimensional input
case $(d=3)$. Suppose the number of neurons in the middle layer is one and its
activity is z, for simplicity. Assume $x_{1}=1$, $x_{2}=2$ and $x_{3}=3$ and
$w^{(1)}_{1}=w^{(1)}_{2}=w^{(1)}_{3}=0$. Then
$(x_{i}+w^{(1)}_{1})<(x_{i}+w^{(1)}_{2})<(x_{i}+w^{(1)}_{3})$ and
$i_{\textrm{max}}=3$ and $i_{\textrm{min}}=1$. Therefore,

$\frac{\partial Q}{\partial
w^{(1)}_{i}}=\left\\{\begin{array}[]{ll}-(y-y^{\textrm{true}})w^{(2)}&(i=1)\\\
0&(i=2)\\\ (y-y^{\textrm{true}})w^{(2)}&(i=3).\end{array}\right.$ (9)

In this case we have $z=3-1=2$ If, furtheremore, $w^{(2)}=1$, then,
$y=w^{(2)}z=2$ and

$\frac{\partial Q}{\partial
w^{(1)}_{i}}=\left\\{\begin{array}[]{ll}-(2-y^{\textrm{true}}),&(i=1)\\\
0,&(i=2)\\\ 2-y^{\textrm{true}}.&(i=3),\end{array}\right.$ (10)

$w^{(1)}_{i}$ can be, for example, updated by the SGD rule: $\Delta
w^{(1)}_{i}=-\eta\frac{\partial Q}{\partial w^{(1)}_{i}}$, where $\eta>0$ is a
learning rate. Then, $w^{(1)}_{3}$ increases (and $w^{(1)}_{1}$ decreases) if
$2>y^{\textrm{true}}$.

######  Remark 18.

It is interesting that only two of $w^{(1)}_{i}$ are modified while the others
remains. Note that $\Delta w^{(1)}$ is orthogonal to the one vector,
$\mathbf{1}:=(1,1,\ldots,1)\in\mathbb{R}^{d}$. It is interesting to elucidate
how this learning rule works as a dynamical system.

##  5 TensorFlow2 Codes for Tropical Neural Networks

In order to boost computing with GPUs, we implement tropical neural networks
in TensorFlow 2 [9]. As is the case for the classical neural networks, the
auto-differential is the key for the GPU implementation of tropical neural
networks. In order to guarantee fast auto-differentials, all the calculations
must be implemented only with the math functions in TensorFlow2 such as
top_k($v$,$d$), which returns the maximum and the minimum of a vector $v$.

In practice, it is essential to create a user-friendly class for the tropical
embedding as the first layer of the tropical neural networks, that is scalable
for big data. The following code defines a hand-made class called TropEmbed(),
which enables us to easily implement the tropical neural networks in the
Keras/Tensorflow style.

⬇

class TropEmbed(Layer):

def __init__(self, units=2, input_dim=3):

super(TropEmbed, self).__init__()

self.w = self.add_weight(shape=(units, input_dim), \

initializer="random_normal")

self.units = units

self.input_dim = input_dim

def call(self, x):

x_reshaped = tf.reshape(x,[-1, 1, self.input_dim])

Bcast = repeat_elements(x_reshaped, self.units, 1)

val, i = tf.math.top_k(Bcast + self.w, self.input_dim)

return val[:,:,0] - val[:,:,-1]

# usage

model = Sequential([TropEmbed(10, d), Dense(1)])

The codes for TropEmbed() class and for reproducing all the figures in this
paper are available at https://github.com/keiji-miura/TropicalNN.

##  6 Weight Initialization Based on Extreme Value Statistics

Weight initializations are important for avoiding the divergence and
vanishment of neural activities after propagating many layers. For the
classical neural networks, Xavier’s and He’s initializations are famous [12,
16]. Here we consider a tropical analogue.

######  Definition 19 (Generalized Hilbert Projective Metric).

For any points
$v:=(v_{1},\ldots,v_{d}),\,w:=(w_{1},\ldots,w_{d})\in\mathbb{R}^{d}\\!/\mathbb{R}{\bf
1}$, the tropical distance (also known as tropical metric) $d_{\rm tr}$
between $v$ and $w$ is defined as:

$d_{\rm
tr}(v,w):=\max_{i\in\\{1,\ldots,d\\}}\bigl{\\{}v_{i}-w_{i}\bigr{\\}}-\min_{i\in\\{1,\ldots,d\\}}\bigl{\\{}v_{i}-w_{i}\bigr{\\}}.$

######  Lemma 20.

Suppose $x_{i},w_{i}\sim N(0,1)$ for $i=1,\ldots,d$. Then the expectation and
variance of $d_{\rm tr}(x,-w)$ can be approximated by
$2\sqrt{2}(a_{d}\gamma+b_{d})$ and $\frac{\pi^{2}}{3\log d}$, respectively,
where $a_{d}=\frac{1}{\sqrt{2\log d}}$ and $b_{d}=\sqrt{2\log
d}-\frac{\log\log d+\log(4\pi)}{2\sqrt{2\log d}}$.

###### Proof.

As $x_{i}+w_{i}\sim N(0,2)$,
$Z:=\frac{\max\\{x+w\\}/\sqrt{2}-b_{d}}{a_{d}}\sim Gumbel(0,1)$ as
$d\to\infty$. Therefore, $\mathrm{Ex}[d_{\rm
tr}(x,-w)]=\mathrm{Ex}[2\max\\{x+w\\}]\xrightarrow[d\to\infty]{}2\sqrt{2}(a_{d}\mathrm{Ex}[Z]+b_{d})$.
$\mathrm{Var}[d_{\rm
tr}(x,-w)]=\mathrm{Var}[\max\\{x+w\\}-\min\\{x+w\\}]=2\mathrm{Var}[\max\\{x+w\\}]+2\mathrm{Cov}[\max\\{x+w\\},-\min\\{x+w\\}]\xrightarrow[d\to\infty]{}2\times
2a_{d}^{2}\mathrm{Var}[Z]=2\times 2a_{d}^{2}\frac{\pi^{2}}{6}$ where
$\mathrm{Cov}[\max\\{x+w\\},-\min\\{x+w\\}]\xrightarrow[d\to\infty]{}0$ was
assumed. ∎

Here we confirm that the above scaling holds actually by numerical
calculations.

Figure 2: histogram of simulated $d_{\rm tr}(x,-w)$ for the same situation as
in Lemma 20 with $d=10000$. For the histogram, $100000$ samples of $d_{\rm
tr}(x,-w)$ are used. The simulated mean is $10.893$ while the theoretical
prediction is $10.954$. The simulated std is $0.604$ while the theoretica
prediction is $0.598$. The mean and std are the same as predicted from the
theory.

One way for better weight initialization is to choose the scale of $w$ so that
the variance of the neural activity in the embedding layer becomes $1$.

######  Lemma 21.

Suppose $x_{i}\sim N(0,1)$ and $w_{i}\sim N(0,\frac{6\log d}{\pi^{2}}-1)$ for
$i=1,\ldots,d$. Then the expectation and variance of $d_{\rm tr}(x,-w)$ can be
approximated by $2\sqrt{\frac{6\log d}{\pi^{2}}}(a_{d}\gamma+b_{d})$ and $1$,
respectively.

###### Proof.

As $x_{i}+w_{i}\sim N(0,\frac{6\log d}{\pi^{2}})$,
$Z:=\frac{\max\\{x+w\\}/\sqrt{\frac{6\log d}{\pi^{2}}}-b_{d}}{a_{d}}\sim
Gumbel(0,1)$. Therefore, $\mathrm{Ex}[d_{\rm
tr}(x,-w)]=\mathrm{Ex}[2\max\\{x+w\\}]\xrightarrow[d\to\infty]{}2\sqrt{\frac{6\log
d}{\pi^{2}}}(a_{d}\mathrm{Ex}[Z]+b_{d})$. $\mathrm{Var}[d_{\rm
tr}(x,-w)]\xrightarrow[d\to\infty]{}2\times\frac{6\log
d}{\pi^{2}}a_{d}^{2}\mathrm{Var}[Z]=2\frac{6\log
d}{\pi^{2}}a_{d}^{2}\frac{\pi^{2}}{6}=1$. ∎

To control the standard deviation of the weights, you can customize an
initializer (instead of simply specifying initializer $=$ ”random_normal”) in
the definition of TropEmbed class.

⬇

> ini = tf.keras.initializers.RandomNormal(mean=0., stddev=1.)

> self.w = self.add_weight(shape=(units, input_dim), \

initializer=ini)

However, as the weght initialization should be done together with the data
preprocessing, in this paper we entirely use the default value of stddev=0.05
for ”random_normal” for simplicity.

⬇

> self.w = self.add_weight(shape=(units, input_dim), \

initializer="random_normal")

##  7 Computational Experiments

In this section, we apply tropical neural networks with one hidden layer (the
tropical embedded layer) to simulated data as well as empirical data. Then
later we compare its performance against neural networks with one hidden layer
with ReLU activator.

###  7.1 Small simulated data

First we illustrate our tropical neural networks with one hidden layer with 16
neurons and with one output Sigmoid function using a small example. First we
generate two dimensional $16+16$ random points from the Gaussian distributions
with means $(0.5,-0.5)$ and $(-0.5,0.5)$ with the unit covariance matrix. Then
these points are randomly translated by $(c,c)$ where $c$ is a Gaussian random
variable whose standard deviation is 4. The left and right figures in Figure 3
show the actual test labels and the predicted probabilities of the test data
by the tropical neural networks with one hidden layer with 16 neurons.

Figure 3: Predicted probabilities by the tropical neural networks on a small
example.

###  7.2 High-dimensional simulated data

Second we demonstrate that a tropical neural network with one hidden layer
with 8 neurons and with one output sigmoid function works against the curse of
dimensionality, where the most of the variables in this high dimensional data
are rather just noises [42]. We generate $d$ dimensional $16+16$ random points
from the Gaussian distributions with means $(0.5,-0.5,0,\ldots,0)$ and
$(-0.5,0.5,0,\ldots,0)$ with the unit covariance matrix. Then these points are
randomly translated by $(c,c,c,\ldots,c)$ where $c$ is a Gaussian random
variable whose standard deviation is 6. The result demonstrates that the
tropical neural networks work robustly against the curse of dimensionality.

Figure 4: Application of the tropical neural networks to a high-dimensional
example. The test accuracy averaged over 100 trials is plotted. The tropical
neural networks work robustly against the curse of dimensionality.

###  7.3 simulated data generated from the multi-species coalescent model

In this subsection we apply the tropical neural networks to a sample of
phylogenetic trees generated under the multi-species coalescent model.

A phylogenetic tree is a weighted tree whose leaves are labeled with
$[m]:=\\{1,2,\ldots,m\\}$, where $m$ is the number of leaves, and whose
internal nodes are unlabeled. A weight on each edge of a phylogenetic tree is
considered as a distance from a node to another node on the tree and in
evolutionary biology, a weight on an edge can be considered as a product of an
evolutionary time and mutationa rate [32]. On this paper we consider a rooted
phylogenetic tree with a leaf label set $[m]$. An equidistant tree with $m$
leaves is called equidistant tree if the total weights on the unique path from
its room to each leaf is the same for all leaf in $[m]$. Under the multi-
species coalescent model which can be used to analyze gene trees, which are
phylogenetic trees reconstructed from genes in a genome, we assume that gene
trees are all equidistant trees. Therefore in this paper we assume that all
phylogenetic trees are equidistant trees.

To conduct a statistical analysis on a set of phylogenetic trees, we consider
a space of phylogenetic trees with fixed $[m]$. A space of phylogenetic trees
on $[m]$ is a set of all possible phylogenetic trees with $[m]$ and it is
well-known that it is not Euclidean [33]. It is also well-known that the space
of all possible equidistant trees on $[m]$ with the tropical metric under the
max-plus algebra is a subspace of the tropical projective space [3, 39]. In
order to define the space of equidistant trees, first we define ultrametrics.
COnsider a map $u:[m]\times[m]\to\mathbb{R}$ such that $u(i,j)=u(j,i)$ and
$u(i,i)=0$. This map is called a dissimilarity map on $[m]$. If a
dissimilarity map $u$ satisfies that

$\max\\{u(i,j),u(i,k),u(j,k)\\}$

achieve at least twice, then we call $u$ an ultrametric.

######  Example 22.

Suppose $m=3$, i.e., $[3]=\\{1,2,3\\}$ and suppose

$u(1,2)=u(2,1)=1,u(1,3)=u(3,1)=1,u(2,3)=u(3,2)=0.5$

and $u(i,i)=0$ for all $i=1,2,3$. Since

$\max\\{u(1,2),u(1,3),u(2,3)\\}=1$

and it achieves twice, i.e., $u(1,2)=u(1,3)=1$. Thus, $u$ is an ultrametric.

Consider dissimilarity maps $u_{T}$ on a phylogenetic tree $T$ on $[m]$ such
that $u(i,j)$ is the total weights on the unique path from a leaf $i$ to a
leaf $j$ for all $i,j\in[m]$. Then we have the following theorem:

######  Theorem 23 ([7]).

Consider an equidistant tree $T$ on $[m]$. Then $u_{T}$ realizes an
equidistant tree $T$ on $[m]$ if and only if a dissimilarity map $u_{T}$ is
ultrametric.

######  Example 24.

Suppose we have an ultrametric from Example 22. An equidistant tree whose
dissimilarity maps are ultrametric in Example 22 is a rooted phylogenetic tree
with leaves $[3]=\\{1,2,3\\}$ shown in Figure 5.

Figure 5: An equidistant tree with the dissimilarity maps which are
ultrametric shown in Example 22.

Therefore, we consider the space of all possible ultrametrics on $[m]$ as the
space of equidistant trees on $[m]$. Then we have the following theorem:

######  Theorem 25 ([3]).

The space of ultrametrics on $[m]$ is the tropicalization of the linear
subspace defined by linear equations

$x_{ij}-x_{ik}+x_{jk}=0$

for $1\leq i\leq j\leq k\leq m$ by replacing sum with max operation and
replacing multiplication by a classical summation.

The space of ultrametrics on $[m]$ is a subspace of the tropical projective
space $(\mathbb{R}^{e}\cup\\{-\infty\\})/\mathbb{R}{\bf 1}$ where
$e=\binom{m}{2}$. Therefore, we apply our method, tropical neural networks, to
simulated data generated from the multi-species coalescent model using a
software Mesquite [21].

The multi-species coalescent model has two parameters: species depth and
effective population size. Here we fix the effective population size
$N_{e}=100000$ and we vary

$R=\frac{SD}{N_{e}}$

where $SD$ is the species depth. We generate species trees using the Yule
process. Then we use the multi-species coalescent model to generate gene trees
with a given species tree. In this experiment, for each $R$, we generate two
different set of $1000$ gene trees: In each set of gene tree, we have a
different species tree so that each set of gene trees is different from the
other. We conduct experiments with $R=0.25,0.5,1,2,5,10$.

Note that when we have a small ratio $R$, then gene trees become more like
random trees since the species tree constrains less on tree topologies of gene
trees. Thus it is more difficult when we have a small $R$ while if we have a
large $R$, then it is easier to classify since the species tree constrains
more on tree topologies of gene trees.

Figure 6: ROC curves for neural networks with ReLU and tropical neural
networks with one hidden layer. We conduct experiments with
$R=0.25,0.5,1,2,5,10$.

In this experiment, we set one hidden layer for each neural network: neural
network with ReLU activators and neural network with tropical activators. We
set the Sigmoid function in the output node in both neural networks. In each
neural network, we set $1000$ neurons in the hidden layer.

Figure 6 shows ROC curves for neural networks with ReLU and tropical neural
networks. In general tropical neural networks perform better than neural
networks with ReLU activator function.

###  7.4 Influenza data

In this subsection we apply our method to genomic data for 1089 full length
sequences of hemagglutinin (HA) for influenza A H3N2 from 1993 to 2017
collected in the state of New York obtained from the GI-SAID EpiFlu database
(www.gisaid.org). These collected data were aligned using muscle developed by
[10] with the default settings. Then we apply the neighbor-joining method with
the p-distance [31] to reconstruct a tree from each sequenced data. Each year
corresponds to the first season. We also apply KDETrees [38] to remove
outliers and a sample size of each year is about 20,000 trees.

We apply tropical neural networks and neural networks with ReLU with one
hidden layer with 10 neurons to all pairs of different years to see if they
are significantly different one year to the other. Heatmaps of accuracy rates
with the probability threshold $0.5$ and AUC values are shown in Figure 7.
Again, tropical neural networks outperform classical neural networks.

Figure 7: Heat maps for (top) classification rates with threshold $0.5$ and
(bottom) AUC values for classical neural networks with ReLU (left) and
tropical neural networks (right).

###  7.5 Tropical Neural Network as a Generalization of Tropical Logistic
Regression for Classification of Gene Trees

A tropical neural network can be interpreted as a generalization of a tropical
logistic regression. The tropical logistic regression model [2] is developed
for binary classification on gene trees and it has been proven to have higher
predictive power than classical logistic regression to classify phylogenetic
trees. It assumes that if $X$ is the covariate (gene tree), the binary
response random variable is $Y\sim\text{Bernoulli}(p(X))$, with

$p(x)=S\left(\lambda_{0}d_{\rm tr}(x,\omega_{0})-\lambda_{1}d_{\rm
tr}(x,\omega_{1})+C\right),$

where $m$ is the number of leaves in each phylogenetic tree in the given
sample, $e:=\binom{m}{2}$, $S$ is the sigmoid function,
$\omega_{0},\,\omega_{1}\in\mathbb{R}^{e}/\mathbb{R}{\bf 1}$, and
$\lambda_{0},\lambda_{1},C\in\mathbb{R}$ with $\lambda_{0}\lambda_{1}\geq 0$.
Note that this model is a special case for Eq. (6), with the sigmoid function
as the link, and two neurons in the hidden layer whose weights are
$w_{0}^{(2)}=\lambda_{0},w_{1}^{(2)}=-\lambda_{1}$ and
$w_{0}^{(1)}=-\omega_{0},w_{1}^{(1)}=-\omega_{1}$. Therefore, tropical
logistic regression is almost identical to a tropical neural network
consisting of one tropical embedding layer with two neurons and a classical
layer, with the additional assumption that $w_{0}^{(2)}w_{1}^{(2)}\leq 0$.

The one-species model described in [2] can be considered to be a neural
network with $e$ neurons in the input layer, no hidden layers and a unique
output neuron. The activation function is the logistic function and the inner
product used is tropical defined as

$\langle x,-\omega\rangle:=d_{\rm tr}(x,\omega)-C,$ (11)

where $C$ can be considered to be a bias variable, similarly to the intercept
variable in classical models. Tropical logistic regression returns the sigmoid
of the tropical inner product. We define the tropical generalised linear model
as an extension to tropical logistic regression, where instead of the sigmoid
function, we may use a different link/activation function. If there are
multiple outputs (multivariate generalized linear model (GLM)) and if we treat
the output layer as the new input layer and iterate this $L$ times, then we
have an $L$-layer neural network. In the same way that classical neural
networks are a stack/recursive application of classical multivariate GLMs,
tropical neural networks can be a stack of tropical multivariate GLMs.
Effectively, all is identical to classical networks, but instead of applying
classical inner products, we apply tropical inner products as defined in Eq.
(11). The $i$-th neuron of the $l$-th layer is defined as $x^{(l)}_{i}$ and
computed through the recursive formula,

$x^{(l)}_{i}=d_{\rm
tr}\left(x^{(l-1)}_{i},\omega^{(l)}_{i}\right)-C^{(l)}_{i},$ (12)

where
$\Omega^{(l)}=(\omega^{(l)}_{1},\omega^{(l)}_{2},\dots,\omega^{(l)}_{N_{l}})\in\mathbb{R}^{N_{l-1}\times
N_{l}}$ is the weight matrix between layer $(l-1)$ and $l$ for the number
$N_{s}$ of neurons in layer $s$, and $C^{(l)}\in\mathbb{R}^{N_{l}}$. By
assuming that all neurons share the same bias variable $c=C^{(l)}_{i}$ for all
$i\in[N_{l}]$, Eq. (12) reduces to Eq. (1), since vectors are defined up to an
additive constant vector $(c,\dots,c)$ in the tropical projective torus. When
the last tropical embedding layer connects to the first classical layer, the
constant bias vector is incorporated in the bias term of the classical layer.
Hence, tropical bias terms are redundant and not considered in the development
of tropical neural networks. Thus, the tropical neural network which we
propose in this paper follows naturally as an extension of the tropical
logistic regression model.

##  8 Summary and Discussion

In this paper, we first developed a tropical embedding layer. We used the
tropical embedding layer as the first layer of the classical neural networks
to keep the invariance in the tropical projective space. To check if this
tropical neural network has enough flexibility, we next proved that this
tropical neural network is a universal approximator. After we derived a
backpropagation rule for the tropical neural networks, we provided TensorFlow
2 codes for implementing a tropical neural network in the same fashion as the
classical one, where the weights initialization problem is considered
according to the extreme value statistics. Finally we showed some applications
as examples.

The tropical neural networks with the tropical metric worked better than the
classical neural networks when the input data are phylogenetic trees which is
included in the tropical projective torus. This is partly because only the
tropical neural network can keep the invariance of the input vector under the
one-vector which is innate in the tropical projective torus.

One of the nice properties of tropical neural networks is its tractability and
interpretability in analysis. The tropical embedding can be interpreted as
taking the tropical distance to a point in the space of the tropical
projective torus. The activities of neurons in the tropical neural networks
with the randomized weights and inputs can be analyzed by using the extreme
value statistics. The backpropagation rule of the tropical neural networks can
be derived and interpreted rather easily.

The TensorFlow 2 codes for the Python class for tropical embedding was
provided in the paper. This makes it possible to implement a tropical neural
network in the same familiar fashion as the classical one. This facilitates,
for example, to compare tropical and classical neural networks for the same
data by using a common code.

Recent work shows that neural networks are vulnerable against adversarial
attacks (for example, [5, 25, 34, 22]). However, our initial computational
experiments on image data from computer vision show that tropical neural
networks are robust against gradient based methods, such as the Fast Gradient
Sign Method [14] and Ensemble Adversarial Training [36]. It is interesting to
investigate why tropical neural networks are robust against such attacks. In
addition, it is interesting to develop adversarial attacks toward tropical
neural networks.

## References

  * [1] M Alfarra, A Bibi, H Hammoud, M Gaafar, and B Ghanem.  On the decision boundaries of neural networks: A tropical geometry perspective.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2022\.
  * [2] Georgios Aliatimis, Ruriko Yoshida, Burak Boyaci, and James A Grant.  Tropical logistic regression model on space of phylogenetic trees.  arXiv preprint arXiv:2306.08796, 2023.
  * [3] F. Ardila and C. J. Klivans.  The bergman complex of a matroid and phylogenetic trees. journal of combinatorial theory.  Series B, 96(1):38–49, 2006.
  * [4] R. Arora, A. Basu, P. Mianjy, and A. Mukherjee.  Understanding deep neural networks with rectified linear units.  ICLR, 2018.
  * [5] Battista Biggio, Igino Corona, Davide Maiorca, Blaine Nelson, Nedim Šrndić, Pavel Laskov, Giorgio Giacinto, and Fabio Roli.  Evasion attacks against machine learning at test time.  In Hendrik Blockeel, Kristian Kersting, Siegfried Nijssen, and Filip Železný, editors, Machine Learning and Knowledge Discovery in Databases, pages 387–402, Berlin, Heidelberg, 2013. Springer Berlin Heidelberg.
  * [6] L.J. Billera, S.P. Holmes, and K. Vogtmann.  Geometry of the space of phylogenetic trees.  Adv Appl Math, 27(4):733–767, 2001.
  * [7] P. Buneman.  A note on the metric properties of trees.  J. Combinatorial Theory Ser. B., 17:48–50, 1974.
  * [8] Ovidiu Calin.  Deep Learning Architectures: A Mathematical Approach.  Springer, 2020.
  * [9] Francois Chollet.  Deep Learning with Python 2nd Edition.  Manning, 2021.
  * [10] R.C. Edgar.  Muscle: multiple sequence alignment with high accuracy and high throughput.  Nucleic Acids Research, 32:1792–1797, 2004.
  * [11] Martin Ford.  Architects of Intelligence: The truth about AI from the people building it.  Packt Publishing, 2018.
  * [12] Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, PMLR, 9:249–256, 2010.
  * [13] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  The MIT Press, 2016.
  * [14] Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and harnessing adversarial examples.  arXiv preprint arXiv:1412.6572, 2014.
  * [15] Philipp Grohs and Gitta Kutyniok.  Mathematical Aspects of Deep Learning.  Cambridge Univ Press, 2023.
  * [16] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016.
  * [17] Michael Joswig.  Essentials of Tropical Combinatorics.  Springer, New York, NY, 2021.
  * [18] B. Lin and R. Yoshida.  Tropical Fermat-Weber points.  SIAM Discrete Math., page https://arxiv.org/abs/1604.04674, 2018\.
  * [19] Zhou Lu, Hongming Pu, Feicheng Wang, Zhiqiang Hu, and Liwei Wang.  The expressive power of neural networks: A view from the width.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [20] D. Maclagan and B. Sturmfels.  Introduction to Tropical Geometry, volume 161 of Graduate Studies in Mathematics.  Graduate Studies in Mathematics, 161, American Mathematical Society, Providence, RI, 2015.
  * [21] W. P. Maddison and D.R. Maddison.  Mesquite: a modular system for evolutionary analysis. version 2.72, 2009\.  Available at http://mesquiteproject.org.
  * [22] Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu.  Towards deep learning models resistant to adversarial attacks.  In 6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings. OpenReview.net, 2018.
  * [23] A. Monod, B. Lin, Q. Kang, and R. Yoshida.  Tropical foundations for probability & statistics on phylogenetic tree space, 2019.
  * [24] Guido Montúfar, Yue Ren, and Leon Zhang.  Sharp bounds for the number of regions of maxout networks and vertices of minkowski sums.  SIAM J. APPL. ALGEBRA GEOMETRY, 6(4):618–649, 2022.
  * [25] Anh Mai Nguyen, Jason Yosinski, and Jeff Clune.  Deep neural networks are easily fooled: High confidence predictions for unrecognizable images.  CoRR, abs/1412.1897, 2014.
  * [26] T. M. W. Nye.  Principal components analysis in the space of phylogenetic trees.  Ann. Stat., 39(5):2716–2739, 2011.
  * [27] T. M. W. Nye, X. Tang, G. Weyenberg, and R. Yoshida.  Principal component analysis and the locus of the fréchet mean in the space of phylogenetic trees.  Biometrika, 104:901–922, 2017.
  * [28] Aytuğ Onan.  Sentiment analysis on product reviews based on weighted word embeddings and deep neural networks.  Concurrency and Computation: Practice and Experience, 33(23):e5909, 2021.
  * [29] Robert Page, Ruriko Yoshida, and Leon Zhang.  Tropical principal component analysis on the space of phylogenetic trees.  Bioinformatics, 36(17):4590–4598, 2020.
  * [30] Jürgen Richter-Gebert, Bernd Sturmfels, and Thorsten Theobald.  First steps in tropical geometry, 2003.
  * [31] N. Saitou and M. Nei.  The neighbor-joining method: a new method for reconstructing phylogenetic trees.  Mol Biol Evol., 4:406–425, 1987.
  * [32] C. Semple and M. Steel.  Phylogenetics, volume 24 of Oxford Lecture Series in Mathematics and its Applications.  Oxford University Press, 2003, 2003.
  * [33] D. Speyer and B. Sturmfels.  Tropical mathematics.  Mathematics Magazine, 82:163–173, 2009.
  * [34] Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian J. Goodfellow, and Rob Fergus.  Intriguing properties of neural networks.  In Yoshua Bengio and Yann LeCun, editors, 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings, 2014.
  * [35] Xiaoxian Tang, Houjie Wang, and Ruriko Yoshida.  Tropical support vector machine and its applications to phylogenomics, 2020.
  * [36] Florian Tramèr, Alexey Kurakin, Nicolas Papernot, Ian J. Goodfellow, Dan Boneh, and Patrick D. McDaniel.  Ensemble adversarial training: Attacks and defenses.  In 6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings. OpenReview.net, 2018.
  * [37] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Ł ukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [38] G. Weyenberg, R. Yoshida, and D. Howe.  Normalizing kernels in the Billera-Holmes-Vogtmann treespace.  IEEE ACM T. Comput. Bi., page doi:10.1109/TCBB.2016.2565475, 2016\.
  * [39] R. Yoshida, L. Zhang, and X. Zhang.  Tropical principal component analysis and its application to phylogenetics.  Bulletin of Mathematical Biology, 81:568–597, 2019.
  * [40] Ruriko Yoshida, David Barnhill, Keiji Miura, and Daniel Howe.  Tropical density estimation of phylogenetic trees, 2022.
  * [41] Ruriko Yoshida, Keiji Miura, and David Barnhill.  Hit and run sampling from tropically convex sets, 2022.
  * [42] Ruriko Yoshida, Misaki Takamori, Hideyuki Matsumoto, and Keiji Miura.  Tropical support vector machines: Evaluations and extension to function spaces.  Neural Networks, 157:77–89, 2023.
  * [43] L Zhang, G Naitzat, and LH Lim.  Tropical geometry of deep neural networks.  International Conference on Machine Learning, pages 5824–5832, 2018\.