# Territorial differences in the spread of COVID-19 in European regions and US
counties

Fabrizio Natale111Corresponding author, Email<EMAIL_ADDRESS>222European Commission, Joint Research Centre, Via Enrico Fermi 2749, 21027
Ispra (VA), Italy Stefano Maria Iacus††footnotemark:  Alessandra Conte
††footnotemark:  Spyridon Spyratos††footnotemark:  Francesco
Sermi††footnotemark:

###### Abstract

This article explores the territorial differences in the onset and spread of
COVID-19 and the excess mortality associated with the pandemic, across the
European NUTS3 regions and US counties. Both in Europe and in the US, the
pandemic arrived earlier and recorded higher Rt values in urban regions than
in intermediate and rural ones. A similar gap is also found in the data on
excess mortality. In the weeks during the first phase of the pandemic, urban
regions in EU countries experienced excess mortality of up to 68pp more than
rural ones. We show that, during the initial days of the pandemic, territorial
differences in Rt by the degree of urbanisation can be largely explained by
the level of internal, inbound and outbound mobility. The differences in the
spread of COVID-19 by rural-urban typology and the role of mobility are less
clear during the second wave. This could be linked to the fact that the
infection is widespread across territories, to changes in mobility patterns
during the summer period as well as to the different containment measures
which reverse the causality between mobility and Rt.

Keywords— Mobile Positioning Data, Business to Government, COVID-19, Mobility

## Introduction

The current COVID-19 pandemic is creating severe social and economic
consequences, with some places experiencing disproportionately high levels of
mortality and economic losses. Urban regions, and particularly large cities,
have been severely affected by the spread of the pandemic in its early stages.
The evolving public discussion on the territorial impact of the pandemic
requires a greater understanding of the way the pandemic is affecting regions
that are diversely vulnerable and will require different recovery plans.
Analysis of the role of population density on the virus spread has led to
mixed results [1, 2, 3]. In addition to density, particular attention has also
been given to some factors related to the urban organisation that would make
some places vulnerable to infection in the first phase. In particular, these
are the connectivity of cities as hubs of national and international transport
systems [4, 5, 6], the structure of employment and industry [7, 8], as well as
overcrowded living conditions in some home environments [9, 10]. Overall, this
analysis is aimed at a deeper understanding of the links between COVID-19,
urban-rural typologies, territorial conditions, and mobility, which is
critical for designing effective public health policy responses. We first
explore the heterogeneity of COVID-19 patterns in its onset, spread, and
associated excess mortality by comparing the results by the level of
urbanisation of European regions and counties in the US. The results show that
the pandemic started earlier in urban regions than in intermediate and rural
areas. Urban regions had the highest Rt values in both Europe and the US
during the first wave, whereas rural counties were more affected than urban
counties in the second wave. Analysis of excess mortality, calculated using
Eurostat statistics and obtained from the difference between reported
fatalities and a baseline model based on historical data between 2011 and
2019, also shows a large gap by urbanisation level during the first wave, with
a median excess mortality up to 73% for urban regions, 18% for intermediate
regions, and 11% for rural regions. In a second phase, we empirically examine
the impact of domestic mobility on virus spread. We model population mobility
in European regions through a harmonised mobility index derived from mobile
phone data. We examine the geographical distribution of mobility changes
through regression models for the weeks in the first and second virus waves.
Our results show that, on the one hand, higher mobility explains most of the
variation in values in the weekly Rt during the first wave, with internal,
inbound, and outbound mobility positively affecting Rt. The effect of internal
mobility, in particular, is more pronounced than that of the degree of
urbanisation, and remains significant even when population and population
density are taken into account. On the other hand, the same regression models
replicated for the second wave weeks show a negative role of mobility on the
local spread of the virus, as well as a higher prevalence of the infection in
rural regions compared to large cities during the second wave. The paper is
organised as follows. The data section describes the data and methods used in
the analyses. In the results and discussion section, we present how the
COVID-19 pandemic was spread in rural, intermediate, and urban regions during
the first and the second wave, and finally, the conclusions are outlined in
the final section.

## Data and methods

In this section we present the data sources and the methods we used to assess
the spread of the COVID-19 pandemic in rural, intermediate and urban regions.
We calculated the reproductive number (Rt) as an indicator to assess how fast
the virus spread across different types of geographical areas. We estimated
the excess mortality to monitor in quantitative terms the evolution and
impacts of COVID-19 pandemic. We used fully anonymised and regionally
aggregated mobility data to get insights about the different regional mobility
patters. Finally, we fitted a linear regression model to study the
relationship between mobility and Rt during the first and the second wave.

### Rt

Rt is the main real-time indicator used to assess the evolution of the
pandemic, design containment measures and monitor their effectiveness. A time-
dependent reproduction number, Rt, was calculated with the R package R0
following a likelihood-based estimation procedure [11] with parameters for
generation time distribution reported in [12]. and on the basis of daily data
on confirmed COVID-19 cases at regional level downloaded through the ’COVID19’
R package. The daily Rt values were aggregated to weekly averages.

### Excess mortality

The baseline for mortality was calculated with Generalised Additive Models
fitted independently for each region. These models include a seasonal
component to account for the increase in mortality during the winter months
linked to influenza outbreaks, and a linear time trend to account for long-
term changes in mortality due to demographic dynamics. The excess mortality is
measured as difference between the reported data in 2020 and the estimated
baseline for all occurrences exceeding the lower or upper 95% confidence
intervals of the estimated baseline. The weekly mortality are from Eurostat
(demormweek3) and cover 900 regions in 26 EU MS and the UK with time series
spanning from 2001 to 2015.

### Mobility

In this study we used fully anonymised and aggregated mobility data shared
with the European Commission (EC) by European Mobile Network Operators (MNOs).
These mobility data comply with the ‘Guidelines on the use of location data
and contact tracing tools in the context of the COVID-19 outbreak’ by the
European Data Protection Board [13]. The mobility data are in the form of
Origin-Destination Matrix (ODM) [14, 15] and they provide valuable insights
into mobility patterns across geographical areas. The data has been used to
derive mobility insights and build tools to inform better targeted containment
measures, in a Mobility Visualisation Platform, available to the Member States
[16].
Given the high variation in the spatial and temporal aggregation across
countries and operators, the original ODMs are harmonised at standardised
spatial and temporal granularity to the derived Mobility Indicators. For more
information about the Mobility Indicators we direct the interested reader to
[17]. In this study, we further aggregated the Mobility Indicators at a weekly
temporal granularity, and we normalised them in both per capital and absolute
terms to enable a better cross-country comparison on the number of recorded
movements. This further normalisation was performed by setting for each
country and type of mobility (internal, inbound, outbound) the value of one
for the NUTS3 regions with the higher mobility over the reference time period
February to December 2020, and the value of zero to the NUTS3 regions with the
lowest mobility over the same time period.

### Regression

The fitting of the regression models was constrained by the necessity of
having regional data on COVID-19 cases for the calculation of Rt and mobility
indicators for the same periods. Data on population was obtained from Eurostat
(demorpjangrp3 and demord3dens). Overall the regressions are based on around
3500 observations in 654 regions for the first wave, and 10500 observations in
551 regions for the second wave. The regression models include country fixed
effects to control for country-specific and time-invariant factors that may
influence the spread of the disease.

## Results and Discussion

### The COVID-19 pandemic started earlier in urban regions

Figure 1 describes the pandemic onset in the NUTS3 regions of some European
countries and counties in the US, clustered by the three levels of
urbanisation (for the EU we use Eurostat NUTS3 rural-urban typologies and for
the US we use Rural-Urban Continuum Codes reduced to 3 classes). We measure
the pandemic onset in each region by the number of days between the
registration of the first 20 confirmed cases of Coronavirus disease and the
beginning of the year 2020. In both Europe and the United States, urban
regions are more vulnerable to the pandemic’s onset. The pandemic started
earlier in most urban regions, while we observe a later onset in intermediate
and rural regions in the first wave of the virus.

Figure 1: Onset of the pandemic across regions by rural-urban typology.

### The infection has spread faster in urban regions during the first wave

Figure 2 displays the Rt values for the first and second waves of the
pandemic. Rt is calculated from daily confirmed cases in 807 NUTS3 regions in
the UK, Netherlands, Germany, Italy, Spain, France, Czech Republic and Austria
(left) and in 3100 counties in the US (right). The indicator is averaged
across regions grouped by rural-urban type and aggregated by days since the
first reported case in each region (upper Figure), and weeks since the start
of the second wave of the pandemic (lower Figure). Looking first at the upper
figure, we observe that urban regions in Europe and the United States recorded
higher Rt values than those found in intermediate and rural regions at the
start of the pandemic. This indicates that the disease spread faster in urban
regions and that containment was more difficult in more densely populated
areas. Approximately 56 days after the start of the pandemic, we find a
general decline in the Rt and a reduction in the differences in Rt between the
three groups of regions. At the start of the pandemic, the rural-urban divide
in Rt values is more pronounced in the US counties. However, even in this
case, the disparity in the pandemic spread by level of urbanisation has
narrowed among the three regional groups, with the Rt index close to 1 at the
end of the first wave. The lower part of Figure 2 shows the median Rt values
across regions and counties in the weeks following the summer period, when the
pandemic began to spread in a second wave of infections. In the European
regions, we observe an initially higher Rt in the urban regions and increasing
and higher values in the intermediate and rural regions as the second wave
progresses. In contrast, in the US, rural and intermediate counties are the
most vulnerable to virus spread for most of the weeks during the second wave,
with a slight change in trend in the last weeks of the period.

Figure 2: Median Rt values in European regions and American counties by rural-
urban typology. The upper part presents the daily Rt values from the first
recorded case in each region during the first wave. The lower part represents
the weekly Rt values during the second wave.

### The excess mortality linked to COVID-19 is higher in the European urban
regions in the first wave

