#  SurfFlow: high-throughput surface energy calculations for arbitrary
crystals

Firat Yalcin  Computational Materials Physics, University of Vienna,
Kolingasse 14-16, 1090, Vienna, Austria  Michael Wolloch  Vasp Software GmbH,
Sensengasse 8/12, 1090, Vienna, Austria

###### Abstract

We introduce SurfFlow, an open-source high-throughput workflow package
designed for automated first-principles calculations of surface energies in
arbitrary crystals. Our package offers a comprehensive solution capable of
handling multi-element crystals, nonstoichiometric compositions, and
asymmetric slabs, for all potential terminations. To streamline the
computational process, SurfFlow employs an efficient pre-screening method that
discards surfaces with suspected high surface energy before conducting
resource-intensive density functional theory computations. The results
generated are seamlessly compiled into an optimade-compliant database,
ensuring easy access and compatibility. Additionally, a user-friendly web
interface facilitates workflow submission and management, provides result
visualization, and enables the examination of Wulff shapes. SurfFlow
represents a valuable tool for researchers looking to explore surface energies
and their implications in a diverse range of systems.

Keywords— high-throughput, surface energies, density functional theory, Wulff
construction

##  1 Introduction

The surface energy of a material determines its stability and influences its
adhesion properties, catalytic ability, and ability to form thin films and
interfaces. It is defined as the work needed to create a new surface from a
bulk crystal. It is a fundamental property that can be used to explain the
stability of surface facets, the Wulff shape [1], which is the equilibrium
crystal shape of the material, and phenomena such as surface reconstruction
[2, 3], segregation [4, 5, 6], and catalysis [7, 8, 9]. Catalytic properties
are closely tied to surface energy, and nanomaterials with high surface
energies show exceptional properties for electrocatalysis, photocatalysis, and
gas sensor applications [10]. The surface energies of a pair of materials are
also an important descriptor for the adhesion energy of these materials, a
property of great importance for the construction of solid-state batteries
[11, 12].

Experimental data on surface energies are scarce because of the technical
difficulty of the measurements. The available data are primarily limited to
specific facets of elemental crystals. Experimental measurements of surface
energies of solids can be based e.g. on cleavage [13, 14], interfacial energy
of small elastic smooth spheres [15], or contact angle measurements of various
liquids [16, 17]. The main problem with cleavage-based methods is plastic
deformation at the crack tip, while the second approach strictly works only
for amorphous solids. Contact angle measurements, on the other hand, are
associated with several difficulties and limitations that affect the accuracy
of results, such as surface contamination, roughness, and assumptions in the
underlying mathematical models.

In the quest to find new and highly functional materials, expensive and time-
consuming experimental studies have been supplemented more and more by large-
scale computational high-throughput (HT) screenings in recent years, mainly
employing density functional theory (DFT). Due to the difficulties with
conducting experiments to determine surface energies, and their crucial
importance for a range of applications (e.g. Refs. [18, 19]), a couple of HT
tools to tackle this problem have been published in recent years [20, 21, 22,
23]. Calculating surface properties via ab-initio atomistic modeling
alleviates the problem of isolating specific facets and provides exceptional
control of the system parameters.

Tran et al. generated a database of the surface energies of elemental crystals
[20], in which they provide surface energies of more than 100 polymorphs of
approximately 70 elements. Yang et al. provide an open-source code [21] to
generate organic surfaces from bulk molecular crystals, which was very
recently updated to better interface with DFT and neural network interatomic
potentials, among other improvements [24]. Brlec et al. provide a python
library [22] that automates the surface cleaving and the processing of raw DFT
outputs to extract materials properties such as surface energy. Furthermore,
there has been an interest in predicting surface properties from unrelaxed
surfaces, skipping DFT relaxations, and using a neural network model to learn
and predict cleavage energies and Wulff constructions [23].

Although these studies are effective in calculating surface energies for
specific systems, they are hampered by some constraints: the ability to handle
multi-element crystals is very limited, and/or they are not fully automated,
requiring the user to manually initiate calculations and run post-processing
scripts. This hinders the high-throughput screening approach as manual
handling of such a large number of systems quickly becomes impractical.

The main problem that so far eluded a comprehensive HT treatment of surface
energies of multi-element crystals is to correctly and automatically handle
asymmetric and/or nonstoichiometric slabs, which are unavoidable for most
systems, while still keeping the computational effort tractable. This is
further complicated by the sheer amount of symmetrically in-equivalent surface
directions (Miller indices), which might have several unique terminations
each. Looking at surfaces with Miller indices of only up to a maximal index of
3, a rather simple crystal with only two distinct elements might have on the
order of 100 unique surface configurations, while only a small fraction of
those will show favorable surface energy and actually contribute to the Wulff
shape.

In the present work, we aim to solve these issues and present a fully
automatic approach to calculate surface energies and Wulff shapes for
arbitrary crystals. Our code package efficiently handles multiple terminations
of asymmetric and nonstoichiometric slabs. Additionally, it can filter out
polar surfaces111Polar surfaces usually have higher surface energies compared
to non-polar ones and usually do not contribute to the Wulff shape. The
surface energy of polar surfaces can be calculated with SurfFlow if wanted,
however, as discussed in section S3 of the SM., and predict low energy ones,
as well as provide fully automatic calculation handling, error correction,
database operations, and output visualization.

The package is developed in Python 3 and builds on well-known and commonly
used tools for material science and HT computations like pymatgen [25],
atomate [26], and Fireworks [27]. It is available from the Python package
index (PyPI) or a GitLab repository as open source software and uses the
Vienna Ab-initio Simulation Package, vasp [28, 29, 30], to perform the DFT
calculations. See Sec. 3 for the settings used. While the methods and
algorithms used in our workflow are not individually new, we believe that no
other freely available Python package offers a comparable richness of
features, ease of use, and efficiency.

In the following sections, we will present our methodology, highlight our
measures to optimize slab size and pre-screen for low-energy surfaces using
bond valence information, and benchmark some results against available
experimental and calculated data. A detailed flowchart and explanation of the
workflow architecture can be found in S5 of the SM.

##  2 Results & Discussion

###  2.1 High-throughput handling of arbitrary surfaces

To model a surface with a DFT code, which generally employs periodic boundary
conditions in all directions, requires the construction of a slab, a thin
slice of material separated by a vacuum region from its duplicates in one
direction. The simplest case of a surface energy ($\gamma$) calculation is for
a symmetric and stoichiometric slab, meaning that both surfaces of the slab
are equivalent and the system contains an integer number of formula units of
the primitive bulk cell of the material. Then the surface energy $\gamma$ is,

$\gamma=\frac{1}{2A}\left[E_{\mathrm{slab}}-NE_{bulk}\right]\quad,$ (1)

where $A$ is the area of the exposed surface, $E_{\mathrm{slab}}$ the total
energy of the relaxed slab, $E_{bulk}$ the energy of the bulk unit222It has
been shown that using a bulk cell with the same lateral symmetry helps to
converge the surface energy quickly with respect to slab thickness [31]., and
$N$ the number of formula units in the slab. For all other cases of symmetry
and stoichiometry, a more generalized expression can be given

$\gamma_{1}=\frac{1}{A}\left[E_{\mathrm{slab}}-\sum_{i}N_{i}\mu_{i}(p,T)\right]-\gamma_{2}\quad,$
(2)

where asymmetry leads to two distinct surface energies given by $\gamma_{1}$
and $\gamma_{2}$, and the bulk energy is replaced by a term involving chemical
potentials $\mu_{i}$, which quantifies the contributions of the missing atoms
in nonstoichiometric slabs. Because the chemical potential of a species is a
function of pressure and temperature, instead of a single value for the
surface energy, we can talk about a region of stability for the surface that
is defined by the environment.

While it is in principle possible to perform reference calculations to
determine the region of stability with respect to chemical potentials, this
gets progressively more difficult as the number of species in the system
increases. It is also not always trivial to choose a bulk or gas-phase
reference [32, 33]. We have thus decided to compute our surface energies with
respect to a complete vacuum and not make an attempt to include chemical
potentials.

There have been a number of ways presented in the literature to
computationally decouple surfaces in asymmetric cases, without relying on
chemical potentials. Notable examples are the wedge method [34, 35], energy
density method [36], surface passivation method [34, 37], twinned slab method
[38], and methods based on combinations of unrelaxed and relaxed surface
energies [39, 40]. All these approaches intend to isolate one side of the slab
by eliminating the contribution to the free energy of the other side.

In this work, we have chosen to employ the method by Tian et al., which
calculates the surface energy as the combination of cleavage (very similar to
the unrelaxed surface energy given in [39]) and relaxation energies in a way
that is highly transferable and able to deal with both asymmetric and
nonstoichiometric slabs [40]. This approach is also similar to that of Eglitis
and Vanderbilt [41] with some improvements in dealing with asymmetric
surfaces. The approach is based on the notion that a slab is first cleaved
from a bulk, and then the surfaces relax into their final shape. The cleavage
energy is defined as

$E_{\mathrm{cleavage}}=\frac{1}{2A}\left[E_{\mathrm{slab}}^{\mathrm{unrelax}}-NE_{\mathrm{bulk}}\right]\quad,$
(3)

where $E_{\mathrm{slab}}^{\mathrm{unrelax}}$ is the total energy of the
unrelaxed slab, $E_{\mathrm{bulk}}$ is the energy of the bulk reference
structure, and $N$ is the number of formula units in the slab. (We should note
that this definition of the cleavage energy differs from its more common
interpretation as the average surface energy of an asymmetric slab.)
Relaxation energy for symmetric slabs is simply given as

$E_{\mathrm{relaxation}}=\frac{1}{2A}\left[E_{\mathrm{slab}}^{\mathrm{relax}}-E_{\mathrm{slab}}^{\mathrm{unrelax}}\right]\quad,$
(4)

where $E_{\mathrm{slab}}^{\mathrm{relax}}$ is the total energy of the fully
relaxed slab.

