longtable

# Jointly Estimating Subnational Mortality for Multiple Populations

Ameer Dharamshi1,2
Monica Alexander1
Celeste Winant3
Magali Barbieri3,4

1 University of Toronto
2 University of Washington
3 University of California, Berkeley
4 Institut National d’Études Démographiques

(October 04, 2023)

###### Abstract

Understanding patterns in mortality across subpopulations is essential for
local health policy decision making. One of the key challenges of subnational
mortality rate estimation is the presence of small populations and zero or
near zero death counts. When studying differences between subpopulations, this
challenge is compounded as the small populations are further divided along
socioeconomic or demographic lines. In this paper, we build on principal
component-based Bayesian hierarchical approaches for subnational mortality
rate estimation to model correlations across subpopulations. The principal
components identify structural differences between subpopulations, and
coefficient and error models track the correlations between subpopulations
over time. We illustrate the use of the model in a simulation study as well as
on county-level sex-specific US mortality data. We find that results from the
model are reasonable and that it successfully extracts meaningful patterns in
US sex-specific mortality. Additionally, we show that ancillary correlation
parameters are a useful tool for studying the convergence and divergence of
mortality patterns over time.

##  1 Introduction

Reliable mortality estimates are crucial to support health policy design,
implementation, and monitoring. Historically, developments in both mortality
modelling and understanding mortality differences across populations have
focused on national-level estimates for cross-country comparisons. However, as
demonstrated by a growing body of evidence (Ezzati 2008; Kindig and Cheng
2013; Murray 2006; Sehgal et al. 2022), patterns in more granular subnational
level mortality rates are needed to facilitate local level decisions. This
insight has led to a recent line of work focused on studying disparities in
mortality outcomes within countries at the subnational level. In the United
States, for example, researchers are analyzing differences in life expectancy
across state jurisdictions (Woolf and Schoomaker 2019; Montez et al. 2020;
Harper, Riddell, and King 2021). At the international level, the United
Nations has recently moved to producing estimates of key global health
indicators, such as the under-five mortality rate, at subnational levels
alongside the longstanding national level estimates (Godwin and Wakefield
2021). In addition, in recent years it has become increasingly clear that
disparities in mortality outcomes along demographic dimensions such as sex or
gender, race or ethnicity, or socioeconomic status, are growing (Hendi 2015;
Crimmins and Zhang 2019; Gutin and Hummer 2021; Masters and Aron 2021). As an
example, in the case of child mortality risks, Bhutta (2016) argued that
subnational determinants explain a greater portion of mortality than country-
level characteristics. In response, estimates of subpopulation mortality risks
are needed to identify and understand the mortality patterns of vulnerable
groups, and track the effects of policy responses.

In this paper, we focus on sex differences in mortality. Differences in
mortality patterns between males and females have long been documented. Women
have had longer life expectancy than men in every population and time period
where and when female discrimination has not been prevalent. The mortality gap
is particularly large at young adult ages, when mortality is mostly determined
by external causes (particularly accidents and homicides), but it remains high
throughout all working ages. It closes progressively after age 65 to 70 years
but never disappears completely. Sex differences in mortality reached a peak
in most high-income countries around 1975-1985 but have declined since, in
large part because of the convergence in male and female behavior (e.g.,
smoking) (Gjonça, Tomassini, and Vaupel 1999). In these countries, the
difference in life expectancy between men and women ranges from a low of
around 3 years in Northern Europe up to nearly 10 years in Eastern Europe. In
the United States, the sex difference in life expectancy increased throughout
the 20th century to reach a maximum at 7.74 years in 1975 before declining to
a low of 4.74 years in 2012. It has increased again since and currently (2021)
reaches 5.75 years (Human Mortality Database 2023). The difference varies
across states, ranging from 3.7 years (in Utah) to 7.2 years (in the District
of Columbia) (United States Mortality Database 2023). Given that the sex gap
in life expectancy is highly location-specific, it is of interest to further
study the behaviour of sex-specific mortality risks at finer geographic
resolutions such as counties. Hence, there is a need for methods that reliably
estimate differential mortality risks between males and females (and between
other demographic subgroups) in small populations.

While the need for subnational estimates of mortality disaggregated by
subpopulation is clear, there are several challenges to obtaining reliable
estimates. Although raw death rates can usually be calculated, the natural
variations lead to uncertain estimates. It is well understood that subnational
mortality estimates are more complex to construct than national level
estimates given the large random fluctuations associated with small numbers
(Stevenson and Olson 1993). This problem is further compounded when examining
mortality across sub-groups within sub-national areas, where inferring
underlying mortality risks from raw death rates is challenged by the erratic
patterns over age.

Producing reliable estimates of mortality rates and measuring the associated
uncertainty requires statistical models that take the natural stochasticity in
the data into account. In recent years, an increasing number of studies have
demonstrated the usefulness of the Bayesian approach in the construction of
such models (see for instance Alexander, Zagheni, and Barbieri (2017); Khana
and Warner (2018); Song and Luan (2022)). Building on previous work, in this
paper, we propose a general Bayesian model framework to estimate age-specific
mortality rates at the subnational level jointly for multiple populations. The
model incorporates characteristic shapes of mortality age schedules within a
Bayesian hierarchical framework, which allows information on mortality
patterns to be shared across populations. We extend previous approaches by
accounting for correlation in mortality experiences across subpopulations,
rather than assuming subpopulations are independent. The model produces
estimates and uncertainty for mortality rates for each sex. In addition, it
estimates higher-level parameters. These parameters summarise trends in
different dimensions of mortality over time and show how similar or dissimilar
the trends are across groups. We illustrate the model with an application to
estimating sex-specific mortality by county in the United States. While we
focus on sex-specific mortality, the modelling framework is generalizable to
population groups defined by other characteristics (such as race/ethnicity).

The remainder of the paper is structured as follows. We first give a brief
overview of recent developments in subnational mortality estimation. We begin
our methods discussion by demonstrating that principal components derived from
a set of reference mortality curves capture structural differences in
subpopulation specific mortality patterns, then follow with a formal statement
of the model. Results from the model along with validation exercises using
simulated and real sex-specific US county-level mortality data are then
provided. Finally, we conclude with a discussion of our findings and identify
directions for future research.

##  2 Background

A large body of research exists on small area estimation issues, and recently
demographers have increasingly taken advantage of computational advances which
make fitting complex statistical models to small-scale mortality data
feasible. In particular, there has been a notable increase in the use of
Bayesian methods in demographic estimation, particularly subnational
estimation. Bayesian methods are particularly suited to demographic contexts
as they provide a useful framework to incorporate different data sources in
the same model, account for various types of uncertainty, and allow for
information exchange across time and space (Bijak and Bryant 2016).

One area of previous work has focused on estimating aggregate indicators at
the subnational level, such as child mortality and life expectancy (Mercer et
al. 2015; Ševčíková and Raftery 2021). Models on aggregate indicators
generally involve temporal and spatial smoothing, allowing for information in
mortality trends to be shared across these dimensions. In some cases, models
rely on covariates (such as education or income) to stabilize mortality rate
estimates from noisy data (Wang et al. 2013; Arias et al. 2018; Murray 2006).

In addition to aggregate mortality indicators, research has focused on
producing estimates of age-specific mortality rates, which can form the basis
of estimating subnational life tables. Recent advances in this area build on
classical demographic approaches of model life tables and relational models,
which identify key patterns in mortality over age across a wide range of
populations, and allow patterns to be shifted based on a reduced set of
parameters (Coale, Demeny, and Vaughan 1983). For example, TOPALS models
consist of a standard age schedule and population-specific deviations away
from that standard, which are smoothed using linear splines (de Beer 2012).
TOPALS-type models have been used to produce subnational estimates of
migration and mortality in varying data quality contexts (Gonzaga and
Schmertmann 2016; Schmertmann and Gonzaga 2018; Dyrting 2020). TOPALS models
do, however, require selecting a single standard mortality curve. This is
challenging when studying multiple populations, as each may have a distinct
mortality profile.

A related modelling approach derives a set of ‘principal components’ from
reference mortality curves which are then used as the basis of a regression
framework. This allows a large set of plausible mortality curves to be
estimated using a reduced set of parameters. This approach is in the spirit of
the Lee-Carter model for mortality forecasting and related work (Lee and
Carter 1992). For example, Clark (2019) uses this approach to formulate a new
set of model age schedules in data-sparse contexts. In particular, this paper
extends an earlier paper by Alexander, Zagheni, and Barbieri (2017) which
introduces a Bayesian hierarchical framework that builds on principal
components derived from national mortality schedules for use at the
subnational level. In more recent work, Alexander and Alkema (2022) have used
a principal components approach to model mortality in the broader context of
subnational population estimation using a cohort component projection
framework.

In general, most existing mortality models inherently assume subpopulations
are independent. They model each population separately, then combine or
interpret after modelling. However, ideally we would model all subpopulations
within one framework as mortality is generally correlated over groups.
Existing joint models focus on modelling mortality by sex as it is well
understood that changes in male and female mortality are correlated (Noymer
and Van 2014). At the national level, sex-specific life expectancy estimates
have been produced by the UN using gap-based approaches such as in Raftery,
Lalic, and Gerland (2014). Under this approach, female life expectancy is
estimated with standard one-sex methods. Estimates of male life expectancy are
then constructed by modelling the gap between female and male life expectancy.
In the present case, we are interested in estimates of all age-specific
mortality rates, not just life expectancy, and a purely gap-based approach to
this problem at smaller scales may not be appropriate due to the small counts
in both areas and age groups. Additionally, gap-based approaches require the
selection of an anchor population such as females in the case of sex-specific
modelling. This works in situations where the composition within a population
is roughly balanced, but could be problematic when considering variables such
as race or ethnicity that vary dramatically across jurisdictions, particularly
if the anchor population is near absent in certain regions.

At the subnational level, Rau and Schmertmann (2020) jointly model age- and
sex-specific mortality in regions across Germany using a Bayesian TOPALS
approach. With regards to sex differences, they obtain a typical pattern of
age-specific sex differences across aggregated regions of Germany, and then
use this information as a basis for a prior for differences at smaller scales,
penalizing large deviations away from the observed differences by sex. In this
paper, we take a different approach, and flip the viewpoint from thinking
about sex differences, to thinking about covariation in mortality by sex (or
any other population subgroups), and explicitly allow for group mortality
rates to move together.

##  3 Methods

To estimate subnational mortality rates while also extracting patterns across
key subpopulations, we propose a Bayesian hierarchical model that builds on a
principal component-based approach, and incorporates structures that capture
subpopulation interactions. Before defining the model equations, we first
discuss why principal components are an attractive modelling strategy in this
context.

###  3.1 Principal components models

The use of principal components is motivated by the fact that age-specific
mortality rates tend to display strong regularities across different
populations. This means that systematic variation in mortality rates is well
captured by a reduced set of parameters which can be modelled such that
information is shared across space and time.

