††thanks: These two authors contributed equally to this work.††thanks: These
two authors contributed equally to this work.

# Making Chiral Topological Superconductivities from Non-topological
Superconductivities Through the Twist

Yu-Bo Liu  School of Physics, Beijing Institute of Technology, Beijing 100081,
China  Jing Zhou  Department of Science, Chongqing University of Posts and
Telecommunications, Chongqing 400065, China  Yongyou Zhang  School of Physics,
Beijing Institute of Technology, Beijing 100081, China  Wei-Qiang Chen
Department of Physics and Shenzhen Institute for Quantum Science and
Engineering, Southern University of Science and Technology, Shenzhen 518055,
China  Shenzhen Key Laboratory of Advanced Quantum Functional Materials and
Devices, Southern University of Science and Technology, Shenzhen 518055, China
Fan Yang<EMAIL_ADDRESS>School of Physics, Beijing Institute of
Technology, Beijing 100081, China

###### Abstract

In this paper, we propose a general scheme to realize chiral TSCs through the
“twistronics”. Suppose we have a $D_{n}$-symmetric monolayer superconductor,
which carries non-topological SC with pairing angular momentum $L=n/2$. Here
we propose that we can obtain chiral TSC with the same $L$, by stacking two
such monolayers with the largest twist angle $\pi/n$, forming a Moireless
quasi-crystal (QC) structure, dubbed as the twist-bilayer QC (TB-QC) here. The
chiral TSC in the TB-QC is driven by the interlay Josephson coupling between
the pairing order parameters of the two layers. An argument based on the
universal Ginzburg-Landau theory is provided to understand this proposal. One
known example which fits our proposal is the $d+id$\- chiral TSC in the
45°-twisted bilayer cuprates. Here, based on the microscopic framework
developed previously to treat with the electron-electron interactions in the
TB-QC, we demonstrate the application of our proposal to a new example, i.e.
the $f+if$\- chiral TSC obtained by twisting two properly-doped honeycomb-
Hubbard-model monolayers by the angle 30°. This example is related to the
newly synthesized 30°- twisted bilayer graphene.

##  I Introduction

The origin and physical properties of topological quantum states is an
important research area, which has captured great research interests in recent
years Kane2010 ; Shoucheng2011 . Among these topological states, the chiral
topological superconductivities (TSCs) Read2000  are particularly interesting
since such quantum states are characterized by nonzero Chern numbers and
resultant Majorana zero modes in the vortex core or on the boundary
Volovik1999 , which carry non-Abelian statistics Ivanov2001  and can be used
in the design of fault-tolerant quantum computation Kitaev2003 ; Nayak2008 .
The chiral TSCs on a 2D lattice are usually generated by the nontrivial
rotational symmetries of the lattice. On an n-fold rotation symmetric 2D
lattice, the pairing symmetries can be classified according to the irreducible
representations (IRRPs) of the point group Sigrist1991 . In the special case
when the leading pairing symmetry belongs to the 2D IRRP, the two degenerate
pairing gap functions would usually be $1:\pm i$ mixed to lower the free
energy below the superconducting $T_{c}$ Cheng2010 , leading to the chiral
TSC. For example, the square lattice possessing $D_{4}$ symmetry can host the
$p+ip$-wave TSC with pairing angular momentum $L=1$ Yao2015 ; Meng2015 , while
the triangle-, the honeycomb- and the Kagome- lattices which possess the
$D_{6}$ symmetry can host $p+ip$\- or $d+id$\- wave TSC carrying pairing
angular momentum $L=1$ Ma2014 ; Chen2015 ; Zhang2015  or $L=2$ Doniach2007 ;
Gonzalez2008 ; Honerkamp2008 ; Pathak2010 ; McChesney2010 ; Nandkishore2012 ;
Wang2012 ; Kiesel2012 ; Honerkamp2014 ; Liufeng2013 ; ZhangLiDa2015 .

The material realization of the chiral TSC has long been a challenging
problem. The Sr2RuO4 used to be a promising candidate of the $p+ip$ chiral TSC
Maeno2003  on the square lattice, but now more and more experimental evidences
don’t support such a point of view. The quarter-doped graphene Doniach2007 ;
Gonzalez2008 ; Honerkamp2008 ; Pathak2010 ; McChesney2010 ; Nandkishore2012 ;
Wang2012 ; Kiesel2012 ; Honerkamp2014  and the properly-doped magic-angle
twisted bilayer graphene Xu2018 ; YangFan2018 ; WuFeng20181 ; Isobe2018 ;
You2019 ; Gonzalez2019  were proposed to host the $d+id$ chiral TSC, but
unambiguous experiment evidences are still lack. It’s interesting to ask the
question whether we can design chiral TSC from existing materials through some
engineering approach. Fortunately, the recently emergent “twistronics” gives
us the hope. Recently, it was proposed Yu_Bo2023 ; cuprates_QC ; JPHu2018 ;
cuprates_QC2 ; cuprates_QC3  that through twisting two cuprate monolayers by
the angles near 45°Zhu2021 ; Zhao2021 , one can obtain the $d+id$ chiral TSC
through the interlayer Josephson coupling (IJC). Here the cuprates monolayer
with $D_{4}$-symmetric square-lattice structure is already superconducting,
which hosts the $d$-wave SC carrying pairing angular momentum $L=2$. But the
$d$-wave pairing on the square lattice belongs to the 1D $B$\- IRRP with non-
degenerate real gap function, which is non topological. However, when two such
non-topological superconducting monolayes are stacked with the proper twist
angle, the chiral TSC can be achieved! Note that the special angle 45° here is
actually the largest possible twist angle between the two four-fold symmetric
monolayer. It’s interesting to generalize such a constructive proposal to more
lattices with different symmetries.

In this paper, we propose a general scheme to realize chiral TSCs through the
“twistronics”. Suppose we have a $D_{n}$-symmetric monolayer superconductor,
which carries non-topological SC with pairing angular momentum $L=n/2$, e.g.
the $d$-wave SC for square lattice or $f$-wave one for the triangle, honeycomb
or Kagome lattice. Here we propose that we can obtain chiral TSC with the same
pairing angular momentum, i.e. the $d+id$\- or $f+if$\- TSC, by stacking two
such monolayers with the largest twist angle $\pi/n$ between them, forming a
Moireless quasi-crystal (QC) structure, dubbed as the twist-bilayer QC (TB-QC)
here. The chiral TSC in the TB-QC is driven by the interlay Josephson coupling
(IJC) between the pairing order parameters (ODPs) of the two layers. An
argument based on the universal Ginzburg-Landau (G-L) theory is provided to
understand this proposal. Then based on the microscopic framework developed
previously to treat with the electron-electron (e-e) interactions in the TB-
QC, we demonstrate the application of our proposal with a new example
different from the known one, i.e. the $f+if$\- chiral TSC obtained by
twisting two properly-doped honeycomb-Hubbard-model monolayers by the angle
30°. This example is related to the newly synthesized 30°- twisted bilayer
graphene Ahn2018 ; Yao2018 ; Pezzini2020 ; Yan2019 ; Deng2020 .

The remaining part of the paper is organized as follow. In Sec. II, we provide
the G-L theory based analysis, which considers what pairing state would be
obtained in a TB-QC when each of its monolayer hosts a pairing state with the
largest pairing angular momentum for the lattice. In Sec.III, we provide an
example to demonstrate the conclusion achieved in Sec. II, i.e. a TB-QC with
each of its monolayer described by a Hubbard model on the honeycomb lattice.
In Sec. IV, a conclusion is arrived after some discussions.

##  II The G-L theory

We start from the classification of pairing symmetries on a 2D lattice
according to the IRRPs of its $D_{n}$ ($n$ is even hereafter) point group
Yu_Bo2023 . The $D_{n}$ point group has four 1D IRRPs and
$\left(\frac{n}{2}-1\right)$ 2D ones (labeled as $E_{L}$
($L\in\left[1,\frac{n}{2}-1\right]$)). For each 2D IRRP $E_{L}$, the two
degenerate basis gap functions would generally be mixed as $1:\pm i$ to lower
the free energy. The resultant gap function $\Delta^{(\pm)}_{L}(\mathbf{k})$
transform as $\Delta^{(\pm)}_{L}(\mathbf{k})\to e^{\mp
iL\Delta\phi}\Delta^{(\pm)}_{L}(\mathbf{k})$ under a $\Delta\phi=2\pi/n$
rotation, corresponding to a TSC with pairing angular momentum
$L\leq\frac{n}{2}-1$, and pairing chirality “$+$” or “$-$”. The four 1D IRRPs
correspond to the non-topological $A_{1,2}$ pairing symmetry with $L=0$ and
$B_{1,2}$ one with $L=\frac{n}{2}$. Here the label $L$ denotes the pairing
angular momentum. Clearly, for a $D_{n}$-symmetric lattice, the largest $L$ is
$\frac{n}{2}$, and the pairing with this $L$ is non topological. While for
$n=4$ the largest $L=2$ corresponds to the $d$-wave pairing, for $n=6$ the
largest $L=3$ corresponds to the $f$-wave one.

Figure 1: (Color online) Schematic illustration of the SC induced by IJC in a
TB-QC formed by two $D_{n}$-symmetric monolayers carrying SC with pairing
angular momentum $L=\frac{n}{2}$. We take $n=6$ for example. The color denotes
the gap sign on the Fermi surfaces.

Let’s take two $D_{n}$-symmetric monolayers, and stack them by the twist angle
$\frac{\pi}{n}$ to form a TB-QC, as shown in Fig. 1 for $n=6$. Obviously, the
point group is $D_{nd}$, isomorphic to $D_{2n}$. There is an additional
symmetry generator in the TB-QC absent its monolayer, i.e. the $C^{1}_{2n}$
rotation accompanied by a succeeding layer exchange, renamed as
$\tilde{C}^{1}_{2n}$ here.

Suppose that driven by some pairing mechanism, the monolayer $\mu=\text{t/b}$
(top/bottom) can host a pairing state with the largest pairing angular
momentum $L=n/2$. The pairing gap function in the $\mu$ layer is

$\Delta^{(\mu)}(\mathbf{k})=\psi_{\mu}\Gamma^{(\mu)}(\mathbf{k}).$ (1)

Here $\Gamma^{(\mu)}(\mathbf{k})$ is the normalized real form factor, and
$\psi_{\mu}$ is the “complex pairing amplitude”. Prominently, the
$\Gamma^{(\mu)}(\mathbf{k})$ for $L=n/2$ changes sign with every $C_{n}^{1}$
rotation, due to the following relation

$e^{iL\Delta\phi}=e^{i\pi}=-1.$ (2)

The geometry shown in Fig. 1 dictates

$\Gamma^{(\text{b})}(\mathbf{k})=\hat{P}_{\frac{\pi}{n}}\Gamma^{(\text{t})}(\mathbf{k}),~{}~{}\hat{P}_{\frac{2\pi}{n}}\Gamma^{(\mu)}(\mathbf{k})=-\Gamma^{(\mu)}(\mathbf{k}).$
(3)

Here $\hat{P}_{\phi}$ indicates the rotation by the angle $\phi$. As the
interlayer coupling in the TB-QC is weak Moon2019 ; Park2019 ; Yuan2020 ;
Yu_Bo2023 , we can only consider the dominant intralayer pairing. However, the
two intralayer pairing ODPs can couple through the IJC describing the combined
hopping of a Cooper pair between the two layers cuprates_QC ; JPHu2018 ;
cuprates_QC2 ; cuprates_QC3 . We shall investigate the ground state induced by
this IJC, based on a G-L theory analysis.

