[1]Johannes H. Jensen These authors contributed equally to this work.

[2]Anders Strømberg These authors contributed equally to this work.

[1]Department of Computer Science, Norwegian University of Science and
Technology, Trondheim, Norway

[2]Department of Electronic Systems, Norwegian University of Science and
Technology

3]ALBA Synchrotron Light Facility, Carrer de la Llum 2 – 26, Cerdanyola del
Vallés, 08290, Barcelona, Spain

# Clocked dynamics in artificial spin ice

<EMAIL_ADDRESS><EMAIL_ADDRESS>Ida Breivik  Arthur Penty
Miguel Angel Niño  Muhammad Waqas Khaliq  Michael Foerster  Gunnar Tufte  Erik
Folven  *  *  [

###### Abstract

Artificial spin ice (ASI) are nanomagnetic metamaterials exhibiting a wide
range of emergent properties, which have recently shown promise for
neuromorphic computing. However, the lack of efficient protocols to control
the state evolution of these metamaterials has been limiting progress. To
overcome this barrier, we introduce _astroid clocking_ , a global field
protocol offering discrete, gradual evolution of spin states. The method
exploits the intrinsic switching astroids and dipolar interactions of the
nanomagnets to selectively address ASI spins in sequence. We demonstrate,
experimentally and in simulations, how astroid clocking of pinwheel ASI allows
ferromagnetic domains to be gradually grown or reversed at will. More complex
dynamics arise when the clock protocol allows both growth and reversal to
occur simultaneously. Astroid clocking offers unprecedented control and
understanding of ASI dynamics in both time and space, extending what is
possible in nanomagnetic metamaterials.

###### keywords:

Artificial Spin Ice, Magnetic Metamaterial, Pinwheel Artificial Spin Ice,
Permalloy, Nanofabrication, Nanomagnets, Coupled Nanomagnetic Ensembles,
Astroid Clocking, Unconventional Computation, Material computation, flatspin

##  1 Controlling artificial spin ice

Artificial spin ice (ASI) are systems of coupled nanomagnets arranged on a
two-dimensional lattice. The nanomagnets are elongated, giving them two stable
magnetization directions, thus behaving as artificial spins. Dipolar
interactions give rise to a rich variety of emergent behavior, as determined
by the ASI geometry[1, 2, 3]. As this behavior can be probed directly, ASIs
have attracted considerable interest as model systems for the study of
fundamental physics[4, 5]. More recently, ASIs have shown promise as
substrates for computation[6, 7, 8, 9, 10, 11, 12].

External fields are the primary method used to perturb ASIs in a controlled
manner. Various global field protocols have been employed. For example, a
cycled in-plane field is often used to characterize magnetization reversal[1,
13, 14, 15, 16, 17, 18, 19, 20, 21]. Another approach is to use a rotating
field with slowly decreasing amplitude to effectively anneal the ASI to a low
energy state[22, 23, 24, 25, 26, 27, 28, 29]. While there are variations of
these simple field protocols, more complex protocols are largely unexplored.

These approaches use field strength to modulate ASI behavior, which will
typically result in uncontrolled avalanches of activity[30]. An in-plane field
will advance ASI state primarily when the strength of the field is increased
beyond the coercivity of the array, and is highly dependent on field
resolution[25, 29, 19]. Consequently, the discrete spin flip dynamics in the
ASI are sudden and hard to control.

Here, we introduce a new field protocol scheme called _astroid clocking_ ,
which produces fundamentally different spin flip dynamics. Astroid clocking of
an ASI results in a step-wise, gradual evolution of spin states. This offers
unprecedented control and understanding of the dynamical process in both time
and space. By exploiting the shape and orientation of the nanomagnet switching
astroids and their dipolar coupling, specific field angles are employed to
selectively address different parts of the ensemble. A clock protocol pulses
fields at these angles in an alternate fashion, driving the intrinsic dynamics
of the ASI. Distinctively, the clock pulses maintain a constant field
amplitude.

In the context of nanomagnetic logic, Nomura et al. [31] demonstrated how the
shape of two overlapping Stoner-Wohlfarth astroids can be exploited to
preferentially switch nanomagnets in a 1D shift register. Astroid clocking
extends and generalizes this concept to 2D nanomagnet arrays and non-
elliptical nanomagnets with different astroid shapes. We show how astroid
clocking reveals the intrinsic dynamics of coupled nanomagnetic systems.

For this study we consider the pinwheel ASI system, but stress that astroid
clocking is readily applicable to other coupled nanomagnetic systems as well.
We demonstrate and analyse how ferromagnetic domains in pinwheel ASI can be
gradually grown and reversed at will using astroid clocking. Different clock
protocols are explored, giving rise to distinct properties of the spin flip
dynamics.

##  2 Astroid clocking

Figure 1: Astroid clocking of pinwheel ASI. a, a small $4\times 4$ pinwheel
ASI, formed by two interleaved sublattices $L_{a}$ (solid outline) and $L_{b}$
(dashed outline) with magnets oriented at $+45\degree$ and $-45\degree$,
respectively. Colors correspond to magnetization direction as indicated by the
white center arrows. The magnetic state shows the four possible ferromagnetic
domains of pinwheel ASI, where the net magnetization forms a counter-clockwise
magnetic flux closure pattern. b, switching astroids of the magnets in
sublattice $L_{a}$ (solid lines) and $L_{b}$ (dashed lines), along with the
four clock fields, $\mathbf{H}_{A}$, $\mathbf{H}_{B}$, $\mathbf{H}_{a}$, and
$\mathbf{H}_{b}$. The astroid edges are colored according to the magnet state
which is promoted when fields cross the edge. Similarly, the colored regions
correspond to fields that _exclusively_ promote a magnet state within a
sublattice. Astroid axes are normalized with respect to the hard axis
switching threshold, $h_{k}$.

Pinwheel ASI[2, 32, 19] consists of nanomagnets arranged on two interleaved
square sublattices, as shown in Fig. 1a. In this study, the magnets in the two
sublattices $L_{a}$ and $L_{b}$ are rotated $+45\degree$ and $-45\degree$ with
respect to the array edges. The sublattice and magnetization of the magnets
are indicated by their color: magnets in sublattice $L_{a}$ are orange or
blue, while magnets in sublattice $L_{b}$ are pink or green. For brevity, we
will refer to magnet state by these four colors.

Pinwheel ASI favors a ferromagnetic ordering, with emergent domains of
coherent magnetization. Fig. 1a shows the four possible domain directions:
rightwards (orange/pink), leftwards (blue/green) and so on. The ferromagnetic
domains are separated by domain walls, which are slightly less energetically
favorable[32].

The switching threshold of a nanomagnet depends on the field angle, and can be
approximated by the Stoner-Wohlfarth astroid. Fig. 1b shows the switching
astroids for the two orientations of stadium-shaped magnets in pinwheel
ASI[33]. A magnet will switch state if the total field acting on it lies
outside the astroid boundary, _and_ the field is directed against the current
magnetization.

Nanomagnet shape largely determines the shape of the astroid. Stadium-shaped
nanomagnets, commonly used in ASI, have a switching astroid with 2-fold
rotational symmetry[33]. This is in contrast to classical Stoner-Wohlfarth
astroids that display 4-fold rotational symmetry derived for elliptical
nanomagnets[34].

Switching astroids that break the 4-fold rotational symmetry, can be exploited
to selectively address nanomagnets that are rotated $90\degree$ with respect
to each other. If the total field lies within the shaded regions in Fig. 1b,
_only_ the nanomagnets in the corresponding sublattice will be able to switch.
A field in the orange/blue shaded regions will address only the magnets in
sublattice $L_{a}$, while a field in the pink/green regions will address only
magnets in $L_{b}$. Furthermore, each region promotes a specific magnet state
within each sublattice, e.g., a field in the blue shaded region promotes blue
magnets by switching orange magnets.

In this study, we define two _bipolar clocks_ A and B along the $+22\degree$
and $-22\degree$ axes, respectively. As shown in Fig. 1b, each clock consists
of a positive and negative _clock field_ of magnitude $H$ along the clock
axis. The four arrows in Fig. 1b are colored according to the magnet states
they promote, e.g., the $\mathbf{H}_{A}$ field only promotes orange magnets.

The dipolar fields $\mathbf{h}_{\text{dip}}$ from neighboring magnets may
either promote or prevent switching. If the dipolar fields are directed out of
(into) the astroid, they effectively promote (prevent) switching. A clock
field can thus selectively address a subset of a sublattice, depending on the
state of the ensemble. The clock angles $\pm 22\degree$ are selected to allow
the dipolar fields to have a large influence on switching, using a field
strength $H$ close to the switching threshold. However, a precise angle is not
crucial and the method tolerates a wide range of clock angles. Our system
tolerates clock angles in the range
$10\text{\,}\mathrm{\SIUnitSymbolDegree}35\text{\,}\mathrm{\SIUnitSymbolDegree}$,
and field strengths accurate to $3\text{\,}\mathrm{mT}4\text{\,}\mathrm{mT}$.

Figure 2: Clock diagram of astroid clocking. Clock protocols are defined by
sequences of clock pulses. The clock diagram shows $AB$ clocking (alternating
pulses of the positive clock fields $\mathbf{H}_{A}$ and $\mathbf{H}_{B}$)
followed by $ab$ clocking (alternating pulses of the negative clock fields
$\mathbf{H}_{a}$ and $\mathbf{H}_{b}$).

Fig. 2 illustrates astroid clocking, where a _clock pulse_ is defined as
ramping a clock field from zero to $H$ and down to zero again. The ramping
speed is much slower than the timescale of nanomagnetic switching. A _clock
protocol_ is a specific sequence of clock pulses. For example, $AB$ clocking
consists of repeated alternating clock pulses of $A$ and $B$. We define a
_clock cycle_ as a single sequence of the clock pulses in a protocol, e.g., an
$aAbB$ clock cycle is the sequence of four pulses $(a,A,b,B)$. A _unipolar_
clock protocol exclusively employs one polarity of each clock, while a
_bipolar_ clock protocol employs both polarities.

##  3 Unipolar clocking

First, we explore the spin flip dynamics of pinwheel ASI when subject to the
unipolar clock protocols $AB$ and $ab$. The $50\times 50$ pinwheel ASI (5100
magnets) is initialized with a small rightwards (orange/pink) domain in the
center of an otherwise leftwards polarized (blue/green) array. Fig. 3 (1)
shows a closeup of the initial state.

Figure 3: Simulation of unipolar astroid clocking of pinwheel ASI in flatspin.
Each snapshot shows a zoomed-in view of the $50\times 50$ nanomagnet system,
at different points during a clock protocol. (1) shows the initial state, a
small orange/pink (rightwards) domain in the center of an otherwise polarized
blue/green (leftwards) array. (2-8) show the state during $AB$ clocking,
resulting in gradual domain growth. (9-12) show the subsequent states during
$ab$ clocking, resulting in gradual domain reversal. Magnets that change state
between snapshots are highlighted with a solid black outline.

Fig. 3 (2-8) shows the state evolution of the array subject to $AB$ clocking,
obtained from flatspin simulations (see flatspin simulations). As expected,
the $A$ pulse selectively switches magnets in sublattice $L_{a}$ from blue to
orange, while the $B$ pulse selectively switches magnets in sublattice $L_{b}$
from green to pink.

Interestingly, the particular magnets that switch are the ones along the
domain wall. As a result, the inner (leftwards) domain grows gradually over
time, with only a thin layer of the domain advancing after each clock pulse.
The growth is _monotonic_ and _step-wise_ , driven by the clock pulses.

A curious property is that the domain grows mainly in the horizontal
direction. In Fig. 3 (2-8), the magnets along the vertical domain walls are
the only ones to switch. If the growing domain reaches the edges of the array,
the direction of growth changes and becomes vertical, eventually filling the
entire array (Fig. 7).

Inverting the clock pulses ($ab$ clocking), will instead grow the outer
(blue/green) domain and consequently shrink (reverse) the inner (orange/pink)
domain. As can be seen in Fig. 3 (9-12), domain reversal from (8) proceeds in
both vertical and horizontal directions, resulting in reversal of the inner
domain in fewer clock cycles compared to growth. Hence there is an apparent
_asymmetry_ in the direction of domain growth and reversal.

##  4 Growth and reversal mechanism

To understand the mechanism behind the domain growth and reversal, we consider
the larger domain shown in Fig. 4c, subject to $\mathbf{H}_{A}$. In Fig. 4a,
we plot the relative locations of all the magnets within their respective
switching astroids. Each dot represents the total field
$\mathbf{h}_{i}=\mathbf{H}_{A}+\mathbf{h}_{\text{dip}}^{(i)}$ experienced by a
magnet $i$ in its _local frame of reference_. There are four clusters of dots
within the astroids, corresponding to the four magnet colors, where only the
blue magnets are close to switching.

Figure 4: Astroid clusters and astroid (black curve) for the pinwheel system
shown in c, when subject to the clock field $\mathbf{H}_{A}$. a, astroid
cluster plot where each dot represents the total field
$\mathbf{h}_{i}=\mathbf{H}_{A}+\mathbf{h}_{\text{dip}}^{(i)}$ experienced by a
magnet $i$, projected onto its parallel ($h_{\parallel}$) and perpendicular
($h_{\perp}$) axis. The colors in the plot correspond to magnet state. Note
that the astroid plot shows location relative to each magnet’s own switching
threshold, e.g., orange magnets are far from switching as they are aligned
with $\mathbf{H}_{A}$. b, closeup of the blue cluster, revealing a sub-group
of blue magnets that lie outside the switching astroid and are eligible for
switching. These magnets are highlighted in c, and are all found to lie along
the vertical and $+45\degree$ domain walls. The data is obtained from flatspin
simulations.

The internal structure of each astroid cluster is a result of the nanomagnet
dipolar coupling, and a direct consequence of the ASI geometry. In the absence
of dipolar fields, each cluster collapses into a single point. The dipolar
fields add complex structure to the clusters, with sub-groups corresponding to
different subsets of magnets within the ASI. For a detailed analysis of
neighbor contributions, see section S1.

The inset shown in Fig. 4b reveals the structure of the blue cluster. Notice
there are a few blue dots that lie outside the astroid, corresponding to
magnets that are eligible for switching, which are highlighted in Fig. 4c.
Evidently, the switchable magnets all lie along the vertical and $+45\degree$
domain walls.

When a magnet switches, its location within the astroid jumps to the cluster
of opposite spin, e.g., a blue magnet switches to the orange state. In
addition, neighboring magnets will see a change in the dipolar fields, causing
movement within their respective clusters. In this way, the switching of a
magnet may enable future switching in neighboring magnets, either during the
current or a future clock pulse.

The observed horizontal domain growth can now be explained from the internal
structure of the astroid clusters. We have seen that magnets along certain
types of domain walls can be selectively switched under an applied clock
field. Switching the blue (highlighted in Fig. 4c) magnets along these domain
walls reverses their dipolar fields, which affects the structure of the green
cluster. Consequently, green magnets that are part of the domain walls will
approach the switching astroid. When the $\mathbf{H}_{B}$ field is
subsequently applied, these magnets will be outside the astroid and hence
switch. As this cycle repeats, the result is an apparent horizontal domain
growth emerging from the dipolar interactions and clock fields.

During domain reversal, both horizontal and vertical domain walls take part in
the process. As a result, reversal requires fewer clock cycles compared to
growth. During reversal, the switchable magnets lie along both the horizontal,
vertical and $-45\degree$ domain walls (Fig. 8).

We find that the horizontal domain wall movement, particular to reversal, is
dependent on the curvature of the reversing domain. If the horizontal domain
wall is surrounded by a blue/green domain on three sides, there is a stronger
dipolar “push” towards the astroid edge. As such, domain shape plays a crucial
role in the reversal process.

When a domain grows to reach the edges of the array, there is an apparent
transition from horizontal to vertical growth (Fig. 7). We find that vertical
growth proceeds by avalanches along the domain wall, starting at the bottom-
left and top-right corners of the domain, close to the array edges.

##  5 Experimental growth and reversal

Next, we demonstrate astroid clocking of pinwheel ASI experimentally. Samples
are imaged with x-ray magnetic circular dichroism photoemission electron
microscopy (XMCD-PEEM), with an in-situ vector magnet to perform astroid
clocking. See Methods for details.

After polarizing all magnets in the leftwards direction (bright contrast), we
perform steps of $AB$ clocking, imaging in-between each clock cycle. Fig. 5a
shows total magnetization of the ensemble, obtained from the XMCD-PEEM images,
which increases in a stable, monotonic fashion. Selected experiment snapshots
are shown in Fig. 5b. Snapshots (1-3) show that domains nucleate at the
vertical edges then predominantly expand horizontally.

Figure 5: Results of growth and reversal with unipolar clock protocols, and
control experiment. a, total magnetization of the ensembles subject to the
different clock protocols. The timeline indicates clock time, labeled by the
clock protocol. During $AB$ clocking, the ensembles undergo growth and hence
an increase in magnetization. The second phase, $ab$ clocking, quickly
reverses domains and total magnetization. The control experiment, consisting
of separate $A$ and $B$ clock sequences, show no development of the domains.
b, magnetic image snapshots (experiemental XMCD-PEEM images and flatspin
simulated XMCD-PEEM contrast images) of the ensembles at the specified points
in time. The depicted ensembles are approximately 12.5x12.5µ ($50\times 50$
pinwheel ASI, 5100 magnets). All XMCD-PEEM images are available in Figs. 9 and
10. Videos of the experiment and simulation are provided in
LABEL:vid:unipolar-peem and LABEL:vid:unipolar-sim.

Domain formation at the vertical array edges can be explained by the dipolar
field-driven mechanism behind $AB$ clocking. While domain nucleation along
horizontal edges is possible, continued growth primarily occurs in the
horizontal direction, preventing further expansion of horizontal edge
nucleated domains.

In any physical ASI system, the nanomagnets will exhibit a range of intrinsic
switching thresholds, a disorder, due to imperfections and microscopic
variations of material composition. Disorder affects both domain shape and
growth dynamics, as evident in our experimental results. Compared to the
idealized simulations, domains appear more organic, with distinct features
such as jagged edges, slanted domain walls, and sporadic holes. In terms of
dynamics, some domain borders get stuck for several clock cycles, while others
advance more than one step during a single cycle (see Fig. 9).

By introducing disorder to the simulations (see flatspin simulations), we
obtain results that more closely resemble the experiment. The magnetization
curve and snapshots from simulations with disorder are included in Fig. 5.
Notice how the simulated snapshots show organic-looking domains that resemble
the domains of the experiment.

After growth, we apply the reversal clock protocol, $ab$ clocking. For each
$ab$ clock cycle, the magnetization reduces sharply, with domains shrinking
more rapidly compared to the increase during growth. Comparing snapshot 3 and
4 of Fig. 5b, it is clear that the domains shrink in both vertical and
horizontal directions.

Next, we conduct a control experiment to verify that simply repeating a clock
pulse $A$ or $B$ does not result in domain growth. After re-initializing the
system, we apply several pulses of $A$, then several pulses of $B$, imaging
after each pulse. As seen in Fig. 10 and the last part of Fig. 5a, only the
first application of $A$ or $B$ results in growth. Growth progresses only when
the type of clock pulse is changed, which confirms that the alternating
pattern of $A$ and $B$ is what drives the observed domain growth.

These experiments affirm the viability of astroid clocking in the face of
experimental sensitivities (as low as $<1\text{\,}\mathrm{mT}$ from Fig. 4)
and potential impediments such as fabrication imperfections, temperature
effects, and material degradation. While unstable individual magnets and
inaccuracies in the image analysis induce some noise, it is negligible
compared to the effect of astroid clocking. Experimental astroid clocking is
surprisingly robust, demonstrating that it is possible to precisely control
the spin flip dynamics of ASIs using global fields.

##  6 Bipolar clocking

In bipolar clocking, each clock may be pulsed in both polarities. We consider
two clock protocols illustrated in Fig. 11, namely $aAbB$ and its inverse,
$AaBb$ clocking. In contrast to unipolar clocking, the magnetic fields in
these bipolar clock protocols are balanced, i.e., the sum of all clock fields
is zero. One might then expect that this results in a net zero magnetization
change.

On the contrary, bipolar clocking also results in domain growth and reversal,
and a net change in magnetization. Fig. 6a plots the total magnetization of
pinwheel ASI subject to bipolar clocking. As can be seen, $aAbB$ clocking
results in net domain growth, while $AaBb$ clocking results in domain
reversal.

Figure 6: Results of growth and reversal with bipolar clocking, and control
experiment. a, total magnetization of the ensembles subject to the different
bipolar clock protocols. The timeline indicates clock time, labeled by the
clock protocol. During the first phase, $AaBb$ clocking, the ensembles undergo
domain growth and increase in magnetization. The controls, $aA$ clocking and
$bB$ clocking, show no net growth. Further growth ($aAbB$ clocking) and
reversal ($AaBb$ clocking) occur after the controls. b, magnetic image
snapshots of the experimental ensemble, and zoomed in views of the flatspin
simulated ensemble, at the specified points in time. The growing domains
change morphology during the clock protocol. All XMCD-PEEM images are
available in LABEL:fig:all-DACB-snapshots. Videos of the experiment and
simulation are provided in LABEL:vid:bipolar-peem and LABEL:vid:bipolar-sim.

In contrast to unipolar clocking, bipolar clocking can also induce
morphological changes to the growing domains. As a result of the bipolarity of
the clock pulses, domains are now able to both grow and shrink within the same
clock cycle. In the experiment snapshots of Fig. 6b, we observe growth from
(1) to (2), followed by a clear change in domain morphology from (2) to (3),
and further growth between (3) to (4). In simulations, we can observe the
step-wise details of simultaneous growth and morphology changes, as shown in
the zoomed in snapshots. Inverting the clock protocol ($AaBb$ clocking)
results in domain reversal.

The deciding factor for growth or reversal is the polarity of the last clock
pulse at the transition between the two clocks. Each clock in $aAbB$ clocking,
for example, ends on the positive polarity at the transition ($aA$ and $bB$),
resulting in growth of the rightwards (orange/pink) domains.

Within a bipolar clock cycle, there is an apparent competition between growth
and reversal. Some domain wall configurations result in net domain growth
(others in net reversal), in a “one step back, two steps forward” process (see
section S2). In this way, a domain may grow horizontally and reverse
vertically, thereby gradually changing shape over time (see Fig. 12). While
the balance between growth and reversal can be delicate, there is a clear
trend for the clock protocols explored here, namely growth for $aAbB$, and
reversal for $AaBb$.

Compared to unipolar clocking, the dynamics in bipolar clocked pinwheel ASIs
are more varied and complex. While there is a gradual net domain growth, the
activity can intermittently spike and linger, depending on the particular
state of the ensemble (see LABEL:vid:bipolar-peem and LABEL:vid:bipolar-sim).
Bipolar clocking hence unlocks a wide variety of complex dynamic behavior in
pinwheel ASI, while at the same time offering considerable control by choice
of clock protocol.

##  7 Conclusions

We have introduced astroid clocking, a scheme for field-driven evolution in
nanomagnetic metamaterials. The method exploits the shape and orientation of
the nanomagnet switching astroids and dipolar coupling to selectively address
subsets of the nanomagnets. Pulsing specific fields in sequence results in
clocked dynamics that are both gradual and discrete in time. Considerable
control of the dynamics is available through choice of clock protocol.

This work demonstrates how astroid clocking can be used to control the growth
and reversal of ferromagnetic domains in pinwheel ASI. In this system,
unipolar clocking results in monotonic domain growth or reversal, while
bipolar clocking adds more complex dynamics that include changes to domain
morphology.

The principles of astroid clocking are not limited to pinwheel ASI, and are
applicable to a range of coupled nanomagnetic systems. Exploring the clocked
dynamics of established and future nanomagnetic metamaterials is an exciting
research direction. The space of possible clock protocols remains vast.

Astroid clocking offers unprecedented control and understanding of ASI
dynamics in both time and space. The method enables new directions in ASI
research and paves the way for novel devices based on nanomagnetic
metamaterials.

##  8 Extended data figures

Figure 7: Unipolar $AB$ clocking of an orange/pink (rightwards) domain as it
reaches the edge of the array. There is an apparent transition from horizontal
to vertical domain growth (5-6). Vertical growth proceeds by avalanches of
spin flips, starting at the bottom-left and top-right corners of the domain at
the array edge. Magnets that change state between snapshots are highlighted by
a solid black outline.  Figure 8: a-b, astroid clusters during reversal, when
the pinwheel system shown in c is subject to the negative clock field
$\mathbf{H}_{b}$. b, astroid clusters during reversal have a different
structure compared to growth. Switchable magnets outside the astroid are
highlighted in c. During reversal, the switchable magnets are along both the
horizontal, vertical and $-45\degree$ domain walls. Switchable magnets along
the horizontal domain wall is attributed to the curvature of the inner domain.
Figure 9: XMCD-PEEM images of all steps from the relevant unipolar clock
protocol series. Time starts at $t=0$, and is incremented by 1 for each clock
step, with clock pulses indicated by the labels. The black (rightwards)
domains grow with application of $AB$ clocking, and quickly reverses with $ab$
clocking. Red circle highlights: The short, vertical domain wall terminating
the black domain in the center region of snapshot $t=20$ exemplifies both
avalanching domain growth and a stuck domain wall. In snapshot $t=21$ the top
part of the domain wall has progressed in an avalanche to form a finger
extension of the domain, while the bottom part of the domain wall remains as
before.  Figure 10: XMCD-PEEM images of the control experiment. The system is
reinitialized at $t=51$ (following from Fig. 9), and $t$ is incremented by 1
for each clock step, with clock pulses indicated by the labels. The first $A$
clock pulse promotes dark (rightwards) magnets, equivalent to half a clock
cycle, while subsequent applications of $A$ incurs no further change. When the
clock pulse is changed to $B$, dark (rightwards) magnets are again promoted,
equivalent to the second half of an $AB$ clock cycle. Furthermore, additional
$B$ clock pulses incurs no change in the state.  Figure 11: Clock diagram of
bipolar $aAbB$ clocking followed by its inverse, $AaBb$ clocking. Bipolar
clocking employs both positive and negative clock pulses.  Figure 12: Bipolar
$aAbB$ clocking of pinwheel ASI. Each snapshot shows a zoomed-in view of a
$50\times 50$ system, at different points during a clock protocol. (1) shows
the initial state, an orange/pink (rightwards) domain in the center of an
otherwise polarized blue/green (leftwards) array. (2-12) show the state during
$aAbB$ clocking, with simultaneous domain growth (horizontally) and reversal
(vertically). As a result the domain gradually changes morphology over time.
Magnets that change state between snapshots are highlighted by a solid black
outline.

## Methods

### Sample fabrication details

The samples are arrays of permalloy nanomagnets fabricated in pinwheel ASI
geometries on a silicon substrate. The resist mixture, 1:2 CSAR 62:anisole, is
spin-coated onto the substrate at
$4000\text{\,}\mathrm{r}\mathrm{p}\mathrm{m}$, achieving a thickness of $\sim
100\text{\,}\mathrm{nm}$. Following coating, samples are soft baked at
$150\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for 1 minute. The desired
patterns, arrays of 220 x 80 stadium shaped nanomagnets in $30\times 30$ and
$50\times 50$ pinwheel geometries, are then exposed using the Elionix ELS-G100
EBL system. Samples are post-exposure baked at
$150\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for 1 minute. The patterned resist
is developed using AR600-546 for 1 minute, rinsed with isopropanol, and
nitrogen dried. Permalloy (Ni_0.79Fe_0.21) is deposited to a thickness of
$25\text{\,}\mathrm{nm}$ via electron beam evaporation using a Pfeiffer Vacuum
Classic 500 system, and capped with a $2\text{\,}\mathrm{nm}$ aluminium layer.
Finally, the samples undergo ultrasound-assisted lift-off using a dedicated
stripper (AR600-71), leaving behind the patterned permalloy nanomagnets. Post-
fabrication, the precision and quality of the fabricated nanomagnet arrays are
inspected using Scanning Electron Microscopy (SEM). This SEM inspection
confirmed that the permalloy nanomagnets are properly formed, free-standing,
and without significant defects.

### XMCD-PEEM and clocking procedure

Experimentally realized clocking of fabricated ASIs is carried out under
magnetic microscopy inspection. We use a photoemission electron microscope
with x-ray magnetic circular dichroism (XMCD-PEEM) for magnetic contrast to
observe single magnet states of the ASI ensembles[35]. An in-plane, bi-axial
quadrupole magnet with two pairs of coils and a split 2D-yoke provides astroid
clocking fields[36]. The signal at the Fe L_3 edge is exploited for
ferromagnetic XMCD contrast.

The orientation of the ASI ensembles, applied magnetic fields, and XMCD
contrast is carefully selected. Samples are mounted with top and bottom
ensemble edges parallel to the synchrotron light, with each nanomagnetic
element oriented $\pm 45\text{\,}\mathrm{\SIUnitSymbolDegree}$ to the light.
This orientation guarantees balanced magnetic contrast for nanomagnets of both
sublattices $L_{a}$ and $L_{b}$. The in-plane field direction is given
relative to the incoming x-ray illumination, with angle values increasing
counter-clockwise. Consequently, the field directions and ensemble orientation
align with the illustration in Fig. 1, with an added light axis (providing
magnetic contrast) parallel to the $h_{x}$-axis.

The general experimental procedure is to initialize the ASI system, then apply
clock protocols interspersed with magnetic imaging. We initialize the system
by applying a strong, polarizing magnetic field ($72\text{\,}\mathrm{mT}$
along $180\text{\,}\mathrm{\SIUnitSymbolDegree}$), followed by two smaller
fields, ($18\text{\,}\mathrm{mT}$ along
$0\text{\,}\mathrm{\SIUnitSymbolDegree}$ and $3.5\text{\,}\mathrm{mT}$ along
$180\text{\,}\mathrm{\SIUnitSymbolDegree}$) to demagnetize the yoke. For the
bipolar clocking, however, the initial field strength is
$82\text{\,}\mathrm{mT}$. The difference in field strength is due to observed
differences in the ensemble coercivity. Successful initialization is confirmed
by imaging a fully polarized ensemble (fully bright contrast (leftwards), as
in snapshot $t=0$ of Fig. 9) and the absence of remaining image translation in
the PEEM (indicating a demagnetized yoke).

After initialization, we perform steps of the clock protocols by alternating
the application of clock pulses $A$, $B$, $a$ or $b$. Each step of a clock
protocol comprises at least one clock pulse (ramping the applied field to
$\mathbf{H}_{i}$, holding the max field value, ramping down to zero applied
field), and a magnetic contrast image acquisition. The value of $H$ that
defines the $\mathbf{H}_{i}$ magnitudes is $62\text{\,}\mathrm{mT}$ for the
unipolar clocking, and $75\text{\,}\mathrm{mT}$ for the bipolar clocking.
After applying the first cycle of a clock protocol, before imaging, we shift
the image, using the electron microscope optics, to re-center the ensemble,
compensating for a small remanent magnetization in the yoke. We carry out
multiple cycles, each consisting of applying clock pulses and capturing an
image, while maintaining the same image shift throughout.

In addition to the growth and reversal protocols, we conduct a control
experiment by applying repeated clock pulses of $A$ and $B$ separately.

### flatspin simulations

Numerical simulations were done using flatspin, a large-scale ASI simulator
[33]. flatspin approximates each nanomagnet as a point dipole with position
$\mathbf{r}_{i}$ and orientation $\theta_{i}$. Each dipole then has two
possible magnetization directions along $\theta_{i}$, i.e., a binary macrospin
$s_{i}\in\\{-1,+1\\}$.

Each spin $i$ is influenced by a total field
$\mathbf{h}_{i}=\mathbf{h}_{\text{dip}}^{(i)}+\mathbf{h}_{\text{ext}}^{(i)}+\mathbf{h}_{\text{th}}^{(i)},$
where $\mathbf{h}_{\text{dip}}^{(i)}$ is the total dipolar field from
neighboring magnets, $\mathbf{h}_{\text{ext}}^{(i)}$ is a global or local
external field, and $\mathbf{h}_{\text{th}}^{(i)}$ is a stochastic magnetic
field representing thermal fluctuations in each magnetic element. The total
dipolar field is given by the magnetic dipole-dipole interaction,

$\mathbf{h}_{\text{dip}}^{(i)}=\alpha\sum_{j\neq
i}\frac{3\mathbf{r}_{ij}(\mathbf{m}_{j}\cdot\mathbf{r}_{ij})}{\lvert\mathbf{r}_{ij}\rvert^{5}}-\frac{\mathbf{m}_{j}}{\lvert\mathbf{r}_{ij}\rvert^{3}},$
(1)

where $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$ is the distance vector
from spin $i$ to $j$, and $\alpha$ scales the dipolar coupling strength
between spins. The coupling strength $\alpha$ is given by
$\alpha=\frac{\mu_{0}M}{4\pi a^{3}}$, where $a$ is the lattice spacing, $M$ is
the net magnetic moment of a single magnet, and $\mu_{0}$ is the vacuum
permeability.

Nanomagnet switching (magnetization reversal) occurs if the total field is
directed against the current magnetization $\mathbf{m}_{i}$ and the magnitude
of the field exceeds the coercive field $h_{\text{c}}$. flatspin employs a
generalized Stoner-Wohlfarth model, where $h_{\text{c}}$ depends on the angle
of the total field $\mathbf{h}_{i}$ with respect to the magnet orientation.
Associated with each magnet is a switching astroid, which describes
$h_{\text{c}}$ in terms of the parallel (easy axis) and perpendicular (hard
axis) component of the total field, $\mathbf{h}_{\parallel}$ and
$\mathbf{h}_{\perp}$. The shape of the switching astroid is described by the
equation

$\left(\frac{h_{\parallel}}{bh_{k}}\right)^{2/\gamma}+\left(\frac{h_{\perp}}{ch_{k}}\right)^{2/\beta}=1,$
(2)

where $h_{k}$ denotes the coercive field along the hard axis. The parameters
$b$, $c$, $\beta$, and $\gamma$ adjust the shape of the astroid: $b$ and $c$
define the height and width, respectively, while $\beta$ and $\gamma$ adjust
the curvature of the astroid at the easy and hard axis, respectively. Astroid
parameters are typically tuned to obtain a shape that agrees with results from
micromagnetic simulations.

Fabrication imperfections are modelled as variation in the coercive fields
$h_{k}^{(i)}$, which are sampled from a normal distribution
$\mathcal{N}(h_{k},\sigma)$, where $\sigma=k_{\text{disorder}}\cdot h_{k}$ and
$k_{\text{disorder}}$ is a user-defined parameter.

Dynamics are modeled using a deterministic single spin flip strategy. At each
simulation step, the total magnetic field $\mathbf{h}_{i}$ is calculated.
Next, we obtain a list of spins that may flip, according to the switching
astroid. Finally, the spin which is _furthest outside its switching astroid_
is flipped. The dipolar fields are recalculated after every spin flip, and the
above process is repeated until there are no more flippable spins. This
relaxation process is performed with constant external and thermal fields.

In this work, a global external field is used
($\mathbf{h}_{\text{ext}}^{(i)}=\mathbf{h}_{\text{ext}}$), and thermal
fluctuations are assumed to be negligible ($\mathbf{h}_{\text{th}}^{(i)}=0$).

The coupling strength $\alpha=0.0013$ was estimated to match the experimental
results from the $50\times 50$ fabricated pinwheel sample (see Sample
fabrication details). The value of $\alpha=0.0013$ is lower than predicted by
theory ($\alpha\approx 0.0025$), which is likely due to demagnetizing
oxidation of the permalloy. A partially oxidized nanomagnet will have a
reduced magnetic moment and a smaller effective size as the surface layer is
no longer ferromagnetic. The smaller $30\times 30$ sample used in Fig. 6 had a
slightly larger magnet spacing and $\alpha=0.0012$ was used in this case.

For the simulation studies, a field strength $H=$76.5\text{\,}\mathrm{mT}$$
and no disorder was used. Simulations accompanying the experimental results
used a slightly lower field strength of $H=$75.8\text{\,}\mathrm{mT}$$ for
Fig. 5 and $H=$75.9\text{\,}\mathrm{mT}$$ for Fig. 6.

Switching parameters were estimated from micromagnetic simulations of a
220x80x25 stadium magnet using mumax[37], namely
$h_{k}=$0.2\text{\,}\mathrm{T}$$, $b=0.38$, $c=1$, $\beta=1.3$, and
$\gamma=3.6$ Other parameters include $k_{\text{disorder}}=4\%$ and a neighbor
distance of $10$.

## Supplementary information

###  S1 Neighborhood interactions

Here we analyse what type of neighbor interactions causes switching to occur
selectively along the vertical and $+45\degree$ domain walls. We consider five
different prototype cases shown in LABEL:fig:astroid-distancec: a uniform
blue/green (leftwards) domain, and two domains separated by horizontal,
vertical, and $\pm 45\degree$ domain walls (DWs). Within each prototype case,
the subject of study is the highlighted blue magnet in the center. The circled
insets in the figure show only a limited neighborhood in the center of a
larger $50\times 50$ system which is initialized according to each prototype
case.

Figure 1: supfigure]fig:astroid-distanceMinimum distance to the astroid edge
as the neighborhood is increased, for the highlighted blue magnet in the
center of each scenario in c. In all cases, the clock field $\mathbf{H}_{A}$
is applied. a, distance to the astroid for the highlighted magnet, as the
neighborhood is increased when calculating the dipolar fields. b, shows a
trace of the position within the astroid as the neighborhood is increased.
Note that the scenarios all start at the same point (no neighbors), then
diverge.