Two key assumptions in this method allow it to treat asymmetric and
nonstoichiometric slabs. First, it is assumed that the cleavage energy
(sometimes referred to as the unrelaxed surface energy in literature) is equal
for complementary terminations (i.e. sequential layers in the infinite bulk).
This assumption follows from the idea that complementary terminations are
created simultaneously by cleaving the bulk in a single plane333In reference
[40], this is claimed to hold only for systems where the constituent species
have similar electronegativities. Based on some theoretical arguments and test
calculations we believe that this is only true for extreme differences in
electronegativity, but can be considered negligible for most systems. See the
SM, section 4, for details.. Thus, one can calculate the contribution of the
cleavage energy to the surface energy on each side of the slab separately,
provided that the slab is stoichiometric, as Eqn. 3 is valid only for
stoichiometric slabs.

The second assumption is that freezing half of the slab is enough to decouple
two surfaces of the slab and simulate a semi-infinite slab which is fulfilled
if the slab is thick enough that the middle layers are bulk-like. The
relaxation energy for asymmetric slabs can then be expressed separately for
different terminations $T$

$E_{\mathrm{relaxation}}(T)=\frac{1}{A}\left[E_{\mathrm{slab}}^{\mathrm{relax}}(T)-E_{\mathrm{slab}}^{\mathrm{unrelax}}\right]\quad.$
(5)

where this time we divide by $A$ since only one side of the slab is relaxed.

These assumptions together allow SurfFlow to compute the surface energy
$\gamma$ as

$\gamma=E_{\mathrm{cleavage}}+E_{\mathrm{relaxation}}\quad,$ (6)

where $E_{\mathrm{relaxation}}$ is either equation 4 or equation 5, depending
on the symmetry of the slab. Of course, equation 6 is dependent on the
termination $T$ in the case of asymmetric slabs.

While it is possible in theory to handle all cases of symmetry and
stoichiometry using this approach, in practice we have to impose further
constraints on the systems to be able to automate the process of calculating
surface energies. First, our computational framework deals strictly with bulk-
terminated surfaces, which means that adatoms or surface vacancies cannot be
considered.

Furthermore, we limit the nonstoichiometric systems to symmetric slabs. This
is done by modifying the bottom surfaces of slabs but is not possible for all
slabs, thus this option is turned off by default and we generate equivalent
stoichiometric, but asymmetric, slabs instead. This is possible in all cases
and all possible terminations of a bulk are covered by this approach. For
these systems, it suffices to perform two relaxation calculations (top and
bottom halves) and two static calculations (static slab and bulk reference) to
calculate both surface energies $\gamma(T_{1})$ and $\gamma(T_{2})$. This case
is depicted for an example system (KCl) in Fig. 1. For a more detailed
explanation of the symmetric but nonstoichiometric case, see the SM, section
2.

With this method [40], we are able to deal with all terminations using only a
few calculations, avoid large system sizes (as for the wedge and twinned-slab
methods), or face difficulties in generalizability (as for energy density and
passivation methods).

Figure 1: Surface energy calculation scheme for an asymmetric and
stoichiometric KCl (mp-23193) (111) slab (a) with complementary terminations
using the method by Tian et al. [40]. Calculations required are relaxations of
the initial slab with only the (b) top, and (c) bottom halves relaxed, (c)
static initial slab, and (d) static oriented unit cell.

###  2.2 Optimizing performance by slab resizing and predicting low energy
surfaces

While simple DFT calculations with a GGA or metaGGA functional can no longer
be considered expensive for systems up to around 100 atoms, slabs have usually
low symmetry and especially higher index surfaces often need slabs with large
lateral extensions. Even if a single calculation is relatively cheap, usually
there are many unique directions with several possible terminations each, even
for relatively simple crystals. Only surfaces with comparatively low energy
are formed in most experimental circumstances, so those are of considerably
more interest compared to unstable ones. It is therefore prudent to minimize
the computational load by optimizing slabs and pre-screening surfaces to
determine if it is likely that they might occur in experiments and/or
contribute to the Wulff shape.

Optimization of slab size is a simpler task but still not trivial. While
pymatgen provides great tools for generating slabs and minimizing their
lateral dimensions, the minimal thickness is set by the size of the oriented
unit cell (OUC). For high-index surfaces, this OUC might become quite large
due to the necessity for periodic boundary conditions. For the slabs, the
periodic boundary condition in the direction of the surface normal is broken
by the vacuum however, SurfFlow reduces the thickness to match the intended
layer count as accurately as possible. However, the original terminations must
be preserved, which SurfFlow can achieve by clustering atomic sites into
layers and removing appropriate chunks of them. This procedure is described in
detail in S5.2 of the supplemental material (SM).

Filtering out surfaces that are unlikely to contribute to the Wulff shape
because they should have very high surface energy is a much more complex task.
In the first step, we attempt to filter out polar surfaces which usually have
higher surface energies than their nonpolar counterparts of the same crystal
[42, 43, 44]. The polarity of the surfaces is identified by pymatgen by
guessing the most likely oxidation states and calculating the dipole moment in
the direction of the surface normal. More about polar slabs can be found in
section 3 of the SM.

Next, SurfFlow attempts to rank the remaining surfaces by quantifying the
bonds broken during slab generation before actually performing any DFT
calculations. Then, the workflow proceeds with $N$ candidates predicted to
have low surface energies. These will be used for the Wulff shape generation.
(Note that high-index and high energy surfaces can be beneficial for catalysis
applications [19], so this option can be turned off when submitting a workflow
to include all possible slabs compatible with the slab generation parameters.)
In the following paragraphs, we will present the broken bond model used by
SurfFlow in some more detail and benchmark the performance on a comprehensive
set of materials with different structures, constituents, and bonding types.

Broken bond models have been used to predict surface energies of certain
materials in the past [45, 46, 47]. However, the performance of the approach
depends heavily on the material and the types of bonds present, which leads to
better performance for some systems than for others. SurfFlow adapts this
method to apply to a wide range of materials by counting and weighing broken
bonds. This is especially important for multi-element crystals where bond
strength can vary significantly. The weights assigned to each bond should
correspond to its strength, as the energy needed to break them is the main
contributor to the surface energy.

SurfFlow makes use of bond valences that have been shown to approximate the
bond energy by Etxebarria et al. [48]. The valence of a bond between species
$i$ and $j$ is given as

$S_{ij}=exp((R_{0}-R_{ij})/b)\quad,$ (7)

where $R_{0}$ is the optimal (for these bonding partners) and $R_{ij}$ is the
realized bond length in the investigated material. The parameter $b$ measures
the softness of the bond.

Eqn. 7 is the expression most widely used for bond valence, and tables of
values for $R_{0}$ and $b$ for numerous bonds are available in the literature.
To be independent of necessarily incomplete tabulated data, SurfFlow uses the
sum of the covalent radii of species forming the bond as the ideal bond
length, and the frequently used 0.37 Å as $b$ parameter [49]444Some studies
suggest that significantly different $b$ values should be used depending on
bond type [49], but this is intractable for an HT approach aimed at
generality. However, we allow the default parameter to be changed in the
setting file.. We define the total bond valence sum (BVS) of a given site $i$
as a sum over all neighbors,

$S_{i}=\sum_{j\neq i}exp((R_{0}-R_{ij})/b)\quad.$ (8)

We consider all neighbors up to 1.2 times the largest bond length of the bulk,
however, due to the exponential decay of the bond valence, mostly nearest
neighbors contribute. The total bond valence sum of broken bonds $S$ is then
simply a sum of all the partial sums over all surface sites,

$S=\sum_{i}S_{i}\quad.$ (9)

A similar, but slightly more complex approach was recently used to estimate
the relative surface energies of homoatomic transition metal crystals with
good success [50]. However, for our more diverse data set, the method
presented here was found to be slightly more reliable.

No method based on bond breaking can differentiate between in-equivalent but
complementary terminations of an asymmetric slab, because both surfaces
originate from a single cleavage. However, the bond valence sum correlates
generally well with DFT surface energies, with a median Pearson correlation
value of 0.87. Thus, in most cases, the method can be used as an excellent
preliminary filter to weed out facets and terminations predicted to have high
surface energies.

###  2.3 Benchmarking the bond-valence model for predicting low energy
surfaces

We have tested our bond valence model for 36 materials (See table 1 in the
SM), encompassing different bonding characteristics, crystal structures, and
constituting elements, from monoatomic to three-component systems. All unique
non-polar surfaces up to a maximal Miller index (MMI) of 3 have been
considered for each material, 653 surface energies in total555For a couple of
systems we found only 4 non-polar unique surfaces, and have decided to include
polar ones with $\mathrm{MMI}\leq 2$ as well for those. Also, slabs with more
than 100 sites were disregarded to conserve computational resources. See S1 in
the SM for more details.. The median correlation coefficient across all test
systems is $r_{\mathrm{med}}=0.870$, while two systems with slightly negative
correlation reduce the average to $r_{\mathrm{avg}}=0.803$. However, 30 of the
36 materials tested show a correlation larger than 0.7, as can be seen in the
histogram in Fig. 2. The four worst performing systems are bcc Li (mp-135,
$r=-0.087$), $\gamma$-TiAl (mp-1953, $r=-0.070$), as well as bcc and hcp Fe
(mp-13, $r=0.573$; mp-136, $r=0.550$). The by far worst correlation values
(essentially 0) are observed for bcc Li and $\gamma$-TiAl. Both materials have
several surfaces that are extremely close to one another in energy. For Li, a
very soft metal, where all surface energies are low, and also the BVS values
are very close together, this is not very surprising. We thus recommend not
relying on BVS filtering when calculating surface energies of very weakly
bound materials. For TiAl, however, the situation is different and might be
rooted in difficulties in evaluating Ti $d-d$ interactions correctly in
intermetallic systems. Sang and coworkers measured the charge density
distribution of $\gamma$-TiAl with highly accurate quantitative convergent
beam electron diffraction [51]. While the qualitative agreement with full
potential DFT calculations was good, there were considerable deviations for
PAW-based calculations as we conduct them here. A more detailed discussion
about all systems with $p<0.6$ can be found in S11 of the SM.