To generate principal components, we compute the singular value decomposition
of a collection of regional log-mortality curves. The collection of log-
mortality curves, $\mathbf{X}$, is a $N\times A$ matrix where $N$ is the
number of region-subpopulation-years under consideration and $A$ is the number
of age groups. For US sex-specific mortality, $N=6066$ (60 years of state-
level sex-specific data) and $A=19$ (ages <1, 1-4, 5-9, …, 75-79, 80-84, 85+).
Note that different age groups, such as one-year age groups, could also be
considered. The singular value decomposition of $\mathbf{X}$ is then:

$\mathbf{X}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{\prime}$ (1)

where $\mathbf{U}$ is the $6066\times 19$ matrix of left singular values,
$\boldsymbol{\Sigma}$ is a $19\times 19$ diagonal matrix of scaling factors,
and $\mathbf{V}$ is the $19\times 19$ matrix of right singular values, which
we term ‘principal components’. These extract key patterns in mortality over
age. As the $\mathbf{X}$ matrix contains data from both male and female
populations, the patterns captured by the principal components are shared
across sexes.

The first four principal component curves are plotted in Figure 1. The first
principal component explains the most variation across mortality curves with
each successive component explaining less variation. The first four principal
components explain over 99% of the variation in the state-level data. Note
that the first principal component is inverted for clarity as its coefficients
are exclusively negative.

Figure 1: First four state-level log-mortality principal components.

In practice, choosing the number of components is a balance between
incorporating components that only pick up on systematic patterns, while still
allowing for enough flexibility in the model. The set of principal components
that offer useful information on sex-specific mortality can be identified
using the $\mathbf{U}$ matrix. The left singular values in each row of the
$\mathbf{U}$ matrix represent the contribution of each principal component to
the corresponding mortality curve in the $\mathbf{X}$ matrix. If the $i$th
principal component contains material subpopulation-specific differences, we
would expect the distributions of the values in the $i$th column of
$\mathbf{U}$ to differ between subpopulations. We thus separate the rows of
$\mathbf{U}$ by sex and examine the resulting distributions. Figure 2 displays
the densities of the coefficients by sex for the first eight principal
components. We choose eight as the singular values suggest there is limited
additional information beyond eight.

Figure 2: Distribution of observed state-level left singular values by sex and
principal component.

The first four principal components demonstrate clear differences in the sex-
specific distributions. In particular, the location of the distributions are
noticeably different suggesting structural differences between sexes. The
distributions of the remaining principal components, including those not
presented in Figure 2, are largely similar. Beyond a visual inspection, simple
t-tests suggest that aside from the fifth, all of the first eight principal
components have location differences. Balancing these findings with the
practical consideration that each additional principal component will lead to
a significant computational burden during model estimation, we suggest that at
least three but ideally four principal components should be used in modelling
sex-specific mortality in the US.

Considering Figure 1 and Figure 2 together offers insight into patterns of
interest. The first principal component has the characteristic ‘J’ shape of
log mortality. The lower female coefficient values in the PC1 panel of Figure
2 indicate that female mortality is generally lower than male mortality.
Similarly, the PC3 panel implies that women experience lower mortality through
the middle-age adult years. In addition, the distribution of coefficients for
the second and fourth principal components indicate higher young adult
mortality in the male population than the female population which is
consistent with the “accident hump” described in the broader mortality
literature, as males are more likely to suffer mortality due to risky
behaviour (Heligman and Pollard 1980).

There is a substantial body of literature documenting the large and systematic
excess mortality of men compared with women (see reviews of this literature in
Gjonça, Tomassini, and Vaupel (1999); Barford et al. (2006); Beltrán-Sánchez,
Finch, and Crimmins (2015)). As previously mentioned, the male disadvantage in
the risk of dying is high at all adult ages up to about 65 to 70 years,
depending on the population, and it is particularly acute in young adults. The
finding that the principal components are recovering these well understood
sex-specific differences motivates their use as the basis for a subpopulation
mortality model. Using principal components as a basis for a statistical model
can be robust in small population settings as the components allow plausible
mortality rates to be estimated even in the presence of high variability.

###  3.2 Model summary

We begin by defining $y_{a,s,c,t}$ as the observed number of deaths in age
group $a$, subpopulation $s$, county $c$, and year $t$. Then, we assume that

$y_{a,s,c,t}|\lambda_{a,s,c,t}\sim\text{Poisson}\left(P_{a,s,c,t}\cdot\lambda_{a,s,c,t}\right),$
(2)

where $P_{a,s,c,t}$ is the population corresponding to age group $a$,
subpopulation $s$, county $c$, and year $t$, and $\lambda_{a,s,c,t}$ is the
mortality rate to be estimated for age group $a$, subpopulation $s$, county
$c$, and year $t$.

The mortality rates $\lambda_{a,s,c,t}$ are estimated on the log scale as
follows:

$\log\left(\lambda_{a,s,c,t}\right)=\sum_{i=1}^{P}\left(\beta_{i,s,c,t}\cdot
Y_{i,a}\right)+\gamma_{a,s,c,t},$ (3)

where $Y_{i}$ is the $i$th principal component, $P$ is the number of principal
components, $\beta$ are the estimated coefficients, and $\gamma$ is an
overdispersion term. The number of principal components, $P$ can be selected
based on their contributions and differences between groups. In the case of
sex-specific mortality in US counties, we use the $P=4$ principal components
plotted in Figure 1.

Equation (3) can be intuitively understood as constructing a log-mortality
curve for each age-sex-county-year as a linear combination of the four
principal components plus some additional age-specific variation. Different
linear combinations of the principal components (that is, different estimated
values of the $\beta$ coefficients) lead to different plausible log-mortality
curves. The overdispersion term accounts for the possibility that deaths will
be overdispersed relative to the patterns captured by the principal
components.

####  3.2.1 Core model

For many counties, population and death counts are small and highly variable,
which would lead to uncertain estimates of $\beta$ if each county and year
were estimated independently. As such, we propose a hierarchical model for
$\beta$ based on the nesting structure of counties within states.
Specifically, we will assume that counties within each state are more likely
to share similar mortality patterns than counties in different states, and
allow high-data counties to share information with low-data counties in the
same state. A similar relationship holds across subgroups. Within a county,
while subgroups of the population are expected to experience specific drivers
of mortality due to their unique experience, they are also expected to jointly
experience more general drivers, whether those are due to policy, the
environment, or other factors. For instance, regulations regarding maternity
leaves have an impact on women’s health, but natural catastrophes like heat
waves or hurricanes will affect everyone in the same area (albeit to various
extents). This structure suggests that $\beta$ coefficients, and thus by
extension log-mortality curves, for all subgroups within a county should be
modelled jointly, thereby allowing for distinct mortality estimates in each
group while also exploiting the dependence between groups.

To capture geographic and subgroup dependence, we propose to model the vector
of $\beta$ coefficients for all $S$ groups within each county jointly as
multivariate normal with a common state-level mean vector and covariance
matrix:

$\displaystyle\begin{pmatrix}\beta_{i,1,c,t}\\\ \dots\\\
\beta_{i,S,c,t}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\mu_{\beta_{i,1,t}}\\\ \dots\\\
\mu_{\beta_{i,S,t}}\end{pmatrix}+\begin{pmatrix}\omega_{i,1,c,t}\\\ \dots\\\
\omega_{i,S,c,t}\end{pmatrix},\quad i=1,\dots,P$ (4)
$\displaystyle\left.\begin{pmatrix}\omega_{i,1,c,t}\\\ \dots\\\
\omega_{i,S,c,t}\end{pmatrix}\right|\sigma_{\beta_{i,t}},L_{i,t}^{(\beta)}$
$\displaystyle\sim\mathcal{N}\left(\mathbf{0}_{S},\sigma_{\beta_{i,t}}\mathbf{1}_{S}L^{(\beta)}_{i,t}L_{i,t}^{(\beta)\top}\mathbf{1}_{S}\sigma_{\beta_{i,t}}\right)$
(5) $\displaystyle L^{(\beta)}_{i,t}L_{i,t}^{(\beta)\top}$
$\displaystyle\sim\text{LKJ}(1)$ (6) $\displaystyle\sigma_{\beta_{i,t}}$
$\displaystyle\sim\mathcal{N}\left(0,1\right)$ (7)

where $\mu_{\beta_{i,s,t}}$ is the state-level coefficient for the $i$th
principal component, subpopulation $s$, and year $t$, and $\omega_{i,s,c,t}$
is the county deviation for the $i$th principal component, subpopulation $s$,
county $c$, and year $t$.

We discuss the two components of $\beta$ in turn, starting with the state-
level mean vector, $\mu_{\beta}$. Information-rich observations from county-
subgroups with large populations are the primary contributors to estimates of
$\mu_{\beta}$ for the corresponding subpopulation. Estimates of $\beta$ for
small county-subgroups with less informative observed death counts are then
partially informed by the high-data counties through the shared state-level
means. This pooling effect stabilizes estimates of $\beta$ for small
populations by pulling the estimates towards the state mean. It is important
to note that this effect occurs at the county-subgroup level and not uniformly
within a county: a large population county with an uneven subgroup composition
can experience stronger pooling effects in its smaller subgroups. This is not
particularly important for the application to sex-specific mortality as the
populations are roughly balanced, but may be critical in other applications
involving subgroups defined by other demographic variables.

The $\omega$ vector models the specific deviations in the $\beta$ vector for
each county from the state-level mean vector. By jointly modelling the
deviations for all subpopulations as in (4), the model captures patterns in
the dispersion of $\beta$ by principal component in
$L_{i,t}^{(\beta)}L_{i,t}^{(\beta)\top}$. This enables information sharing
across subpopulations. In counties where there is an imbalance in the number
of non-zero death observations across subpopulations, jointly generating
$\beta$ coefficients allows higher data groups to support lower data groups.
An example of such a county is given in Figure 12 in Appendix B.

We assign the uninformative LKJ(1) prior to the subpopulation correlation
Cholesky factors to allow the observed data to fully determine the correlation
structure (Lewandowski, Kurowicka, and Joe 2009).

By modelling correlations in principal component coefficients between sexes,
it is possible for subgroups to experience joint variation in certain parts of
their mortality curves but not in others. For example, the coefficients
corresponding to the first principal component (ie. baseline mortality) could
experience strong correlation, leading to a common trend in baseline mortality
levels, but the coefficients on the second principal component (ie. the
accident hump) may not, leading to differences by sex in young adult
mortality.

An illustration of the correlated principal component coefficients is provided
in Figure 3. Each plot captures the patterns of one of the four principal
component coefficients for California in 2017. The red point represents the
state-level value, the black points represent each county, and the blue
contours are the distribution of $\omega$ centered at the state means. High
correlation in the first principal component is expected since when the
general conditions in a region improve, both sexes are positively impacted,
leading to highly correlated changes in baseline mortality. Similarly, the low
correlation in the second principal component is consistent with males
experiencing a relatively elevated level of mortality in young adult ages
(i.e. the accident hump) whereas females tend not to. High correlations are
recovered for the third and fourth principal components, suggesting that
changes in young adult and middle-age adult mortality experiences are similar
across sexes.

Figure 3: Median posterior county-level principal component coefficients
($\beta$) and the corresponding state-level mean ($\mu_{\beta}$) for
California in 2017. The blue contour plots indicate the recovered correlation
patters between male and female population.