LABEL:fig:astroid-distancea plots the distance to the astroid for the center
magnet, as the number of neighbors are increased when calculating the dipolar
fields. In other words, we compute the total dipolar field from all magnets
within a radius of the Nth nearest neighbor (NN). After adding the total
dipolar field to the external clock field $H_{A}$, the shortest distance to
the astroid is calculated. We define astroid distance as positive outside the
astroid and negative inside.

Astroid distance is plotted for each of the five prototype cases in
LABEL:fig:astroid-distancec. With zero neighbors, and hence no dipolar fields,
all five cases start at the same point outside the astroid. As the first NNs
are included, the cases split into four: the uniform domain and the
$-45\degree$ domain wall enter the astroid. In other words, the dipolar fields
from the first NNs stabilize and prevent switching in these two cases.
Including also the second NNs causes the horizontal domain wall to enter the
astroid. Horizontal domain walls are hence stabilized by 2nd NN interactions.
For the horizontal and $-45\degree$ domain walls, astroid distance does not
change significantly as the neighborhood is increased further. For the uniform
domain, however, astroid distance increases further as the number of NNs are
increased, with significant stabilizing interactions also beyond 9NNs.

Next, we consider the two cases where switching _does_ occur, namely the
vertical and $+45\degree$ domain walls. Somewhat curious, the astroid distance
for the vertical domain wall appears to stay nearly constant across all NNs.
The $+45\degree$ domain walls travel further outside the astroid due to 1st NN
interactions, then the 3rd NN interactions bring it closer to the astroid
again, after which it remains at a near-constant distance.

