# The suppression of magnetism and the development of superconductivity within
the collapsed tetragonal phase of Ca0.67Sr0.33Fe2As2 at high pressure

J. R. Jeffries  Condensed Matter and Materials Division, Lawrence Livermore
National Laboratory, Livermore, CA 94550, USA  N. P. Butch  Condensed Matter
and Materials Division, Lawrence Livermore National Laboratory, Livermore, CA
94550, USA  K. Kirshenbaum  Center for Nanophysics and Advanced Materials,
Department of Physics, University of Maryland, College Park, MD 20742, USA  S.
R. Saha  Center for Nanophysics and Advanced Materials, Department of Physics,
University of Maryland, College Park, MD 20742, USA  S. T. Weir  Condensed
Matter and Materials Division, Lawrence Livermore National Laboratory,
Livermore, CA 94550, USA  Y. K. Vohra  Department of Physics, University of
Alabama at Birmingham, Birmingham, Alabama 35294, USA  J. Paglione  Center for
Nanophysics and Advanced Materials, Department of Physics, University of
Maryland, College Park, MD 20742, USA

###### Abstract

Structural and electronic characterization of (Ca0.67Sr0.33)Fe2As2 has been
performed as a function of pressure up to 12 GPa using conventional and
designer diamond anvil cells. The compound (Ca0.67Sr0.33)Fe2As2 behaves
intermediate between its end members—CaFe2As2 and SrFe2As$2$—displaying a
suppression of magnetism and the onset of superconductivity. Like other
members of the AEFe2As2 family, (Ca0.67Sr0.33)Fe2As2 undergoes a pressure-
induced isostructural volume collapse, which we associate with the development
of As-As bonding across the mirror plane of the structure. This collapsed
tetragonal phase abruptly cuts off the magnetic state, giving rise to
superconductivity with a maximum $T_{c}$=22.2 K. The maximum $T_{c}$ of the
superconducting phase is not strongly correlated with any structural
parameter, but its proximity to the abrupt suppression of magnetism as well as
the volume collapse transition suggests that magnetic interactions and
structural inhomogeneity may play a role in its development. The pressure-
dependent evolution of the ordered states and crystal structures in
(Ca,Sr)Fe2As2 provides an avenue to understand the generic behavior of the
other members of the AEFe2As2 family.

superconductivity, magnetism, x-ray diffraction, pressure

###### pacs:

74.62.Fj, 74.70.Xa, 75.50.Ee, 61.50.Ks

##  I Introduction

Since the first reports of superconductivity with a critical temperature
$T_{c}$=26 K in fluorine-doped LaFeAsO,Kamihara2008  researchers have rapidly
expanded the number of Fe-based superconductors,Ishida2009  raised the $T_{c}$
to about 55 K,Ren2008  and identified five different, but related, crystal
structures in which these Fe-based superconductors crystallize.Paglione2010 ;
Johnston2010  Like the cuprate superconductors,Lee2006  the different Fe-based
superconductors display many common themes in both the electronic and
structural properties: the presence of corrugated Fe-pnictogen or Fe-chalcogen
layers within a tetragonal unit cell, and the occurrence of antiferromagnetic
order in the undoped or ambient-pressure compounds.Lumsden2010  The ubiquity
of these common elements makes these systems fertile playgrounds to explore
the interplay between magnetism, structure, and superconductivity.

One of the archetypal Fe-based superconductor structures is the “122”
structure: AEFe2X2 (ThCr2Si2-type), with AE an alkaline earth element (Ca, Sr,
Ba) an alkali metal (K, Rb, Cs) or Eu and X a pnictogen element.Paglione2010 ;
Johnston2010 ; Lumsden2010 ; Gooch2010 ; Jeevan2008  Variants of the 122
structure have been widely studied owing to the availability of a wide range
of chemical substitutions on different crystallographic sites (e.g., Co for
Fe, K for Ba, P for As, etc.) as well as their tendency to form macroscopic,
high-purity crystals. The parent compounds within the 122 systems are
paramagnetic metals at room temperature, but at low temperatures each member
of the 122 family exhibits a concomitant structural and magnetic transition.
Despite their different structural/magnetic transition temperatures—spanning a
range greater than 100 K—and chemical compositions, each of the 122 parent
compounds displays the same low-temperature structural and magnetic phases.
The tetragonal I4/mmm space group stable at room temperature undergoes a
distortion that leads to a low-temperature orthorhombic (Fmmm space group)
crystal structure, where the basal plane of the orthorhombic unit cell is
rotated by 45∘ with respect to that of the tetragonal unit cell.Rotter2008 ;
Huang2008  At ambient pressure or without doping, spin-density-wave (SDW),
antiferromagnetic (AFM) order occurs simultaneous with the tetragonal-
orthorhombic structural transition. The AFM state is characterized by a (101)
wavevector (note: the magnetic and orthorhombic unit cells are identical),
yielding Fe moments directed along the orthorhombic $a$-axis that are
antiferromagnetically arranged along $a$ and $c$ (between Fe layers) and
ferromagnetically coupled along $b$. In contrast to the wide range of AFM
ordering temperatures ($T_{N}$), the ordered moment within the 122 family
varies only slightly between 0.80 and 1.01 $\mu_{B}$.Goldman2008 ; Kaneko2008
; Su2009

With applied pressure or doping, both the structural and AFM transitions are
generally suppressed. In the case of Co or K doping in BaFe2As2, the nominally
concomitant structural and magnetic transitions separate from one another,
with the structural transition preceding the magnetic transition upon
cooling.Ni2008 ; Chu2009 ; Pratt2009 ; Urbano2010  Near the suppression of the
AFM state, with either doping or pressure, superconductivity arises with
critical temperatures ranging from roughly 9-47 K. Han2009 ; Saha2012  A
notable exception to this general observation is the lack of superconductivity
in CaFe2As2 under highly hydrostatic pressure conditions.Yu2009  In addition
to superconductivity, an isostructural volume collapse to a collapsed
tetragonal phase is seen as a function of pressure as well as for a small
minority of chemical substitutions.Saha2012 ; Goldman2009 ; Uhoya2010a ;
Uhoya2010b ; Uhoya2011 ; Mittal2011 ; Danura2011  The proximity of the
superconducting state with both structural and magnetic instabilities has
prompted suggestions that the maximum $T_{c}$ in the 122 family of compounds
could be controlled by structural parameters, magnetic interactions, or
both.Ishida2009 ; Paglione2010 ; Johnston2010 ; Lumsden2010 ; Yildrim2009
Each of these factors could have ramifications on the pairing symmetry of the
superconducting state itself,Mazin2009  and, as such, exploring the
relationships between superconductivity, magnetism, and structural
instabilities is an important component of understanding the unconventional,
high-temperature superconductivity seen in the ferropnictide compounds.

In this article we report a pressure-dependent structural and electrical
transport study of (Ca0.67Sr0.33)Fe2As2. The isoelectronic substitution of Sr
for Ca in this pseudobinary alloy expands the ambient-pressure lattice volume
and rapidly increases $T_{N}$ close to that of SrFe2As2 with with the addition
of approximately 30% Sr.Kirshenbaum2012  The effects of the larger volume and
the enhanced $T_{N}$ are to expand the phase space occupied by the AFM state
to higher temperatures and higher pressures relative to that of pure CaFe2As2,
thus pushing the destruction of magnetism to higher pressures and allowing for
a larger region of study under pressure.

##  II Experimental Details

Single crystals of (Ca0.67Sr0.33)Fe2As2 were synthesized with a flux-growth
technique previously described.Saha2009  The samples were verified with x-ray
diffraction to crystallize in the I4/mmm ThCr2Si2-type crystal structure with
ambient-pressure lattice constants $a$=3.907 Å and $c$=11.988 Å.

Pressure-dependent electrical transport measurements were performed using two
pressure cells: (i) a hydrostatic clamp cell employing n-pentane:isoamyl
alcohol as a pressure-transmitting medium was used up to 1.1 GPa; and (ii) a
designer diamond anvil cell (DAC) loaded with quasihydrostatic solid steatite
as a pressure-transmitting medium was used for pressures above 1.76 GPa. The
designer DAC was composed of a 300-$\mu$m culet, 8-probe designer diamond
anvilWeir2000 ; Patterson2000 ; Jackson2006  paired with a matching standard
diamond anvil. In order to facilitate electrical contact with the sample,
tungsten contact pads were lithographically deposited onto the microprobes
exposed at the culet of the designer diamond anvil. A non-magnetic MP35N
gasket was pre-indented to a thickness of 40 $\mu$m and a 130-$\mu$m hole was
drilled in the center of the indentation by means of an electric discharge
machine (EDM). A small, thin crystallite (approximately 70 $\times$ 70
$\times$ x 20 $\mu$m) was placed on the culet of the designer diamond anvil in
contact with the tungsten contact pads. The pressure was calibrated using the
shift in the R1 fluorescence line of ruby.Mao1986 ; Vos1991  The ruby R2
fluorescence line remained distinguishable from the R1 line to pressures in
excess of 7 GPa, implying a nearly hydrostatic environment below 7 GPa.
Temperature-dependent, electrical resistance measurements were performed in a
commercial cryostat.

For x-ray diffraction measurements, the DAC was composed of a pair of opposed
diamond anvils with 700-$\mu$m culets and a nickel gasket. The gasket was pre-
indented to a thickness of 65 $\mu$m and a 250-$\mu$m hole was drilled in the
center of the indentation with an EDM. The (Ca0.67Sr0.33)Fe2As2 crystals were
crushed in a mortar and pestle and loaded into the sample space along with a
few small ruby chips for initial pressure calibration and fine Cu powder (3-6
$\mu$m, Alfa Aesar) for in situ x-ray pressure calibration. A 4:1
methanol:ethanol mixture served as the pressure-transmitting medium.

Room-temperature, angle-dispersive x-ray diffraction (ADXD) experiments were
performed at the HPCAT beamline 16 BM-D of the Advanced Photon Source at
Argonne National Laboratory. A 5x10 $\mu$m, 29.2 keV (${\lambda}_{inc}$=0.4246
Å) incident x-ray beam, calibrated with CeO2, was used. The diffracted x-rays
were detected with a Mar345 image plate; exposure times ranged from 300-600
seconds. 2D diffraction patterns were collapsed to 1D intensity versus
2$\Theta$ plots using the program FIT2D.Hammersley1996  Pressure-dependent
lattice parameters were extracted by indexing the positions of the Bragg
reflections using the EXPGUI/GSAS package.Larson1994 ; Toby2001

##  III Results

###  III.1 Crystal structure

A typical, powder x-ray diffraction pattern for (Ca0.67Sr0.33)Fe2As2, taken at
a pressure of 1.98 GPa in the DAC, is shown in Fig. 1. The Bragg reflections
corresponding to the tetragonal I4/mmm structure of (Ca0.67Sr0.33)Fe2As2 are
indicated by the red tickmarks below the data (thin, black crosses), for which
the background has been subtracted.The green stars represent the positions of
the Bragg peaks of the Cu pressure marker. The diffraction pattern is well
described, as indicated by the absence of additional peaks in the pattern, by
including only a combination of (Ca0.67Sr0.33)Fe2As2 and Cu. The
(Ca0.67Sr0.33)Fe2As2 specimen displays a preferred orientation, with the
crystallites of the powder tending to form small platelets aligned with the
$c$-axis parallel to surface of the culet of the diamond anvil (i.e., parallel
to the incident x-ray beam). This preferred orientation results in a relative
decrease in the intensity of the (00l) reflections and an increase in the
intensity of the (hk0) reflections. Nonetheless, a full refinement (red line
through data) of the diffraction pattern results in a good fit to the data,
allowing for determination of the lattice parameters as well as the
$z$-coordinate of the As atoms.

Figure 1: (color online) An example x-ray diffraction pattern acquired at 1.98
GPa in a DAC. The refinement runs through the data points as a red line, and
the residual of the refinement is shown below the pattern as a light blue
line. Bragg reflections of the (Ca0.67Sr0.33)Fe2As2 sample are shown as red
tickmarks, while Bragg peaks from the Cu pressure marker are indicated by the
green stars.

The structural parameters extracted from refinements of the x-ray diffraction
data under pressure are shown in Fig. 2 up to 12 GPa. Ambient-pressure values
are from [Saha2011, ]. With increasing pressure, the $c$-axis of the
tetragonal unit cell monotonically decreases, but with a steeper slope between
roughly 2 and 6 GPa. The $a$-axis, on the other hand, increases with pressure
within the same 2-6 GPa range, followed by a more conventional compression for
pressures in excess of 6 GPa. The unit cell volume and the $c/a$ ratio (Fig.
2b) naturally reflect the pressure dependences of the lattice parameters, with
both quantities exhibiting an increased slope centered around 4 GPa.