Figure 3 shows the trend in the excess mortality for the European regions
during the year 2020. The increase in weekly mortality compared to past trends
is used as an indirect measure to monitor the evolution of COVID-19. This
indicator has the downside of including fatalities not necessarily linked to
COVID-19, such as those caused by the saturation of hospital capacity, but has
the advantage of being less influenced by the underestimation of the real
infection rate due to asymptomatic cases or differences in testing strategies
over time and regions [18]. The bars in Figure 3 show the weekly total excess
mortality calculated from Eurostat statistics for most EU countries and the
UK. The excess mortality is obtained from the difference between the reported
fatalities and a modelled baseline estimated from historical data between 2011
and 2019. The number of weekly fatalities attributable to COVID-19 peaked at
the beginning of April, with about 41 400 deaths in excess compared to the
baseline.333This peak represents 21 600 more cases than the excess mortality
recorded in the same countries during the second week of January 2017,
corresponding to a particularly severe year for the seasonal flu. The lines in
the figure show the median excess mortality in the NUTS3 regions classified
according to their degree of urbanisation. At the peak of the pandemic in
third week of April, the median excess mortality in urban regions reached its
peak with an excess mortality of 73%, which was 58 pp higher than in
intermediate regions and 68 pp higher than in rural regions in the same week.
In the second wave of the pandemic, the disparities among regions appear less
pronounced. There is also a reverse in the trend of excess mortality, with
rural and intermediate regions having higher rates, 38% and 32% respectively,
than urban regions with an excess mortality rate of 26%.

Figure 3: Total excess mortality (bars) and median excess mortality by rural-
urban typology (lines) on a weekly basis.

### Mobility is higher in urban regions

One possible explanation for the higher Rt and excess mortality in urban
regions is that in cities the infection can spread more rapidly given the
higher population density, larger use of public transportation and higher
number of social interactions. The intensity of social interaction is
reflected in mobility indicators which can be calculated from social media and
mobile phones data. In fact, the relation between intensity of social
contacts, mobility and infection is at the basis of mobility restriction that
most governments have put in place to contain the pandemic. We analyse the
patterns of mobility within, from and toward European regions with anonymised
and aggregated mobile indicators derived from mobile phone data as described
in the Data and methods Section. Figure 4 shows the median patterns of weekly
mobility of 1033 NUTS3 regions in 22 EU countries, grouped by rural,
intermediate and urban typology, in absolute (upper chart) and per capita
terms (lower chart). The trends in the two charts in Figure 4 reflect the
implementation of generalised lockdown until April, the reopening during the
summer period and the new restrictions on mobility after summer. The upper
chart in Figure 4, reflects that the mobility in absolute terms in urban areas
is higher compared to the intermediate and rural ones, mainly due to their
different population size. The lower chart in Figure 4 describes that during
the first wave, and independently from the implementation of the restriction
measures, the level of per capita mobility was higher in urban regions in
respect of intermediate and rural ones. During the second wave, the per capita
mobility is almost equal across all areas, indicating substantial reduction of
mobility in urban and intermediate regions at the beginning of summer and an
increase in rural regions. This shift in mobility patterns is exemplified in
Figure 5 showing the weekly relative changes in mobility for each Italian
region (rows) in respect of the levels recorded during the last week of
February (first column). In this case, regions are sorted on the basis of
their proximity to the sea or mountains to better appreciate the mobility
linked to domestic tourism. In May, after the lifting of lockdown, all Italian
regions recorded an increase of mobility to the levels of February. However,
during summer, in coastal and mountain regions mobility increased to higher
values than at the beginning of the year. The highest increase was recorded in
the second week of August in the renowned region of Olbia in Sardinia (+373%).
With the re-opening of schools in September, the level of mobility started
again to increase uniformly across all regions.

Figure 4: Median absolute and per capita internal mobility across regions
grouped by rural-urban typology Figure 5: Changes in internal mobility for
regions in Italy in respect of the level in the last week of February

### The higher mobility in urban regions may explain great part of the
territorial gaps in Rt during the first wave

To support our intuition about the heterogeneity in the spread of Covid during
the first and second waves, we examine the effect of different mobility
patterns on the spread of the pandemic through a series of regression models
using the Rt values recorded in each European region as the dependent
variable. Table 1 shows the results of the ordinary least squares regressions
on the first wave of infection considering the Rt values in the first 28 days
after the start of the pandemic. Table 2 presents results for the second wave
on the Rt values in the weeks after August and Figure 6 provides estimates of
the relation between Rt and mobility using different temporal shifts between
the two variables. All specifications include country fixed effects to account
for differences in transmission resulting from invariant country
characteristics. We estimate a significant relationship between the effective
reproduction number, Rt, and the levels of urbanisation (Column 1 in Table 1).
During the first wave of the pandemic, Rt values are lower in rural and
intermediate regions than in the urban regions used as reference. In Columns
2-4 we include the mobility controls separately, i.e. internal, incoming and
outgoing mobility, given the correlation between these measures within
countries. In all specifications that include controls separately, each
mobility indicator is positively correlated with Rt values. Mobility is also
analysed using a per capita mobility indicator (Column 5). The positive and
significant coefficient of the per capita mobility confirms a pattern of Rt
that increases as the internal mobility measured on the total population
increases. The demographic controls of the total population and density,
presented in Columns 6 and 7, also have a positive effect on Rt in the first
wave. Finally, in Columns 8 and 9, we simultaneously estimate the effect of
the internal mobility, the level of urbanisation of the regions and the
population density. The coefficient of the internal mobility remains
significant and positive even when we include the other control variables.
However, the variables indicating the level of urbanisation lose their
significance while the variables on population size and density change sign.
Internal mobility appears to be a critical determinant of the rate of COVID-19
cases during the first wave, positively influencing the spread of the virus
possibly through increased social interactions. These results confirm that
great part of the territorial characteristics influencing the higher
epidemiological risk at the onset of the pandemic in urban regions can be
explained by the role of mobility. Table 2 examines the relationship between
Rt and different mobility patterns in the European regions in a similar way to
Table 1 but with data for the second wave (from August). The estimates show an
inversion of sign in respect the first wave with a positive association
between the virus spread and the rural and intermediate regions compared to
large cities. These results may reflect a behavioural response as well as more
severe containment measures in the most severely affected areas. In the second
wave, different mobility patterns are associated with decreased Rt values, as
presented in Column 2-5. The results thus indicate that the relationship
between mobility and the regional virus transmission has changed over time and
that shifts in mobility are used to control the epidemic. However, these
changes were not sufficient to prevent a second wave of infection in most of
the regions. The population demographic variable is significant and negative
on the virus spread. The models that simultaneously estimate the effect of
internal mobility and different regional characteristics confirm a negative
relationship between this mobility pattern and the virus transmission, as well
as a higher prevalence of the infection in rural regions compared to large
cities during the second wave. When interpreting the relationship between the
virus transmission and the local mobility during the second wave, it is
important to consider that mobility changes may respond endogenously to Rt
values, as mobility may be related to either individual or government-led
changes in preventive behaviour. These estimates may therefore be
underestimating the link between mobility and infections during the second
wave and a further investigation through an instrumental variables strategy
for our specifications would be preferable to address these endogeneity
concerns. The Figure 6 show the regression coefficients with shifts of
internal mobility of 3 weeks before and after the reference week of Rt. During
the first wave the relation of mobility on Rt is positive and peaks during the
same week. During the second wave the relation is negative and constantly
decreasing. The fact that the relation during the first wave is becoming
clearer towards the reference week is indicating that mobility is having an
effect on Rt. On the contrary during the second wave there is an inversion of
causality and mobility rather than influencing seem to be reacting to changes
in Rt moving in opposite directions. Intuitively this is in line with the
consideration that during the advanced stages of the pandemic mobility is
highly conditioned by restriction measures and lockdowns which are put in
place in correspondence with increases in Rt.

Figure 6: Regression coefficients with shifts of internal mobility of 3 weeks before and after the reference week of Rt. Table 1: Regression on Rt during the first wave (28 days since onset) | Rt first wave
---|---
| (1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9)
Intermediate |  $-$0.12∗∗ |  |  |  |  |  |  | 0.05 |  $-$0.04
| (0.05) |  |  |  |  |  |  | (0.05) | (0.06)
Rural |  $-$0.16∗∗∗ |  |  |  |  |  |  | 0.07 |  $-$0.06
| (0.06) |  |  |  |  |  |  | (0.06) | (0.07)
Internal |  |  1.84∗∗∗ |  |  |  |  |  |  1.99∗∗∗ |  4.54∗∗∗
|  | (0.22) |  |  |  |  |  | (0.25) | (0.40)
Inbound |  |  |  2.20∗∗∗ |  |  |  |  |  |
|  |  | (0.17) |  |  |  |  |  |
Outbound |  |  |  |  2.42∗∗∗ |  |  |  |  |
|  |  |  | (0.21) |  |  |  |  |
Internal pca |  |  |  |  |  0.27∗∗∗ |  |  |  |
|  |  |  |  | (0.02) |  |  |  |
Population |  |  |  |  |  |  0.05∗∗ |  |  |  $-$0.27∗∗∗
|  |  |  |  |  | (0.02) |  |  | (0.03)
Population density |  |  |  |  |  |  |  0.04∗∗∗ |  |  $-$0.05∗
|  |  |  |  |  |  | (0.02) |  | (0.02)
Observations | 3,572 | 3,572 | 3,570 | 3,571 | 3,572 | 3,572 | 4,550 | 3,572 | 3,572
R2 | 0.02 | 0.03 | 0.06 | 0.05 | 0.07 | 0.02 | 0.01 | 0.03 | 0.05
Adjusted R2 | 0.01 | 0.03 | 0.05 | 0.05 | 0.07 | 0.01 | 0.01 | 0.03 | 0.05
Note: |  ∗p$<$0.1; ∗∗p$<$0.05; ∗∗∗p$<$0.01
Table 2: Regression on Rt during the second wave (after August) | Rt second wave
---|---
| (1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9)
Intermediate |  0.004∗∗∗ |  |  |  |  |  |  | 0.001 | 0.002
| (0.001) |  |  |  |  |  |  | (0.002) | (0.002)
Rural |  0.01∗∗∗ |  |  |  |  |  |  |  0.01∗∗∗ |  0.01∗∗∗
| (0.002) |  |  |  |  |  |  | (0.002) | (0.002)
Internal |  |  $-$0.05∗∗∗ |  |  |  |  |  |  $-$0.03∗∗∗ |  $-$0.07∗∗∗
|  | (0.01) |  |  |  |  |  | (0.01) | (0.01)
Inbound |  |  |  $-$0.04∗∗∗ |  |  |  |  |  |
|  |  | (0.005) |  |  |  |  |  |
Outbound |  |  |  |  $-$0.04∗∗∗ |  |  |  |  |
|  |  |  | (0.01) |  |  |  |  |
Internal pca |  |  |  |  |  $-$0.004∗∗∗ |  |  |  |
|  |  |  |  | (0.001) |  |  |  |
Population |  |  |  |  |  |  $-$0.003∗∗∗ |  |  |  0.004∗∗∗
|  |  |  |  |  | (0.001) |  |  | (0.001)
Population density |  |  |  |  |  |  |  $-$0.0004 |  | 0.0002
|  |  |  |  |  |  | (0.0005) |  | (0.001)
Observations | 10,542 | 10,542 | 10,534 | 10,538 | 10,542 | 10,542 | 16,204 | 10,542 | 10,542
R2 | 0.04 | 0.04 | 0.04 | 0.04 | 0.04 | 0.03 | 0.05 | 0.04 | 0.04
Adjusted R2 | 0.04 | 0.04 | 0.04 | 0.04 | 0.04 | 0.03 | 0.05 | 0.04 | 0.04
Note: |  ∗p$<$0.1; ∗∗p$<$0.05; ∗∗∗p$<$0.01