In the inset of Fig. 2 we plot examples of the correlation between $S$ and
$\gamma$ for four systems, colored according to the bins into which they are
sorted.

To test the bond valence sum filtering approach on Wulff shapes, we use the
$N$ surfaces with the lowest bond valence sums to construct Wulff shapes. All
of the surface energies are calculated with DFT, however. We then compare the
mean absolute errors (MAE) in the area fractions of these Wulff shapes with
respect to the full DFT results up to the MMI 3. The area fractions describe
which percentage of the Wulff shape each facet contributes.

Figure 2: Main plot: Histogram sorting the investigated materials into bins
according to the Pearson correlation numbers of the bond valence sum with the
DFT surface energies. Inset: Surface energy $\gamma$ vs bond valence sum $S$
for four exemplary materials representing different success situations of the
BVS estimation of the surface energy. Colors indicate the bins of the main
plot. The shaded regions represent the 95% confidence interval for the
regression estimate.

In Fig. 4 we plot Wulff shapes for the four systems already presented in Fig.
2 calculated by DFT from the lowest $N$ BVS surfaces. These give examples of
Wulff shape evolution if more and more surfaces are considered for different
classes of BVS/$\gamma$ correlation.

In the case of CoSi2 (Fig. 4(a)), where the correlation is really good, we see
that for 3 computed surfaces, the (331) facet is quite prominently featured
alongside the dominant (111) facet and a sliver of (110). The total MAE of the
area fractions is already very low at 1.8%. Adding four more surface energy
calculations reduces the MAE a bit more by including (211), (320), and (321)
facets, where the latter two disappear in the final shape in favor of the
(311) and (310) facets that are both considerably higher in energy at $N=11$.

Figure 3: Coverage required to reach MAE threshold of 0.05 for area fraction
compared to the full Wulff shape (coverage 100%). The inset visualizes how
many $N$ lowest BVS predictions must be calculated with DFT for all materials
to find the lowest energy surface.

For LiPd (Fig. 4(b)) we have a total of 41 surfaces to consider. Again, a high
index surface appears early with the (301) facet, which shows up even at $N=3$
and is ever so slightly lower in energy than (101). At $N=9$, the (320) and
(001) facets also appear, but (320) ends up not contributing to the final
shape at $N=38$, being replaced by the (101) surface. The (103) facet shown
with considerable area fraction at $N=38$ does not appear until quite late at
$N=32$, due to its relatively high surface energy. We see that with many
surfaces close in energy, even systems with decent correlations between BVS
and $\gamma$ might not produce the correct Wulff shape unless many surfaces
are considered. MAEs are therefore higher here, at 9.4% for $N=9$ and 7.2% for
$N=14$. However, calculating the Wulff shape by computing all unique surfaces
with Miller indices smaller than 1 or 2 (as commonly done) results in even
higher errors of 14 and 13.2%, and at costs of 7 and 13 calculations,
respectively.

For Al3Pt2 (Fig. 4(c)), where the correlation is in the bottom six of our test
set, low-index facets dominate the Wulff shape. Here, MAE starts off
relatively high at $N=3$ with 11.2% and as we add more surfaces, drops to 6.9%
at $N=13$ mostly due to the (100) facet appearing which also contributes to
the final shape. At $N=27$, we see contributions from (111) and (2-12) which
remain in the final Wulff shape, and at this value of $N$, the MAE drops to
3.6%. Calculating up to MMI 1 or 2, however, results in MAE of 5.2 and 0.0%,
although at relatively high costs of 11 and 23 calculations, respectively. For
this system, this ends up being the better approach, providing higher accuracy
at fewer calculations.

For $\gamma$-TiAl (Fig. 4(d)), where we have no correlation and do not expect
this method to work well, we expect similar problems. Indeed SurfFlow misses
the lowest energy surface, (110), which also contributes to the Wulff shape if
we set $N<15$. However, the second lowest energy surface, (101) is immediately
captured, as is the high energy (001) facet. The MAEs for $N=3$ and $N=9$ are
decent at 9.1 and 6.8%, and at $N=15$, we have the correct Wulff shape, even
though there are 29 distinct surfaces for $\mathrm{MMI}=3$. However, computing
up to MMI 1 or 2 is more advantageous here as well, with MAE of 3.2 and 2.4%,
at costs of 5 and 12 calculations respectively.

Overall, we can confidently say that utilizing the BVS is better than just
relying on low-index surfaces, although low surface energy alone does not
determine that a facet contributes to the Wulff shape. The average MAE
(compared to the full Wulff shape for MMI 3) for all materials is 17.1% if an
MMI of 1 is used and 12.5% if MMI is 2. Using the same number of calculations
per material as for MMI 1, or 2, but using the BVS predictions to select them,
the average MAE goes down to 9.0 and 2.7%, respectively.

Fig. 3 shows a histogram of materials with respect to the coverage666In this
context, coverage pertains to the ratio between the number of surfaces under
consideration for a given material and the total number of surfaces available
for that material. required to reach 5% of MAE when compared to the full Wulff
shape at 100% coverage. The results show that a 40% coverage is sufficient to
reach this threshold for the vast majority of materials studied ($\sim 75\%$),
while only 4 of the 36 materials need more than 60%.

Finally, the lowest-energy orientation of a crystal is also of considerable
interest. In the study of 36 materials, the BVS ranking correctly predicts the
surface with the lowest overall energy in 67% of the cases, and in over 85%,
it is within the lowest 3 predictions (see inset in Fig. 3). If finding the
lowest energy surface for a set of materials is the primary concern,
SurfFlow’s BVS filtering is thus an excellent tool to save computation time.

(a) CoSi2 (mp-2379)

(b) LiPd (mp-2744)

(c) Al3Pt2 (mp-10905)

(d) TiAl (mp-1953)

Figure 4: Wulff shapes of 4 materials with varying correlations between $S$
and $\gamma$ for different numbers of calculated Miller-index/surface shift
combinations. The facet colors correspond to the surface energy, from dark
blue (lowest) to dark red (highest).

###  2.4 Benchmarking the SurfFlow workflow

As mentioned in section 1, it is hard to measure surface energies
experimentally. Thus, it is not easy to benchmark our workflow using hard
experimental data as a reference.

We have decided to first validate our workflow using some monoatomic systems
with data from the materials project as reference [20]. The agreement is
excellent other than for some outliers which hint at some problems in the data
of reference [20] since our results match other previous studies. A detailed
discussion of these results can be found in S10 of the SM.

A more rigorous test of the important capability to treat multi-component
systems is to compare our calculated Wulff shapes with the experimentally
determined and previously calculated nanoparticle shape of the anatase and
rutile phases of TiO2 (mp-390 and mp-2657). These oxide materials potentially
feature polar surfaces and reconstructions, and thus present a particularly
hard challenge for SurfFlow’s algorithms. They are also well-studied systems,
both experimentally and especially theoretically (see e.g. the reviews by
Diebold and Liu et al. [52, 53]).

For these calculations we have employed the defaults of SurfFlow (see S5.1 in
the SM), disregarded polar surfaces, and utilized the bond valence sum for
pre-screening, only calculating the lowest 10 BVS terminations with DFT.

For anatase, SurfFlow predicts the (101) facet to be dominant (99%) in the
Wulff shape and also having the lowest surface energy of all facets at 0.41
J/m2 which is in good agreement with previous results of 0.44 J/m2 [54, 55].
The only other surface contributing to the Wulff shape is a small (001) facet,
which has a considerably higher surface energy of 0.95 J/m2 (0.90 J/m2 [54]).
Nanoparticles with this shape were already grown more than 20 years ago by
Penn and Banfield [56]. Further calculated low energy surfaces are the (100)
and (201) facets, at 0.51 J/m2 (0.53 [54]) and 0.63 J/m2, respectively.
Indeed, TiO2 anatase nanoparticles with exposed (201) facets have been grown
previously [57, 58]. The other high-index surfaces showing relatively low
surface energy we calculated are (112), (310), and (103). Another termination
of the (201) facet with low BVS, but featuring an undercoordinated oxygen atom
did not converge, so only 9 surface energies are reported in table 1. The
correlation between the BVS and the DFT computed surface energy is very high
at $r_{\mathrm{anatase}}=0.96$.

anatase | rutile
---|---
hkl | $\gamma$ | AF | hkl | $\gamma$ | AF
101 | 0.41 | 98.9 | 110 | 0.29 | 82.1
100 | 0.51 | - | 100 | 0.58 | -
201 | 0.63 | - | 311 | 0.75 | -
112 | 0.71 | - | 201 | 0.79 | 8.5
310 | 0.74 | - | 211 | 0.82 | -
210 | 0.89 | - | 331 | 0.88 | -
103 | 0.91 | - | 101 | 0.93 | -
001 | 0.95 | 1.1 | 332 | 0.97 | 9.4
110 | 1.00 | - | 320 | 1.25 | -
Table 1: Surface energies $\gamma$ in [J/m2] and Wulff shape area fractions
(AF) in % for anatase and rutile TiO2.

For rutile, (110) is found to be the dominant facet, occupying 82% of the
total area of the Wulff shape, with a surface energy of 0.29 J/m2. Two more
facets, (201) at 0.80 J/m2, and (101) at 0.97 J/m2 each contribute about 9% to
the Wulff shape. In literature, only low index surfaces (MMI=1) are
calculated, and the Wulff shapes computed feature the (110), (101), (100), and
(001) facets [59, 60]. We do calculate a low surface energy of 0.58 J/m2 for
the (100) facet, but it does not contribute to the Wulff shape, while the
(001) surface is not calculated because its BVS sum is not within the lowest
10 surfaces.