The $\beta$ specification is intended to be flexible. If there is only a small
number of years of data available or there is prior information suggesting
that correlation structures are time-invariant, one could reasonably omit the
time index or share correlations across a short time interval. Similarly, if
certain correlations are known in advance, the correlation matrices themselves
could be constrained accordingly.

####  3.2.2 Temporal smoothing

State-level means are smoothed over time by penalizing the second-order
differences to produce gradual changes. Smoothing occurs at the state-level as
opposed to the county-level to allow the needed flexibility for counties to
experience irregular mortality patterns driven by localized events such as a
natural disaster, or a pandemic.

$\displaystyle\mu_{\beta_{i,s,t}}|\mu_{\beta_{i,s,t-1}},\mu_{\beta_{i,s,t-2}},\sigma_{\mu_{\beta_{i}}}$
$\displaystyle\sim\mathcal{N}\left(2\cdot\mu_{\beta_{i,s,t-1}}-\mu_{\beta_{i,s,t-2}},\sigma_{\mu_{\beta_{i}}}\right)$
(8) $\displaystyle\sigma_{\mu_{\beta_{i}}}$
$\displaystyle\sim\mathcal{LN}\left(-1.5,0.5\right).$ (9)

####  3.2.3 Overdispersion term

Finally, the $\gamma_{a,s,c,t}$ term that allows for overdispersion of the
log-mortality rate is modelled similarly to $\omega_{i,s,c,t}$ in a
$S$-dimensional multivariate normal setup:

$\displaystyle\left.\begin{pmatrix}\gamma_{a,1,c,t}\\\ \dots\\\
\gamma_{a,S,c,t}\end{pmatrix}\right|\sigma_{a},L_{a,t}^{(\gamma)}$
$\displaystyle\sim\mathcal{N}\left(\mathbf{0}_{S},\sigma_{a}\mathbf{1}_{S}L^{(\gamma)}_{a,t}L_{a,t}^{(\gamma)\top}\mathbf{1}_{S}\sigma_{a}\right)$
(10) $\displaystyle\sigma_{a}$
$\displaystyle\sim\mathcal{N}\left(0,0.25\right)$ (11) $\displaystyle
L^{(\gamma)}_{a,t}L_{a,t}^{(\gamma)\top}$ $\displaystyle\sim\text{LKJ}(1).$
(12)

Relationships between subpopulation $\gamma$’s are captured using the age-year
correlation matrix $L^{(\gamma)}_{a,t}L_{a,t}^{(\gamma)\top}$. An intuitive
way of understanding this component of the model is that the principal
components produce the expected mortality derived from aggregate and local
patterns while $\gamma$ captures additional age-specific deviations. These
deviations from the expectation are often correlated, though this correlation
may differ across ages and over time. Given that the log-mortality values for
different age groups occur in substantially different parts of the log curve,
we estimate a separate scaling factor $\sigma$ for each age group.

###  3.3 Computation

The model described here was fit using a Bayesian framework. Posterior samples
are drawn using the No-U-turn sampling (NUTS) Hamiltonian Monte Carlo
algorithm (Hoffman and Gelman 2014; Neal 2011) implemented in the Stan R
package (Stan Development Team 2021). We execute 4 chains with 500 iterations
of burn-in and 2 500 iterations of samples. Convergence was diagnosed using
trace plots, effective sample sizes, and the Gelman and Rubin diagnostic
(Gelman and Rubin 1992).

###  3.4 Simulation study

We conducted a simulation study to test the ability of the model to
simultaneously estimate mortality rates and extract patterns in subgroup
mortality. We generated 10 years of population data for 25 simulated counties
of various sizes composed of 5 population subgroups ranging in size from 10%
to 50% of the total population. We then generated deaths for all subgroups
using log-mortality rates constructed as the linear combination of a pair of
standard curves with coefficients subject to a variety of correlation
patterns. The goal of the simulation study was to recover the true
correlations and age-specific mortality rates when fitting the proposed model
to the simulated data. Specific details on the data generating process are
provided in Appendix A.

The model is run on the simulated data using the standard curves as the
principal components. To validate the results, we compare the estimated
correlation matrices and log-mortality rates to the corresponding true values.
Figure 4b plots estimated correlation matrices extracted from the model.
Comparing corresponding facets against Figure 4a, it appears that the model
has successfully identified and extracted the patterns observed in the data.

(a) True correlation matrices (b) Estimated correlation matrices

Figure 4: Simulation study true and estimated posterior median correlation
matrices.

Table 1 presents the coverage values for correlation and log-mortality rate
parameters at the 80%, 90%, and 95% nominal levels. Coverage is defined for
correlation and log-mortality rates as
$\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{l_{i}\leq\theta_{i}\leq u_{i}}$ where
$n$ is the number of parameters, $\theta_{i}$ refers to an individual
parameter, and $l_{i}$ and $u_{i}$ are respectively the lower and upper bounds
of the posterior credible intervals. For the correlation matrices, we compute
entrywise coverage values for the off-diagonal entries of the lower triangles
of the correlation matrices.

Table 1: Simulation study coverage values for correlation and log-mortality rates. | Coverage (80%) | Coverage (90%) | Coverage (95%)
---|---|---|---
Correlations | 0.78 | 0.90 | 0.94
Log-mortality rates | 0.83 | 0.92 | 0.96

For both sets of parameters, the model’s coverage values are in line with the
nominal values at all levels, suggesting that the model is well calibrated.

##  4 Application to age- and sex-specific mortality rates at the county level
in the United States

We applied the model to estimate sex- and age-specific mortality at the county
level in the US over the period 1982-2019. In the US, relatively high quality
data are collected through vital registration systems. The National Center for
Health Statistics (NCHS) publishes detailed mortality records for each year
included in our series from which we can compile the mortality tabulations.
For purposes of the present analysis, we obtained access to the necessary
protected data with individual death records from 1989 to 2019 from the NCHS
under a Data User Agreement as the publicly available data do not include
geographic information. The United States Census Bureau publishes mid-year
population estimates by county of residence, year, sex, age. These data are
available to the public through the Census Bureau. To comply with our data
agreement, counties are anonymized in all figures.