LABEL:fig:astroid-distanceb shows a trace of the location within the astroid
as the NNs are increased. For the vertical domain wall (red line), there is
indeed movement due to dipolar interactions, but the movement is exclusively
_parallel_ to the astroid edge. Hence, the astroid distance in this case
remains constant. For the $+45\degree$ domain wall (purple line), the movement
is purely in the perpendicular ($h_{\perp}$) direction for the 1st NN
interactions, then purely parallel ($h_{\parallel}$) from the 3rd NN fields.

An even more detailed picture is provided in LABEL:fig:astroid-distance-
neighbors, where each neighbor magnet is colored according to the contribution
of its dipolar field. Specifically, a magnet is colored red (blue) if its
dipolar field pushes the center magnet further out of (into) the astroid. The
shade of red (blue) represents how much the dipolar field contributes to
promote (prevent) switching of the center magnet. A magnet is colored white if
its dipolar field has no contribution on the resulting astroid distance.

Figure 2: supfigure]fig:astroid-distance-neighborsNeighborhood influence with
respect to center magnet. Each scenario depicts the magnetization state (top)
and the corresponding influence of the neighbors (bottom). Stronger red
signifies that the magnet is biasing the center magnet towards switching, and
stronger blue signifies that the magnet is biasing the center magnet away from
switching.