Figure 2: (color online) Structural parameters extracted from refinements of
x-ray diffraction patterns under pressure: (a) tetragonal lattice parameters
$c$ (red circles, left axis) and $a$ (blue squares, right axis), (b) unit cell
volume (green diamonds, left axis) and $c/a$ ratio (orange triangles, right
axis). The inset shows the refined z-coordinate of the As site as a function
of pressure. In all cases, lines are guides to the eye.

These pressure-dependent evolution of the structural parameters shown in Fig.
2 indicate the presence of an isostructural volume collapse, identical to that
seen in the other pure, alkaline-earth 122 compounds: CaFe2As2, SrFe2As2, and
BaFe2As2. The structural parameters all exhibit inflection points near 4 GPa,
providing a consistent estimate for the volume-collapse transition pressure
(vertical, grey bar in Fig. 2) in (Ca0.67Sr0.33)Fe2As2. Above 4 GPa,
(Ca0.67Sr0.33)Fe2As2 is in the collapsed tetragonal phase. The $z$-coordinate
of the As atoms, a free parameter within the ThCr2Si2 structure, also exhibits
an anomaly near the volume-collapse transition as seen in the inset of Fig.
2b. The $z$-coordinate increases slightly at low pressures before exhibiting a
significant increase near 4 GPa. Above 4 GPa, the $z$-coordinate decreases
before increasing and recovering toward the general trend (dashed line) seen
at low pressure, suggesting a correlation between the As atoms and onset of
the collapsed tetragonal phase.

###  III.2 Electrical transport

The electrical resistivity $\rho$ as a function of temperature for selected
pressures is presented in Fig. 3. The electrical resistivity data have been
normalized such that the ambient pressure value of $\rho$(300 K) is equal to
one. In the ambient-pressure curve, the concomitant magnetic and structural
transition is evident as a pronounced jump near 200 K. With applied pressure,
$\rho$(300 K) decreases and the magnetic/structural transition is smoothly
suppressed, disappearing between 1.10 and 1.76 GPa. There is no evidence
suggesting a splitting of the structural and magnetic transitions. Within the
magnetic state at 0.87 and 1.10 GPa, there is a rapid reduction in resistivity
just below 20 K; this behavior is reminiscent of the strain-induced
superconductivity observed in pure SrFe2As2 crystals,Saha2009  although a full
resistive transition is lacking here within the magnetic state of
(Ca0.67Sr0.33)Fe2As2.

Figure 3: (color online) Normalized electrical resistivity as a function of
temperature for selected pressures (denoted in GPa unless otherwise
specified). The magnetic/structural transition and the onset of
superconductivity are visible in (a). The downward arrows ($T_{ct}$) indicate
inflection points in the electrical resistivity curves (see text). The
evolution of the superconducting transition with pressure is highlighted in
(b).

At 1.76 GPa, the lowest measured pressure where the signature of the
magnetic/structural transition is no longer visible, the electrical
resistivity displays an inflection point (downward arrows in Fig. 3a) near 80
K and full resistive superconducting transition at $T_{c}$=22.2 K. Higher
pressures reduce $T_{c}$, as clearly seen in Fig. 3b, but substantially
increase the temperature of the inflection point. Previous experiments on
CaFe2As2 revealed a similar occurrence and pressure-dependent behavior of this
inflection point.Torikachvili2008a

The features extracted from the pressure-dependent electrical resistivity
measurements are collected in the phase diagram of Fig. 4. The closed, red
squares represent the onset of magnetic order ($T_{N}$) and its associated
structural transition, while the closed, blue circles reveal the evolution of
$T_{c}$ with pressure. The open, blue circles at 0.87 and 1.10 GPa indicate
incomplete transitions possibly associated with strain-induced filamentary
superconductivity. The inflection point, seen for $P>$1.76 GPa, is shown as
open, green squares. The pressure dependence of the inflection point in the
electrical resistivity intersects with the volume collapse transition pressure
at room temperature (described in III.1), leading to the conclusion that the
inflection point seen in the temperature-dependent electrical resistivity
signifies the onset of the collapsed tetragonal phase, a conclusion consistent
with previous pressure-dependent and Rh-substitution studies of
CaFe2As2Torikachvili2008a ; Torikachvili2008b ; Yu2009 ; Canfield2009 ;
Danura2011  We thus denote this feature (i.e., the inflection point in the
electrical resistivity) as $T_{ct}$. Bulk superconductivity, as inferred from
a complete resistive transition, is seemingly limited to the collapsed
tetragonal phase.

Figure 4: (color online) Phase diagram of (Ca0.67Sr0.33)Fe2As2 showing the
suppression of magnetism ($T_{N}$ \- red squares), the development of
superconductivity ($T_{c}$ \- blue circles), and the progression of the volume
collapse transition ($T_{ct}$ \- green, crossed squares) with pressure. The
room-temperature value of $T_{ct}$ is determined from x-ray diffraction data,
all other data points are from electrical transport measurements. The open,
blue circles at lower pressures represent incomplete superconducting
transitions. Lines and shaded regions are guides to the eye.

##  IV Discussion

Given the strong link between the appearances of superconductivity and the
collapsed tetragonal phase in the 122 ferropnictide family of superconductors,
it is naturally important to explore what driving mechanisms or correlations
may be responsible for each phenomena.

###  IV.1 Isostructural volume collapse

At room temperature, the isostructural volume collapse in (Ca0.67Sr0.33)Fe2As2
occurs near 4 GPa (III.1); the volume collapse transition shifts to lower
pressures with reduced temperature (III.2). The evolution of the position of
the As atoms within the unit cell (given by the $z$-coordinate in Fig. 2) upon
passing through the volume collapse transition suggests that the As atoms are
involved in this transition.

Figure 5: (color online) Pressure dependence of the As-As distance ($d_{As-
As}$) across the mirror plane of the crystal structure shown in the inset.
Lines through the data are guides to the eye.

Figure 5 shows the interlayer As-As spacing across the mirror plane of the
unit cell, $d_{As-As}$, as a function of pressure at room temperature. From
ambient pressure, the As-As spacing decreases continuously, and nearly
linearly, with applied pressure, reaching a value of $d_{As-As}$=3.06 Å at 3.2
GPa. Between 3.2 and 3.8 GPa, $d_{As-As}$ abruptly decreases to a value of
2.94 Å, a 4% reduction occurring over 0.6 GPa. Further pressure causes a
continuous, monotonic decrease in $d_{As-As}$ up to the highest pressure
measured. The onset of the collapsed tetragonal phase in (Ca0.67Sr0.33)Fe2As2,
therefore, is signified by a collapse in the As-As separation across the
mirror plane of the unit cell. From Fig. 5, the midpoint of this collapse
occurs when $d_{As-As}$=3.0 Å.

The tendency for a collapse across the mirror plane of the ThCr2Si2-type
structure has been discussed previously. Hoffman and Zheng formulated this
collapse for BaMn2P2,Hoffman1985  but the effect can be generalized to other
compounds with this structure. For the purpose of discussion we refer to a
general formula AB2X2. The basic description of the collapse put forth by
Hoffman and Zheng is predicated on the chemistry of the B2X${}_{2}^{-2}$
layer, which yields a schematic density of states (Fig. 6a) with X-X bonding
and anti-bonding $p$-states separated by the $d$-states arising from the B
atom. As the atomic number of B increases within a row of the Periodic Table,
the Fermi level shifts downward, leaving the anti-bonding $p$-states of the
density of states unfilled, resulting in the development of an X-X bond across
the mirror plane of the structure. X-X bonding across the width of the unit
cell does not occur because that dimension is fixed by the $a$ lattice
parameter, which is at least partly set by the size of the A cation and
typically larger than 3.5 Å.Just1996  There is thus a chemical route to
describing the uncollapsed and collapsed tetragonal phases of the
ThCr2Si2-type structure, which sheds light on the mechanism under pressure.

Figure 6: (color online) (a) Schematic density of states of the
B2As${}_{2}^{-2}$ layers for a hypothetical AB2As2 compound. The Fermi level
lowers with increasing d-band occupancy, depleting the As-As anti-bonding
states, creating an As-As mirror plane bond, and collapsing the structure. (b)
Schematic unit cell energy (adapted from [Hoffman1985, ]) versus $d_{As-As}$
for a hypothetical AB2As2 specimen; pressure or d-electron element
substitution should push $d_{As-As}$ leftward, providing a driving force for
the collapsed phase.

Full electronic structure calculations by Hoffman and Zheng indicate that for
BaMn2P2 there is a maximum structural energy (Fig. 6b) when the P-P distance
is about 2.7 Å.Hoffman1985  In BaMn2P2, the P-P bond length is thus shifted
roughly 0.5 Å above the bare P-P bond length.Pyykko2009  If a similar value of
the X-X bond-length shift occurs in (Ca0.67Sr0.33)Fe2As2, then the As-As bond
length would shift from the bare As-As bond length of 2.4 Å to about 2.9 Å, in
excellent agreement with the value $d_{As-As}$=3.0 Å defining the volume
collapse transition pressure.

The onset of such As-As bonding would naturally be directed along the $c$-axis
of the unit cell (inset, Fig. 5), and would tend to pull the previously weakly
connected FeAs cages toward one another, accounting for the contraction of the
crystallographic $c$-axis upon entering the collapsed tetragonal phase. By
conservation, the increase in bonding between mirror plane As atoms would
likely reduce the Fe-As bond strength, which would relax the FeAs cages, alter
the Fe-As bond angles, and increase the $a$-axis of the unit cell as seen
experimentally.

In addition to the structural consequences of the development of this new As-
As bond within the structure, there are likely electronic structure effects.
Band structure calculations conclude that the transition into the collapsed
tetragonal phase results in a downward shift of the bands relative to the
uncollapsed phase and a reduction in the density of states at the Fermi
level.Yildrim2009 ; Goldman2009  Indeed, magnetotransport measurements in
rare-earth doped CaFe2As2 show a dramatic reduction in the magnitude of the
Hall coefficient upon cooling through the volume collapse transition.Saha2012
In addition, first-principles calculations for CaFe2As2 show that the strength
of the Fe-As bonds, the As-As mirror plane bonds, and the Fe spin-state, and
henceforth magnetism, are strongly coupled.Yildrim2009  Thus, the
disappearance of magnetic order with the onset of the collapsed tetragonal
phase may be an unsurprising consequence of the electronic structure mandated
by the collapsed tetragonal phase.

The onset of As-As interlayer bonding, as indicated by the contraction of
$d_{As-As}$, is not unique to (Ca0.67Sr0.33)Fe2As2. In fact, the other pure
alkaline earth 122 ferropnictide superconductors as well as some of their
rare-earth-doped counterpartsSaha2012  display identical behavior in $d_{As-
As}$, albeit at different pressures. Figure 7 shows $d_{As-As}$ as a function
of pressure for members of the (AE)Fe2As2 family. The horizontal line
represents the onset of As-As bonding at $d_{As-As}$=3.0 Å. Each compound has
been shown to undergo the volume collapse transition, and, accordingly, each
displays an abrupt contraction of the As-As separation. As the size of the
alkaline earth element increases, the unit cell volume of the crystal
structure increases, and $d_{As-As}$ increases. A natural consequence of this
unit cell volume expansion is that a larger pressure is required to achieve
sufficient lattice compression to invoke As-As bonding across the mirror plane
of the unit cell, and the volume collapse transition concordantly shifts to
higher pressures with increasing atomic radius of the alkaline earth element
(inset, Fig. 7).

Figure 7: (color online) Comparison of the pressure dependence of $d_{As-As}$
for the alkaline earth (AE)Fe2As2 compounds. The horizontal dashed line
indicates $d_{As-As}=$3.0 Å. Data for CaFe2As2, SrFe2As2, and BaFe2As2 are
from references [Goldman2009, ], [Uhoya2011, ], and [Mittal2011, ],
respectively. Inset: volume collapse transition pressure, $P_{ct}$, as a
function of unit cell volume at $P$=0. Lines through data points are guides to
the eye.

###  IV.2 As-Fe-As bond angles