## Conclusion

In this article we analysed the territorial differences in the onset and
spread of COVID-19 and the associated excess mortality, across the European
NUTS3 regions and US counties during the first and second COVID-19 wave.
During the first wave, the COVID-19 pandemic arrived earlier, recorded higher
Rt values and had a higher impact in terms of excess mortality in urban
regions compared to the intermediate and the rural ones. In the first wave,
mobility influenced the spread of COVID-19, since the higher mobility of urban
regions is explaining entirely the differences between the three groups of
regions. The fact that these effects are more difficult to recognise in later
stages of the pandemic can be tentatively explained by the widespread of the
infection, the implementation of restriction measures which invert the
causality between mobility and Rt, often applied on a territorial basis, and
the more complex mobility patterns experienced during the summer period.

Our findings are in line with previous studies identifying the role of
mobility on virus spread in the early stages of the pandemic. To our
knowledge, our research is unique in providing a broad geographic coverage and
a high level of geographical detail, and in examining the role of regional
mobility for the spread of COVID-19 through a unique data set on regional
mobility derived from mobile phone data. In terms of policy implication, our
research contributes to a better understanding of territorial characteristics
of the spread of COVID-19, which is critical for designing effective public
health policy responses, often decided at regional level.

## Acknowledgments

The authors acknowledge the support of European MNOs (among which 3 Group -
part of CK Hutchison, A1 Telekom Austria Group, Altice Portugal, Deutsche
Telekom, Orange, Proximus, TIM Telecom Italia, Tele2, Telefonica, Telenor,
Telia Company and Vodafone) in providing access to aggregate and anonymised
data. The authors would also like to acknowledge the GSMA444GSMA is the GSM
Association of Mobile Network Operators., colleagues from Eurostat555Eurostat
is the Statistical Office of the European Union. and ECDC666ECDC: European
Centre for Disease Prevention and Control. An agency of the European Union.
for their input in drafting the data request.
Finally, the authors would also like to acknowledge the support from JRC
colleagues, and in particular the E3 Unit, for setting up a secure environment
to host and process of the data provided by MNOs, as well as the E6 Unit (the
“Dynamic Data Hub team”) for their valuable support in setting up the
database.

## References

  * [1] Andrew J Stier, Marc G Berman, and Luis M A Bettencourt.  Covid-19 attack rate increases with city size, 2020.
  * [2] Haroldo V. Ribeiro, Andre S. Sunahara, Jack Sutton, Matjaž Perc, and Quentin S. Hanley.  City size and the spreading of covid-19 in brazil.  PLOS ONE, 15(9):1–12, 09 2020.
  * [3] Samuel Heroy.  Metropolitan-scale COVID-19 outbreaks: how similar are they?
  * [4] Mattia Mazzoli, David Mateo, Alberto Hernando, and Sandro Meloni.  Effects of mobility and multi-seeding on the propagation of the covid-19 in spain, 2020.
  * [5] Michiel Gerritse.  Cities and COVID-19 infections: Population density, transmission speeds and sheltering responses, 2020.
  * [6] Matteo Chinazzi, Jessica T. Davis, Marco Ajelli, Corrado Gioannini, Maria Litvinova, Stefano Merler, Ana Pastore y Piontti, Kunpeng Mu, Luca Rossi, Kaiyuan Sun, Cécile Viboud, Xinyue Xiong, Hongjie Yu, M. Elizabeth Halloran, Ira M. Longini, and Alessandro Vespignani.  The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak.
  * [7] Milena Almagro and Angelo Orane-Hutchinson.  The determinants of the differential exposure to COVID-19 in new york city and their evolution over time.
  * [8] Mauro Agnoletti, Simone Manganelli, and Francesco Piras.  Covid-19 and rural landscape: The case of italy.  Landscape and Urban Planning, 204:103955, 2020.
  * [9] Chrstian Bayer and Moritz Kuhn.  Intergenerational ties and case fatality rates, March 2020.
  * [10] Marianna Belloc, Paolo Buoananno, and Francesco Drago.  Cross-country correlation analysis for research on COVID-19 | VOX, CEPR Policy Portal, 2020.
  * [11] J. Wallinga.  Different epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures.
  * [12] Zhanwei Du, Xiaoke Xu, Ye Wu, Lin Wang, Benjamin J. Cowling, and Lauren Ancel Meyers.  Serial interval of COVID-19 among publicly reported confirmed cases.
  * [13] EDPB.  Guidelines 04/2020 on the use of location data and contact tracing tools in the context of the covid-19 outbreak, 2020.  Available at https://edpb.europa.eu/our-work-tools/our-documents/linee-guida/guidelines-042020-use-location-data-and-contact-tracing_en.
  * [14] M. Mamei, N. Bicocchi, M. Lippi, S. Mariani, and F. Zambonelli.  Evaluating origin–destination matrices obtained from cdr data.  Sensors, 19:1440, 2019.
  * [15] M. Fekih, T. Bellemans, Z. Smoreda, P. Bonnel, A. Furno, and S. Galland.  A data-driven approach for origin–destination matrix construction from cellular network signalling data: a case study of lyon region (france), 2020\.
  * [16] European Commission.  Communication from the commission to the european parliament and the council staying safe from covid-19 during winter, com/2020/786, 2020.
https://eur-lex.europa.eu/legal-content/EN/ALL/?uri=COM:2020:786:FIN.

  * [17] Carlos Santamaria, Francesco Sermi, Spyridon Spyratos, Stefano Maria Iacus, Alessandro Annunziato, Dario Tarchi, and Michele Vespe.  Measuring the impact of covid-19 confinement measures on human mobility using mobile positioning data. a european regional analysis.  Safety Science, 132:104925, 2020.
  * [18] Krzysztof Bartoszek, Emanuele Guidotti, Stefano Maria Iacus, and Marcin Okrój.  Are official confirmed cases and fatalities counts good enough to study the covid-19 pandemic dynamics? a critical assessment through the case of italy.  Nonlinear Dynamics, 101(3):1951–1979, Aug 2020.


∎

11institutetext:  University of Wroclaw, Plac M. Borna 9, 50-204 Wroclaw,
Poland
11email<EMAIL_ADDRESS>
# QCD phase transition drives supernova explosion of a very massive star

Dependence on metallicity of the progenitor star

Tobias Fischer

(Received: 28 May 2021 / Accepted: 17 August 2021
©The Author(s) 2021
Communicated by Laura Tolos )

###### Abstract

The nature of core-collapse supernova (SN) explosions is yet incompletely
understood. The present article revisits the scenario in which the release of
latent heat due to a first-order phase transition, from normal nuclear matter
to the quark–gluon plasma, liberates the necessary energy to explain observed
SN explosions. Here, the role of the metallicity of the stellar progenitor is
investigated, comparing a solar metallicity and a low-metallicity case, both
having a zero-age main sequence (ZAMS) mass of 75 M⊙. It is found that low-
metallicity models belong exclusively to the failed SN branch, featuring the
formation of black holes without explosions. It excludes this class of massive
star explosions as possible site for the nucleosynthesis of heavy elements at
extremely low metallicity, usually associated with the early universe.

##  1 Introduction

All stars more massive than about 8 times the mass of the sun (M⊙) end their
life as a core-collapse SN. These events belong to the most powerful outbursts
in the modern universe, where most of the energy is carried in form of
neutrinos. SN explosions are associated with the revival of the stalled bounce
shock. The latter forms when the collapsing stellar of a massive star reaches
central densities slightly in excess of normal nuclear matter density
($\rho_{\rm sat}\simeq 2.5\times 10^{14}$ g cm-3, equivalent $0.15$ fm-3 in
nuclear units), when the short-range repulsive nuclear force balances gravity
such that the collapsing core bounces back. The extreme conditions reached at
the SN interior, supersaturation densities and temperatures on the order of
10–50 MeV as well as large neutron excess given by the electron abundance of
$Y_{e}=0.01-0.3$, make the core-collapse SN ideal high-energy nuclear and
particle physics laboratories Fischer:2017 .

Several SN explosions scenarios, i.e. the liberation of energy from the
nascent compact central hot and dense object—the proto-neutron star (PNS)—to
the bounce shock have been proposed; the magneto-rotational mechanism
LeBlanc:1970kg  and the presently considered standard neutrino-heating
mechanism Bethe:1985ux . The latter, however, requires multi-dimensional
simulations where the neutrino heating efficiency increases in the presence of
convection and hydrodynamics instabilities, to obtain massive star explosions
for iron-core progenitors more massive than about 9 M⊙ Janka:2007hi ;
Mirizzi:2016 . The exceptions are the slightly less massive progenitors in the
zero-age main sequence mass range of 8–9 M⊙ Nomoto:1987 ; Kitaura:2006 ;
Melson:2015 , for which neutrino-driven SN explosions can be obtained in
spherical symmetric simulations.

The present article focuses on an additional SN explosion mechanism. It is
related to the energy, more precisely latent heat release, due to a first-
order phase transition from normal nuclear matter to the quark–gluon plasma at
high density Takahara:1988yd ; Gentile:1993ma ; Sagert:2008ka . First
principle calculations of Quantum Chromodynamics (QCD), the theory of strong
interactions with quarks and gluons as the degrees of freedom, are only
possible at vanishing baryon densities, where the QCD equations are solved
numerically being mapped onto a space-time lattice. These computations predict
a smooth cross-over transition at a pseudo-critical temperature in the range
of 150–160 MeV Bazavov:2014 ; Borsanyi:2014 ; Bazavov:2019 . At finite and
high baryon density, phenomenological quark matter models have commonly been
employed in astrophysical studies of compact stellar objects.