As can be seen in LABEL:fig:astroid-distance-neighbors, the neighborhood in
the uniform domain is dominated by magnets that prevent switching (colored
blue), with the highest contribution from the first NNs along the hard axis of
the center magnet. The same subset of the NNs are also the primary stabilizing
force of the $-45\degree$ DW. For the horizontal DW case, the dipolar fields
from the first NNs cancel out, and it is the second NNs that prevent
switching.

For the vertial DW, there is an apparent symmetry between neighbors that
prevent and promote switching. As a result the vertical DW is not stabilized
and hence easily switched. We saw earlier how this is because the dipolar
fields are directed parallel to the astroid edge. The $+45\degree$ DW is the
least stable, where 3/4 of the first NNs promote switching (colored red).

###  S2 Growth and reversal in bipolar clocking

During bipolar clocking, domain growth and reversal in a single clock cycle
can be observed for several domain wall configurations. LABEL:fig:aAbB-
clocking-dw shows the time evolution of different types of domain walls,
subject to $aAbB$ clocking. A straightforward example of simultaneous growth
and reversal can be seen in LABEL:fig:aAbB-clocking-dwd, which shows a
$+45\degree$ domain wall. Notice that the first clock pulse $a$ moves the
domain wall one step towards the left, and hence a reversal of the orange/pink
domain. However, the subsequent $A$ pulse immediately undoes this change _and_
moves the domain wall another step towards the right, advancing the domain
wall a total of two layers of the sublattice $L_{a}$ (orange magnets). Next,
the $b$ pulse has no effect, since the pink magnets along the domain wall are
stabilized by the dipolar fields from their neighbors. Finally, the $B$ pulse
moves the domain another step towards the right, flipping the next layer of
magnets from sublattice $L_{b}$ (from green to pink). As can be seen, the
result is an apparent growth of the orange/pink domain by a single layer along
the domain wall. The other domain wall cases in LABEL:fig:aAbB-clocking-dw
also show simultaneous growth and reversal, but are not discussed in further
detail.

Figure 3: supfigure]fig:aAbB-clocking-dwBipolar $aAbB$ clocking of four types
of domain walls in pinwheel ASI: a, horizontal DW, b, vertical DW, c,
$-45\degree$ DW and d, $+45\degree$ DW. Each domain wall is initialized to
fill the whole $50\times 50$ system from edge to edge. Each snapshot shows a
zoomed-in view of the system, at different points during a clock protocol. (1)
shows the initial state. (2-12) show the state during $aAbB$ clocking. Magnets
that change state between snapshots are highlighted by a solid black outline.

There is an apparent competition between growth and reversal. For the
$+45\degree$ domain wall discussed earlier, the competition seems to favor
growth. However, the situation strongly depends on the particular shape of the
domain. Fig. 12 (main text) shows the time evolution of a hexagonal domain
subject to $aAbB$ clocking. As can be seen, the domain both grows horizontally
and reverses vertically, and hence gradually changes shape over time. Since
vertical domain reversal depends on the curvature of the domain, the process
will stop when the domain grows too wide. The domain will continue to grow
horizontally, as horizontal domain growth is not dependent on curvature. As a
result, domain growth seems to out-compete reversal in this case. The end
result is an apparent tendency towards horizontally elongated domains.

Figure 4: supfigure]fig:all-DACB-snapshotsXMCD-PEEM images of all steps from
the bipolar clock protocol series. The time starts at $t=0$, and is
incremented by 1 for each image, with clock pulses indicated by the labels.
The black (rightwards) domains grow and change shape as the $aAbB$ protocol is
applied. There are two control series where $aA$ and $bB$ are applied, where
no change occurs. Note that there is missing data for $t=33$, but the ensemble
was still subjected to the clock pulses. At $t=39$ we image after each single
clock pulse. From $t=42$ the reverse protocol $BbAa$ is applied, and the black
(rightwards) domains shrink. Note that during the reversal protocol we still
image after a final $A$ or $B$ pulse, in order to keep a constant image shift.

###  S3 Videos

Video 1: supvideo]vid:unipolar-peemVideo showing all XMCD-PEEM images from the
unipolar clock protocol series. Each frame is labeled by the clock pulse(s)
preceding it. The series consists of growth, reversal, re-initialization,
control experiment and a second phase of growth.  Video 2:
supvideo]vid:unipolar-simVideo showing flatspin simulation of the unipolar
clock protocol series. Each frame depicts the clock pulse preceding it. The
series consists of growth, reversal, re-initialization and control experiment.
Video 3: supvideo]vid:bipolar-peemVideo showing all XMCD-PEEM images from the
bipolar clock protocol series. Each frame is labeled by the clock pulse(s)
preceding it. The series consists of growth, control experiment and reversal.
Video 4: supvideo]vid:bipolar-simVideo showing flatspin simulation of the
bipolar clock protocol series. Each frame depicts the clock pulse preceding
it. The series consists of growth, control experiment and reversal.

## Acknowledgments

These experiments were performed at the CIRCE beamline at ALBA Synchrotron
with the collaboration of ALBA staff.

This work was funded in part by the Norwegian Research Council through the
IKTPLUSS project SOCRATES (Grant no. 270961) and the TEKNOKONVERGENS project
SPrINTER (Grant No. 331821), and in part by the EU FET-Open RIA project
SpinENGINE (Grant No. 861618).

The Research Council of Norway is acknowledged for the support to the
Norwegian Micro- and Nano-Fabrication Facility, NorFab, project number 295864.
Simulations were executed on the NTNU EPIC compute cluster[38].

MAN, MF and MWK acknowledge funding from MCIN through grant number
PID2021-122980OB-C54 and MWK also acknowledges support through Marie
Sklodowska-Curie grant agreement No. 754397 (DOC-FAM) from EU Horizon 2020.