The symmetry-allowed free energy $F$ as function of $\psi_{\text{t/b}}$ can be
decomposed into the monolayers $F_{0}(\left|\psi_{\mu}\right|^{2})$ term and
the IJC $F_{J}$ term as cuprates_QC ; JPHu2018 ; cuprates_QC2 ; cuprates_QC3

$\displaystyle F\left(\psi_{\text{t}},\psi_{\text{b}}\right)$ $\displaystyle=$
$\displaystyle
F_{0}(\left|\psi_{\text{t}}\right|^{2})+F_{0}(\left|\psi_{\text{b}}\right|^{2})+F_{J}\left(\psi_{\text{t}},\psi_{\text{b}}\right),$
(4)

Up to the first-order IJC, the $F_{J}$ term should take the following
U(1)-gauge symmetry allowed form,

$\displaystyle F^{(1)}_{J}\left(\psi_{\text{t}},\psi_{\text{b}}\right)$
$\displaystyle=$
$\displaystyle-A\left(e^{i\theta}\psi_{\text{t}}\psi_{\text{b}}^{*}+c.c\right).$
(5)

The invariance of the free energy $F$ under the time-reversal (TR) operation:
$\psi_{\text{t/b}}\to\psi^{*}_{\text{t/b}}$ dictates $\theta=0$, leading to

$\displaystyle F^{(1)}_{J}\left(\psi_{\text{t}},\psi_{\text{b}}\right)$
$\displaystyle=$
$\displaystyle-A\left(\psi_{\text{t}}\psi_{\text{b}}^{*}+c.c\right).$ (6)

Note that the TB-QC possesses an additional symmetry absent in each of its
monolayer, i.e. the $\tilde{C}^{1}_{2n}$ symmetry. Under the
$\tilde{C}^{1}_{2n}$ operation, the gap function on the $\mu$ layer changes
from
$\Delta^{(\mu)}(\mathbf{k})=\psi_{\mu}\Gamma^{\left(\mu\right)}(\mathbf{k})$
to
$\tilde{\Delta}^{(\mu)}(\mathbf{k})=\psi_{\bar{\mu}}\hat{P}_{\frac{\pi}{n}}\Gamma^{\left(\bar{\mu}\right)}(\mathbf{k})$
which, under Eq. (3), can be rewritten as
$\tilde{\psi}_{\mu}\Gamma^{\left(\mu\right)}(\mathbf{k})$ with

$\tilde{\psi_{\text{b}}}=\psi_{\text{t}},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\tilde{\psi_{\text{t}}}=-\psi_{\text{b}}.$
(7)

The invariance of $F$ under $\tilde{C}^{1}_{2n}$ requires $A=0$, indicating
that the first-order IJC should vanish in the TB-QC.

Therefore, the following U(1)-gauge and TR symmetries allowed second-order IJC
should be considered,

$\displaystyle
F_{J}\left(\psi_{\text{t}},\psi_{\text{b}}\right)=A_{0}\left(\psi_{\text{t}}^{2}\psi_{\text{b}}^{2*}+{\rm
c.c.}\right)+O\left(\psi^{6}\right).$ (8)

Eq. (8) is minimized at $\psi_{b}=\pm i\psi_{t}$ for $A_{0}>0$ or
$\psi_{b}=\pm\psi_{t}$ for $A_{0}<0$. The form case is usually energetically
favored as the $1:i$ mixing manner of the ODPs from the two layers leads to
complex gap function which is fully gapped. In such a case, under the
$\tilde{C}^{1}_{2n}$, the pairing gap function
$\left(\psi_{b},\psi_{t}\right)$ would be changed to
$\left(\tilde{\psi_{b}},\tilde{\psi_{t}}\right)=\left(\psi_{t},-\psi_{b}\right)=\mp
i\left(\psi_{b},\psi_{t}\right)=e^{\mp i\frac{2\pi L}{2n}}$ with
$L=\frac{n}{2}$, suggesting a chiral TSC belonging to the $E_{n/2}$ IRRP of
the $D_{nd}$ or $D_{2n}$ point group.

To summarize, taking a TB-QC formed by two $D_{n}$-symmetric monolayer, when
each monolayer hosts a non-topological pairing state with pairing angular
momentum $L=n/2$, the TB-QC would most probably host a chiral TSC with the
same $L$, driven by the IJC.

While the sign of the above coefficient $A_{0}$ of the second-order IJC cannot
be determined by the G-L theory itself, it should be determined by the
microscopic calculations. Previous microscopic calculations favor $A_{0}>0$
for the 45°-twisted bilayer cuprates Yu_Bo2023 ; cuprates_QC . In the next
section, we shall do a microscopic calculation to determine the sign of the
coefficient $A_{0}$ for a TB-QC made from two monolayers described by the
honeycomb lattice Hubbard model which hosts the $f$-wave SC. It would be shown
that the $f+if$ chiral TSC would be obtained for this TB-QC.

##  III Microscopic calculations

Figure 2: (Color online) Results for the monolayer honeycomb lattice. (a) Band
structure along the high-symmetry lines. (b) FS for $\delta=0.17$ hole doping.
(c) The $\lambda\sim\delta$ relations for various leading pairing symmetries
at $U=0.3t$. The dashed line marks the doping level $\delta=0.17$. (d) The
distribution of the gap function for the $f$-wave SC on the FS.

In this section, we provide an example to demonstrate the universal conclusion
obtained on the above, which is different from the known one, i.e. the
45°twisted bilayer cuprates. Here we choose the honeycomb lattice as an
example for $n=6$, and study the Hubbard model. We shall first study the
pairing symmetry on a monolayer by the random-phase-approximation (RPA)
approach. It will be seen that in the phase diagram obtained, there exists a
doping regime in which the $f$-wave SC is the leading pairing symmetry. Then
we study the 30°-twisted bilayer system of the lattice. Through the
microscopic framework developed previously Yu_Bo2023 , we shall obtain the
band structure of this TB-QC. Using the RPA approach, we shall further study
the pairing state of the bilayer for a typical doping point at which the
monolayer hosts the $f$-wave pairing symmetry. Consequently, our results
unambiguously yield the $f+if$ chiral TSC state for this TB-QC.

###  III.1 The $f$-wave pairing for the monolayer

Let’s consider the following monolayer honeycomb-lattice Hubbard model with
nearest-neighbor (NN) and next-nearest-neighbor (NNN) hopping terms,

$\displaystyle H=$ $\displaystyle-$ $\displaystyle
t\sum_{\left\langle\mathbf{i,j}\right\rangle\sigma}c^{\dagger}_{\mathbf{i}\sigma}c_{\mathbf{j}\sigma}+h.c.-t^{\prime}\sum_{\left\langle\left\langle\mathbf{i,j}\right\rangle\right\rangle\sigma}c^{\dagger}_{\mathbf{i}\sigma}c_{\mathbf{j}\sigma}+h.c.$
(9) $\displaystyle+$ $\displaystyle
U\sum_{\mathbf{i}}n_{\mathbf{i}\uparrow}n_{\mathbf{i}\downarrow}$

The corresponding hopping parameters $t,t^{\prime}$ satisfy $t^{\prime}=0.1t$.
For the Hubbard interaction parameter $U$, we have taken a typical $U=0.3t$
friendly for the RPA calculations. The result for larger $U$ below the
critical one for the spin density wave instability is qualitatively the same.

The band structure of the model in the absence of $U$ is shown in Fig. 2(a)
along the high-symmetry lines, and the Fermi surface (FS) for a typical hole
doping $\delta=0.17$ is shown in Fig. 2(b). Obviously, Fig. 2(b) illustrates a
sextuple symmetric FS including two hole pockets centering around the
K-points. Using the standard multi-orbital RPA approach RPA1 ; RPA2 ; RPA3 ;
RPA4 , we obtain the largest pairing eigenvalue $\lambda$ and the
corresponding pairing eigenvector for each pairing symmetry for a given doping
level. Here $\lambda$ is related to the $T_{c}$ via $T_{c}\propto
e^{-\frac{1}{\lambda}}$, and therefore the pairing symmetry with the largest
$\lambda$ is the leading one. The relative pairing gap function of the leading
pairing symmetry is given by the normalized pairing eigenvector.

The doping $\delta$ dependence of the largest pairing eigenvalue $\lambda$ for
various leading pairing symmetries including the non-degenerate $s$-wave, the
degenerate $\left(p_{x},p_{y}\right)$-wave, the degenerate
$\left(d_{x^{2}-y^{2}},d_{xy}\right)$-wave and the non-degenerate $f$-wave are
shown in Fig. 2(c). Clearly, the $f$-wave pairing occupies the doping regime
$\delta\in(0.1,0.17)$ to be the leading pairing symmetry, which is
qualitatively consistent with Ref. Kiesel2012 . The distribution of the gap
function for the obtained $f$-wave pairing is shown in Fig. 2(d) on the FS.
This gap function changes sign with every 60°- rotation, and it hosts six
nodes along the Brillouin-Zone (BZ) diagonal direction.

###  III.2 The $f+if$-wave pairing for the TB-QC

Then let’s stack two monolayers described by Eq. (9) by the twist angle 30°,
to form a TB-QC. The total Hamiltonian now reads

$H=-\sum_{\mathbf{ij}\sigma}t_{\mathbf{ij}}c^{\dagger}_{\mathbf{i}\sigma}c_{\mathbf{j}\sigma}+U\sum_{\mathbf{i}}n_{\mathbf{i}\uparrow}n_{\mathbf{i}\downarrow}$
(10)

Here the index $\mathbf{i}$ labels all the sites belonging to both layers, and
$t_{\mathbf{ij}}$ represents the hopping integral between the sites
$\mathbf{i}$ and $\mathbf{j}$ which can locate in either the same or different
layers. For $\mathbf{i}$ and $\mathbf{j}$ locating within the same layer, the
formula of $t_{\mathbf{ij}}$ has been given by Eq. (9). The formula of
$t_{\mathbf{ij}}$ for the interlayer hopping is given as Moon2019

$t_{\mathbf{ij}}=t_{\mathbf{ij}\pi}\left[1-\left(\frac{\mathbf{R}_{\mathbf{ij}}\cdot\mathbf{e}_{\mathbf{z}}}{R}\right)^{2}\right]+t_{\mathbf{ij}\sigma}\left(\frac{\mathbf{R}_{\mathbf{ij}}\cdot\mathbf{e}_{\mathbf{z}}}{R}\right)^{2},$
(11)

with

$t_{\mathbf{ij}\pi}=t_{\pi}e^{-\left(R_{\mathbf{ij}}-a\right)/r_{0}},\quad
t_{\mathbf{ij}\sigma}=t_{\sigma}e^{-(R_{\mathbf{ij}}-d)/r_{0}}.$

Here, $R_{\mathbf{ij}}$ is the length of the 3D vector
$\mathbf{R}_{\mathbf{ij}}$ pointing from $\mathbf{i}$ to $\mathbf{j}$, and
$\mathbf{e}_{\mathbf{z}}$ is the unit vector perpendicular to the layer. The
parameters $a\approx 0.142$ nm, $d\approx 0.335$ nm, $r_{0}\approx 0.0453$ nm,
$t_{\pi}=t\approx 2.7$ eV and $t_{\sigma}\approx-0.48$ eV denote the lattice
constant, interlayer spacing, normalization distance, in-plane hoping and
vertical hoping, respectively. For these parameters, we have adopted the
corresponding parameters for the 30°-twisted bilayer graphene Moon2019 .