SurfFlow’s result for the (110) facet is in line with a previous PBE result of
0.31 J/m2 [54], while Perron et al. report a Perdew & Wang GGA value of 0.48
J/m2 [61], Bredow et al. find 0.63 J/m2 for the same functional and Jiang and
coworkers compute an even larger value of 0.74 J/m2 with PBE [60].

This broad spread of values is due to well-known and large oscillations of the
electronic properties of this surface with respect to the odd/even number of
titanium layers [62], which makes the surface energy very hard to converge.
Such oscillations are not found for the (100) surfaces. Our (110) result was
computed for 4 titanium layers according to SurfFlow’s defaults, which yields
a particularly low value of the surface energy, crowding out the (100) surface
from our Wulff shape. High-resolution scanning tunneling microscopy images
suggest that the (100) facet should however be present [60].

Again, one of the facets, now (320), did not converge, leaving us with 9
computed surface energies reported in table 1. Note that of the 9 computed
surface energies, 6 have $\mathrm{MMI}>1$, and the correlation between BVS and
computed $\gamma$ is decent at $r_{\mathrm{rutile}}=0.80$, indicating that
those surfaces indeed have low surface energy compared to not calculated lower
index surfaces.

In conclusion, we have developed SurfFlow, a Python package that efficiently
computes surface energies and Wulff shapes for arbitrary crystals in a high-
throughput manner. By employing symmetric considerations, optimizing slab
thickness, and implementing a pre-screening approach based on the bond valence
sum, we successfully reduce the computational cost associated with determining
low-energy surface orientations and terminations. For all but very weakly
bound materials we find a good to excellent correlation between the bond
valence sum and the surface energy. Utilizing this correlation, we have shown
that we can reliably construct more accurate Wulff shapes with fewer
calculations than by just computing low-index surfaces. This paper presents
the algorithms and approximations utilized as well as the workflow
architecture, web tools, and a comprehensive database that includes hundreds
of surface energies for crystals with diverse chemistries and structures.
Notably, SurfFlow offers features such as streamlined job submission,
automatic error correction, post-processing capabilities, and integration with
an optimade compliant MongoDB database.

##  3 Methods

We use the Vienna Ab-Initio Simulation Package (VASP, version 6.3.1, [28, 29,
30]) with the potpaw.54 potential set [63]. The potential mapping is
equivalent to the one defined in the MPRelaxSet of the Materials Project,
except for tungsten, where we use W_sv instead of the deprecated W_pv. The
user can easily overwrite this potential mapping while submitting workflows.
Our package allows for easy changes to computational and structural parameters
such as plane wave energy cutoff, k-point density, and slab thickness. We have
nevertheless taken care to set default values that promise both accurate
results and reasonable computation time. Those are e.g. using 400 eV as the
plane wave cutoff, 5.0 Å-1 for the k-point density, and slab thicknesses of at
least 8 layers or 10 Å, whatever is greater. All data presented here are
computed with those defaults, or in some cases of older data, with higher
values. Smearing parameters change for the relaxations depending on material
parameters (especially bandgap) and might be corrected during runs by
custodian, but for total energies, we always include an extra static run using
the tetrahedron method with Blöchl corrections. Our workflow is optimized to
use the PBE or SCAN functionals, but LDA or other GGAs, as well as van der
Waals corrections, can also be selected. The results calculated for this paper
are using the PBE functional. More information about the choice of functional
and pseudopotential set is available in section 8 of the SM.

## Data availability

All data described in this paper is available in an optimade compliant
database accessible via a web app at: https://surfflow-db.onrender.com/. For
more advanced users that would like to connect their MongoDB apps to browse
the database, we also provide the necessary
<EMAIL_ADDRESS>.
Additionally, we include a table with most of the results in S12 of the SM.

## Code availability

The SurfFlow code can be installed directly from PyPI with the same name, and
the source code can be found on GitHub.

## Acknowledgements

This research was funded by the Austrian Science Fund (FWF) [P 32711]. The
authors thank M. Reticcioli for fruitful discussions.

## Author contributions

Firat Yalcin: Methodology, Software, Validation, Formal analysis,
Investigation, Data curation, Writing – Original Draft, Writing – Review &
Editing, Visualization; Michael Wolloch: Conceptualization, Software,
Resources, Writing - Review & Editing, Supervision, Project administration,
Funding acquisition

## Competing interest

F.Y. declares no competing interest. M.W. is a part-time employee of the VASP
software GmbH.

## References

  * [1] G. Wulff, Zur frage der geschwindigkeit des wachstums und der auflösung der kristallformen, Zeitschrift für Kristallographie 34 (5) (1901) 449–530.
  * [2] G. A. Somorjai, Surface Reconstruction and Catalysis, Annual Review of Physical Chemistry 45 (1) (1994) 721–751.  doi:10.1146/annurev.pc.45.100194.003445.
URL https://doi.org/10.1146/annurev.pc.45.100194.003445

  * [3] K. Oura, V. Lifshits, A. Saranin, A. Zotov, M. Katayama, Surface science: An introduction, Springer Science & Business Media, 2013, [Online; accessed 2023-01-25].
  * [4] P. Wynblatt, R. C. Ku, Surface energy and solute strain energy effects in surface segregation, Surface Science 65 (2) (1977) 511–531.  doi:10.1016/0039-6028(77)90462-9.
  * [5] M. P. Seah, Quantitative prediction of surface segregation, Journal of Catalysis 57 (3) (1979) 450–457.  doi:10.1016/0021-9517(79)90011-3.
  * [6] M. Polak, L. Rubinovich, The interplay of surface segregation and atomic order in alloys, Surface Science Reports 38 (4) (2000) 127–194.  doi:https://doi.org/10.1016/S0167-5729(99)00010-2.
URL https://www.sciencedirect.com/science/article/pii/S0167572999000102

  * [7] B. Hammer, J. K. Nørskov, Theoretical surface science and catalysis—calculations and concepts, Advances in Catalysis 45 (C) (2000) 71–129.  doi:10.1016/S0360-0564(02)45013-4.
  * [8] V. Stamenkovic, B. S. Mun, K. J. Mayrhofer, P. N. Ross, N. M. Markovic, J. Rossmeisl, J. Greeley, J. K. Nørskov, Changing the activity of electrocatalysts for oxygen reduction by tuning the surface electronic structure, Angewandte Chemie - International Edition 45 (18) (2006) 2897–2901.  doi:10.1002/anie.200504386.
  * [9] J. K. Nørskov, T. Bligaard, J. Rossmeisl, C. H. Christensen, Towards the computational design of solid catalysts, Nature Chemistry 1 (1) (2009) 37–46.  doi:10.1038/nchem.121.
URL https://doi.org/10.1038/nchem.121

  * [10] Z.-Y. Zhou, N. Tian, J.-T. Li, I. Broadwell, S.-G. Sun, Nanomaterials of high surface energy with exceptional properties in catalysis and energy storage, Chem. Soc. Rev. 40 (2011) 4167–4185.  doi:10.1039/C0CS00176G.
URL http://dx.doi.org/10.1039/C0CS00176G

  * [11] I. D. Seymour, E. Quérel, R. H. Brugge, F. M. Pesci, A. Aguadero, Understanding and engineering interfacial adhesion in solid-state batteries with metallic anodes, ChemSusChem 16 (12) (2023) e202202215.  arXiv:https://chemistry-europe.onlinelibrary.wiley.com/doi/pdf/10.1002/cssc.202202215, doi:https://doi.org/10.1002/cssc.202202215.
URL https://chemistry-
europe.onlinelibrary.wiley.com/doi/abs/10.1002/cssc.202202215

  * [12] P. Restuccia, G. Losi, O. Chehaimi, M. Marsili, M. C. Righi, High-throughput first-principles prediction of interfacial adhesion energies in metal-on-metal contacts, ACS Applied Materials & Interfaces 15 (15) (2023) 19624–19633, pMID: 37015021.  doi:10.1021/acsami.3c00662.
URL https://doi.org/10.1021/acsami.3c00662

  * [13] J. W. Obreimoff, The splitting strength of mica, Proceedings of The Royal Society A: Mathematical, Physical and Engineering Sciences 127 (1930) 290–297.
  * [14] J. J. Gilman, Direct measurements of the surface energies of crystals, Journal of Applied Physics 31 (12) (1960) 2208–2218.  arXiv:https://doi.org/10.1063/1.1735524, doi:10.1063/1.1735524.
URL https://doi.org/10.1063/1.1735524

  * [15] K. Kendall, N. McN.Alford, J. D. Birchall, A new method for measuring the surface energy of solids, Nature 325 (6107) (1987) 794–796.  doi:10.1038/325794a0.
URL https://doi.org/10.1038/325794a0

  * [16] D. Y. Kwok, A. W. Neumann, Contact angle measurement and contact angle interpretation, Advances in Colloid and Interface Science 81 (3) (1999) 167–249.  doi:https://doi.org/10.1016/S0001-8686(98)00087-6.
URL https://www.sciencedirect.com/science/article/pii/S0001868698000876

  * [17] A. Kozbial, Z. Li, C. Conaway, R. McGinley, S. Dhingra, V. Vahdat, F. Zhou, B. D’Urso, H. Liu, L. Li, Study on the Surface Energy of Graphene by Contact Angle Measurements, Langmuir 30 (28) (2014) 8598–8606.  doi:10.1021/la5018328.
URL https://doi.org/10.1021/la5018328

  * [18] J. A. Williams, H. R. Le, Tribology and mems, Journal of Physics D: Applied Physics 39 (12) (2006) R201.  doi:10.1088/0022-3727/39/12/R01.
URL https://dx.doi.org/10.1088/0022-3727/39/12/R01

  * [19] C. Xiao, B.-A. Lu, P. Xue, N. Tian, Z.-Y. Zhou, X. Lin, W.-F. Lin, S.-G. Sun, High-index-facet- and high-surface-energy nanocrystals of metals and metal oxides as highly efficient catalysts, Joule 4 (12) (2020) 2562–2598.  doi:10.1016/j.joule.2020.10.002.