The classical but simplistic quark-bag model EOS Farhi:1984qu  has been the
foundation in the core-collapse SN study of Ref. Sagert:2008ka . It is based
on the equation of state (EOS) of non-interacting quarks, shifted by a
constant bag pressure. One of the most striking problems for this class of
quark-matter hybrid EOS in applications to astrophysical studies, featuring
hot and dense matter, is the difficulty to obtain maximum neutron star masses
compatible with the current observational constraint of about 2 M⊙
Antoniadis:2013 ; Fonseca:2016 ; Cromartie:2020NatAs ; Fonseca:2021 . This
caveat has been overcome recently, with the introduction of repulsive quark
matter interactions Benic:2015 ; Klaehn:2015 ; Kaltenborn:2017 , providing
sufficient stiffness with increasing density. The extension to finite
temperatures and arbitrary isospin asymmetry is seemingly trivial for the
class of quark-bag model EOS Klaehn:2017 . It represents a challenge for
microscopic quark-matter EOS, such as the Nambu–Jona-Lasinio models NJL:1961 ;
Buballa:2003qv ; Klaehn:2013  as well as Schwinger–Dyson models Roberts:2000aa
; Chen:2008zr ; Chen:2011my ; Bashir:2012fs ; Cloet:2013jya ; Chen:2015mda .
These phenomenological approaches are applicable to the non-perturbative
regime of QCD, however, they are built on a certain truncation of QCD. Most
importantly, they are lacking confinement. On the other hand, perturbative QCD
is valid only near the limit of asymptotic freedom, where quarks are no more
strongly coupled Kurkela:2009gj ; Kurkela:2014vha ; Kurkela:2020NatPh .

The recently developed phenomenological string-flip quark-matter EOS is the
foundation of the present work. It features not only a first-order phase
transition from normal nuclear matter, it is also built on a mechanism for
(de)confinement Kaltenborn:2017 . It has been extended to arbitrary isospin
asymmetry and finite temperatures Fischer:2018 ; Bastian:2021 . With model
parameters selected such that the conditions for the hadron–quark phase
transition to occur at the interior of massive compact stars, i.e. densities
for the onset of quark-matter on the order of 2–5 $\times\rho_{\rm sat}$, one
is led to the conclusions that the associated SN explosion mechanism is likely
to operate for the high-mass end of core-collapse SN progenitors, associated
with the ZAMS mass of $>30$ M⊙ Fischer:2018 , whereas less massive progenitor
stars would be subject to other SN explosion scenarios. Note that also the
failed SN branch, in which a black hole forms instead of an explosion, is
considered within the class of hadron–quark phase transition SN. It depends on
the compactness, more precisely, on the mass enclosed inside the PNS at the
onset of the hadron–quark phase transition. If the latter exceeds the maximum
mass of the hadron–quark hybrid EOS, then an explosion is likely to fail.

The string-flip class of hadron–quark hybrid EOS has also been the foundation
for studies of binary neutron star mergers Bauswein:2019PhRvL ;
Bauswein:2020PhRvL ; Blacker:2020 . One of the key goals has been the
identification of a first-order phase transition through possibly observable
signatures, in the neutrino signal from SN and gravitational waves
OConnor:2020 , the latter also in the context of binary neutron star mergers
Most:2019PhRvL .

One important aspect of the SN explosion scenario driven by the QCE phase
transition is the role of metallicity111Astronomers define as metals all
elements heavier than helium., which has not yet been addressed. It strongly
affects the evolution of massive stars through mass loss at the star’s surface
occurring during the different nuclear burning stages, and hence the core
structure at the onset of core collapse. This article reports the SN explosion
of a very massive ZAMS progenitor star of 75 M⊙ of solar metallicity, driven
by the hadron–quark phase transition, in comparison with those launched from a
75 M⊙ of extremely low metallicity.

The manuscript is organised as follows: In Sect. 2 the hadron–quark matter
hybrid EOS is revisited. Section 3 reviews briefly the SN model employed in
this study together with a discussion of the stellar progenitors. The SN
simulations are discussed in Sect. 4, and the manuscript closes with the
summary and conclusions in Sect. 5.

##  2 Dense matter equation of state with phase transition

For the present study the reference hadronic DD2F EOS is implemented. It is a
nuclear relativistic mean-field (RMF) model with density-dependent
meson–nucleon couplings Typel:2005 ; Typel:2009sy , which has been modified
from the original DD2 RMF parametrization in order to be consistent with the
flow constraint Danielewicz:2002Sci . As a consequence, DD2F is slightly
softer than DD2 at densities in excess of saturation density. At low
densities, below the saturation density, and at temperatures below about 15–20
MeV, nuclear clusters are taken into account based on the nuclear statistical
equilibrium model of Ref. Hempel:2009mc , with several 1000 nuclear species
with tabulated and partly calculated nuclear masses. It implements the
dissolving of all clusters through a geometric excluded volume approach, which
can be adjusted to calculations of in-medium nuclear properties in order to
overcome the limitations of the statistical model Fischer:2020c .

Table 1: DD2F-RDF 1.2 hybrid EOS properties, onset density for the phase transition ($\rho_{\rm onset}$), the density for reaching the pure quark-matter phase ($\rho_{\rm final}$), as well as inset neutron star mass ($M_{\rm onset}$) and maximum mass ($M_{\rm max}$) together with the central density of the maximum mass configuration ($\rho_{\rm central}|_{M_{\rm max}}$), for $T=0$ and a constant entropy per particle of $s=3~{}k_{\rm B}$. DD2F-RDF 1.2 |  |  | $\rho_{\rm onset}$ |  |  | $\rho_{\rm final}$ |  |  | $M_{\rm onset}$ |  |  | $M_{\rm max}^{\rm RDF~{}1.2}$ |  |  | $\left.\rho_{\rm central}\right|_{M_{\rm max}}$ |  |  | $M_{\rm max}^{\rm DD2F}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  | $[\rho_{\rm sat}]$ |  |  | $[\rho_{\rm sat}]$ |  |  | $[{\rm M}_{\odot}]$ |  |  | $[{\rm M}_{\odot}]$ |  |  | $[\rho_{\rm sat}]$ |  |  | $[\rho_{\rm sat}]$
$T=0$ |  |  | 3.05 |  |  | 3.60 |  |  | 1.37 |  |  | 2.16 |  |  | 7.14 |  |  | 2.09
$s=3~{}k_{\rm B}$ |  |  | 1.15 |  |  | 3.25 |  |  | 1.44 |  |  | 2.07 |  |  | 7.22 |  |  | 2.17

The DD2F hadronic EOS has been extended to take into account a phase
transition to the quark–gluon plasma. Therefore, the string-flip EOS for quark
matter is employed Kaltenborn:2017 . It is based on the relativistic density
functional (RDF) approach to quark matter. The present work uses the
parametrization denoted as RDF 1.2 from the comprehensive catalog of Ref.
Bastian:2021 222All EOS of the present study are available at the CompOSE data
base under https://compose.obspm.fr Typel:2014 .. For the phase-transition
construction between the DD2F hadronic EOS and the string-flip EOS, Maxwell’s
condition of phase equilibrium has been applied, i.e kinetic equilibrium for
the pressures with respect to the baryon chemical potential, for fixed
temperature, $T$, and charge chemical potential, $\mu_{Q}$. The latter is
defined as the difference between proton and neutron chemical potentials,
$\mu_{p}$ and $\mu_{n}$, $\mu_{Q}=\mu_{p}-\mu_{n}$ in the hadronic phase, and
as the difference between up and down chemical potentials, $\mu_{u}$ and
$\mu_{d}$, $\mu_{Q}=\mu_{u}-\mu_{d}$ in the quark matter phase. Maxwell’s
construction results in a jump condition for all thermodynamic quantities,
from the pure hadronic phase to the pure quark matter phase. However, for
practical purposes in applications of astrophysical simulations, data are
needed in the coexistence region in between. Therefore, a quark-volume
fraction has been employed, with a linear dependence on the baryon density.
Further details can be found in Ref. Bastian:2021 .

(a) $T=0$

(b) $s=3~{}k_{\rm B}$

Figure 1: EOS under the condition of $\beta$-equilibrium, showing the pressure
$P$, the proton abundance $Y_{p}$ and the temperature $T$ (only for finite
entropy), with respect to the restmass density $\rho$, for $T=0$ in graph (a)
and constant entropy per particle of $s=3~{}k_{\rm B}$ in graph (b), as well
as the corresponding mass–radius relations, comparing the hadronic DD2F
reference model (red dached lines) and SF2 (solid gray lines). The thin
vertical solid and dash-dotted lines mark the location for the onset of the
phase transition and when reaching the pure quark-matter phase, respectively,
whereas the horizontal line in the mass–radius relation graphs marks the
corresponding onset mass, together with the present maximum neutron star
constraint of $2.08\pm 0.07$ M⊙ from the measurement of the pulsar PSR
J0740+6620 Antoniadis:2013 ; Fonseca:2021 .

The class of string-flip EOS is build on a confining mechanism through
divergent particle masses towards vanishing density. It still remains to be
shown in how far it matches the (dynamical) restoration of chiral symmetry. In
addition, in order to obtain sufficient stiffness at high densities, repulsive
interactions are taken into account. The original string-flip EOS includes
besides linear terms higher-order repulsive interactions, which are omitted
here in the RDF 1.2 version, such that the problem of causality—when the speed
of sound approaches the speed of light—does not occur in the density range
encountered at the interior of compact stellar objects. The string-flip model
has been extended to include an isospin mean field, i.e. the vector–isovector
interaction channel, which are comparable to the $\rho$-meson interactions in
nuclear RMF-type models. The coupling can be adjusted to a smooth behavior of
the symmetry energy across the phase transition boundary (a detailed
discussion of this aspect is given in Sect. II C3. of Ref. Bastian:2021 ).

Selected properties of the resulting DD2F-RMF 1.2 hybrid EOS are listed in
Table 1, for zero temperature and finite entropy per particle of $s=3~{}k_{\rm
B}$. Particularly interesting is the reduction of the onset density for the
hadron–quark phase transition, $\rho_{\rm onset}$, with increasing entropy,
viz. temperature, while the final density for reaching the pure quark matter
phase, $\rho_{\rm final}$, remains nearly constant for the temperatures of
relevance here on the order of $T=0-50$ MeV. These aspects will become
relevant in the SN simulations which will be discussed below.