Early studies of the iron-bearing oxypnictide superconductors noted an
empirical correlation between the As-Fe-As bond angles and the maximum
observed $T_{c}$.Ishida2009 ; Lee2008  Within the corrugated FeAs layers of
the crystal structure, there are two As-Fe-As bond angles: the two-fold or
intralayer angle, denoted as $\alpha$; and the four-fold or interlayer angle,
denoted as $\beta$ (see the inset of Fig. 8).Johnston2010  Due to the crystal
structure, these two bond angles move oppositely (as $\alpha$ increases,
$\beta$ decreases), but if the As atoms surrounding the Fe atoms are perfectly
tetrahedrally coordinated, then $\alpha$=$\beta$=109.47∘. In reference
[Lee2008, ], it was found that as the As-Fe-As bond angles of LnFeAsO1-y
approached this ideal tetrahedral angle, $T_{c}$ reached a maximum value near
55 K. Since then, this general empirical relationship has been noted in nearly
all families of ferropnictide superconductors.Paglione2010 ; Johnston2010
However, like many “rules” in condensed-matter physics, there are exceptions,
notably CsFe2As2 with a low $T_{c}$=2.6 K and As-Fe-As bond angles of 109.58∘
and 109.38∘.Gooch2010

Figure 8: (color online) (a) The room-temperature values of the two-fold,
$\alpha$, and four-fold, $\beta$, As-Fe-As bond angles as a function of
pressure. The ideal tetrahedral angle (109.47∘) is marked by the horizontal,
green, dashed line. The inset defines $\alpha$ and $\beta$ within the
corrugated FeAs component of the crystal structure. (b) The pressure
dependence of $T_{c}$ on the same pressure axis; symbols identical to Fig. 4.
Lines are guides to the eye.

Figure 8a displays the pressure dependence of the $\alpha$ and $\beta$ As-Fe-
As bond angles at room temperature. At ambient pressure, the corrugated FeAs
layers exhibit coordination close to the ideal tetrahedral configuration
($\alpha$=109.93∘, $\beta$=109.24∘). Applied pressure drives the structure
away from this ideal tetrahedral coordination, and, above the volume collapse
transition, the bond angles settle into relatively pressure-independent values
significantly disparate from the ideal tetrahedral condition.

The evolution of $T_{c}$ with pressure is reproduced in Fig. 8b for comparison
with the bond-angle evolution. While $T_{c}$ seems to be correlated with the
bond angles, with $T_{c}$ decreasing as the bond angles deviate from that of
the ideal tetrahedron, it should be emphasized that a one-to-one
correspondence is likely too simple of an explanation. The bond angle data in
Fig. 8a was determined at room temperature, while the determination of $T_{c}$
is clearly at low temperature. From the phase diagram in Fig. 4, it can be
seen that the superconducting phase occurs within the collapsed tetragonal
phase. If the relatively constant bond angles seen in the collapsed tetragonal
phase at room temperature are representative of that phase even at low
temperatures, then one might expect that the maximum in $T_{c}$ would occur
with bond angles $\alpha{\approx}$116∘ and $\beta{\approx}$106∘, distinctly
deviating from the ideal tetrahedral angle. Unfortunately, no low-temperature
structural data were acquired in this study, but low-temperature structural
characterization of pure and rare-earth doped CaFe2As2, which still exhibit
superconductivity, indicate that the bond angles tend away from the ideal
tetrahedral angle upon cooling.Kreyssig2008 ; Saha2012  Furthermore, the phase
diagram of Fig. 4 reveals that the superconducting state occurs in proximity
to the destruction of magnetic order, its associated structural transition,
and the occurrence of an isostructural volume collapse, certainly suggesting
that the appearance of superconductivity may be correlated with factors other
than structural parameters at room temperature.

###  IV.3 Structural and electronic phase diagrams

The structural and electronic phase diagram of (Ca0.67Sr0.33)Fe2As2
interpolates very well with the phase diagrams of its end member compounds as
well as the related BaFe2As2 compound. These phase diagram are shown together
in Fig. 9, highlighting the qualitative similarities within the AEFe2As2
system. Each compound obeys some general behavioral rules. At ambient-
pressure, each compound undergoes a structural/magnetic transition ($T_{N}$)
at sub-ambient temperatures. $T_{N}$ is suppressed with applied pressure and
abruptly disappears well above $T=0$. Superconductivity develops around this
discontinuous destruction of magnetism and persists as a several-GPa-wide dome
or half-dome in P-T space with a maximum $T_{c}$ occurring close to the
pressure at which $T_{N}$ abruptly vanishes. Each compound exhibits a
pressure-induced volume collapse at room temperature, and the volume collapse
transition $T_{ct}$ occurs with positive slope in P-T space (i.e.,
$dT_{c}/dP>0$) where it has been measured. The pressure axis of Fig. 9 does
not extend far enough to include the volume collapse transition in BaFe2As2,
which, from the data shown in Fig. 7, occurs near 26 GPa.Mittal2011

In CaFe2As2 and (Ca0.67Sr0.33)Fe2As2, measurements of the $T_{ct}(P)$ line
strongly suggest that the volume collapse itself is likely responsible for the
abrupt destruction of magnetic order. This further implies that magnetism is
limited to the orthorhombic phase, and that the collapsed tetragonal phase
does not support magnetic order. While the destruction of magnetism may be
linked to the onset of the collapsed tetragonal phase, whether that
destruction is driven by a reduction in the Fe moments, an altering of some
exchange coupling, or a more subtle change in the electronic structure is an
open question likely requiring both theoretical and experimental input to
reach a conclusion.

Figure 9: (color online) Comparison of the temperature-pressure phase diagrams
for (AE)Fe2As2. Closed circles represent $T_{N}$, open diamonds represent
$T_{c}$, and crossed or diagonal boxes represent $T_{ct}$. The volume collapse
transition as determined by electrical transport and x-ray diffraction are
denoted by ET and XRD, respectively. Electrical transport data for CaFe2As2,
SrFe2As2, and BaFe2As2 are from [Torikachvili2008a, ; Torikachvili2008b, ;
Colombier2009, ]; structural data are from [Goldman2009, ; Uhoya2011, ;
Mittal2011, ]. Lines and shaded regions are guides to the eye.

Unlike the volume collapse transition, the facts about the occurrence of
superconductivity in the AEFe2As2 systems are less clear. With the use of more
hydrostatic pressure conditions, researchers have found that the
superconducting dome, as defined by complete resistive transitions or
susceptibility data, in the parent compounds is generally excluded from the
magnetic, orthorhombic phase.Alireza2009 ; Colombier2009  The results on
(Ca0.67Sr0.33)Fe2As2 are consistent with this finding even with the less
hydrostatic steatite pressure-transmitting medium used in this study. However,
it is imperative to note that a study using He as the pressure-transmitting
medium revealed the absence of superconductivity in both the orthorhombic
(magnetic) and the collapsed tetragonal phases of CaFe2As2.Yu2009  A simple,
quantitative shift in the pressure at which superconductivity appears as a
function of hydrostaticity would not be entirely surprising, but the
qualitative difference (i.e., a complete lack of superconductivity) as a
function of pressure media creates a conundrum regarding the appearance of
superconductivity in the AEFe2As2 systems.

While the generic phase diagrams of chemical substitution and doping look
strikingly similar, there are subtle differences that question the roles of
structural and magnetic instabilities and suggest that a one-to-one
correspondence between pressure and substitution is too simple. Studies of
electron- and hole-doped BaFe2As2 have revealed a splitting of the nominally
coupled paramagnetic-antiferromagnetic and tetragonal-orthorhombic phase
transitions with increasing dopant content.Ni2008 ; Chu2009 ; Pratt2009 ;
Urbano2010  Superconductivity is seen to develop within the antiferromagnetic,
orthorhombic phase with both Co and K substitution in BaFe2As2,Ni2008  but the
optimal $T_{c}$ is achieved within the paramagnetic, uncollapsed tetragonal
phase. Though superconductivity develops within the antiferromagnetic state,
the suppression of that state, at least to some degree, seems to be a
necessary ingredient for superconductivity. Unlike doping of the BaFe2As2 end
member, Rh-doping into CaFe2As2 does not reveal a significant splitting of the
paramagnetic-antiferromagnetic and tetragonal-orthorhombic phase transitions,
but $T_{N}$ is nonetheless suppressed with increasing doping.Danura2011
Superconductivity is seen only in the paramagnetic, uncollapsed tetragonal
phase, and the onset of a doping-induced collapsed tetragonal phase destroys
superconductivity. Recently, a doping study using rare earth (RE) elements in
Ca1-xRExFe2As2 has shown that the components of electron doping and chemical
pressure can be effectively separated, and that each of these components plays
a different role in manifesting superconductivity.Saha2012  Within the
Ca1-xRExFe2As2 series, superconductivity occurs in either the uncollapsed or
collapsed tetragonal phases, but not within the antiferromagnetic,
orthorhombic phase.

Assimilating the doping- and pressure-dependent phase diagrams of the 122
systems in order to better understand the nature of the high-temperature
superconductivity seen therein is a challenging problem. Given the
overwhelming evidence that electron- or hole-doping plays an important role in
the development of superconductivity within the chemically substituted
AEFe2As2 systems and given the occurrence of strain-induced superconductivity
in SrFe2As2,Saha2009  it is tempting to posit that the small superconducting
dome seen under pressure may be a product of some effective doping induced by
non-hydrostatic pressure conditions, possibly resulting not only from the
pressure-transmitting medium but from the volume collapse transition itself.
Alternatively, without high-fidelity, low-temperature structural data under
pressure, it is difficult to exclude structural phase inhomogeneity (e.g.,
coexistence of the collapsed and uncollapsed tetragonal phases) as a possible
cause of the proximity of superconductivity to the observed structural
instabilities. More work demarcating the $T_{ct}(P)$ lines in SrFe2As2 and
BaFe2As2 under pressure may help to illuminate any possible connections
between pressure-induced superconductivity and the isostructural volume
collapse.

##  V Conclusions

The compound (Ca0.67Sr0.33)Fe2As2 under pressure behaves intermediate between
its two end member parent compounds. With applied pressure, the concomitant
structural/magnetic transition is suppressed with no evidence favoring a
splitting of the two nominally coupled transitions. The AFM, orthorhombic
phase is abruptly cut off by an isostructural volume collapse resulting in a
collapsed tetragonal phase. The volume collapse transition is driven by the
development of As-As bonding across the mirror plane of the crystal structure,
contracting the $c$-axis of the unit cell and likely affecting the Fe-As
bonding and potentially the magnetic state of the Fe atoms. The collapsed
tetragonal phase supports superconductivity with a maximum $T_{c}$ near 22 K.
There is no obvious structural parameter that defines the magnitude of
$T_{c}$, but the proximity of the superconducting phase to the suppression of
magnetism as well as the onset of the collapsed tetragonal phase suggests that
magnetic interactions and/or structural inhomogeneity may both play a role in
the development of pressure-induced superconductivity in these systems.

##  VI Acknowledgments

We are grateful to Z. Jenei and K. Visbeck for assistance with cell
preparations. JRJ and STW are supported by the Science Campaign at Lawrence
Livermore National Laboratory. Portions of this work were performed under LDRD
(11-LW-003). Lawrence Livermore National Laboratory is operated by Lawrence
Livermore National Security, LLC, for the U.S. Department of Energy, National
Nuclear Security Administration under Contract DE-AC52-07NA27344. Portions of
this work were performed at HPCAT (Sector 16), Advanced Photon Source (APS),
Argonne National Laboratory. HPCAT is supported by CIW, CDAC, UNLV and LLNL
through funding from DOE-NNSA, DOE-BES and NSF. Use of the Advanced Photon
Source, an Office of Science User Facility operated for the U.S. Department of
Energy (DOE) Office of Science by Argonne National Laboratory, was supported
by the U.S. DOE under Contract No. DE-AC02-06CH11357. Beamtime was provided
through the Carnegie-DOE Alliance Center (CDAC). This work was partially
supported by AFOSR-MURI Grant No. FA9550-09-1-0603. YKV acknowledges support
from DOE-NNSA Grant No. DE-FG52-10NA29660.