URL https://doi.org/10.1016/j.joule.2020.10.002

  * [20] R. Tran, Z. Xu, B. Radhakrishnan, D. Winston, W. Sun, K. A. Persson, S. P. Ong, Surface energies of elemental crystals, Scientific Data 3 (1) (2016) 160080.  doi:10.1038/sdata.2016.80.
URL https://doi.org/10.1038/sdata.2016.80

  * [21] S. Yang, I. Bier, W. Wen, J. Zhan, S. Moayedpour, N. Marom, Ogre: A Python package for molecular crystal surface generation with applications to surface energy and crystal habit prediction, Journal of Chemical Physics 152 (24) (2020).  doi:10.1063/5.0010615.
URL https://doi.org/10.1063/5.0010615

  * [22] K. Brlec, D. Davies, D. Scanlon, Surfaxe: Systematic surface calculations, Journal of Open Source Software 6 (61) (2021) 3171.  doi:10.21105/joss.03171.
  * [23] A. Palizhati, W. Zhong, K. Tran, Z. W. Ulissi, Predicting intermetallic surface energies with high-throughput dft and convolutional neural networks, ChemRxiv (2019) 1–18doi:10.26434/chemrxiv.8709566.
  * [24] S. Moayedpour, I. Bier, W. Wen, D. Dardzinski, O. Isayev, N. Marom, Structure prediction of epitaxial organic interfaces with ogre, demonstrated for tetracyanoquinodimethane (tcnq) on tetrathiafulvalene (ttf), The Journal of Physical Chemistry C 127 (21) (2023) 10398–10410.  doi:10.1021/acs.jpcc.3c02384.
URL https://doi.org/10.1021/acs.jpcc.3c02384

  * [25] S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, G. Ceder, Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis, Computational Materials Science 68 (2013) 314–319.  doi:10.1016/j.commatsci.2012.10.028.
URL https://www.sciencedirect.com/science/article/pii/S0927025612006295

  * [26] K. Mathew, J. H. Montoya, A. Faghaninia, S. Dwarakanath, M. Aykol, H. Tang, I. heng Chu, T. Smidt, B. Bocklund, M. Horton, J. Dagdelen, B. Wood, Z. K. Liu, J. Neaton, S. P. Ong, K. Persson, A. Jain, Atomate: A high-level interface to generate, execute, and analyze computational materials science workflows, Computational Materials Science 139 (2017) 140–152.  doi:10.1016/j.commatsci.2017.07.030.
URL https://www.sciencedirect.com/science/article/pii/S0927025617303919

  * [27] A. Jain, S. P. Ong, W. Chen, B. Medasani, X. Qu, M. Kocher, M. Brafman, G. Petretto, G.-M. Rignanese, G. Hautier, D. Gunter, K. A. Persson, FireWorks: a dynamic workflow system designed for high-throughput applications, Concurrency and Computation: Practice and Experience 27 (17) (2015) 5037–5059.  doi:https://doi.org/10.1002/cpe.3505.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/cpe.3505

  * [28] G. Kresse, J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47 (1993) 558–561.  doi:10.1103/PhysRevB.47.558.
URL https://link.aps.org/doi/10.1103/PhysRevB.47.558

  * [29] G. Kresse, J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6 (1) (1996) 15–50.  doi:10.1016/0927-0256(96)00008-0.
  * [30] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Physical Review B - Condensed Matter and Materials Physics 54 (16) (1996) 11169–11186.  doi:10.1103/PhysRevB.54.11169.
  * [31] W. Sun, G. Ceder, Efficient creation and convergence of surface slabs, Surface Science 617 (2013) 53–59.  doi:10.1016/j.susc.2013.05.016.
URL https://www.sciencedirect.com/science/article/pii/S003960281300160X

  * [32] T. Jia, Z. Zeng, H. Paudel, D. J. Senor, Y. Duan, First-principles study of the surface properties of $\gamma$-LiAlO2: Stability and tritium adsorption, Journal of Nuclear Materials 522 (2019) 1–10.  doi:10.1016/j.jnucmat.2019.05.007.
  * [33] T. Jia, D. J. Senor, Y. Duan, First-principles study of the surface properties of LiAl5O8: Stability and tritiated water formation, Journal of Nuclear Materials 555 (2021) 153111.  doi:10.1016/j.jnucmat.2021.153111.
  * [34] J. Zhang, Y. Zhang, K. Tse, B. Deng, H. Xu, J. Zhu, New approaches for calculating absolute surface energies of wurtzite (0001)/(000 1): A study of ZnO and GaN, Journal of Applied Physics 119 (20) (2016) 0–8.  doi:10.1063/1.4952395.
  * [35] C. Ma, W. Jin, X. Duan, X. Ma, H. Han, Z. Zhang, J. Yu, Y. Wu, From the absolute surface energy to the stabilization mechanism of high index polar surface in wurtzite structure: The case of ZnO, Journal of Alloys and Compounds 772 (2019) 482–488.  doi:10.1016/j.jallcom.2018.09.194.
URL https://doi.org/10.1016/j.jallcom.2018.09.194

  * [36] N. Chetty, R. M. Martin, First-principles energy density and its applications to selected polar surfaces, Physical Review B 45 (11) (1992) 6074–6088.  doi:10.1103/PhysRevB.45.6074.
  * [37] J. W. Kaminski, P. Kratzer, C. Ratsch, Towards a standardized setup for surface energy calculations, Physical Review B 95 (8) (2017) 1–11.  doi:10.1103/PhysRevB.95.085408.
  * [38] M. Bruno, S. Ghignone, A new computational strategy to calculate the surface energy of a dipolar crystal surface, CrystEngComm 23 (27) (2021) 4791–4798.  doi:10.1039/d1ce00403d.
  * [39] E. Heifets, R. I. Eglitis, E. A. Kotomin, J. Maier, G. Borstel, _Ab Initio_ modeling of surface structure for SrTiO 3 perovskite crystals, Physical Review B 64 (23) (2001) 235417.  doi:10.1103/PhysRevB.64.235417.
  * [40] X. Tian, T. Wang, L. Fan, Y. Wang, H. Lu, Y. Mu, A DFT based method for calculating the surface energies of asymmetric MoP facets, Applied Surface Science 427 (2018) 357–362.  doi:10.1016/j.apsusc.2017.08.172.
URL
http://dx.doi.org/10.1016/j.apsusc.2017.08.172https://doi.org/10.1016/j.apsusc.2017.08.172

  * [41] R. I. Eglitis, D. Vanderbilt, _Ab Initio_ calculations of Ba Ti O 3 and Pb Ti O 3 (001) and (011) surface structures, Physical Review B 76 (15) (2007) 155439.  doi:10.1103/PhysRevB.76.155439.
  * [42] C. Noguera, Polar oxide surfaces, Journal of Physics Condensed Matter 12 (31) (2000).  doi:10.1088/0953-8984/12/31/201.
  * [43] J. Goniakowski, F. Finocchi, C. Noguera, Polarity of oxide surfaces and nanostructures, Reports on Progress in Physics 71 (1) (2008).  doi:10.1088/0034-4885/71/1/016501.
  * [44] C. E. Dreyer, A. Janotti, C. G. Van De Walle, Absolute surface energies of polar and nonpolar planes of GaN, Physical Review B - Condensed Matter and Materials Physics 89 (8) (2014) 1–4.  doi:10.1103/PhysRevB.89.081305.
  * [45] M. Methfessel, D. Hennig, M. Scheffler, Calculated surface energies of the 4 d transition metals: A study of bond-cutting models, Applied Physics A Solids and Surfaces 55 (5) (1992) 442–448.  doi:10.1007/BF00348331.
  * [46] I. Galanakis, N. Papanikolaou, P. H. Dederichs, Applicability of the broken-bond rule to the surface energy of the fcc metals, Surface Science 511 (1-3) (2002) 1–12.  arXiv:0110236, doi:10.1016/S0039-6028(02)01547-9.
  * [47] Z. Y. Gao, W. Sun, Y. H. Hu, Mineral cleavage nature and surface energy: Anisotropic surface broken bonds consideration, Transactions of Nonferrous Metals Society of China (English Edition) 24 (9) (2014) 2930–2937.  doi:10.1016/S1003-6326(14)63428-2.
URL http://dx.doi.org/10.1016/S1003-6326(14)63428-2

  * [48] I. Etxebarria, J. M. Perez-Mato, A. García, P. Blaha, K. Schwarz, J. Rodriguez-Carvajal, Comparison of empirical bond-valence and first-principles energy calculations for a complex structural instability, Physical Review B - Condensed Matter and Materials Physics 72 (17) (2005) 1–8.  doi:10.1103/PhysRevB.72.174108.
  * [49] I. D. Brown, Recent developments in the methods and applications of the bond valence model, Chemical Reviews 109 (12) (2009) 6858–6919.  doi:10.1021/cr900053k.
  * [50] H. Ma, Y. Jiao, W. Guo, X. Liu, Y. Li, X.-D. Wen, Predicting crystal morphology using a geometric descriptor: A comparative study of elemental crystals with high-throughput dft calculations, The Journal of Physical Chemistry C 124 (29) (2020) 15920–15927.  doi:10.1021/acs.jpcc.0c03537.
URL https://doi.org/10.1021/acs.jpcc.0c03537

  * [51] X. Sang, A. Kulovits, G. Wang, J. Wiezorek, High precision electronic charge density determination for l10-ordered $\gamma$-tial by quantitative convergent beam electron diffraction, Philosophical Magazine 92 (35) (2012) 4408–4424.  doi:10.1080/14786435.2012.709324.
URL https://doi.org/10.1080/14786435.2012.709324

  * [52] U. Diebold, The surface science of titanium dioxide, Surface Science Reports 48 (5) (2003) 53–229.  doi:https://doi.org/10.1016/S0167-5729(02)00100-0.