The model is applied to each state individually using the first four principal
components plotted in Figure 1. Note that estimates based on this model are
published as part of the United States Mortality DataBase (USMDB)
(https://usa.mortality.org/).

As an illustration of the model outputs, in Figure 5 we plot the estimated
sex-specific log-mortality curves in 1982 and 2019 for a subset of US counties
along with associated uncertainty intervals. The observed log-mortality rates
for age groups with non-zero deaths are presented as points, and the female
and male county-level estimates and 95% credible intervals are presented in
red and blue respectively. The selected counties represent settings with
small, medium, and large populations (of approximately 13,000, 40,000, and
10,000,000 people, respectively). As expected, smaller counties such as County
1 have greater uncertainty as compared to larger ones such as County 3 due to
the higher levels of noise caused by low or zero death counts. For County 3,
modelled estimates follow the data exactly.

Figure 5: Examples of US county-level sex-specific log-mortality plots by
county-year-sex for three counties of varying sizes. Black dots represent
observed values and coloured curves and regions indicate posterior medians and
95% credible intervals respectively.

###  4.1 State-level patterns in sex-specific mortality

In addition to producing mortality rates for subnational areas, the proposed
model specification contains a number of parameters that are useful in
observing broad patterns in mortality at the state-level. Specifically,
changes in $\mu_{\beta}$ describe broad, state-level trends in mortality
patterns over time, and the correlation matrices offer insight into county-
level trends.

In Figures 6 and 7, we plot posterior medians and 95% credible intervals of
the estimated state-level coefficients for the first and second principal
components for male and female populations in all states and years. States are
organized roughly along their relative geographic location. The most notable
finding is that in the plot for the first principal component, the gap in the
baseline mortality coefficients between men and women is shrinking over time.
This suggests a convergence in mortality across the sexes (Seligman,
Greenberg, and Tuljapurkar 2016). Geographically, we see that baseline
mortality patterns differ substantially by state, suggesting evidence for a
divergence in mortality across US states due to stagnation or regression in
some states and continuous improvement in others (Fenelon 2013). For example,
states such as California and New Jersey experienced significant declines in
baseline mortality in both sexes. Others such as Louisiana experienced limited
declines and have stagnated in the 2010s. Some states such as Ohio saw a
deterioration in baseline mortality in the 2010s. This trend reversal occurs
across the eastern half of the US.

Figure 6: Posterior medians and 95% credible intervals for the state-level
coefficients for the first principal component ($\mu_{\beta}$). Figure 7:
Posterior medians and 95% credible intervals for the second principal
component ($\mu_{\beta}$).

For the second principal component, declines are experienced across the US in
both sexes with female values transitioning from positive to negative. As the
second principal component contains an inverted accident hump shape, this
suggests increasing young adult mortality relative to the baseline. Our
findings suggest that this trend accelerated in the 2010s for males in New
Jersey, West Virginia, Ohio, and more broadly across the Northeast and
Midwest. This is consistent with the sharp increase in opioid overdose deaths
in these regions of the US (Alexander, Kiang, and Barbieri 2018).

For the final two principal components (shown in Appendix B) there are again
clear geographic clusters of patterns over time. For the third principal
component, which has relatively higher mortality in mid adult ages, trends are
generally declining, but stagnating more in the South, and there is evidence
of a trend reversal in Kentucky and West Virginia. Unlike the rest of the
states, Alaska does not show any significant gap between the two sexes in the
third principal component, suggesting that the deviation to baseline mortality
associated with this component contributes similarly to male and female
mortality. However, it is important to recognize that the population in Alaska
is small, leading to large uncertainty in the extracted trends.

For the fourth principal component, a transition from positive to negative
coefficient values is found across the US. The interpretation of this trend is
similar to that of the second principal component. An inverted fourth
principal component implies relatively higher early adult mortality, which is
again consistent with the opioid crisis.

Figure 8 plots the posterior values of the between-sex principal component
coefficient correlations extracted from
$L^{(\beta)}_{i,t}L_{i,t}^{(\beta)\top}$ for a subset of states over time.
Plots for all states are given in Appendix B. Each series includes the
associated uncertainty intervals along with a horizontal line at zero. Across
states, the first principal component has high correlation between the two
sexes, reinforcing the notion that baseline mortality for both sexes move
together. In contrast, there is a limited relationship in the second principal
component. This is consistent with the idea that the accident hump is a
predominantly male phenomenon, which would imply weak correlation between
sexes. The third and fourth principal component correlations deviate by state.
For Alaska, the relationships are weak at best, for California and Texas,
there are strong or growing correlations, and New Jersey has declining
correlations.

Figure 8: Time-series of posterior principal component correlations for seven
states.

In Figure 9, we plot the posterior medians and 95% credible intervals for the
between sexes correlations in the overdispersion terms captured by
$L_{a,t}^{(\gamma)}L_{a,t}^{(\gamma)\top}$ for California as an example. We
find that patterns vary substantially by age group. For the 85+ category,
consistent high correlation is found. However, in the youth age groups, there
does not appear to be any meaningful correlation. The most interesting
patterns are those in the mid adult age groups where transitions from limited
correlation to strong positive correlation are found.

Figure 9: Posterior medians and 95% credible intervals for $\gamma$
correlation over time by age group for California.

###  4.2 Validation

To formally evaluate the model, we perform out-of-sample validation exercises
comparing the present model with the model proposed by Alexander, Zagheni, and
Barbieri (2017). By comparing these two models, we can focus on the
contributions of the between-subpopulation correlation matrices introduced.
The validation exercises focus on five states (Alaska, California, Louisiana,
New Jersey, and Texas) that capture the diversity in population size and
number of counties across the US.

For each state, we leave out 20% of the observed data for each county in the
years 1982-2019 as an out-of-sample dataset, then execute the model on the
remaining in-sample dataset. We then generate a distribution of deaths for all
left-out observations using the posterior log-mortality rates for the
corresponding county-age-sex-years. Finally, we calculate coverage at the 80%,
90%, and 95% nominal levels for the death uncertainty intervals as well as
mean squared errors (MSE) and mean absolute deviations (MAD) between the
median death estimates and the observed deaths. In Table 2, the results of
this exercise by state are presented. Note that the model introduced in this
paper is labeled as the “joint” model and the version in Alexander, Zagheni,
and Barbieri (2017) with independent sex modelling is denoted “independent”.

Table 2: Out-of-Sample Coverage and Errors State | Model | Cov80 | Cov90 | Cov95 | MAD | MSE
---|---|---|---|---|---|---
Alaska | independent | 0.884 | 0.943 | 0.973 | 1.931 | 14.447
Alaska | joint | 0.887 | 0.946 | 0.974 | 1.899 | 12.614
California | independent | 0.843 | 0.920 | 0.961 | 9.747 | 2004.650
California | joint | 0.846 | 0.923 | 0.963 | 8.568 | 1194.914
Louisiana | independent | 0.873 | 0.939 | 0.970 | 3.054 | 40.795
Louisiana | joint | 0.874 | 0.940 | 0.972 | 2.921 | 30.257
New Jersey | independent | 0.843 | 0.925 | 0.963 | 8.388 | 298.082
New Jersey | joint | 0.850 | 0.930 | 0.964 | 7.838 | 238.458
Texas | independent | 0.878 | 0.940 | 0.970 | 3.232 | 92.449
Texas | joint | 0.879 | 0.941 | 0.970 | 3.030 | 61.618

The joint model consistently outperforms the independent model on the out-of-
sample set across metrics. The level of outperformance is most pronounced in
the larger states, notably California. The outperformance of the joint model
on the out-of-sample set suggests that it may have uses beyond the US county
context for estimating subnational mortality by subpopulation in jurisdictions
without complete data.

##  5 Discussion

In this paper we extended principal component-based methods to jointly
estimate subnational mortality across subpopulations. This approach leverages
the inherent structural mortality patterns associated with the individual
principal components and extracts correlations between groups to offer insight
into the joint movements of mortality trends across groups. The model centers
on a regression-based framework with four principal components that allow core
patterns in age-specific mortality to be captured. The principal component
coefficients are modelled hierarchically to allow for information exchange
within states. County-specific effects on each component are assumed to be
correlated across subgroups within a county. Our approach is validated using
both a simulation study and with out-of-sample exercises using real
subnational mortality data. We find that the proposed model is well calibrated
and that its errors are smaller than those of the independent model in out-of-
sample exercises.

We illustrate the model through estimating US county-level sex-specific
mortality rates. An investigation into the parameters of the model highlights
general trends in US mortality as well as state specific patterns. Notably,
while movements in baseline mortality tend to manifest in both sexes
simultaneously, specific features of the mortality curve such as elevated
young adult mortality appear independent by sex. These results can offer
direction to those working to reduce mortality pressures. For example, when
correlations are high, aggregate policies may be effective. However, in
jurisdictions where correlations are low or zero, more targeted policies may
be necessary. State-level parameter estimates also highlighted clear
geographic clustering of mortality patterns over time. Notably, there is
evidence of stagnation in parts of the country, with increasing early adult
mortality particularly in the eastern states, and a stagnating improvement in
middle-age mortality in the southern states.

We identify two directions for future extensions to this work. The first is to
study how the model performs in real-world circumstances similar to the
simulation study where there are more than two subpopulations and the
composition of the population is not approximately balanced across groups. For
example, race/ethnicity-based mortality data in US counties. Many counties are
dominated by one subgroup, and thus the frequency of low or zero death data is
significantly higher than with aggregate or sex-specific data. The second
direction is to apply the model in contexts where there are differences in
mortality data collection or death registration coverage. In many countries,
complete mortality data such as the US county data used here are not available
and death registration coverage can vary substantially by sex, geographical
area, or time (Peralta et al. 2019; Basu and Adair 2021). Typically, male
coverage is higher than female coverage. By studying relationships between
sexes in the years and regions with higher quality data, one may be able to
support estimation of female mortality rates by using the male data and the
modelled correlations.

##  6 References

reAlexander, Monica, and Leontine Alkema. 2022. “A Bayesian Cohort Component
Projection Model to Estimate Women of Reproductive Age at the Subnational
Level in Data-Sparse Settings.” _Demography_ 59 (5): 1713–37.
https://doi.org/10.1215/00703370-10216406.

preAlexander, Monica, Mathew V Kiang, and Magali Barbieri. 2018. “Trends in
Black and White Opioid Mortality in the United States, 1979–2015.”
_Epidemiology (Cambridge, Mass.)_ 29 (5): 707.

preAlexander, Monica, Emilio Zagheni, and Magali Barbieri. 2017. “A Flexible
Bayesian Model for Estimating Subnational Mortality.” _Demography_ 54\.
https://doi.org/10.1007/s13524-017-0618-7.

preArias, Elizabeth, Loraine A. Escobedo, Jocelyn Kennedy, Chunxia Fu, and
Jodi Cisewki. 2018. “U.s. Small-Area Life Expectancy Estimates Project:
Methodology and Results Summary.” _Vital and Health Statistics Series 2_ 181\.
https://www.cdc.gov/nchs/data/series/sr_02/sr02_181.pdf.

preBarford, Anna, Danny Dorling, George Davey Smith, and Mary Shaw. 2006.
“Life Expectancy: Women Now on Top Everywhere.” _Bmj_. British Medical Journal
Publishing Group.

preBasu, J. K., and T. Adair. 2021. “Have inequalities in completeness of
death registration between states in India narrowed during two decades of
civil registration system strengthening?” _Int J Equity Health_ 20 (1): 195.

preBeltrán-Sánchez, Hiram, Caleb E. Finch, and Eileen M. Crimmins. 2015.
“Twentieth Century Surge of Excess Adult Male Mortality.” _Proceedings of the
National Academy of Sciences_ 112 (29): 8993–98.
https://doi.org/10.1073/pnas.1421942112.

preBhutta, Zulfiqar A. 2016. “Mapping the geography of child mortality: a key
step in addressing disparities.” _The Lancet Global Health_ 4 (12): E877–78.
https://doi.org/10.1016/S2214-109X(16)30264-9.

preBijak, J., and J. Bryant. 2016. “Bayesian demography 250 years after
Bayes.” _Population Studies_ 70 (1): 1–19.
https://doi.org/10.1080/00324728.2015.1122826.

preClark, Samuel J. 2019. “A General Age-Specific Mortality Model with an
Example Indexed by Child Mortality or Both Child and Adult Mortality.”
_Demography_ 56 (3). https://doi.org/10.1007/s13524-019-00785-3.

preCoale, Ansley J., Paul Demeny, and Barbara Vaughan, eds. 1983. Second
Edition. Academic Press. https://doi.org/10.1016/C2013-0-07295-7.

preCrimmins, Eileen M, and Yuan S Zhang. 2019. “Aging Populations, Mortality,
and Life Expectancy.” _Annual Review of Sociology_ 45: 69–89.

prede Beer, Joop. 2012. “Smoothing and Projecting Age-Specific Probabilities
of Death by TOPALS.” _Demographic Research_ 27: 543–92.
http://www.jstor.org/stable/26349934.

preDyrting, Sigurd. 2020. “Smoothing migration intensities with P-TOPALS.”
_Demographic Research_ 43 (55): 1607–50.
https://doi.org/10.4054/DemRes.2020.43.55.

preEzzati, Ari B AND Kulkarni, Majid AND Friedman. 2008. “The Reversal of
Fortunes: Trends in County Mortality and Cross-County Mortality Disparities in
the United States.” _PLOS Medicine_ 5 (4): 1–12.
https://doi.org/10.1371/journal.pmed.0050066.

preFenelon, Andrew. 2013. “Geographic Divergence in Mortality in the United
States.” _Population and Development Review_ 39 (4): 611–34.
https://doi.org/10.1111/j.1728-4457.2013.00630.x.

preGelman, Andrew, and Donald B. Rubin. 1992. “Inference from Iterative
Simulation Using Multiple Sequences.” _Statistical Science_ 7 (4): 457–72.
http://www.jstor.org/stable/2246093.

preGjonça, Arjan, Cecilia Tomassini, and James Vaupel. 1999. “Male-Female
Differences in Mortality in the Developed World.” Rostock: Max Planck
Institute for Demographic Research.

preGodwin, Jessica, and Jon Wakefield. 2021. “Space-Time Modeling of Child
Mortality at the Admin-2 Level in a Low and Middle Income Countries Context.”
_Statistics in Medicine_ 40 (7): 1593–1638.

preGonzaga, Marcos R., and Carl P. Schmertmann. 2016. “Estimation of Mortality
Rates by Age and Sex for Small Areas with TOPALS Regression: An Application
for Brazil in 2010.” _Revista Brasileira De Estudos De População_ 33 (3):
629–52. https://doi.org/10.20947/S0102-30982016c0009.

preGutin, Iliya, and Robert A Hummer. 2021. “Social Inequality and the Future
of US Life Expectancy.” _Annual Review of Sociology_ 47: 501–20.

preHarper, Sam, Corinne A Riddell, and Nicholas B King. 2021. “Declining Life
Expectancy in the United States: Missing the Trees for the Forest.” _Annual
Review of Public Health_ 42: 381–403.

preHeligman, L., and J. H. Pollard. 1980. “The Age Pattern of Mortality.”
_Journal of the Institute of Actuaries_ 107 (1): 49–80.
https://doi.org/10.1017/S0020268100040257.

preHendi, Arun S. 2015. “Trends in US Life Expectancy Gradients: The Role of
Changing Educational Composition.” _International Journal of Epidemiology_ 44
(3): 946–55.

preHoffman, Matthew D., and Andrew Gelman. 2014. “The No-u-Turn Sampler:
Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.” _Journal of
Machine Learning Research_ 15 (47): 1593–623.
http://jmlr.org/papers/v15/hoffman14a.html.

preHuman Mortality Database. 2023. Max Planck Institute for Demographic
Research (Germany), University of California, Berkeley (USA),; French
Institute for Demographic Studies (France). http://www.mortality.org (data
downloaded on August 25, 2023).

preKhana, Lauren M. Rossen, Diba, and Margaret Warner. 2018. “A Bayesian
Spatial and Temporal Modeling Approach to Mapping Geographic Variation in
Mortality Rates for Subnational Areas with r-INLA.” _Journal of Data Science_
16 (1): 147–82.