## References

  *   * Ladak et al. [2010] Ladak, S., Read, D.E., Perkins, G.K., Cohen, L.F., Branford, W.R.: Direct observation of magnetic monopole defects in an artificial spin-ice system. Nature Physics 6(5), 359–363 (2010) https://doi.org/10.1038/nphys1628
  * Gliga et al. [2017] Gliga, S., Hrkac, G., Donnelly, C., Büchi, J., Kleibert, A., Cui, J., Farhan, A., Kirk, E., Chopdekar, R.V., Masaki, Y., Bingham, N.S., Scholl, A., Stamps, R.L., Heyderman, L.J.: Emergent dynamic chirality in a thermally driven artificial spin ratchet. Nature Materials 16(October) (2017) https://doi.org/10.1038/nmat5007
  * Sendetskyi et al. [2019] Sendetskyi, O., Scagnoli, V., Leo, N., Anghinolfi, L., Alberca, A., Lüning, J., Staub, U., Derlet, P.M., Heyderman, L.J.: Continuous magnetic phase transition in artificial square ice. Physical Review B 99(21), 214430 (2019) https://doi.org/10.1103/PhysRevB.99.214430
  * Heyderman and Stamps [2013] Heyderman, L.J., Stamps, R.L.: Artificial ferroic systems: Novel functionality from structure, interactions and dynamics. Journal of Physics: Condensed Matter 25(36), 363201 (2013) https://doi.org/10.1088/0953-8984/25/36/363201
  * Skjærvø et al. [2020] Skjærvø, S.H., Marrows, C.H., Stamps, R.L., Heyderman, L.J.: Advances in artificial spin ice. Nature Reviews Physics 2(1), 13–28 (2020) https://doi.org/10.1038/s42254-019-0118-3
  * Jensen et al. [2018] Jensen, J.H., Folven, E., Tufte, G.: Computation in artificial spin ice. In: ALIFE 2018: The 2018 Conference on Artificial Life, pp. 15–22. MIT Press, Tokyo, Japan (2018). https://doi.org/10.1162/isal_a_00011
  * Arava et al. [2018] Arava, H., Leo, N., Schildknecht, D., Cui, J., Vijayakumar, J., Derlet, P., Kleibert, A., Heyderman, L.: Engineering Relaxation Pathways in Building Blocks of Artificial Spin Ice for Computation, 1–6 (2018) arxiv:1812.06936
  * Gypens et al. [2018] Gypens, P., Leliaert, J., Van Waeyenberge, B.: Balanced Magnetic Logic Gates in a Kagome Spin Ice. Physical Review Applied 9(3), 034004 (2018) https://doi.org/10.1103/PhysRevApplied.9.034004
  * Jensen and Tufte [2020] Jensen, J.H., Tufte, G.: Reservoir Computing in Artificial Spin Ice. In: ALIFE 2020: The 2020 Conference on Artificial Life, pp. 376–383. MIT Press, Online (2020). https://doi.org/10.1162/isal_a_00268
  * Hon et al. [2021] Hon, K., Kuwabiraki, Y., Goto, M., Nakatani, R., Suzuki, Y., Nomura, H.: Numerical simulation of artificial spin ice for reservoir computing. Applied Physics Express 14(3), 033001 (2021) https://doi.org/10.35848/1882-0786/abdcd8
  * Kaffash et al. [2021] Kaffash, M.T., Lendinez, S., Jungfleisch, M.B.: Nanomagnonics with artificial spin ice. Physics Letters A 402, 127364 (2021) https://doi.org/10.1016/j.physleta.2021.127364
  * Gartside et al. [2022] Gartside, J.C., Stenning, K.D., Vanstone, A., Holder, H.H., Arroo, D.M., Dion, T., Caravelli, F., Kurebayashi, H., Branford, W.R.: Reconfigurable training and reservoir computing in an artificial spin-vortex ice via spin-wave fingerprinting. Nature Nanotechnology 17(5), 460–469 (2022) https://doi.org/10.1038/s41565-022-01091-7
  * Schumann et al. [2010] Schumann, A., Sothmann, B., Szary, P., Zabel, H.: Charge ordering of magnetic dipoles in artificial honeycomb patterns. Applied Physics Letters 97(2), 022509 (2010) https://doi.org/%****␣main.bbl␣Line␣275␣****10.1063/1.3463482
  * Mengotti et al. [2011] Mengotti, E., Heyderman, L.J., Rodríguez, A.F., Nolting, F., Hügli, R.V., Braun, H.B.: Real-space observation of emergent magnetic monopoles and associated Dirac strings in artificial kagome spin ice. Nature Physics 7(1), 68–74 (2011) https://doi.org/10.1038/nphys1794
  * Morgan et al. [2011] Morgan, J.P., Stein, A., Langridge, S., Marrows, C.H.: Magnetic reversal of an artificial square ice: Dipolar correlation and charge ordering. New Journal of Physics 13(10), 105002 (2011) https://doi.org/10.1088/1367-2630/13/10/105002
  * Pollard et al. [2012] Pollard, S.D., Volkov, V., Zhu, Y.: Propagation of magnetic charge monopoles and Dirac flux strings in an artificial spin-ice lattice. Physical Review B 85(18), 180402 (2012) https://doi.org/10.1103/PhysRevB.85.180402
  * León [2013] León, A.: Heavy and light monopoles in magnetic reversion in artificial spin ice. Current Applied Physics 13(9), 2014–2018 (2013) https://doi.org/10.1016/j.cap.2013.08.010
  * Gilbert et al. [2015] Gilbert, I., Chern, G.-W., Fore, B., Lao, Y., Zhang, S., Nisoli, C., Schiffer, P.: Direct visualization of memory effects in artificial spin ice. Physical Review B 92(10), 104417 (2015) https://doi.org/10.1103/PhysRevB.92.104417 arxiv:1508.06330
  * Li et al. [2019] Li, Y., Paterson, G.W., Macauley, G.M., Nascimento, F.S., Ferguson, C., Morley, S.A., Rosamond, M.C., Linfield, E.H., MacLaren, D.A., Macêdo, R., Marrows, C.H., McVitie, S., Stamps, R.L.: Superferromagnetism and Domain-Wall Topologies in Artificial “Pinwheel” Spin Ice. ACS Nano 13(2), 2213–2222 (2019) https://doi.org/10.1021/acsnano.8b08884
  * Bingham et al. [2022] Bingham, N.S., Zhang, X., Ramberger, J., Heinonen, O., Leighton, C., Schiffer, P.: Collective Ferromagnetism of Artificial Square Spin Ice. Physical Review Letters 129(6), 067201 (2022) https://doi.org/10.1103/PhysRevLett.129.067201 . Accessed 2023-06-11
  * Begum Popy et al. [2022] Begum Popy, R., Frank, J., Stamps, R.L.: Magnetic field driven dynamics in twisted bilayer artificial spin ice at superlattice angles. Journal of Applied Physics 132(13), 133902 (2022) https://doi.org/10.1063/5.0118078 . Accessed 2023-06-11
  * Wang et al. [2006] Wang, R.F., Nisoli, C., Freitas, R.S., Li, J., McConville, W., Cooley, B.J., Lund, M.S., Samarth, N., Leighton, C., Crespi, V.H., Schiffer, P.: Artificial ’spin ice’ in a geometrically frustrated lattice of nanoscale ferromagnetic islands. Nature 439(7074), 303–306 (2006) https://doi.org/10.1038/nature04447 arxiv:cond-mat/0601429
  * Wang et al. [2007] Wang, R.F., Li, J., McConville, W., Nisoli, C., Ke, X., Freeland, J.W., Rose, V., Grimsditch, M., Lammert, P., Crespi, V.H., Schiffer, P.: Demagnetization protocols for frustrated interacting nanomagnet arrays. Journal of Applied Physics 101(9), 09–104 (2007) https://doi.org/10.1063/1.2712528
  * Nisoli et al. [2007] Nisoli, C., Wang, R., Li, J., McConville, W.F., Lammert, P.E., Schiffer, P., Crespi, V.H.: Ground state lost but degeneracy found: The effective thermodynamics of artificial spin ice. Physical Review Letters 98(21), 1–4 (2007) https://doi.org/10.1103/PhysRevLett.98.217203 arxiv:cond-mat/0702320
  * Ke et al. [2008] Ke, X., Li, J., Nisoli, C., Lammert, P.E., McConville, W., Wang, R.F., Crespi, V.H., Schiffer, P.: Energy Minimization and ac Demagnetization in a Nanomagnet Array. Physical Review Letters 101(3), 037205 (2008) https://doi.org/10.1103/PhysRevLett.101.037205
  * Qi et al. [2008] Qi, Y., Brintlinger, T., Cumings, J.: Direct observation of the ice rule in an artificial kagome spin ice. Physical Review B 77(9), 094418 (2008) https://doi.org/10.1103/PhysRevB.77.094418 arxiv:0802.0034
  * Nisoli et al. [2010] Nisoli, C., Li, J., Ke, X., Garand, D., Schiffer, P., Crespi, V.H.: Effective Temperature in an Interacting Vertex System: Theory and Experiment on Artificial Spin Ice. Physical Review Letters 105(4), 047205 (2010) https://doi.org/10.1103/PhysRevLett.105.047205
  * Phatak et al. [2011] Phatak, C., Petford-Long, A.K., Heinonen, O., Tanase, M., De Graef, M.: Nanoscale structure of the magnetic induction at monopole defects in artificial spin-ice lattices. Physical Review B 83(17), 174431 (2011) https://doi.org/10.1103/PhysRevB.83.174431
  * Perrin et al. [2016] Perrin, Y., Canals, B., Rougemaille, N.: Extensive degeneracy, Coulomb phase and magnetic monopoles in artificial square ice. Nature 540(7633), 410–413 (2016) https://doi.org/10.1038/nature20155 . Number: 7633 Publisher: Nature Publishing Group. Accessed 2023-06-11
  * Bingham et al. [2021] Bingham, N.S., Rooke, S., Park, J., Simon, A., Zhu, W., Zhang, X., Batley, J., Watts, J.D., Leighton, C., Dahmen, K.A., Schiffer, P.: Experimental Realization of the 1D Random Field Ising Model. Physical Review Letters 127(20), 207203 (2021) https://doi.org/10.1103/PhysRevLett.127.207203
  * Nomura et al. [2017] Nomura, H., Yoshioka, N., Miura, S., Nakatani, R.: Controlling operation timing and data flow direction between nanomagnet logic elements with spatially uniform clock fields. Applied Physics Express 10(12), 1–4 (2017) https://doi.org/%****␣main.bbl␣Line␣625␣****10.7567/APEX.10.123004
  * Macêdo et al. [2018] Macêdo, R., Macauley, G.M., Nascimento, F.S., Stamps, R.L.: Apparent ferromagnetism in the pinwheel artificial spin ice. Physical Review B 98(1), 014437 (2018) https://doi.org/10.1103/PhysRevB.98.014437
  * Jensen et al. [2022] Jensen, J.H., Strømberg, A., Lykkebø, O.R., Penty, A., Leliaert, J., Själander, M., Folven, E., Tufte, G.: Flatspin: A large-scale artificial spin ice simulator. Physical Review B 106(6), 064408 (2022) https://doi.org/10.1103/PhysRevB.106.064408
  * Tannous and Gieraltowski [2008] Tannous, C., Gieraltowski, J.: The Stoner–Wohlfarth model of ferromagnetism. European Journal of Physics 29(3), 475–487 (2008) https://doi.org/%****␣main.bbl␣Line␣675␣****10.1088/0143-0807/29/3/008 . Publisher: IOP Publishing. Accessed 2019-04-09
  * Aballe et al. [2015] Aballe, L., Foerster, M., Pellegrin, E., Nicolas, J., Ferrer, S.: The ALBA spectroscopic LEEM-PEEM experimental station: layout and performance. Journal of Synchrotron Radiation 22(3), 745–752 (2015) https://doi.org/10.1107/S1600577515003537 . Publisher: International Union of Crystallography. Accessed 2023-06-08
  * Foerster et al. [2016] Foerster, M., Prat, J., Massana, V., Gonzalez, N., Fontsere, A., Molas, B., Matilla, O., Pellegrin, E., Aballe, L.: Custom sample environments at the ALBA XPEEM. Ultramicroscopy 171, 63–69 (2016) https://doi.org/10.1016/j.ultramic.2016.08.016 . Accessed 2023-06-08
  * Vansteenkiste et al. [2014] Vansteenkiste, A., Leliaert, J., Dvornik, M., Helsen, M., Garcia-Sanchez, F., Van Waeyenberge, B.: The design and verification of MuMax3. AIP Advances 4(10), 107133 (2014) https://doi.org/10.1063/1.4899186 arxiv:1406.7635
  * Själander et al. [2019] Själander, M., Jahre, M., Tufte, G., Reissmann, N.: EPIC: An energy-efficient, high-performance GPGPU computing research infrastructure. arXiv:1912.05848 [cs.DC] (2019)


# JefiGPU: Jefimenko’s Equations on GPU

Jun-Jie Zhang  Jian-Nan Chen  Guo-Liang Peng  Tai-Jiao Du  Hai-Yan Xie
Beijing Institute of Technology, Beijing 100081, China  Northwest Institute of
Nuclear Technology, Xi’an 710024, China

###### Abstract

We have implemented a GPU version of the Jefimenko’s equations — JefiGPU.
Given the proper distributions of the source terms $\rho$ (charge density) and
$\mathbf{J}$ (current density) in the source volume, the algorithm gives the
electromagnetic fields in the observational region (not necessarily overlaps
the vicinity of the sources). To verify the accuracy of the GPU
implementation, we have compared the obtained results with that of the
theoretical ones. Our results show that the deviations of the GPU results from
the theoretical ones are around 5%. Meanwhile, we have also compared the
performance of the GPU implementation with a CPU version. The simulation
results indicate that the GPU code is significantly faster than the CPU
version. Finally, we have studied the parameter dependence of the execution
time and memory consumption on one NVIDIA Tesla V100 card. Our code can be
consistently coupled to RBG (Relativistic Boltzmann equations on GPUs) and
many other GPU-based algorithms in physics.

###### keywords:

Jefimenko’s equations; GPU; Heavy-Ion collisions; numba; ray;cupy.

††journal: Computer Physics Communications

PROGRAM SUMMARY/NEW VERSION PROGRAM SUMMARY

Program Title: JefiGPU
CPC Library link to program files: (to be added by Technical Editor)
Developer’s repository link: (if available)
Code Ocean capsule: (to be added by Technical Editor)
Licensing provisions(please choose one): Apache-2.0
Programming language: Python
Supplementary material:
Journal reference of previous version:*
Does the new version supersede the previous version?:*
Reasons for the new version:*
Summary of revisions:*
Nature of problem(approx. 50-250 words):Jefimenko’s equations are numerically
stable with given sources $\rho$ (charge density) and $\mathbf{J}$ (current
density). Providing proper sources, Jefimenko’s equations can give the
electromagnetic fields without boundary conditions. However, the relevant
integrations of the Jefimenko’s equations involve a 4D space-time which is
both memory and time consuming. With the help of the state-of-art GPU
technique, these integrations can be evaluated with acceptable accuracy and
execution time. In the present work, we have implemented a GPU version of the
Jefimenko’s equations and tested the accuracy and the performance of the code.
Our work have demonstrated a significant improvement of the performance
compared with a similar CPU implementation.
Solution method(approx. 50-250 words):The GPU code mainly deals with the
implementation of the integrations in Jefimenko’s equations. Each CUDA kernel
corresponds to a spatial point of $\mathbf{E}$ and $\mathbf{B}$. The
integrations involving the source region and the retarded time are performed
in each CUDA kernel. Ray, Numba and Cupy packages are used to manipulate the
scaling of GPU allocation and evaluation.
Additional comments including restrictions and unusual features (approx.
50-250 words):

##  1 Introduction

Electromagnetic force, as one of the four basic interactions, plays an
important role in almost all fields in nature. The governing equations of the
classical electromagnetic force are the well-known Maxwell’s equations. To
solve the Maxwell’s equations numerically, various algorithms including open-
source[1, 2, 3] and commercial software[4, 5, 6, 7] have been designed for
both general and specific usages. Among the many approaches to solve the
Maxwell’s equations, the Jefimenko’s equations[8, 9], known as the general
solutions of the Maxwell’s equations, are the ones that require only the
spatial distribution of the source terms $\rho$ (charge density) and
$\mathbf{J}$ (current density) to give the relevant electromagnetic (EM)
fields.