URL https://www.sciencedirect.com/science/article/pii/S0167572902001000

  * [53] G. Liu, H. G. Yang, J. Pan, Y. Q. Yang, G. Q. M. Lu, H.-M. Cheng, Titanium dioxide crystals with tailored facets, Chemical Reviews 114 (19) (2014) 9559–9612.  doi:10.1021/cr400621z.
URL https://doi.org/10.1021/cr400621z

  * [54] M. Lazzeri, A. Vittadini, A. Selloni, Erratum: Structure and energetics of stoichiometric ${\mathrm{tio}}_{2}$ anatase surfaces [phys. rev. b 63, 155409 (2001)], Phys. Rev. B 65 (2002) 119901.  doi:10.1103/PhysRevB.65.119901.
URL https://link.aps.org/doi/10.1103/PhysRevB.65.119901

  * [55] X.-Q. Gong, A. Selloni, M. Batzill, U. Diebold, Steps on anatase tio2(101), Nature Materials 5 (8) (2006) 665–670.  doi:10.1038/nmat1695.
URL https://doi.org/10.1038/nmat1695

  * [56] R. Penn, J. F. Banfield, Morphology development and crystal growth in nanocrystalline aggregates under hydrothermal conditions: insights from titania, Geochimica et Cosmochimica Acta 63 (10) (1999) 1549–1557.  doi:https://doi.org/10.1016/S0016-7037(99)00037-X.
URL https://www.sciencedirect.com/science/article/pii/S001670379900037X

  * [57] H. B. Wu, J. S. Chen, X. W. D. Lou, H. H. Hng, Asymmetric anatase tio2 nanocrystals with exposed high-index facets and their excellent lithium storage properties, Nanoscale 3 (2011) 4082–4084.  doi:10.1039/C1NR10854A.
URL http://dx.doi.org/10.1039/C1NR10854A

  * [58] L. Wu, H. B. Jiang, F. Tian, Z. Chen, C. Sun, H. G. Yang, Ti0.89si0.11o2 single crystals bound by high-index 201 facets showing enhanced visible-light photocatalytic hydrogen evolution, Chem. Commun. 49 (2013) 2016–2018.  doi:10.1039/C3CC38105F.
URL http://dx.doi.org/10.1039/C3CC38105F

  * [59] M. Ramamoorthy, D. Vanderbilt, R. D. King-Smith, First-principles calculations of the energetics of stoichiometric ${\mathrm{tio}}_{2}$ surfaces, Phys. Rev. B 49 (1994) 16721–16727.  doi:10.1103/PhysRevB.49.16721.
URL https://link.aps.org/doi/10.1103/PhysRevB.49.16721

  * [60] F. Jiang, L. Yang, D. Zhou, G. He, J. Zhou, F. Wang, Z.-G. Chen, First-principles atomistic wulff constructions for an equilibrium rutile tio2 shape modeling, Applied Surface Science 436 (2018) 989–994.  doi:https://doi.org/10.1016/j.apsusc.2017.12.050.
URL https://www.sciencedirect.com/science/article/pii/S0169433217336371

  * [61] H. Perron, C. Domain, J. Roques, R. Drot, E. Simoni, H. Catalette, Optimisation of accurate rutile tio2 (110), (100), (101) and (001) surface models from periodic dft calculations, Theoretical Chemistry Accounts 117 (4) (2007) 565–574.  doi:10.1007/s00214-006-0189-y.
URL https://doi.org/10.1007/s00214-006-0189-y

  * [62] T. Bredow, L. Giordano, F. Cinquini, G. Pacchioni, Electronic properties of rutile $\mathrm{Ti}{\mathrm{o}}_{2}$ ultrathin films: Odd-even oscillations with the number of layers, Phys. Rev. B 70 (2004) 035419.  doi:10.1103/PhysRevB.70.035419.
URL https://link.aps.org/doi/10.1103/PhysRevB.70.035419

  * [63] G. Kresse, D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59 (1999) 1758–1775.  doi:10.1103/PhysRevB.59.1758.
URL https://link.aps.org/doi/10.1103/PhysRevB.59.1758


[a]Laura Olivera-Nieto

# Standardized Formats for Gamma-Ray Analysis Applied to HAWC Observatory Data

, Vikas Joshi  , Harm Schoorlemmer  and Axel Donath

###### Abstract

A wide range of data formats and proprietary software have traditionally been
used in $\gamma$-ray astronomy, usually developed for a single specific
mission or experiment. However, in recent years there has been an increasing
effort towards making astronomical data open and easily accessible. Within the
$\gamma$-ray community this has translated to the creation of a common data
format across different $\gamma$-ray observatories: the "gamma-astro-data-
format" (GADF). Based on a similar premise, open-source analysis packages,
such as Gammapy, are being developed and aim to provide a single, robust tool
which suits the needs of many experiments at once. In this contribution we
show that data from the High-Altitude Water Cherenkov (HAWC) observatory can
be made compatible with the GADF and present the first GADF-based production
of event lists and instrument response functions for a ground-based wide-field
instrument. We use these data products to reproduce with excellent agreement
the published HAWC Crab spectrum using Gammapy. Having a common data format
and analysis tools facilitates joint analysis between different experiments
and effective data sharing. This will be especially important for next-
generation instruments, such as the proposed Southern Wide-field Gamma-ray
Observatory (SWGO) and the planned Cherenkov Telescope Array (CTA).

##  1 Introduction

Historically, a variety of instrument-specific and largely proprietary tools
and data formats have been used in $\gamma$-ray astronomy, which hinders
effective data-sharing and reproducibility. However, in recent years there has
been a shift towards making data more accessible and easier to share in the
context of joint analysis. A big driver of this trend has been the upcoming
Cherenkov Telescope Array (CTA) [1], the data of which will become public
after a short proprietary period. Motivated by this development, there has
been an advent of openly developed analysis tools, such as Gammapy [2], which
are able to replace the existing instrument-specific packages by offering a
single, common tool. Gammapy is a community developed Python package for
$\gamma$-ray astronomy selected to be part of the CTA science tools. It has
been successfully used and validated for analysis of Imaging Atmospheric
Cherenkov Telescope (IACT) data [3] and joint analysis with data from the
Fermi Large Area Telescope [4].

In parallel to these efforts, a common data format across different
observatories, the gamma-astro-data-format (GDAF)111https://gamma-astro-data-
formats.readthedocs.io/en/latest/ [5] has been developed. The scope of this
standard is to cover all high-level data products from telescopes, starting at
the level of event lists and instrument response functions (IRFs). This format
relies on file storage by the Flexible Image Transport System (FITS) format
[6], which is widely used by the whole astronomical community. It builds on
existing standards such as
OGIP222https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/ofwg_intro.html and
expands them to tailor the specific needs of the $\gamma$-ray community. The
focus of these efforts has largely been IACT data, ignoring another type of
$\gamma$-ray experiment: wide-field ground-based observatories.

The High-Altitude Water Cherenkov (HAWC) is a $\gamma$-ray detector that
consists of 300 water Cherenkov detectors, each outfitted with four
photomultiplier tubes (PMTs). Its wide field of view covers two-thirds of the
sky uninterruptedly, allowing for constant monitoring and deep observations of
the $\gamma$-ray sky. It has been in operation since November 2014.

It is worth noting that the HAWC Accelerated Likelihood (HAL) [7] framework
and the Multi-Mission Maximum Likelihood framework (3ML) [8], the packages
primarily used by the HAWC observatory, are also open-source, but spouse a
different philosophy to that of the packages described above. Instead of
replacing the existing frameworks from different observatories by a single,
common tool, packages like 3ML provide a common framework in which the
instrument-specific tools (such as HAL) interface. This approach has the
advantage to include not only $\gamma$-ray instruments but also
multiwavelength and multimessenger observations. However, it has the
disadvantage to rely on instrument-specific analysis software, which, in the
case of $\gamma$-ray observations, performs very similar tasks between
different instruments.

##  2 HAWC data and IRFs in the GADF

Events recorded by the HAWC observatory are binned by size, that is, on the
fraction of PMTs from the array that were triggered by the event. A total of 9
event size bins are defined, as detailed in Table 2. The event size only
weakly correlates with energy [9]. In order to estimate the energy on an
event-by-event basis, more advanced algorithms have been developed. The ground
parameter (GP) algorithm is based on the charge density deposited at ground by
the shower. The neural network (NN) algorithm estimates energies with an
artificial neural network that takes as input several quantities computed
during the event reconstruction. A detailed overview of both algorithms can be
found in [10]. Energy bins are typically defined beforehand in HAWC analysis
and referred to with letter names, as described in Table 2. The combination of
event sizes and energy bins leads to a 2-dimensional bin scheme, with 108
analysis bins resulting from the combination of each event size bin 1 to 9
with the 12 energy bins. However, only a subset of these bins are populated
with enough event statistics, for example, low energy events are very unlikely
to have large event sizes.

Bin number | Low edge | High edge
---|---|---
1 | 0.067 | 0.105
2 | 0.105 | 0.162
3 | 0.162 | 0.247
4 | 0.247 | 0.356
5 | 0.356 | 0.485
6 | 0.485 | 0.618
7 | 0.618 | 0.740
8 | 0.740 | 0.840
9 | 0.840 | 1.00

Table 1: Event size ("nhit") bins. Bins are defined from the fraction of PMTs
triggered by each event.

Bin | Low edge (TeV) | High edge (TeV)
---|---|---
a | 0.316 | 0.562
b | 0.562 | 1.00
c | 1.00 | 1.78
d | 1.78 | 3.16
e | 3.16 | 5.62
f | 5.62 | 10.0
g | 10.0 | 17.8
h | 17.8 | 31.6
i | 31.6 | 56.2
j | 56.2 | 100
k | 100 | 177
l | 177 | 316

Table 2: The energy bins. Note that the first two bins are not used in the
analysis as the estimate is highly biased [10].

###  2.1 Events and good time intervals