preKindig, David A., and Erika R. Cheng. 2013. “Even as Mortality Fell in Most
US Counties, Female Mortality Nonetheless Rose in 42.8 Percent of Counties
from 1992 to 2006.” _Health Affairs_ 32 (3): 451–58.
https://doi.org/10.1377/hlthaff.2011.0892.

preLee, Ronald D., and Lawrence R. Carter. 1992. “Modeling and Forecasting u.
S. Mortality.” _Journal of the American Statistical Association_ 87 (419):
659–71. http://www.jstor.org/stable/2290201.

preLewandowski, Daniel, Dorota Kurowicka, and Harry Joe. 2009. “Generating
Random Correlation Matrices Based on Vines and Extended Onion Method.”
_Journal of Multivariate Analysis_ 100 (9): 1989–2001.
https://doi.org/https://doi.org/10.1016/j.jmva.2009.04.008.

preMasters, Ryan, and Laudan Aron. 2021. “Changes in US Life Expectancy in the
Wake of COVID-19: Differences by Race/Ethnicity and Relative to Other High-
Income Countries.”

preMercer, Laina D., Jon Wakefield, Athena Pantazis, Angelina M. Lutambi,
Honorati Masanja, and Samuel Clark. 2015. “Space-Time Smoothing of Complex
Survey Data: Small Area Estimation for Child Mortality.” _The Annals of
Applied Statistics_ 9 (December): 1889–1905.
https://doi.org/10.1214/15-AOAS872.

preMontez, Jennifer Karas, Jason Beckfield, Julene Kemp Cooney, Jacob M
Grumbach, Mark D Hayward, Huseyin Zeyd Koytak, Steven H Woolf, and Anna
Zajacova. 2020. “US State Policies, Politics, and Life Expectancy.” _The
Milbank Quarterly_ 98 (3): 668–99.

preMurray, Sandeep C AND Michaud, Christopher J. L AND Kulkarni. 2006. “Eight
Americas: Investigating Mortality Disparities Across Races, Counties, and
Race-Counties in the United States.” _PLOS Medicine_ 3 (9): 1–12.
https://doi.org/10.1371/journal.pmed.0030260.

preNeal, Radford M. 2011. “MCMC Using Hamiltonian Dynamics.” In _Handbook of
Markov Chain Monte Carlo_. Chapman; Hall/CRC. https://doi.org/10.1201/b10905.

preNoymer, Andrew, and Viola Van. 2014. “Divergence without decoupling: Male
and female life expectancy usually co-move.” _Demographic Research_ 31 (51):
1503–24. https://doi.org/10.4054/DemRes.2014.31.51.

prePeralta, A., J. Benach, C. Borrell, V. Espinel-Flores, L. Cash-Gibson, B.
L. Queiroz, and M. Marí-Dell’Olmo. 2019. “Evaluation of the mortality registry
in Ecuador (2001-2013) - social and geographical inequalities in completeness
and quality.” _Popul Health Metr_ 17 (1): 3.

preRaftery, Adrian E., Nevena Lalic, and Patrick Gerland. 2014. “Joint
probabilistic projection of female and male life expectancy.” _Demographic
Research_ 30 (27): 795–822. https://doi.org/10.4054/DemRes.2014.30.27.

preRau, Roland, and Carl P. Schmertmann. 2020. “Lebenserwartung auf Kreisebene
in Deutschland.” _Dtsch Arztebl International_ 117 (29-30): 493–99.
https://doi.org/10.3238/arztebl.2020.0493.

preSchmertmann, Carl P., and Marcos R. Gonzaga. 2018. “Bayesian Estimation of
Age-Specific Mortality and Life Expectancy for Small Areas With Defective
Vital Records.” _Demography_ 55 (4): 1363–88.
https://doi.org/10.1007/s13524-018-0695-2.

preSehgal, Neil Jay, Dahai Yue, Elle Pope, Ren Hao Wang, and Dylan H. Roby.
2022. “The Association Between COVID-19 Mortality and the County-Level
Partisan Divide in the United States.” _Health Affairs_ 41 (6): 853–63.
https://doi.org/10.1377/hlthaff.2022.00085.

preSeligman, Benjamin, Gabi Greenberg, and Shripad Tuljapurkar. 2016.
“Convergence in male and female life expectancy: Direction, age pattern, and
causes.” _Demographic Research_ 34 (38): 1063–74.
https://doi.org/10.4054/DemRes.2016.34.38.

preŠevčíková, Hana, and Adrian E. Raftery. 2021. “Probabilistic Projection of
Subnational Life Expectancy.” _Journal of Official Statistics_ 37 (3):
591–610. https://doi.org/doi:10.2478/jos-2021-0027.

preSong, Insang, and Hui Luan. 2022. “The Spatially and Temporally Varying
Association Between Mental Illness and Substance Use Mortality and
Unemployment: A Bayesian Analysis in the Contiguous United States, 2001–2014.”
_Applied Geography_ 140: 102664.
https://doi.org/https://doi.org/10.1016/j.apgeog.2022.102664.

preStan Development Team. 2021. “RStan: The R Interface to Stan.” https://mc-
stan.org/.

preStevenson, John M., and David R. Olson. 1993. “Methods for Analysing
County-Level Mortality Rates.” _Statistics in Medicine_ 12 (3-4): 393–401.
https://doi.org/https://doi.org/10.1002/sim.4780120320.

preUnited States Mortality Database. 2023. University of California, Berkeley
(USA). http://usa.mortality.org (data downloaded on August 25, 2023).

preWang, Haidong, Austin E. Schumacher, Carly E. Levitz, Ali H. Mokdad, and
Christopher JL Murray. 2013. “Left Behind: Widening Disparities for Males and
Females in US County Life Expectancy, 1985–2010.” _Population Health Metrics_
11 (8). https://doi.org/10.1186/1478-7954-11-8.

preWoolf, Steven H, and Heidi Schoomaker. 2019. “Life Expectancy and Mortality
Rates in the United States, 1959-2017.” _Jama_ 322 (20): 1996–2016.

p

##  Appendix A Simulation study details

To illustrate the model’s ability to estimate mortality rates and extract
mortality patterns across subgroups, we simulate population and death data for
populations composed of five subgroups of varying sizes. The simulation study
data includes population and mortality data for 10 years, 25 subnational
areas, and 5 population subgroups. For each year-area-age-subgroup, we
construct the population data as follows:

  1. 1.

The total population in year 1 in each area is set as 100,000 times the area
index, which ranges from 1-25.

  2. 2.

The total population in each area in subsequent years is computed by
increasing the total population in each area by 1% annually.

  3. 3.

Total populations in each year-area are allocated to the 19 age groups using a
jittered version of the age distribution of Los Angeles County in California.

  4. 4.

Each year-area-age population is then disaggregated into the 5 subgroups using
the following proportions: 50% in Group A, 20% in Group B, and 10% in Groups
C, D, E.

To generate corresponding mortality data for the populations in each year-
area-age-subgroup, we first construct true log-mortality curves for each year-
area-subgroup. The mortality curves are constructed using a linear combination
of two standard curves: the first representing a baseline mortality curve, and
the second containing an accident hump. Plots of the standard curves are given
in Figure 10.

Figure 10: Simulation study mortality curve components

The log-mortality curves are then constructed with the following steps:

  1. 1.

The baseline mortality curve is multiplied by a random coefficients sampled
from a 5-dimensional multivariate normal distribution with mean
$\mathbf{1}_{5}$ and a year-dependent covariance matrix. The covariance matrix
describes the dependence between the 5 subgroups. It is assigned a constant
variance of $0.1^{2}$ and one of the correlation structures plotted in Figure
4a. In years 1-3, the independent correlation matrix is used. In years 4-6,
the exchangeable correlation matrix (ie. constant off-diagonal values) is
used. In years 7-10, the unstructured correlation matrix based on state-level
observed race/ethnicity-based US mortality is used.

  2. 2.

An accident hump is added to the mortality curve by adding the accident hump
standard curve multiplied by random coefficients sampled from a 5-dimensional
multivariate normal distribution with mean $\mathbf{0}_{5}$, and covariance
matrix with a constant variance of $0.5^{2}$ and one of the correlation
matrices plotted in Figure 4a. The matrices used in each year mirror those
used for the baseline mortality standard.

Finally, observed deaths for each year-area-age-subgroup are generated from
the log-mortality curves using a Poisson likelihood with rate equal to the
corresponding population multiplied by the exponential of the corresponding
log-mortality rate.

Examples of simulated data in two subnational areas are given in Figure 11.
Points correspond to log-mortality observations (ie. the log of observed
deaths divided by population). Note that zero death observations are encoded
as -10 on the log-scale.

Figure 11: Examples of simulated log-mortality observations.

##  Appendix B Additional results

Figure 12 plots two county-years where the amount of non-zero death data
varies by sex. The observed log-mortality rates for age groups with non-zero
deaths are presented as points, the estimated state means for each sex are
denoted by the grey lines, and the male and female county-level estimates and
95% credible intervals are presented in blue and red respectively. In both
county-years, the majority of age groups had zero female deaths in the year
whereas the male data is nearly complete. As opposed to defaulting to the
female state-level means, estimates for female mortality are also informed by
the behaviour of male mortality relative to the male state-level means.

Figure 12: Examples of US counties with limited female but sufficient male
data.

Figures 13 and 14 plot maps of the posterior medians and 95% credible
intervals of the estimated state-level coefficients for the third and fourth
principal components for male and female populations for all years.

Figures 15, 16, 17, and 18 plot the posterior medians and 95% credible
intervals for principal component coefficient correlations between sexes over
time. As expected, strong correlation in the first principal component is
found in most states whereas there is negligible correlation in the second
principal component. Patterns in the third and fourth principal component
coefficient correlations vary regionally with some evidence for stronger
correlations in the fourth principal component coefficients in the Great Lakes
region and surrounding states.

Figure 13: Posterior medians and 95% credible intervals for the third
principal component ($\mu_{\beta}$). Figure 14: Posterior medians and 95%
credible intervals for the fourth principal component ($\mu_{\beta}$). Figure
15: Time-series of posterior medians and 95% credible intervals for the first
principal component correlations. Figure 16: Time-series of posterior medians
and 95% credible intervals for the second principal component correlations.
Figure 17: Time-series of posterior medians and 95% credible intervals for the
third principal component correlations. Figure 18: Time-series of posterior
medians and 95% credible intervals for the fourth principal component
correlations.


# Dissecting physics of carbon ordering in bcc iron

Sam Waseda<EMAIL_ADDRESS>Tilmann Hickel  Jörg Neugebauer  Max-Planck-
Institut für Eisenforschung, Max-Planck-Straße 1, D-40237 Düsseldorf, Germany
Julien Morthomas  Patrice Chantrenne  Michel Perez  Univ Lyon, INSA Lyon,
UCBL, MATEIS, UMR5510, 69621 Villeurbanne, France

###### Abstract