The benefits of the Jefimenko’s equations are two folds. First, they contain
two integrals of $\rho$ and $\mathbf{J}$. Combined with the continuum relation
$\nabla\cdot\mathbf{J}+\partial\rho/\partial t=0$, these integrals are
equivalent to the Maxwell’s equations. Integrals are much more stable than
differences. Compared with the FTDT (Finite-Difference Time-Domain)[10, 11, 1]
and FDFD (Finite-Difference Frequency- Domain)[12] methods, Jefimenko’s
equations do not have the trouble of numerical blow up (this is only true when
we use the concept of charge screening[13]). Second, the Jefimenko’s equations
do not need the boundary conditions — as long as one provides the proper
distributions of the sources. The $\mathbf{E}$ and $\mathbf{B}$ fields can be
obtained by performing an integration in the vicinity of the sources. This is
convenient especially when we calculate the EM fields in a fixed box, since in
FEM (Finite Element Method)[14, 15, 16], FDTD and FDFD one has to carefully
design the boundary conditions[17, 18, 19] to absorb the reflections of the EM
fields.

However, the Jefimenko’s equations also encounter certain difficulties. The
integral is performed over the entire source volume space with retarded time.
This requires all the information of the sources in the relevant space-time,
and is therefore memory demanding. Meanwhile, the denominator of the
integrated involves a divergence at $\mathbf{r}=\mathbf{r}^{\prime}$ (see Esq.
(1) ~ 3). The divergence is somewhat troublesome in numerical evaluations.
Thus, the refinement’s equations are usually used to solve far-field EM
fields[20, 21]. Except for the numerical difficulties, the source term $\rho$
may not be easily available since it is experimentally hard to measure the
charge density distributions[22]. To provide a distribution of the charge and
current densities, one may need to solve the transport equations (e.g., the
Boltzmann equation [23]) which is another challenge so far.

In this regard, we have applied the Jefimenko’s equation on the state-of-art
GPU based on our previous work ZMCintegral[24, 25] and RBG[23]. Our present
code mainly contributes to the field of Heavy-Ion collisions[26, 27, 28, 29]
where the electromagnetic fields in the early time significantly influence the
evolution of the quark-gluon matter. Our hope is to combine the Jefimenko’s
equations to the relativistic Boltzmann equation[23] and eventually answer the
question of “early thermalization”[30, 31, 32]. Due to the high clock rate,
high instruction per cycle[33], and multiple cores of the GPUs, more and more
numerical algorithms that were hard in the past are now possible[24, 34]. In
the present work, we use the Python package Ray[35], Numba[36] and Cupy[37] to
perform the jefimenko’s equations on GPU. Our work has demonstrated
significant performance improvement compared with the usual CPU
implementation. The results are stable and within an acceptable margin of
error.

The paper is organized in the following structure. In Sec. 2, we introduce the
expressions of the Jefimenko’s equations and the details of the GPU
implementation. In Sec. 3, we compare the GPU results with the theoretical
results. The corresponding execution time of the GPU code is also compared
with a similar CPU implementation. In Sec. 4, we have investigated the
parameter dependence of the GPU code on one Tesla V100 card, and presented the
time consumption of HToD (Host To Device), DToH (Device To Host) and the
evaluation. The conclusion and outlook of JefiGPU is in Sec. 5. Through out
the paper we use the Rationalized-Lorentz-Heaviside-QCD unit.

##  2 GPU Implementation

###  2.1 Jefimenko’s equations

The Jefimenko’s equation can be directly derived from the retarded potential
solution of Maxwell’s equations[22]. In Rationalized-Lorentz-Heaviside-QCD
unit, the Jefimenko’s equations take the following form

$\displaystyle\mathbf{E}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\int\left[\frac{\mathbf{r}-\mathbf{r}^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{3}}\rho(\mathbf{r}^{\prime},t_{r})\right.$
(1)
$\displaystyle+\frac{\mathbf{r}-\mathbf{r}^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{2}}\frac{\partial\rho(\mathbf{r}^{\prime},t_{r})}{\partial
t}$
$\displaystyle\left.-\frac{1}{|\mathbf{r}-\mathbf{r}^{\prime}|}\frac{\partial\mathbf{J}(\mathbf{r}^{\prime},t_{r})}{\partial
t}\right]d^{3}\mathbf{r}^{\prime}$ $\displaystyle\mathbf{B}(\mathbf{r},t)$
$\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\int\left[\frac{\mathbf{r}-\mathbf{r}^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{3}}\times\mathbf{J}(\mathbf{r}^{\prime},t_{r})\right.$
(2)
$\displaystyle+\left.\frac{\mathbf{r}-\mathbf{r}^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{2}}\times\frac{\partial\mathbf{J}(\mathbf{r}^{\prime},t_{r})}{\partial
t}\right]d^{3}\mathbf{r}^{\prime}$ $\displaystyle t_{r}$ $\displaystyle=$
$\displaystyle t-|\mathbf{r}-\mathbf{r}^{\prime}|,$ (3)

where $\mathbf{E}$ and $\mathbf{B}$ are the electric and magnetic field at
space-time point $(\mathbf{r},t)$, $\rho$ and $\mathbf{J}$ are the charge and
current densities, and $t_{r}$ is the retarded time. Once the space-time
distributions of $\rho$ and $\mathbf{J}$ are given, one can use Eqs. (1) ~ 3
to obtain the field distributions in space-time.

The numerical difficulties of the Jefimenko’s equations lie in the following
two aspects. First, to obtain the numerical value of $\mathbf{E}$ and
$\mathbf{B}$ at each point in the observational region (not necessarily
overlaps the source region), an integration involving a 4-Dimensional (4D)
space-time is required. Supposing the observational region has
$n_{x,\text{o}}\times n_{y,\text{o}}\times n_{z,\text{o}}=100^{3}\sim 10^{6}$
points, we need to perform $10^{6}$ integrations at each time step. If each
integration involves a summation of $n_{x,\text{s}}\times n_{y,\text{s}}\times
n_{z,\text{s}}=100^{3}\sim 10^{6}$ points, the task will be both time and
memory consuming. Second, Jefimenko’s equations encounter a divergence at
$\mathbf{r}=\mathbf{r}^{\prime}$ when the observational region overlaps the
source region. The near-source divergence has its profound physical origins
and the usual treatment is to adopt the concept of coarse screening[13], which
is constrained by the grid size. With coarse screening, the $\mathbf{E}$ and
$\mathbf{B}$ fields in a specific grid can only be generated by the charge and
density sources from other grids. If one is working in the far-field
calculations, then the divergence will not be a problem.

###  2.2 GPU implementation of Jefimenko’s equations

We define two spatial regions of sizes $n_{x,\text{s}}\times
n_{y,\text{s}}\times n_{z,\text{s}}$ and $n_{x,\text{o}}\times
n_{y,\text{o}}\times n_{z,\text{o}}$, where the subscripts o and s denote the
physical quantities in the source and observational regions, receptively. The
infinitesimal differences in these two regions are denoted as
$dx_{\text{s}}\times dy_{\text{s}}\times dz_{\text{s}}$ and
$dx_{\text{o}}\times dy_{\text{o}}\times dz_{\text{o}}$. In the current work,
we only use one NVIDIA V100 card to test the code, and the functionality of
multi-GPU on clusters is also supported via the Ray package.

For numerical convenience, we discretize Eqs. (1) ~ 3 as

$\displaystyle\mathbf{E}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\Sigma_{i,j,k}\left[\frac{\mathbf{r}-\mathbf{r}_{i,j,k}^{\prime}}{|\mathbf{r}-\mathbf{r}_{i,j,k}^{\prime}|^{3}}\rho(\mathbf{r}_{i,j,k}^{\prime},t_{r})\right.$
(4)
$\displaystyle+\frac{\mathbf{r}-\mathbf{r}_{i,j,k}^{\prime}}{|\mathbf{r}-\mathbf{r}_{i,j,k}^{\prime}|^{2}}\frac{\rho(\mathbf{r}_{i,j,k}^{\prime},t_{r})-\rho(\mathbf{r}_{i,j,k}^{\prime},t_{r}-dt)}{dt}$
$\displaystyle\left.-\frac{1}{|\mathbf{r}-\mathbf{r}_{i,j,k}^{\prime}|}\frac{\mathbf{J}(\mathbf{r}_{i,j,k}^{\prime},t_{r})-\mathbf{J}(\mathbf{r}_{i,j,k}^{\prime},t_{r}-dt)}{dt}\right]d\Omega^{\prime}$
$\displaystyle\mathbf{B}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\Sigma_{i,j,k}\left[\frac{\mathbf{r}-\mathbf{r}_{i,j,k}^{\prime}}{|\mathbf{r}-\mathbf{r}_{i,j,k}^{\prime}|^{3}}\times\mathbf{J}(\mathbf{r}_{i,j,k}^{\prime},t_{r})\right.$
(5)
$\displaystyle+\left.\frac{\mathbf{r}-\mathbf{r}_{i,j,k}^{\prime}}{|\mathbf{r}-\mathbf{r}_{i,j,k}^{\prime}|^{2}}\times\frac{\mathbf{J}(\mathbf{r}_{i,j,k}^{\prime},t_{r})-\mathbf{J}(\mathbf{r}_{i,j,k}^{\prime},t_{r}-dt)}{dt}\right]d\Omega^{\prime}$
$\displaystyle t_{r}$ $\displaystyle=$ $\displaystyle
t-|\mathbf{r}-\mathbf{r}_{i,j,k}^{\prime}|,$ (6)

where index $i\in\\{1,2,...,n_{x}\\}$, $j\in\\{1,2,...,n_{y}\\}$,
$k\in\\{1,2,...,n_{z}\\}$ and volume element
$d\Omega^{\prime}=dx^{\prime}dy^{\prime}dz^{\prime}$.

Figure 1: Schematic of the GPU implementation of Jefimenko’s equations. Each
CUDA kernel corresponds to the evaluations of $\mathbf{E}(\mathbf{r}_{i},t)$
and $\mathbf{B}(\mathbf{r}_{i},t)$ at space-time point $(\mathbf{r}_{i},t)$
where $i$ is the index of the spatial grid. The integration at point
$(\mathbf{r}_{i},t)$ is performed in each CUDA kernel via Eqs. (4) ~ (6).

The single GPU implementation is illustrated in Fig. 1. In each CUDA kernel we
perform an integration to obtain $\mathbf{E}(\mathbf{r}_{i},t)$ and
$\mathbf{B}(\mathbf{r}_{i},t)$ at space-time point $(\mathbf{r}_{i},t)$ where
$i$ is the index of the spatial grid. At each call of the integration method
(jargon in Python) we claim an empty space for $\mathbf{E}(\mathbf{r},t)$ and
$\mathbf{B}(\mathbf{r},t)$ in the global GPU memory. The calculated EM fields
via Eqs (4) ~ (6) will be stored in this empty space. After the evaluation of
the integral (i.e., summation of Eqs. (4) ~ (6)), new values of
$\mathbf{E}(\mathbf{r},t)$ and $\mathbf{B}(\mathbf{r},t)$ at time snapshot $t$
will be obtained (i.e., stored in the pre-allocated empty space) and
transferred to the host memory. Since $\mathbf{E}(\mathbf{r},t)$ and
$\mathbf{B}(\mathbf{r},t)$ are frequently transferred to the host at each
time, we store them in the global GPU memory for fast data transferring. The
numerical integration of Eqs. (4) ~ (6) in each CUDA kernel requires a full
access of $\rho(\mathbf{r}^{\prime},t_{r})$ and
$\mathbf{J}(\mathbf{r}^{\prime},t_{r})$ at all relevant space-time points (see
Sec. 4 for details), thus $\rho$ and $\mathbf{J}$ should be stored in the
global GPU memory. Meanwhile, the GPU memory corresponds to $\rho$ and
$\mathbf{J}$ will not be released after each call of the integration method.
Other temporary variables such as $\mathbf{r}^{\prime}$ and $t_{r}$ are
defined in local GPU memories.

The scaling of JefiGPU on GPU clusters is straight forward — the Ray package
provides a succinct API to manipulate multiple GPU cards and cluster nodes.
There are three approaches to conduct the GPU scaling, as is illustrated in
Fig. 2. We can divide the regions into several rectangular regions of S-O
(Source region and Observational region) pairs. Each pair can be performed on
one GPU card. In the code, we provide the users with a Python API to
manipulate these procedures.

Figure 2: Three approaches for GPU scaling. The source and observational
regions are divided into several sub-regions. The sub-regions will form many
S-O pairs. Each S-O pair can be performed on one GPU card using the algorithm
illustrated in Fig. 1.

##  3 Comparisons with the theoretical results and the CPU based algorithm