## References

  * (1) Y. Kamihara, T. Watanabe, M. Hirano, and H. Hosono, J. Am. Chem. Soc. 130, 3296 (2008).
  * (2) K. Ishida, Y. Nakai, and H. Hosono, J. Phys. Soc. Jpn. 78, 062001 (2009).
  * (3) Z. -A. Ren, W. Lu, J. Yang, W. Yi, X. -L. Shen, Z. -C. Li, G. -C. Che, X. -L. Dong, L. -L. Sun, F. Zhou, and Z. -X. Zhao, Chin. Phys. Lett. 25, 2215 (2008).
  * (4) J. Paglione and R. L. Greene, Nat. Phys. 6, 645 (2010).
  * (5) D. C. Johnston, Adv. Phys. 59, 803 (2010).
  * (6) P. A. Lee, N. Nagaosa, and X. -G. Wen, Rev. Mod. Phys. 78, 17 (2006).
  * (7) M. D. Lumsden and A. D. Christianson, J. Phys.: Condens. Matter 22, 203203 (2010).
  * (8) M. Gooch, B. Lv, K. Sasmal, J. H. Tapp, A. J. Tang, A. M. Guloy, B. Lorenz, and C. W. Chu, Physica C 470, S276 (2010).
  * (9) H. S. Jeevan, Z. Hossain, D. Kasinathan, H. Rosner, C. Geibel, and P. Gegenwart, Phys. Rev. B 78, 052502 (2008).
  * (10) M. Rotter, M. Tegel, D. Johrendt, I. Schellenberg, W. Hermes, and R. Pöttgen, Phys. Rev. B 78, 020503(R) (2008).
  * (11) Q. Huang, Y. Qiu, W. Bao, M. A. Green, J. W. Lynn, Y. C. Gasparovic, T. Wu, G. Wu, and X. H. Chen, Phys. Rev. Lett. 101, 257003 (2008).
  * (12) A. I. Goldman, D. N. Argyriou, B. Ouladdiaf, T. Chatterji, A. Kreyssig, S. Nandi, N. Ni, S. L. Bud’ko, P. C. Canfield, and R. J. McQueeney, Phys. Rev. B 78, 100506(R) (2008).
  * (13) K. Kaneko, A. Hoser, N. Caroca-Canales, A. Jesche, C. Krellner, O. Stockert, and C. Geibel, Phys Rev. B 78, 212502 (2008).
  * (14) Y. Su, P. Link, A. Schneidewind, Th. Wolf, P. Adelmann, Y. Xiao, M. Meven, R. Mittal, M. Rotter, D. Johrendt, Th. Brueckel, and M. Loewenhaupt, Phys. Rev. B 79, 064504 (2009).
  * (15) N. Ni, M. E. Tillman, J. -Q. Yan, A. Kracher, S. T. Hannahs, S. L. Bud’ko, and P. C. Canfield, Phys. Rev. B 78, 214515 (2008).
  * (16) J. -H. Chu, J. G. Analytis, C. Kucharczyk, and I. R. Fisher, Phys. Rev. B 79, 014506 (2009).
  * (17) D. K. Pratt, W. Tian, A. Kreyssig, J. L. Zaretsky, S. Nandi, N. Ni, S. L. Bud’ko, P. C. Canfield, A. I. Goldman, and R. J. McQueeney, Phys. Rev. Lett. 103, 087001 (2009).
  * (18) R. R. Urbano, E. L. Green, W. G. Moulton, A. P. Reyes, P. L. Kuhns, E. M. Bittar, C. Adriano, T. M. Garitezi, L. Bufaiçal, and P. G. Pagliuso, Phys. Rev. Lett. 105, 107001 (2010).
  * (19) F. Han, X. Zhu, P. Cheng, G. Mu, Y. Jia, L. Fang, Y. Wang, H. Luo, B. Zeng, B. Shen, L. Shan, C. Ren, H. -H. Wen, Phys. Rev. B 80, 024506 (2009).
  * (20) S. R. Saha, N. P. Butch, T. Drye, J. Magill, S. Ziemak , K. Kirshenbaum, P. Y. Zavalij, J. W. Lynn, and J. Paglione, Phys. Rev. B 85, 024525 (2012).
  * (21) W. Yu, A. A. Aczel, T. J. Williams, S. L. Bud’ko, N. Ni, P. C. Canfield, and G. M. Luke, Phys. Rev. B 79, 020511(R) (2009).
  * (22) A. I. Goldman, A. Kreyssig, K. Prokes̆, D. K. Pratt, D. N. Argyriou, J. W. Lynn, S. Nandi, S. A. J. Kimber, Y. Chen, Y. B. Lee, G. Samolyuk, J. B. Leão, S. J. Poulton, S. L. Bud’ko, N. Ni, P. C. Canfield, B. N. Harmon, and R. J. McQueeney, Phys. Rev. B 79, 0.24513 (2009).
  * (23) W. Uhoya, A.Stemshorn, G. Tsoi, Y. K. Vohra, A. S. Sefat, B. C. Sales, K. M. Hope, and S. T. Weir, Phys. Rev. B 82, 144118 (2010).
  * (24) W. Uhoya, G. Tsoi, Y. K. Vohra, M. A. McGuire, A. S. Sefat, B. C. Sales, D. Mandrus, and S. T. Weir, J. Phys.: Condens. Matter 22, 292202 (2010).
  * (25) W. O. Uhoya, J. M. Montgomery, G. M. Tsoi, Y. K. Vohra, M. A. McGuire, A. S. Sefat, B. C. Sales, and S. T. Weir, J. Phys.: Condens. Matt. 23, 122201 (2011).
  * (26) R. Mittal, S. K. Mishra, S. L. Chaplot, S. V. Ovsyannikov, E. Greenberg, D. M. Trots, L. Dubrovinsky, Y. Su, Th. Brueckel, S. Matsuishi, H. Hosono, G. Garabino, Phys. Rev. B 83, 054503 (2011).
  * (27) M. Danura, K. Kudo, Y. Oshiro, S. Araki, T. C. Kobayashi, M. Nohara, J. Phys. Soc. Jpn. 80, 103701 (2011).
  * (28) T. Yildrim, Phys. Rev. Lett. 102, 037003 (2009).
  * (29) I. I. Mazin and J. Schmalian, Physica C 469, 614 (2009).
  * (30) K. Kirshenbaum, N. P. Butch, S. R. Saha, P. Y. Zavalij, B. G. Ueland, J. W. Lynn, and J. Paglione, arXiv (2012).
  * (31) S. R. Saha, N. P. Butch, K. Kirshenbaum, J. Paglione, and P. Zavalij, Phys. Rev. Lett. 103, 037005 (2009).
  * (32) S. T. Weir, J. Akella, C. Aracne-Ruddle, Y. Vohra, and S. A. Catledge, Appl. Phys. Lett. 77, 3400 (2000).
  * (33) J. R. Patterson, S. A. Catledge, Y. K. Vohra, J. Akella, and S. T. Weir, Phys. Rev. Lett. 85, 5364 (2000).
  * (34) D. D. Jackson, J. R. Jeffries, W. Qiu, J. D. Griffith, S. McCall, C. Aracne, M. Fluss, M. B. Maple, S. T. Weir, and Y. K. Vohra, Phys. Rev. B 74, 174401 (2006).
  * (35) H. K. Mao, J. Xu, and P. M. Bell, J. Geophys. Res. 91, 4673 (1986).
  * (36) W. L. Vos and J. A. Schouten, J. Appl. Phys. 69, 6744 (1991).
  * (37) A. Hammersley, S. Svensson, M. Hanfland, A. Fitch, and D. Häusermann, High Press. Res. 14, 235 (1996).
  * (38) A. C. Larson and R. B. Von Dreele, Los Alamos National Laboratory Report LAUR 86-748 (1994).
  * (39) B. H. Toby, J. Appl. Crystallogr. 34, 210 (2001).
  * (40) S. R. Saha, K. Kirshenbaum, N. P. Butch, J. Paglione, and P. Zavalij, J. Phys. Conf. Ser. 273, 012104 (2011).
  * (41) M. Torikachvili, S. L. Bud’ko, N. Ni, and P. C. Canfield, Phys. Rev. Lett. 101, 057006 (2008).
  * (42) M. Torikachvili, S. L. Bud’ko, N. Ni, and P. C. Canfield, Phys. Rev. B 78, 104527 (2008).
  * (43) P. C. Canfield, S. L. Bud’ko, N. Ni, A. Kreyssig, A. I. Golman, R. J. McQueeney, M. S. Torikachvili, D. N. Argyriou, G. Luke, and W. Yu, Physica C 469, 404 (2009).
  * (44) R. Hoffman and C. Zheng, J. Phys. Chem. 89, 4175 (1985).
  * (45) G. Just and P. Paufler, J. Alloy Compd. 232, 1 (1996).
  * (46) P. Pyykkö and M. Atsumi, Chem. Eur. J. 15, 186 (2009).
  * (47) C.-H. Lee, A. Iyo, H. Eisaki, H. Kito, M. T. Fernandez-Diaz, T. Ito, K. Kihou, H. Matsuhata, M. Braden, and K. Yamada, J. Phys. Soc. Jpn. 77, 083704 (2008).
  * (48) A. Kreyssig, M. A. Green, Y. Lee, G. D. Samolyuk, P. Zajdel, J. W. Lynn, S. L. Bud’ko, M. S. Torikachvili, N. Ni, S. Nandi, J. B. Leao, S. J. Poulton, D. N. Argyriou, B. N. Harmon, R. J. McQueeney, P. C. Canfield, and A. I. Goldman, Phys. Rev. B 78, 184517 (2008).
  * (49) P. L. Alireza, Y. T. C. Ko, J. Gillett, C. M. Petrone, J. M. Cole, G. G. Lonzarich, and S. E. Sebastian, J. Phys.: Condens. Matter 21, 012208 (2009).
  * (50) E. Colombier, S. L. Bud’ko, N. Ni, and P. C. Canfield, Phys. Rev. B 79, 224518 (2009).


# CORTEX: Large-Scale Brain Simulator Utilizing Indegree Sub-Graph
Decomposition on Fugaku Supercomputer

Tianxiang Lyu  Graduate School of Medicine
Juntendo University, Japan
Email<EMAIL_ADDRESS>Mitsuhisa Sato  Faculty of Health Data Science
Graduate School of Medicine
Juntendo University
RIKEN Center of Computational Science, Japan
Email<EMAIL_ADDRESS>Shigeki Aoki  Faculty of Health Data Science
Graduate School of Medicine
Juntendo University, Japan
Email<EMAIL_ADDRESS>Ryutaro Himeno  Faculty of Health Data Science
Graduate School of Medicine
Juntendo University, Japan
Email<EMAIL_ADDRESS>Zhe Sun  Faculty of Health Data Science
Graduate School of Medicine
Juntendo University
RIKEN Center for Advanced Photonics, Japan
Email<EMAIL_ADDRESS>
###### Abstract

We introduce CORTEX, an algorithmic framework designed for large-scale brain
simulation. Leveraging the computational capacity of the Fugaku Supercomputer,
CORTEX maximizes available problem size and processing performance. Our
primary innovation, Indegree Sub-Graph Decomposition, along with a suite of
parallel algorithms, facilitates efficient domain decomposition by segmenting
the global graph structure into smaller, identically structured sub-graphs.
This segmentation allows for parallel processing of synaptic interactions
without inter-process dependencies, effectively eliminating data racing at the
thread level without necessitating mutexes or atomic operations. Additionally,
this strategy enhances the overlap of communication and computation. Benchmark
tests conducted on spiking neural networks, characterized by biological
parameters, have demonstrated significant enhancements in both problem size
and simulation performance, surpassing the capabilities of the current leading
open-source solution, the NEST Simulator. Our work offers a powerful new tool
for the field of neuromorphic computing and understanding brain function.

###### Index Terms:

Fugaku Supercomputer, Large Scale Brain Simulation, Spiking Neural Networks,
Indegree Sub-Graph Decomposition, HPC Application

##  I BACKGROUND AND MOTIVATION