The DD2F-RDF 1.2 hadron–quark hybrid EOS in $\beta$-equilibrium is illustrated
in Fig. 1a for $T=0$ and in Fig. 1b for finite entropy per particle of
$s=3\,k_{\rm B}$, in comparison to the hadronic DD2F reference EOS. Note that
the pressure, $P$, contains all contributions, hadrons and/or quarks as well
as electrons, positrons and photons. The last two are added following Ref.
Timmes:1999 . The vertical lines mark $\rho_{\rm onset}$ (solid lines) and
$\rho_{\rm final}$ (dash-dotted lines), where it is interesting to note that
while for $T=0$ the pressure is nearly constant in the hadron–quark mixed
region, for finite entropy case there is a finite slope. The latter aspect is
attributed to the fact that there are different temperatures for the phase
transition, $T(\rho_{\rm onset})$, and for reaching the pure quark-matter
phase, $T(\rho_{\rm final})$. In order to preserve the entropy per particle,
the temperature decreases throughout the mixed phase, for entropies less than
about $s\simeq 6-7~{}k_{\rm B}$, while the temperature increases for higher
entropies. This feature is already illustrated in Fig. 1 in Ref. Fischer:2018
(see also the solid green lines in Fig. 7 which will be discussed below).

Figure 2: Comparison of the progenitors with a ZAMS masses of 75 M⊙ but
different metallicities, from Ref. Woosley:2002zz , showing the density $\rho$
and entropy per particle $s$ (left panels) and the abundance of helium denoted
as He (dotted lines), carbon and oxygen denoted as C+O (thin dash-dotted
lines) as well as silicon and sulphur denoted as Si+S (dashed lines) and iron-
group nuclei collectively denoted as ’Fe’ (thick solid lines), with respect to
the mass coordinate. In addition, for the low-metallicity progenitor, 56Ni is
shown as thin solid line. (a) Solar metallicity
(b) $Z=10^{-4}\times Z_{\odot}$

It is also interesting to note that the aforementioned Maxwell construction
for the phase transition results in gradients for various quantities, e.g.,
the proton fraction, in general the charge fraction across the coexistence
region, in order to conserve the charge chemical potential (see the middle
panels in Fig. 1). This behavior is imposed by the linear dependence on
density in the coexistence region. For a real Maxwell construction, both
pressure and charge fraction would jump from $\rho_{\rm onset}$ to $\rho_{\rm
final}$. Note that in quark matter we construct proton and neutron abundances,
$Y_{p}$ and $Y_{n}$, from the abundances of up- and down-quark, $Y_{u}$ and
$Y_{d}$, as follows: $Y_{p}=2Y_{u}+Y_{d}$ and $Y_{n}=Y_{u}+2Y_{d}$.

The corresponding mass–radius relations are shown the right panels of Fig. 1.
For the DD2F-RDF 1.2 parameters, at $T=0$, the maximum mass of the hybrid star
exceeds the maximum mass of the DD2F reference hadronic EOS; the exact values
are given in Table 1. The onset mass for the hybrid branch, i.e. the third
family of compact stellar objects, is at $M_{\rm onset}=1.37$ M⊙. These
hydrostatic properties change when considering hot configurations, in
particular at entropy of $s=3~{}k_{\rm B}$, which represents a canonical value
at the PNS interior. Note that it has been found that the universal relations
YagiYunes:2013Sci ; YagiYunes:2013PhRvD ; Maselli:2013  hold only for constant
but low values of the entropy Raduta:2020 ; Sedrakian:2021 , since the entropy
is a conserved quantity. Since the latter is also the case for SN, it is
suitable to consider the situation of constant entropy, rather than constant
temperature, in order to understand what happens in simulations of core-
collapse SN regarding thermal effects. The first observation is that the
coexistence region extends due to the different behaviors of hadronic and
quark matter EOS with increasing temperature (see Fig. 1b), besides the
significantly larger gradients between $\rho_{\rm onset}$ (vertical solid
lines) and $\rho_{\rm final}$ (vertical dash-dotted lines). Secondly, the
radii of the hydrostatic configurations become significantly larger, on the
order of 20–40 km. Secondly, while the maximum mass of the DD2F reference
hadronic EOS increases, the maximum mass of the DD2F-RDF 1.2 hybrid
configuration reduces. Moreover, an extended region of gravitational
instability develops after the onset of quark matter at the transition from
the stable hadronic to the stable quark matter branch Hempel:2016 . This twin
phenomenon has long been studied in the context of cold neutron stars, i.e.
the existence of two compact stars with the same gravitational mass but
different radii and hence belonging to different families of compact stars
(cf. Ref. Benic:2015  and references therein). It is interesting to note that
the twin phenomenon is only present in case of finite, and in particular high
temperatures, illustrated by the temperature jump around $T\simeq 50$ MeV at
the phase transition (middle panel in Fig. 1b).

Figure 3: Post-bounce evolution of the central and maximum temperatures in
graph (a), $T_{\rm central}$ (thick lines) and $T_{\rm max}$ (thin lines), the
central density $\rho_{\rm central}$ in graph (b) and the PNS mass $M_{\rm
PNS}$ in graph (c) as well as in graph (d) the PNS radius $R_{\rm PNS}$ (thick
lines) and the SN shock radius $R_{\rm shock}$ (thin lines), including the
second shock wave due to the hadron–quark phase transition (red dash-dotted
lines), as well as in graph (e) the diagnostic energy $E_{\rm diagnostic}$,
comparing the SN simulations of the 75 M⊙ progenitors of solar metallicity
(solid gray lines) and $10^{-4}$ of solar metallicity (blue dashed lines). The
thin horizontal lines in graph (b) mark the onset densities for the
hadron–quark phase transition for the corresponding temperatures and a
represented proton fraction of $Y_{p}=0.3$.

##  3 Supernova model and stellar progenitors

The present numerical study of core-collapse SN is based on AGILE-BOLTZTRAN
(see Ref. Liebendoerfer:2004  and references therein), which is a general
relativistic neutrino-radiation hydrodynamics model in spherical symmetry
featuring accurate 6-species Boltzmann neutrino transport Fischer:2020d . The
set of non-muonic weak interactions considered here can be found in Table 1 of
Ref. Fischer:2020a , with the references given accordingly, and the muonic
weak processes are introduced and discussed in Refs. Guo:2020 ; Fischer:2020d
. AGILE-BOLTZTRAN has an adaptive baryon mass mesh refinement
Liebendoerfer:2002 ; Fischer:2009af . In the present simulations 207 radial
grid points are used, and the neutrino phase space is discretized following
Refs. Bruenn:1985en ; Mezzacappa:1993gm , with 6 neutrino lateral angle bins,
$\cos\theta\in\\{-1,+1\\}$, and 36 neutrino energy bins, from 0.5–300 MeV.

For the present investigation of the phase transition in core-collapse SN as
explosion mechanism, one of the yet-incompletely addressed questions is the
dependence of the stellar progenitor model. In particular, what is the maximum
progenitor mass that the hadron–quark phase transition is able to drive the
explosion, and how does it depend on metallicity, denoted as $Z$? Therefore,
the most massive stellar progenitor is selected from the comprehensive stellar
evolution series of Ref. Woosley:2002zz , with a ZAMS mass of 75 M⊙, while a
stellar model will ZAMS mass of 50 M⊙ has been explored previously
Umeda:2007wk ; Fischer:2018 .