As the QC structure doesn’t possess translation symmetry, the traditional
band-structure theory cannot apply to the electronic structure of this
material. However, due to the large twist angle, the interlayer hybridization
is weak, and the perturbational-band theory Ahn2018 ; Yao2018 ; Moon2013 ;
Koshino2015 ; Yu_Bo2023  is suitable to treat with the electronic structure.
To involve the e-e interaction, we adopted the following revised
perturbational-band theory Yu_Bo2023  developed previously by some of the
authors of this paper.

Concretely, we decompose the tight-binding part of the Hamiltonian into the
zeroth-order intralayer hopping term $H_{0}$ and perturbational interlayer
tunneling term $H^{\prime}$. We first diagonalize $H_{0}$ in the
$\mathbf{k}$-space to obtain its eigen state
$\left|\mathbf{k}\alpha^{(\text{t/b})}\right\rangle$ and eigen energy
$\varepsilon_{\mathbf{k}\alpha}^{\text{t/b}}$. The $H^{\prime}$ can be written
as a hybridization form between top-layer state
$\left|\mathbf{k}\alpha^{(\text{t})}\right\rangle$ and bottom-layer state
$\left|\mathbf{q}\beta^{(\text{b})}\right\rangle$. Consequently, for a given
$\left|\mathbf{k}\alpha^{(\text{t})}\right\rangle$ state from the top layer,
only a few isolate $\left|\mathbf{q}\beta^{(\text{b})}\right\rangle$ state
from the bottom layer can couple with it, justifying the perturbational
treatment. Gathering all the $\left|\mathbf{q}\beta^{(\text{b})}\right\rangle$
related to $\left|\mathbf{k}\alpha^{(\text{t})}\right\rangle$, we can
calculate the perturbation-corrected eigen state and eigen energy brought
about by the $H^{\prime}$ term numerically, which are labeled as
$|\widetilde{\mathbf{k}\alpha^{(\text{t})}}\rangle$ and
$\tilde{\varepsilon}^{\text{t}}_{\mathbf{k}\alpha}$. Similarly, we get
$|\widetilde{\mathbf{q}\beta^{(\text{b})}}\rangle$ and
$\tilde{\varepsilon}^{\text{b}\beta}_{\mathbf{q}}$. We have checked that
different $|\widetilde{\mathbf{k}\alpha^{(\mu)}}\rangle$ thus obtained are
almost mutually orthogonal, qualifying
$\\{|\widetilde{\mathbf{k}\alpha^{(\mu)}}\rangle\\}$ as a good set of single-
particle bases to facilitate the succeeding studies involving e-e interaction.

The obtained band structure for the TB-QC is shown in Fig. 3(a), in comparison
with the two uncoupled band structures from the two separate monlayers, and
the FS for the $\delta=0.17$ hole doping is shown in Fig. 3(b). The most
prominent feature at the low hole-doping regime lies in that the two uncoupled
bands from the two separate layers cross at the $X$ point (or more generally
on the $\Gamma$-$X$ line) and strongly hybridize there, after which the two
bands are split into the lower band and the higher band. For the hole doping
level $\delta=0.17$ studied here, only the higher band crosses the Fermi
level, leading to a dodecagonal-symmetric FS, as shown in Fig. 3(b).

Figure 3: (Color online) Results for the $30\degree$-twisted bilayer honeycomb
lattice. (a) Band structure along the high-symmetry lines: solid (dashed)
lines for the coupled bilayer (two uncoupled monolayers). (b) FSs for
$\delta=0.17$ hole doping. The colors in (a) and (b) represent layer
component. (c) The $\lambda\sim\frac{U}{t}$ relations for various leading
pairing symmetries at $\delta=0.17$ hole doping, (d) Mixing-phase-angle
$\theta$ dependence of the energy for the degenerate $f$-wave pairings for
$U=t$.

When the Hubbard interaction is considered, we adopt the standard multi-
orbital RPA approach by using the above perturbation-corrected band structure.
Considering only intraband pairing between opposite momenta, we get the
effective BCS Hamiltonian. Under the mean-field (MF) treatment, we can obtain
the following linearized gap equation near the superconducting $T_{c}$ RPA1 ,

$-\frac{1}{(2\pi)^{2}}\sum_{\nu\beta}\oint
dq_{\parallel}\frac{V^{\mu\nu}_{\alpha\beta}(\mathbf{k},\mathbf{q})}{v^{\nu\beta}_{F}(\mathbf{q})}\Delta_{\nu\beta}(\mathbf{q})=\lambda\Delta_{\mu\alpha}(\mathbf{k}).$
(12)

where $V^{\mu\nu}_{\alpha\beta}(\mathbf{k},\mathbf{q})$ is the effective
pairing interaction given in Ref. Yu_Bo2023 . This equation is solved to yield
the largest pairing eigenvalue $\lambda$ and corresponding eigenvector
$\Delta_{\mu\alpha}(\mathbf{k})$. The former and latter determine the $T_{c}$
and the gap function respectively.

The $U/t$-dependence of the largest pairing eigenvalue $\lambda$ for various
pairing symmetries are shown in Fig. 3(c). The doping level is fixed at
$\delta=0.17$ hole doping, at which the $f$-wave pairing is the leading
pairing symmetry for the monolayer. Due to the classification according to the
IRRPs of the $D_{6d}$ point group, there can be non-degenerate $s$-, $i$-,
$i^{\prime}$-, $i*i^{\prime}$\- wave pairing symmetries and degenerate
$\left(p_{x},p_{y}\right)$-, $\left(d_{x^{2}-y^{2}},d_{xy}\right)$-,
$\left(f_{x^{3}-3xy^{2}},f_{3x^{2}y-y^{3}}\right)$-,
$\left(g_{x^{4}+y^{4}-6x^{2}y^{2}},g_{x^{3}y-xy^{3}}\right)$\- and
$\left(h_{x^{5}-10x^{3}y^{2}+5xy^{4}},h_{5x^{4}y-10x^{2}y^{3}+y^{5}}\right)$\-
wave pairing symmetries for this TB-QC. Here we only show the several leading
pairing symmetries.

Figure 4: (Color online) Distributions of the obtained pairing gap functions
on the FS of the 30°-twisted bilayer honeycomb lattice. The distributions of
the gap functions of the obtained $f_{x^{3}-3xy^{2}}$\- and
$f_{3x^{2}y-y^{3}}$-wave pairings are shown in (a) and (b) respectively. The
distributions of the amplitude and the phase of the gap function of the
$f+if$-wave pairing state obtained are shown in (c) and (d) respectively. The
hole doping level is $\delta=0.17$ and the interaction parameter is $U=t$.

From Fig. 3(c), the degenerate
$\left(f_{x^{3}-3xy^{2}},f_{3x^{2}y-y^{3}}\right)$-wave pairing symmetry is
the leading pairing symmetry for all the $U/t$ values shown. Note that
although the leading pairing symmetry for the monolayer system and the
30°-twisted bilayer one is both the $f$-wave carrying pairing angular momentum
3, it is non-degenerate for the former case and degenerate for the latter
case. The reason lies in that the point group has been enlarged from $D_{6}$
for the former case to $D_{6d}$ (isomorphic to $D_{12}$) for the latter case,
and the $f$-wave belongs to the $B_{2}$ and $E_{6}$ IRRPs for the two cases
respectively. The distributions of the pairing gap functions for the
degenerate $f_{x^{3}-3xy^{2}}$\- and $f_{3x^{2}y-y^{3}}$-wave pairing
symmetries are shown on the FSs in Fig. 4(a) and (b), respectively. While both
$f$-wave pairing gap functions change sign with every 60°-rotation, they
possess different symmetric and antisymmetric axes, as well as different nodal
lines. Clearly, the two gap functions are mutually related by 30°-rotation.

The two degenerate pairing components of the $f$-wave pairing possess the same
$T_{c}$, and would be mixed below $T_{c}$. We mix them as $1:e^{i\theta}$ to
minimize the ground-state energy. Here we have set $U=t$ instead of $U=0.3t$
so that the condensation energy is obviously larger than the machine accuracy.
Consequently, the $E(\theta)$ function shown in Fig. 3(d) is minimized at
$\theta=\pm\pi/2$, leading to the $f_{x^{3}-3xy^{2}}\pm if_{3x^{2}y-y^{3}}$
($f+if$ for abbreviation) -wave pairing state, consistent with the G-L theory.
The distributions of the amplitude and phase of the gap function for the
obtained $f+if$-wave pairing state are shown in Fig. 4(c) and (d),
respectively.

Figure 5: (Color online) Topological properties of the obtained $f+if$-wave
chiral TSC on the 30°-twisted bilayer honeycomb lattice. (a) The real-space
distribution of the squared amplitude of the wave function of a typical
Majorana zero-energy state. (b) The real-space distribution of the spontaneous
super current. The doping level and the interaction parameter are the same as
those in Fig. 4.

###  III.3 Topological Properties

The revised perturbational-band theory based microscopic framework adopted
here possesses clear advantages over the real-space approaches in the study of
the topological properties of TSCs in the TB-QC. In the weak-pairing limit
which applies to most superconductors including the TB-QCs studied here, the
Chern number for a fully-gapped pairing state is determined by the winding
number of the pairing phase around the FS Qi2010 ; Alicea2012 . As shown in
Fig. 4(c), the gap function of the $f+if$-wave pairing state is fully gapped,
which provides condition for the realization of TSC. Further more, Fig. 4(d)
shows that the distribution of the pairing phase repeats three times for each
run around the FS, leading to the winding number 3. Consequently, the Chern
number is 3. Therefore, we have obtained here the chiral $f+if$-wave TSC with
nontrivial high Chern number 3.

Usually, a chiral TSC is accompanied with Majorana chiral Fermion states and
spontaneous chiral super current on the boundary. To study such topological
properties, we take an open boundary condition which respects the dodecagonal
symmetry of the TB-QC, as shown in Fig. 5. Then, setting the obtained
$f+if$-wave pairing gap function as the input of the BCS MF Hamiltonian, we
diagonalize the Hamiltonian and obtain the spectrum and wave functions for the
Bogoliubov quasi-particles. Consequently, the real-space distribution of the
squared amplitude of the wave function of a typical zero-energy Bogoliubov
quasi-particle is shown in Fig. 5 (a). Obviously, the Majorana chiral Fermion
mainly distributes on the boundary of the TB-QC. Meanwhile, we calculate the
real-space distribution of the current. Consequently, there exists spontaneous
chiral super current which propagates along the boundary of the TB-QC, as
shown in Fig. 5 (b).

##  IV Discussion and Conclusion

In conclusion, we propose a general scheme to realize chiral TSCs from non-
topological SCs through the “twistronics”. Briefly, taking an n-fold symmetric
monolayer which hosts non-topological SC with pairing angular momentum
$L=n/2$, the corresponding TB-QC would host chiral TSC with the same pairing
angular momentum. Besides the known example of acquiring $d+id$ TSC in the
45°-twisted bilayer cuprates, here we propose another example in which we can
get the $f+if$\- TSC in a TB-QC made of two $D_{6}$ symmetric monolayers which
carry the $f$-wave pairing.

