# Field-free Josephson diode effect in a $d$-wave superconductor
heterostructure

Hamed Vakili  Department of Physics and Astronomy and Nebraska Center for
Materials and Nanoscience, University of Nebraska, Lincoln, Nebraska 68588,
USA  Moaz Ali  Department of Physics and Astronomy and Nebraska Center for
Materials and Nanoscience, University of Nebraska, Lincoln, Nebraska 68588,
USA  Alexey A. Kovalev  Department of Physics and Astronomy and Nebraska
Center for Materials and Nanoscience, University of Nebraska, Lincoln,
Nebraska 68588, USA

(August 28, 2024)

###### Abstract

We study superconductor/normal region/superconductor (S$|$N$|$S) Josephson
junction formed using superconductors with $d$, $d+id^{\prime}$, and $d+is$
superconducting pairings. We show that the quality factor of the Josephson
diode effect and its sign can be substantially tuned by the external magnetic
field, gate voltage, and the length of the junction for all three types of
pairings. We also identify the conditions under which the anomalous Josephson
and Josephson diode effects can appear in the junction by analyzing
appropriate symmetries. In particular, by breaking a $\pi$-rotation symmetry,
we show how a large field-free Josephson diode effect can be realized even in
the absence of spin-orbit coupling. We also study the role of edge states
appearing in the case of chiral superconductor with $d+id^{\prime}$ pairing.
Our results demonstrate that the Josephson diode effect in a planar geometry
can be used as a signature of unconventional superconducting pairings.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