Zener ordering is a phenomenon that octahedral interstitial atoms such as
carbon occupy the same sublattice inside bcc matrix such as iron. The original
formulation relies on a mean field theory, which is still most in use today.
We employ multiple methods, such as Molecular Dynamics, Metropolis Monte
Carlo, Mean Field Theory with chemical interactions and finite temperature
effects to show that the Zener ordering for iron carbon systems is governed by
local chemical interactions and finite temperature effects and less of mean
field nature as described originally by Zener.

††preprint: APS/123-QED

##  I Introduction

Iron-carbon alloys, commonly known as steels, serve as the backbone of modern
industrial infrastructure. Understanding the effects of different carbon
concentrations and alloying elements on the mechanical properties of steel is
essential for producing materials that are strong, ductile, and corrosion-
resistant. This knowledge enables engineers to design structures that
withstand immense pressures, varying temperatures, and environmental
challenges, ensuring safety and longevity in critical applications like
bridges, skyscrapers, and pipelines.

In spite of the extensive utilization of steel, which has spanned both a
substantial temporal continuum and a significant volumetric scale within human
history, full understanding of the energetics as well as the kinetics of
carbon inside steels has yet not been attained [1]. This deficiency of
comprehension arises from the intricate interplay of chemistry, elasticity,
and magnetism within the context of iron-carbon as well as carbon-carbon
interactions in steels. In ferrite, the anisotropic nature of the three types
of sublattice makes certain sublattice alignments of carbon atoms elastically
favourable, while it is unfavourable in terms of entropy. Such an alignment is
today commonly known as Zener ordering. The understanding of the
thermodynamics of Zener ordering will give us insight into the stability of
martensite, the kinetics of carbon atoms at the austenite-ferrite transition,
as well as the interplay with defect segregation [1].

In its original form, the theory of Zener [2] described the ordering
phenomenon on the basis of a mean-field theory (MFT), counting the number of
carbon atoms in each of the three sublattices without considering direct
interactions between the carbon atoms as well as the temperature effect beyond
configurational entropy. Khachaturyan, in his Microscopic Elasticity Theory
(MET) [3, 4, 5] included pairwise carbon-carbon interactions, which are
repulsive in nature. Waseda _et al._ [6] compared the ordering phenomenon
described by MFT, Metropolis Monte Carlo (MMC) and Molecular Dynamics (MD) and
found that short-range carbon-carbon interactions hinder Zener ordering. The
combination of the two opposing effects led to the unexpected congruence
between the results of MFT and MD at higher temperature [7, 8, 9, 10]. In
addition, their results [6] showed pronounced importance of precipitate-like
phases, such as Fe16C2. This precipitate structure is commonly observed in
systems like iron-nitrogen, but has not been detected for iron-carbon
experimentally. The result is also supported by other computational studies
that emphasized the importance of short range carbon-carbon interactions
leading to the formation of carbides [11]. These results, however, failed to
show the full analysis of finite temperature effects beyond configurational
entropy, which was speculated to be the reason why in other works a good
agreement between MD and MFT was observed [9, 7]. Moreover, according to the
MFT it is expected that there is an energy barrier between the ordered and
disorderes states, meaning the free energy minimization scheme of MMC might
converge to the local minimum and not the true free energy minimum.

In this study, we utilized the metadynamics method [12] to understand both
long- and short-range interactions of carbon atoms for temperatures ranging
from very high values employed in previous studies [9, 7], to below room
temperature, which cannot be handled by the conventional MD, as well as
varying the carbon content, in order to understand the order-disorder
transition. The interpretation of the results is substantiated by the
combination of these other methods: MFT alone, MFT with short range carbon-
carbon interactions (MFT+CC), MFT+CC with temperature dependent thermodynamic
parameters (MFT+CC+FT), conjugated with previous results [6] employing MMC and
classical MD. In order to develop the logic, we first start with the simplest
method, namely MFT, and add physical complexity to progressively converge
towards the outcomes encompassing the complete intricacies afforded by
metadynamics.

##  II Theory

###  II.1 Metadynamics

Metadynamics is an innovative approach in computational science that enhances
the sampling of complex molecular systems. It achieves this by dynamically
adding an artificial bias to the system’s energy landscape during MD
simulations. As the simulation progresses, metadynamics continuously
introduces small Gaussian potentials, or “hills”, onto the energy surface at
specific collective variables $Z$ of interest. These hills counteract the
system’s tendency to remain trapped in local energy minima, encouraging it to
explore different regions of the energy landscape. Over time, these added
potentials create a “history” that guides the system towards previously
unvisited configurations.

For the collective variable $Z$, we introduce the artificial bias $B(Z)$ via

$\displaystyle B(Z)=w\sum_{i}e^{-\frac{(Z-Z_{i})^{2}}{2\sigma^{2}}},$ (1)

where $Z_{i}$ are the measured collective variables in the previous steps, $w$
is the potential height and $\sigma$ is the smearing parameter, determining
the flatness of the bias potential. Therefore, the bias force $f_{a}$ for each
atom $a$ is calculated by

$\displaystyle f_{a}=-\frac{\partial B}{\partial Z}\frac{\partial Z}{\partial
x_{a}}.$ (2)

By gradually reducing energy barriers and promoting transitions between
different states, metadynamics allows for the exploration of rare and critical
events that might otherwise be inaccessible within the limited timescales of
traditional simulations.

The method has been heavily employed for protein folding [13, 14, 15] in
biophysics, but also in materials science for example to investigate phase
transformations [16].

###  II.2 Order parameter

In a system of $n_{\mathrm{C}}$ carbon atoms, we define the order parameter
$Z$ via

$\displaystyle Z=$
$\displaystyle\sqrt{\frac{3}{2}\frac{\sum_{i}n_{i}^{2}}{n_{\mathrm{C}}^{2}}-\frac{1}{2}},$
(3) $\displaystyle\text{where}\quad n_{i}=$
$\displaystyle\frac{\sum_{ap}\cos^{2}\left(\sum_{k}M_{pik}x_{ak}\right)}{\sum_{apj}\cos^{2}\left(\sum_{k}M_{pjk}x_{ak}\right)}.$
(4)

$n_{i}$ are to be interpreted as the number of carbon atoms occupying the sub-
lattice type $i$ ($i=x,y,z$) and $x_{ak}$ is the $k$-coordinate of the carbon
atom $a$ ($a=1,\dots,n_{\mathrm{C}}$). A site type is determined by its
orientation, i.e. the direction along which the nearest iron atoms lie. We
define the rotation matrices $M_{pjk}$ by

$\displaystyle M_{0jk}=\frac{\pi}{a_{0}}\begin{pmatrix}0&1&1\\\ 1&0&1\\\
1&1&0\end{pmatrix},\quad M_{1jk}=\frac{\pi}{a_{0}}\begin{pmatrix}0&1&-1\\\
-1&0&1\\\ 1&-1&0\end{pmatrix}$ (5)

where $a_{0}$ is the lattice parameter. This definition of $Z$ allows us to
calculate the derivative with respect to the position of each atom. Therefore
by using the order parameter $Z$ for the collective variable for metadynamics,
we can calculate the bias forces following eq. 2. This definition of the order
parameter is equivalent to the one used in [6], where the difference comes
from the fact that in their formulation the carbon atoms were placed only in
perfect octahedral sites and there was only one preferential site (i.e.
$n_{z}>n_{x}=n_{y}$). This choice of the order parameter intrinsically
captures the periodicity and the symmetry of the sublattices and therefore
allows us to efficiently obtain the biased energy and forces on the fly. In
the following, we consider the system to be ordered if $Z>0.5$ and otherwise
disordered.

###  II.3 Mean field theory (MFT) and extension of pair-wise chemical
interactions (MFT+CC)

Following MFT, the potential energy $U$ in a system of volume $V$ is given by

$\displaystyle
U_{\mathrm{elast}}=-\frac{\Delta\lambda^{2}(1+\nu)}{2VE}\sum_{i}n_{i}^{2},$
(6)

where $\Delta\lambda=\lambda^{\prime}-\lambda$ is the difference of dipole
moments [17, 18, 19, 20] along the direction of orientation (i.e. direction of
the nearest neighboring iron atoms represented by $\lambda$) and along the two
other directions represented by $\lambda^{\prime}$, $\nu$ is the Poisson’s
ratio and $E$ is the Young’s modulus. In our model, we also append the short
range carbon-carbon interactions by

$\displaystyle U_{\mathrm{chem}}=\frac{1}{2}\sum_{ij}\phi_{i}K_{ij}\phi_{j}$
(7)

where $\phi_{i}$ is the occupancy interstitial sites by C of the site $i$
($i=1,\dots,N_{\mathrm{sites}}=3N_{\mathrm{Fe}}$) and $K_{ij}$ is the chemical
interaction tensor, which is defined by

$\displaystyle
K_{ij}=\kappa_{ij}+c\kappa_{ij}^{\prime}+c^{2}\kappa_{ij}^{\prime\prime}.$ (8)

where $\kappa_{ij}$, $\kappa_{ij}^{\prime}$ and $\kappa_{ij}^{\prime\prime}$
are fitting parameters. In this context, we consider $\kappa_{ij}$ as the main
pairwise interaction parameters (simple pairwise) and $\kappa_{ij}^{\prime}$
and $\kappa_{ij}^{\prime\prime}$ as correction terms (effective manybody). The
entropy of the system $S$ is given by:

$\displaystyle
S=-k_{\mathrm{B}}\sum_{i}(\phi_{i}\ln\phi_{i}+(1-\phi_{i})\ln(1-\phi_{i}))$
(9)

The free energy $F$ is given by:

$\displaystyle F=U-TS$ (10)

where $U$ is only $U_{\mathrm{elast}}$ in the framework of the mean field
model, and $U_{\mathrm{elast}}+U_{\mathrm{chem}}$ to account for the short
range chemical interactions. Within the framework of the mean field model, the
free energy can be reduced to the function of a single variable $z$ and can be
solved analytically. In the classical Zener model, no temperature dependence
on the thermodynamic quantities are considered. However, as shown previously
[6], the temperature dependence is proposed to play a significant role in the
ordering behaviour. Therefore we look into the temperature dependence of the
dipole tensor, atomic volume, Young’s modulus and Poisson’s ratio. On the
other hand, the carbon concentration dependence of these thermodynamic
parameters was previously found to play no significant role [21] and therefore
we omit it in this study.

###  II.4 Simulation setup

The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [22],
as implemented in the integrated development environment (IDE) pyiron [23] was
employed to conduct MD simulations, and for the measurement of thermodynamic
parameters. The modified Becquart interatomic potential [24, 25] was chosen to
accurately represent iron-carbon interactions. This potential has been
employed previously for the investigation of Zener ordering [7, 6] and is
known to reproduce qualitatively the ordering behaviour measured
experimentally [26]. The simulation box contained 2000 iron atoms, with carbon
concentrations ranging from 1% to 11%. Temperatures spanned 100 K to 1000 K,
covering a wide thermal range, where the box shape was allowed to vary to keep
the box free of pressure.

To explore the impact of initial configurations, 20 simulation boxes were
utilized per temperature. Ten of these boxes began with ordered atomic
arrangements, while the other ten had disordered configurations, where it was
made sure that no two atoms are within a distance of one lattice parameter in
order to avoid an unrealistic starting condition. Simulations were conducted
until the convergence of the free energy histogram was reached, but at least
for 1 ns to capture thermal effects accurately.