Radial profiles, as a function of the enclosed mass, for selected quantities
are shown in Fig. 2a for the solar metallicity model. The current value of the
solar surface metallicity is $Z_{\odot}=0.0187\pm 0.0009$, Asplund:2009 ;
Asplund:2021 . This features a rather compact stellar core at the onset of
stellar collapse, with an iron-core mass of about 1.5 M⊙ and with central
densities in excess of $5\times 10^{9}$ g cm-3, as well as a narrow
silicon–sulfur shell above the iron core (dashed line in the right panel of
Fig. 2a. Typical for solar-metallicity stars is a large mass-loss rate, which
results in a sharp transition from the stellar core, i.e. carbon-oxygen shell
(dash-dotted line in the right panel of Fig. 2a) and the transition to the
helium-rich envelope, with a sharp density drop over many orders of magnitude.
The resulting total enclosed mass is about 6.36 M⊙ at the onset of stellar
core collapse Woosley:2002zz , i.e. about 68.64 M⊙ has been lost during the
previous nuclear burning stages, mostly during hydrogen and helium burning
Woosley:2002zz . The later collision of the SN ejecta with the material
ejected pre-collapse may power the brightest SN explosions known as super-
luminous SN Fischer:2018 .

The situation is very different for stars with lower metallicity. Figure 2b
shows the same radial profiles for the same 75 M⊙ ZAMS mass but at a
metallicity of $10^{-4}$ of the solar value, from the same stellar evolution
series Woosley:2002zz . This stellar model has experienced nearly zero mass
loss. It it related with the photon opacity, which scales with the number of
protons and is hence substantially higher for ions. The lack of efficient mass
loss leaves a nearly as massive progenitor star, with a total mass at the
onset of stellar collapse of 75 M⊙ Woosley:2002zz , featuring a central
density of about $10^{9}$ g cm-3 and a higher central entropy per baryon of
about $s\simeq 2~{}k_{\rm B}$, in comparison to the solar metalicty progenitor
with $s\simeq 1.5~{}k_{\rm B}$ (see Fig. 2). Note also that the nascent iron
core contains nearly 3 M⊙ (the abundance of 56Ni is shown as thin solid line
in the right panel of Fig. 2b).

##  4 Supernova simulations

In this section, results will be discussed from simulations of core-collapse
SN, launched from the stellar progenitors introduced in Sect. 3, with a ZAMS
mass of 75 M⊙ of solar metallicity and $10^{-4}$ of the solar metallicity
Woosley:2002zz . Note that the SN simulations of the solar metallicity
progenitor contains the entire remaining mass of 6.36 M⊙. Both SN simulations
employ the DD2F-RDF-1.2 hadron–quark hybrid EOS, reviewed in Sect. 2.

Figure 4: Post-bounce evolution of the mass accretion rate, $\dot{M}$, sampled
at the PNS surface, comparing the SN simulations of the 75 M⊙ progenitors of
solar metallicity and $10^{-4}$ solar metallicity. The inlay shows a zoom in
the later post-bounce evolution for the solar metallicity model.

Both SN simulations proceed qualitatively similarly through the stellar core
collapse, featuring slightly higher central temperatures for the low-
metallicity model as well as lower central density, than the solar-metallicity
one, as illustrated in graphs (a) and (b) in Fig. 3. However, of particular
interest is the post-bounce evolution, during which differences grow
substantially between the two models. These, in turn, can be related to the
mass accretion rates, shown in Fig 4, sampled at the PNS surface. The mass
accretion rate, $\dot{M}=4\pi r^{2}\rho v$, with rest-mass density $\rho$ and
velocity $v$ at the PNS surface, is significantly larger for the low-
metallicity case as a direct consequence of the higher density of the far more
massive iron core and of the extended silicon-sulphur layer above the iron
core (see the left panels in Fig 2). The mass enclosed inside the PNS grows
very rapidly shortly after core bounce, reaching shortly 2 M⊙, as illustrated
in Fig. 3c. It leads to a rapid rise of the central density and of the central
temperature, and in particular to maximum temperatures at the PNS interior,
denoted as $T_{\rm max}$, in excess of 30–50 MeV. The very massive PNS gives
rise to high neutrino luminosities, on the order of $L_{\nu}\simeq 1-2\times
10^{53}$ erg s-1 (see Fig. 5) during the evolution after the
$\nu_{e}$-deleptonization burst has been launched, about 50 ms post-bounce.
This is nearly one order of magnitude higher than the one obtained in normal
SN simulations after the high-density silicon-sulfur layer has fallen onto the
SN shock (c.f. Refs. Lentz:2015 ; Mueller:2015 ; Mirizzi:2016 ; Mueller:2017 ;
Kotake:2018 ; Burrows:2019MNRAS ; Bollig:2020arXiv201010506B  and references
therein). Exactly this phenomenon has not been the case for the low metallicty
model here, for the entire time of the post-bounce evolution the neutrino
luminosities are powered by extraordinary high mass accretion rates, being
approximated as accretion luminosity for $\nu$ and $\bar{\nu}_{e}$ as follows,

$\displaystyle L_{\nu_{e}}\propto 2\times 10^{53}~{}{\rm
erg~{}s}^{-1}\left(\frac{M_{\rm PNS}}{2~{}{\rm
M}_{\odot}}\right)\left(\frac{\dot{M}}{4~{}{\rm
M}_{\odot}{/s}}\right)\left(\frac{10^{2}~{}{\rm km}}{R_{\nu}}\right)~{}.$

The values for the PNS mass, $M_{\rm PNS}$, is larger than for the solar
metallicity case by about one third, and the mass accretion rate is larger my
more than a factor of two (see Figs. 3 and 4). The high neutrino luminosities
are the reflection of enhanced neutrino heating rates, which, in turn, result
in larger neutrinosphere radii and also a systematically larger shock radius,
as illustrated in Fig. 3d. Such phenomenon has been reported based on the SN
simulation of a zero-metallicity stellar progenitor with ZAMS mass of 70 M⊙
Kotake:2018MNRAS , which, however, eventually lead to a failed SN explosion
and the subsequent black hole formation at about 253 ms post bounce for this
low-metallicity model. It is associated with the sudden rise of the maximum
temperature and central density (see graphs (a) and (b) in Fig. 3) and the
abrupt cease of the mass accretion rate (see Fig. 4).

Here, for the low-metallicity 75 M⊙ progenitor, the high mass accretion rate
during the early post-bounce evolution results in a critical behavior when the
conditions are reached for the hadron–quark phase transition. These are
illustrated in Fig. 3 (b) through thin horizontal lines, showing the critical
density for the onset of the phase transition, $\rho_{\rm
onset}=2\times\rho_{\rm sat}$, for the corresponding temperature, shown in
Fig. 3 (a), and a representative proton fraction of $Y_{p}=0.3$. Once the
phase transition takes place, the central quark-matter volume fraction starts
to rise rapidly for this model and the PNS becomes gravitationally unstable
due to the reduced pressure gradient in the co-existence region, between the
two stable hadronic and quark matter phases, as was discussed in Sect. 2. The
PNS collapse proceeds supersonically. It halts when a sufficient mass at the
PNS interior is converted into the quark matter phase, where the EOS stiffens
again. As a consequence, a second shock wave forms which propagates to
increasingly larger radii where it turns into an accretion front. However,
before the accretion front can accelerate along the decreasing density
gradient towards the PNS surface, the PNS collapses behind this accretion
front. This is related to the fact that the enclosed PNS mass of $M\simeq 2.6$
M⊙ (see graph (c) of Fig. 3) exceeds the maximum mass of the DD2F-RDF-1.2
hybrid EOS (see Table 1). This phenomenon has already been reported in Ref.
Fischer:2018  for the most compact stellar progenitor with ZAMS mass of 25 M⊙
explored in that study.

Figure 5: Evolution of the neutrino luminosities (top panel) and average
energies (bottom panel) for all neutrino flavors, $\nu_{e}$, $\bar{\nu}_{e}$
and $nu_{x}$ collectively denoted for all heavy lepton flavor neutrinos and
$\bar{\nu}_{x}$ for their antineutrinos respectively, sampled in the co-moving
frame of reference at 500 km, for the SN simulation of the low-metallicity 75
M⊙ progenitor with $Z=10^{-4}\times Z_{\odot}$.

Figure 6: The same as Fig. 6 but for the SN simulations for the solar-
metallicity 75 M⊙ progenitor.

The SN post-bounce evolution proceeds differently for the progenitor of solar
metallicity, featuring a substantially lower mass accretion during the entire
post-bounce evolution (solid line in Fig. 4). It results in an extended post-
bounce mass accretion phase which lasts for nearly 1.0 second, until the
central density reaches a value of about twice saturation density, in
comparison to the low metallicity case when the same central density is
reached already at about 253 ms post bounce (see graph (b) in Fig. 3). The low
mass accretion rate for the solar metallicity case, with $\dot{M}\leq 0.1$ M⊙
s-1 already after about 500 ms post bounce, leads to a long timescale for the
PNS mass to grow. In fact the PNS mass never reaches the maximum mass of the
DD2F-RDF-1.2 hybrid EOS, it stays around $M~{}{\rm PNS}\simeq 1.8$ M⊙ (see
graph (c) in Fig. 3). Moreover, the SN shock contracts below $R_{\rm
shock}<100-80$ km once the silicon-sulfur layer has been fallen onto the
shock, and the neutrino luminosities decrease to values on the order of
$L_{\nu_{x}}\simeq 1\times 10^{51}$ erg s-1 ($\nu_{x}$ denote collectively all
heavy lepton flavors, here $\nu_{x}=\nu_{\mu}$) and
$L_{(\bar{\nu}_{e})\nu_{e}}\simeq 2\times 10^{51}$ erg s-1 (see Fig. 6).

Figure 7: Evolution of the SN simulation of the 75 M⊙ progenitor, of solar
metallicity, in the $T-\rho$ phase diagram where the color-coding is due to
the electron fraction, $Y_{e}$. Indications for the phase transition are set
as onset of the quark matter phase, $\rho_{\rm onset}$, and reaching the pure
quark matter phase, $\rho_{\rm final}$. The evolution of the central
temperature and density are also shown (thick solid black line), together with
curves of constant entropy per particle for selected values of $s=1-7~{}k_{\rm
B}$ at fixed proton fraction of $Y_{p}=0.3$ (solid green lines).

Subsequently, the conditions for the hadron–quark phase transition are reached
only after about 1.795 s post bounce for the solar-metallicity case,
corresponding to an onset density for the phase transition of $\rho_{\rm
onset}=5.6\times 10^{14}$ g cm-3, indicated by the thin horizontal line in
graph (b) of Fig. 3, for the central temperature of $T_{\rm central}\simeq 30$
MeV, shown in graph (a). Note that also for this model, maximum temperatures
reach as high as $T_{\rm max}\simeq 40-50$ MeV. However, unlike in the case
for the low-metallicity model discussed above, the second shock wave, which is
associated with the PNS collapse and transition into the hybrid branch,
expands to continuously larger radii along the decreasing density at the PNS
surface while the PNS remains in a metastable stable. The enclosed mass is
well below the maximum mass of the hybrid EOS. During the PNS collapse, the
central density rises slightly above $\rho_{\rm central}=5\times\rho_{\rm
sat}$, and settles back to a value of about $\rho_{\rm
central}=4.2\times\rho_{\rm sat}$ once the second shock wave accelerates. The
latter quickly reaches the location of the bounce shock, at a radius of about
50 km (see graph (d) in Fig. 3), taking over the bounce shock and expanding
continuously to increasingly larger radii. It reaches a radius of 1000 km
about 10 ms after its formation (thin red dash-dotted line in graph (d) of
Fig. 3). It defines the onset of the SN explosion triggered by the PNS
collapse due to the hadron–quark phase transition, indicated by the sudden
rise of the diagnostic energy, $E_{\rm diagnostic}$, shown in graph (e) of
Fig. 3,

$E_{\rm diagnostic}(t)=-\int_{\rm M}^{a_{\rm cut}}da\,E_{\rm
specific}(t,a)~{},$ (2)

which is the baryon mass integration, $da$, of the total specific energy,
$E_{\rm specific}$, integrated from the stellar surface, M, towards the
center. The specific energy is given as follows,

$E_{\rm specific}(t,a)=\Gamma
e+\frac{2}{\Gamma+1}\left(\frac{u^{2}}{2}-\frac{m}{r}\right)~{},$ (3)

with metric function, $\Gamma=\sqrt{1+m/r-u}$, specific internal energy, $e$,
velocity, $u$, and the gravitational mass, $m$. The integration is performed
until the mass cut, $a_{\rm cut}$. The latter defines the mass coordinate
inside of which material is not ejected, which is determined at the end of the
simulation time of about 12 s post bounce. Note that it includes the neutrino-
driven wind phase. Note further that the diagnostic energy contains the
gravitational binding energy of the stellar envelop, since the entire
progenitor mass is included in the SN simulations.

The diagnostic energy evolves during the SN simulations (see graph (e) in Fig.
3). It is often used as indicator for a successful SN explosion when $E_{\rm
diagnostic}>0$, referring to the SN explosion energy which is the kinetic
energy of the ejecta which should match the asymptotic value obtained for
$E_{\rm diagnostic}$ (see the inlay of graph (e) in Fig. 3).

The entire SN evolution for the solar-metallicity model is illustrated in the
temperature-density phase diagram, Fig. 7, with the color coding of the
electron fraction. The evolution of the central values is given as thick black
line. It is interesting to note that the central evolution across the phase
co-existence, marked by the blue thick solid (onset density) and thin dash-
dotted lines (reaching the pure quark matter phase), proceeds along curves of
constant entropy per particle of about $s=1.5-3~{}k_{\rm B}$, which indicates
the adiabatic collapse of the PNS. As the entropy is conserved, i.e. neutrinos
cannot leave the hot and dense PNS interior, the temperature decreases during
the phase transition slightly. This phenomenon has already been realized in
Ref. Fischer:2018 . Consequently, when considering the stability of hybrid
matter PNS, one should take into account the maximum mass of compact stars at
finite entropy per particle (see Table 1). It is interesting to note that the
SN evolution in the phase diagram is similar as those of binary neutron star
mergers, which consider the hadron quark phase transition Bauswein:2019PhRvL ;
Most:2019PhRvL ; Blacker:2020 , however, reaching somewhat higher
temperatures.

The rapid SN shock expansion of the second shock across the neutrinospheres
releases an outburst of neutrinos of all flavors, however, dominated by
$\bar{\nu}_{e}$ and heavy lepton flavor neutrinos (see Fig. 6). This burst
lasts for only a few tens of a second. The associated luminosities and average
neutrino energies rise to values of several times $10^{53}$ erg s-1 and 40
MeV, respectively. It has been pointed out that this non-standard neutrino
signal—not present in neutrino-driven SN simulation models—is a smoking-gun
signature for a sufficiently strong first-order phase transition. It is
observable at the present generation of neutrino detectors, e.g., Super-
Kamiokande Fischer:2018 .

##  5 Summary and conclusions

The present article reports about results of two core-collapse SN simulations
of very massive progenitor stars with the same ZAMS mass of 75 M⊙ but
different metallicities. This important aspect was not explored in previous
studies. The metallicity has important consequences for the mass loss during
the stellar evolution; high(low) metallicity results in high(low) mass loss
rates and consequently in high(low) compactness of the nascent iron core at
the onset of stellar core collapse. This, in turn, results in substantial
differences during the SN post-bounce evolution, featuring low(high) mass
accretion rates. This study considers an EOS with a first-order phase
transition from normal nuclear matter to the quark–gluon plasma at high baryon
density. In the subsequent SN simulations, there are two competing time-
scales: (i) the growth of the central density and temperature to reach the
critical conditions for the onset of the phase transition, and (ii) the growth
of the mass enclosed inside the PNS. Both are determined by the mass accretion
rate of the still collapsing stellar core onto the bounce shock and by the
compression behavior of the PNS, which, in turn is determined by the high-
density behavior of the nuclear EOS. Here, the previous finding are confirmed
that a high mass accretion rate may result in a massive PNS at the onset of
the hadron–quark phase transition, exceeding the maximum mass of the EOS. Only
when the timescale for the PNS compression, i.e. the timescale associated with
the rising central density and temperature, exceeding the growth timescale of
the PNS mass, the nascent PNS will be in a meta-stable state after the phase
transition. The formation of a second shock wave determines the onset of the
SN explosion. In the present work, this is the case for the progenitor of
solar metallicity only. Contrary, for the low-metallicity model, the enclosed
mass exceeds the maximum hybrid EOS mass and hence a black hole forms. The SN
explosion of the solar-metallicity model confirms the unique observable
signature in the neutrino signal.

These findings have important consequences for the dependence on metallicity
of this class of massive star explosion, driven by a sufficiently strong
hadron–quark phase transition. While the first generation of stars had zero
metal content, the metal enrichment is considered to rise continuously during
the evolving universe. Such considerations lead to the age–metallicity
relation Hirai:2015 . Here, a metallicity of $10^{-4}$ of the solar value is
considered as representative example of extremely metal poor stars
Umeda:2003ut ; Tominaga:2007pb . Observations of metal poor stars with surface
abundances enriched with heavy $r$-process elements, such as strontium and
europium, have been observed frequently Sneden:2008 ; Roederer:2012ApJS ;
Roederer:2012ApJ ; Frebel:2018 ; Cowan:2021RvMP . It is yet incompletely
understood which astrophysical sites can account for these observations,
especially at metalicities of $Z=10^{-4}\times Z_{\odot}$ and below.
Specifically, whether binary neutron star mergers can account therefore is
presently unclear Argast:2004A&A ; Wehmeyer:2015 ; Wehmeyer:2019 . On the
other hand, even though SN explosions could have enriched the galaxy at low
metallicity, canonical, neutrino-driven core-collapse SN explosions cannot
account for the production of heavy neutron-capture elements heavier than
molybdenum, with atomic number of $Z=42$ MartinezPinedo:2014 . The reason for
this is the absence of sufficient neutron rich ejecta. It points to rare and,
perhaps speculative sites for the $r$ process associated with massive star
explosions, such as magneto-rotational driven jet-like explosions
Winteler:2012 ; Moesta:2014 ; Moesta:2015 ; Kotake:2020 . In addition, SN
driven by the hadron–quark phase transition have been suggested as possible
but rare $r$-process site Fischer:2020b , with an abundance pattern featuring
an overproduction of nuclei in the vicinity of the second $r$-process peak in
the region with nuclear mass numbers $140<A<180$. It was suggested that this
site may operate at low metallicity and hence account for the observed metal-
poor star data. However, the findings of the present study points to the
absence of this class of massive star explosions driven by the hadron-qaurk
phase transition at a metallicity of $Z\leq 10^{-4}\times Z_{\odot}$.

Zero- and low-metalicty stellar evolution calculations predict qualitatively
the same evolutionary pattern for massive stars with ZAMS masses $>30$ M⊙
Woosley:2002zz ; Umeda:2003ut ; Umeda:2007wk ; Tominaga:2007pb , namely an
extended high-density silicon-sulfur layer above the iron core due to
negligible mass loss. Consequently, all these massive stars would belong to
the failed SN branch with the failed SN shock revival and hence the formation
of black holes. It excludes SN explosions of massive stars driven by the
hadron–quark phase transition as $r$-process site operating at low
metallicity. At which value of the increasing metallicity during the galactic
evolution they might have started to contribute remains to be explored.

Acknowledgements The author acknowledges support from the Polish National
Science Center (NCN) under Grant No. 2019/ 33/B/ST9/03059. The SN simulations
were performed at the Wroclaw Center for Scientific Computing and Networking
(WCSS) in Wroclaw (Poland).

Data Availability Statement This manuscript has associated data in a data
repository. The data associated with the equation of state and the supernova
simulations are available from the corresponding author on reasonable request.

Open Access This article is licensed under a Creative Commons Attribution 4.0
International License, which permits use, sharing, adaptation, distribution
and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the
Creative Commons licence, and indicate if changes were made. The images or
other third party material in this article are included in the article’s
Creative Commons licence, unless indicated otherwise in a credit line to the
material. If material is not included in the article’s Creative Commons
licence and your intended use is not permitted by statutory regulation or
exceeds the permitted use, you will need to obtain permission directly from
the copyright holder. To view a copy of this licence, visit
http://creativecommons.org/licenses/by/4.0.

## References

  * (1) T. Fischer, N.-U. Bastian, D. Blaschke, M. Cierniak, M. Hempel, T. Klähn, G. Martínez-Pinedo, W. G. Newton, G. Röpke, S. Typel, Publ. Astron. Soc. Austr. 34, e067 (2017).
  * (2) J. M. LeBlanc, J. R. Wilson, Astrophys. J. 161, 541 (1970)
  * (3) H. A. Bethe, R. Wilson, James, Astrophys. J. 295, 14 (1985)
  * (4) H.-T. Janka, K. Langanke, A. Marek, G. Martínez-Pinedo, Phys. Rep. 422, 38 (2007).
  * (5) A. Mirizzi, I. Tamborra, H.-T. Janka, N. Saviano, K. Scholberg, R. Bollig, L. Hüdepohl, S. Chakraborty, Nuovo Cimento Rivista Serie 39, 1 (2016)
  * (6) K. Nomoto, Astrophys. J. 322, 206 (1987)
  * (7) F. Kitaura, H.-T. Janka, W. Hillebrandt, Astron. Astrophys. 450, 345 (2006)
  * (8) T. Melson, H.-T. Janka, A. Marek, Astrophys. J. 801, L24 (2015)
  * (9) M. Takahara, K. Sato, Prog. Theor. Phys. 80, 861 (1988)
  * (10) N. Gentile, M. Aufderheide, G. Mathews, F. Swesty, G. Fuller, Astrophys. J. 414, 701 (1993)
  * (11) I. Sagert, T. Fischer, M. Hempel, G. Pagliara, J. Schaffner-Bielich, et al., Phys. Rev. Lett. 102, 081101 (2009)
  * (12) A. Bazavov, et al., HotQCD Collaboration, Phys. Rev. D 90, 094503 (2014)
  * (13) S. Borsányi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, K. K. Szabó, Phys. Lett. B 730, 99 (2014)
  * (14) A. Bazavov, et al., HotQCD Collaboration, Phys. Lett. B 795, 15 (2019)
  * (15) E. Farhi, R. Jaffe, Phys. Rev. D 30,2379 (1984)
  * (16) J. Antoniadis, et al., Science 340, 448 (2013)
  * (17) E. Fonseca, et al., Astrophys. J. 832, 167 (2016)
  * (18) H. T. Cromartie, et al., Nat. Astron. 4, 72 (2020)
  * (19) E. Fonseca, et al., Astrophys. J. 915, L12 (2021)
  * (20) S. Benić, D. Blaschke, D. E. Alvarez-Castillo, T. Fischer, S. Typel, Astron. Astrophys. 577, A40 (2015)
  * (21) T. Klähn, T. Fischer, Astrophys. J. 810, 134 (2015)
  * (22) M. A. R. Kaltenborn, N.-U. F. Bastian, D. B. Blaschke, Phys. Rev. D 96, 056024 (2017)
  * (23) T. Klähn, T. Fischer, M. Hempel, Astrophys. J. 836, 89 (2017)
  * (24) Y. Nambu, G. Jona-Lasinio, Phys. Rev. 122, 345 (1961)
  * (25) M. Buballa, Phys. Rept. 407, 205 (2005)
  * (26) T. Klähn, R. Łastowiecki, D. Blaschke, Phys. Rev. D 88, 085001 (2013)
  * (27) C. D. Roberts, S. M. Schmidt, Prog. Part. Nucl. Phys. 45, S1 (2000)
  * (28) H. Chen, W. Yuan, L. Chang, Y.-X. Liu, T. Kl’́ahn, et al., Phys. Rev. D 78, 116015 (2008)
  * (29) H. Chen, M. Baldo, G. Burgio, H.-J. Schulze, Phys. Rev. D 84, 105023 (2011)
  * (30) A. Bashir, L. Chang, I. C. Cloet, B. El-Bennich, Y.-X. Liu, et al., Commun. Theor. Phys. 58, 79 (2012)
  * (31) I. C. Cloet, C. D. Roberts, “Explanation and Prediction of Observables using Continuum Strong QCD,” Prog. Part. Nucl. Phys. 77, 1 (2014)
  * (32) H. Chen, J. B. Wei, M. Baldo, G. Burgio, H. J. Schulze, Phys. Rev. D 91, 105002 (2015)
  * (33) A. Kurkela, P. Romatschke, A. Vuorinen, Phys. Rev. D 81, 105021 (2010)
  * (34) A. Kurkela, E. S. Fraga, J. Schaffner-Bielich, A. Vuorinen, Astrophys. J. 789, 127 (2014)
  * (35) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, A. Vuorinen, Nat. Phys. 16, 907 (2020)
  * (36) T. Fischer, N.-U. F. Bastian, M.-R. Wu, P. Baklanov, E. Sorokina, S. Blinnikov, S. Typel, T. Klähn, D. B. Blaschke, Nat. Astron. 2, 980 (2018)
  * (37) N.-U. F. Bastian, Phys. Rev. D 103, 023001 (2021)
  * (38) A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, K. Chatziioannou, J. A. Clark, T. Fischer, and M. Oertel, “Identifying a First-Order Phase Transition in Neutron-Star Mergers through Gravitational Waves,” Phys. Rev. Lett., vol. 122, p. 061102, Feb. 2019.
  * (39) A. Bauswein, S. Blacker, V. Vijayan, N. Stergioulas, K. Chatziioannou, J. A. Clark, N.-U. F. Bastian, D. B. Blaschke, M. Cierniak, T. Fischer, Phys. Rev. Lett. 125, 141103 (2020)
  * (40) S. Blacker, N.-U. F. Bastian, A. Bauswein, D. B. Blaschke, T. Fischer, M. Oertel, T. Soultanis, S. Typel, Phys. Rev. D 102, 123023 (2020)
  * (41) S. Zha, E. P. O’Connor, M.-c. Chu, L.-M. Lin, S. M. Couch, Phys. Rev. Lett. 125, 051102 (2020)
  * (42) S. Typel, Phys. Rev. C 71, 064301 (2005)
  * (43) S. Typel, G. Röpke, T. Klähn, D. Blaschke, H. Wolter, Phys. Rev. C 81, 015803 (2010)
  * (44) P. Danielewicz, R. Lacey, W. G. Lynch, Science 298, 1592 (2002)
  * (45) M. Hempel, J. Schaffner-Bielich, Nucl. Phys. A 837, 210 (2010)
  * (46) T. Fischer, S. Typel, G. Röpke, N.-U. F. Bastian, G. Martínez-Pinedo, Phys. Rev. C 102, 055807 (2020)
  * (47) S. Typel, Phys. Rev. C 89, 064321 (2014)
  * (48) F. X. Timmes, D. Arnett, Astrophys. J. Suppl. 125, 277 (1999)
  * (49) K. Yagi, N. Yunes, Science 341, 365 (2013)
  * (50) K. Yagi, N. Yunes, Phys. Rev. D 88, 023009 (2013)
  * (51) A. Maselli, V. Cardoso, V. Ferrari, L. Gualtieri, P. Pani, Phys. Rev. D 88, 023007 (2013)
  * (52) A. R. Raduta, M. Oertel, A. Sedrakian, Mon. Not. Roy. Astron. Soc. 499, 914 (2020)
  * (53) S. Khadkikar, A. R. Raduta, M. Oertel, A. Sedrakian, Phys. Rev. C 103, 055811 (2021)
  * (54) M. Hempel, O. Heinimann, A. Yudin, I. Iosilevskiy, M. Liebendörfer, F.-K. Thielemann, Phys. Rev. D 94, 103001 (2016)
  * (55) M. Liebendörfer, O. Messer, A. Mezzacappa, S. Bruenn, C. Cardall, et al. Astrophys. J. Suppl. 150, 263 (2004)
  * (56) T. Fischer, G. Guo, G. Martínez-Pinedo, M. Liebendörfer, A. Mezzacappa, Phys. Rev. D 102, 123001 (2020)
  * (57) T. Fischer, G. Guo, A. A. Dzhioev, G. Martínez-Pinedo, M.-R. Wu, A. Lohs, Y.-Z. Qian, Phys. Rev. C 101, 025804 (2020)
  * (58) G. Guo, G. Martínez-Pinedo, A. Lohs, T. Fischer, Phys. Rev. D 102, 023037 (2020)
  * (59) M. Liebendörfer, S. Rosswog, F.-K. Thielemann, Astrophys. J. Suppl. 141, 229 (2002)
  * (60) T. Fischer, S. Whitehouse, A. Mezzacappa, F.-K. Thielemann, M. Liebendörfer, Astron. Astrophys. 517, A80 (2010)
  * (61) S. W. Bruenn, Astrophys. J. Suppl. 58, 771 (1985)
  * (62) A. Mezzacappa, S. Bruenn, Astrophys. J. 405, 637 (1993)
  * (63) S. Woosley, A. Heger, T. Weaver, Rev/ Mod. Phys. 74, 1015 (2002)
  * (64) H. Umeda, K. Nomoto, Astrophys. J. 673, 1014 (2008)
  * (65) A. M. Serenelli, S. Basu, J. W. Ferguson, M. Asplund, Astrophys. J. 705, L123 (2009)
  * (66) M. Asplund, A. M. Amarsi, N. Grevesse, arXiv e-prints, astro-ph/2105.01661 (2021)
  * (67) E. J. Lentz, S. W. Bruenn, W. R. Hix, A. Mezzacappa, O. E. B. Messer, E. Endeve, J. M. Blondin, J. A. Harris, P. Marronetti, K. N. Yakunin, Astrophys. J. 807, L31 (2015)
  * (68) B. Müller, Mon. Not. Roy. Astron. Soc. 453, 287 (2015)
  * (69) H. Andresen, B. Müller, E. Müller, H. T. Janka, Mon. Not. Roy. Astron. Soc. 468, 2032 (2017)
  * (70) K. Kotake, T. Takiwaki, T. Fischer, K. Nakamura, G. Martínez-Pinedo, Astrophys. J. 853, 170 (2018)
  * (71) A. Burrows, D. Radice, D. Vartanyan, Mon. Not. Roy. Astron. Soc. 485, 3153 (2019)
  * (72) R. Bollig, N. Yadav, D. Kresse, H. T. Janka, B. Mueller, A. Heger, Astrophys. J. 915, 28 (2021)
  * (73) T. Kuroda, K. Kotake, T. Takiwaki, F.-K. Thielemann, Mon. Not. Roy. Astron. Soc. 477, L80 (2018)
  * (74) E. R. Most, L. J. Papenfort, V. Dexheimer, M. Hanauske, S. Schramm, H. Stöcker, L. Rezzolla, Phys. Rev. Lett. 122, 061101 (2019)
  * (75) Y. Hirai, Y. Ishimaru, T. R. Saitoh, M. S. Fujii, J. Hidaka, T. Kajino, Astrophys. J. 814, 41 (2015)
  * (76) H. Umeda, K. Nomoto, Astrophys. J. 619, 427 (2005)
  * (77) N. Tominaga, H. Umeda, K. Nomoto, Astrophys. J. 660, 516 (2007)
  * (78) C. Sneden, J. J. Cowan, R. Gallino, Ann. Rev. Astron. Astrophys. 46, 241 (2008)
  * (79) I. U. Roederer, J. E. Lawler, J. S. Sobeck, T. C. Beers, J. J. Cowan, A. Frebel, I. I. Ivans, H. Schatz, C. Sneden, I. B. Thompson, Astrophys. J. Suppl. 203, 27 (2012)
  * (80) I. U. Roederer, J. E. Lawler, Astrophys. J. 750, 76 (2012)
  * (81) A. Frebel, Ann. Rev. Nucl. Part. Sci. 68, 237 (2018)
  * (82) J. J. Cowan, C. Sneden, J. E. Lawler, A. Aprahamian, M. Wiescher, K. Langanke, G. Martínez-Pinedo, F.-K. Thielemann, Rev/ Mod. Phys. 93, 015002 (2021)
  * (83) D. Argast, M. Samland, F. K. Thielemann, Y. Z. Qian, Astron. Astrophys. 416, 997 (2004)
  * (84) B. Wehmeyer, M. Pignatari, F. K. Thielemann, Mon. Not. Roy. Astron. Soc. 452, 1970 (2015)
  * (85) B. Wehmeyer, C. Fröhlich, B. Côté, M. Pignatari, F. K. Thielemann, Mon. Not. Roy. Astron. Soc. 487, 1745 (2019)
  * (86) G. Martínez-Pinedo, T. Fischer, L. Huther, J. Phys. G Nucl. Part. Phys. 41, 044008 (2014)
  * (87) C. Winteler, R. Käppeli, A. Perego, A. Arcones, N. Vasset, N. Nishimura, M. Liebendörfer, F. K. Thielemann, Astrophys. J. 750, L22 (2012)
  * (88) P. Mösta, S. Richers, C. D. Ott, R. Haas, A. L. Piro, K. Boydstun, E. Abdikamalov, C. Reisswig, E. Schnetter, Astrophys. J. 785, L29 (2014)
  * (89) P. Mösta, C. D. Ott, D. Radice, L. F. Roberts, E. Schnetter, R. Haas, Nature 528, 376 (2015)
  * (90) T. Kuroda, A. Arcones, T. Takiwaki, K. Kotake, Astrophys. J. 896, 102 (2020)
  * (91) T. Fischer, M.-R. Wu, B. Wehmeyer, N.-U. F. Bastian, G. Martínez-Pinedo, F.-K. Thielemann, Astrophys. J. 894, 9 (2020)