For each of these analysis bins, gamma-hadron separation cuts are optimized
[10] and applied to the reconstructed data. Additional direction corrections
are applied and the event coordinates are converted to the J2000 epoch. The
resulting event lists are stored in FITS [6] files with headers and columns
compliant with the GADF. Additionally, other columns are stored that contain
information pertinent to the characteristics of wide-field arrays, such as
core location in the array, direction in local coordinates and event size.

Due to the large sky coverage and high duty cycle, it is common among wide-
field observatories to make sky-maps as their primary data product directly
from reconstructed events. However, producing event lists as an intermediate
step has the advantage of adding more flexibility. It allows to easily select
a subset of the whole dataset for a given analysis, which simplifies the study
of time-dependent signals. It also facilitates extensive systematic checks.

Good Time Intervals (GTIs) are defined as the time intervals during which the
detector is on and taking data continuously. These intervals are then used to
build the exposure information. Detector downtime is caused by a variety of
factors, ranging from hardware issues to meteorological conditions. These
interruptions are not uniformly distributed over time, with certain
meteorological events being more likely during particular times of the year.
This leads to fluctuations in the exposure as a function of hour angle, or,
equivalently, right ascension (R.A.). For HAWC data collected between June
2015 and June 2019, we compute the number of transits - i.e., sidereal days -
for which the detector was taking data as a function of the R.A. of zenith,
shown by the green line in Figure 1.

Figure 1: Number of transits during which data was recorded as a function of
R.A.

It can be useful, for example for background modeling, to remove this R.A.
dependency on the exposure, to "flatten" it. To do this, one can remove GTIs
from the time selection until the resulting exposure is no longer a function
of R.A. An example of this process can be seen in Figure 1 This exposure
"flattening" leads to a loss of data of around 1-2%, maintaining the overall
detector efficiency above 90%.

###  2.2 Intrument response functions

The IRFs describe the combined detection abilities and precision of an
instrument data-taking and reconstruction procedure. They are computed by
simulating a point source emitting $\gamma$-rays following a given energy
spectrum, usually $\propto$E-2. These events are processed with the detector
simulation and reconstruction procedure, see [9] for more details. The
reconstructed events are binned as described in Section 2 and gamma-hadron
separation is applied. This process yields information on the number of events
successfully identified as $\gamma$-rays as well as the accuracy and precision
of the energy and direction assigned to them. In the GADF, this information is
split into three components respectively, the effective area, the energy
dispersion matrix, and the point-spread function (PSF). Note that this
framework neglects the correlation between the different IRF components. This
is currently mostly sufficient for IACT analysis and also present in the
standard HAWC framework. However, this will be re-addressed for CTA and thus,
likely for the GADF as well in future.

####  2.2.1 Effective area

The effective area of a detector is the combination of its detection
efficiency with the observable area. Typically in the computation of HAWC IRFs
the simulated spectrum is convolved with a source transit over the
observatory, leading to a comparison of events detected to events thrown per
square meter and transit. Analysis bins described in Section 2 include a cut
on zenith angle at 45° [10]. We can compute the number of hours that a source
at each declination spends at zenith angles lower than 45°, that is, the
duration of a transit as a function of declination, to recover the effective
area in units of m2. Figure 2 shows these curves for all the available
declinations resulting from analysis bins defined using each of the energy
estimation schemes described in Section 2. The effective is greatest for
declinations close to 19°, the terrestrial longitude of the HAWC location.
Sources at this declination pass through the local zenith as they transit the
sky over the observatory. The effective area per transit, that is, the
effective exposure for one transit, is also a useful quantity, especially when
considering the long-term study of sources. Given its declination dependence,
we fill a Gammapy Map333https://docs.gammapy.org/0.18.2/maps/index.html with
the one-transit effective exposure. Then, for a given data range selection,
the number of transits like the one shown in Figure 1 is combined with this
map to produce the effective exposure map used for the analysis.

Figure 2: Effective area after background rejection for bins defined using
each of the energy estimators described in Section 2.

###  2.3 Point-Spread Function

The point-spread function is a measure of the precision achieved in the event
direction reconstruction. It is computed as the spatial probability
distribution of events produced by a point source. It is assumed to be
radially symmetric, and so can be stored as a function of offset from the
source location. Gammapy provides several options to store and use this
information in the irf444https://docs.gammapy.org/0.18.2/irf/psf.html#irf-psf
class. Like the other IRFs, the HAWC PSF depends on declination. For this
reason, we fill a PSFMap with the PSF radial profile information at each sky
position. From this, we can compute the 68% containment radius, which we
compare to the published values in [10] to validate the procedure. This
comparison can be seen in Figure 3, showing good agreement.

Figure 3: Comparison of the PSF containment at the Crab location for analysis
bins using the energy estimators described in Section 2.

###  2.4 Energy dispersion matrix

The energy dispersion is a measure of the accuracy and precision achieved in
the event energy estimation. It is computed as the probability of an event
with a given simulated energy (Etrue) to be reconstructed with a different
energy (Ereco). We use the
EDispMap555https://docs.gammapy.org/0.18.2/api/gammapy.irf.EDispMap.html#gammapy.irf.EDispMap
class in Gammapy, which is a dedicated 4-dimensional sky-map. At each sky
position, it contains the probability matrix that quantifies the energy
dispersion. An example of such matrices at the Crab location is shown in
Figure 4 for each of the event size bins defined in Table 2 and Ereco energy
bins defined in Table 2. As can be seen, the energy resolution improves with
increasing event size bin.

Figure 4: Energy dispersion matrix at the Crab location for each of the energy
estimators described in Section 2

##  3 Validation analysis with Gammapy: the Crab Nebula

The Crab Nebula is one of the brightest sources in the $\gamma$-ray sky. For
this reason, it is typically used for calibration and as a reference analysis.
The HAWC Collaboration published a measurement of the Crab spectrum extending
up to 100 TeV using the energy estimators described in [10]. Using the same
data range and background estimation method as in that work, and the IRFs as
described above, we can reproduce that result using Gammapy. We fit a combined
3-dimensional (spatial and spectral) model. For the spatial model we use a
point-source assumption, and a log-parabola for the spectral shape. The result
of this fit is shown in Figure 5, for both energy estimators. For both
figures, the top panel compares the best-fit spectrum obtained with Gammapy
with the published one, showing excellent agreement. In the bottom panel, the
flux points computed with Gammapy are compared to the model in [10].

Figure 5:  Best-fit Crab spectrum obtained with Gammapy compared with [10] for
both energy estimators described in the text. The residual shows the
comparison of the flux points with the model in [10].

##  4 Conclusion

The data from ground-based, wide-field observatories, and in this case, from
the HAWC observatory, can also be made compatible with the GDAF and thus can
be analyzed using the related open-source tools with minor adjustments. We
find excellent agreement with the results published in [10] using an analysis
tool that is built with a different philosophy and structure. This is a
powerful check of both the scientific result and both tools involved, as well
as the production of HAWC data and IRFs in the GDAF-compliant format.

Having a common data format and analysis tools facilitates joint analysis
between different experiments and effective data sharing. This synergy between
experiments is particularly relevant given the complimentary nature of
pointing and wide-field instruments. This will be specially important for
future observatories like SWGO [11].

The lifetime of observatories is finite, and one of the concerns at the end of
the operation is to ensure that the archival data is available and easy to use
for future studies and reproducibility of results. Having data in a format
that is common to other observatories and which can be analyzed with a general
use tool is a great advantage in this regard.

Gammapy has recently been selected as the official CTA Science tool. This
ensures that it will be maintained and used by the overall $\gamma$-ray
community, a much larger developer and user base than any of the other
collaboration-specific tools individually.

## Acknowledgments

We acknowledge the support from: the US National Science Foundation (NSF); the
US Department of Energy Office of High-Energy Physics; the Laboratory Directed
Research and Development (LDRD) program of Los Alamos National Laboratory;
Consejo Nacional de Ciencia y Tecnología (CONACyT), México, grants 271051,
232656, 260378, 179588, 254964, 258865, 243290, 132197, A1-S-46288,
A1-S-22784, cátedras 873, 1563, 341, 323, Red HAWC, México; DGAPA-UNAM grants
IG101320, IN111716-3, IN111419, IA102019, IN110621, IN110521; VIEP-BUAP; PIFI
2012, 2013, PROFOCIE 2014, 2015; the University of Wisconsin Alumni Research
Foundation; the Institute of Geophysics, Planetary Physics, and Signatures at
Los Alamos National Laboratory; Polish Science Centre grant,
DEC-2017/27/B/ST9/02272; Coordinación de la Investigación Científica de la
Universidad Michoacana; Royal Society - Newton Advanced Fellowship 180385;
Generalitat Valenciana, grant CIDEGENT/2018/034; Chulalongkorn University’s
CUniverse (CUAASC) grant; Coordinación General Académica e Innovación (CGAI-
UdeG), PRODEP-SEP UDG-CA-499; Institute of Cosmic Ray Research (ICRR),
University of Tokyo, H.F. acknowledges support by NASA under award number
80GSFC21M0002. We also acknowledge the significant contributions over many
years of Stefan Westerhoff, Gaurang Yodh and Arnulfo Zepeda Dominguez, all
deceased members of the HAWC collaboration. Thanks to Scott Delay, Luciano
Díaz and Eduardo Murrieta for technical support.