To see the performance of the GPU based algorithm, we will compare it with the
corresponding CPU implementation. The hardware set-up of the two cases are
summarized in Tab. 1. The CPU code is performed via Python with the similar
algorithm stated in Fig. 1. The theoretical results are obtained in
Mathematica with the built-in function NIntegrate using the GlobalAdaptive
method. To avoid near-source divergence in the numerical integrations of Eqs.
(1) ~ 3 in NIntegrate, we introduce a screening volume of size $\delta^{3}$.
Thus both $\rho$ and $\mathbf{J}$ give zero contribution to $\mathbf{E}$ and
$\mathbf{B}$ when $|i-i^{\prime}|\leq\delta$, where $i\in\\{x,y,z\\}$.

Table 1: Hardware set-up of the two conditions. | Hardware of CPU | Hardware of GPU
---|---|---
GPU | Intel(R) Xeon(R) Silver 4110 | NVIDIA Tesla V100
Implementation |<EMAIL_ADDRESS>GHz with 32 cores | with 32Gb memory
CPU | Intel(R) Xeon(R) Silver 4110 | None
Implementation |<EMAIL_ADDRESS>GHz with 1 core |

###  3.1 Constant sources

We will firstly use the constant sources to see whether the GPU-Code gives
stable and acceptable results, and then compare the execution time of the GPU-
Code to the CPU-Code. The sources of $\rho$ and $\mathbf{J}$ take the
following form

$\displaystyle\rho(x,y,z,t)$ $\displaystyle=$
$\displaystyle\begin{cases}0&\text{if }t<0\\\ 1&\text{if }t\geq 0\end{cases}$
$\displaystyle\mathbf{J}(x,y,z,t)$ $\displaystyle=$
$\displaystyle\begin{cases}\mathbf{0}&\text{if }t<0\\\ \mathbf{1}&\text{if
}t\geq 0\end{cases},$ (7)

where $\mathbf{0}=(0,0,0)$ and $\mathbf{1}=(1,1,1)$.

Figs. 3~5 show the results of the obtained electric and magnetic fields in the
XOY plane. The integration domain (also the domain of the observational
region) is of size $[-3\text{ GeV}^{-1},3\text{$\text{ GeV}^{-1}$}]^{3}$. The
theoretical results are obtained via NIntegrate in Mathematica with screening
length $\delta=0.2\text{ $\text{ GeV}^{-1}$}$, which is half the grid length
in the GPU-Code. In the GPU implementation we have chosen
$n_{x,\text{s}}=n_{y,\text{s}}=n_{z,\text{s}}=n_{x,\text{o}}=n_{y,\text{o}}=n_{z,\text{o}}=15$
such that each 3D grid is of volume size $(0.4\text{ $\text{ GeV}^{-1}$})^{3}$
.

|  |
---|---|---
|  |
Figure 3: Electric and magnetic field in the XOY plane with constant $\rho$
and $\mathbf{J}$. The upper panel gives the theoretical results via
Mathematica with built-in function NIntegrate, while the lower panel gives the
results obtained by the GPU code.

Fig. 3 compares $\mathbf{E}$ and $\mathbf{B}$ at time snapshots 2.5
$\text{GeV}^{-1}$, 5 $\text{GeV}^{-1}$ and 7.5 $\text{GeV}^{-1}$. The upper
panel shows the theoretical results with screening volume
$\delta^{3}=(0.2\text{ $\text{GeV}^{-1}$})^{3}$, and the lower panel gives the
GPU results. We can see that the two results follow the similar patterns at
all time snapshots. Since the sources are constant, we expect a saturation of
$\mathbf{E}$ and $\mathbf{B}$ at later times (e.g., when $t\geq 5\text{$\text{
GeV}^{-1}$}$ the fields remain almost unchanged). The near saturation at
around 5 $\text{GeV}^{-1}$ is reasonable since the longest transmission length
of the field is $6\text{ $\text{GeV}^{-1}$}\times\sqrt[3]{3}\sim
8.7\text{$\text{ GeV}^{-1}$}$,. Thus the system must reach a constant field
distribution after $8.7\text{$\text{ GeV}^{-1}$}$.

Fig. 4 and 5 give the comparison of
$|\mathbf{E}_{\text{Theo}}-\mathbf{E}_{\text{GPU}}|$ with
$|\mathbf{E}_{\text{Theo}}|$ and
$|\mathbf{B}_{\text{Theo}}-\mathbf{B}_{\text{GPU}}|$ with
$|\mathbf{B}_{\text{Theo}}|$. In the XOY plane, the norms are taken as

$\displaystyle|\mathbf{E}_{\text{Theo}}-\mathbf{E}_{\text{GPU}}|$
$\displaystyle=$
$\displaystyle\left[(E_{x,\text{Theo}}-E_{x,\text{GPU}})^{2}\right.$
$\displaystyle+\left.(E_{y,\text{Theo}}-E_{y,\text{GPU}})^{2}\right]^{1/2}$
$\displaystyle|\mathbf{E}_{\text{Theo}}|$ $\displaystyle=$
$\displaystyle\sqrt{E_{x,\text{Theo}}^{2}+E_{y,\text{Theo}}^{2}}$
$\displaystyle|\mathbf{B}_{\text{Theo}}-\mathbf{B}_{\text{GPU}}|$
$\displaystyle=$
$\displaystyle\left[(B_{x,\text{Theo}}-B_{x,\text{GPU}})^{2}\right.$
$\displaystyle+\left.(B_{y,\text{Theo}}-B_{y,\text{GPU}})^{2}\right]^{1/2}$
$\displaystyle|\mathbf{B}_{\text{Theo}}|$ $\displaystyle=$
$\displaystyle\sqrt{B_{x,\text{Theo}}^{2}+B_{y,\text{Theo}}^{2}},$ (8)

where we neglect the $z-$component of the fields for numerical convenience.
From the legend bars in Fig. 4 and 5, we can see that the maximum deviation of
the GPU-Code from the theoretical results is within $0.015/0.5\sim 3$ for both
$\mathbf{E}$ and $\mathbf{B}$.

|  |
---|---|---
|  |
Figure 4: Norm of $\mathbf{E}_{\text{Theo}}-\mathbf{E}_{\text{GPU}}$ and $\mathbf{E}_{\text{Theo}}$ with constant sources. The norm of the two vectors are defined in the XOY plane, i.e., $|\mathbf{A}|\equiv\sqrt{A_{x}^{2}+A_{y}^{2}}$. |  |
---|---|---
|  |
Figure 5: Norm of $\mathbf{B}_{\text{Theo}}-\mathbf{B}_{\text{GPU}}$ and
$\mathbf{B}_{\text{Theo}}$ with constant sources.

###  3.2 Sinusoidal sources

In the sinusoidal case, we choose two different domains for the observational
($[-3\text{ GeV}^{-1},3\text{$\text{ GeV}^{-1}$}]^{2}\times[10\text{
GeV}^{-1},16\text{ GeV}^{-1}]$, where $[10\text{ GeV}^{-1},16\text{
GeV}^{-1}]$ is the interval in $z$ direction) and source regions ($[-3\text{
GeV}^{-1},3\text{$\text{ GeV}^{-1}$}]^{3}$). The sources of sinusoidal $\rho$
and $\mathbf{J}$ take the following form

$\displaystyle\rho(x,y,z,t)$ $\displaystyle=$
$\displaystyle\begin{cases}0&\text{if }t<0\\\
3(\text{cos}(t)-1)\text{cos}(x+y+z)&\text{if }t\geq 0\end{cases}$
$\displaystyle\mathbf{J}(x,y,z,t)$ $\displaystyle=$
$\displaystyle\begin{cases}\mathbf{0}&\text{if }t<0\\\
(j_{x},j_{y},j_{z})&\text{if }t\geq 0\end{cases},$ (9)

where $j_{x}=j_{y}=j_{z}=\text{sin}(x+y+z)\text{sin}(t)$. Eq. (9) satisfies
the continuum relation $\nabla\cdot\mathbf{J}+\partial\rho/\partial t=0$.

Fig. 6 confirms that the patterns of the EM fields obtained from the
theoretical and the GPU code are similar. Different from the case of the
constant sources, there is no saturation of the EM fields. Instead, we find
the periodic features of the EM fields which are consistent with the
sinusoidal sources. Fig. 7 and 8 show a maximum deviation from the theoretical
results of $0.004/0.035\sim 11.2$ for $\mathbf{E}$ and $0.005/0.035\sim 14.3$
for $\mathbf{B}$, respectively.

|  |
---|---|---
|  |
Figure 6: Electric and magnetic field in the XOY plane with sinusoidal $\rho$ and $\mathbf{J}$.  |  |
---|---|---
|  |
Figure 7: Norm of $\mathbf{E}_{\text{Theo}}-\mathbf{E}_{\text{GPU}}$ and $\mathbf{E}_{\text{Theo}}$ with sinusoidal sources. |  |
---|---|---
|  |
Figure 8: Norm of $\mathbf{B}_{\text{Theo}}-\mathbf{B}_{\text{GPU}}$ and
$\mathbf{B}_{\text{Theo}}$ with sinusoidal sources.

###  3.3 Comparison of the execution time

For a comparison of the execution time, we have performed a similar algorithm
depicted in Fig. 1 on a CPU apparatus whose hardware set-up is listed in Tab.
1. We have chosen the same condition as in Sec. 3.1. Tab. 2 shows the
execution time of both the CPU and GPU implementations. One can see that the
use of GPU can significantly enhance the performance by a factor of up to a
thousand. The detailed time consumption of HToD (Host to Device) and DToH
(Device to Host) will be demonstrated in Sec. 4.

Table 2: Comparison of the Execution time in GPU and CPU implementations. To save time, the CPU code only evaluates the EM field for one time-step, while the GPU code runs for totally 410 time snapshots ranging from 0 $\text{GeV}^{-1}$ to 20.5 $\text{GeV}^{-1}$. The two time-info in the table correspond to the constant/sinusoidal sources (e.g., 8.696/8.111 s means that 8.696 s is the execution time for the constant source and 8.111 s for the sinusoidal source).The last row gives the effective execution time where both implementations are assumed to be executed for the same time snapshots. Apparatus | GPU | CPU
---|---|---
Execution time | 8.696/8.111 s | 18.012/18.540 s
Evaluated time snapshots | 410 | 1
Effective execution time | 8.696/8.111 s | 2.051/2.112 h

##  4 Parameter dependence of the algorithm on one Tesla V100 card

There are two parameters that mainly affect the performance of the algorithm,
i.e., “lts” short for length_of_time_snapshots and “ntg” short for
number_of_total_grids. In Eqs. (4) ~ (6), the summations of the source terms
$\rho$ and $\mathbf{J}$ are at the retarded time $t_{r}$. Therefore we need to
store $\rho$ and $\mathbf{J}$ in the GPU memory at all relevant times, i.e.,
$\rho_{\text{GPU}}\mapsto[\rho_{0}(\mathbf{r}),\rho_{1}(\mathbf{r}),\cdots\rho_{\text{lts}}(\mathbf{r})]$
and
$\mathbf{J}_{\text{GPU}}\mapsto[\mathbf{J}_{0}(\mathbf{r}),\mathbf{J}_{1}(\mathbf{r}),\cdots\mathbf{J}_{\text{lts}}(\mathbf{r})]$.
If the maximum distance between the source point $\mathbf{r}_{\text{s}}$ and
the observational point $\mathbf{r}_{\text{o}}$is
$d_{\text{m}}=\text{max}\\{|\mathbf{r}_{\text{s}}-\mathbf{r}_{\text{o}}|\\}$,
and the time step used in the calculation is $dt$, then we can determine the
value of length_of_time_snapshots by $\text{\textquotedblleft
lts\textquotedblright}=\left\lceil d_{\text{m}}/dt\right\rceil$.
number_of_total_grids counts the total grid sizes of both the observational
and source region, and is determined by the pair $\text{ \textquotedblleft
ntg\textquotedblright}=[n_{x,\text{s}}\times n_{y,\text{s}}\times
n_{z,\text{s}},n_{x,\text{o}}\times n_{y,\text{o}}\times n_{z,\text{o}}]$. The
total amount of float numbers on one GPU card can be evaluated approximately
as

$\displaystyle N$ $\displaystyle=$ $\displaystyle(\left\lceil
d_{\text{m}}/dt\right\rceil+1)\times n_{x,\text{s}}\times n_{y,\text{s}}\times
n_{z,\text{s}}\times 4$ (10) $\displaystyle+n_{x,\text{o}}\times
n_{y,\text{o}}\times n_{z,\text{o}}\times 6,$

where number “4” corresponds to $\rho,J_{x},J_{y},J_{z}$ and number “6”
corresponds to $E_{x},E_{y},E_{z},B_{x},B_{y},B_{z}$. We can use $N$ to
evaluate the occupied GPU memory in the execution. Apart from the memory
consumption, the execution time is mainly constrained by “ntg” since the
summation and loop are on the spatial grids.