Recent years have seen an unprecedented accumulation of data at various levels
of brain organisation, thanks in large part to extensive biological research
exemplified by works such as [1], with a detailed mapping of the brain’s
cytoarchitecture. There will be a sure demand for integrating all of them into
a unified model, providing a unified view on how brain works as a whole[3].
However, the fundamental barriers are scale and complexity[4], making it hard
to achieve a general comprehension by experimental or theoretical methods[2,
5], dealing with 86 billion neurons in the human brain and the interactions
between molecules, neurons, microcircuits, and brain areas. To transcend these
barriers, brain simulation was developed to provide a new way to gain an
overview of all the forests of datasets, with reorganizing and integrating
them in the context of the whole brain. Optimistically, it seems possible to
not only fill the gaps that experiments and theoretical analysis will never
manage to fill, but also find ways through the forests by considering the
ecosystem of the brain architecture[5].
Enabled by the computational power of leading-edge supercomputers, which has
been steadily and exponentially increasing for the past 20 years[7], brain
simulation towards human scale is one of the most ambitious scientific
challenges in the 21st century[9], playing an indispensable role in the
investigation of the multi-scale brain[5]. The goal of brain simulation is to
achieve a dense digital reconstruction of brain dynamics from experimental
datasets and the most fundamental principles we can isolate to understand,
linking the multiple layers from cells, circuits, areas to brain function and
behavior, in order to make progress systematic and understanding tractable[5].
However, due to the sparsity of brain architecture in both spatial and
temporal domains, many existing algorithms and libraries[18] can’t be directly
applied in brain simulation for optimizations. What’s worse, in some
simulators, mutexes or atomic operations[12, 13] are introduced, in order to
avoid data racing at the thread level, but resulting in poor performance. At
last, the communication overhead might be unbearable in large-scale
simulation, because the ratio of computation and communication in brain
simulation is much lower than primary HPC applications.
Fortunately, most instances of interest are far from the worst-case scenario
and the performance of simulation can have orders of magnitude improvements
compared with the naive approach, by focusing on the unique properties of
brain architecture. Based on directed graph[17], the main technical
contribution of our paper is Indegree Sub-Graph Decomposition, along with a
suite of parallel algorithms, which can not only perform domain decomposition
of entire brain architecture among all processes with maximized problem size,
but also avoid data racing at the thread level without any mutex or atomic
operation. Having achieved impressive improvements in various aspects, CORTEX
represents as a powerful brain simulation framework on leading-edge
supercomputers[14].

###  I-A Brain Architecture and Spiking Neural Networks

In this section, a brief introduction about brain architecture consists of
neurons and synapses fundamentally. More detailed reviews on neuron models
have been published elsewhere[11]. A neuron is a basic unit in the brain,
connected with others through synapses as Fig. 1 shows. A ”spike”, also known
as a neuronal spike or action potential, is the fundamental unit of
information transfer in the nervous system. When a neuron is sufficiently
stimulated, it generates a rapid, transient change in voltage. This phenomenon
is referred to as a spike. All these neurons and synapses are forming a whole,
which can be regarded as spiking neural networks (SNNs).
In brain simulation, the states inside neurons will be updated at each time
step, according to the integral of its dynamics equations from modeling
methods. In this paper, the leaky-integrated and fire (LIF) model[11] is used
as the modeling method of neurons, with ordinary differential equations shown
at the following:

$\tau\frac{du_{i}(t)}{dt}=-u_{i}(t)+u_{rest}+R(I_{syn,i}(t)+I_{ext,i}(t))$ (1)
$\rm{if}\,\mathit{u(t)}=\theta\Longrightarrow\mathit{u\rightarrow
u_{r}\quad\rm{within}\quad\mathit{(t,t+t_{RP}]}}$ (2)

where $u_{i}$ is membrane potential, $u_{rest}$ is the resting potential,
$\tau$ is membrane time constant, $R$ is membrane resistance, $I_{syn,i}$ is
synaptic current, $I_{ext,i}$ is external current, $t$ is time, $\theta$ is
spike threshold, $u_{r}$ is reset value of membrane potential, and $t_{RP}$ is
refractory period. As shown in Fig. 1, membrane potential will rapidly rise
and then fall into $u_{r}$ during $(t,t+t_{RP}]$, which is called a spike. The
synaptic current[11] can be described as:

$I_{syn,i}(t)=\sum_{j}\sum_{f}\delta(t-t^{f}_{j})W_{ji}(t)g_{syn}(u_{i}(t)-E_{syn})$
(3)

where $W_{ji}$ is a synaptic weight from neuron $j$ to $i$, $g_{syn}$ is a
time-dependent function of synaptic conductance, $t^{f}_{j}$ is the arrival
time of spike $f$ from neuron $j$, and $E_{syn}$ is reversal potential of a
synaptic ion channel. An exponential function was used as $g_{syn}$.
Different from other simulators, most ideas of optimizations in CORTEX are not
from practical coding or test, but from graph abstraction of SNNs with
theoretical traceable analysis, guiding the parallelization from the process
level to the thread level. In the section of OVERVIEW OF CORTEX, graph
abstraction of SNNs will be described in detail.

Figure 1: Schematic Representation of A Neuron and Its Spike: Illustrating the
basic elements of (SNNs), where neurons are connected through synapses,
forming a complex network. The ”spike” serves as the primary unit of
information processing in SNNs.

###  I-B Sparsity in Spatial and Temporal Domains

In brain simulation, the computing intensity is much lower than that in
primary HPC applications, because of its sparsity in both spatial and temporal
domains, which presents the main challenges. First, it’s easy to realize that
the synaptic interactions inside a brain can be really sparse, with only 160
trillion synapses for 86 billion neurons in our human brain. In other words,
each neuron only connects with thousands of others, which are much fewer than
the neurons themselves. Therefore, the computation among these sparse
interactions must be hard to be accelerated by mainstream computing
architectures designed with increasingly higher ratios of FLOPS$/$Byte.
Second, another type of sparsity is shown in temporal domain, related to
varied synaptic delays. Different from other graph computations, synaptic
interactions from pre-synaptic neurons to post-synaptic neurons will not take
effect immediately but after a few time steps, where the delays of different
synaptic interactions can be not exactly the same. In general, these varied
synaptic delays must be the main reason for low computing efficiency.

###  I-C Ratio of Computation and Communication

As a typical sparse problem, the ratio of computation and communication has a
significant impact on general performance. Also, we are keen to discuss the
reason behind the chosen LIF neuron model[21, 22, 23]. In the section on
VERIFICATION AND EVALUATION, the selected modeling method has been generally
applied in brain science investigation[27, 28, 29, 30], because of its
lightweight computing intensity, but precise approximation in neural dynamics,
which is only slightly different from Hodgkin–Huxley neuron model[31]
incorporating more biological details with much higher computing intensity.
Simulations[32] using neuron models with very high computing intensity show
absolutely better results in scalability, which can be regarded as ”good”
cases for massive parallel computing. In our consideration, these ”good” cases
exhibits only the upper bound performance of the simulator, which is too
trivial to demonstrate the contribution in optimizations. Therefore, only
”bad” cases, with low ratio of computation to communication, are driving the
design and implementation of brain simulation towards higher actual
performance with maximal problem size, facing the sparsity of brain
architecture in reality.

###  I-D Current State of The Art

Various types of large-scale brain simulation have been proposed with levels
of description and scales according to their purposes. A large-scale
simulation of the rodent primary somatosensory cortex, with detailed
morphological features of neurons, has been reported by the Blue Brain
Project[10]. In this case, the study focused on conditions of synchronization
and influences of morphological features on network activities, simulated on
the supercomputer Piz Daint at EPFL. Large-scale simulations of the cerebral
cortex have been performed by various research teams. On Blue Gene P[33], a
cat-scale corticothalamic circuit with 1.6 billion neurons was presented by
IBM, using the C2 simulator. Subsequently, using the NEST Simulator [34, 35],
Kunkel and her colleagues performed non-spatial cortical models with 1.8
billion neurons and ten trillion of synapses on the K computer. Forked from
NEST Simulator, MONET[37] has successfully simulated a layered sheet-type of a
model cerebral cortex with five billion neurons on the K computer by Igarashi
and his colleagues, also a human-scale cerebellar model with 68 billion
neurons[6]. However, the cases we mentioned above shouldn’t be a direct
benchmark or baseline, because there might be some tricks in modeling methods
or accuracy compression on floating point numbers with single or half
precision. Up till now, human-scale whole brain simulation with cerebral
cortex, cerebellum, and thalamus has never been conducted to date. Above all,
NEST Simulator is the SOTA[15] available with open-source code. Therefore, the
comparison will be shown between CORTEX and NEST Simulator, in the section of
VERIFICATION AND EVALUATION.

###  I-E System and Environment of Fugaku

The supercomputer Fugaku[14] is a homogeneous, CPU-based distributed machine
with 158,976 nodes. Based on the ARM v8.2A architecture, the CPU of the
system, A64FX, has 32 GB HBM2 memory and 48 compute cores, with total memory
bandwidth 1024 GB/s. 12 compute cores form one Core Memory Group (CMG). On the
ALUs side, there is a couple of 512-bit wide SIMD units inside each compute
core, where the whole CPU performs 3.072 TF/s for double-precision, 6.144 TF/s
for single-precision, and 12.288 TF/s for half-precision. As to cache, Each
compute core contains 64 KB of L1 data cache, and 32 MB L2 cache is shared by
all. For interconnect, low latency and high throughput are achieved by a
6-dimensional mesh/torus network, with a 6.8 GB/s link bandwidth and 40.8 GB/s
of injection bandwidth per compute node.
Using a hybrid parallelized strategy, an MPI process is assigned to each CMG,
with an OpenMP thread allocated to each compute core, to achieve maximal
bandwidth and minimal latency in this NUMA architecture. On the compiler side,
the Fujitsu C/C++ compiler 4.10 is used in Trad mode with automatic
vectorization and software pipelining.

##  II OVERVIEW OF CORTEX

In this paper, CORTEX, our highly optimized framework for brain simulation on
Fugaku supercomputer is proposed, focusing on biological data driven brain
architecture, based on the graph abstraction. In this section, the concepts
and principles will be described first, showing the basic ideas from which we
start.

###  II-A Graph Abstraction of SNNs

(a) Graph Abstraction

(b) Basic Elements

Figure 2: Graph Abstraction of SNNs

In a graph to represent the brain architecture, each vertex can be a neuron,
and each edge is a synaptic interaction from, or input to, a neuron (vertex).
The neural dynamics are computed on neurons, where synaptic interactions act
on target neurons along edges. These basic abstractions of SNNs are shown in
Fig. 2, and we are going to describe how they form a brain architecture in
simulation.

####  II-A1 Indegree and Outdegree Sub-graph

A directed graph is an ordered pair $G=(V,E)$, comprising:

  * •

$V$, a set of vertices, which can be regarded as neurons;

  * •

$E\subseteq\\{(x,y)|(x,y)\in V^{2}$ and $x\neq y\\}$, a set of edges, as
synaptic interactions, which are ordered pairs of vertices, where the
condition: $x\neq y$ can be ignored in some SNNs.

Also, there can be a triplet:

${}^{*}\\!S=(^{*}\\!V^{pre},^{*}\\!V^{post},^{*}\\!E)$ (4)

where $*\in\\{in,out\\}$, defined as an indegree or outdegree format of graph
$G$, comprising:

  * •

${}^{*}\\!V^{pre}\subseteq V$, the set of pre-vertices;

  * •

${}^{*}\\!V^{post}\subseteq V$, the set of post-vertices.

  * •

${}^{*}\\!E\subseteq\\{(x,y)|x\in V^{pre}$ and $y\in V^{post}\\}$, the set of
edges, with a bijection $f:E\rightarrow^{*}\\!E$

However, if the scale of simulation becomes large enough, the entire graph
should be stored in distributed memory with domain decomposition. Assumed that
$\tilde{V}$ is a subset of $V$, the indegree sub-graph based on $\tilde{V}$
can be well defined as:

${}^{in}\\!S(\tilde{V})=(^{in}\\!\tilde{V}^{pre},\tilde{V},^{in}\\!\tilde{E})$
(5)

where ${}^{in}\\!\tilde{V}^{pre}=\\{x|(x,y)\in^{in}\\!\tilde{E}\\}$,
${}^{in}\\!\tilde{E}=\\{(x,y)|(x,y)\in^{in}\\!E$ and $y\in\tilde{V}\\}$.
Similarly, the outdegree sub-graph should be:

${}^{out}\\!S(\tilde{V})=(\tilde{V},^{out}\\!\tilde{V}^{post},^{out}\\!\tilde{E})$
(6)

with ${}^{out}\\!\tilde{V}^{post}=\\{y|(x,y)\in^{out}\\!\tilde{E}\\}$,
${}^{out}\\!\tilde{E}=\\{(x,y)|(x,y)\in^{out}\\!E$ and $x\in\tilde{V}\\}$.
When ${}^{*}\\!S$ is defined without specification of $\tilde{V}$, definition
(4) should be adopted. When ${}^{*}\\!S(\tilde{V})$ has been specified with
$\tilde{V}$, as a function from sets of vertices to sub-graphs, we use the
definition in (5) or (6), without misunderstanding. Then, it is
straightforward to define binary operation $\circledast$ on the following:

${}^{*}\\!S_{a}\circledast^{*}\\!S_{b}=(^{*}\\!V^{pre}_{a}\odot^{*}\\!V^{pre}_{b},^{*}\\!V^{post}_{a}\odot^{*}\\!V^{post}_{b},^{*}\\!E_{a}\odot^{*}\\!E_{b})$
(7)

where $(\circledast,\odot)\in\\{(\barwedge,\cap),(\veebar,\cup)\\}$.
Obviously, each of them is commutative and associated, further, there is a
homomorphism:

${}^{*}\\!S(V_{a})\circledast^{*}\\!S(V_{b})=^{*}\\!S(V_{a}\odot V_{b})$ (8)

An illustration can be shown as Fig. 3.

(a) $G=(V,E)$

(b) $\tilde{V}\subseteq V$

(c) Pre-Vertices and Edges of the Indegree Sub-graph

(d) Post-Vertices and Edges of the Indegree Sub-graph

(e) Pre-Vertices and Edges of the Outdegree Sub-graph

(f) Post-Vertices and Edges of the Outdegree Sub-graph

Figure 3: A Sub-graph of Indegree and Outdegree Format

####  II-A2 Synaptic Interactions on Sub-graphs

Considering synaptic interactions, the question comes to which kind of sub-
graphs we should choose for parallel implementation. First, there can be a
well partition of $V$, into $n$ parts:

$\\{V_{1},V_{2},V_{3},...,V_{n}\\}$ (9)

where $V_{i}\cap V_{j}=\varnothing$ with $i\neq j$, and
$\bigcup^{n}_{1}V_{i}=V$. Also, there are $n$ sub-graphs as:

$\\{^{*}\\!S(V_{1}),^{*}\\!S(V_{2}),^{*}\\!S(V_{3}),...,^{*}\\!S(V_{n})\\}$
(10)

where each one is corresponded to one process.
In one time step, there are very few pre-synaptic neurons in spiking state,
with their corresponding synaptic interactions and post-synaptic neurons. So,
there can be a set of all spiking pre-synaptic neurons ${}^{*}\\!V^{pre}_{s}$,
and the corresponding post-synaptic neurons ${}^{*}\\!V^{post}_{s}$ with
synaptic interactions ${}^{*}\\!E_{s}$, forming a spiking graph
${}^{*}\\!S_{s}=(^{*}\\!V^{pre}_{s},^{*}\\!V^{post}_{s},^{*}\\!E_{s})$. For a
sub-graph ${}^{*}\\!S(V_{i})$ , its spiking sub-graph can be defined as:

${}^{*}\\!S_{s}(V_{i})=^{*}\\!S(V_{i})\barwedge^{*}\\!S_{s}$ (11)

Then, as we mentioned before, the synaptic interactions will effect at post-
synaptic neurons, where the data instances of them should be writable. In
parallel implementation, we want that all write operations can be performed
without any dependency between spiking sub-graphs. Based on what we mentioned
before, dependencies between 2 arbitrary sub-graphs can be expressed as:

${}^{*}\\!S_{a}\barwedge^{*}\\!S_{b}$ (12)

Because of the homomorphism in (8), with the commutative and associated
operation, we have:

${}^{*}\\!S_{s}(V_{i})\barwedge^{*}\\!S_{s}(V_{j})=^{*}\\!S(V_{i}\cap
V_{j})\barwedge^{*}\\!S_{s}$ (13)

Obviously, for indegree sub-graphs, we have:

${}^{in}\\!S(V_{i}\cap
V_{j})=(^{in}\\!V^{pre}_{i}\cap^{in}\\!V^{pre}_{j},\varnothing,\varnothing)$
(14)

while in outdegree sub-graphs:

${}^{out}\\!S(V_{i}\cap
V_{j})=(\varnothing,^{out}\\!V^{post}_{i}\cap^{out}\\!V^{post}_{j},\varnothing)$
(15)

What we mentioned in (14) or (15) can indicate which kind of data should be
synchronized for each modification, as shown in Fig. 4 and 5. Therefore,
indegree sub-graphs should be the only choice.

(a) Indegree Spiking Graph

(b) Indegree Sub-graph Tangerine

(c) Indegree Sub-graph Grape

(d) Indegree Sub-graph Ice

Figure 4: Synaptic Interactions on the Indegree Sub-graph: In Fig.4(a),
spiking neurons with active synaptic interactions have been highlight with red
border (No. 1 and No.6). At each time step, once the spiking graph has been
well defined, synaptic interactions can take effect on corresponding post-
synaptic neurons without dependency between different indegree sub-graphs.

(a) Outdegree Spiking Graph

(b) Outdegree Sub-graph Tangerine

(c) Outdegree Sub-graph Grape

(d) Outdegree Sub-graph Ice

Figure 5: Synaptic Interactions on the Outdegree Sub-graph: In Fig. 5(a),
spiking neurons with active synaptic interactions have been highlight with red
border (No. 1 and No.6). At this time step, there are 2 synaptic interactions
taking effect on neuron 9 in 2 outdegree sub-graphs respectively. Between each
synaptic interaction with its nonlinear dynamics, the states of all post-
synaptic neurons with the same ID should be synchronized among all sub-graphs.

###  II-B Simulation Procedure

(a) Connectome Data

(b) Problem Description

(c) Domain Decomposition

(d) Simulation

Figure 6: Brain Simulation Procedure: fig (a) is an illustration of brain
connectome data extracted from Diffusion magnetic resonance imaging. In fig
(d), the neurons in color red are the spiking neurons, whose synaptic
interactions with color red too will effect at the post-synaptic neurons after
corresponding synaptic delay.

As to run a brain simulation by CORTEX on Fugaku supercomputer, it might
contain several steps. Fig. 6 describes 4 main steps of brain simulation.
Different from other primary HPC applications like Computational Fluid
Dynamics (CFD) or Molecular Dynamics (MD), only a few parts of synaptic
interactions (edges) will be active in one time step. As shown in Fig. 6(d),
the synaptic interactions from spiking pre-synaptic neurons in red are taking
effect on their post-synaptic neurons. After that, all neurons (vertices) will
update their variables like membrane potential by neural dynamics, generating
new spiking pre-synaptic neurons, to start a new loop.

##  III MAJOR INNOVATIVE CONTRIBUTIONS

Based on our graph abstraction of SNNs, the main idea for optimized
parallelization has been clear. Then, it comes to implementation, facing the
sparsity in both spatial and temporal domains. Efforts on storage, computing
and networking are adopted for not only application-available problem size but
also general performance of simulation. In brief, our key innovations are
summarized as follows:

  * •

A customized domain decomposition method comprising of 2 steps: Area-Processes
Mapping and Multisection Division with Sampling Method, is introduced to
achieve maximal problem size with high performance.

  * •

A elaborate multi-threading parallel scheme without mutex or atomic operation
when computing synaptic interactions on edges, which is the hotspot of whole
simulation.

  * •

An optimized spikes broadcast method with dedicated thread for communication,
in order to overlap communication and computation as much as possible.

###  III-A Domain Decomposition

Fortunately, with the homomorphism in (8), the domain decomposition of graph
into sub-graphs can be transferred to the partition on vertices. Then, the
problem arises to generate appropriate partitions of vertices with their
corresponding indegree sub-graphs to achieve best performance in simulation,
which is a typical NP-hardness[19]. Therefore, the solution should be
generally derived using heuristics and approximation algorithms.

####  III-A1 Varied Density of Synaptic Interactions

Before domain decomposition, the properties of biological data driven SNNs
should be taken into consideration, owing to the varied density of both
synaptic interactions and neurons.

Figure 7: An example of marmoset brain atlas from Brain/MINDS Project[20],
where each color patch represents a distinct brain area.

As shown in Fig. 7, with the Atlas, a brain architecture can be separated into
areas. Each area is related to different macro behavior or brain function. For
example, the primary vision cortex V1 area is related to vision and the
primary motor cortex M1 area is related to motor control. As we can find from
biological data, the density of synaptic interactions within one area is much
higher than that between areas.

####  III-A2 Area-Processes Mapping

Intuitively, a principle for decomposition should be utilizing the areas
described in the Atlas. First, all neurons are partitioned into $n$ parts as
shown in (9), according to areas described in the Atlas. Then, indegree sub-
graphs of each area can be well generated by the partition with (5). After
that, memory consumption of each sub-graph can be estimated, making it easy to
determine how many processes should be mapped to this area. Finally, as to the
name, Area-Processes Mapping is to assign several specific processes to each
area.
Additionally, compared with the naive approach: Random Equivalent Mapping as
shown in 8(a), the advantage of Area-Processes Mapping is to reduce the number
of pre-vertices (pre-synaptic neurons) in indegree sub-graphs. As shown in
Fig. 9, in Random Equivalent Mapping, the pre-vertices ${}^{in}\\!V^{pre}_{i}$
of sub-graph ${}^{in}\\!S(V_{i})=(^{in}\\!V^{pre}_{i},V_{i},^{in}\\!E_{i})$
are randomly selected from $V$. In the worst condition, it can be
${}^{in}\\!V^{pre}_{i}=V$, because the post-vertices in $V_{i}$ are possible
to be pointed by edges from arbitrary pre-vertices in $V$. Then, regarding
data instances, almost all pre-vertices in $V$ should be stored in each
process. Therefore, in Random Equivalent Mapping, the memory consumption
becomes unacceptable in large-scale simulation, with up to billions of neurons
(vertices).
As shown in Fig. 10, in Area-Processes Mapping, a sub-graph $S^{in}(V_{i})$
can be further decomposed into
${}^{in}\\!S^{l}(V_{i})=(V_{i},V_{i},^{in}\\!E^{l}_{i})$ and
${}^{in}\\!S^{r}(V_{i})=(^{in}\\!V^{r}_{i},V_{i},^{in}\\!E^{r}_{i})$, named
local indegree sub-graph and remote indegree sub-graph respectively. And we
have:

${}^{in}\\!S(V_{i})=^{in}\\!S^{l}(V_{i})\veebar^{in}\\!S^{r}(V_{i})=(^{in}\\!V_{i}\cup^{in}\\!V^{r}_{i},V_{i},^{in}\\!E^{l}_{i}\cup^{in}\\!E^{r}_{i})$
(16)

where $V^{r}_{i}\cap V_{i}=\varnothing$, and so $E^{l}_{i}\cap
E^{r}_{i}=\varnothing$.
As we have mentioned above, remote pre-synaptic neurons and remote synaptic
interactions of each area are much less than the local ones as shown in Fig.
8(b), expressed as $n(^{in}\\!V^{r}_{i})\ll n(V_{i})$ and
$n(^{in}\\!E^{r}_{i})\ll n(^{in}\\!E^{l}_{i})$. In other words, most edges are
in ${}^{in}\\!S^{l}(V_{i})$, whose number of post-vertices is fixed to
$n(V_{i})$. Obviously, in large-scale simulations, $n(V_{i})\ll n(V)$, and we
still get $n(V_{i})+n(^{in}\\!V^{r}_{i})\ll n(V)$.
Therefore, with Area-Processes Mapping, the memory consumption of pre-synaptic
neurons and post-synaptic neurons can be kept in a low level, making more
memory available to synaptic interactions, which consumes the most memory.
Utilizing varied densities of synaptic interactions between and inside areas,
we are advancing the available problem size to the next level.

(a) Random Equivalent Mapping

(b) Area-processes Mapping

Figure 8: Mapping Methods for Domain Decomposition: Remote synaptic
interactions with red is much less than local synaptic interactions with
black.

(a) Pre-synaptic Neurons

(b) Post-synaptic Neurons

Figure 9: Random Equivalent Mapping: There are 29 pre-synaptic neurons in
process grape, while the number of post-synaptic neurons is 12, which is a
quarter of total post-synaptic neurons.

(a) Pre-synaptic Neurons

(b) Post-synaptic Neurons

Figure 10: Area-processes Mapping: There are 16 pre-synaptic neurons in
process grape, which has been reduced obviously compared to Random Equivalent
Mapping. And the number of post-synaptic neurons is 12, remaining unchanged.

####  III-A3 Multisection Division with Sampling Method

Area-Processes mapping is our first step in domain decomposition, because each
area might contain millions of neurons and billions of synaptic interactions,
which are still hard to be stored in one process. Fortunately, in the
definition of indegree sub-graphs, edges are bound to post-synaptic neurons,
making it possible to perform a decomposition on post-synaptic neurons only to
generate more indegree sub-graphs. In CORTEX, a stable implementation from
FDPS[24, 25], called Multisection Division with Sampling Method[26], is
introduced. This method can generate divisions in grid cells for a roughly
equal number of points (post-synaptic neurons) with their coordinates in
euclidean space, even in a non-uniform distribution. After the post-synaptic
neurons has been well defined by the generated divisions, an indegree sub-
graph can be generated. As shown in Fig. 11, there are several steps to
generate such divisions and apply them into the original distribution of post-
synaptic neurons with load balance.

(a) Original Distribution

(b) Sample Distribution