For metadynamics histogram parameters, we employed the smearing value of
$\sigma=0.01$. The histogram $B(Z)$ was updated every 100 steps with the value
of $w=0.1$ meV. The histogram was stored in the form of a 1-D mesh grid with
1,000 points, meaning within one $\sigma$, there were 10 grid points. At each
integration step, the metadynamics forces on the carbon atoms are calculated
according to eq. 2. In order to prevent the collective motion of the system,
the counteracting forces are added to the iron atoms. In particular, the sum
of all metadynamics forces divided by the number of iron atoms is added to the
force of each of the iron atoms.

In parallel, we also ran conventional MD simulations for the temperature range
of 600 K to 1,400 K for 1 ns and calculated the order parameter based on the
final configuration.

The temperature dependence of the thermodynamic parameters was measured by
running MD calculations for various temperatures for 100 ps. The necessary
values (pressures, energy etc.) were measured by taking the average of all
steps after the initial 10 ps.

##  III Results

###  III.1 Temperature effects on thermodynamic parameters

In order to include the vibrational entropy into MFT, we first studied the
temperature effects on the following thermodynamic parameters: difference of
dipole moments $\Delta\lambda$, Young’s modulus, thermal expansion and
Poisson’s ratio.

Figure 1: Temperature dependence of thermodynamic parameters normalized to the
value at 300 K. Experimental results were taken from [27] for the Poisson’s
ratio, [28] for the Young’s modulus and [29] for the atomic volume.

Fig. 1 displays the variation of thermodynamic parameters as a function of the
temperature. The comparison with the experimental results shows a good
agreement. It can be clearly seen that the dipole moment hardly changes with
the temperature. There is some slight change for the Poisson’s ratio, but as
the change enters with the factor of $1+\nu$ in eq. 6, this temperature change
as well as the deviations from experiment can be safely ignored, as it lies
around a value of 0.3. The temperature dependence of the Young’s modulus, on
the other hand, has an important impact on the hardness of the material, and
we can expect the elastic energy to have a higher weight as the temperature
increases compared to the case where the 0 K Young’s modulus is used for the
entire temperature spectrum. This effect is somewhat weakened by the increase
in volume with the increasing temperature.

Figure 2: Free energy as a function of order parameter $Z$ for 7.5 at.%C at
900 K with temperature dependence for Poisson’s ratio $\nu$, difference of
dipole moments $\Delta\lambda$, Young’s modulus $E$ and volume $V$
individually, as well as no temperature dependence (i.e., values measured for
0 K) and all of them together.

Fig. 2 shows the free energy variation as a function of the order parameter
for 7.5 at.%C at 900 K. For the four parameters that have been investigated,
we can see that only the Young’s modulus has a considerable effect on the free
energy. This figure also shows that the softening of iron becomes so important
that the system becomes ordered by considering it, as the free energy minimum
falls in the region $Z>0.5$ in this case.

###  III.2 Short range chemical interactions

Next, in order to understand the short-range carbon-carbon interactions, we
inserted two carbon atoms in a box with 2,000 iron atoms for all possible
configurations, for which a structure optimization was performed to obtain the
minimum energy. The interaction between two carbon atoms $i$ and $j$ is
calculated via

$\displaystyle\kappa_{ij}=E(C_{i},C_{j})-2E(C)+E(\varnothing)+\Delta E_{ij},$
(11)

where $E(C_{i},C_{j})$ is the total energy of the box containing two carbon
atoms at positions $i$ and $j$, $E(C)$ is the total energy of the box with one
carbon atom and $E(\varnothing)$ is the box energy without carbon atoms. The
additional value $\Delta E_{ij}$ is calculated by (cf. eq. 6)

$\displaystyle\Delta
E_{ij}=\begin{cases}\frac{\Delta\lambda^{2}(1+\nu)}{VE},&\text{if same
orientation}\\\ 0,&\text{otherwise}\end{cases}$ (12)

$\Delta E_{ij}$ is needed to counteract double counting of long-range elastic
interactions. Two carbon atoms share the same orientation, when their
respective nearest iron atoms lie along the same direction.

Figure 3: Chemical interactions between two carbon atoms as a function of the
distance. Distances corresponding to the Fe16C2 configuration are marked with
circles additionally. Positive values are repulsive interactions.

As previously also shown [30, 31], the direct carbon-carbon interactions are
mostly repulsive, as shown in fig. 3. There are a few exceptional points, most
notable for the ones which would correspond to the Fe16C2 configuration. As a
matter of fact, the interaction between two carbon atoms is most attractive
for this potential, when two carbon atoms sit at a distance of one lattice
parameter in each direction, which is the shortest distance between two carbon
atoms in the Fe16C2 configuration.

In order to compute $\kappa_{ij}^{\prime}$ and $\kappa_{ij}^{\prime\prime}$ in
eq. 8, we distributed carbon atoms via simple Monte Carlo using $\kappa_{ij}$
for various temperatures and carbon contents, measured the energies and fitted
them using the Ridge regression by minimizing the following cost function:

$\displaystyle\mathcal{L}=\sum_{k}(U(\phi_{k})-\sum_{ij}\phi_{k,i}K_{ij}\phi_{k,j})^{2}+\lambda\sum_{ij}(\kappa_{ij}^{\prime
2}+\kappa_{ij}^{\prime\prime 2}),$ (13)

where

$\displaystyle
U(\phi)=E(\phi)-n_{\mathrm{C}}E(C)+(n_{\mathrm{C}}-1)E(\varnothing)$ (14)

and $\phi_{k}$ is the occupancy vector (i.e. filled with 0 and 1, where 0 is
empty and 1 is filled), $\lambda$ is the regularization parameter of the Ridge
regression, and $E(\phi)$ is the box energy of the occupancy $\phi$. The cost
function $\mathcal{L}$ is minimized with respect to $\kappa_{ij}^{\prime}$ and
$\kappa_{ij}^{\prime\prime}$. In order to see the quality of fitting, we
created a data set of 2,000 points, where we set 25 % of the points aside for
the test set. We estimated the value of $\lambda$ using the Leave-One-Out
(LOO) using the remaining 75 % of the data points.

Figure 4: Energy prediction using the pairwise interactions without correction
(simple pairwise, without $\kappa_{ij}^{\prime}$ and
$\kappa_{ij}^{\prime\prime}$) and with correction (effective manybody, with
$\kappa_{ij}^{\prime}$ and $\kappa_{ij}^{\prime\prime}$) for the test set

As it turned out, the superposition of the pairwise interactions with only
$\kappa_{ij}$ for boxes containing more than 2 carbon atoms the energy was
generally underestimated (_cf._ fig. 4). Despite the lack of explicit many-
body interactions, the energies with more than 2 atoms can be correctly
reproduced with the correction terms.

MFT vs. MFT+CC vs. MMC
---
Figure 5: Comparison of MMC results and MFT+CC results (one with simple
pairwise interactions, the other with the effective many-body interactions) at
8 at.%C. For a comparison, the MFT and MMC results [6] are also shown.

Using these chemical interactions, we obtained the order-disorder transition
behaviour at 8 %at.C as shown in fig. 5 for MFT+CC. The results are compared
against the classical MFT and the results from MMC [6], which intrinsically
has long and short range carbon-carbon interactions, but does not have the
temperature effects beyond configurational entropy. As it was speculated (cf.
[6]), the addition of direct carbon-carbon interactions decreases the
transition temperature at high carbon concentration in MFT+CC compared to MFT,
and it shows a significantly better agreement with the MMC results, which also
contain direct carbon-carbon interactions. The simple pairwise interactions
gave ambiguous behaviour, staying somewhere between ordering and disordering
at lower temperatures, and only slowly converging to the disordered state
around the same temperature as MFT. We have observed that there are
configurations for two carbon atoms which deliver particularly low energies.
These values got corrected with effective many-body interactions.

MFT vs. MFT+CC vs. MMC @ 1 at.%C
---
Figure 6: Comparison of MMC results [6], MFT+CC and MFT at 1 at.%C

Previous results [6] had indicated the presence of the Fe16C2 structure, which
was predicted to be the reason why MMC showed a higher transition temperature
than MFT. In our study, we could observe the same behaviour at 1 %at.C (fig.
6). In order to better understand the sublattice distribution, we plotted the
site radial distribution in fig. 7, which we defined by:

Figure 7: Site radial distribution as defined in eq. 15 for MFT+CC and for 1
%at.C at 175 K

$\displaystyle\rho(r_{k})=\sum_{i>j}\frac{\sqrt{\phi(r_{i})\phi(r_{j})}}{n_{k}|\phi|}\delta(r_{k}-|r_{i}-r_{j}|)$
(15)

where $n_{k}$ is the number of site pairs for the shell $k$ and $|\phi|$ is
the average occupation. We could indeed observe that at the transition
temperature predicted by MMC, there is a strong preference for the sites
corresponding to Fe16C2. As a matter of fact, we can also observe that for all
other pairs, there is no strong preference, regardless of whether the site
types are the same or not, indicating without Fe16C2, the system could very
well be disordered.

Figure 8: Comparison of MD results and MFT+CC+FT at 8 %at.C for the order
parameter as a function of the temperature

We now add the temperature dependence of thermodynamic parameters in MFT+CC to
make MFT+CC+FT and compare the results against classical MD results in fig 8,
which intrinsically have both long and short range carbon-carbon interactions
and all temperature effects. Firstly, we can observe the significant increase
of the transition temperature in both MFT+CC+FT and MD, compared to the
previous results for MFT+CC and MMC in fig. 5. Quite crucially, the simple
pairwise interactions deliver a much lower order parameter in the area, in
which the carbon atoms are ordered according to the MD results. This behaviour
is corrected in the effective many-body interactions, following the MD results
almost perfectly.

###  III.3 Metadynamics and comparison of all methods

Finally, we compare the foregoing results with the ones obtained from
metadynamics, which are comparable to MD results, but can explore the free
energy surface along the order parameter for any temperature and any carbon
content.

MFT vs. Metadynamics @ 8 at.%C
---
Figure 9: Free energy measured with MFT and metadynamics with 8 at.%C,
relative to the respective minimum.

Fig. 9 shows the free energy surface for 8.0 at.%C across the transition
temperature for metadynamics, and also for MFT for a comparison. In both cases
we can see a clear first order behaviour with qualitatively similar values,
except for the lower end of the minimum, which is $z=0$ for MFT and somewhere
around $z=0.25$ for metadynamics at higher temperature (1100 K for MFT and
1200 K for Metadynamics). It could be related to the fact that in MFT there
are only discrete octahedral sites to take, while in metadynamics any position
of any C atoms has its contribution (cf. eq. 4), so that thermal vibrations of
C atoms lead to a slight deviation from $z=0$. Similarly, Metadynamics has a
strong increase at higher order parameter compared to MFT.

Figure 10: Phase diagram of order-disorder transition for Metadynamics
(background and phase boundary interpolated by a magenta line) in comparison
with other methods discussed in previous sections