## References

  * [1] CTA-Consortium, _Science with the Cherenkov Telescope Array_ , WORLD SCIENTIFIC (2019), 10.1142/10986, [https://www.worldscientific.com/doi/pdf/10.1142/10986].
  * [2] C. Deil et al., _Gammapy - A prototype for the CTA science tools_ , in _35th International Cosmic Ray Conference (ICRC2017)_ , vol. 301 of _International Cosmic Ray Conference_ , p. 766, Jan., 2017 [1709.01751].
  * [3] L. Mohrmann et al., _Validation of open-source science tools and background model construction in $\gamma$-ray astronomy_, _Astronomy and Astrophysics_ 632 (2019) A72 [1910.08088].
  * [4] C. Nigro et al., _Towards open and reproducible multi-instrument analysis in gamma-ray astronomy_ , _Astronomy and Astrophysics_ 625 (2019) A10 [1903.06621].
  * [5] C. Deil et al., _Open high-level data formats and software for gamma-ray astronomy_ , in _6th International Symposium on High Energy Gamma-Ray Astronomy_ , vol. 1792 of _American Institute of Physics Conference Series_ , p. 070006, Jan., 2017, DOI [1610.01884].
  * [6] D.C. Wells et al., _FITS - a Flexible Image Transport System_ , _Astronomy and Astrophysics, Supplement_ 44 (1981) 363.
  * [7] A.U. Abeysekara et al., _Characterizing gamma-ray sources with HAL (HAWC Accelerated likelihood) and 3ML_ , in _Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021)_ , vol. 395, p. 828, 2021, DOI.
  * [8] G. Vianello et al., _The Multi-Mission Maximum Likelihood framework (3ML)_ , _arXiv e-prints_ (2015) arXiv:1507.08343 [1507.08343].
  * [9] A.U. Abeysekara et al., _Observation of the Crab Nebula with the HAWC Gamma-Ray Observatory_ , _Astrophysical Journal_ 843 (2017) 39 [1701.01778].
  * [10] A.U. Abeysekara et al., _Measurement of the Crab Nebula Spectrum Past 100 TeV with HAWC_ , _Astrophysical Journal_ 881 (2019) 134 [1905.12518].
  * [11] A. Albert et al., _Science Case for a Wide Field-of-View Very-High-Energy Gamma-Ray Observatory in the Southern Hemisphere_ , _arXiv e-prints_ (2019) arXiv:1902.08429 [1902.08429].

## Full Authors List: HAWC Collaboration

A.U. Abeysekara48, A. Albert21, R. Alfaro14, C. Alvarez41, J.D. Álvarez40,
J.R. Angeles Camacho14, J.C. Arteaga-Velázquez40, K. P. Arunbabu17, D. Avila
Rojas14, H.A. Ayala Solares28, R. Babu25, V. Baghmanyan15, A.S. Barber48, J.
Becerra Gonzalez11, E. Belmont-Moreno14, S.Y. BenZvi29, D. Berley39, C.
Brisbois39, K.S. Caballero-Mora41, T. Capistrán12, A. Carramiñana18, S.
Casanova15, O. Chaparro-Amaro3, U. Cotti40, J. Cotzomi8, S. Coutiño de León18,
E. De la Fuente46, C. de León40, L. Diaz-Cruz8, R. Diaz Hernandez18, J.C.
Díaz-Vélez46, B.L. Dingus21, M. Durocher21, M.A. DuVernois45, R.W.
Ellsworth39, K. Engel39, C. Espinoza14, K.L. Fan39, K. Fang45, M. Fernández
Alonso28, B. Fick25, H. Fleischhack51,11,52, J.L. Flores46, N.I. Fraija12, D.
Garcia14, J.A. García-González20, J. L. García-Luna46, G. García-Torales46, F.
Garfias12, G. Giacinti22, H. Goksu22, M.M. González12, J.A. Goodman39, J.P.
Harding21, S. Hernandez14, I. Herzog25, J. Hinton22, B. Hona48, D. Huang25, F.
Hueyotl-Zahuantitla41, C.M. Hui23, B. Humensky39, P. Hüntemeyer25, A.
Iriarte12, A. Jardin-Blicq22,49,50, H. Jhee43, V. Joshi7, D. Kieda48, G J.
Kunde21, S. Kunwar22, A. Lara17, J. Lee43, W.H. Lee12, D. Lennarz9, H. León
Vargas14, J. Linnemann24, A.L. Longinotti12, R. López-Coto19, G. Luis-Raya44,
J. Lundeen24, K. Malone21, V. Marandon22, O. Martinez8, I. Martinez-
Castellanos39, H. Martínez-Huerta38, J. Martínez-Castro3, J.A.J. Matthews42,
J. McEnery11, P. Miranda-Romagnoli34, J.A. Morales-Soto40, E. Moreno8, M.
Mostafá28, A. Nayerhoda15, L. Nellen13, M. Newbold48, M.U. Nisa24, R. Noriega-
Papaqui34, L. Olivera-Nieto22, N. Omodei32, A. Peisker24, Y. Pérez Araujo12,
E.G. Pérez-Pérez44, C.D. Rho43, C. Rivière39, D. Rosa-Gonzalez18, E. Ruiz-
Velasco22, J. Ryan26, H. Salazar8, F. Salesa Greus15,53, A. Sandoval14, M.
Schneider39, H. Schoorlemmer22, J. Serna-Franco14, G. Sinnis21, A.J. Smith39,
R.W. Springer48, P. Surajbali22, I. Taboada9, M. Tanner28, K. Tollefson24, I.
Torres18, R. Torres-Escobedo30, R. Turner25, F. Ureña-Mena18, L. Villaseñor8,
X. Wang25, I.J. Watson43, T. Weisgarber45, F. Werner22, E. Willox39, J.
Wood23, G.B. Yodh35, A. Zepeda4, H. Zhou30

1Barnard College, New York, NY, USA, 2Department of Chemistry and Physics,
California University of Pennsylvania, California, PA, USA, 3Centro de
Investigación en Computación, Instituto Politécnico Nacional, Ciudad de
México, México, 4Physics Department, Centro de Investigación y de Estudios
Avanzados del IPN, Ciudad de México, México, 5Colorado State University,
Physics Dept., Fort Collins, CO, USA, 6DCI-UDG, Leon, Gto, México, 7Erlangen
Centre for Astroparticle Physics, Friedrich Alexander Universität, Erlangen,
BY, Germany, 8Facultad de Ciencias Físico Matemáticas, Benemérita Universidad
Autónoma de Puebla, Puebla, México, 9School of Physics and Center for
Relativistic Astrophysics, Georgia Institute of Technology, Atlanta, GA, USA,
10School of Physics Astronomy and Computational Sciences, George Mason
University, Fairfax, VA, USA, 11NASA Goddard Space Flight Center, Greenbelt,
MD, USA, 12Instituto de Astronomía, Universidad Nacional Autónoma de México,
Ciudad de México, México, 13Instituto de Ciencias Nucleares, Universidad
Nacional Autónoma de México, Ciudad de México, México, 14Instituto de Física,
Universidad Nacional Autónoma de México, Ciudad de México, México, 15Institute
of Nuclear Physics, Polish Academy of Sciences, Krakow, Poland, 16Instituto de
Física de São Carlos, Universidade de São Paulo, São Carlos, SP, Brasil,
17Instituto de Geofísica, Universidad Nacional Autónoma de México, Ciudad de
México, México, 18Instituto Nacional de Astrofísica, Óptica y Electrónica,
Tonantzintla, Puebla, México, 19INFN Padova, Padova, Italy, 20Tecnologico de
Monterrey, Escuela de Ingeniería y Ciencias, Ave. Eugenio Garza Sada 2501,
Monterrey, N.L., 64849, México, 21Physics Division, Los Alamos National
Laboratory, Los Alamos, NM, USA, 22Max-Planck Institute for Nuclear Physics,
Heidelberg, Germany, 23NASA Marshall Space Flight Center, Astrophysics Office,
Huntsville, AL, USA, 24Department of Physics and Astronomy, Michigan State
University, East Lansing, MI, USA, 25Department of Physics, Michigan
Technological University, Houghton, MI, USA, 26Space Science Center,
University of New Hampshire, Durham, NH, USA, 27The Ohio State University at
Lima, Lima, OH, USA, 28Department of Physics, Pennsylvania State University,
University Park, PA, USA, 29Department of Physics and Astronomy, University of
Rochester, Rochester, NY, USA, 30Tsung-Dao Lee Institute and School of Physics
and Astronomy, Shanghai Jiao Tong University, Shanghai, China, 31Sungkyunkwan
University, Gyeonggi, Rep. of Korea, 32Stanford University, Stanford, CA, USA,
33Department of Physics and Astronomy, University of Alabama, Tuscaloosa, AL,
USA, 34Universidad Autónoma del Estado de Hidalgo, Pachuca, Hgo., México,
35Department of Physics and Astronomy, University of California, Irvine,
Irvine, CA, USA, 36Santa Cruz Institute for Particle Physics, University of
California, Santa Cruz, Santa Cruz, CA, USA, 37Universidad de Costa Rica, San
José , Costa Rica, 38Department of Physics and Mathematics, Universidad de
Monterrey, San Pedro Garza García, N.L., México, 39Department of Physics,
University of Maryland, College Park, MD, USA, 40Instituto de Física y
Matemáticas, Universidad Michoacana de San Nicolás de Hidalgo, Morelia,
Michoacán, México, 41FCFM-MCTP, Universidad Autónoma de Chiapas, Tuxtla
Gutiérrez, Chiapas, México, 42Department of Physics and Astronomy, University
of New Mexico, Albuquerque, NM, USA, 43University of Seoul, Seoul, Rep. of
Korea, 44Universidad Politécnica de Pachuca, Pachuca, Hgo, México,
45Department of Physics, University of Wisconsin-Madison, Madison, WI, USA,
46CUCEI, CUCEA, Universidad de Guadalajara, Guadalajara, Jalisco, México,
47Universität Würzburg, Institute for Theoretical Physics and Astrophysics,
Würzburg, Germany, 48Department of Physics and Astronomy, University of Utah,
Salt Lake City, UT, USA, 49Department of Physics, Faculty of Science,
Chulalongkorn University, Pathumwan, Bangkok 10330, Thailand, 50National
Astronomical Research Institute of Thailand (Public Organization), Don Kaeo,
MaeRim, Chiang Mai 50180, Thailand, 51Department of Physics, Catholic
University of America, Washington, DC, USA, 52Center for Research and
Exploration in Space Science and Technology, NASA/GSFC, Greenbelt, MD, USA,
53Instituto de Física Corpuscular, CSIC, Universitat de València, Paterna,
Valencia, Spain