(c) Division Generation

(d) Division Application

Figure 11: Multisection Division with Sampling Method

####  III-A4 Load Balance

From we mentioned above, the memory consumption of each process can be
$O(n_{pre}+n_{post}+n_{edges})$, where edges consume most of the memory. From
homomorphism in (8), the partition of edges can be transferred to the
partition on vertices, which is easy to achieve load balance by 2 steps we
mentioned before. First, in Area-Processes Mapping, memory consumption of each
area in the adopted Atlas can be estimated. Then, each area will be mapped to
several processes according to the estimation. After that, each area will be
further decomposed into several sub-graphs by Multisection Division with
Sampling Method, assuming that all properties are homogeneous inside an area.

###  III-B Multi-Threading Parallelization

(a) Increasing Index Order

(b) Increasing Delay Order

(c) Spiking Sub-graph

Figure 12: Data Instance of Synaptic Interactions: Deeper in color means
higher delay of this synaptic interaction.

In general, variables in synaptic interactions and post-synaptic neurons
should be able to be read and written, while pre-synaptic neurons can be read-
only and shared by all threads. Therefore, if a synaptic interaction and its
post-synaptic neuron can be accessed by different threads, mutexes or atomic
operations[12, 13] must be introduced, in order to avoid data racing,
resulting in poor computing throughput. To transcend these barriers, an
elaborate parallelization scheme has been implemented, which can handle
massive processing of synaptic interactions at the thread level, without any
mutex or atomic operation.

####  III-B1 Avoid Data Racing

To avoid data racing, each synaptic interaction and its post-synaptic neuron
should be accessed by a specific thread only. Therefore, the first thing is to
partition the indegree sub-graph again and assign each synaptic interaction
and its post-synaptic neuron to one thread in correspondence. Based on
homomorphism in (8) with (9) and (10), sub-graph generations can be
implemented by dividing all post-synaptic neurons with roughly equal numbers.
After that, we map each vertex and edge to a specific thread as shown in Fig.
13.
Intuitively, Fig. 14(b) is an equivalent representation of Fig. 14(a) in
sparse matrix, where the columns are pre-synaptic neurons, the rows are post-
synaptic neurons and elements are synaptic interactions. Each elements in
sparse matrix (edges in graph) is assigned to a specific thread, as shown by
its color, which is the same thread assigned to its post-synaptic neurons. In
computation, each thread will only access to the synaptic interactions and
post-synaptic neurons assigned with itself, where pre-synaptic neurons can be
shared by all threads, because they are read-only during the computation of
synaptic interactions, as we mentioned above.

(a) Vertices Mapping

(b) Edges Mapping

Figure 13: Vertices and Edges Mapping at the Thread Level

(a) Directed Graph

(b) Sparse Matrix

Figure 14: Sparse Matrix of Synaptic Interactions

As shown in Fig. 12(a), there is the data instance of a sub-graph in memory,
including pre-synaptic neurons, synaptic interactions and post-synaptic
neurons. Except color ’aqua’ for shared accessing to all threads, each color
in each element represents a specific thread with accessing permission to this
element. For synaptic interactions, deeper in color means higher delay, and
the number is the pointer to specific post-synaptic neuron. Another reason for
using such data format is that a spiking sub-graph of (11) can be easily
generated. Because the only necessary information is the spiking pre-synaptic
neurons, the corresponding edges and post-synaptic neurons can be conveniently
found, as shown in Fig. 12(c).

####  III-B2 Varied Synaptic Delays

As we mentioned above, once a pre-synaptic neuron generates a spike, its
synaptic interactions will not take effect on the post-synaptic neurons
immediately, but after a specific delay. However, the synaptic delay varies
between neurons, leading to poor efficiency for checking whether a synaptic
delay falls into the range of current time step. What’s worse, in one
simulator[12, 13], parallel computing is performed on temporal domain. That is
an principled error, because there must be some dependencies between two time
steps.
Fortunately, an innovative scheduling scheme is implemented to deal with
varied synaptic delays. First, synaptic interactions will be reordered
according to their delays and the corresponding threads as shown in Fig.
12(b), where deeper in color means higher delay of this synaptic interaction.
Then, the synaptic interactions will be performed from lowest synaptic delays
to highest synaptic delays according to the current time step. As shown in
Fig. 15, the elements with red border are the currently accessing ones in this
time step, and the elements in color ’aluminium’ are the ones which have been
accessed in the past. By this simple but clear arrangement, parallel
computation of synaptic interactions at the thread level can be well executed
with varied synaptic delays.

(a) Delay Step 1

(b) Delay Step 2

(c) Delay Step 3

Figure 15: Synaptic Interactions with Varied Delays: The elements with red
border are current accessing elements in this time step, and the elements in
color ’aluminium’ are the elements which have been accessed in the past

###  III-C Overlapping of Communication and Computation

As a typical sparse problem, another bottleneck is the inter-process
communication, especially in large-scale simulation. Therefore, a practical
way of general improvement is to overlap communication and computation. In
this section, we will describe why the overlapping is possible, and our
optimized implementation.

####  III-C1 Spikes Broadcast and Buffer

Figure 16: Overlapping of Communication and Computation in One Time Step

As we mentioned in (11), (13) and (14), in indegree sub-graphs, the only
dependencies are pre-synaptic neurons. And also, the spiking graph
${}^{in}\\!S_{s}$ can be defined by spiking pre-synaptic neurons
${}^{in}\\!V^{pre}_{s}$ only, where the post synaptic neurons
${}^{in}\\!E_{s}$ and synaptic interactions ${}^{in}\\!V_{s}$ can be easily
found out in (4). Therefore, the goal of communication is to let all processes
know which pre-synaptic neurons generate spikes in each time step, reducing
lots of requirement for data transfer between processes. This procedure can be
called as Spikes Broadcast, where each process broadcast their spiking pre-
synaptic neurons to others.
Additionally, as we mentioned before, because of varied synaptic delays, the
computation of synaptic interactions from one spiking pre-synaptic neurons
might spans several time steps. So, a buffer is required to place these
spiking pre-synaptic neurons, until their synaptic interactions are all
finished. However, these previous spiking pre-synaptic neurons in buffer make
the overlapping become possible. In each process, the synaptic interactions
from them can be computed before new spiking pre-synaptic neurons arriving
from other processes, to overlap communication and computation, as shown in
Fig. 16.

####  III-C2 Dedicated Thread for Communication

In order to overlap communication and computation as much as possible, a
dedicated thread for communication has been set up in each process, while
others remain unchanged for computation. Different from this implementation
[40] by using OpenMP Tasks for overlapping, the synaptic interactions of
CORTEX are based on the indegree sub-graph of each thread, where the threads
for computation should be specific during the entire simulation. As shown in
Fig. 17, the dataflow is similar to a circulatory system. Neural dynamics with
synaptic interactions will generate spikes from post-synaptic neurons. Then,
the communicating thread broadcasts these spikes to other processes. In the
destination, these spikes will be received and placed in a buffer as spiking
pre-synaptic neurons. After that, all computing threads will fetch these
spiking pre-synaptic neurons to perform synaptic interactions and neural
dynamics without dependencies at both process and thread level. And thus,
there will be new spikes generated from post-synaptic neurons, to start a new
loop.

Figure 17: Dedicated Thread for Communication

##  IV VERIFICATION AND EVALUATION

Because of various modeling methods and biological parameters, it is hard to
demonstrate which one is better in performance or scalability without
simulating a same brain model. However, it is impossible to perform an apple-
to-apple comparison among all simulators, because some of them might not
support some kinds of modeling methods, and there should be lots of work to
build the same simulation on all existing platforms. So, to perform a fair and
meaningful comparison, cases with same modeling method will be simulated by
our simulator: CORTEX and NEST Simulator, the SOTA[15] with available open-
source code. Two cases will be in comparison as follows:

  * •

Verification: A random balanced network with STDP mechanism in nonlinear
synaptic dynamics, provided by NEST[38].

  * •

Evaluations: A Multi-scale spiking network model of marmoset in paxinos
structural connectome, built from the connectome data at the website:
https://connectome.marmosetbrainmapping.org/Paxinos_SC_matrix, cell density:
https://www.marmosetbrain.org/cell_density, interareal distance:
https://analytics.marmosetbrain.org/static/data/marmoset_brain_connectivity_1_0_interareal_distance_matrix.txt
and the internal architecture is derived from[30]:
https://github.com/nest/nest-
simulator/tree/master/pynest/examples/Potjans_2014

###  IV-A Verification

Each aspect we mentioned above has its own meaning. The first case is a
balanced random network in which the excitatory-excitatory neurons exhibit
STDP with multiplicative depression and power-law potentiation. The number of
incoming synaptic interactions per neuron is fixed and independent of network
size. By this case, we want to show that CORTEX is able to support nonlinear
synaptic interactions, including complex computation with varied data
structures like queue, stack and so on, and still free from data racing
without any mutex or atomic operation. Also, the results of thread mapping
will be check in this case. Simply, if a edge or post-vertex is accessed by
different threads, Abort will be called by CORTEX.
As described in NEST[38], the neurons’ firing rates of this network should be
lower than 10 Hz, where the results might differ due to the chaotic nature of
the dynamics. On different simulators or runtime platforms, minor differences
in the computation of neural and synaptic interactions can lead to completely
different spike sequences after a short time[39].
Although this case is named as ”HPC Benchmark” by NEST, we don’t use it as a
real performance benchmark, because this case just represents a simple
theoretical SNN without enough biological properties from brain architecture
in reality, especially the varied density and delays of synaptic interactions,
upon which many of the optimization efforts we mentioned above are focused.

###  IV-B Evaluation

Using LIF neuron model[21, 22], the marmoset cerebral cortex simulation is a
typical large-scale SNNs built from biological data, which can show the
superiority of CORTEX in terms of not only application-available problem size
but also computing performance. In our belief, this superiority of CORTEX is
not limited to a specific case as we mentioned above, almost every brain
architecture driven by biological data can enjoy CORTEX. Unfortunately, up to
now, many existing cases[27, 28, 29] from other simulators can’t run on the
Fugaku supercomputer due to too many technical problems with system
environment and software dependencies. Also, another reason is that the
dynamic response and computing intensity of this brain architecture remain
unchanged with different scaling problem sizes. Therefore, it is suitable to
be as a baseline, to demonstrate the performance between different problem
sizes.
Even so, the source code of CORTEX will be available on Github, and everyone
is welcome to build any brain architecture on CORTEX and to make comparisons
with other simulators on leading-edge supercomputers, where the neuron and
synapse templates along with connecting definitions will be updated to support
a wide range of modeling methods in future.
The performance is measured using the average time spent on the benchmark test
as we mentioned above. The comparison between CORTEX and NEST Simulator
includes memory consumption and simulation time, exhibiting the general
performance in multi-scaling problem size, as shown in Fig. 18. The normalized
problem size of 1 contains 1 million neurons with 3.8 billion synaptic
interactions, and there are 4 processes per node. Also, all variables for
numerical computing are in IEEE 754 64-bit floating point format without any
compression on accuracy.

###  IV-C Results and Visualization

(a) Memory Usage

(b) Simulation Time

Figure 18: Performance Result: Memory usage is the maximal value of
consumption in each node. The normalized problem size 1 contains 1 million
neurons with 3.7 billion synaptic interactions. All variables for numerical
computing are in IEEE 754 64-bit floating point format without any compression
on accuracy.

Because of the chosen LIF model [21, 22], whose computing intensity is much
lower than that of other models, the scalability might not be as good as shown
in other publications[32] using Hodgkin-Huxley model[31] with much higher
computing intensity. Also, although it is possible to perform larger
simulation on CORTEX, the gap of memory consumption and computing performance
between CORTEX and NEST Simulator within 1536 processes (384 nodes) has been
huge enough.
Two raster results of area V1 are shown in Fig. 19, which are similar to each
other with slight differences, because the computing procedures and random
generations between these two simulators are not exactly the same.

(a) CORTEX Simulator

(b) NEST Simulator

Figure 19: Rasters Plots of Cortical Activities

##  V DISCUSSION

Despite introducing slightly greater complexities to our algorithmic
framework, a key point we wish to highlight is the impact of these innovations
on scaling simulation performance and capacity. By maintaining numerical
integrity, we overcome the limitations imposed by sparsity. This approach
significantly advances our ability to support simulations of the entire human
brain, marking a substantial improvement over previous methodologies. Our
successful application of CORTEX on the Fugaku Supercomputer demonstrates its
potential to handle increasingly larger problems, moving us closer to the
ambitious goal of a full human-scale brain simulation, the pinnacle challenge
in our field.