Note that the TB parameters adopted in our model are the same as those for the
30°-twisted bilayer graphene. However, the largest $U/t$ adopted in Fig. 3(c)
is only about $1$, less than the real value of about $2\sim 3$, because the
realistic $U/t$ has gone beyond the range that can be treated in the RPA
approach. In Ref. Kiesel2012 , the functional renormalization group based
study reveal that the $f$-wave pairing can indeed takes place in the low hole-
doped graphene with realistic interaction parameters. Then from our universal
G-L theory, the $f+if$-wave chiral TSC would indeed be realized in the low
hole-doped 30°-twisted bilayer graphene.

## Acknowledgements

This work is supported by the National Natural Science Foundation of China
under the Grant Nos.12074031, 12234016, 11674025, 12074037, and 12141402.
W.-Q. Chen is supported by the Science, Technology and Innovation Commission
of Shenzhen Municipality (No.ZDSYS20190902092905285), Guangdong Basic and
Applied Basic Research Foundation under Grant No.2020B1515120100 and Center
for Computational Science and Engineering of Southern University of Science
and Technology, and the National Key R and D Program of China (Grants No.
2022YFA1403700).

## Note

During the preparation of this manuscript, we became aware of a recent
preprintZhou2022  which suggested the possibility of chiral $f+if$ pairing in
the context of maximally twisted double-layer spin-triplet valley-singlet
superconductors.

## References

  * (1) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * (2) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * (3) N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
  * (4) G. E. Volovik, JETP Lett. 70, 609 (1999).
  * (5) D. A. Ivanov, Phys. Rev. Lett. 86, 268 (2001).
  * (6) A. Y. Kitaev, Ann. Phys. 303, 2 (2003).
  * (7) C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. D. Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * (8) M. Sigrist and K. Ueda, Rev. Mod. Phys. 63,239 (1991).
  * (9) M. Cheng, K. Sun, V. Galitski, and S. D. Sarma, Phys. Rev. B 81, 024504 (2010).
  * (10) H. Yao and F. Yang, Phys. Rev. B 92, 035132 (2015).
  * (11) Z.-Y. Meng, F. Yang, K.-S. Chen, H. Yao, and H.-Y. Kee, Phys. Rev. B 91, 184509 (2015).
  * (12) T.-X. Ma,F. Yang, H. Yao, and H.-Q. Lin, Phys. Rev. B 90, 245114 (2014).
  * (13) X. Chen, Y.-G. Yao, H. Yao, F. Yang,and J. Ni, Phys. Rev. B 92, 174503 (2015).
  * (14) L.-D. Zhang, F. Yang, and Y.-G. Yao, Phys. Rev. B 92, 104504 (2015).
  * (15) A. M. Black-Schaffer, and S. Doniach, Phys. Rev. B 75, 134512 (2007).
  * (16) J. González, Phys. Rev. B 78, 205431 (2008).
  * (17) C. Honerkamp, Phys. Rev. Lett. 100, 146404 (2008).
  * (18) S. Pathak, V. B. Shenoy, and G. Baskaran, Phys. Rev. B 81, 085431 (2010).
  * (19) J. L. McChesney, A. Bostwick, T. Ohta, T. Seyller, K. Horn, J. González, and E. Rotenberg, Phys. Rev. Lett. 104, 136803 (2010).
  * (20) R. Nandkishore, L. S. Levitov, and A. V. Chubukov, Nat. Phys. 8, 158 (2012).
  * (21) W.-S. Wang, Y.-Y. Xiang, Q.-H. Wang, F. Wang, F. Yang, and D.-H. Lee, Phys. Rev. B 85, 035414 (2012).
  * (22) M. L. Kiesel, C. Platt, W. Hanke, D. A. Abanin, R. Thomale, Phys. Rev. B 86, 020507(R) (2012).
  * (23) A. M. Black-Schaffer and C. Honerkamp, J. Phys. Condens. Matter 26, 423201 (2014).
  * (24) F. Liu, C.-C. Liu, K.-H. Wu, F. Yang, and Y.-G. Yao, Phys. Rev. Lett. 111, 066804 (2013).
  * (25) L.-D. Zhang, F. Yang and Y.-G. Yao, Sci. Rep. 5 : 8203 (2015).
  * (26) A. P. Mackenzie and Y. Maeno, Rev. Mod. Phys. 75, 657 (2003).
  * (27) C.-K. Xu and L. Balents, Phys. Rev. Lett. 121, 087001 (2018).
  * (28) C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Phys. Rev. Lett. 121, 217001 (2018).
  * (29) F.-C. Wu, A. H. MacDonald, and I. Martin, Phys. Rev. Lett. 121, 257001 (2018).
  * (30) H. Isobe, N. F. Q. Yuan, and L. Fu, Phys. Rev. X 8, 041041 (2018).
  * (31) Y.-Z. You and A. Vishwanath, NPJ Quan. Mater. 4, 1 (2019).
  * (32) J. González and T. Stauber, Phys. Rev. Lett. 122, 026801 (2019).
  * (33) Y.-B. Liu, Y.-Y. Zhang, W.-Q. Chen, and F. Yang, Phys. Rev. B 107, 014501 (2023).
  * (34) O. Can, T. Tummuru, R. P. Day, I. Elfimov, A. Damascelli, and M. Franz, Nat. Phys. 17, 519(2021).
  * (35) Z.-S. Yang, S.-S. Qin, Q. Zhang, C. Fang, and J.-P. Hu, Rev. B 98, 104515 (2018).
  * (36) A. Mercado, S. Sahoo, and M. Franz, Phys. Rev. Lett. 128, 137002 (2022).
  * (37) T. Tummuru, S. Plugge, and M. Franz, Phys. Rev. B 105, 064501 (2022).
  * (38) Y.-Y. Zhu, M.-H. Liao, Q.-H. Zhang, H.-Y. Xie, F.-Q. Meng, Y.-W. Liu, Z.-H. Bai, S.-H. Ji, J. Zhang, K.-L. Jiang, R.-D. Zhong, J. Schneeloch, G.-D. Gu, L. Gu, X.-C. Ma, D. Zhang, and Q.-K. Xue, Phys. Rev. X 11, 031011 (2021).
  * (39) S. Y. F. Zhao, N. Poccia, X.-M. Cui, P. A. Volkov, H. Yoo, R. Engelke, Y. Ronen, R.-D. Zhong, G.-D. Gu, S. Plugge, T. Tummuru, M. Franz, J. H. Pixley, P. Kim, arXiv: 2108.13455.
  * (40) S. J. Ahn, P. Moon, T.-H. Kim, H.-W. Kim, H.-C. Shin, E. H. Kim, H. W. Cha, S.-J. Kahng, P. Kim, M. Koshino, Y.-W. Son, C.-W. Yang, J. R. Ahn, Science 361, 782 (2018).
  * (41) W. Yao, E.-Y. Wang, C.-H. Bao, Y.-O. Zhang, K.-N. Zhang, K.-J. Bao, C.-K. Chan, C.-Y. Chen, J. Avila, M. C. Asensio, J.-Y. Zhu, and S.-Y. Zhou, PNAS 115, 6928 (2018).
  * (42) C. Yan, D.-L. Ma, J.-B. Qiao, H.-Y. Zhong, L. Yang, S.-Y. Li, Z.-Q. Fu, Y. Zhang and L. He, 2D Mater. 6, 045041 (2019).
  * (43) S. Pezzini, V. Mis̆eikis, G. Piccinini, S. Forti, S. Pace, R. Engelke, F. Rossella, K. Watanabe, T. Taniguchi, P. Kim and C. Coletti, Nano Lett. 20, 3313 (2020).
  * (44) B. Deng, B.-B. Wang, N. Li, R.-T. Li, Y.-N. Wang, J.-L. Tang, Q. Fu, Z. Tian, P. Gao, J.-M. Xue and H.-L. Peng, ACS Nano 14, 1656 (2020).
  * (45) P. Moon, M. Koshino, and Y.-W. Son, Phys. Rev. B 99, 165430 (2019).
  * (46) M. J. Park, H. S. Kim, and S. B. Lee, Phys. Rev. B 99, 245401(2019).
  * (47) G.-D. Yu, Z.-W. Wu, Z. Zhan, M. I. Katsnelson, and S.-J. Yuan, Phys. Rev. B 102, 115123 (2020).
  * (48) K. Kuroki, S. Onari, R. Arita, H. Usui, Y. Tanaka, H. Kontani and H. Aoki, Phys. Rev. Lett. 101, 087004 (2008).
  * (49) S. Graser, T. A. Maier, P. J. Hirschfeld and D. J. Scalapino, New Journal of Physics 11, 025016 (2009).
  * (50) T. A. Maier, S. Graser, P. J. Hirschfeld, and D. J. Scalapino, Phys. Rev. B 83, 100515(R) (2011).
  * (51) X.-X. Wu, F. Yang, C.-C. Le, H. Fan, and J.-P. Hu, Phys. Rev. B 92, 104511 (2015).
  * (52) P. Moon, M. Koshino, Phys. Rev. B 87, 205404 (2013).
  * (53) M. Koshino, New J. Phys. 17, 015014 (2015).
  * (54) X.-L. Qi, T. L. Hughes and S.-C. Zhang, Phys.Rev.B 82, 184516 (2010).
  * (55) J. Alicea, Rep. Prog. Phys. 75, 076501 (2012).
  * (56) B. T. Zhou, S. Egan, D. Kush, and M. Franz, arXiv: 2206.05599.


# Towards Programmable Memory Controller for Tensor Decomposition