Table 3: Occupied GPU memory at various “ntg” and “lts”. | “ntg” | $[160^{3},160^{3}]$ | $[100^{3},100^{3}]$ | $[70^{3},70^{3}]$ | $[60^{3},60^{3}]$ | $[20^{3},20^{3}]$
---|---|---|---|---|---|---
“lts” |  |  |  |  |  |
9000 |  |  |  |  | 30025MiB | 1461MiB
6000 |  |  |  | 31775MiB | 20137MiB | 1093MiB
2000 |  |  | 30899MiB | 10839MiB | 6953MiB | 605MiB
500 |  | 31707MiB | 8011MiB | 2991MiB | 2009MiB | 421MiB

In Tab. 3, we have given the GPU memory consumption with different values of
“ntg” and “lts”. The memory consumption is roughly consistent with Eq. 10. For
specific tasks, users can choose the proper parameter values according to the
GPU condition. Tab. 4 shows the time consumption of HToD, DToH and GPU
execution for one time-step. As a comparison with the CPU implementation, we
have also conducted a similar calculation with $\text{\textquotedblleft
ntg\textquotedblright}=[20^{3},20^{3}]$ and $\text{\textquotedblleft
lts\textquotedblright}=500$. The corresponding execution time for one time-
step in the CPU calculation is 102.25 s, which is around 100 times slower than
the GPU implementation. However, this does not mean that the GPU version is
only 100 faster than the CPU version. From Eqs. (4) ~ (6), we see that the
execution time increases exponentially with “ntg”. Therefore, one should
expect a significant performance improvement with larger “ntg” in the GPU
implementation. In our CPU version, we have only used one CPU core, therefore
the execution time should be extremely large for $\text{\textquotedblleft
ntg\textquotedblright}\sim[70^{3},70^{3}]$, and is impossible to perform.

Table 4: Time consumption with various values of “ntg” and “lts”. The upper
and lower panels give the time consumption of [DToH, HToD] and the GPU
execution at each time step.

| “lts” | 9000 | 6000 | 2000 | 500
---|---|---|---|---|---
“ntg” |  |  |  |  |
$[160^{3},160^{3}]$ |  |  |  |  | [0.060,0.042] s
$[100^{3},100^{3}]$ |  |  |  | [0.027,0.008] s | [0.027,0.009] s
$[70^{3},70^{3}]$ |  |  | [0.009,0.008] s | [0.009, 0.004] s | [0.009, 0.005] s
$[60^{3},60^{3}]$ |  | [0.005,0.001] s | [0.005,0.004] s | [0.006,0.003] s | [0.002,0.004] s
$[20^{3},20^{3}]$ |  | [0.0004,0.016] s | [0.0004,0.002] s | [0.0004,0.002] s | [0.0004,0.002] s

| “lts” | 9000 | 6000 | 2000 | 500
---|---|---|---|---|---
“ntg” |  |  |  |  |
$[160^{3},160^{3}]$ |  |  |  |  | 521.970 s
$[100^{3},100^{3}]$ |  |  |  | 32.583 s | 32.565 s
$[70^{3},70^{3}]$ |  |  | 4.845 s | 4.813 s | 4.858 s
$[60^{3},60^{3}]$ |  | 2.706 s | 2.668 s | 2.679 s | 2.720 s
$[20^{3},20^{3}]$ |  | 1.006 s | 1.008 s | 1.067 s | 1.001 s

##  5 Conclusion

In the current work, we have provided an implementation of the Jefimenko’s
equations on GPU based on our previous works[23, 24, 25]. Our code gives
stable and accurate results compared with the theoretical calculations. To see
the enhancement of the GPU code, we have also performed a similar algorithm in
with 1 CPU core available. A significant improvement of the execution time has
been observed in the GPU implementation. Finally, we give the parameter
dependence of the GPU code on one NVIDIA Tesla V100 card. The simulation
results in Tab. 4 shows that the execution time of the GPU implementation is
at least 100 times faster compared with single CPU.

In the API we have provided users with scaling functionalities. With specified
domains of both the observational and source regions, the code can be easily
used in multi-GPU clusters. However, the current version of the code only
supports multi-GPU manipulations via the Ray package. For an automatic use of
these scaling, we still require a heuristic module to allocate the GPUs.
Meanwhile, the current version does not support expanding systems like the
condition in heavy-ion collisions. We will consider all these factors in the
future.

## References

  * [1] C. Warren, A. Giannopoulos, A. Gray, I. Giannakis, A. Patterson, L. Wetter, A. Hamrah, A CUDA-based GPU engine for gprMax: Open source FDTD electromagnetic simulation software, Computer Physics Communications 237 (2019) 208–218.  doi:10.1016/j.cpc.2018.11.007.
  * [2] C. Warren, A. Giannopoulos, I. Giannakis, gprMax: Open source software to simulate electromagnetic wave propagation for ground penetrating radar, Computer Physics Communications 209 (2016) 163–170.  doi:10.1016/j.cpc.2016.08.020.
  * [3] A. Fedeli, C. Montecucco, G. L. Gragnani, Open-source software for electromagnetic scattering simulation: The case of antenna design, Electronics 8 (12) (2019) 1506.  doi:10.3390/electronics8121506.
  * [4] G. Yoon, J. Rho, MAXIM: Metasurfaces-oriented electromagnetic wave simulation software with intuitive graphical user interfaces, Computer Physics Communications (2021) 107846doi:10.1016/j.cpc.2021.107846.
  * [5] [link].
URL https://www.emworks.com/whitepapers

  * [6] P. F. Lopez, C. Arcambal, D. Baudry, S. Verdeyme, B. Mazari, Simple electromagnetic modeling procedure: From near-field measurements to commercial electromagnetic simulation tool, IEEE Transactions on Instrumentation and Measurement 59 (12) (2010) 3111–3121.  doi:10.1109/tim.2010.2063070.
  * [7] [link].
URL https://www.remcom.com/

  * [8] O. Jefimenko, Electricity and magnetism : an introduction to the theory of electric and magnetic fields, Electret Scientific Co, Star City, W. Va, 1989.
  * [9] D. Griffiths, Introduction to electrodynamics, Prentice Hall, Upper Saddle River, N.J, 1999.
  * [10] K. Yee, Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media, IEEE Transactions on Antennas and Propagation 14 (3) (1966) 302–307.  doi:10.1109/tap.1966.1138693.
  * [11] A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. Joannopoulos, S. G. Johnson, Meep: A flexible free-software package for electromagnetic simulations by the FDTD method, Computer Physics Communications 181 (3) (2010) 687–702.  doi:10.1016/j.cpc.2009.11.008.
  * [12] N. J. Champagne, J. G. Berryman, H. Buettner, FDFD: A 3d finite-difference frequency-domain code for electromagnetic induction tomography, Journal of Computational Physics 170 (2) (2001) 830–848.  doi:10.1006/jcph.2001.6765.
  * [13] M. Peskin, An introduction to quantum field theory, CRC Press, Boca Raton, FL, 2018\.
  * [14] R. Otin, J. Verpoorte, H. Schippers, R. Isanta, A finite element tool for the electromagnetic analysis of braided cable shields, Computer Physics Communications 191 (2015) 209–220.  doi:10.1016/j.cpc.2015.02.007.
  * [15] Jin, Finite Element Electromagnetic, John Wiley and Sons, 2014.
  * [16] G. Y. Delisle, K. L. Wu, J. Litva, Coupled finite element and boundary element method in electromagnetics, Computer Physics Communications 68 (1-3) (1991) 255–278.  doi:10.1016/0010-4655(91)90203-w.
  * [17] D. Komatitsch, J. Tromp, A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation, Geophysical Journal International 154 (1) (2003) 146–153.  doi:10.1046/j.1365-246x.2003.01950.x.
  * [18] W. Tang, S. xu Wang, Z. ru Wang, Perfectly matched layer absorbing boundary condition of finite element prestack reverse time migration in rugged topography, Procedia Engineering 37 (2012) 304–308.  doi:10.1016/j.proeng.2012.04.244.
  * [19] B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Mathematics of Computation 31 (139) (1977) 629–629.  doi:10.1090/s0025-5718-1977-0436612-4.
  * [20] A. Baev, Y. Kuznetsov, A. Gorbunova, M. Konovalyuk, J. A. Russer, Modeling of near-field to far-field propagator based on the jefimenko's equations, in: 2019 International Conference on Electromagnetics in Advanced Applications (ICEAA), IEEE, 2019.  doi:10.1109/iceaa.2019.8879083.
  * [21] Y. Kuznetsov, A. Baev, M. Konovalyuk, A. Gorbunova, J. A. Russer, M. Haider, Time-domain stochastic electromagnetic field propagator based on jefimenko's equations, in: 2018 Baltic URSI Symposium (URSI), IEEE, 2018.  doi:10.23919/ursi.2018.8406708.
  * [22] X.-M. Shao, Generalization of the lightning electromagnetic equations of uman, McLain, and krider based on jefimenko equations, Journal of Geophysical Research: Atmospheres 121 (7) (2016) 3363–3371.  doi:10.1002/2015jd024717.
  * [23] J.-J. Zhang, H.-Z. Wu, S. Pu, G.-Y. Qin, Q. Wang, Towards a full solution of the relativistic boltzmann equation for quark-gluon matter on GPUs, Physical Review D 102 (7) (2020) 074011.  doi:10.1103/physrevd.102.074011.
  * [24] H.-Z. Wu, J.-J. Zhang, L.-G. Pang, Q. Wang, Zmcintegral: a package for multi-dimensional monte carlo integration on multi-gpus, Comput. Phys. Commun. 248 (2020) 106962.  arXiv:1902.07916, doi:10.1016/j.cpc.2019.106962.
  * [25] J.-J. Zhang, H.-Z. Wu, Zmcintegral-v5: Support for integrations with the scanning of large parameter space on multi-gpus, Comput. Phys. Commun. 251 (2020) 107240.  arXiv:2020.107240, doi:10.1016/j.cpc.2020.107240.
  * [26] K. Yagi, T. Hatsuda, Y. Miake, Quark-Gluon Plasma, Cambridge University Press, 2008\.
  * [27] N. Andersson, G. L. Comer, Relativistic fluid dynamics: Physics for many different scales, Living Reviews in Relativity 10 (1) (jan 2007).  doi:10.12942/lrr-2007-1.
  * [28] J. Zhao, K. Zhou, S. Chen, P. Zhuang, Heavy flavors under extreme conditions in high energy nuclear collisions, Progress in Particle and Nuclear Physics 114 (2020) 103801.  doi:10.1016/j.ppnp.2020.103801.
  * [29] K. Fukushima, Evolution to the quark-gluon plasma, Rept. Prog. Phys. 80 (2) (2016) 022301.  arXiv:1603.02340, doi:10.1088/1361-6633/80/2/022301.
  * [30] J. Berges, M. P. Heller, A. Mazeliauskas, R. Venugopalan, Thermalization in QCD: theoretical approaches, phenomenological applications, and interdisciplinary connections (5 2020).  arXiv:2005.12299.
  * [31] R. Ryblewski, W. Florkowski, Early anisotropic hydrodynamics and the rhic early-thermalization and hbt puzzles, Phys. Rev. C 82 (2010) 024903.  arXiv:1004.1594, doi:10.1103/PhysRevC.82.024903.
  * [32] U. W. Heinz, P. F. Kolb, Two rhic puzzles: Early thermalization and the hbt problem, in: 18th Winter Workshop on Nuclear Dynamics (WWND 2002) Nassau, Bahamas, January 20-22, 2002, 2002, pp. 205–216.  arXiv:hep-ph/0204061.
  * [33] M. Gorelick, High Performance Python, O’Reilly Media, Inc, USA, 2014.
  * [34] J.-j. Zhang, R.-h. Fang, Q. Wang, X.-N. Wang, A microscopic description for polarization in particle scatterings, Phys. Rev. C 100 (6) (2019) 064904.  arXiv:1904.09152, doi:10.1103/PhysRevC.100.064904.
  * [35] P. Moritz, R. Nishihara, S. Wang, A. Tumanov, R. Liaw, E. Liang, M. Elibol, Z. Yang, W. Paul, M. I. Jordan, I. Stoica, Ray: A distributed framework for emerging ai applications, in: 13th USENIX Symposium on Operating Systems Design and Implementation (OSDI 18), USENIX Association, Carlsbad, CA, 2018, pp. 561–577.
URL https://www.usenix.org/conference/osdi18/presentation/moritz

  * [36] S. K. Lam, A. Pitrou, S. Seibert, Numba, in: Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC - LLVM, ACM Press, 2015.  doi:10.1145/2833157.2833162.
  * [37] R. Okuta, Y. Unno, D. Nishino, S. Hido, C. Loomis, Cupy: A numpy-compatible library for nvidia gpu calculations, in: Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS), 2017.