####  V-1 Domain Specific Optimizations

On the practical side, there are still many ways to improve CORTEX at various
aspects, especially with domain specific optimizations. In future, CORTEX will
be further implemented with optimizations on specific architectures, which are
not limited in homogeneous CPU-based systems. Admittedly, implementation on
heterogeneous system like CUDA GPU[42] or MN-Core[43] accelerators should be
more challenging, directly facing the sparsity of brain architecture in both
spatial and temporal domains. Therefore, we expect that further improvements
on both the algorithmic and engineering considerations can trigger a
breakthrough in future.

####  V-2 Further Improvement on Communication

In brain simulation, there are lots of irregular spike broadcasts in each time
step. Obviously, MPI_Bcast are not suitable for irregular large amounts of
communication. Brain Simulation Broadcast (BSB) is a broadcast acceleration
library specifically designed for this communication pattern, which can
automatically packs/unpacks spikes into/from messages and adaptively routes
the messages among processes to decrease the number of small messages in the
physical network. To decrease the communication latency, the interfaces of BSB
are designed in a producer-consumer model, and its implementation utilizes
Unified Notifiable RMA (UNR) library to achieve synchronization-free
communication. We will optimize the communication of CORTEX with BSB in the
next update.

####  V-3 Towards Whole Human-Scale Brain Simulation

As we mentioned at the outset, human-scale brain simulation is one of the most
ambitious scientific challenges of the 21st century, where the whole brain
includes the cerebellum and basal ganglia, with synaptic plasticity, whose
properties are much more different from the architecture of cerebrum cortex
only. Within this brain architecture, minor differences in parameters and
mechanisms can lead to completely different firing rates, synaptic density and
mechanisms, related to differences in computing intensity and memory
consumption. Obviously, with such a heterogeneous brain architecture, load
balance of both computation and storage is more challenging to achieve at the
same time.

##  VI ACKNOWLEDGMENTS

This work was supported by Program for Promoting Researches on the
Supercomputer Fugaku (hp200139). We are very grateful to Dr. Jun Igarashi and
Prof. Kenji Doya.

## References

  * [1] Amunts, Katrin, et al. ”Julich-Brain: A 3D probabilistic atlas of the human brain’s cytoarchitecture.” Science 369.6506 (2020): 988-992.
  * [2] Frackowiak, R., and Markram, H. (2015). The future of human cerebral cartography: a novel approach. Philos. Trans. R. Soc. Lond. B Biol. Sci. 370:20140171.
  * [3] Churchland, P. S. (1986). Neurophilosophy: Toward a Unified Science of the Mind-Brain. Cambridge, MA: MIT Press.
  * [4] DeFelipe, J. (2015). The anatomical problem posed by brain complexity and size: a potential solution. Front. Neuroanat. 9:104.
  * [5] Fan, Xue, and Henry Markram. “A Brief History of Simulation Neuroscience.” Frontiers in Neuroinformatics 13 (May 7, 2019): 32.
  * [6] Yamazaki, Tadashi, Jun Igarashi, and Hiroshi Yamaura. “Human-Scale Brain Simulation via Supercomputer: A Case Study on the Cerebellum.” Neuroscience 462 (May 2021): 235–46.
  * [7] November. (2023). TOP500 supercomputer sites. Available online at: https://www.top500.org/lists/top500/2023/11/ [Accessed March 9, 2024].
  * [8] Billeh, Yazan N., Binghuang Cai, Sergey L. Gratiy, Kael Dai, Ramakrishnan Iyer, Nathan W. Gouwens, Reza Abbasi-Asl, et al. “Systematic Integration of Structural and Functional Data into Multi-Scale Models of Mouse Primary Visual Cortex.” Neuron 106, no. 3 (May 2020): 388-403.e18.
  * [9] Kandel ER, Schwartz JH, Jessell TM, Mack S, editors. Principles of neural science. McGraw-Hill Medical.
  * [10] H. Markram, E. Muller, S. Ramaswamy, M. W. Reimann, M. Abdellah, C. A. Sanchez, A. Ailamaki, L. Alonso-Nanclares, N. Antille, S. Arsever et al., “Reconstruction and simulation of neocortical microcircuitry,” Cell, vol. 163, no. 2, pp. 456–492, 2015.
  * [11] Gerstner, Wulfram, Werner M. Kistler, Richard Naud, and Liam Paninski. Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition. 1st ed. Cambridge University Press, 2014.
  * [12] P. Qu, H. Lin, M. Pang, X. Liu, W. Zheng and Y. Zhang, ”ENLARGE: An Efficient SNN Simulation Framework on GPU Clusters,” in IEEE Transactions on Parallel and Distributed Systems, vol. 34, no. 9, pp. 2529-2540, Sept. 2023.
  * [13] Qu, Peng, Youhui Zhang, Xiang Fei, and Weimin Zheng. “High Performance Simulation of Spiking Neural Network on GPGPUs.” IEEE Transactions on Parallel and Distributed Systems 31, no. 11 (November 1, 2020): 2510–23.
  * [14] Fugaku Supercomputer — RIKEN Center of Computational Science: https://www.r-ccs.riken.jp/en/fugaku/ [Accessed March 9, 2024].
  * [15] Kuniyoshi, Yusuke, Rin Kuriyama, Shu Omura, Carlos Enrique Gutierrez, Zhe Sun, Benedikt Feldotto, Ugo Albanese, et al. “Embodied Bidirectional Simulation of a Spiking Cortico-Basal Ganglia-Cerebellar-Thalamic Brain Model and a Mouse Musculoskeletal Body Model Distributed across Computers Including the Supercomputer Fugaku.” Frontiers in Neurorobotics 17 (October 5, 2023): 1269848.
  * [16] Izhikevich, Eugene M., and Gerald M. Edelman. “Large-Scale Model of Mammalian Thalamocortical Systems.” Proceedings of the National Academy of Sciences 105, no. 9 (March 4, 2008): 3593–98.
  * [17] Deo, Narsingh (1974). Graph Theory with Applications to Engineering and Computer Science. Englewood, New Jersey: Prentice-Hall. ISBN 0-13-363473-6.
  * [18] KarypisG, Kumar V (1998) A fast and high quality multi level scheme for partitioning irregular graphs. SIAM J Sci Comput 20(1):359–392
  * [19] Feldmann, Andreas Emil; Foschini, Luca (2012). ”Balanced Partitions of Trees and Applications”. Proceedings of the 29th International Symposium on Theoretical Aspects of Computer Science: 100–111.
  * [20] “The Brain/MINDS Population-Based Marmoset Brain Atlas,”
https://dataportal.brainminds.jp/3d-atlas-viewer-bma2023

  * [21] Rotter, Stefan, and Markus Diesmann. “Exact Digital Simulation of Time-Invariant Linear Systems with Applications to Neuronal Modeling.” Biological Cybernetics 81, no. 5–6 (November 5, 1999): 381–402.
  * [22] Brette, Romain, and Wulfram Gerstner. “Adaptive Exponential Integrate-and-Fire Model as an Effective Description of Neuronal Activity.” Journal of Neurophysiology 94, no. 5 (November 2005): 3637–42.
  * [23] Morrison, Abigail, Ad Aertsen, and Markus Diesmann. “Spike-Timing-Dependent Plasticity in Balanced Random Networks.” Neural Computation 19, no. 6 (June 2007): 1437–67.
  * [24] Iwasawa, Masaki, Ataru Tanikawa, Natsuki Hosono, Keigo Nitadori, Takayuki Muranushi, and Junichiro Makino. “FDPS: A Novel Framework for Developing High-Performance Particle Simulation Codes for Distributed-Memory Systems.” In Proceedings of the 5th International Workshop on Domain-Specific Languages and High-Level Frameworks for High Performance Computing, 1–10. Austin Texas: ACM, 2015.
  * [25] Iwasawa, Masaki, Ataru Tanikawa, Natsuki Hosono, Keigo Nitadori, Takayuki Muranushi, and Junichiro Makino. “Implementation and Performance of FDPS: A Framework Developing Parallel Particle Simulation Codes.” Publications of the Astronomical Society of Japan 68, no. 4 (August 1, 2016): 54.
  * [26] Ishiyama, Tomoaki, Keigo Nitadori, and Junichiro Makino. “4.45 Pflops Astrophysical N-Body Simulation on K Computer – The Gravitational Trillion-Body Problem.” In 2012 International Conference for High Performance Computing, Networking, Storage and Analysis, 1–10. Salt Lake City, UT: IEEE, 2012.
  * [27] Schmidt, Maximilian, Rembrandt Bakker, Claus C. Hilgetag, Markus Diesmann, and Sacha J. Van Albada. “Multi-Scale Account of the Network Structure of Macaque Visual Cortex.” Brain Structure and Function 223, no. 3 (April 2018): 1409–35.
  * [28] Schuecker J, Schmidt M, van Albada SJ, Diesmann M and Helias M (2017) Fundamental Activity Constraints Lead to Specific Interpretations of the Connectome. PLOS Computational Biology, 13(2): e1005179.
  * [29] A multi-scale layer-resolved spiking network model of resting-state dynamics in macaque cortex. PLOS Computational Biology, 14(9): e1006359.
  * [30] Potjans, Tobias C., and Markus Diesmann. “The Cell-Type Specific Cortical Microcircuit: Relating Structure and Activity in a Full-Scale Spiking Network Model.” Cerebral Cortex 24, no. 3 (March 2014): 785–806.
  * [31] Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544.
  * [32] Lu, Wenlian, Longbin Zeng, Xin Du, Wenyong Zhang, Shitong Xiang, Huarui Wang, Jiexiang Wang, et al. “Digital Twin Brain: A Simulation and Assimilation Platform for Whole Human Brain.” arXiv, August 2, 2023.
  * [33] R. Ananthanarayanan, S. K. Esser, H. D. Simon, and D. S. Modha, “The cat is out of the bag: cortical simulations with 109 neurons, 1013 synapses,” in Proceedings of the conference on high performance computing networking, Storage Anal, vol. 63, pp. 1–12, 2009.
  * [34] M. Helias, S. Kunkel, G. Masumoto, J. Igarashi, J. M. Eppler, S. Ishii, T. Fukai, A. Morrison, and M. Diesmann, “Supercomputers ready for use as discovery machines for neuroscience,” Frontiers in Neuroinformatics, vol. 6, p. 26, 2012.
  * [35] S. Kunkel, M. Schmidt, J. M. Eppler, H. E. Plesser, G. Masumoto, J. Igarashi, S. Ishii, T. Fukai, A. Morrison, M. Diesmann, and M. Helias, “Spiking network simulation code for petascale computers,” Frontiers in Neuroinformatics, vol. 8, p. 78, 2014.
  * [36] J. Jordan, T. Ippen, M. Helias, I. Kitayama, M. Sato, J. Igarashi, M. Diesmann, and S. Kunkel, “Extremely scalable spiking neuronal network simulation code: from laptops to exascale computers,” Frontiers in Neuroinformatics, vol. 12, p. 2, 2018.
  * [37] J. Igarashi, H. Yamaura, and T. Yamazaki, “Large-scale simulation of a layered cortical sheet of spiking network model using a tile partitioning method,” Frontiers in Neuroinformatics, vol. 13, p. 71, 2019.
  * [38] “Random balanced network HPC benchmark,” NEST Simulator
https://github.com/nest/nest-simulator/blob/master/pynest/examples/hpc_
benchmark.py

  * [39] “NEST Simulator Documentation,” NEST Simulator
https://nest-simulator.readthedocs.io/en/stable/faqs/faqs.html#faqs-precise-
neurons

  * [40] Valero-Lara, Pedro, Raül Sirvent, Antonio J. Peña, and Jesús Labarta. “MPI+OpenMP Tasking Scalability for Multi-Morphology Simulations of the Human Brain.” Parallel Computing 84 (May 2019): 50–61.
  * [41] Y. Ajima, T. Kawashima, T. Okamoto, N. Shida, K. Hirai, T. Shimizu, S. Hiramoto, Y. Ikeda, T. Yoshikawa, K. Uchida, and T. Inoue, “The tofu interconnect D,” in Proceedings of the IEEE International Conference on Cluster Computing ICCC 2018-September, (Belfast: IEEE), pp. 646–654, 2018.
  * [42] NVIDIA H200 Tensor Core GPU — NVIDIA Corporation: https://www.nvidia.com/en-us/data-center/h200/
  * [43] MN-Core™ Series — Preferred Networks, Inc.: https://projects.preferred.jp/mn-core/en/ [Accessed March 9, 2024].