Sasindu Wijeratne∗, Ta-Yang Wang∗, Rajgopal Kannan†, Viktor Prasanna∗
∗ University of Southern California, Los Angeles, USA
†US Army Research Lab, Los Angeles, USA
Email: {kangaram<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Tensor decomposition has become an essential tool in many data science
applications. Sparse Matricized Tensor Times Khatri-Rao Product (MTTKRP) is
the pivotal kernel in tensor decomposition algorithms that decompose higher-
order real-world large tensors into multiple matrices. Accelerating MTTKRP can
speed up the tensor decomposition process immensely. Sparse MTTKRP is a
challenging kernel to accelerate due to its irregular memory access
characteristics. Implementing accelerators on Field Programmable Gate Array
(FPGA) for kernels such as MTTKRP is attractive due to the energy efficiency
and the inherent parallelism of FPGA. This paper explores the opportunities,
key challenges, and an approach for designing a custom memory controller on
FPGA for MTTKRP while exploring the parameter space of such a custom memory
controller.

##  1 INTRODUCTION

Recent advances in collecting and analyzing large datasets have led to the
information being naturally represented as higher-order tensors. Tensor
Decomposition transforms input tensors to a reduced latent space which can
then be leveraged to learn salient features of the underlying data
distribution. Tensor Decomposition has been successfully employed in many
fields, including machine learning, signal processing, and network analysis
[Mondelli and Montanari, 2019, Cheng et al., 2020, Wen et al., 2020].
Canonical Polyadic Decomposition (CPD) is the most popular means of
decomposing a tensor to a low-rank tensor decomposition model. It has become
the standard tool for unsupervised multiway data analysis. The Matricized
Tensor Times Khatri-Rao product (MTTKRP) kernel is known to be the
computationally intensive kernel in CPD. Since the real-world tensors are
sparse, specialized hardware accelerators are becoming common means of
improving compute efficiency of sparse tensor computations. But external
memory access time has become the bottleneck due to irregular data access
patterns in sparse MTTKRP operation.

Since real-world tensors are sparse, specialized hardware accelerators are
attractive for improving the compute efficiency of sparse tensor computations.
But external memory access time is the bottleneck due to irregular data access
patterns in sparse MTTKRP operation.

Field Programmable Gate Arrays (FPGAs) are an attractive platform to
accelerate CPD due to the vast inherent parallelism and energy efficiency
FPGAs can offer. Since sparse MTTKRP is memory bound, improving the sustained
memory bandwidth and latency between the compute units on the FPGA and the
external DRAM memory can significantly reduce the MTTKRP compute time. FPGA
facilitates near memory computing with custom adaptive hardware due to its
reconfigurability and large on-chip BlockRAM memory [Xilinx, 2019]. It enables
the development of memory controllers and compute units specialized for
specific data formats; such customization is not supported on CPU and GPU.

The key contributions of this paper are:

  * •

We investigate possible sparse MTTKRP compute patterns and possible pitfalls
while adapting sparse MTTKRP computation to FPGA.

  * •

We scrutinize the importance of a FPGA-based memory controller design to
reduce the total memory access time of MTTKRP. Since MTTKRP on FPGA is a
memory-bound operation, it leads to significant acceleration in total MTTKRP
compute time.

  * •

We explore possible hardware solutions for Memory Controller design with
memory modules (e.g., DMA controller and cache) that can use to reduce the
overall memory access time.

The rest of the paper is organized as follows: Section 2 focuses on the
background of tensor decomposition and spMTTKRP. Section 3 and Section 4
investigate the compute patterns and memory access patterns of spMTTKRP.
Section 5 discusses the properties of configurable Memory Controller design.
Finally, we discuss the work in progress in Section 6.

##  2 BACKGROUND

###  2.1 Tensor Decomposition

Canonical Polyadic Decomposition (CPD) decomposes a tensor into a sum of 1D
tensors [Kolda and Bader, 2009]. For example, it approximates a 3D tensor
$\mathcal{X}\in\mathbb{R}^{I_{0}\times I_{1}\times I_{2}}$ as

$\mathcal{X}\approx\sum_{r=1}^{R}\lambda_{r}\cdot\mathbf{a}_{r}\otimes\mathbf{b}_{r}\otimes\mathbf{c}_{r}=:[\\![{\lambda};\mathbf{A},\mathbf{B},\mathbf{C}]\\!],$

where $R\in\mathbb{Z}_{+}$ is the rank, $\mathbf{a}_{r}\in\mathbb{R}^{I_{0}}$,
$\mathbf{b}_{r}\in\mathbb{R}^{I_{1}}$, and
$\mathbf{c}_{r}\in\mathbb{R}^{I_{2}}$ for $r=1,\ldots,R$. The components of
the above summation can be expressed as factor matrices, i.e.,
$\mathbf{A}=[\mathbf{a}_{1},\ldots,\mathbf{a}_{R}]\in\mathbb{R}^{I_{0}\times
R}$ and similar to $\mathbf{B}$ and $\mathbf{C}$. We normalize these vectors
to unit length and store the norms in
$\lambda=[\lambda_{1},\ldots,\lambda_{R}]\in\mathbb{R}^{R}$. Since the problem
is non-convex and has no closed-form solution, existing methods for this
optimization problem rely on iterative schemes.

The Alternating Least Squares (ALS) algorithm is the most popular method for
computing the CPD. Algorithm 1 shows a common formulation of ALS for 3D
tensors. In each iteration, each factor matrix is updated by fixing the other
two; e.g., $\mathbf{A}\leftarrow\mathcal{X}_{(0)}(\mathbf{B}\odot\mathbf{C})$.
This Matricized Tensor-Times Khatri-Rao product (MTTKRP) is the most expensive
kernel of ALS.

1Input: A tensor $\mathcal{X}\in\mathbb{R}^{I_{0}\times I_{1}\times I_{2}}$,
the rank $R\in\mathbb{Z}_{+}$

2 Output: CP decomposition
$[\\![{\lambda};\mathbf{A},\mathbf{B},\mathbf{C}]\\!]$,
$\lambda\in\mathbb{R}^{R}$, $\mathbf{A}\in\mathbb{R}^{I_{0}\times R}$,
$\mathbf{B}\in\mathbb{R}^{I_{1}\times R}$,
$\mathbf{C}\in\mathbb{R}^{I_{2}\times R}$

3 while  __stopping criterion not met__ do

4 $\mathbf{A}\leftarrow\mathcal{X}_{(0)}(\mathbf{B}\odot\mathbf{C})$

5 $\mathbf{B}\leftarrow\mathcal{X}_{(1)}(\mathbf{A}\odot\mathbf{C})$

6 $\mathbf{C}\leftarrow\mathcal{X}_{(2)}(\mathbf{A}\odot\mathbf{B})$

7 Normalize $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$ and store the norms as
$\lambda$

8

Algorithm 1 CP-ALS for the 3D tensors Figure 1: An illustration of MTTKRP for
a 3D tensor $\mathcal{X}$ in mode $0$

Figure 1 illustrates the update process of MTTKRP for mode $0$. With a sparse
tensor stored in the coordinate format, sparse MTTKRP (spMTTKRP) for mode $0$
can be performed as shown in Algorithm 1. For each non-zero element $x$ in a
sparse 3D tensor $\mathcal{X}\in\mathbb{R}^{I_{0}\times I_{1}\times I_{2}}$ at
$(i,j,k)$, the $i$th row of $\mathbf{A}$ is updated as follows: the $j$th row
of $\mathbf{B}$ and the $k$th row of $\mathbf{C}$ are fetched, and their
Hadamard product is computed and scaled with the value of $x$. The main
challenge for efficient computation is how to access the factor matrices and
non-zero elements for the spMTTKRP operation. We will use a hypergraph model
for modeling these dependencies in the spMTTKRP operation in Section 3.

Algorithm 2 shows the sequential sparse MTTKRP (spMTTKRP) approach for third-
order tensors in Coordinate (COO) Format, where IndI, indJ, and indK
correspond to the coordinate vectors of each non-zero tensor element. “nnz”
refers to the number of non-zero values inside the tensor.

Input: indI[nnz], indJ[nnz], indK[nnz], vals[nnz], $\mathbf{B}[J][R]$,
$\mathbf{C}[K][R]$

Output: $\mathbf{\tilde{A}}[I][R]$

1 for  _$z=0\emph{ to nnz}-1$_ do

2 $i=$ indI[z]

3 $j=$ indJ[z]

4 $k=$ indK[z]

5

6 for  _$r=0\emph{ to }R-1$_ do

7 $\mathbf{\tilde{A}}[i][r]$ += vals[z] $\cdot$ $\mathbf{B}[j][r]$ $\cdot$
$\mathbf{C}[k][r]$

8

9

10return $\mathbf{\tilde{A}}$

Algorithm 2 Single iteration of COO based spMTTKRP for third order tensors

###  2.2 FPGA Technologies

FPGAs are especially suitable for accelerating memory-bound applications with
irregular data accesses which require custom hardware. The logic cells on
state-of-the-art FPGA devices consist of Look Up Tables (LUTs), multiplexers,
and flip-flops. FPGA devices also have access to a large on-chip memory
(BRAM). High-bandwidth interfaces to external memory can be implemented on
FPGA. Current FPGAs are comprised of multiple Super Logic Regions (SLRs),
where each SLR is connected to a single or several DRAMs using a memory
interface IP.

HBM technology is used in state-of-the-art FPGAs as the high bandwidth
interconnections particularly benefit FPGAs [Kuppannagari et al., 2019]. The
combination of high bandwidth access to large banks of memory from logic
layers makes 3DIC architectures attractive for new approaches to computing,
unconstrained by the memory wall. Cache Coherent Interconnect supports shared
memory and cache coherency between the processor (CPU) and the accelerator.
Both FPGA and the processor have access to shared memory in the form of
external DRAM, while the cache coherency protocol ensures that any
modifications to a local copy of the data in either device are visible to the
other device. Protocols such as CXL [CXL, 2021] and CCIX [CCIX, 2021] develop
to realize coherent memory.

##  3 SPARSE MTTKRP COMPUTE PATTERNS

The spMTTKRP operation for a given tensor can be represented using a
hypergraph. For illustrative purposes, we consider a 3 mode sparse tensor
$\mathcal{X}\in\mathbb{R}^{I_{0}\times I_{1}\times I_{2}}$ where
$(i_{0},i_{1},i_{2})$ denote the coordinates of $x$ in $\mathcal{X}$. Here
$I_{0}$, $I_{1}$, and $I_{2}$ represent the size of each tensor mode. Note
that the following approaches can be applied to tensors with any number of
modes.

For a given tensor $\mathcal{X}$, we can build a hypergraph $H=(V,E)$ with the
vertex set $V$ and the hyperedge set $E$ as follows: vertices correspond to
the tensor indices in all the modes and hyperedges represent its non-zero
elements. For a 3D sparse tensor $\mathcal{X}\in\mathbb{R}^{I_{0}\times
I_{1}\times I_{2}}$ with $M$ non-zero elements, its hypergraph $H=(V,E)$
consists of $|V|=I_{0}+I_{1}+I_{2}$ vertices and $|E|=M$ hyperedges. A
hyperedge $\mathcal{X}(i,j,k)$ connects the three vertices $i$, $j$, and $k$,
which correspond to the indices of rows of the factor matrices. Figure 2 shows
an example of the hypergraph for a sparse tensor.

Figure 2: A hypergraph example of a sparse tensor

Our goal is to determine a mapping of $\mathcal{X}$ into memory for each mode
so that the total time spent on (1) loading tensor data from external memory,
(2) loading input factor matrix data from the external memory, (3) storing
output factor matrix data to the external memory, and (4) element-wise
computation for each non-zero element of the tensor is minimized.

Considering the current works in the literature [Srivastava et al., 2019]
[Nisa et al., 2019] [Li et al., 2018] [Helal et al., 2021], for a given mode,
there are 2 ways to perform sparse MTTKRP.

  * •

Approach 1: Output-mode direction computation

  * •

Approach 2: Input-mode direction computation

Algorithm 3 and Algorithm 4 show the mode $0$ MTTKRP of a tensor with three
modes using Approach 1 and Approach 2, respectively.

1Input: A sparse tensor $\mathcal{X}\in\mathbb{R}^{I_{0}\times I_{1}\times
I_{2}}$, dense factor matrices $\mathbf{{B}}\in\mathbb{R}^{I_{1}\times R}$,
$\mathbf{{C}}\in\mathbb{R}^{I_{2}\times R}$

2 Output: Updated dense factor matrix $\mathbf{A}\in\mathbb{R}^{I_{0}\times
R}$

3

4for  _each $i_{0}$ output factor matrix row in $\mathbf{A}$_ do

5 $\mathbf{A}(i_{0},:)=0$

6 for  _each nonzero element in $\mathcal{X}$ at $(i_{0},i_{1},i_{2})$ with
$i_{0}$ coordinates_ do

7 Load($\mathcal{X}(i_{0},i_{1},i_{2})$)

8 Load($\mathbf{B}(i_{1},:)$)

9 Load($\mathbf{C}(i_{2},:)$)

10 for  _$r=1,\ldots,R$_ do

11
$\mathbf{A}(i_{0},r)+=\mathcal{X}(i_{0},i_{1},i_{2})\times\mathbf{B}(i_{1},r)\times\mathbf{C}(i_{2},r)$

12

13

14 Store($\mathbf{A}(i_{0},:)$)

15

return $\mathbf{{A}}$

Algorithm 3 Approach 1 for mode 0 of a tensor with 3 modes

1Input: A sparse tensor $\mathcal{X}\in\mathbb{R}^{I_{0}\times I_{1}\times
I_{2}}$, dense factor matrices $\mathbf{{B}}\in\mathbb{R}^{I_{1}\times R}$,
$\mathbf{{C}}\in\mathbb{R}^{I_{2}\times R}$

2 Output: Updated dense factor matrix $\mathbf{A}\in\mathbb{R}^{I_{0}\times
R}$

3

4for  _each $i_{1}$ input factor matrix row in $\mathbf{B}$_ do

5 Load($\mathbf{B}(i_{1},:)$)

6 for  _each nonzero element in $\mathcal{X}$ at $(i_{0},i_{1},i_{2})$ with
$i_{i}$ coordinates_ do

7 Load($\mathcal{X}(i_{0},i_{1},i_{2})$)

8 Load($\mathbf{C}(i_{2},:)$)

9 for  _$r=1,\ldots,R$_ do

10
$\mathbf{p_{A}}(i_{0},r)=\mathcal{X}(i_{0},i_{1},i_{2})\times\mathbf{B}(i_{1},r)\times\mathbf{C}(i_{2},r)$

11

12 Store($\mathbf{p_{A}}(i_{0},:)$)

13

14

15 for  _each $i_{0}$ output factor matrix row in $\mathbf{A}$_ do

16 $\mathbf{A}(i_{0},:)=0$

17 for  _each partial element $\mathbf{p_{A}}$ with $i_{0}$ coordinates_ do

18 for  _$r=1,\ldots,R$_ do

19 Load($\mathbf{p_{A}}(i_{0},:)$)

20 $\mathbf{A}(i_{0},r)+=\mathbf{p_{A}}(i_{0},r)$

21

22 Store($\mathbf{A}(i_{0},:)$)

23

24

return $\mathbf{{A}}$

Algorithm 4 Approach 2 for mode 0 of a tensor with 3 modes

We use the hypergraph model of the tensor to describe the different
approaches. The main difference between these two approaches is the hypergraph
traversal order. Hence, we denote the two approaches based on the order of
hyperedge traversal.

In Approach 1, all hyperedges that share the same vertex of the output mode
are accessed consecutively. For each hyperedge, all the input vertices are
traversed to access their corresponding rows of input factor matrices. In
Approach 2, all hyperedges that share the same vertex of one of the input
modes are accessed sequentially. For each vertex, all its incident hyperedges
are iterated consecutively. For each hyperedge, the rest of the input vertices
of the hyperedge are traversed to access rows of the remaining input factor
matrices. It follows the element-wise multiplication and addition.

In Approach 1, since the order of hyperedge depends on the output mode, the
output factor matrix can be calculated without generating intermediate partial
sums (Algorithm 3: line 10). However, in Approach 2, since the hyperedges are
ordered according to the input mode coordinates, we need to store the partial
sums (Algorithm 4: line 9) in the FPGA external memory. It leads to
accumulating the partial sums to generate the output factor matrix (Algorithm
4: line 11-17).

The total computations of both approaches are the same: for a general sparse
tensor with $|T|$ non-zero elements, $N$ modes, and factor matrices with rank
$R$, since every hyperedge will be traversed once, and there are $N-1$
multiplications and one addition for computing MTTKRP, the total computation
per mode is $N\times|T|\times R$. However, their total external memory
accesses are different: both approaches require $|T|$ load operations for all
the hyperedges and the total factor matrix elements transferred per mode is
$(N-1)\times|T|\times R$, which corresponds to accessing input factor matrices
of vertices in the hypergraph model. However, in Approach 2, the partial value
needs to be stored in the memory (Algorithm 4: line 10), which requires an
additional $|T|\times R$ external memory storage. Let $I_{out}$ and $I_{in}$
represent the length of the output mode and the input mode, respectively. Then
the total amount of data transferred is $|T|+(N-1)\times|T|\times
R+I_{out}\times R$ for Approach 1 and $|T|+N\times|T|\times R+I_{in}\times R$
for Approach 2. Therefore, Approach 1 benefits from avoiding loading and
storing partial sums. Table 1 summarizes the properties of the two approaches.

Table 1: Comparison of the Approaches

Approach | Total Computations | Total external memory accesses | Size of total partial sums
---|---|---|---
1 | $N\times|T|\times R$ | $|T|+(N-1)\times|T|\times R+I_{out}\times R$ | $0$
2 | $N\times|T|\times R$ | $|T|+N\times|T|\times R+I_{in}\times R$ | $|T|\times R$

In the following, we discuss these in detail and identify the memory access
characteristics.

###  3.1 spMTTKRP on FPGA

In this paper, we consider large-scale data decomposition on very large
tensors. Hence the FPGA stores the tensor and the factor matrices inside their
external DRAM memory. Therefore, we need to optimize the FPGA memory
controller to the DRAM technology. In this section, we first explain the DRAM
timing model following the memory access patterns of sparse MTTKRP. Figure 3
shows the conceptual overall design.

Approach 2 is not practical for FPGA due to the large external memory
requirement to store the partial sums during the computation. In the work of
this paper, we focus on Approach 1.

For Approach 1, the tensor is sorted according to the coordinates of the
output mode. Typically, spMTTKRP is calculated for all the modes. To adapt
Approach 1 to compute the factor matrices corresponding to all the modes, (1)
Use multiple copies of the tensor. Each tensor copy is sorted according to the
coordinates of a tensor mode or (2) re-order the tensor in the output
direction before computing spMTTKRP for a mode.

Using multiple copies of a tensor is not a practical solution due to the
limited external memory of the FPGA. Hence, our memory solution focuses on
remapping the tensor in the output direction before computing spMTTKRP for a
mode. It enables to perform spMTTKRP using Approach 1 for each tensor mode.

Algorithm 5 summarizes the Approach 1 with remapping. The algorithm focus on
computing spMTTKRP of mode 1. Initially, we assume the sparse tensor is
ordered according to the coordinates of mode 0 after computing the factor
matrix of mode 0. Before starting the spMTTKRP for mode 1, we remap the
according to the mode 1 coordinates (lines 3 - 6). After remapping, all the
non-zero values with the same output mode coordinates are brought to the
compute unit consecutively (line 9). For each non-zero value, the
corresponding rows of the input factor matrices are brought into the compute
units following the element-wise multiplication and addition. Since the tensor
elements with the same output mode coordinates are brought together, the
processing unit can calculate the output factor matrix without storing the
partial values in FPGA external memory. Once a row of factor matrix is
computed, the value is stored in the external memory. Here, Load/Store
corresponds to loading/storing an element from the external memory. Also,
“$:$” refers to performing an operation for an entire factor matrix row.

The proposed approach introduces several implementation overheads:
Additional external memory accesses:

The remapping required additional external memory load and a store (Algorithm
5: lines 4 and 6). The total access to the external memory is increased by
$2\times|T|$ for a tensor of size $|T|$. The communication overhead per mode
is:

$\dfrac{2\times|T|}{|T|+(N-1)\times|T|\times R+I_{out}\times R}$
$\approx\dfrac{2}{1+(N-1)\times R}$

For a typical scenario (N = 3-5 and R = 16-64), the total external memory
communication only increases by less than $6\%$.
Additional external memory space:

During the remapping process, the remapped data requires an additional space
equal to the size of the tensor ($|T|$) to store the remapped tensor elements
in the memory.

Excessive memory address pointers to store the remapped tensor:

Figure 3: Conceptual overall design

The remapping brings the tensor elements with the same output mode coordinate
together (Algorithm 5: line 5). To achieve this, the memory controller needs
to track the memory location (i.e., memory address) of the next tensor element
with each output coordinate needs to be stored. This required memory address
pointers, which track the memory address to store a tensor element depending
on its output mode coordinate. Algorithm 5 requires such memory pointers
proportionate to the size of the output mode of a given tensor.

The number of address pointers may not fit in the FPGA internal memory for a
large tensor. For example, a tensor with an output mode with 10 million
coordinate values requires 40 MB to store the memory address pointers (i.e.,
32-bit memory addresses are considered). It does not fit in the FPGA on-chip
memory. Hence the address pointers should be stored in the external memory. It
introduces additional external memory access for each tensor element.

Also, the number of tensor elements with the same output mode coordinate value
is different for each output coordinate due to the sparsity of the tensor. It
complicates the memory layout of the tensor.

An ideal memory layout should guarantee: (1) The number of memory address
pointers required for remapping fit in the internal memory of the FPGA, and
(2) Each tensor partition contains the same number of tensor elements.

1Input: A sparse tensor $\mathcal{X}\in\mathbb{R}^{I_{0}\times I_{1}\times
I_{2}}$ sorted in mode 0, dense factor matrices
$\mathbf{{A}}\in\mathbb{R}^{I_{0}\times R}$,
$\mathbf{{C}}\in\mathbb{R}^{I_{2}\times R}$

2 Output: Updated dense factor matrix $\mathbf{B}\in\mathbb{R}^{I_{1}\times
R}$

3

4for  _each non-zero element in $\mathcal{X}$ at $(i_{0},i_{1},i_{2})$ with
$i_{0}$ coordinates_ do

5 Load($\mathcal{X}(i_{0},i_{1},i_{2})$)

6 pos${}_{i_{1}}$ = Find(Memory address of $i_{1}$)

7 Store($\mathcal{X}(i_{0},i_{1},i_{2})$ at memory address pos${}_{i_{1}}$)

8

9

10for  _each $i_{1}$ output factor matrix row in $\mathbf{B}$_ do

11 $\mathbf{B}(i_{1},:)=0$

12 for  _each non-zero element in $\mathcal{X}$ at $(i_{0},i_{1},i_{2})$ with
$i_{1}$ coordinates_ do

13 Load($\mathcal{X}(i_{0},i_{1},i_{2})$)

14 Load($\mathbf{A}(i_{0},:)$)

15 Load($\mathbf{C}(i_{2},:)$)

16 for  _$r=1,\ldots,R$_ do

17
$\mathbf{B}(i_{1},r)+=\mathcal{X}(i_{0},i_{1},i_{2})\times\mathbf{A}(i_{0},r)\times\mathbf{C}(i_{2},r)$

18

19

20 Store($\mathbf{B}(i_{1},:)$)

21

return $\mathbf{{B}}$

Algorithm 5 Approach 1 with remapping for mode 1 of a tensor with 3 modes

##  4 SPARSE MTTKRP MEMORY ACCESS PATTERNS

The proposed sparse MTTKRP computation has 5 main actions: (1) load a non-zero
tensor element, (2) load corresponding factor matrices, (3) perform spMTTKRP
operation, (4) store remapped tensor elements, and (5) store the final output.

The objective of the memory controller is to decrease the total DRAM memory
access time. To identify the opportunities to reduce the memory access time,
we analyze the memory access patterns of the proposed tensor format and memory
layout. The summary of memory access patterns is as follows:

  1. 1.

The tensor elements can be loaded as streaming accesses while remapping and
computing spMTTKRP.

  2. 2.

Each remapped tensor element can be stored element-wise.

  3. 3.

The different rows of each input factor matrices are random accesses.

  4. 4.

Each row of output factor matrix can be stored in streaming memory access.

Accessing the data in bulk (i.e., a large chunk of data stored sequentially)
can reduce the total memory access time. It is due to the characteristics of
the DRAM. DMA (Direct Memory Access) is the standard method to perform bulk
memory transfers. Further, the random accesses can be performed as element-
wise memory accesses through a cache to explore the temporal and spatial
locality of the accesses. It can improve the total access time.

Thus memory transfer types are as follows:

  1. 1.

Cache transfers: Supports random memory accesses. Load/store individual
requests in minimum latency. Access patterns with high spatial and temporal
locality are transferred using cache lines.

  2. 2.

DMA stream transfers: Supports streaming accesses. Load/store operations on
all requested data with minimum latency from memory.

  3. 3.

DMA element-wise transfers: Transfer the requested data element-wise. This
method is used with data with no spatial and temporal locality.

##  5 TOWARDS CONFIGURABLE MEMORY CONTROLLER

Figure 4: Proposed Memory Controller

To support the memory accesses, we propose a programmable memory controller as
shown in Figure 4. It consists of a Cache Engine, a tensor remapper, and a DMA
Engine. We evaluate the impact of using caches and DMAs as intermediate
buffering techniques to reduce the total execution time of sparse MTTKRP.

The modules inside the memory controller (e.g., Cache Engine, tensor remapper,
and DMA Engine) can be developed as configurable hardware. These modules are
programmable during the FPGA synthesis time. For example, the Cache Engine can
have a different number of cache lines and associativity that can be
configured during synthesis time. Also, the resource utilization of each
module heavily depends on the configuration. On the other hand, the FPGA
contains limited on-chip resources. Hence the FPGA resources should be
distributed among different modules optimally, such that the overall memory
access time is minimized (see Section 5.3).

###  5.1 Memory Controller Architecture

####  5.1.1 Cache Engine

The Cache Engine can be used to satisfy a single memory request with minimum
latency. The Cache Engine is used to explore the spatial and temporal locality
of requested data. We intend to use Cache Engine to explore the locality of
the input factor matrices. When a tensor computation requests a row of the
factor matrix, the memory controller first look-ups the Cache Engine. If the
requested factor matrix row is already in the Cache Engine due to prior
requests, the factor matrix row is forwarded to the tensor computation from
the Cache Engine. Otherwise, the tensor row is loaded to the Cache Engine from
the FPGA external memory. Then a copy of the matrix row is forwarded to the
computation while storing it in the cache.

####  5.1.2 DMA Engine

The DMA Engine can process bulk transfers between the compute units inside
FPGA and FPGA external memory. A DMA Engine can have several DMA buffers
inside.

The main advantages of having a DMA Engine are: (a) DMA requests can request
more than one element at once, unlike the Cache Engine, and reduce the input
traffic of the memory controller, (b) Using a DMA Engine to access data
without polluting the cache inside the Cache Engine, and (c) DMA transfers can
utilize the external memory bandwidth for bulk transfers.

####  5.1.3 Tensor Remapper

Tensor remapper includes a DMA buffer and the proposed remapping logic
discussed in Section 3. It loads each partition of the tensor as a bulk
transfer similar to the DMA Engine. After, it stores each tensor element
depending on the output mode coordinate value in an element-wise fashion.

Required memory consistency of the memory controller:

The suggested memory controller above should have a weak memory consistency
model with the following properties:

  * •

Consistency of DMA Engine, Cache Engine and Tensor Remapper: They process its
requests based on a first-in-first-out basis.

  * •

Consistency between Cache Engine, Tensor Remapper and DMA Engine:  The first-
in first-served basis is maintained. Since the same memory location is not
accessed by the Cache Engine, Tensor Remapper and DMA Engine at the same time,
weak consistency is maintained.

###  5.2 Programmable Parameters

The Cache Engine and DMA Engine use on-chip FPGA memory (i.e., BRAM and URAM).
These resources need to be shared among the modules optimally to achieve
significant improvements in memory access time. The resource requirement of
the Cache Engine and DMA Engine depends on their configurable parameters
mentioned below.

####  5.2.1 Memory Controller Parameters

Cache Engine parameters include cache line width, number of cache lines, and
associativity of the cache.

The design parameters of the DMA Engine are: the number of DMAs, the number of
DMA buffers per DMA, and the size of DMA buffers.

The design parameters of the Tensor Remapper include: (1) size of the DMA
buffer, (2) width of a tensor element, and (3) the maximum number of address
pointers Tensor Remapper can track.

Table 2: Characteristics of sparse tensors in FROSTT Repository

Metric | Value
---|---
Length of a tensor mode |  $17$-$39$ M
Width of a matrix $(R)$ |  $8-32$ (Typical = 16)
Number of non-zeros |  $3$-$144$ M
Number of modes |  $3$, $4$, $5$
Tensor size |  $\leq 2.25$ GB
Size of a factor matrix |  $<4.9$ GB

###  5.3 Exploring the Design Space

The tensor datasets can have different characteristics depending on the domain
from which the dataset is extracted. Table 2 shows the characteristics of the
tensors in The Formidable Repository of Open Sparse Tensors and Tools (FROSTT)
[Smith et al., 2017]. It is commonly used in the high-performance computing
community to benchmark custom accelerator designs for sparse MTTKRP.

Tensor datasets from separate domains of applications have different
characteristics such as sparsity, size of the modes, and the number of modes.
Hence, the datasets extracted from various applications show the least memory
access time with different configurations of the memory controller. Hence,
performance estimator software is required to estimate the optimal
configurable parameters for datasets of a domain. We introduce the features of
a Performance Model Simulator (PMS) software to estimate the total execution
time of spMTTKRP for a given dataset. It can use with multiple datasets from
the same domain to estimate the average execution time ($t_{avg}$) for a
selected domain. Also, it should estimate the total FPGA on-chip memory
requirement for a given set of programmable parameters to make sure the memory
controller fits in the FPGA device. We will explore the possible inputs
required for a PMS concerning: (1) available FPGA resources (i.e., total
BRAMs, and URAMs of the selected FPGA and data width of memory interface), (2)
size of data structures (e.g., size of an input tensor element, size of an
input factor matrix element, and rank of the input factor matrices), and (3)
Parameters of the memory controller (i.e., DMA buffer sizes, number of cache
lines, associativity of cache, and number of factor matrices shared by a
cache).

A module-by-module (e.g., Cache Engine and DMA Engine) exhaustive parameter
search can be proposed to identify the optimal parameters for the memory
controller.

##  6 DISCUSSION

In this paper, we investigated the characteristics of a custom memory
controller that can reduce the total memory access time of sparse MTTKRP on
FPGAs. Sparse MTTKRP is a memory-bound operation. It has 2 types of memory
access patterns that can be optimized to reduce the total memory access time.
A memory controller design that can be configured during compile/synthesis
time depending on the application and targeted hardware is required.

We are developing a configurable memory controller and a memory layout for
sparse tensors to reduce the total memory access time of sparse MTTKRP
operation.

Since synthesizing a FPGA can take a long time, optimizing the memory
controller parameters for a given application can be a time-consuming process.
Hence, we are developing a Performance Model Simulator (PMS) software to
identify the optimal parameters for a given application on a selected FPGA.

## ACKNOWLEDGEMENTS

This work was supported by the U.S. National Science Foundation (NSF) under
grants NSF SaTC # 2104264 and PPoSS- 2119816.

## REFERENCES

  * CCIX, 2021 CCIX (2021).  Cache Coherent Interconnect for Accelerators (CCIX).  https://www.ccixconsortium.com/.
  * Cheng et al., 2020 Cheng, Z., Li, B., Fan, Y., and Bao, Y. (2020).  A novel rank selection scheme in tensor ring decomposition based on reinforcement learning for deep neural networks.  In ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3292–3296. IEEE.
  * CXL, 2021 CXL (2021).  Compute Express Link (CXL).  https://www.computeexpresslink.org/.
  * Helal et al., 2021 Helal, A. E., Laukemann, J., Checconi, F., Tithi, J. J., Ranadive, T., Petrini, F., and Choi, J. (2021).  Alto: Adaptive linearized storage of sparse tensors.  In Proceedings of the ACM International Conference on Supercomputing, ICS ’21, page 404–416, New York, NY, USA. Association for Computing Machinery.
  * Kolda and Bader, 2009 Kolda, T. G. and Bader, B. W. (2009).  Tensor decompositions and applications.  SIAM review, 51(3):455–500.
  * Kuppannagari et al., 2019 Kuppannagari, S. R., Rajat, R., Kannan, R., Dasu, A., and Prasanna, V. (2019).  Ip cores for graph kernels on fpgas.  In 2019 IEEE High Performance Extreme Computing Conference (HPEC), pages 1–7.
  * Li et al., 2018 Li, J., Sun, J., and Vuduc, R. (2018).  Hicoo: Hierarchical storage of sparse tensors.  In Proceedings of the International Conference for High Performance Computing, Networking, Storage, and Analysis, SC ’18. IEEE Press.
  * Mondelli and Montanari, 2019 Mondelli, M. and Montanari, A. (2019).  On the connection between learning two-layer neural networks and tensor decomposition.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1051–1060. PMLR.
  * Nisa et al., 2019 Nisa, I., Li, J., Sukumaran-Rajam, A., Vuduc, R., and Sadayappan, P. (2019).  Load-balanced sparse mttkrp on gpus.  In 2019 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 123–133.
  * Smith et al., 2017 Smith, S., Choi, J. W., Li, J., Vuduc, R., Park, J., Liu, X., and Karypis, G. (2017).  FROSTT: The formidable repository of open sparse tensors and tools.
  * Srivastava et al., 2019 Srivastava, N., Rong, H., Barua, P., Feng, G., Cao, H., Zhang, Z., Albonesi, D., Sarkar, V., Chen, W., Petersen, P., Lowney, G., Herr, A., Hughes, C., Mattson, T., and Dubey, P. (2019).  T2s-tensor: Productively generating high-performance spatial hardware for dense tensor computations.  In 2019 IEEE 27th Annual International Symposium on Field-Programmable Custom Computing Machines (FCCM), pages 181–189.
  * Wen et al., 2020 Wen, F., So, H. C., and Wymeersch, H. (2020).  Tensor decomposition-based beamspace esprit algorithm for multidimensional harmonic retrieval.  In ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 4572–4576. IEEE.
  * Xilinx, 2019 Xilinx (2019).  Alveo u250 data center accelerator card.  https://www.xilinx.com/products/boards-and-kits/alveo/u250.html.


[c,d]Anabella Araudo

# Using PIC and PIC-MHD to investigate cosmic ray acceleration in mildly
relativistic shocks

Artem Bohdan  Allard Jan van Marle  Fabien Casse  Alexandre Marcowith

###### Abstract

Astrophysical shocks create cosmic rays by accelerating charged particles to
relativistic speeds. However, the relative contribution of various types of
shocks to the cosmic ray spectrum is still the subject of ongoing debate.
Numerical studies have shown that in the non-relativistic regime, oblique
shocks are capable of accelerating cosmic rays, depending on the Alfvénic Mach
number of the shock. We now seek to extend this study into the mildly
relativistic regime. In this case, dependence of the ion reflection rate on
the shock obliquity is different compared to the nonrelativistic regime.
Faster relativistic shocks are perpendicular for the majority of shock
obliquity angles therefore their ability to initialize efficient DSA is
limited. We define the ion injection rate using fully kinetic PIC simulation
where we follow the formation of the shock and determine the fraction of ions
that gets involved into formation of the shock precursor in the mildly
relativistic regime covering a Lorentz factor range from 1 to 3. Then, with
this result, we use a combined PIC-MHD method to model the large-scale
evolution of the shock with the ion injection recipe dependent on the local
shock obliquity. This methodology accounts for the influence of the self-
generated or pre-existing upstream turbulence on the shock obliquity which
allows study substantially larger and longer simulations compared to classical
hybrid techniques.

##  1 Introduction

The origin and acceleration mechanisms of cosmic rays have been the subject of
intense research in astrophysics. It is widely accepted that the bulk of
galactic cosmic rays are produced in the aftermath of supernova explosions,
where particles gain energy through Diffusive Shock Acceleration (DSA) [1].
However, the origin of Ultra-High-Energy Cosmic Rays (UHECRs) remains an
enigma, and the quest to understand their sources continues to captivate the
scientific community. One intriguing candidate for accelerating particles to
extreme energies is found in the mildly relativistic shocks within the
backflows of Active Galactic Nuclei jets [2]. These shocks have the potential
to accelerate particles up to energies of $10^{18}$ eV and beyond.

In this study, we focus on mildly relativistic shocks with Lorentz factors
ranging from 1 to 3, exploring their potential as sources of high-energy
particles. The efficiency of DSA in these shocks depends heavily on their
specific parameters, making them intriguing candidates for understanding the
origin of UHECRs. By delving into the intricate details of particle
acceleration within these shocks, we aim to shed light on the mechanisms
responsible for producing cosmic rays with such extreme energies.

To gain insight into the acceleration processes, we need to combine two
different methods, similarly how it was done for nonrelativistic oblique
supernova remnant shocks [3]. Firstly, we need to use the particle-in-cell
(PIC) method to model the shock structure in order to determine the fraction
of particles that gets reflected into the upstream medium. Once this has been
established, we can use this injection rate to run further simulations with
the PIC-magnetohydrodynamic (MHD) method to determine whether the injection
rate suffices to trigger the non-resonant streaming instability [4, 5].

##  2 PIC simulations

PIC shock simulations are performed using an optimized fully-relativistic
electromagnetic 2D PIC code with MPI-openMP hybrid parallelization developed
from the TRISTAN code [6, 7]. Shocks are initialized using the reflecting-wall
setup (Fig. 1(a)). Interaction of the upstream electron-ion plasma flow with
the reflecting wall results in a shock which propagates towards the upstream
plasma flow in positive $x$ direction. The large-scale magnetic field, $\bf
B_{0}$, is carried by the upstream plasma. The shock obliquity angle,
$\theta_{B}$, is defined as the angle between $\bf B_{0}$ and the shock normal
vector. Here we utilise the _out-of-plane_ magnetic field configuration,
therefore $\varphi=90\deg$, where $\varphi$ the angle between the upstream
magnetic field and y-axis which lies within the simulation plane. The upstream
plasma velocity is chosen to have the resulting Lorentz factor of the shock of
$\gamma_{\rm sh}=1.86$. This sets the Alfvénic and sonic Mach numbers as
$M_{\rm A}\approx 20$ and $M_{\rm s}\approx 22$, which are defined in the
upstream reference frame. The plasma beta value, which denotes the thermal-to-
magnetic energy density ratio in the upstream region is $\beta=1$. The ion-to-
electron mass ratio is set to $m_{i}/m_{e}=25$. To ensure adequate resolution
of all important timescales in our simulation, we use a simulation timestep of
$\delta t=1/16\leavevmode\nobreak\ \omega_{\mathrm{pe}}^{-1}$ and a resolution
of $\lambda_{\mathrm{se}}=c/\omega_{\mathrm{pe}}=8\Delta$, where $\Delta$
represents the size of one cell in our simulation grid. The number of
particles per cell per species is $N_{\mathrm{ppc}}=40$.

We conducted a comprehensive study with 6 PIC simulations, varying the
obliquity angle $\theta_{B}$, as presented in Table 1. The Lorentz factor of
the shock is fixed at $\gamma_{\rm sh}=1.86$. Depending on the value of
$\theta_{B}$, the shock can be categorized as perpendicular (superluminal
case), oblique (where extra energy is required for a particle to escape the
shock), or parallel (where no extra energy is needed for a particle to escape
the shock), with the critical transitions occurring at $\theta_{B}=27.8$
(parallel-to-oblique) and $\theta_{B}=32$ (oblique-to-perpendicular), measured
in the upstream reference frame.

To study the shock evolution, we followed the simulations for a duration of
$220\,\Omega_{\mathrm{i}}$, during which the electromagnetic shock precursor
formed and the ion injection stabilized. The parameters of the stabilized ion
injection are summarized in Table 1. As expected, the number density of shock-
reflected particles remains stable for parallel shocks (PIC01-04). For oblique
shocks (PIC05), the number density decreases, and no injection occurs in the
superluminal case (PIC06). The normalized energy of injected particles,
$U_{\rm inj}/U_{0}$, is found to be minimal for simulations PIC01-02, as
particles do not require extra energy to outrun the shock in these cases. In
contrast, simulations PIC03-04 show an increase in $U_{\rm inj}/U_{0}$,
indicating that Shock Drift Acceleration (SDA) is active at these obliquities,
despite the shock remaining parallel, and ions not requiring extra energy to
escape the shock.

Figure 1: Simulation setups. Name | $\theta_{B}$ | $\theta_{B}$ | $N_{\rm inj}/N_{\rm 0}$ | $U_{\rm inj}/U_{\rm sh}$ | $U_{\rm inj}/U_{\rm sh}$
---|---|---|---|---|---
| SHRF | UFR |  | URF | SHRF
PIC01 | 0 | 0 | 0.025 | 3.2 | 0.4
PIC02 | 17.7 | 6.5 | 0.019 | 3.7 | 0.5
PIC03 | 33.4 | 13.3 | 0.021 | 7.3 | 1.8
PIC04 | 46.3 | 20.5 | 0.023 | 12.6 | 4.3
PIC05 | 56.7 | 28.6 | 0.013 | 10.9 | 4.6
PIC06 | 61.1 | 32.9 | 0 | - | -
Table 1: PIC Simulation Parameters. $\theta_{B}$ is the obliquity angle
measures in the upstream reference frame (URF) and the shock reference frame
(SHRF). For the shock Lorentz factor of 1.86, the parallel-to-oblique
transition at $\theta_{B}=27.8$, the oblique-to-perpendicular transition at
$\theta_{B}=32$, where $\theta_{B}$ is measured in URF. $U_{\rm inj}/U_{\rm
sh}=(\gamma_{\rm inj}-1)/(\gamma_{\rm sh}-1)$, where $\gamma_{\rm inj}$ is the
average energy of injected particles in the respective reference frame and
$(\gamma_{\rm sh}-1)mc^{2}$ is the energy of the upstream particles in the
shock reference frame.

##  3 PIC-MHD simulation

To investigate the large scale evolution of the shock, we perform a 2.5D PIC-
MHD simulation, using a relativistic MHD version of the code described in [8],
which is based on the MPI-AMRVAC code [9]. This code combines aspects of both
traditional MHD and PIC in that it treats the thermal plasma as a fluid
through MHD and the non-thermal charged particles through the PIC method.
These two components interact with each other through the equation of motion
for the particles and a modified version of Ohm’s law for the fluid [10]. We
set up our simulations in the rest-frame of the shock (Fig. 1(b)), starting
from the analytical solution of the Rankine-Hugoniot conditions. Once the
simulation has started, we begin to introduce the non-thermal particles at the
shock front according to the injection rate given in Table 1.

For the PIC-MHD simulation, we start from a standing shock with a Lorentz
factor of 1.86, and an obliquity of 30 degrees in the upstream rest-frame. At
the shock, we inject non-thermal protons at a rate of 0.5 percent of the
number of ions passing through the shock (hence compatible with PIC05
results). The non-thermal particles have a Lorentz factor of 3 times that of
the pre-shock gas and are injected with an isotropic velocity distribution in
the rest-frame of the shock. The simulation space has a box-size of
180$\times$30 $R_{L}$ with $R_{L}$ the Larmour radius defined by the velocity
of the injected particles and the upstream magnetic field. In order to assure
adequate resolution of the particle motion, we use a fixed grid with
individual cells of 0.15 $R_{L}$.

Figure 2 shows the result of the PIC-MHD simulation, displaying the magnetic
field relative to the initial upstream field, the relative density of non-
thermal to thermal particles, and the thermal gas density relative to the
initial upstream density at $t\,=\,100$ and $500\,R_{L}/c$. As can be seen,
magnetic field shows both upstream and downstream distortion, which enables
the reflection of non-thermal particles across the shock. However, the
upstream disturbance is relatively small. As a result, we would expect that
diffusive shock acceleration will be relatively inefficient, since many
particles will escape without being reflected. This situation is further
enhanced by the size of the perturbed pre-shock region.

Figure 2: PIC-MHD results. From top to bottom, the figures show the magnetic
field relative to the initial upstream field, the ratio of non-thermal to
thermal particles, and the thermal gas density relative to the initial
upstream density at $t\,=\,100\ R_{L}/c$ (left) and $500\ R_{L}/c$ (right).
The bottom panel also includes the magnetic field lines.

The spectral energy distributions resulting from the PIC-MHD simulation are
shown in Fig, 3. Initially, the spectrum is broad, but lacks an exponential
high-energy tail, indicating that shock-drift is the primary acceleration
mechanism. A high-energies the tail starts to appear over time, signalling
that diffusive shock acceleration has started. However, a longer simulation is
required in order to determine whether diffusive shock acceleration can become
efficient in this particular situation.

Figure 3: Spectral energy distribution over time for the PIC-MHD simulation.

##  4 Conclusions

  * •

Utilizing a fully-relativistic electromagnetic 2D3V PIC code, we perform a
series of shock simulations with varying obliquity angles ($\theta_{B}$) to
study the shock dynamics and particle acceleration processes. Stable shock
reflected particle densities are observed for parallel shocks, whereas it
decreases for oblique shocks and it is absent in the superluminal case.
Results are consistent with theoretical predictions for the selected Lorentz
factor of the shock.

  * •

The PIC-MHD simulation reveals that the magnetic field exhibits both upstream
and downstream turbulence, enabling the reflection of non-thermal particles
across the shock. However, the upstream disturbance is relatively small,
leading to the expectation that diffusive shock acceleration will be
relatively inefficient and longer simulations are required to determine
whether diffusive shock acceleration will become efficient in this specific
scenario.

## Acknowledgments

A. B. was supported by the German Research Foundation (DFG) as part of the
Excellence Strategy of the federal and state governments - EXC 2094 -
390783311. Computations were performed on the HPC system Raven at the Max
Planck Computing and Data Facility.

## References

  * Drury [1983] Drury, L.O’.C., 1983, Rep. Prog. Phys., 46, 973
  * Bell et al. [2019] Bell, A. R., Matthews, J. H., Blundell, K. M., et al. 2019, MNRAS, 487, 4571
  * van Marle et al. [2022] van Marle, A. J., Bohdan, A., Morris, P. J., et al. 2022, ApJ, 929, 7.
  * Bell [2004] Bell, A. R. 2004, MNRAS, 353, 550
  * Bell [2005] Bell, A. R. 2005, MNRAS, 358, 181
  * Buneman, O. [1993] Buneman, O., Computer Space Plasma Physics: Simulation Techniques and Software Eds.: H. Matsumoto & Y. Omura, Tokyo: Terra Scientific, 1993
  * Niemiec et al. [2008] Niemiec, J., Pohl, M., Stroman, T., et al. 2008, ApJ, 684, 1174
  * van Marle et al. [2018] van Marle, A. J., Casse, F., & Marcowith, A. 2018, MNRAS, 473, 3394
  * van der Holst et al. [2008] van der Holst, B., Keppens, R., & meliani, Z. 2008, CPC, 179, 617
  * Bai et al. [2015] Bai, X-N., Caprioli, D., Sironi, L., & Spitkovsky, A. 2015, ApJ, 809, 55