The superconducting diode (SDE) and Josephson diode effects (JDE) have
attracted considerable attention in recent years due to many potential
applications in superconducting devices and quantum computing [1, 2, 3, 4, 5,
6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. The effects are characterized by non-
reciprocal behavior in the critical current, in particular, the current can be
dissipationless in one direction and dissipative in the opposite direction. In
the absence of both time-reversal and inversion symmetry, Josephson junctions
(JJ) can exhibit the JDE, characterized by non-reciprocal behavior in the
critical supercurrent ($|I_{c}|$) [16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
26]. This means that the critical supercurrent in the “$+$” direction
($|I^{+}_{c}|$) differs from that in the opposite “$-$” direction
($|I^{-}_{c}|$). This nonreciprocity in supercurrents holds promise for
applications in superconducting electronics [27, 28].

The $d$-wave superconductivity is very common and appears in cuprate
superconductors. The bulk of traditional $d$-wave superconductors is gapless.
Recently, chiral $d$-wave superconductivity has been proposed in magic angle
bilayer twisted graphene [29, 30], in Bi/Ni bilayers [31], and in Sn/Si(111)
[31]. It has been predicted that twisted bilayer cuprates can realize $d+is$
or $d+id^{\prime}$ pairings without gapless states in the bulk [32, 33];
however, further research is necessary to confirm these predictions [34, 35].
With the $d+id^{\prime}$ pairing, the bulk states are gapped, and instead,
chiral edge states appear. Proximity induced superconductivity has been
suggested in two-dimensional electron gas (2DEG) materials such as Bi2O2Se
proximized with cuprates [36]. Furthermore, it has been demonstrated that a
twisted bilayer cuprate exhibits time-reversal symmetry (TRS) breaking and JDE
at specific twist angles [37, 38, 34].

In this work, we study superconductor/normal region/superconductor (S$|$N$|$S)
Josephson junction formed with $d$-wave superconductors. We consider $d$,
$d+id^{\prime}$, and $d+is$ superconducting pairings that have been suggested
in exfoliated layers of cuprate superconductors. We find that for the case of
$d$-wave pairing, the finite JDE can only occur in the presence of spin-orbit
coupling (SOC) and magnetic field. The $d+id^{\prime}$ and $d+is$ pairings
break TRS. As a result, in the presence of additional asymmetry (see Fig. 1)
the TRS breaking results in JDE even in the absence of SOC and magnetic field.
According to our results for $d+id^{\prime}$ pairing, the chiral edge modes
show a clear contribution to JDE, which suggests that interference effects can
be used to control JDE. We further confirm this by studying the distribution
of the Josephson current in the junction and by identifying the Andreev bound
states (ABS) from the local density of states (LDOS). Our proposal realizes a
large, field-free JDE and demonstrates that JDE in a planar geometry can be
used as a signature of unconventional superconducting pairings.

Figure 1: (a), (b) The schematics of the JJ used and the pairing lobes shown
for $d+id^{\prime}$ and $d+is$ superconductor pairing in the leads,
respectively. The shaded area of the normal region is the gated portion. For
both $d+id^{\prime}$ and $d+is$ pairing, the left and right planes are rotated
so that the pairing lobes have different angles, which breaks the necessary
symmetries. (c) parameters used are: $\Delta_{0}=0.01t$,
$\Delta_{0}^{\prime}=0.001t$, $k_{B}T=0.01\Delta_{0}$,
$L_{y}=20\xi\leavevmode\nobreak\ (d+id^{\prime}),10\xi\leavevmode\nobreak\
(d+is)$, $W_{J}=\xi/2$, $\mu=0.1\leavevmode\nobreak\
(d+id^{\prime}),0.2t\leavevmode\nobreak\ (d+is)$.
$I_{0}=e\Delta_{0}L_{y}/\hbar\xi$. $\theta_{L}=-\theta_{R}=0.15\pi$ for
$d+id^{\prime}$ and $\theta_{L}=0\pi,\theta_{R}=0.5\pi$ for $d+is$ pairing.
(d) The schematic of I-V plot for the JJ. Upward (downward) sweep measurement
determines $I_{c}^{+}$ ($I_{c}^{-}$).

##  II Model and methods

We consider S$|$N$|$S Josephson junction (JJ) formed using 2DEG with induced
$d$-wave, $d+id^{\prime}$, or $d+is$ pairing, see Fig. 1. An electrical gate
$V_{g}$ is applied to the normal region as shown in Fig. 1, which allows for
the tunability of JJ.

The Bogoliubov–de Gennes (BdG) Hamiltonian of the system, written in the Nambu
basis
($\psi_{\uparrow},\psi_{\downarrow},\psi^{\dagger}_{\downarrow},-\psi^{\dagger}_{\uparrow}$)
in the continuous limit reads [39]

$\displaystyle H$
$\displaystyle=\left(\frac{\hbar^{2}\bm{k}^{2}}{2m^{*}}-(\mu+V)+\alpha_{SO}\left(\hat{\bm{\sigma}}\times\bm{k}\right)_{z}\right)\hat{\tau}_{z}$
$\displaystyle+\bm{h}.\bm{\sigma}+\hat{\tau}_{x}\text{Re}\Delta(\bm{k},x)-\hat{\tau}_{y}\text{Im}\Delta(\bm{k},x).$
(1)

Here, $m^{*}$ is the effective electron mass, $\mu$ is the chemical potential,
$\alpha$ is the Rashba spin-orbit coupling strength, and $\bm{h}$ describes
the Zeeman energy, e.g., for an external magnetic field. The gate voltage is
defined as $V=V_{G}\;\Theta(W_{G}/2-|x|)$. The pair potential
$\Delta(\bm{k},x)$ is parametrized for $d+id^{\prime}$ pairing as

$\displaystyle\Delta(\bm{k},x)$
$\displaystyle=2e^{i\phi(x)}\Bigl{[}\tilde{\Delta}(x)[\cos{2\theta}\left(\cos{k_{x}}-\cos{k_{y}}\right)$
$\displaystyle+\sin{2\theta}\sin{k_{x}}\sin{k_{y}}]+$ $\displaystyle
i\tilde{\Delta}^{\prime}(x)[\cos{(2\theta+\pi/2)}\left(\cos{k_{x}}-\cos{k_{y}}\right)$
$\displaystyle+\sin{(2\theta+\pi/2)}\sin{k_{x}}\sin{k_{y}}]\Bigr{]},$ (2)

and for $d+is$ pairing as

$\displaystyle\Delta(\bm{k},x)$
$\displaystyle=2e^{i\phi(x)}\Bigl{[}\tilde{\Delta}(x)[\cos{2\theta}\left(\cos{k_{x}}-\cos{k_{y}}\right)$
$\displaystyle+\sin{2\theta}\sin{k_{x}}\sin{k_{y}}]-i\tilde{\Delta}^{\prime}[\cos{k_{x}}+\cos{k_{y}}]\Bigr{]},$
(3)

where $\theta$ describes the angle between the $x$-axis and the lobe direction
of the pair potential. We denote $\theta_{L}$ ($\theta_{R}$) for the left
(right) superconductor. The terms
$\tilde{\Delta}(x),\tilde{\Delta}^{\prime}(x)$ are only nonzero in the regions
covered by superconductors where $\tilde{\Delta}(x)=\Delta_{0}$,
$\tilde{\Delta}^{\prime}(x)=\Delta^{\prime}_{0}$. For the region covered by
the left superconductor, we take $\phi(x)=\phi_{L}$, and for the region
covered by the right superconductor, we take $\phi(x)=\phi_{R}$. We discretize
the Hamiltonian (1) and use the tight binding approach where
$t=\frac{\hbar^{2}}{2m^{*}a^{2}}$ is the nearest neighbor hopping and $a$ is
the lattice spacing. We use dimensionless units measuring energy in units of
$t$ and length in units of $a$. The SC coherence length is defined as
$\xi=\frac{\hbar^{2}k_{F}}{m^{*}\Delta_{0}}=2a{\frac{\sqrt{\mu
t}}{\Delta_{0}}}$. Typical material parameters we use in numerical
calculations are as follows: $a=2.5nm$, $m^{*}=0.14m_{e}$,
$\alpha_{SO}=0.1ta=0.1eV{\textup{\AA}}$, $\Delta_{0}=0.01t=0.4meV$, and
$\mu=0.1t=4meV$ [40].

We calculate the conserved current inside the junction (see Fig. 8(a)) which
does not require the self-consistent calculation. The self-consistent
procedure would be required inside superconductor to ensure the current
conservation. To calculate the equilibrium Josephson current, we use the
Matsubara Green’s function $G(i\omega_{n})$ obtained using the analytical
continuation of the retarded Green’s function. The local current between the
sites $\mathbf{n}$ and $\mathbf{m}$ can be expressed using the standard
approach [41, 42, 43, 44, 45, 46]:

$\displaystyle
I_{\mathbf{n},\mathbf{m}}=2\frac{ek_{B}T}{\hbar}\sum_{\mathbf{\omega_{\mathbf{n}}}}\mathrm{Tr}\left[Im\left(H_{\mathbf{n}\mathbf{m}}G_{\mathbf{m}\mathbf{n}}-H_{\mathbf{m}\mathbf{n}}G_{\mathbf{n}\mathbf{m}}\right)\right],$
(4)

where $G_{\mathbf{n}\mathbf{m}}$ is the Green’s function and
$H_{\mathbf{n}\mathbf{m}}$ is the Hamiltonian submatrix calculated between the
sites $\mathbf{n}$ and $\mathbf{m}$, with the lattice sites taken inside the
normal region. The summation is performed over the fermionic Matsubara
frequencies $\omega_{n}=(2n+1)\pi k_{B}T$ and the Green’s functions
$G_{\mathbf{n}\mathbf{m}}$ also accounts for the self energy of the leads. The
total current is obtained by summing all contributions along the cut in
$y$-direction, $I=\sum_{\mathbf{n}}I_{\mathbf{n},\mathbf{n}+\mathbf{e}_{x}}$,
where $\mathbf{e}_{x}$ is the unit vector along the $x$-axis (see Fig. 8(a)).
The diode effect is then characterized by the quality factor
$Q=(I^{+}_{c}-|I^{-}_{c}|)/(I^{+}_{c}+|I^{-}_{c}|)$, where $I^{+}_{c}$ and
$I^{-}_{c}$ are the critical current calculated by finding the maximum and
minimum of current over $\phi=\phi_{L}-\phi_{R}\in\\{0,2\pi\\}$. As an
example, in Fig. 1(c) we show the current-phase relation (CPR) for
$d+id^{\prime}$ and $d+is$ pairings. The figure shows the nonvanishing
anomalous Josephson effect at $\phi=0$ provided that the two pairing lobes are
rotated with respect to each other. Furthermore, CPR in the figure results in
the quality factor $Q$ equal to $0.2$ for $d+id^{\prime}$ pairing and $0.14$
for $d+is$ pairing.

We also calculate the local density of states (LDOS) using the expression:

$\displaystyle\rho(\omega,n)=-\frac{1}{\pi}\mathrm{Tr_{n}}(\mathrm{Im}[G(\omega+i\eta)]),$
(5)

where $i\eta$ is the band broadening. The location of LDOS peaks allows us to
determine energies of the Andreev bound states (ABS).

SC Pairing | $\mathcal{R}_{x}(\pi)$ | $\mathcal{R}_{y}(\pi)$ | $\mathcal{T}$ | Symmetries to break or
---|---|---|---|---
|  |  |  | terms to add for JDE
$d_{x^{2}-y^{2}}$ | $\checkmark$ | $\checkmark$ | $\checkmark$ |  $h\hat{\sigma}_{y}$ \+ Rashba SOC
$d_{xy}$ | $\checkmark^{\star}$ | $\checkmark^{\star}$ | $\checkmark$ |  $h\hat{\sigma}_{y}$ \+ Rashba SOC
$d_{x^{2}-y^{2}}+d_{xy}$ | $\times$ | $\times$ | $\checkmark$ |  $h\hat{\sigma}_{y}$ \+ Rashba SOC
$d_{x^{2}-y^{2}}+id_{xy}$ | $\times$ | $\times$ | $\times$ | $\mathcal{R}_{x}(\pi)\mathcal{R}_{y}(\pi)\&\mathcal{R}_{x}(\pi)\mathcal{T}$
$d_{x^{2}-y^{2}}+is$ | $\checkmark$ | $\checkmark$ | $\times$ | $\mathcal{R}_{y}(\pi)$
Table 1: The symmetry of each SC pairing type and corresponding JJ needed to
have JDE.⋆ With a $\pi$ gauge transformation of SC pairing. Figure 2: (a) The
pictorial symmetry analysis of the JJ with $d+id$ pairing in the leads.
$I_{c}$ is the current inside the normal region, shown by the black arrow. To
recover the correct chirality and phase sign, a rotation around y-axis
followed by a rotation around x-axis is needed (inversion with respect to the
center of JJ). (b) With $d+is$ pairing in the leads, $\mathcal{R}_{x}(\pi)$
would recover $I(\phi)=-I(-\phi)$ only if left and right side satisfy
$\beta_{L}=\beta_{R}$. (c) For $d_{x^{2}-y^{2}}$ pairing in the leads, we need
non-zero magnetic and SOC terms. The SOC term prevents the spin rotation to
recover the correct spin-momentum coupling.

##  III Symmetry analysis

To understand how the anomalous Josephson effect and JDE can be achieved, we
begin with a symmetry analysis of the system, followed by numerical
simulations of the diode effect while varying various parameters. In order to
get a non-reciprocal transport in a Josephson junction, it is required to have
$I(\phi)\neq-I(-\phi)$. If however there is a symmetry $\tilde{T}$ that
reverses $I$ and $\phi$ then the anomalous Josephson effect and JDE (by
extension of the above argument to applied voltage $V$) vanish. To identify
such a symmetry, we start with analysis of systems with $d$-wave pairings. In
Table 1, we list whether TRS $\mathcal{T}$ and $\pi$-rotation symmetries with
respect to $x$ and $y$ axis ($\mathcal{R}_{x}(\pi),\mathcal{R}_{y}(\pi)$) are
satisfied for the system in Fig. 1 for simple alignment of the pairing lobes.
In the last column, we list modifications that are necessary to break the
symmetry or symmetries prohibiting the anomalous Josephson effect. For
example, in case of $d_{x^{2}-y^{2}}+d_{xy}$ pairing, TRS is present, and
$\mathcal{T}e^{i\phi}\Delta\mathcal{T}^{-1}=e^{-i\phi}\Delta$. This means that
to achieve the nonreciprocity, we need to break both TRS and spin rotational
symmetry. The choice of the direction of the magnetic field and the presence
of Rashba SOC guarantee that we cannot use $\mathcal{R}_{x}(\pi)$ and
$\mathcal{R}_{y}(\pi)$ symmetries for the transformation that reverses both
$I$ and $\phi$. With $d_{x^{2}-y^{2}}+id_{xy}$ pairing, all of the three
symmetries ($\mathcal{T}$, $\mathcal{R}_{x}(\pi)$, $\mathcal{R}_{y}(\pi)$) are
separately broken. However, the combined symmetries
$\mathcal{R}_{x}(\pi)\mathcal{R}_{y}(\pi)$ and
$\mathcal{R}_{x}(\pi)\mathcal{T}$ still need to be broken for the presence of
the anomalous Josephson effect. This can be achieved by using JJ with an
asymmetry (e.g. induced by gating) or by rotating the pairing lobes of the
left and right superconductors with respect to each other. Finally, for the
$d_{x^{2}-y^{2}}+is$ pairing, the rotational symmetries are present. Since
$\mathcal{R}_{y}(\pi)$ reverses both $I$ and $\phi$ it has to be broken for
the presence of the anomalous Josephson effect. As before, this can be
achieved by using JJ with an asymmetry (e.g. induced by gating) or by rotating
the pairing lobes of the left and right superconductors with respect to each
other. The latter approach is explored in details in our numerical
simulations. Figure 2 shows the pictorial demonstration of the symmetry
analysis. In particular, Fig. 2(a) shows that for $d+id^{\prime}$ pairing the
orientations of the pairing lobes in the left and right SCs (i.e.
$\theta_{L}\neq\theta_{R}$) can break
$\mathcal{R}_{x}(\pi)\mathcal{R}_{y}(\pi)$ and
$\mathcal{R}_{x}(\pi)\mathcal{T}$ symmetries. In Fig. 2(b), we show that for
$d+is$ pairing, the orientations of the pairing lobes (i.e.
$\theta_{L}\neq\theta_{R}$) break $\mathcal{R}_{y}$ symmetry provided that
$\theta_{L}\neq-\theta_{R}$. In Fig. 2(c), we assume $d_{x^{2}-y^{2}}+d_{xy}$
pairing. Here, a Zeeman field along the $y$ axis is necessary to break TRS. To
further break the spin rotation symmetries, the Rashba SOC term has to be
added. It is worth mentioning that the anomalous Josephson effect has been
predicted for $s$-wave pairing without the presence of SOC by introducing
appropriate magnetic textures [47].

##  IV Leads with $d$-wave pairing

Here, we consider a ballistic JJ with superconducting leads featuring the pure
$d$-wave pairing. To treat systems with large $L_{y}$, we assume here periodic
boundaries along the $y$ axis and numerically find solutions characterized by
the transverse wave vector $k_{y}$. As a result, we solve numerically the
effective one-dimensional scattering problem along the direction of the
junction, afterwards summing the contributions from each $k_{y}$. With this
approach, we are able to reproduce analytical results in Ref. [48] where the
Zeeman and Rashba terms were only present in the junction region.

For an arbitrary orientation of the pairing lobes, the pairing potential can
take a general form of $d_{x^{2}-y^{2}}+d_{xy}$. For nonvanishing JDE, as
shown before, we need the Zeeman term to break the time-reversal symmetry, and
the Rashba SOC to break the spin rotation symmetry. Figure 3 shows the quality
factor $Q$ plotted as a function of the pairing lobe orientation angles
$\theta_{L}$ and $\theta_{R}$. Variations in the lobe orientation can be
achieved by adjusting the crystal orientation of SC leads. In Fig. 3(a), the
Zeeman field is applied in the $x$ direction. The figure shows regions of
large $Q$ throughout the phase space. In Fig. 3(b), the Zeeman field is
applied in the $y$ direction. For this direction of the Zeeman field, even
with no rotation of lobes, i.e., for $\theta_{L}=0$ and $\theta_{R}=0$, we
observe large values of $Q$. Finally, we observe relatively weak JDE for the
Zeeman field along the $z$ axis (not shown in the figure) in contrast to a JJ
based on the topological insulator [9].

Figure 3: Q is plotted as $\theta_{L}$ and $\theta_{R}$ is varied. The Zeeman
energy is in the x, y direction for (a), (b), respectively. Parameters used:
$\mu=0.1\>t$, $\alpha_{SO}=0.1\>ta$, $\Delta_{0}=0.01\>t$,
$\Delta^{\prime}_{0}=0\>t$, $k_{B}T=0.01\>\Delta_{0}$, $W_{J}=\xi/2$,
$|h|=0.1\>\Delta_{0}$. (a) $h=(h_{x},0,0)$. At the off diagonal values
($\theta_{L}=-\theta_{R}$) Q is zero. (b) $h=(0,h_{y},0)$. The calculations
are performed using a Fourier transform along the y-axis.  Figure 4: (a) The
plot of Q as $W_{J}$ is increased for two cases of Zeeman energy in the x and
y directions. (b) The plot of Q as the gate potential $V_{G}$ is varied.
Parameters used: $\mu=0.1\>t$, $\alpha_{SO}=0.1\>ta$, $\Delta_{0}=0.01\>t$,
$\Delta^{\prime}_{0}=0$, $k_{B}T=0.01\>\Delta_{0}$,
$\theta_{L}=\theta_{R}=-0.15\>\pi$ for $h=(h_{x},0,0)$ and
$\theta_{L}=-\theta_{R}=-0.10\>\pi$ for $h=(0,h_{y},0)$ plots.
$h_{x}=0.2\>\Delta_{0}$, $h_{y}=0.1\>\Delta_{0}$ is used for (a) and (b)
plots. $W_{J}=\xi/2$ in (b) and (c) plots. $W_{G}=\xi/8$ in (b). (c) The plot
of Q as $|h|$ is varied in $x$ and $y$ directions. Figure 5: $I_{max}$ and
$I_{min}$ are plotted as $W_{J}$ is varied for h in $x$ and $y$ direction in
(a) and (b), respectively. (e) and (f) are plots of $I_{max}$ and $I_{min}$ as
$|h|$ is varied. The parameters used are the same as Fig. 4.

Several parameters can be adjusted to fine-tune the quality factor of JJ.
Fist, we vary the junction width $W_{J}$. As shown in Fig. 4(a), the quality
factor shows the large variation and even sign change as $W_{J}$ is increased.
This behavior stems from contributions of the higher order harmonics to the
current-phase relation (CPR). Next, we also explore how $Q$ depends on the
gate voltage $V_{G}$ and the Zeeman field $\mathbf{h}$. Figure 4(b)
demonstrates that $Q$ can be effectively tuned by the gate voltage within the
junction region for the Zeeman field along the $y$ axis. The variation of $Q$
as a function of the Zeeman field along the $x$ and $y$ directions is shown in
Fig. 4(c). In Fig. 4(c), we observe large variations as the Zeeman field is
increased, with eventual suppression of JDE at large Zeeman fields. Even
though the dependence of the quality factor $Q$ on various parameters is
complicated, some trends can be identified. This can be useful in fine tuning
the diode effect for applications such as Josephson transistor.

This is further elucidated in Fig. 5 showing the evolution of the maximum,
$I_{max}$, and minimum, $I_{min}$, Josephson currents. Figure 5 shows the
evolution of $I_{max}$ and $I_{min}$ for the calculations in Fig. 4, for the
Zeeman field $h$ along the $x$ axis in (a), (c), and (e) plots, and for the
Zeeman field $h$ along the $y$ axis in (b), (d), and (f) plots. The high
tunability of Josephson current is seen in Fig. 5(b) where a dip in $I_{min}$
corresponds to the peak observed in $Q$. In comparison, Fig. 5(a) shows a
monotonic decay in both $I_{max}$ and $I_{min}$ corresponding to the smoother
change in $Q$, as can be seen in Fig. 4(a). In the long junction limit
$W_{J}\gg\xi$, we observe that $Q$ tends to zero. By looking at the variation
of $I_{max}$ and $I_{min}$ as the gate voltage $V_{G}$ is varied in Fig. 5(c),
we observe only minor change. However, Fig. 5(d) demonstrates high tunability
by gate voltage as we observe a dip in $I_{min}$ and a peak in $I_{max}$ at
around $V_{G}=0$. In this case, the gating also substantially modifies $Q$,
even leading to a sign change. Figures 5(e) and (f) show large $I_{max}$ and
$I_{main}$ at small $h$. Away from the region $h=0$ the magnetic field
suppresses the supercurrent. We observe that the Zeeman field skews $I_{max}$
and $I_{min}$ in the opposite directions, which correspond to the sign change
of $Q$ seen in Fig. 4(c). Overall, Figs. 4 and 5 demonstrate high tunability
of JDE for leads with the $d$-wave pairing.

##  V Leads with $d\pm id^{\prime}$ pairing

Now we consider a ballistic JJ with superconducting leads featuring the $d\pm
id^{\prime}$ pairing. For this pairing, the TRS is broken. To satisfy the
conditions of JDE, we still need to break rotational symmetries (see Table 1).
This can be achieved by rotating the pairing lobes parametrized by angle
$\theta$ in Eq. (II). Figure 6(a) shows the quality factor $Q$ as a function
of $\theta_{L}$ and $\theta_{R}$ without any Zeeman field or SOC for a
$d+id^{\prime}|$N$|d+id^{\prime}$ JJ. We observe a possibility of large JDE
even in the absence of the Zeeman field and SOC. As the system may have the
Rashba SOC, we also repeated the calculations with the inclusion of Rashba
SOC, finding only minor differences. Figure 6(b) shows the quality factor $Q$
as a function of $\theta_{L}$ and $\theta_{R}$ without any Zeeman field or SOC
for a $d-id^{\prime}|$N$|d+id^{\prime}$ JJ. Compared to the previous case, we
observe significantly smaller values of $Q$. For this case no JDE is present
on the off diagonal, i.e., $\theta_{L}=-\theta_{R}$. We further look at how Q
evolves as $W_{J}$ is varied. Figure 7(a) shows $Q$ for the case of a
$d+id^{\prime}|$N$|d+id^{\prime}$ JJ. The short junction limit $W_{J}\ll\xi$
shows the largest $Q$, although the decrease in $Q$ as $W_{J}$ is increased is
not significant. Figure 7(b) is the calculations for a
$d-id^{\prime}|$N$|d+id^{\prime}$ JJ. As before, we observe much smaller $Q$
for superconductors with opposite chiralities.

In the case of $d+id^{\prime}$ pairing, the superconductors exhibit chiral
edge states, with the relative sign of the $d$ and $d^{\prime}$ components
determining the chirality of the system. To investigate how edge currents
contribute to JDE, we examine the current distribution along the transverse
direction for a JJ of finite width $L_{y}$. Figure 8(a) illustrates the
schematics of the current distribution calculation. We consider a
$d+id^{\prime}|$N$|d+id^{\prime}$ JJ in Figs. 8(b) and (c). We observe a
concentration of Josephson current at the top and bottom edges, consistent
with the effect of edge states. As expected, in Fig. 8(c) we observe vanishing
JDE. In Fig. 8(b), the current in the bulk of the channel exhibits a
significantly stronger anomalous Josephson effect compared to edges, and the
edge part of the current mostly follows a $\sin{\phi}$ CPR behavior resulting
in relatively weak contribution to JDE. Figures 8(d) and (e) show the current
distribution for a $d-id^{\prime}|$N$|d+id^{\prime}$ JJ. In Fig. 8(e), we
observe vanishing JDE. In Fig. 8(d), both the bulk and the edge modes seem to
contribute to JDE; however, overall JDE is much weaker in this case. If we
compare the current distribution calculations for $d\pm id^{\prime}$ in Figs.
8(b–e) with similar calculations for $d\pm is$ case in Figs. 8(f, g), we
observe a clear difference that can be associated with the edge states present
in $d\pm id^{\prime}$ systems. Figures 8(f, g) will be discussed in more
detail in the next section.

Figure 6: Q phase space as $\theta_{L}$, $\theta_{R}$ are varied for a
$d+id^{\prime}|$N$|d+id^{\prime}$ JJ in (a) and a
$d+id^{\prime}|$N$|d+id^{\prime}$ JJ in (b). (a) Along the
$\theta_{R}=\theta_{L}$ line, Q is zero. (b) In this case, on the off-diagonal
line ($\theta_{R}=-\theta_{L}$), Q is zero. Parameters used: $\mu=0.1\>t$,
$\Delta_{0}=0.01\>t$, $\Delta^{\prime}_{0}=0.001\>t$,
$k_{B}T=0.01\>\Delta_{0}$, and $W_{J}=\xi/2$. The calculations are performed
using a Fourier transform along the y-axis.  Figure 7: The plot of Q as
$W_{J}$ is increased for various spin-orbit coupling strength $\alpha_{SO}$
for the two cases of the same chirality ($d+id^{\prime}|$N$|d+id^{\prime}$)
and opposite chirality ($d-id^{\prime}|$N$|d+id^{\prime}$) in (a) and (b),
respectively. Parameters used: $\mu=0.1\>t$, $\Delta_{0}=0.01\>t$,
$\Delta^{\prime}_{0}=0.001\>t$, $k_{B}T=0.01\>\Delta_{0}$.
$\theta_{L}=-\theta_{R}=-0.15\>\pi$ used for (a) and $\theta_{L}=0$,
$\theta_{R}=0.2\>\pi$ used for (b).  Figure 8: (a) Schematics of the tight
binding of JJ with SC leads at the left and right sides. The current
distribution along y as $\phi$ is varied for a
$d+id^{\prime}|$N$|d+id^{\prime}$ JJ with $\theta_{L}=-\theta_{R}=0.15\>\pi$
in (b) and $\theta_{L}=\theta_{R}=0.15\>\pi$ in (c). The JJ in (d) and (e) is
the $d-id^{\prime}|$N$|d+id^{\prime}$ type. $\theta_{L}=0$,
$\theta_{R}=0.2\>\pi$ is used in (d) and $\theta_{L}=-\theta_{R}=0.15\>\pi$
used in (e). Current density - phase relation with $d+is$ pairing shown in (f)
with $\theta_{L}=0$, $\theta_{R}=0.5\>\pi$ and in (g) with
$\theta_{L}=\theta_{R}=0.15\>\pi$. The dark red of the colormap corresponds to
the minimum of the current distribution, and the dark blue corresponds to the
maximum of the current distribution in each plot. Parameters used:
$\mu=0.1\>t$, $\Delta_{0}=0.01\>t$, $\Delta^{\prime}_{0}=0.001\>t$,
$k_{B}T=0.01\>\Delta_{0}$, $W_{J}=\xi/2$.

Next, we look at the LDOS to identify the behavior of ABS in the vicinity of
edges. Fig. 9(a) shows the LDOS at the top edge of the channel for a system
with leads featuring $d+id^{\prime}$ pairing. We observe a clear asymmetry as
$\phi$ is changed which signifies the nonreciprocity of the current density.
In Fig. 9(b), we show the LDOS calculated in the middle of the channel, which
also reveals the asymmetry and nonreciprocity. However, the pattern of the
LDOS is clearly different compared to the LDOS at the top edge. This agrees
with the behavior of the current density in Fig. 8 attributed to the presence
of edge states. The LDOS at the bottom edge in Fig. 9(c) resembles the LDOS of
the top edge. Above behavior is in contrast with Fig. 10 calculated for $d+is$
pairing where we observe similar patterns for the top, bottom, and middle
parts, and no signatures of edge states.

Figure 9: LDOS at the (a) top edge, (b) middle, and (c) bottom edge of the
channel with $d+id^{\prime}$ pairing in the leads. The peaks of LDOS
correspond to the ABS energies. Most of the bound states are present at both
the top and bottom. However, a few bound states only appear at one end. The
bound states in the middle of the channel are significantly different from the
two edges. $\mu=0.1\>t$, $\Delta_{0}=0.01\>t$, $\Delta^{\prime}_{0}=0.001\>t$,
$\theta_{L}=-\theta_{R}=0.15\>\pi$, $L_{y}=40\>\xi$ and $W_{J}=\xi/2$.  Figure
10: LDOS at the (a) top edge, (b) middle, and (c) bottom edge of the channel
with $d+is$ pairing in the leads. Parameters used: $\mu=0.1\>t$,
$\Delta_{0}=0.01\>t$, $\Delta^{\prime}_{0}=0.001\>t$, $\theta_{L}=0$,
$\theta_{R}=0.5\>\pi$, $L_{y}=20\>\xi$ and $W_{J}=\xi/2$.

To separate the edge contributions from the bulk contributions, we look at the
behavior of $Q$ as the transverse size $L_{y}$ is increased. Note that to make
the calculation feasible, we had to choose unrealistically large value of the
pairing potential $\Delta_{0}$. Figure 11(a) shows that in case of
$d+id^{\prime}$ pairing, a JJ with a relatively narrow channel
($L_{y}\sim\xi$) exhibits interplay of the edge and bulk contributions where
$Q$ reaches a maximum. As $L_{y}$ is increased further, the bulk current
contribution becomes dominant, which results in the sign reversal of $Q$. We
eventually get close to the asymptotic value corresponding to
$L_{y}\rightarrow\infty$ where $Q$ is calculated using the Fourier transform
approach that ignores the edge contributions. In comparison, for a similar
calculation with leads featuring the $d+is$ pairing, the asymptotic limit of
the channel width is reached very quickly ($L_{y}\approx 2\xi$), as shown in
Fig. 11(b).

Figure 11: Q is plotted as $L_{y}$ is increased for $d+id$ (a) and
$d+is^{\prime}$ (b) pairing in the leads. The peak in Q for $d+id^{\prime}$
pairing corresponds to the edge-dominated transport. Two sets of parameters
used: $\mu=0.1\>t$, $\Delta_{0}=0.1\>t$, $\Delta^{\prime}_{0}=0.01\>t$, and
$\mu=0.2\>t$, $\Delta_{0}=0.1\>t$, $\Delta^{\prime}_{0}=0.01\>t$.
$W_{J}=\xi/2$ for both sets.

Similar to the JJ considered in the previous section, the systems featuring
$d+id^{\prime}$ pairing can be tuned by a gate potential in the normal region.
This is demonstrated in Fig. 13 where $Q$ is calculated as a function of the
gate potential $V_{G}$. The strong dependence on the gate voltage is observed
for $\theta_{L}$ and $\theta_{R}$ that are in the vicinity of points where $Q$
changes sign in Fig. 6(a), as can be seen by comparing different curves in
Figs. 13(a) and (b). We can achieve tunability and even the sign change, thus
replicating the transistor functionality.

Figure 12: (a) Phase diagram of Q as $\theta_{L},\theta_{R}$ are varied in a
$d+is|$N$|d+is$ JJ. Along the $\theta_{R}=\pm\theta_{L}$ lines, Q is zero. A
sign change also occurs at these lines. (b) Plot of Q for various
$\alpha_{SO}$ as $W_{J}$ is increased with $\theta_{L}=0$,
$\theta_{R}=0.5\>\pi$. Parameters used: $\mu=0.1\>t$, $\Delta_{0}=0.01\>t$,
$\Delta^{\prime}_{0}=0.001\>t$, $k_{B}T=0.01\>\Delta_{0}$, and $W_{J}=\xi/2$
for (a). The calculations are performed using a Fourier transform along the
y-axis

##  VI Leads with $d\pm is$ pairing

Finally, we consider a ballistic JJ with superconducting leads featuring the
$d\pm is$ pairing. Figure 12 (a) shows the quality factor $Q$ as a function of
the orientation angles $\theta_{L}$ and $\theta_{R}$. For a $d+is|$N$|d+is$
JJ, we observe in Fig. 12(a) that there is no JDE for orientations described
by $\theta_{L}=\theta_{R}$ or $\theta_{L}=-\theta_{R}$, which is consistent
with our symmetry analysis. The vanishing JDE for orientations described by
$\theta_{L}=-\theta_{R}$ is in contrast with the case of $d+id^{\prime}$
pairing in Fig. 6(a), and can be used as a signature to distinguish between
the two types of pairings experimentally. The quality factror a for
$d-is|$N$|d+is$ JJ can be obtained from Fig. 12(a) by a transformation
$\theta_{L}\rightarrow\theta_{L}+\pi/2$ or
$\theta_{R}\rightarrow\theta_{R}+\pi/2$. This follows from Eq. (II) as the
transformation $\theta\rightarrow\theta+\pi/2$ changes $d-is$ to $-d-is$. In
Fig. 6(b), we study the variation of $Q$ with the junction length $W_{J}$ and
find that the quality factor decreases by a factor of two for a longer
junction. Furthermore, the presence of SOC only leads to minor changes. In
Figs 8(f, g), we examine the current distribution along the transverse
direction for a JJ of finite width $L_{y}$. In Fig. 8(f) for a $d+is|$N$|d+is$
JJ, we observe a uniform distribution of current along the $y$ direction for
orientation angles $\theta_{L}=0$ and $\theta_{R}=0.5\pi$. In Fig. 8(g) for
orientation angles $\theta_{L}=\theta_{R}=0.15\pi$, we observe only a minor
difference between the edge and the bulk and a vanishing diode effect. To
further identify the differences between $d+id^{\prime}$ and $d+is$ pairings,
we study the width dependence of JDE in Fig. 11. For the $d+is$ case in Fig.
11(b), we observe that JDE quickly reaches the asymptotic value, which is
consistent with the absence of edge contributions for $d\pm is$ pairing. In
Fig. 13(b), we study how a JJ with $d+is$ pairing can be tuned by a gate
potential in the normal region. As in the previously considered cases, we
observe that the qualuity factor $Q$ can be tuned by the gate voltage $V_{G}$,
especially for orientation angles $\theta_{L}$ and $\theta_{R}$ in Fig. 12(b)
corresponding to the sign change of $Q$.

Figure 13: The plot of Q as $V_{G}$ is varied for $d+id^{\prime}$ (a) and
$d+is$ (b) pairing in the leads. Parameters used: $\mu=0.1\>t$,
$\Delta_{0}=0.01\>t$, $\Delta^{\prime}_{0}=0.001\>t$,
$k_{B}T=0.01\>\Delta_{0}$, $W_{G}=\xi/4$, and $W_{J}=\xi/2$.

##  VII Conclusions

In this work, we have studied JDE in S$|$N$|$S Josephson junction formed using
2DEG and superconductors with $d$, $d+id^{\prime}$, and $d+is$ superconducting
pairings. We have shown that the quality factor and its sign $Q$ can be
substantially tuned by the external magnetic field, gate voltage, and the
length of the junction for all three types of pairings. We have demonstrated
that due to the intrinsic TRS breaking in $d+is$ and $d+id^{\prime}$
superconductors, we can realize field-free JDE, where even a relatively small
TRS-breaking superconducting pairing component can result in a large JDE. In
particular, by breaking certain rotational symmetries a large JDE is possible
in the absence of the Zeeman field and SOC. By calculating LDOS and the
current distribution inside the normal region, we have also studied the role
of edge states associated with $d+id^{\prime}$ pairing. We see clear
signatures of the edge states in systems with $d+id^{\prime}$ pairing where
the edge currents can have a significant contributions to JDE. We find that
JDE vanishes in systems with $d+is$ and $d+id^{\prime}$ pairings under
different symmetry constraints. This leads us to the conclusion that JDE can
be used to differentiate between different superconducting pairings. Our
proposal can be useful for superconductor rectification and realizations of
transistor-like functionalities, e.g., also using magnetic textures [49, 50],
in superconducting devices operating at low temperatures.

###### Acknowledgements.

This work was supported by the National Science Foundation/EPSCoR RII Track-1:
Emergent Quantum Materials and Technologies (EQUATE), Award OIA-2044049.

## References

  * Ando  _et al._ [2020] F. Ando, Y. Miyasaka, T. Li, J. Ishizuka, T. Arakawa, Y. Shiota, T. Moriyama, Y. Yanase, and T. Ono, Nature 584, 373–376 (2020).
  * Miyasaka  _et al._ [2021] Y. Miyasaka, R. Kawarazaki, H. Narita, F. Ando, Y. Ikeda, R. Hisatomi, A. Daido, Y. Shiota, T. Moriyama, Y. Yanase, and T. Ono, Appl. Phys. Express 14, 073003 (2021).
  * Zhang  _et al._ [2020] E. Zhang, X. Xu, Y.-C. Zou, L. Ai, X. Dong, C. Huang, P. Leng, S. Liu, Y. Zhang, Z. Jia, X. Peng, M. Zhao, Y. Yang, Z. Li, H. Guo, S. J. Haigh, N. Nagaosa, J. Shen, and F. Xiu, Nat Commun 11, 5634 (2020).
  * Wakatsuki  _et al._ [2017] R. Wakatsuki, Y. Saito, S. Hoshino, Y. M. Itahashi, T. Ideue, M. Ezawa, Y. Iwasa, and N. Nagaosa, Sci. Adv. 3, e1602390 (2017).
  * Lyu  _et al._ [2021] Y.-Y. Lyu, J. Jiang, Y.-L. Wang, Z.-L. Xiao, S. Dong, Q.-H. Chen, M. V. Milošević, H. Wang, R. Divan, J. E. Pearson, P. Wu, F. M. Peeters, and W.-K. Kwok, Nat Commun 12, 2703 (2021).
  * Narita  _et al._ [2022] H. Narita, J. Ishizuka, R. Kawarazaki, D. Kan, Y. Shiota, T. Moriyama, Y. Shimakawa, A. V. Ognev, A. S. Samardak, Y. Yanase, and T. Ono, Nat. Nanotechnol. 17, 823–828 (2022).
  * Strambini  _et al._ [2022] E. Strambini, M. Spies, N. Ligato, S. Ilić, M. Rouco, C. González-Orellana, M. Ilyn, C. Rogero, F. S. Bergeret, J. S. Moodera, P. Virtanen, T. T. Heikkilä, and F. Giazotto, Nat Commun 13, 2431 (2022).
  * He  _et al._ [2022a] J. J. He, Y. Tanaka, and N. Nagaosa, New J. Phys. 24, 053014 (2022a).
  * Tanaka  _et al._ [2022] Y. Tanaka, B. Lu, and N. Nagaosa, Phys. Rev. B 106, 214524 (2022).
  * Wu  _et al._ [2022] H. Wu, Y. Wang, Y. Xu, P. K. Sivakumar, C. Pasco, U. Filippozzi, S. S. P. Parkin, Y.-J. Zeng, T. McQueen, and M. N. Ali, Nature 604, 653–656 (2022).
  * Nadeem  _et al._ [2023] M. Nadeem, M. S. Fuhrer, and X. Wang, Nat Rev Phys 5, 558–577 (2023).
  * Bocquillon  _et al._ [2016] E. Bocquillon, R. S. Deacon, J. Wiedenmann, P. Leubner, T. M. Klapwijk, C. Brüne, K. Ishibashi, H. Buhmann, and L. W. Molenkamp, Nat. Nanotechnol. 12, 137–143 (2016).
  * Pal  _et al._ [2022] B. Pal, A. Chakraborty, P. K. Sivakumar, M. Davydova, A. K. Gopi, A. K. Pandeya, J. A. Krieger, Y. Zhang, M. Date, S. Ju, N. Yuan, N. B. M. Schröter, L. Fu, and S. S. P. Parkin, Nat. Phys. 18, 1228–1233 (2022).
  * Baumgartner  _et al._ [2022] C. Baumgartner, L. Fuchs, A. Costa, J. Picó-Cortés, S. Reinhardt, S. Gronin, G. C. Gardner, T. Lindemann, M. J. Manfra, P. E. Faria Junior, D. Kochan, J. Fabian, N. Paradiso, and C. Strunk, J. Phys.: Condens. Matter 34, 154005 (2022).
  * Gupta  _et al._ [2023] M. Gupta, G. V. Graziano, M. Pendharkar, J. T. Dong, C. P. Dempsey, C. Palmstrøm, and V. S. Pribiag, Nat Commun 14, 3078 (2023).
  * Zhang  _et al._ [2022] Y. Zhang, Y. Gu, P. Li, J. Hu, and K. Jiang, Phys. Rev. X 12, 041013 (2022).
  * Misaki and Nagaosa [2021] K. Misaki and N. Nagaosa, Phys. Rev. B 103, 245302 (2021).
  * Kopasov  _et al._ [2021] A. A. Kopasov, A. G. Kutlin, and A. S. Mel’nikov, Phys. Rev. B 103, 144520 (2021).
  * He  _et al._ [2022b] J. J. He, Y. Tanaka, and N. Nagaosa, New J. Phys. 24, 053014 (2022b).
  * Daido  _et al._ [2022] A. Daido, Y. Ikeda, and Y. Yanase, Phys. Rev. Lett. 128, 037001 (2022).
  * Yuan and Fu [2022] N. F. Q. Yuan and L. Fu, Proc Natl Acad Sci 119, e2119548119 (2022).
  * Ilić and Bergeret [2022] S. Ilić and F. S. Bergeret, Phys. Rev. Lett. 128, 177001 (2022).
  * Halterman  _et al._ [2022] K. Halterman, M. Alidoust, R. Smith, and S. Starr, Phys. Rev. B 105, 104508 (2022).
  * Karabassov  _et al._ [2022] T. Karabassov, I. V. Bobkova, A. A. Golubov, and A. S. Vasenko, Phys. Rev. B 106, 224509 (2022).
  * Souto  _et al._ [2022] R. S. Souto, M. Leijnse, and C. Schrade, Phys. Rev. Lett. 129, 267702 (2022).
  * Fominov and Mikhailov [2022] Y. V. Fominov and D. S. Mikhailov, Phys. Rev. B 106, 134514 (2022).
  * Braginski [2018] A. I. Braginski, J Supercond Nov Magn 32, 23–44 (2018).
  * Rymarz  _et al._ [2021] M. Rymarz, S. Bosco, A. Ciani, and D. P. DiVincenzo, Phys. Rev. X 11, 011032 (2021).
  * Oh  _et al._ [2021] M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, X. Liu, K. Watanabe, T. Taniguchi, and A. Yazdani, Nature 600, 240–245 (2021).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43–50 (2018).
  * Gong  _et al._ [2017] X. Gong, M. Kargarian, A. Stern, D. Yue, H. Zhou, X. Jin, V. M. Galitski, V. M. Yakovenko, and J. Xia, Sci. Adv. 3, e1602579 (2017).
  * Sigrist [1998] M. Sigrist, Progress of Theoretical Physics 99, 899–929 (1998).
  * Can  _et al._ [2021] O. Can, T. Tummuru, R. P. Day, I. Elfimov, A. Damascelli, and M. Franz, Nat. Phys. 17, 519–524 (2021).
  * Zhao  _et al._ [2023] S. Y. F. Zhao, X. Cui, P. A. Volkov, H. Yoo, S. Lee, J. A. Gardener, A. J. Akey, R. Engelke, Y. Ronen, R. Zhong, G. Gu, S. Plugge, T. Tummuru, M. Kim, M. Franz, J. H. Pixley, N. Poccia, and P. Kim, Science 382, 1422–1427 (2023).
  * Wang  _et al._ [2023] H. Wang, Y. Zhu, Z. Bai, Z. Wang, S. Hu, H.-Y. Xie, X. Hu, J. Cui, M. Huang, J. Chen, Y. Ding, L. Zhao, X. Li, Q. Zhang, L. Gu, X. J. Zhou, J. Zhu, D. Zhang, and Q.-K. Xue, Nat Commun 14, 5201 (2023).
  * Margalit  _et al._ [2022] G. Margalit, B. Yan, M. Franz, and Y. Oreg, Phys. Rev. B 106, 205424 (2022).
  * Tummuru  _et al._ [2022] T. Tummuru, S. Plugge, and M. Franz, Phys. Rev. B 105, 064501 (2022).
  * Song  _et al._ [2022] X.-Y. Song, Y.-H. Zhang, and A. Vishwanath, Phys. Rev. B 105, L201102 (2022).
  * Vakili  _et al._ [2023] H. Vakili, M. Ali, M. Elekhtiar, and A. A. Kovalev, Phys. Rev. B 108, L140506 (2023).
  * Wu  _et al._ [2017] J. Wu, H. Yuan, M. Meng, C. Chen, Y. Sun, Z. Chen, W. Dang, C. Tan, Y. Liu, J. Yin, Y. Zhou, S. Huang, H. Q. Xu, Y. Cui, H. Y. Hwang, Z. Liu, Y. Chen, B. Yan, and H. Peng, Nat. Nanotechnol. 12, 530 (2017).
  * Asano [2001] Y. Asano, Phys. Rev. B 63, 052512 (2001).
  * Lu  _et al._ [2023] B. Lu, S. Ikegaya, P. Burset, Y. Tanaka, and N. Nagaosa, Phys. Rev. Lett. 131, 096001 (2023).
  * Ostroukh  _et al._ [2016] V. P. Ostroukh, B. Baxevanis, A. R. Akhmerov, and C. W. J. Beenakker, Phys. Rev. B 94, 094514 (2016).
  * Zuo  _et al._ [2017] K. Zuo, V. Mourik, D. B. Szombati, B. Nijholt, D. J. van Woerkom, A. Geresdi, J. Chen, V. P. Ostroukh, A. R. Akhmerov, S. R. Plissard, D. Car, E. P. A. M. Bakkers, D. I. Pikulin, L. P. Kouwenhoven, and S. M. Frolov, Phys. Rev. Lett. 119, 187704 (2017).
  * de Vries  _et al._ [2018] F. K. de Vries, T. Timmerman, V. P. Ostroukh, J. van Veen, A. J. A. Beukman, F. Qu, M. Wimmer, B.-M. Nguyen, A. A. Kiselev, W. Yi, M. Sokolich, M. J. Manfra, C. M. Marcus, and L. P. Kouwenhoven, Phys. Rev. Lett. 120, 047702 (2018).
  * Baumgartner  _et al._ [2021] C. Baumgartner, L. Fuchs, A. Costa, S. Reinhardt, S. Gronin, G. C. Gardner, T. Lindemann, M. J. Manfra, P. E. Faria Junior, D. Kochan, J. Fabian, N. Paradiso, and C. Strunk, Nat. Nanotechnol. 17, 39–44 (2021).
  * Legg  _et al._ [2023] H. F. Legg, K. Laubscher, D. Loss, and J. Klinovaja, Phys. Rev. B 108, 214520 (2023).
  * Du  _et al._ [2023] B.-H. Du, M. Yang, and L.-B. Hu, Chinese Phys. B 32, 077201 (2023).
  * Hess  _et al._ [2023] R. Hess, H. F. Legg, D. Loss, and J. Klinovaja, Phys. Rev. B 108, 174516 (2023).
  * Güngördü and Kovalev [2022] U. Güngördü and A. A. Kovalev, J. Appl. Phys. 132, 041101 (2022).


# Performance evaluation of accelerated complex multiple-precision LU
decomposition

Tomonori Kouya
Shizuoka Institute of Science and Technology
OrcidID:0000-0003-0178-5519

###### Abstract

The direct method is one of the most important algorithms for solving linear
systems of equations, with LU decomposition comprising a significant portion
of its computation time. This study explores strategies to accelerate complex
LU decomposition using multiple-precision floating-point arithmetic of the
multiple-component type. Specifically, we explore the potential efficiency
gains using a combination of SIMDization and the 3M method for complex matrix
multiplication. Our benchmark tests compare this approach with the direct
method implementation in MPLAPACK, focusing on computation time and numerical
errors.

##  1 Introduction

Linear systems of equations are frequently encountered in various scientific
computations, therefore various solvers have been developed to address the
computation environment and mathematical properties of these equations. The
direct method, including LU decomposition, is considered standard because it
can solve equations with a limited number of algebraic operations and maintain
numerical stability through the use of a pivoting strategy. The HPL benchmark
test[8] employs the direct method to select the top 500 supercomputers
concurrently.

In recent years, there has been much research on mixed-precision computing,
where the number of precision digits (the length of mantissa) of the floating-
point (FP) numbers used can be tailored to the requirements of the problem.
This approach allows for improving computation performance by combining
various kinds of hardware-oriented IEEE754 FP arithmetic, such as binary16
(11-bit precision), binary32 (24-bit), and binary64 (53-bit). Nevertheless,
inadequate accuracy in calculation results would pose a significant
disadvantage. There is an increasing demand for “reproducible computing" to
minimize environmental dependence of the calculation results. For good-
condition problems, speed-up should be pursued by using fast binary16, 32, and
64 FP arithmetic, which have shorter calculation digits. Conversely, for bad-
condition problems, multiple-precision FP operations exceeding binary64 should
be employed to ensure accuracy and reproducibility. In particular, for stable
solutions of nonlinear problems based on linear computation, it is necessary
to minimize rounding errors in the intermediate results. This, utilizing a
standard multiple-precision FP arithmetic library known for its user-friendly
interface and high performance, such as QD[1] and MPFR[9] becomes essential.
MPLAPACK/MPBLAS[7], serving as a de-facto standard multiple-precision linear
computation environment incorporating these libraries, becomes essential. We
have developed three high-performance multiple-precision real basic linear
calculations of double-double (DD, 106 bits), triple-double (TD, 159 bits),
and quadruple-double (QD, 212 bits) using AVX2[5]. The results have been
compiled into BNCmatmul[4] and released as open-source software.

Based on these results, this paper discusses the performance evaluation of
complex multiplicative precision LU decomposition, particularly focusing on
multi-component DD, TD, and QD precision, while also employing the 3M method,
a type of complex multiplication method. Benchmark tests compared the
acceleration of complex arithmetic using AVX2, the LU decomposition using the
Strassen algorithm and the Ozaki scheme, and parallelization using OpenMP.

##  2 Acceleration of complex basic linear computation

As mentioned above, we have already implemented real basic linear computation
and LU decomposition with SIMD operations such as AVX2 and a library of
multiple-precision linear computation supporting parallelization with OpenMP.
Additionally, we have developed multi-component fixed-precision operations
supporting DD, TD, and QD precision, which are faster than MPLAPACK/MPBLAS,
and arbitrary precision floating-point operations based on MPFR. Building upon
this foundation, we have implemented complex BLAS[6]. It is worth noting that
while we have published a collection of results on arbitrary-precision complex
LU decomposition[6], we have yet to parallelize it. In the following, we
explain the strategies for improving the performance of complex linear
computation, which serves as the basis for complex LU decomposition.

###  2.1 The 3M method for complex linear computation

For any $a,b\in\mathbb{C}$, standard complex multiplication $ab\in\mathbb{C}$
is calculated using the 4M (multiplication) method, which involves four real
multiplications:

$\begin{split}ab&:=(\mathrm{Re}(a)\mathrm{Re}(b)-\mathrm{Im}(a)\mathrm{Im}(b))+\left\\{\mathrm{Im}(a)\mathrm{Re}(b)\right.\\\
&\left.+\mathrm{Re}(a)\mathrm{Im}(b)\right\\}\cdot\mathrm{i},\ \end{split}$
(1)

In contrast, the 3M method is executed as follows:

$\begin{split}t_{1}&:=\mathrm{Re}(a)\mathrm{Re}(b),\
t_{2}:=\mathrm{Im}(a)\mathrm{Im}(b)\in\mathbb{R}\\\
ab&:=(t_{1}-t_{2})+\left\\{(\mathrm{Re}(a)+\mathrm{Im}(a))\right.\\\
&\left.(\mathrm{Re}(b)+\mathrm{Im}(b))-t_{1}-t_{2}\right\\}\cdot\mathrm{i}\in\mathbb{C},\end{split}$
(2)

The 3M method retains the same computation in the real part as the 4M method
but reduces the number of real multiplications in the imaginary part, hence
its name. We cannot ascertain the efficiency of the 3M method without
benchmark testing of its implementations because it includes three additional
additions compared to the 4M method. Nevertheless, The 3M method has proven to
be valuable for tasks such as matrix multiplications[10] and arbitrary-
precision complex arithmetic[3][2], particularly in scenarios with extensive
computation times for multiplication.

We implement multiple-precision matrix multiplication based on the 3M method.
For any $A\in\mathbb{C}^{m\times n}$ and $B\in\mathbb{C}^{n\times l}$:

$\begin{split}T_{1}&:=\mathrm{Re}(A)\mathrm{Re}(B),\
T_{2}:=\mathrm{Im}(A)\mathrm{Im}(B)\in\mathbb{R}^{m\times l}\\\
AB&:=(T_{1}-T_{2})+\left\\{(\mathrm{Re}(A)+\mathrm{Im}(A))\right.\\\
&\left.(\mathrm{Re}(B)+\mathrm{Im}(B))-T_{1}-T_{2}\right\\}\cdot\mathrm{i}\in\mathbb{C}^{m\times
l},\end{split}$ (3)

This approach allows for complex matrix multiplication based on existing high-
performance real matrix multiplication techniques. While the efficiency of
this implementation remains uncertain due to the additional three additions
and subtractions of matrices, benchmark tests are necessary to evaluate the
performance of our multiple-precision complex matrix multiplication.

###  2.2 Acceleration using AVX2 SIMDization

As described above, our complex basic linear computation, which includes the
3M method is constructed upon real multiple-precision basic linear subprograms
(BLAS) developed with SIMDization such as AVX2 in our library. Complex vectors
and matrices are allocated as separate arrays for real and imaginary parts in
multi-component-type multiple-precision FP data types. On the other hand, for
arbitrary-precision FP arithmetic, we adopt the same data type (a set of real
and imaginary parts) as MPC[3] for a single element. MPC, based on MPFR,
employs the 3M method for element-wise complex multiplication, ensuring high
speed simply by using it. The data structure for these complex vectors and
matrices is illustrated in Figure1.

Figure 1: Data structure for multiple-precision complex vectors
$\mathbf{v}\in\mathbb{C}^{3}$ and matrices $A\in\mathbb{C}^{3\times 3}$

This approach simplifies the implementation of fast complex matrix
multiplication based on the formula in (3) for multi-component-type multiple-
precision calculations. We have already implemented four matrix multiplication
algorithms: simple matrix multiplication using triple loops, block matrix
multiplication, Strassen matrix multiplication (and the Winograd algorithm),
and the Ozaki scheme. Complex matrix multiplication is achieved by combining
them.

Therefore, SIMDization of complex arithmetic units is not necessary for basic
matrix multiplication alone. However, in order to speed up complex LU
decomposition, we decided to SIMDize the arithmetic, similar to the real LU
decomposition, to speed up the element-by-element forward elimination
calculation. For multi-component complex vectors and matrices, since the real
and imaginary parts are separated, SIMD Load/Store instructions can be used to
retrieve and store the elements. Using this approach, complex DD operations
using AVX2 (`_bncavx2_rcdd_add`, `sub`, `mul`, `div`, `inv`(reciprocal)) can
be implemented as shown in Figure2.

Figure 2: Complex linear computation with AVX2

Our SIMDized complex arithmetic is based on a combination of SIMDized real
arithmetic.

###  2.3 Parallelization with OpenMP

For parallelization, multi-component-type multiple precision complex linear
computation is based on real linear computation units. The real matrix
multiplication based on the Strassen algorithm has been parallelized in
sections, as illustrated in Figure3. The parallelization is conducted section
by section as depicted in. Due to the nested nature of the parallelization,
performance improvement beyond 8 threads is not currently expected. Arbitrary-
precision complex matrix multiplication is implemented using MPC operations,
including parallelization.

Figure 3: The parallelized Strassen algorithm

The Ozaki scheme divides the computation into lengths to avoid introducing
rounding errors in the base short-precision matrix multiplication (xGEMM,
using binary64 DGEMM in this case). The product obtained by xGEMM is then
added with a multiple-precision matrix addition to obtain a high-precision
matrix product. To parallelize this process, it is necessary to use parallel
xGEMMs or perform OpenMP loop parallelization for each xGEMM, as illustrated
in Figure4. In this case, the latter approach was used for parallelization.

Figure 4: The parallelized Ozaki scheme

###  2.4 Benchmark test for complex matrix multiplication

We use the following EPYC computing environment to present the results
obtained through benchmark tests.

CPU and OS

AMD EPYC 9354P 3.7 GHz 32 cores, Ubuntu 20.04.6 LTS

C/C++

Intel Compiler version 2021.10.0

Compiler Option

`-O3 -std=c++11 -fp-model precise`

with AVX2

-axCORE-AVX2 -march=skylake -mtune=skylake -mcpu=skylake

MPLAPACK

2.0.1, GNU MP 6.2.1, MPFR 4.1.0, MPC 1.2.1

We address the problem of computing the matrix product $C:=AB$ for given
complex square matrices $A$ and $B$$\in\mathbb{C}^{n\times n}$. The elements
of $A$ and $B$ used for benchmarking purposes are uniformly distributed random
numbers in the interval $[0,1]$. For the Ozaki scheme, we used a maximum of 6
divisions for DD accuracy, 8 divisions for TD accuracy, and 12 divisions for
QD accuracy to ensure sufficient accuracy.

At first, we show the results of DD-precision matrix multiplication on
Figure5. The left figure shows the computation time in seconds, and the right
figure shows the corresponding speed-up ratio achieved by parallelization with
OpenMP.

The left figure shows that our implementations of the Strassen and Ozaki
schemes achieve up to 100 times performance improvement compared to the Cgemm
routine of MPBLAS for serial calculations. Additionally, the 3M method
exhibits a speed improvement of about 1.3 times over the 4M method, while the
Strassen method achieves a speed-up of 3 to 4 times. The computation time of
the Ozaki scheme is nearly equivalent to that of the AVX2 Strassen scheme due
to the small number of partitions. Furthermore, the block matrix
multiplication based on the 3M method is slightly slower than the Strassen
matrix multiplication.

In OpenMP parallelization, the Strassen matrix multiplication does not exhibit
any acceleration, even when using 32 threads. Similarly, the Ozaki scheme
demonstrates only a modest speed-up of about twice as fast. In contrast, the
block matrix multiplication achieves a significant improvement, with speeds up
to 12 times faster.

Figure 5: DD-prec.: Computation time in seconds (left) and speed-up ratio when
using 32 threads (right)

Next, we show the results obtained by TD-precision matrix multiplication in
Figure6.

The AVX2 SIMDized TD-precision Strassen, block matrix multiplication is
observed to be slower than the non-SIMDized version, although the reason for
this discrepancy is unclear. The 3M method demonstrates a speed improvement of
about 1.3 times faster compared to the 4M method, which is common with DD
precision. The Ozaki scheme achieves speed up to about 8 times faster than the
Strassen 3M method.

In OpenMP parallelization, the Strassen matrix multiplication is observed to
be slower, and the Ozaki scheme is only about 2 times faster, even when
utilizing 32 threads, as in the DD accuracy. The block matrix multiplication
is up to 20 times faster.

Figure 6: TD-prec. : Computation time in second (left) and speed-up ratio when
using 32 threads (right)

Finally, the results of QD-precision matrix multiplication are shown in
Figure7.

Similar to DD and TD precision, the Strassen matrix multiplication achieves
approximately 2 times faster performance with AVX2 optimization and 1.3 times
faster performance with the use of the 3M method.

In OpenMP parallelization, the Strassen matrix multiplication demonstrates
faster performance, while the Ozaki scheme achieves only 2 times speed-up
compared to the DD- and TD-precision schemes, despite utilizing 32 threads.
The block matrix multiplication is 20 times faster than TD precision.

Figure 7: QD-prec. : Computation time in seconds (left) and speed-up ratio
when using 32 threads (right)

These results demonstrate that while the 3M method can speed up serial complex
matrix multiplication, it does not speed up the parallelization in relation to
the number of threads, except in the case of block matrix multiplication.

##  3 Acceleration of multiple-precision complex LU decomposition

The complex LU decomposition is implemented based on a complex basic linear
computation that integrates all the implemented speed-up techniques (3M
method, SIMD, OpenMP parallelization) to achieve high performance. The
following is a concise description of the implementation.

The target $n$-dimensional complex linear system is shown as

$A\mathbf{x}=\mathbf{b},$ (4)

where $A\in\mathbb{C}^{n\times n}$, $\mathbf{x}\in\mathbb{C}^{n}$ and
$\mathbf{b}\in\mathbb{C}^{n}$. We solve (4) to determine the unknown vector
$\mathbf{x}$ given $A$ and $b$.

The current standard implementation of LU decomposition uses fast matrix
multiplication. As shown in Figure8, a certain block size $K$($=\mbox{NDIM}$)
is set in advance, and $A-L_{12}U_{21}$ is performed to compute the $A_{22}$
part. Consequently, the computation performance varies depending on the value
of $K$. The LU decomposition (Cgetrf function) of MPLAPACK used for comparison
also uses complex matrix multiplication (Cgemm function). In all cases, LU
decomposition is computed using a partial pivoting strategy.

Figure 8: LU decomposition with matrix multiplication

In addition, a standard LU decomposition approach (referred to as “Normal LU"
or “Normal"), which eliminates rows sequentially without updating or
performing matrix multiplication, was implemented using complex arithmetic
with AVX2, as illustrated in Figure2, owing to the ease of speeding up the
process with SIMD instructions. AVX2 demonstrated effectiveness in real LU
decomposition, with the Ozaki scheme emerging as the fastest in this case. The
question remains: will the same effect be obtained for complex LU
decomposition ?

##  4 Benchmark test of complex LU decomposition

During the benchmark test, the computation time and maximum error of the
numerical solution of the LU decomposition were measured by varying $K$ while
calculating DD, TD, and QD precision, respectively. The graphs below show the
computation time for normal LU and LU decomposition using matrix
multiplication, as well as the relative error of the numerical solutions
obtained through forward and backward substitutions.

Our focus is on solving the $n$-dimensional complex linear system depicted in
equation (4). For benchmarking purposes, the elements of $A$ are generated as
uniform random numbers within the range $[0,1]$.

The true solution $\mathbf{x}$ is $(\mathbf{x})_{k}=k+k\cdot\mathrm{i}$, and
the constant vector is derived by computing $\mathbf{b}:=A\mathbf{x}$ using
2048-bit precision calculation.

###  4.1 DD-precision computation

First, the results of the DD accuracy calculation are presented. This involved
measuring the relative error and computation time while increasing $1\leq
K\leq n$ by 32.

Examining the relative error (Figure9), the MPLAPACK and AVX2 SIMDized normal
LU decompositions demonstrate superior performance. The normal LU
decomposition and LU decomposition using the Ozaki scheme and Strassen matrix
multiplication yield similar results. Overall, the results also show that the
accuracy tends to decrease when the Strassen is used. For this problem, the
Ozaki scheme with DGEMM does not achieve optimal accuracy with 5 divisions,
necessitating more than 6 divisions. Therefore, results are shown for 5, 6,
and 7 divisions.

Figure 9: DD-precision complex LU decomposition: Max. relative error

Figure10 demonstrates the measured serial computation time (left) and
32-thread parallel computation time (right). The results with normal LU
decomposition and AVX2 acceleration are represented by solid orange and dashed
lines, respectively. The relative error of MPLAPACK is depicted as a solid
black line, although it is not included in the computation time graph due to
the significant difference in scale (130.7 s).

Figure 10: DD-precision complex LU decomposition: Seconds of serial
computation (left), Seconds of 32-thread computation (right)

The AVX2 normal LU decomposition exhibits the fastest serial computation,
while the AVX2 Strassen LU decomposition and Ozaki scheme show similar
computation times, with smaller $K$ values requiring less time. However, they
are not as fast as the AVX2 normal LU decomposition.

In the case of DD-precision calculations, OpenMP parallelization of matrix
multiplication initially increases the computation time with 2 threads, after
which it becomes faster proportionally to the number of threads. Therefore,
with 32 threads, the computation time is approximately the same as for serial
computation, and as a result, the parallelization effect is not enhanced. On
the other hand, as depicted in Figure11, the normal LU decomposition achieves
speed-up in proportion to the number of threads, and the AVX2 version remains
the fastest even with 32 threads.

Figure 11: Parallelization of normal complex LU decomposition: DD-precision

The reason for this phenomenon is believed to be the higher computation
complexity inherent in complex LU decomposition compared to real LU
decomposition, despite the enhanced performance of complex operations achieved
through SIMD. As previously mentioned, while SIMD algorithms are not essential
for complex basic linear computation implementations, they may facilitate
acceleration in algorithms such as the Fast Fourier Transform, which are
required locally as LU decomposition.

###  4.2 QD-precision computing

Next, we present the results of the QD accuracy calculations.

Figure12 displays the relative errors obtained from QD-precision benchmark
tests. As depicted in this figure, the AVX2 SIMDized normal LU decomposition
outperforms MPLAPACK, while the LU decomposition using the Ozaki scheme and
the Strassen matrix multiplication exhibits larger relative errors compared to
the normal LU decomposition.

Figure12

Figure 12: QD-prec. complex LU decomposition: Max. relative error

Figure13 illustrates the computation time obtained from QD-precision benchmark
tests. The graph showcases the measured serial computation time (left) and
32-thread parallel computation time (right) with the standard LU decomposition
and the AVX2 accelerated approach, represented by solid and dashed orange
lines, respectively. For this particular problem, the Ozaki scheme using DGEMM
does not reach the highest QD accuracy with 11 divisions and requires more
than 12 divisions, hence results for 11, 12, and 13 divisions are displayed.
The relative error for MPLAPACK is represented by the solid black line, but
the serial computation time is 595.4 seconds, so it is not included in the
graph of computation time because it is on a different scale.

Figure 13: QD-prec. complex LU decomposition:Seconds of serial computation
(left), Seconds of 32-thread computation (right)

For serial calculations, the AVX2 SIMDized normal LU decomposition generally
exhibits the fastest performance. However, the Ozaki scheme LU decomposition
is slightly faster when $K$ is small. The Strassen LU decomposition is not as
fast as the AVX2 SIMDized normal LU decomposition, even when AVX2 SIMDized.

In QD-precision calculations, the OpenMP parallelization of matrix
multiplication accelerates the computation linearly with the number of
threads. However, the parallelization effect of the normal LU decomposition
surpasses this, with the AVX2 SIMDized normal LU decomposition remaining the
fastest overall, except in certain cases where the AVX2 SIMDized Strassen LU
decomposition outperforms it. As illustrated in Figure14, although the
parallelization efficiency is high, it is slowed down by AVX2.

Figure 14: Parallelization of normal complex LU decomposition: QD-prec.

###  4.3 TD-precision computing

Finally, we show the results of the TD-precision calculation, which represents
an intermediate accuracy level between DD and QD.

In this case, as depicted in Figure15, for TD accuracy, the Ozaki scheme with
DGEMM does not achieve the highest accuracy with 7 divisions, requiring 8 or
more.

Figure 15: TD-precision complex LU decomposition: Max. relative error

The relative error analysis reveals that, similar to the case of QD accuracy,
the AVX2-ized normal LU decomposition performs the best. Furthermore, the
relative error observed in the LU decomposition using the Ozaki scheme and
Strassen matrix multiplication is greater compared to the normal LU
decomposition.

A graph depicting computation time is presented in Figure16, omitting the
32-thread parallelization using TD accuracy due to its significantly increase
in computation time for small $K$.

Figure 16: TD-precision complex LU decomposition: Seconds of serial
computation (left), Seconds of 32-thread computation (right)

For serial calculations, the AVX2 SIMDized normal LU decomposition generally
does not outperform others, with the LU decomposition using the Ozaki scheme
being largely faster for small $K$ values. In TD-precision calculations,
OpenMP parallelization of the matrix multiplication displays straightforward
speed-up proportionate to the number of threads. However, the parallelization
effect of the normal LU decomposition outweighs it. The normal LU
decomposition with AVX2 remains the fastest, except in cases where LU
decomposition with the Ozaki scheme shows better performance. As shown in
Figure17, the parallelization efficiency is considerable; however, the ratio
is slowed down compared with DD and QD normal LU decomposition.

Figure 17: Parallelization of normal complex LU decomposition: TD-precision

###  4.4 Summary of high-performance complex LU decomposition

These benchmark results indicate that the AVX2 SIMDized normal LU
decomposition generally exhibits the fastest performance in terms of both DD
and QD accuracy, with good parallelization. A summary solely focusing on the
computation time is presented in Table 1, highlighting the minimum computation
time achieved for Strassen with AVX2 for DD and QD, and without AVX2 in TD.

Table 1: Minimum computation time of complex LU decomposition in second DD
precision, $n=1024$
---
| MPLAPACK | Normal | Normal AVX2 | Strassen AVX2 | OZ
Serial | 130.7 | 7.9 | 1.5 | 2.1 ($K=32$) | 2.0 ($K=64$)
Parallel(32T) | $-$ | 0.29 | 0.18 | 2.1 ($K=32$) | 1.9 ($K=96$)
TD precision, $n=1024$
Serial | $-$ | 54.0 | 41.8 | 46.7 ($K=384$) | 11.9($K=64$)
Parallel(32T) | $-$ | 7.4 | 4.9 | 8.5 ($K=992$) | 4.3($K=192$)
QD precision, $n=1024$
Serial | 595.4 | 227.6 | 64.7 | 60.1($K=32$) | 40.7($K=64$)
Parallel(32T) | $-$ | 8.3 | 6.1 | 3.6($K=32$) | 9.8($K=352$)

Finally, it was demonstrated that the 32-thread AVX2 normal LU decomposition
outperformed MPLAPACK with serial computation by approximately 726 times for
DD accuracy, and approximately 91 times for QD accuracy. These findings
suggest that the default setting for complex LU decomposition should be the
AVX2 SIMDized normal LU decomposition.

##  5 Conclusion and future works

The implementation of the multiple precision complex basic linear computation,
based on the enhanced multiple-precision real basic linear computation,
demonstrated improved speed by using the 3M method for matrix multiplication.
Additionally, the complex LU decomposition was found to be sufficiently fast
through the utilization of SIMD complex operations and parallelization with
OpenMP. However, it was observed that the speed-up of the LU decomposition is
significant enough to surpass conventional LU decomposition algorithms in DD-
and QD-precision, given the current complex matrix multiplication, both in
serial and parallel computations.

Future work will focus on enhancing the performance of numerical algorithms
for ill-conditioned nonlinear problems by improving the following functions.

  1. 1.

Implementation of a multiple precision sparse matrix and extension of our
implementation to use real and complex sparse matrices.

  2. 2.

Speed up reproducible computations for binary64 and other types of precision
arithmetic.

  3. 3.

Extending fast, multi-precision computation to the Python environment.

This includes seeking further speed-up of real and complex basic linear
computation of dense matrices and LU decomposition. Of particular interest is
confirming the potential speed-up of SIMD-enhanced multi-component-type real
and complex arithmetic for other numerical algorithms.

####  5.0.1 Acknowledgment

This research is supported by the JSPS Grant-in-Aid for Scientific Research
23K11127.

## References

  * [1] Bailey, D.: QD. https://www.davidhbailey.com/dhbsoftware/
  * [2] Bini, D.A., Fiorentino, G.: Design, analysis, and implementation of a multiprecision polynomial rootfinder. Numerical Algorithms 23, 127–173 (2000). https://doi.org/https://doi.org/10.1023/A:1019199917103
  * [3] Enge, A., Théveny, P., Zimmermann, P.: MPC. http://www.multiprecision.org/mpc/
  * [4] Kouya, T.: BNCmatmul. https://github.com/tkouya/bncmatmul
  * [5] Kouya, T.: Acceleration of multiple precision matrix multiplication based on multi-component floating-point arithmetic using avx2. In: Gervasi, O., Murgante, B., Misra, S., Garau, C., Blečić, I., Taniar, D., Apduhan, B.O., Rocha, A.M.A.C., Tarantino, E., Torre, C.M. (eds.) Computational Science and Its Applications – ICCSA 2021. pp. 202–217. Springer International Publishing, Cham (2021)
  * [6] Kouya, T.: Acceleration of complex matrix multiplication using arbitrary precision floating-point arithmetic (2023). https://doi.org/10.48550/arXiv.2307.06072
  * [7] MPLAPACK/MPBLAS: Multiple precision arithmetic LAPACK and BLAS. https://github.com/nakatamaho/mplapack
  * [8] Petitet, A., Whaley, R.C., Dongarra, J., Cleary, A.: Hpl - a portable implementation of the high-performance linpack benchmark for distributed-memory computers
  * [9] Project, M.: The MPFR library. https://www.mpfr.org/
  * [10] Van Zee, F.G., Smith, T.M.: Implementing high-performance complex matrix multiplication via the 3m and 4m methods. ACM Trans. Math. Softw. 44(1) (jul 2017). https://doi.org/10.1145/3086466.