Finally, we show the phase diagram delivered by metadynamics in fig. 10 with
the transition temperature curves from previous sections. We can divide this
diagram in three parts: low carbon range ($<$2 at.%C), intermediate range
($>$2 at.%C & $<$7 at.%C) and high carbon range ($>$7 at.%C). In the low
carbon range, models including carbon-carbon interactions (MFT+CC, MFT+CC+FT,
MMC and metadynamics) show ordering behaviour, which cannot be explained
within the framework of MFT. Then in the intermediate range, there is a
suppression of ordering, which is observed in all models containing short
range carbon-carbon interactions, i.e. all but MFT. Then at high carbon range,
there is a strong increase of the transition temperature, which is observed by
models containing vibrational entropy (MFT+CC+FT, MD and metadynamics), which
is suppressed in MMC and MFT+CC. MFT, on the other hand, goes through the
middle of the phase diagram, by coincidence reaching more or less the same
area as the results from metadynamics, as the combination of all effects
roughly nullifies the overall effect. We found that with metadynamics the
transition temperature at very low C range significantly higher than with MMC
or MFT+CC+FT. It may contain some effects that we did not consider in this
article, but one reason for this is because at lower temperatures the
fluctuations in the energy landscape were a lot larger than the energy
difference between ordered and disordered states, so that it was extremely
difficult to have fully converged results. Therefore we may see a different
result if we increase the number of calculations significantly.

##  IV Conclusion and discussion

In this article, we did a comparison of order-disorder transition behaviour
using different methods, each of which captures certain thermodynamic
properties others do not offer. In the initial phase of this study, our focus
was directed towards a comprehensive examination of the thermodynamic
parameters that constitute the foundational framework of the mean-field theory
(MFT). Notably, this investigation revealed a pronounced modulation in the
Young’s modulus as a function of temperature. This observation implies that
the material exhibits a propensity for softening at elevated temperatures. In
contrast, all other thermodynamic parameters entering MFT showed
small/negligible changes in temperature. Including short-range carbon-carbon
interactions in MFT in form of pairwise interactions (MFT+CC) lowered the
order-disorder transition temperature. This tendency is in agreement with the
previous MMC results, which also contain short-range carbon-carbon
interactions. By including the explicit temperature dependence of the
thermodynamic parameters (MFT+CC+FT), we could show a good agreement with the
MD results in transition temperature and explain the difference between MMC
and MD quantitatively. Finally, we employed metadynamics to explore all areas
of the phase diagram and showed a good agreement with MFT+CC+FT.

The findings of this work have the following important implications: Firstly,
it shows that the order-disorder transition is difficult to observe in
supersaturated iron until very high temperatures are reached, since the direct
carbon-carbon interactions suppress the ordering. This might indeed be the
reason why there are not many experimental results that attest to the actual
order-disorder transition. Secondly, while we did not handle the issue in this
article, we could observe a strong dependence of the transition behaviour on
the short-range carbon-carbon interactions by varying their values slightly.
This is because long-range elastic interactions are a lot weaker than the
short-range interactions, and therefore it is possible that by employing a
different model (e.g. DFT etc.) it forms a carbide rather than converging to a
Zener ordering, which was indeed also been witnessed in previous studies [32].
At the same time, we could also see that both linear elasticity and pairwise
interactions are limited to delivering qualitative understanding when it comes
to intricate interactions of carbon atoms, as they are both outstandingly non-
linear and of many body nature, in addition to the finite temperature effects
shown in this article. Lastly, the combination of these factors, namely the
sensitivity to the temperature dependence and many-body nature of the carbon
interactions is inherently very difficult to address to using conventional
atomistic methods such as DFT, especially if more factors such as other
chemical elements or presence of defects have to be taken into account.

## Supplementary material

All materials used to produce this article can be found on this page:
https://github.com/samwaseda/zener_ordering

## Acknowledgement

We gratefully acknowledge the financial support from the German Research
Foundation (DFG) under grant HI 1300/15-1 within the DFG-ANR project C-TRAM.

## References

  * Zhang  _et al._ [2020] X. Zhang, H. Wang, T. Hickel, J. Rogal, Y. Li, and J. Neugebauer, Mechanism of collective interstitial ordering in Fe-C alloys, Nature Materials 19, 849 (2020).
  * Zener [1948] C. Zener, Theory of strain interaction of solute atoms, Physical Review 74, 639 (1948).
  * Khachaturian and Shatalov [1955] A. Khachaturian and G. Shatalov, K teorii uporyadocheniya atomov ugleroda v kristalle martensita, Fiz. Met. Metalloved. 32, 5 (1955).
  * Kurdjumov and Khachaturyan [1972] G. Kurdjumov and A. Khachaturyan, Phenomena of carbon atom redistribution in martensite, Metallurgical Transactions 3, 1069 (1972).
  * Udyansky  _et al._ [2009] A. Udyansky, J. Von Pezold, V. Bugaev, M. Friák, and J. Neugebauer, Interplay between long-range elastic and short-range chemical interactions in Fe-C martensite formation, Physical Review B 79, 224112 (2009).
  * Waseda  _et al._ [2018] O. Waseda, J. Morthomas, F. Ribeiro, P. Chantrenne, C. W. Sinclair, and M. Perez, Ordering of carbon in highly supersaturated $\alpha$-Fe, Modelling and Simulation in Materials Science and Engineering 27, 015005 (2018).
  * Sinclair  _et al._ [2010] C. Sinclair, M. Perez, R. Veiga, and A. Weck, Ordering of carbon in highly supersaturated $\alpha$-Fe: A molecular dynamics study, Physical Review B: Condensed Matter and Materials Physics (1998-2015) 81, 224204 (2010).
  * Maugis  _et al._ [2018] P. Maugis, S. Chentouf, and D. Connétable, Stress-controlled carbon diffusion channeling in bct-iron: A mean-field theory, Journal of Alloys and Compounds 769, 1121 (2018).
  * Maugis  _et al._ [2017] P. Maugis, F. Danoix, H. Zapolsky, S. Cazottes, and M. Gouné, Temperature hysteresis of the order-disorder transition in carbon-supersaturated $\alpha$-Fe, Physical Review B 96, 214104 (2017).
  * Chirkov  _et al._ [2015] P. Chirkov, A. Mirzoev, and D. Mirzaev, Molecular-dynamics simulations of carbon ordering in bcc Fe and its impact on martensite transition, Materials Today: Proceedings 2, S553 (2015).
  * Naraghi  _et al._ [2014] R. Naraghi, M. Selleby, and J. Ågren, Thermodynamics of stable and metastable structures in Fe-C system, Calphad 46, 148 (2014).
  * Theodoropoulos  _et al._ [2000] C. Theodoropoulos, Y.-H. Qian, and I. G. Kevrekidis, “coarse” stability and bifurcation analysis using time-steppers: A reaction-diffusion example, Proceedings of the National Academy of Sciences 97, 9840 (2000).
  * Barducci  _et al._ [2011] A. Barducci, M. Bonomi, and M. Parrinello, Metadynamics, Wiley Interdisciplinary Reviews: Computational Molecular Science 1, 826 (2011).
  * Granata  _et al._ [2013] D. Granata, C. Camilloni, M. Vendruscolo, and A. Laio, Characterization of the free-energy landscapes of proteins by nmr-guided metadynamics, Proceedings of the National Academy of Sciences 110, 6817 (2013).
  * Mendels  _et al._ [2018] D. Mendels, G. Piccini, Z. F. Brotzakis, Y. I. Yang, and M. Parrinello, Folding a small protein using harmonic linear discriminant analysis, The Journal of chemical physics 149 (2018).
  * Rogal  _et al._ [2019] J. Rogal, E. Schneider, and M. E. Tuckerman, Neural-network-based path collective variables for enhanced sampling of phase transformations, Physical Review Letters 123, 245701 (2019).
  * Siems [1968] R. Siems, Mechanical interactions of point defects, physica status solidi (b) 30, 645 (1968).
  * Teodosiu [2013] C. Teodosiu,  _Elastic models of crystal defects_ (Springer Science & Business Media, 2013).
  * Bacon  _et al._ [1980] D. Bacon, D. Barnett, and R. O. Scattergood, Anisotropic continuum theory of lattice defects, Progress in Materials Science 23, 51 (1980).
  * Leibfried and Breuer [2006] G. Leibfried and N. Breuer,  _Point defects in metals I: introduction to the theory_ , Vol. 81 (Springer, 2006).
  * Janßen  _et al._ [2016] J. Janßen, N. Gunkelmann, and H. M. Urbassek, Influence of c concentration on elastic moduli of $\alpha$’-Fe1-x Cx alloys, Philosophical Magazine 96, 1448 (2016).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, _et al._ , Lammps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Computer Physics Communications 271, 108171 (2022).
  * Janssen  _et al._ [2019] J. Janssen, S. Surendralal, Y. Lysogorskiy, M. Todorova, T. Hickel, R. Drautz, and J. Neugebauer, pyiron: An integrated development environment for computational materials science, Computational Materials Science 163, 24 (2019).
  * Becquart  _et al._ [2007] C. Becquart, J.-M. Raulot, G. Bencteux, C. Domain, M. Perez, S. Garruchet, and H. Nguyen, Atomistic modeling of an Fe system with a small concentration of C, Computational Materials Science 40, 119 (2007).
  * Veiga  _et al._ [2014] R. G. Veiga, C. Becquart, and M. Perez, Comments on “atomistic modeling of an fe system with a small concentration of C”, Computational materials science 82, 118 (2014).
  * Xiao  _et al._ [1995] L. Xiao, Z. Fan, Z. Jinxiu, Z. Mingxing, K. Mokuang, and G. Zhenqi, Lattice-parameter variation with carbon content of martensite. i. x-ray-diffraction experimental study, Physical Review B 52, 9970 (1995).
  * Batista Abreu  _et al._ [2014] J. C. Batista Abreu, L. Vieira, M. H. Abu-Hamd, and B. W. Schafer, development of performance-based fire design for cold-formed steel, Fire Science Reviews 3, 1 (2014).
  * Fukuhara and Sanpei [1993] M. Fukuhara and A. Sanpei, Elastic moduli and internal friction of low carbon and stainless steels as a function of temperature, ISIJ international 33, 508 (1993).
  * Denand  _et al._ [2020] B. Denand, V. Esin, M. Dehmas, G. Geandier, S. Denis, T. Sourmail, and E. Aeby-Gautier, Carbon content evolution in austenite during austenitization studied by in situ synchrotron x-ray diffraction of a hypoeutectoid steel, Materialia 10, 100664 (2020).
  * Ruban [2014] A. V. Ruban, Self-trapping of carbon atoms in $\alpha$’-Fe during the martensitic transformation: A qualitative picture from ab initio calculations, Physical Review B 90, 144106 (2014).
  * Domain  _et al._ [2004] C. Domain, C. Becquart, and J. Foct, Ab initio study of foreign interstitial atom (C, N) interactions with intrinsic point defects in $\alpha$-Fe, Physical Review B 69, 144112 (2004).
  * Huang and Maugis [2023] L. Huang and P. Maugis, Atomistic simulation of the collective carbon motion in body-centered tetragonal iron: A new insight into the martensite ageing, Acta Materialia 249, 118846 (2023).