# Nonlinearity in Oscillatory Flows over Sand Ripple

Ruma Dutta111Electronic address<EMAIL_ADDRESS>Dept of Civil
engineering & Geodetic Science, The Ohio State Univ  Ethan
Kubatko222Electronic Address<EMAIL_ADDRESS>Dept. of Civil Engineering &
Geodetic Science, The Ohio State University

###### Abstract

In this report, We investigated the nonlinear phenomena in the study of flow
dynamics of velocity component. In our studies, we observed nonlinear term in
the vertical component of velocity by vittori et al vit The time series
simulation of vertical component which increases with Reynold stress value. We
developed direct numerical simulation under two dimensional grid system to
study the flow dynamics and vorticity parameter. Flow pattern and flow
dynamics near wavy boundary wall in the victinity of ripple bottom was
readdressed under direct numerical simualtion(DNS) framework. Both horizontal
and vertical component of fluid velocity were studied under pulating force of
flow. Vorticity is calculated under complex framework by taking into higher
order interaction term. we tried to carry out similar simulation with same
particle ejection in the viscous bed using DNS simulation for pulsating flow.
Our focus was to observe particle motion using DNS simulation and study the
particle phase under vortex structures formed here.

##  I Introduction

Many observations were made cencerning a complex bed form pattern in the
victinity of offshore region. Temporal chaos in fluid turbulence is the
symptotatic of spatial chaos. The chaos in turbulence can be studied can be
studied by many numerical approaches.

As widely discussed, the chaotic pattern is detected in turbulence both for
fixed geometric configuration increasing the flow Reynolds number and for
fixed characteristic of oscillatory flow increasing the amplitude of the wall
waviness. From an analysis of the The tide driven current in the offshore
region is the m major source sediment transport in beach areas. The current
deflects toward the crests due to increase in bottom friction with a
decreasing water depth. Hence the cross ridge velocity increases to satisfy
continuity whereas the along ridge velocity decreases owing to increase of
bottom friction. Huthnance was the first to present a mathematical description
using simple flow model. He used a simplified version of depth average shallow
water equation with a power law relationship for sediment transport corrected
for downhill gravitational transport.

and turbulence pattern under complex bed-fluid intearction mechanism. In fact,
vast sea water(both offshore and near shore) is not clear Underneath sea with
complex bed form having ripple is observed occasionally. sleath
Indeed many observations were made concerning a complex bed form pattern in
the victinity of offshore region. The bottom topography is indeed very complex
in nature. These ripples are formed due to oscillatory nature of turbulence
flow over sand bed form in a complex manner where the bottom topography often
takes form of The chaotic phenomena related to turbulence studies were
observed by blond brick or tile pattern depending on the complexity of the
bed-fluid interaction under turbulence flow. Flow becomes more complex near
the shore region where the flow due to oscillatory nature affects the ripple
formation. It is believed that waves are mainly responsible for sediment
transport where currents carry the entrained sediments away. Field
observations focussing tide driven region indicate the presence of symmetrical
and asymmetrical waves with crest almost perpendicular to the direction of
main current and characterized by wavelengths of few hundered metres. In this
region, full understanding of turbulence with oscillatory flow is needed to
understand the sediment flow. Turbulent fluctuations due to oscillatory nature
of flow are usually confined within a thin oscillating boundary layer. This
situation makes very difficult to take experimental data accurately in this
region which makes sometimes the interpretation of data controversial. There
are different numerical approaches to study these problems such as Reynolds
Averaged Naveir Stokes Equation(RANS), ‘Direct Numerical Simulation(DNS). LES
model is based on solving Navier Stokes equation with large eddy function and
small eddies are neglected in this calculation. It well known for its
simplicity and low computer cost but can not give more insight in the complex
phenomena. On the other hand, Reynolds Averaged Navier Stokes Equation alsoo
known as RANS model is based on On the other hand, DNS simulation does not
involve any turbulence model but uses unsteady flow using grid system that are
sufficiently fine to resolve all scales of motion. A first attempt to explain
the mechanism of this flow made by Hara an Mei hara who developed a three
dimensional model investigating the stability of Stokes layer induced by sea
wave. The associated steady state along the ripple surface showed a tendency
to accumulate sediment particle in various pattern. B Blondeaux & Vittori blon
tried to model sand ripple based on brick pattern form on the basis of three
dimensional simulation model. They studied three dimensional vortex structure
in the oscillatory form of flow under two dimensional ripple bed. In both
cases, a unidirectional oscillatory flow is considered and fluid particles are
considered on top of the bottom boundary to oscillate to and fro. Ripple at
sea bed affects the sediment transport rate and cause additional energy
dissipation enhancing mixing in the vicinity of ripple. However the detail
knowledge of the flow structure and the dynamics of vortex structure generated
by flow seperation is not clear in this region. Recently Scandura & Blondeaux
Scan studied by means of numerical simulations, the flow induced by wavy wall
under uniform oscillatory motion. They observed in the simulation that
velocity is periodic under weak flow and vorticity is shed just above the
crest which has a tendency for pitchfork bifurcation above critical value of
velocity.
In our numerical simulation, we tried to readdress the problem of chaos in
velocity and vorticity field and observed similar bifurcation at critical
value of velocity. Our direct numerical simulation is based on finite
difference scheme for oscillatory flow of fluid and we observe development of
bifurcation in the normal and streamline flow component which increases with
$u_{\rm 0}$. The nonlinearity nature of the vertical component of the velocity
also shows similar pattern.

##  II Formulation of the Problem

The problem was formulated in the following way, We consider incompressible
fluid of density $\rho$ and kinematic viscosity $\nu$ induced close to a wavy
wall by a uniform oscillating pressure gradient. We define Cartesian
orthogonal corordinate We start with Navier Stokes equation for an
incompressible fluid flow in rectangular $(x,y,z)$ coordinates. We also
consider wall profile described parametrically by the relationship

$y=-{\frac{h}{2}}[cos({\it k}{\xi})+{\sum_{n=1}}^{N}c_{\rm n}cos({\gamma}_{\rm
n})]$ (1) $x={\xi}+\frac{h}{2}[sin(k{\xi})+{\sum_{n=1}}^{N}c_{\rm
n}sin({\gamma_{n}})]$ (2)

where ${\rm k}=\frac{2{\pi}}{l}$ is the wavenumber of the waviness. $\xi$ is a
dummy variable and $\gamma_{\rm n}=nk{\xi}+{\phi}_{\rm n}$.

$\displaystyle\ \frac{\partial u}{\partial t}+{\bf u}\cdot\bigtriangledown
u=-\frac{1}{\rho}\frac{\partial p}{\partial
x}+{\nu}\bigtriangledown^{2}u+F^{(x)}$ (3) $\displaystyle\ \frac{\partial
w}{\partial t}+{\bf u}\cdot\bigtriangledown w=-\frac{1}{\rho}\frac{\partial
p}{\partial z}+{\nu}\bigtriangledown^{2}w$ (4) $\displaystyle\
\bigtriangledown\cdot{\bf u}=0$ (5)

For Direct numerical simulation algorithm, we choose collocated, nonstaggered
grid system. The algorithm we

$\displaystyle-\frac{1}{\rho}\frac{\partial p}{\partial
y}+{\nu}\bigtriangledown^{2}v$ (6) $\displaystyle\ \frac{\partial w}{\partial
t}+{\bf u}\cdot\bigtriangledown w=-\frac{1}{\rho}\frac{\partial p}{\partial
z}+{\nu}\bigtriangledown^{2}v$ (7) $\displaystyle\ \bigtriangledown\cdot{\bf
u}=0$ (8)

where the field velocities ( u,v,w) are along (x,y,z) directions repectively.

For the sediment particle, the basic equation is controlled by spherical
particle moving under gravity in viscous fluid.

###  II.1 Discussion and Conclusion of the Results

We consider the flow of an incompressible viscous fluid of density $\rho$ and
kinematic viscosity $\nu$ induced close to wavy wall by a uniform oscillating
pressure gradient. The nonlinear fluctuation and periodicity was observed here
for velocity component both in streamline and vertical flow field. The
nonlinearity increases with $u_{\rm o}$ above the threshold value. Whwn the
shear stress experienced by the interface betwen the flowing fluid and the
resting particles is low, the flow is unable to entrain the particles lying on
the bed, which then remains immobile. As the shear stress increases,

##  III Acknowledgement

This work was supported by Naval Research Lab Grant.

## References

  * (1) O.E Landford III, Annual Rev. Fluid Mech., 14, 347 (1982).
  * (2) J. Guckenheimer, Annual Rev. FLuid Mech., 18, 15 (1986).
  * (3) Hara. T & Mei C.C (1990), Centrifugal Instability of an oscillatory flow over periodic ripples Journal of Fluid Mechanics, 217, 1-32.
  * (4) P. Blondeaux & G. Vittori, A route to chaos in an oscillatory flow: Feigenbaum scenario, Phys. Fluids A 3(11), 2492-2495, 1991.
  * (5) Blondeaux P. 1990, Sand Ripples under seawaves Part 1, Ripple Formation, J. Fluid Mechanics,218, 1-17.
  * (6) Three dimensional oscillatory flow over steep ripples; J. Fluid Mechanics, 412, 355-478.
  * (7) P.Scandura, G.Vittori and P. Blondeaux, Bifurcations in the Oscillatory Flow over a Wavy Wall, Mechanics, 37, 305-311, 2002.


# A 120dB Programmable-Range On-Chip Pulse Generator for Characterizing
Ferroelectric Devices ††thanks: This work received funding from the European
Union’s H2020 research and innovation programme under the H2020
BeFerrosynaptic (871737), MeM-Scales (871371) projects.

Shyam Narayanan1, Erika Covi2, Viktor Havel2, Charlotte Frenkel 1, Suzanne
Lancaster 2,
Quang Duong 2, Stefan Slesazeck 2 Thomas Mikolajick 2, Melika Payvand 1,
Giacomo Indiveri 1
1Institute of Neuroinformatics, University of Zurich and ETH Zurich,
Switzerland
2NaMLab gGmbH, Dresden, Germany
Email: {shyam, melika<EMAIL_ADDRESS>
###### Abstract

Novel non-volatile memory devices based on ferroelectric thin films represent
a promising emerging technology that is ideally suited for neuromorphic
applications. The physical switching mechanism in such films is the nucleation
and growth of ferroelectric domains. Since this has a strong dependence on
both pulse width and voltage amplitude, it is important to use precise pulsing
schemes for a thorough characterization of their behavior. In this work, we
present an on-chip 120 dB programmable range pulse generator, that can
generate pulse widths ranging from 10 ns to 10 ms $\mathbf{\pm}$2.5% which
eliminates the RLC bottleneck in the device characterisation setup. We
describe the pulse generator design and show how the pulse width can be tuned
with high accuracy, using Digital to Analog converters. Finally, we present
experimental results measured from the circuit, fabricated using a standard
180 nm CMOS technology.

###### Index Terms:

Pulse Generator, FTJ, FeFET, DAC, memory characterization

##  I Introduction

Thanks to their ultra-low power operation, Complementary Metal-Oxide-
Semiconductor (CMOS)-compatibility, and non-volatility, ferroelectric hafnium
oxide-based devices such as ferroelectric capacitors (FeCAP) [1],
ferroelectric field effect transistors (FeFET) [2], and ferroelectric
tunneling junctions (FTJ) [3] are widely explored for digital storage
conventional memory array architectures [1, 2] as well as in novel bio-
inspired architectures [4]. Besides multi-bit storage in digital memories, the
analog switching properties of these nano-scale devices are important for the
realization of synaptic weight elements [5, 6] or artificial neurons [7] for
adoption in neuromorphic applications [8]. In order to allow proper device
optimization and defining the constraints for circuit design, a profound
understanding of the switching kinetics in such nano-scale devices is needed
[9].

(a)

(b)

Figure 1: Switching transitions: (1(a)) switched polarization as a function of
pulse width for different pulse amplitudes, for a bilayer FTJ device. (1(b))
Time needed to switch 50% of polarization for each voltage amplitude, expected
switching curve for HZO only (black dotted line from [10]) and asymptotic RC-
cutoff (red dotted line)

Figure 1(a) depicts the results of such a switching kinetics measurement for a
bilayer FTJ device featuring a 10 nm ferroelectric hafnium-zirconium-oxide
(HZO) and a 2 nm alumnium oxide (Al2O3) layer sandwiched between two titanium
nitride (TiN) electrodes [3]. The diagram depicts the portion of switched
polarization upon application of different switching pulses. Extracting the
respective switching voltage and time at $50\%$ switched polarization results
in the switching kinetics curve depicted in Fig. 1(b). Extrapolation of this
data yields both the data retention time of the respective device [11] as well
as the required pulse voltages for ns-switching times. The switching kinetics
of a given device will be material- and size-dependent [9], and as such, the
characterization should be performed on individual devices. As an example,
when looking at the voltage-time dependence in Fig. 1(b) and comparing to the
data for an HZO-only sample (extracted from [10]), there is a shift to larger
voltages due to the voltage drop over the Al2O3 layer, and an asymptotic
saturation of switching time between 100 ns and 1 $\mathrm{\mu}$s. This
saturation is related to a RC-delay within the FTJ samples and measurement
set-up that cannot be easily avoided. Similar issues occur for FeFET and FeCAP
single devices, even though switching times down to just 14 ns to 20 ns have
been demonstrated in integrated memory arrays [1, 2]. Obviously, the RC-delay
has a strong impact on the extrapolation of switching kinetics, especially at
shorter time scales, and consequently a pulse generator is required which is
able to perform such a thorough characterization by applying pulses with
varying amplitudes and timewidths directly on-chip. Moreover, for the
characterization of multi-bit and analogue storage a very precise control of
the switching pulses is mandatory [12], especially for the steep switching
gradients at shorter pulse times (see Fig. 1(a)). In order to tackle this
issue, we developed an optimized on-chip pulse-generator circuit, enabling the
application of well-controlled voltage pulses with timings that span over six
orders of magnitude.

##  II System Architecture

Figure 2 shows the architecture of the wide-range programmable on-chip pulse
generator. On the arrival of the input pulses (shown as $CLK$), the
programmable capacitor $C_{DAC}$ is abruptly charged to $V_{DD}$, and once the
pulse is removed, it slowly discharges through the $I_{DAC}$ current. This
results in a ramp voltage on node $V_{RAMP}$ which when compared to a
reference voltage ($V_{REF}$), generates a pulse at the output ($V_{PULSE}$).
The slope of the ramp voltage is determined by the magnitude of the current
$I_{DAC}$ and the value of $C_{DAC}$. Therefore, the width of the pulse at the
output of the comparator can be determined by the following equation:

$I_{DAC}=C_{DAC}\cdot\frac{\Delta V_{REF}}{\Delta t}$ (1)

To have a precise control on the output pulse width ($\Delta t$), the circuit
should be independent of the input pulse width, and employ a high-precision
capacitor $C_{DAC}$, high-precision discharge current $I_{DAC}$, and a low-
offset comparator. In this section we describe each of these blocks.

Figure 2: Architecture of the Pulse Generator

###  II-A Edge detector

To decouple the $V_{PULSE}$ pulse-width from that of the $CLK$ pulses, we
implemented an Edge Detector block which detects the rising edge of the $CLK$
signal and sends a fixed-width pulse to the transistor $M_{PD}$, independent
of the $CLK$ pulse width [13].

###  II-B Wide-Range Current Digital to Analog Converter (IDAC)

The wide-range IDAC block of Fig. 2 represents a 7-bit current output binary
weighted Digital to Analog Converter (DAC). It consists of identical modular
unit cells, combined in a binary weighted fashion to work as p-type current
sources. The unit cells are sized appropriately to provide a precise value of
current. Both p-type transistors where one being the current source and the
other being the gating switch were sized to be of the same width but different
lengths. The primary considerations for the sizing were the output resistance,
matching, dynamic range and leakage depending on the function of each device.
The DAC is biased by an external precise reference current which varies in
three steps to cover the large operating range of the circuit as described in
table I.

One important consideration for the design of the current DAC is to ensure the
current sources are accurate at large output currents, and that the current
source transistors are not out of the saturation. The output impedance of the
$I_{DAC}$ should be high to produce reliable output current with good
accuracy. We decided to limit the systematic error $<1\%$ to help the system
deliver a highly accurate output current. This is especially important,
because as as shown in Figure 2, the $I_{DAC}$ output is injected to the
$V_{RAMP}$ node whose voltage is changing over time, and the $I_{DAC}$ current
should not be susceptible to this change through the drain-source voltage
across the current sources of the $I_{DAC}$. We verified that the systematic
error in the $I_{DAC}$ current is under 0.3% of the desired value, if the
operating headroom is at least 400-450 mV.

TABLE I: $I_{DAC}$ Programmable Range Pulse Range | $\mathbf{I_{BIAS}}$ | $\mathbf{I_{DAC}}$
---|---|---
7.8 ns–1 $\mu$s | 160 $\mu$A | 10 $\mu$A–1.28 mA
780 ns–100 $\mu$s | 1.60 $\mu$A | 100 nA–12.8 $\mu$A
78 $\mu$s–10 ms | 16 nA | 1 nA–128 nA

We simulated the functionality of the circuit with a wide range of currents,
and verified with experimental measurements the correct functionality of the
circuit down to the lowest current setting of 1 nA (see Section IV).
Typically, the circuit’s systematic error is higher at lower currents, when
the p-type transistors operate in weak inversion. However, narrow pulses are
needed to characterize the memristive devices. Therefore, in this application
scenario to generate narrow pulses, the required $I_{DAC}$ current is in the
order of magnitude where the IDAC operates in strong inversion saturation.

We performed a MonteCarlo simulation to quantify the expected accuracy and
quality of this design, as we could not measure through a large set of silicon
samples. In this analysis, we quantified the error for five different current
settings ranging from the minimum to maximum of the 7-bit IDAC code across
different range of operation as described in Table I. The deviation from the
target value of the programmed current was calculated as:

$I_{DAC,Error}[\%]=\Big{(}\frac{I_{DAC,Sim}}{I_{DAC^{*}}}-1\Big{)}\cdot 100,$
(2)

where $I_{DAC,Sim}$ is the actual simulated $I_{DAC}$ output current and
$I_{DAC^{*}}$ is the desired output current programmed with the 7-bit code.
The results are shown in Fig. 3. For this simulation, we limited the lower
current setting of $I_{DAC^{*}}$ to 100 nA as the MonteCarlo simulation in
this technology yields a very large systematic offset (probably arising from
inaccurate sub-threshold mismatch models). It can be observed from Fig. 3 that
the lower current setting of 100 nA shows the largest error, and the highest
current setting of 1.28 mA (required for the narrowest pulse) shows the least
error, indicative of its high accuracy. Overall, the combined $3\sigma$ error
of $I_{DAC}$ is within $\pm$ 6% for all the current ranges tested, without any
calibration process. Extra tuning for higher accuracy is possible by adjusting
the capacitance value (trimming).

Figure 3: MonteCarlo Analysis of $I_{DAC}$ Error

###  II-C $C_{DAC}$

The primary capacitor in this system is an 8-bit binary weighted capacitor DAC
($C_{DAC}$). The $C_{DAC}$ is constructed out of identical unit DAC cells,
carefully laid out to provide a good level of matching among all the elements
in the array. The resolution of the $C_{DAC}$ is 62.5 fF and the total
capacitance we chose for our application is 10 pF. This enables us to adjust
the DAC value and correct/trim out any non-idealities in the design with a
precision of $\pm 0.625\%$. As the desired value of 10 pF (maximum $C_{DAC}$
value = 16pF) is not exactly at the centre of the DAC code range, we have non-
uniform correction range centered around the desired value. However, one could
potentially choose the exact center value of this DAC to be the capacitance
required for pulse generation by either adjusting the $I_{DAC}$ current or the
reference voltage $V_{REF}$ as described in Eq. 1.

###  II-D Hysteretic Comparator

For this application, the primary considerations for a comparator was to have
hysteresis to avoid false triggering due to fluctuations on the node the
comparator monitors. In addition to this the other main considerations were to
have low-offset and high-speed. For the comparator design, we chose a well-
known architecture of a hysteretic comparator [14, 15] as shown in Fig. 4.
Positive feedback is introduced in this circuit to have hysteresis.

Figure 4: Differential to single- ended hysteretic comparator

##  III Implementation and system level simulation

The pulse generator was designed and fabricated in a standard 180 nm
technology node. The final design excluding the probe pads occupy an area of
369 $\mu$m$\times$139 $\mu$m as shown in Fig. 5. The chip microphotograph is
shown in Fig. 6.

Figure 5: Pulse generator layout Figure 6: Left: pulse generator die
photograph. Right: Zoom of the contacts for ferroelectric device BEOL
integration.

The pulse generator can be programmed using internal Serial Peripheral
Interface (SPI) registers. Both $I_{DAC}$ and $C_{DAC}$ are internally wired
to their respective registers which are tightly integrated with rest of the
analog circuitry as shown in Fig. 5. The I/Os are wired up to metal probe
pads. They are powered and clocked externally through the wafer probing
equipment. The FTJs will be fabricated on the same die with Back End Of Line
(BEOL) post-processing integration steps. Therefore, the output of the pulse
generator is directly connected to the FTJs with on chip low parasitic
interconnects and switches.

Table II compares the performance of this work against the implementations
designed for similar application. This design can reliably cover a
programmable range at least 4 orders of magnitude more than [16], the digital
implementation on the same technology node. Although compared to [17], the
area benefit of our work is explainable by the smaller technology node and the
programming range is 1 order of magnitude less, the accuracy of the pulses in
our design across the full programmable range is $\mathbf{\pm}$2.5% which is
not reported in [16] and [17].

TABLE II: Performance comparison

Parameters | This work | [17] | [16]
---|---|---|---
Technology [nm] | 180 | 350 | 180
Implementation | Analog | Mixed | Digital
Silicon Area [mm2] | 0.051 | 0.203 | 0.015
Programmable Range [dB] | 120 | 150 | 42
Output Pulse Accuracy [%] |  $\mathbf{\pm}$2.5 | N.R * | N.R *

  * *

Not Reported

##  IV Hardware setup and electrical characterization

To measure data from the fabricated pulse generator circuit we used a probe-
station and applied constant voltages and currents using a Semiconductor
Parameter Analyzer (see Fig. 7). The SPI and CLK signals of Fig. 7 were
generated using an Arbitrary Waveform Generator.

Figure 7: Schematic of the electrical setup that provides the biases and
signals to operate the pulse generator.

(a)

(b)

Figure 8: (8(a)) Difference between the measured and the simulated pulse time
width using the nominal conditions. The dotted line indicates the ideal
behavior. (8(b)) Error of the measured pulse time width after calibrating the
value of the capacitor.

(a)

(b)

(c)

(d)

Figure 9: Pre (red trace) and Post-calibrated (black dotted trace) pulse width
measurements (9(a)) 10 ms; (9(b)) 100 $\mu$s; (9(c)) 1 $\mu$s (9(d)) 336 ns

Initially, we biased the circuit using the same parameters used in the circuit
simulations to produce pulses of desired widths. However, the fabricated
circuits produced outputs that differ significantly from the circuit
simulations (Fig. 8(a)). The option of precisely calibrating the capacitance
value using the $C_{DAC}$ allows accurate pulse widths, with errors below $\pm
2.5\%$ (Fig. 8(b)). Here, only pulses down to 350 ns width could be measured,
due to the RC parasitic load induced by the setup. This is confirmed in Fig.
9, showing a selection of waveforms generated by the pulse generator under
different conditions. The trend shown in Fig. 8(a) proves that the pulse
generator can work at shorter pulse widths. In normal operating conditions,
the device to be tested is placed right next to the Pulse Generator, whose
photograph, as well as a detail of the contacts for the BEOL-integration of
ferroelectric devices, are shown in Fig. 6. In this configuration, the
integrated Pulse Generator is no longer limited by the parasitic load, and the
generated waveform are expected to be adequate to precisely investigate the
switching kinetics of ferroelectric devices.

##  V Conclusion

We presented an on-chip pulse-generator designed and fabricated in 180 nm
technology. We demonstrated pulse generation across 6 orders of magnitude with
an error of $\pm 2.5\%$, while capable of driving relatively high parasitic
capacitive loads present in the measurement setup. A low-parasitic on-chip
interconnect between the the pulse-generator and the FeFET/FTJ devices placed
on the same silicon ensures a reliable delivery of programming pulses with
high accuracy. Future work includes exploring performance of the two DACs and
their effect on the output pulses. This work is a key enabler for thorough
characterization of the FeFETs/FTJs, otherwise limited to the parasitics of
the measurement set up.

## References

  * [1] J. Okuno, T. Kunihiro, K. Konishi, H. Maemura, Y. Shuto, F. Sugaya, M. Materano, T. Ali, M. Lederer, K. Kuehnel, K. Seidel, U. Schroeder, T. Mikolajick, M. Tsukamoto, and T. Umebayashi, “High-endurance and low-voltage operation of 1t1c feram arrays for nonvolatile memory application,” in _2021 IEEE International Memory Workshop (IMW)_ , 2021, pp. 1–3.
  * [2] S. Beyer, S. Dünkel, M. Trentzsch, J. Müller, A. Hellmich, D. Utess, J. Paul, D. Kleimaier, J. Pellerin, S. Müller, J. Ocker, A. Benoist, H. Zhou, M. Mennenga, M. Schuster, F. Tassan, M. Noack, A. Pourkeramati, F. Müller, M. Lederer, T. Ali, R. Hoffmann, T. Kämpfe, K. Seidel, H. Mulaosmanovic, E. T. Breyer, T. Mikolajick, and S. Slesazeck, “Fefet: A versatile cmos compatible device with game-changing potential,” in _2020 IEEE International Memory Workshop (IMW)_ , 2020, pp. 1–4.
  * [3] B. Max, M. Hoffmann, S. Slesazeck, and T. Mikolajick, “Direct correlation of ferroelectric properties and memory characteristics in ferroelectric tunnel junctions,” _IEEE Journal of the Electron Devices Society_ , vol. 7, pp. 1175–1181, 2019.
  * [4] E. Covi, Q. T. Duong, S. Lancaster, V. Havel, J. Coignus, J. Barbot, O. Richter, P. Klein, E. Chicca, L. Grenouillet, A. Dimoulas, T. Mikolajick, and S. Slesazeck, “Ferroelectric tunneling junctions for edge computing,” in _2021 IEEE International Symposium on Circuits and Systems (ISCAS)_ , 2021, pp. 1–5.
  * [5] H. Mulaosmanovic, J. Ocker, S. Müller, M. Noack, J. Müller, P. Polakowski, T. Mikolajick, and S. Slesazeck, “Novel ferroelectric fet based synapse for neuromorphic systems,” in _2017 Symposium on VLSI Technology_ , 2017, pp. T176–T177.
  * [6] M. Halter, L. Bégon-Lours, V. Bragaglia, M. Sousa, B. J. Offrein, S. Abel, M. Luisier, and J. Fompeyrine, “Back-end, cmos-compatible ferroelectric field-effect transistor for synaptic weights,” _ACS Applied Materials & Interfaces_, vol. 12, no. 15, pp. 17 725–17 732, 2020, pMID: 32192333. [Online]. Available: https://doi.org/10.1021/acsami.0c00877
  * [7] B. Suresh, M. Bertele, E. T. Breyer, P. Klein, H. Mulaosmanovic, T. Mikolajick, S. Slesazeck, and E. Chicca, “Simulation of integrate-and-fire neuron circuits using hfo¡inf¿2¡/inf¿-based ferroelectric field effect transistors,” in _2019 26th IEEE International Conference on Electronics, Circuits and Systems (ICECS)_ , 2019, pp. 229–232.
  * [8] G. Indiveri and S.-C. Liu, “Memory and information processing in neuromorphic systems,” _Proceedings of the IEEE_ , vol. 103, no. 8, pp. 1379–1397, 2015\.
  * [9] H. Mulaosmanovic, J. Ocker, S. Müller, U. Schroeder, J. Müller, P. Polakowski, S. Flachowsky, R. van Bentum, T. Mikolajick, and S. Slesazeck, “Switching kinetics in nanoscale hafnium oxide based ferroelectric field-effect transistors,” _ACS Applied Materials & Interfaces_, vol. 9, no. 4, pp. 3792–3798, 2017, pMID: 28071052. [Online]. Available: https://doi.org/10.1021/acsami.6b13866
  * [10] M. Materano, P. D. Lomenzo, H. Mulaosmanovic, M. Hoffmann, A. Toriumi, T. Mikolajick, and U. Schroeder, “Polarization switching in thin doped hfo2 ferroelectric layers,” _Applied Physics Letters_ , vol. 117, no. 26, p. 262904, 2020.
  * [11] H. Mulaosmanovic, F. Müller, M. Lederer, T. Ali, R. Hoffmann, K. Seidel, H. Zhou, J. Ocker, S. Mueller, S. Dünkel, D. Kleimaier, J. Müller, M. Trentzsch, S. Beyer, E. T. Breyer, T. Mikolajick, and S. Slesazeck, “Interplay between switching and retention in hfo¡sub¿2¡/sub¿-based ferroelectric fets,” _IEEE Transactions on Electron Devices_ , vol. 67, no. 8, pp. 3466–3471, 2020.
  * [12] E. Covi, A. Cabrini, L. Vendrame, L. Bortesi, R. Gastaldi, and G. Torelli, “On-wafer analog pulse generator for fast characterization and parametric test of resistive switching memories,” _IEEE Transactions on Semiconductor Manufacturing_ , vol. 27, no. 2, pp. 134–150, 2014.
  * [13] N. Weste, _CMOS VLSI Design: A Circuits and Systems Perspective_. Addison Wesley, 2011.
  * [14] D. Allstot, “A precision variable-supply cmos comparator,” _IEEE Journal of Solid-State Circuits_ , vol. 17, no. 6, pp. 1080–1087, 1982.
  * [15] P. Allen and D. Holberg, _CMOS Analog Circuit Design_ , ser. The Oxford Series in Electrical and Computer Engineering. OUP USA, 2011.
  * [16] S. Lee and H. Yoo, “Configurable digital pulse generator for neuromorphic devices,” in _2021 International Conference on Electronics, Information, and Communication (ICEIC)_ , 2021, pp. 1–3.
  * [17] H. Erfanijazi, T. Serrano-Gotarredona, and B. Linares-Barranco, “Novel programmable single pulse generator for producing pulse widths in different time scales,” in _2021 International Conference on Content-Based Multimedia Indexing (CBMI)_ , 2021, pp. 1–6.

  *[CMOS]: Complementary Metal-Oxide-Semiconductor
  *[IDAC]: Current Digital to Analog Converter
  *[DAC]: Digital to Analog Converter


# Spin-pumping from a ferromagnetic insulator to an unconventional
superconductor
with interfacial Andreev bound-states

Chi Sun  Center for Quantum Spintronics, Department of Physics, Norwegian
University of Science and Technology, NO-7491 Trondheim, Norway  Jacob Linder
Center for Quantum Spintronics, Department of Physics, Norwegian
University of Science and Technology, NO-7491 Trondheim, Norway

###### Abstract

Spin-pumping from a ferromagnetic insulator into a high-$T_{c}$ superconductor
with a $d$-wave superconducting order parameter has recently been
experimentally observed. Such unconventional superconducting order is known to
produce interfacial bound-states for certain crystallographic orientations.
Here, we present a methodology which can be used to study spin-pumping into
unconventional superconductors, including the role of interfacial bound-
states. As an example, we determine how the crystallographic orientation of
the $d$-wave order parameter relative the interface changes the spin-pumping
effect. We find that the spin-pumping effect is slightly enhanced at low
temperatures for orientations hosting interfacial bound-states compared to
other superconducting states. However, the spin-pumping effect does not show a
coherence peak close to $T_{c}$ for such orientations, and instead remains
smaller than the normal state value for all $T$. For orientations not hosting
interfacial bound-states, we find that the pumped spin current can be
increased to several times relative the normal-state spin current at
frequencies that are small compared to the superconducting gap. Our results
show that the spin-pumping dependency on frequency and temperature changes
qualitatively depending on the crystallographic orientation of unconventional
superconducting order parameters relative the interface.

##  I Introduction

Spin-pumping [1, 2] is an effect where a spin current is injected from a
ferromagnet with a precessing magnetization into a material in contact with
the ferromagnet. Its reciprocal effect is spin-transfer torque [3], and the
two constitute central phenomena in the field of spintronics [4].

The precessional motion of the magnetization is typically induced by an
external static magnetic field that is misaligned with the magnetization. The
motion is then maintained by a transverse microwave field which supplies the
required energy to maintain the precession despite damping in the ferromagnet.
The damping originates from a spin-lattice coupling, which ultimately causes
conversion of energy to dissipated heat. The spin-pumping effect is
experimentally manifested as an increase in the damping due to the angular
momentum loss of the ferromagnet in the form of a spin current. The damping is
measurable as an increase in the linewidth of a ferromagnetic resonance (FMR)
[5] signal.

Spin-pumping by both ferromagnetic metals and ferromagnetic insulators (FIs)
into normal metals has been studied widely in the literature [6]. The injected
spin current increases with the frequency of the magnetization precession, and
creates a spin accumulation in the normal metal that eventually decays as one
moves away from the interface with the magnetic material. When the normal
metal is replaced by a superconductor, the spin-pumping effect is altered
considerably [7, 8, 9, 10, 11, 12]. Considering a conventional superconductor
such as Al or Pb, technically denoted a BCS $s$-wave superconductor, the gap
in the quasiparticle energy spectrum suppresses the spin-pumping effect at low
temperatures. However, when the temperature is close to the superconducting
critical temperature $T_{c}$, the spin pumping effect at low frequencies can
be strongly enhanced in the superconducting state [13, 14] compared to the
normal state due to the large quasiparticle density of states close to the gap
edge of the superconductor.

Spin-pumping is also possible via antiferromagnetic materials and has been
considered for antiferromagnets coupled with normal metals [15, 16] and, more
recently, superconductors [17]. The spin-pumping effect has similar features
as in the ferromagnetic case, with the exception that the FMR frequency which
matches the natural precession frequency of antiferromagnetic materials is
typically orders of magnitude larger than in ferromagnetic materials. This
difference stems from the very large exchange coupling between the magnetic
sublattices in an antiferromagnet.

Recently, spin-pumping has also been considered using unconventional
superconductors [18, 19], such as high-$T_{c}$ cuprates. The superconducting
order parameter in these materials can have a $d$-wave symmetry, meaning that
it is anisotropic in momentum space. Quasiparticles moving in different
directions in the lattice can thus experience both a different magnitude and
different sign of the superconducting order parameter. Experimentally, spin-
pumping in a ferromagnetic metal/$d$-wave superconductor bilayer was recently
reported [20]. The modified damping was measured in YBa2Cu3O7-d/Ni80Fe20
(YBCO/Py) structures, where one found a suppression of the damping upon
cooling the temperature sufficiently to transition into the superconducting
state of YBCO. When the temperature was lowered further, some samples showed
an upturn in the linewidth, corresponding to an increase in the damping and
thus more efficiently injected spin current in YBCO. This was interpreted as a
consequence of the crystallographic dependence of the $d$-wave superconducting
order parameter relative the Py interface: some orientations feature nodes in
the order parameter, where the quasiparticle do not have any gap in their
excitation spectrum. This should give rise to spin-pumping features closer to
the normal metal case. Importantly, some orientations also permit for the
existence of interfacial Andreev bound-states [21, 22]. These states exist at
the normal-state Fermi level (zero energy) when the $d$-wave superconductor is
interfaced with a non-magnetic material, like a normal metal. However, when
the superconductor (SC) is interfaced with a ferromagnet, it has been shown
that the bound-state energies are shifted to finite energies and that their
density of states is strongly suppressed [23]. Therefore, the low-temperature
behavior of the spin-pumping in a FI/$d$-wave SC cannot immediately be
explained by the formation of a large density of zero-energy interfacial
Andreev bound-states that occur when a $d$-wave SC is in contact with a non-
magnetic barrier.

The spin-pumping in FI/$d$-wave SC was considered theoretically [18] for a
fixed crystallographic orientation that does not feature any interfacial
Andreev bound-states. This effect was also studied in $p$-wave
superconductors, but in neither case was the effect of interfacial states
taken into account [19]. Here, we develop a theory for spin-pumping into
unconventional superconductors for arbitrary crystallographic orientations,
which thus includes the effect of interfacial bound-states. In particular, we
consider a $d$-wave superconducting order parameter
$\Delta=\Delta_{0}(T)\cos(2\theta-2\alpha)$ where $\Delta_{0}(T)$ is the
temperature-dependent magnitude of the superconducting order parameter,
$\theta$ is the angle describing the quasiparticle motion, and $\alpha$ is a
parameter that describes different crystallographic orientations of the
superconductor lattice relative the interface to the FI. For $\alpha=0$, the
order parameter has a $d_{x^{2}-y^{2}}$ wave form relative the interface, so
that quasiparticles normally incident toward the FI experience the maximum
value of the gap. For $\alpha=\pi/4$, the order parameter has a $d_{xy}$ wave
form relative the interface, so that quasiparticles normally incident toward
the FI experience a node of the gap. For $\alpha=0$, our results are
qualitatively consistent with the recent predictions of Ref. [18], although
some quantitative differences are present due to the fact that we consider a
specularly reflecting (clean) interface rather than a rough interface. At very
small frequencies compared to the superconducting gap, the spin current is
enhanced in the $d$-wave $\alpha=0$ case compared to the normal state. For
$\alpha=\pi/4$, we find that the spin-pumping is different compared to both
the $s$-wave and $d$-wave $\alpha=0$ cases. We find that at low temperatures
and small frequencies, the spin current is slightly enhanced compared to the
other superconducting pairing symmetries. Moreover, we find that the coherence
peak in the spin current close to $T_{c}$ vanishes in the $\alpha=\pi/4$ case,
so that the spin current is always smaller than that in the normal state of
the system. We discuss the physical origin of all our findings. Our theory can
be applied to other unconventional superconducting materials and accounts for
the presence of interfacial Andreev bound-states, making it a useful tool to
explore spin transport in unconventional superconducting junctions [24, 25,
26, 27, 28, 29].

##  II Theory

Below, we present the methodology used to compute the pumped spin current. It
is based on a BTK (Blonder-Tinkham-Klapwijk)-like [30] scattering approach,
but modified to account for a rotating magnetization which gives rise to
additional time-dependent factors in the wavefunctions. We provide detailed
expressions for the wavefunctions that participate in the scattering
processes, the relevant boundary conditions at the interface between the FI
and the SC, and explain how the pumped spin current can be computed. The FI
region occupies $x<0$ while the SC occupies $x>0$.

###  II.1 Expressions in the SC

Based on the BTK theory [30], the propagation of quasiparticles [electron
(hole)-like quasiparticles with spin up and down] in the SC can be described
by the Bogoliubov-de Gennes (BdG) equation

$\hat{H}_{S}\Psi_{S}=E_{s}\Psi_{S},$ (1)

where $\Psi_{S}$ is the eigenstate with energy eigenvalue $E_{s}$. The
Hamiltonian of the SC is given by

$\hat{H}_{S}=\begin{pmatrix}-{\frac{\triangledown^{2}}{2m_{e}}}-\mu&0&0&\Delta\\\
0&-{\frac{\triangledown^{2}}{2m_{e}}}-\mu&-\Delta&0\\\
0&-\Delta^{*}&{\frac{\triangledown^{2}}{2m_{e}}}+\mu&0\\\
\Delta^{*}&0&0&{\frac{\triangledown^{2}}{2m_{e}}}+\mu\end{pmatrix},$ (2)

where $m_{e}$ is the effective mass of the quasiparticles and $\mu$ is the
chemical potential. The superconducting gap is denoted as
$\Delta=\Delta_{0}(T)g(\theta)$, where $\Delta_{0}(T)$ is the temperature-
dependent gap amplitude and $g(\theta)$ describes the superconducting pair
symmetry. In this work, the conventional BCS temperature dependence
$\Delta_{0}(T)=\Delta_{0}\tanh(1.74\sqrt{T_{c}/T-1})$ is considered. By
solving the BdG equation as an eigenvalue problem, the wavefunction in the SC
can be obtained.

In a $s$-wave superconductor, the superconducting gap is isotropic, i.e.,
$\Delta=\Delta_{0}(T)$ with $g(\theta)=1$. Let us first consider the
wavefunction of an electron-like quasiparticle with spin up that is incoming
towards the interface, which has the form of

$\Psi_{S}=\left[\begin{pmatrix}u_{0}\\\ 0\\\ 0\\\
v_{0}\end{pmatrix}e^{-ik_{x+}x}+a_{1}\begin{pmatrix}u_{0}\\\ 0\\\ 0\\\
v_{0}\end{pmatrix}e^{ik_{x+}x}+b_{1}\begin{pmatrix}v_{0}\\\ 0\\\ 0\\\
u_{0}\end{pmatrix}e^{-ik_{x-}x}\right]e^{-\frac{iE_{s}t}{\hbar}}e^{ik_{y}y}+\left[c_{1}\begin{pmatrix}0\\\
u^{{}^{\prime}}_{0}\\\ -v^{{}^{\prime}}_{0}\\\
0\end{pmatrix}e^{ik_{x+}^{{}^{\prime}}x}+d_{1}\begin{pmatrix}0\\\
v^{{}^{\prime}}_{0}\\\ -u^{{}^{\prime}}_{0}\\\
0\end{pmatrix}e^{-ik_{x-}^{{}^{\prime}}x}\right]e^{-\frac{iE_{s}^{{}^{\prime}}t}{\hbar}}e^{ik_{y}y},$
(3)

where $a_{1}$, $b_{1}$, $c_{1}$ and $d_{1}$ are coefficients describing normal
reflection without spin-flip, Andreev reflection without spin-flip, normal
reflection with spin-flip, and Andreev reflection with spin-flip,
respectively. These coefficients can be determined by applying appropriate
boundary conditions, which we will get back to. To differentiate it from the
incident energy $E_{s}$, the energy after the spin-flip is denoted as
$E_{s}^{{}^{\prime}}$. In terms of $E_{s}$, we have the usual coherence
factors
$u_{0}=\sqrt{\frac{1}{2}(1+\frac{\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)}}{E_{s}})}$
and
$v_{0}=\sqrt{\frac{1}{2}(1-\frac{\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)}}{E_{s}})}$.
The $x$ components of the wave vectors are

$\displaystyle k_{x+}$
$\displaystyle=\frac{\sqrt{2m_{e}(\mu+\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)})-\hbar^{2}k_{y}^{2}}}{\hbar},$
$\displaystyle k_{x-}$
$\displaystyle=\frac{\sqrt{2m_{e}(\mu-\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)})-\hbar^{2}k_{y}^{2}}}{\hbar},$
(4)

which describe the electron-like and hole-like quasiparticles, respectively.
Similarly, $u_{0}^{{}^{\prime}}$, $v_{0}^{{}^{\prime}}$,
$k_{x+}^{{}^{\prime}}$ and $k_{x-}^{{}^{\prime}}$ have the same forms with
respect to $E_{s}^{{}^{\prime}}$. In terms of the Fermi wave vector
$k_{F}=\sqrt{2m_{e}\mu}/\hbar$ in the normal-state, the transverse component
of the wave vector $k_{y}=k_{F}\sin{\theta}$ is conserved across the
interface.

In a $d$-wave superconductor, the superconducting gap is anisotropic, i.e.,
$\Delta=\Delta_{0}(T)g(\theta)$ with $g(\theta)=\cos{(2\theta-2\alpha)}$, as
introduced in Sec. I. The wavefunction due to an incident electron-like
quasiparticle with spin up becomes

$\displaystyle\Psi_{S}$ $\displaystyle=\left[\begin{pmatrix}u_{-}\\\ 0\\\ 0\\\
v_{-}e^{-i\gamma_{-}}\end{pmatrix}e^{-ik_{x+,-}x}+a_{1}\begin{pmatrix}u_{+}\\\
0\\\ 0\\\
v_{+}e^{-i\gamma_{+}}\end{pmatrix}e^{ik_{x+,+}x}+b_{1}\begin{pmatrix}v_{-}e^{i\gamma_{-}}\\\
0\\\ 0\\\
u_{-}\end{pmatrix}e^{-ik_{x-,-}x}\right]e^{-\frac{iE_{s}t}{\hbar}}e^{ik_{y}y}$
(5) $\displaystyle+\left[c_{1}\begin{pmatrix}0\\\ u^{{}^{\prime}}_{+}\\\
-v^{{}^{\prime}}_{+}e^{-i\gamma_{+}}\\\
0\end{pmatrix}e^{ik_{x+,+}^{{}^{\prime}}x}+d_{1}\begin{pmatrix}0\\\
v^{{}^{\prime}}_{-}e^{i\gamma_{-}}\\\ -u^{{}^{\prime}}_{-}\\\
0\end{pmatrix}e^{-ik_{x-,-}^{{}^{\prime}}x}\right]e^{-\frac{iE_{s}^{{}^{\prime}}t}{\hbar}}e^{ik_{y}y},$

in which the coherence factors and wave vectors are modified compared with
those in the $s$-wave case due to the anisotropy. Depending on the
quasiparticle motion, the coherence factors are
$u_{\pm}=\sqrt{\frac{1}{2}(1+\frac{\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)g^{2}(\theta_{\pm})}}{E_{s}})}$
and
$v_{\pm}=\sqrt{\frac{1}{2}(1-\frac{\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)g^{2}(\theta_{\pm})}}{E_{s}})}$
with $\theta_{+}=\theta$ and $\theta_{-}=\pi-\theta$. Similarly, the wave
vectors become

$\displaystyle
k_{x+,\pm}=\frac{\sqrt{2m_{e}(\mu+\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)g^{2}(\theta_{\pm})})-\hbar^{2}k_{y}^{2}}}{\hbar},$
$\displaystyle
k_{x-,\pm}=\frac{\sqrt{2m_{e}(\mu-\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)g^{2}(\theta_{\pm})})-\hbar^{2}k_{y}^{2}}}{\hbar},$
(6)

in which the second $\pm$ sign in the subscript represents the propagation of
the quasiparticles along the $\pm x$ axis. In addition, the factor
$e^{i\gamma_{\pm}}=\frac{g(\theta_{\pm})}{|g(\theta_{\pm})|}$ is introduced.

Note that the energy-dependent wave-vectors and coherence factors appearing in
Eqs. (3) and (5) are only applicable for positive energies, i.e., $E_{s}>0$
and $E_{s}^{{}^{\prime}}>0$. When $E_{s}<0$, the following replacements should
be made: $k_{x\pm}\rightarrow k_{x\mp}$, $u_{0}\rightarrow-v_{0}^{*}$ and
$v_{0}\rightarrow u_{0}^{*}$ for $s$-wave and $k_{x\pm,\pm}\rightarrow
k_{x\mp,\pm}$, $u_{\pm}\rightarrow-v_{\pm}^{*}$ and $v_{\pm}\rightarrow
u_{\pm}^{*}$ for $d$-wave. Similar modifications apply for
$E_{s}^{{}^{\prime}}<0$. A detailed explanation regarding the negative energy
wavefunctions can be found in the Appendix.

###  II.2 Expressions in the FI

In the FI, the Hamiltonian for electron-like and hole-like quasiparticles has
the form

$\hat{H}_{F}=\begin{pmatrix}\underline{H}_{F}&0\\\
0&-\underline{H}_{F}^{*}\end{pmatrix}$ (7)

with

$\underline{H}_{F}=-{\frac{\triangledown^{2}}{2m_{e}}}+U+J\underline{\sigma}\cdot\textbf{M}(t),$
(8)

in which $\underline{\bm{\sigma}}$ denotes the Pauli matrix vector and
$\underline{\ldots}$ is notation for a 2$\times$2 matrix. The potential $U$ is
larger than the Fermi energy $\mu$ in the nearby SC. $J$ decribes the exchange
interaction in the ferromagnet between the localized spin magnetization and
the itinerant electrons. The magnetization is defined as

$\textbf{M}(t)=(m\cos\Omega t,m\sin\Omega t,\sqrt{1-m^{2}}),$ (9)

where $m$ is the magnetization amplitude and $\Omega$ denotes the FMR
frequency for spin-pumping. Here the microwave energy for FMR is considered as
small compared with the exchange energy, i.e., $\hbar\Omega\ll J$. By
employing a wavefunction with the structure
$e^{-\frac{iEt}{\hbar}}(e^{-\frac{i\Omega t}{2}},e^{\frac{i\Omega
t}{2}},e^{\frac{i\Omega t}{2}},e^{-\frac{i\Omega t}{2}})^{T}$ for its time-
dependence , the non-stationary Schrödinger equation can be solved as an
eigenvalue problem. The four eigenpairs are obtained as: $E_{1}=E_{+}$ with
$(a_{+},b_{+},0,0)^{T}$, $E_{2}=E_{-}$ with $(a_{-},b_{-},0,0)^{T}$,
$E_{3}=-E_{+}$ with $(0,0,a_{+},b_{+})^{T}$ and $E_{4}=-E_{-}$ with
$(0,0,a_{-},b_{-})^{T}$. To linear order of the adiabaticity parameter
$\beta=\frac{\hbar\Omega}{2J}$, the eigenenergies are given by

$E_{\pm}=U+\frac{\hbar^{2}(k_{x}^{2}+k_{y}^{2})}{2m_{e}}\pm
J(1-\beta\sqrt{1-m^{2}}).$ (10)

The corresponding eigenstates are described by the coefficients

$a_{+}=-\frac{(1-\beta)(1+\sqrt{1-m^{2}})-m}{\sqrt{[(1-\beta)(1+\sqrt{1-m^{2}})-m]^{2}+[(1+\beta)(1-\sqrt{1-m^{2}})-m]^{2}}},$
(11)
$b_{+}=\frac{(1+\beta)(1-\sqrt{1-m^{2}})-m}{\sqrt{[(1-\beta)(1+\sqrt{1-m^{2}})-m]^{2}+[(1+\beta)(1-\sqrt{1-m^{2}})-m]^{2}}},$
(12)
$a_{-}=-\frac{(1+\beta)(1-\sqrt{1-m^{2}})+m}{\sqrt{[(1-\beta)(1+\sqrt{1-m^{2}})+m]^{2}+[(1+\beta)(1-\sqrt{1-m^{2}})+m]^{2}}},$
(13)
$b_{-}=\frac{(1-\beta)(1+\sqrt{1-m^{2}})+m}{\sqrt{[(1-\beta)(1+\sqrt{1-m^{2}})+m]^{2}+[(1+\beta)(1-\sqrt{1-m^{2}})+m]^{2}}}.$
(14)

Based on the above, the total wavefunction in the FI is constructed as

$\Psi_{F}=\left[m_{1}\begin{pmatrix}a_{+}e^{\frac{-i\Omega t}{2}}\\\
b_{+}e^{\frac{i\Omega t}{2}}\\\ 0\\\
0\end{pmatrix}e^{-ik_{F1}x}e^{\frac{-iE_{1}t}{\hbar}}+n_{1}\begin{pmatrix}a_{-}e^{\frac{-i\Omega
t}{2}}\\\ b_{-}e^{\frac{i\Omega t}{2}}\\\ 0\\\
0\end{pmatrix}e^{-ik_{F2}x}e^{\frac{-iE_{2}t}{\hbar}}+p_{1}\begin{pmatrix}0\\\
0\\\ a_{+}e^{\frac{i\Omega t}{2}}\\\ b_{+}e^{\frac{-i\Omega
t}{2}}\end{pmatrix}e^{-ik_{F3}x}e^{\frac{-iE_{3}t}{\hbar}}+q_{1}\begin{pmatrix}0\\\
0\\\ a_{-}e^{\frac{i\Omega t}{2}}\\\ b_{-}e^{\frac{-i\Omega
t}{2}}\end{pmatrix}e^{-ik_{F4}x}e^{\frac{-iE_{4}t}{\hbar}}\right]e^{ik_{y}y},$
(15)

where $m_{1}$, $n_{1}$, $p_{1}$ and $q_{1}$ are transmission coefficients to
be determined by applying appropriate boundary conditions. In order to match
the time-dependence of the wavefunction components on the FI and SC sides, we
can obtain $E_{s}^{{}^{\prime}}=E_{s}-\hbar\Omega$ and
$E_{1}=E_{2}=E_{3}=E_{4}=E_{s}-\frac{\hbar\Omega}{2}$. In terms of $E_{s}$,
the corresponding $x$ component of the four wave vectors in the FI are
expressed as

$k_{F1}=\frac{\sqrt{2m_{e}[E_{s}-U-J(1-\beta\sqrt{1-m^{2}}+\beta)]-\hbar^{2}k_{y}^{2}}}{\hbar},$
(16)
$k_{F2}=\frac{\sqrt{2m_{e}[E_{s}-U+J(1-\beta\sqrt{1-m^{2}}-\beta)]-\hbar^{2}k_{y}^{2}}}{\hbar},$
(17)
$k_{F3}=\frac{\sqrt{2m_{e}[-E_{s}-U-J(1-\beta\sqrt{1-m^{2}}-\beta)]-\hbar^{2}k_{y}^{2}}}{\hbar}$
(18)

and

$k_{F4}=\frac{\sqrt{2m_{e}[-E_{s}-U+J(1-\beta\sqrt{1-m^{2}}+\beta)]-\hbar^{2}k_{y}^{2}}}{\hbar}.$
(19)

Note that all wave numbers in the FI possess imaginary values since a large
potential $U$ is required to ensure the ferromagnet to be insulating.

###  II.3 Boundary conditions and spin pumping current

The boundary conditions at the FI/SC interface (i.e., $x=0$) are

$\Psi_{F}\big{|}_{x=0}=\Psi_{S}\big{|}_{x=0},$ (20)
$\partial_{x}\Psi_{S}\big{|}_{x=0}-\partial_{x}\Psi_{F}\big{|}_{x=0}=2Zk_{F}\Psi_{F}\big{|}_{x=0},$
(21)

which describe the continuity of the wavefunction and its derivative
appropriate for the $\delta$ function potential at the interface,
respectively. $Z$ is a dimensionless parameter which characterizes the
interface transparency.

Figure 1: (Color online) Temperature-dependence of the pumped spin current
into the superconductor. The spin current is normalized on its normal-state
value and shown for the $s$-wave (left), $d$-wave $\alpha=0$ (middle), and
$d$-wave $\alpha=\pi/4$ (right) cases for several frequencies $\Omega$.

The longitudinal spin-pumping current polarized along the $z$ axis can be
calculated as

$I_{s}=\int_{-W_{l}}^{W_{h}}\int_{-\theta_{c}}^{\theta_{c}}j_{sz}(E_{s},\theta)f_{0}(E_{s})N_{s}\cos{\theta}d\theta
dE_{s},$ (22)

where the quantum mechanical spin current
$j_{sz}=-\frac{\hbar^{2}}{2m_{e}}{\Im\mathrm{m}}[\Psi_{S}^{\dagger}\partial_{x}\mathbf{\hat{\sigma_{z}}}\Psi_{S}]$
with $\mathbf{\hat{\sigma_{z}}}=\begin{pmatrix}\underline{\sigma}_{z}&0\\\
0&\underline{\sigma}_{z}^{*}\end{pmatrix}$ produced by a given incoming
particle. $f_{0}(E_{s})$ denotes the Fermi-Dirac distribution and
$N_{s}={\mathrm{Re}}[\frac{|E_{s}|}{\sqrt{E_{s}^{2}-\Delta^{2}}}]$ is the
density of states in the SC. Consider the specific case of an incoming spin-up
electron-like quasiparticle from the SC. To ensure that the $x$-component of
the wave vectors (before and after spin-flip) in the SC are real and
propagating for normal incidence (i.e., $\theta=0$), the upper and lower
integration limit for $E_{s}$ are determined by
$W_{h}=\sqrt{\mu^{2}+\Delta_{0}^{2}(T)}$ and
$W_{l}=\sqrt{\mu^{2}+\Delta_{0}^{2}(T)}-\hbar\Omega$, in which the
$\hbar\Omega$ shift originates from the energy relation
$E_{s}^{{}^{\prime}}=E_{s}-\hbar\Omega$. Based on this energy range, a
critical angle
$\theta_{c}=\arcsin{\sqrt{1-\frac{\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)}}{\mu}}}$
is utilized for the same reason. However, in practice a much smaller
integration range for $E_{s}$ can be employed to greatly enhance the numerical
integration efficiency for the SC when $T\ll T_{c}$. As explained in the
Appendix, a temperature-dependent integration range may be utilized in
practice. In all plots, we set the barrier parameter to $Z=2$. Varying $Z$
does not give rise to any qualitatively new features, since the ferromagnetic
region is insulating and thus total reflection (either normal or Andreev) of
incoming particles.

Eq. (22) covers the spin-pumping current induced by the electron-like
quasiparticle with spin up (i.e., electron-up) incident. The total spin
current is given by the sum of electron-up and electron-down contributions. In
the Appendix, we present the wavefunctions in the SC induced by electron-down,
hole-up and hole-down incidents. Based on these wavefunctions, it is shown
that the spin-pumping current can also be calculated as the sum of hole-up,
hole-down and additional background contributions (see Appendix for details).

##  III Results and Discussion

###  III.1 Temperature dependence

In Fig. 1, we investigate the temperature dependence of the pumped spin
current into the SC. Due to the existence of interfacial Andreev bound-states
in N/$d$-wave or I/$d$-wave junctions in the $\alpha=\pi/4$ case, one might
expect that the spin-pumping effect would be strongly enhanced compared to the
$s$-wave or $d$-wave $\alpha=0$ case. However, our results do not show any
such strong enhancement. To understand these results, it is instructive to
consider the quantum transport of a ferromagnet/insulator/$d$-wave
superconductor junction as considered in Ref. [23]. With increasing
polarization of the ferromagnet, the density of interfacial bound-states is
reduced and shifted away from the Fermi level. Moreover, it was shown in Ref.
[23] that the interfacial states do not contribute to the spin conductance of
the junction, whereas they contribute strongly to the charge conductance. This
can be understood from the fact that the interfacial states are confined to a
distance of order superconducting coherence length from the interface, and are
thereafter converted into a Cooper pair charge supercurrent propagating into
the bulk of the SC. However, the Cooper pairs in the bulk of the $d$-wave SC
cannot carry a spin current since they are in a singlet state, which means
that the interfacial bound-states cannot contribute to spin transport into the
SC.

Figure 2: (Color online) Frequency-dependence of the pumped spin current into
the superconductor. The spin current is normalized on its normal-state value
and shown for the $s$-wave (left), $d$-wave $\alpha=0$ (middle), and $d$-wave
$\alpha=\pi/4$ (right) cases for several temperatures $T$.

This is consistent with our results. In the $\alpha=\pi/4$ case, there is only
a slight increase in the spin-pumping current at low temperatures and small
frequencies $\hbar\Omega/\Delta_{0}\ll 1$ compared to the $s$-wave and
$d$-wave $\alpha=0$ case. Rather than attributing this to the effect of
interfacial bound-states, this small increase in current should instead stem
from the nodal positioning of the gap relative the interface. Normally
incident particles feel no gap, i.e. $\Delta=0$, in the $\alpha=\pi/4$ case,
so that the spin current is not blocked by the excitation spectrum gap unlike
the $s$-wave and $d$-wave $\alpha=0$ cases.

The second qualitative difference in the $\alpha=\pi/4$ case is the absence of
the coherence peak in the spin-pumping current for temperatures $T\simeq
T_{c}$. Such a peak appears in the $s$-wave at low frequencies
$\hbar\Omega\ll\Delta_{0}$. The reason is that close to $T_{c}$, the
superconducting gap is suppressed, which makes it possible for spin-flip
excitations of quasiparticles with energy $\hbar\Omega$ to occur. Such
excitations can carry a spin current. For energies right above the
superconducting gap, the density of quasiparticle states is strongly enhanced
(the so-called coherence peak). This large number of available states is the
reason for the superconducting spin current exceeding the normal state spin
current at low frequencies, since in the normal metal case there is a range of
energies $\hbar\Omega$ around the Fermi level which contribute to the pumped
spin current. As $\Omega$ increases, the relative increase in the
superconducting spin current vs. normal spin current is lost since the normal
state contribution grows more rapidly. A similar reasoning applies to the
$d$-wave $\alpha=0$ case, although the peaked spin current close to $T_{c}$ is
smaller than in the $s$-wave case due to the nodes in the anisotropic $d$-wave
gap. As we will discuss later, this becomes different for very small
frequencies compared to the superconducting gap, which are experimentally the
most straightforward ones to obtain.

The difference between the two cases described above and the $d$-wave
$\alpha=\pi/4$ case which causes the coherence peak at $T\simeq T_{c}$ to be
absent, is the nodal orientation of the gap relative the interface.
Quasiparticle trajectories close to normal incidence $\theta=0$ will
contribute the most to spin transport. The gap is very small close to normal
incidence and indeed vanishes exactly at $\theta=0$. Therefore, the coherence
peak effect described for the $s$-wave case is present already at low
temperatures (due to the smallness of the gap) for quasiparticle trajectories
close to $\theta=0$. It is not present for trajectories away from normal
incidence, since the gap eventually recovers its full value at $\theta=\pi/4$.
This is the reason for why the spin current for $\alpha=\pi/4$ is slightly
enhanced at low frequencies compared to the other order parameter cases. As
the temperature gets closer to $T_{c}$, however, there is no large pile-up of
available quasiparticle states close to normal incidence that suddenly become
accessible and boost the superconducting spin current, unlike in the $s$-wave
and $d$-wave $\alpha=0$ cases. Therefore, the spin current in the
superconducting case always remains smaller than the normal state value for
the $\alpha=\pi/4$ crystallographic orientation.

###  III.2 Frequency dependence

We now consider the dependence of the ratio between the pumped spin current in
the superconducting and normal-state on frequency in Fig. 2. At $\Omega=0$,
the spin current vanishes for all systems since there is no magnetization
precession and thus no pumped spin current. We note that for temperatures
approaching $T_{c}$, there is an abrupt increase in the the ratio of the spin
current in the superconducting state and normal state vs. frequency as soon as
$\Omega$ becomes non-zero in the $s$-wave and $d$-wave $\alpha=0$ cases. This
means that the spin-pumping is much more efficient in the superconducting
state compared to the normal state. This can be understood from the argument
provided in the previous section: at small frequencies
$\hbar\Omega\ll\Delta_{0}$, the large density of quasiparticles that can carry
spin current with energies close to the gap edge in the SC becomes available
and exceeds the normal-state spin current for small frequencies.
Interestingly, the figure shows that the spin current in the anisotropic
$d$-wave $\alpha=0$ case is several times larger than the normal-state current
for low frequencies $\hbar\Omega/\Delta_{0}\ll 1$ at $T/T_{c}\simeq 0.5.$ This
cannot be explained by the nodal quasiparticle transport in the $d$-wave case
compared to the isotropic $s$-wave gap, since such an effect would make the
spin current closer to its normal-state value. Instead, this enhancement in
the $d$-wave case stems from the pileup of quasiparticle density of states
near energies $E_{s}=\Delta_{0}(T)g(\theta)$. Since the gap magnitude varies
with $\theta$, and is generally smaller or equal to $\Delta_{0}$,
quasiparticles are more easily accessible via thermal excitations for small
frequencies.

In ferromagnets, the FMR resonance frequencies typically lies in the 1 GHz -
100 GHz range. A conventional BCS $s$-wave SC usually has a gap of order 1
meV, corresponding to roughly 250 GHz, whereas the gap magnitude in
high-$T_{c}$ $d$-wave SCs can be ten times as large. It is clear that
$\hbar\Omega/\Delta_{0}\ll 1$ for such parameters, suggesting that the large
enhancement of the spin current in the superconducting state relative the
normal state that takes place exactly for $\hbar\Omega/\Delta_{0}\ll 1$ is
accessible experimentally. Two remarks are in order here. First, it is
possible to access higher frequencies relative the gap,
$\hbar\Omega/\Delta_{0}\sim 1$, by considering spin-pumping into proximity-
induced superconducting regions. In effect, by considering spin-pumping from a
FI into a normal metal where the latter is proximitized by a SC, the normal
metal may effectively be considered as a SC with a smaller gap
$\Delta_{\text{eff}}$ compared to its host SC. Secondly, larger ratios of
$\hbar\Omega/\Delta_{0}$ are also attainable by replacing the FI with an
antiferromagnetic insulator [17]. In antiferromagnets (AFMs), the FMR
resonance frequency is of order THz, and recent theory suggests that the spin-
pumping dependence on temperature and frequency is similar in an AFI/SC as in
an FI/SC bilayer [17].

###  III.3 Crystallographic axis dependence

In Fig. 3, we consider exclusively the $d$-wave case and vary continuously the
misalignment $\alpha$ between the antinodal direction of the gap and the
interface normal. Experimentally, this corresponds to varying the
crystallographic orientation of how the superconducting sample is grown on top
of the FI. We see that the $\alpha=\pi/4$, corresponding to the orientation
which hosts zero-energy interfacial states in the absence of ferromagnetism,
provides the smallest spin current for moderate to high temperatures relative
$T_{c}$. Only at low temperatures $T\ll T_{c}$, the spin current for
$\alpha=\pi/4$ is comparable to other crystallographic orientations. This
result is consistent with the explanation given in Sec. III.1, which shows why
the $\alpha=\pi/4$ direction does not provide any spin current enhancement
relative the normal-state as the temperature approaches $T_{c}$. The fact that
$\alpha=\pi/4$ hosts interfacial bound-states, pinned to zero-energy for any
angles of incidence $\theta$ of the quasiparticles in the absence of
magnetism, but shifted to finite energies in the presence of magnetism, does
not help boosting the spin current. As mentioned previously, these states are
confined to close to the interface region and cannot be converted into a
Cooper pair current carrying the spin in the bulk of the superconductor.

##  IV Summary

Motivated by a recent experiment demonstrating spin-pumping from a
ferromagnetic insulator into a high-$T_{c}$ superconductor with a $d$-wave
superconducting order parameter, we have presented a methodology which can be
used to study spin-pumping into unconventional superconductors, including the
role of interfacial bound-states. To be concrete, we have focused on how the
crystallographic orientation of a $d$-wave superconducting order parameter
relative the interface changes the spin-pumping effect. Such unconventional
superconductivity is known to produce interfacial bound-states for certain
crystallographic orientations. We have shown that the spin-pumping effect is
slightly enhanced at low temperatures for orientations hosting interfacial
bound-states in the absence of ferromagnetism compared to other
superconducting states. However, the spin-pumping effect did not show a
coherence peak close to $T_{c}$ for such orientations, and instead remains
smaller than the normal state value for all $T$. For orientations not hosting
interfacial bound-states, we found that the pumped spin current can be
increased to several times the normal-state spin current at frequencies that
are small compared to the superconducting gap. Our results show that the spin-
pumping dependency on frequency and temperature changes qualitatively
depending on the crystallographic orientation of unconventional
superconducting order parameters relative the interface.

Figure 3: (Color online) Spin current in the superconducting state normalized
against the normal-state value of the spin current vs. crystallographic axis
orientation of the $d$-wave superconducting material relative the interface.
The $\alpha=0$ case corresponds to $d_{x^{2}-y^{2}}$-wave pairing relative the
interface, whereas $\alpha=\pi/4$ corresponds to $d_{xy}$-wave pairing
(featuring zero-energy interface states in the absence of ferromagnetism). We
have set $\hbar\Omega/\Delta_{0}=0.1.$

###### Acknowledgements.

We thank J. A. Ouassou for useful discussions. We acknowledge funding via the
Research Council of Norway Grant numbers 323766, as well as through its
Centres of Excellence funding scheme QuSpin, project number 262633.

##  Appendix A Complete wavefunctions

In the main text, we gave a detailed expression for the wavefunction for
incoming electron-like quasiparticles with spin up (referred to as electron-
up) in the superconductor. Here, we provide the corresponding expressions for
the other possible incoming quasiparticles. For notation simplicity, we use
$e\downarrow$, $h\uparrow$ and $h\downarrow$ to represent the electron-down,
hole-up and hole-down incidents, respectively. In a $s$-wave SC, the
wavefunctions are

$\Psi_{S,e\downarrow}=\left[\begin{pmatrix}0\\\ u_{0}\\\ -v_{0}\\\
0\end{pmatrix}e^{-ik_{x+}x}+a_{2}\begin{pmatrix}0\\\ u_{0}\\\ -v_{0}\\\
0\end{pmatrix}e^{ik_{x+}x}+b_{2}\begin{pmatrix}0\\\ v_{0}\\\ -u_{0}\\\
0\end{pmatrix}e^{-ik_{x-}x}\right]e^{-\frac{iE_{s}t}{\hbar}}e^{ik_{y}y}+\left[c_{2}\begin{pmatrix}u^{{}^{\prime}}_{0}\\\
0\\\ 0\\\
v^{{}^{\prime}}_{0}\end{pmatrix}e^{ik_{x+}^{{}^{\prime}}x}+d_{2}\begin{pmatrix}v^{{}^{\prime}}_{0}\\\
0\\\ 0\\\
u^{{}^{\prime}}_{0}\end{pmatrix}e^{-ik_{x-}^{{}^{\prime}}x}\right]e^{-\frac{iE_{s}^{{}^{\prime}}t}{\hbar}}e^{ik_{y}y},$
(23) $\Psi_{S,h\uparrow}=\left[\begin{pmatrix}0\\\ v_{0}\\\ -u_{0}\\\
0\end{pmatrix}e^{ik_{x-}x}+a_{3}\begin{pmatrix}0\\\ v_{0}\\\ -u_{0}\\\
0\end{pmatrix}e^{-ik_{x-}x}+b_{3}\begin{pmatrix}0\\\ u_{0}\\\ -v_{0}\\\
0\end{pmatrix}e^{ik_{x+}x}\right]e^{-\frac{iE_{s}t}{\hbar}}e^{ik_{y}y}+\left[c_{3}\begin{pmatrix}v^{{}^{\prime}}_{0}\\\
0\\\ 0\\\
u^{{}^{\prime}}_{0}\end{pmatrix}e^{-ik_{x-}^{{}^{\prime}}x}+d_{3}\begin{pmatrix}u^{{}^{\prime}}_{0}\\\
0\\\ 0\\\
v^{{}^{\prime}}_{0}\end{pmatrix}e^{ik_{x+}^{{}^{\prime}}x}\right]e^{-\frac{iE_{s}^{{}^{\prime}}t}{\hbar}}e^{ik_{y}y},$
(24) $\Psi_{S,h\downarrow}=\left[\begin{pmatrix}v_{0}\\\ 0\\\ 0\\\
u_{0}\end{pmatrix}e^{ik_{x-}x}+a_{4}\begin{pmatrix}v_{0}\\\ 0\\\ 0\\\
u_{0}\end{pmatrix}e^{-ik_{x-}x}+b_{4}\begin{pmatrix}u_{0}\\\ 0\\\ 0\\\
v_{0}\end{pmatrix}e^{ik_{x+}x}\right]e^{-\frac{iE_{s}t}{\hbar}}e^{ik_{y}y}+\left[c_{4}\begin{pmatrix}0\\\
v^{{}^{\prime}}_{0}\\\ -u^{{}^{\prime}}_{0}\\\
0\end{pmatrix}e^{-ik_{x-}^{{}^{\prime}}x}+d_{4}\begin{pmatrix}0\\\
u^{{}^{\prime}}_{0}\\\ -v^{{}^{\prime}}_{0}\\\
0\end{pmatrix}e^{ik_{x+}^{{}^{\prime}}x}\right]e^{-\frac{iE_{s}^{{}^{\prime}}t}{\hbar}}e^{ik_{y}y}.$
(25)

In a $d$-wave SC, we have

$\displaystyle\Psi_{S,e\downarrow}$ $\displaystyle=\left[\begin{pmatrix}0\\\
u_{-}\\\ -v_{-}e^{-i\gamma_{-}}\\\
0\end{pmatrix}e^{-ik_{x+,-}x}+a_{2}\begin{pmatrix}0\\\ u_{+}\\\
-v_{+}e^{-i\gamma_{+}}\\\
0\end{pmatrix}e^{ik_{x+,+}x}+b_{2}\begin{pmatrix}0\\\ v_{-}e^{i\gamma_{-}}\\\
-u_{-}\\\
0\end{pmatrix}e^{-ik_{x-,-}x}\right]e^{-\frac{iE_{s}t}{\hbar}}e^{ik_{y}y}$
(26) $\displaystyle+\left[c_{2}\begin{pmatrix}u^{{}^{\prime}}_{+}\\\ 0\\\ 0\\\
v^{{}^{\prime}}_{+}e^{-i\gamma_{+}}\end{pmatrix}e^{ik_{x+,+}^{{}^{\prime}}x}+d_{2}\begin{pmatrix}v^{{}^{\prime}}_{-}e^{i\gamma_{-}}\\\
0\\\ 0\\\
u^{{}^{\prime}}_{-}\end{pmatrix}e^{-ik_{x-,-}^{{}^{\prime}}x}\right]e^{-\frac{iE_{s}^{{}^{\prime}}t}{\hbar}}e^{ik_{y}y},$
$\displaystyle\Psi_{S,h\uparrow}$ $\displaystyle=\left[\begin{pmatrix}0\\\
v_{+}e^{i\gamma_{+}}\\\ -u_{+}\\\
0\end{pmatrix}e^{ik_{x-,+}x}+a_{3}\begin{pmatrix}0\\\ v_{-}e^{i\gamma_{-}}\\\
-u_{-}\\\ 0\end{pmatrix}e^{-ik_{x-,-}x}+b_{3}\begin{pmatrix}0\\\ u_{+}\\\
-v_{+}e^{-i\gamma_{+}}\\\
0\end{pmatrix}e^{ik_{x+,+}x}\right]e^{-\frac{iE_{s}t}{\hbar}}e^{ik_{y}y}$ (27)
$\displaystyle+\left[c_{3}\begin{pmatrix}v^{{}^{\prime}}_{-}e^{i\gamma_{-}}\\\
0\\\ 0\\\
u_{-}^{{}^{\prime}}\end{pmatrix}e^{-ik_{x-,-}^{{}^{\prime}}x}+d_{3}\begin{pmatrix}u_{+}^{{}^{\prime}}\\\
0\\\ 0\\\
v^{{}^{\prime}}_{+}e^{-i\gamma_{+}}\end{pmatrix}e^{ik_{x+,+}^{{}^{\prime}}x}\right]e^{-\frac{iE_{s}^{{}^{\prime}}t}{\hbar}}e^{ik_{y}y},$
$\displaystyle\Psi_{S,h\downarrow}$
$\displaystyle=\left[\begin{pmatrix}v_{+}e^{i\gamma_{+}}\\\ 0\\\ 0\\\
u_{+}\end{pmatrix}e^{ik_{x-,+}x}+a_{4}\begin{pmatrix}v_{-}e^{i\gamma_{-}}\\\
0\\\ 0\\\ u_{-}\end{pmatrix}e^{-ik_{x-,-}x}+b_{4}\begin{pmatrix}u_{+}\\\ 0\\\
0\\\
v_{+}e^{-i\gamma_{+}}\end{pmatrix}e^{ik_{x+,+}x}\right]e^{-\frac{iE_{s}t}{\hbar}}e^{ik_{y}y}$
(28) $\displaystyle+\left[c_{4}\begin{pmatrix}0\\\
v^{{}^{\prime}}_{-}e^{i\gamma_{-}}\\\ -u^{{}^{\prime}}_{-}\\\
0\end{pmatrix}e^{-ik_{x-,-}^{{}^{\prime}}x}+d_{4}\begin{pmatrix}0\\\
u^{{}^{\prime}}_{+}\\\ -v^{{}^{\prime}}_{+}e^{-i\gamma_{+}}\\\
0\end{pmatrix}e^{ik_{x+,+}^{{}^{\prime}}x}\right]e^{-\frac{iE_{s}^{{}^{\prime}}t}{\hbar}}e^{ik_{y}y}.$

Similarly, by matching the time-dependence of each component of the
wavefunctions on both FI and SC sides, we can obtain that
$E_{s}^{{}^{\prime}}=E_{s}-\hbar\Omega$ for $h\downarrow$ incident but
$E_{s}^{{}^{\prime}}=E_{s}+\hbar\Omega$ for $e\downarrow$ and $h\uparrow$
incidents. Note that the above wave functions are only applicable for positive
energies, i.e., $E_{s}>0$ and $E_{s}^{{}^{\prime}}>0$. For negative energy
wavefunctions, the energy-dependent wave vectors and coherence factors should
be modified, as indicated in the main text.

##  Appendix B Negative energy wavefunctions in the SC

The wavefunctions in the SC are obtained by solving the BdG equation
$\hat{H}_{S}\Psi_{S}=E_{s}\Psi_{S}$, in which $E_{s}$ can be either positive
or negative. It is found that if $(\psi_{1},0,0,\psi_{2})^{T}$ is an
eigenstate for $E_{s}=\sqrt{\xi^{2}+\Delta^{2}}$, then
$(\psi_{1},0,0,\psi_{2})^{T}$ is an eigenstate for
$E_{s}=-\sqrt{\xi^{2}+\Delta^{2}}$ with
$\xi=\frac{\hbar^{2}\textbf{k}^{2}}{2m_{e}}-\mu$.

Let us consider the wavefunction in a $s$-wave SC describing an electron-up
incident as an example. For the positive eigenenergy
$E_{s}=\sqrt{\xi^{2}+\Delta^{2}}$, the corresponding solution of the
eigenstate are

$\psi_{1}=\sqrt{\frac{1}{2}(1+\frac{\xi}{E_{s}})},\
\psi_{2}=\sqrt{\frac{1}{2}(1-\frac{\xi}{E_{s}})}.$ (29)

There are two possible choices of $\xi$:

$\xi=\sqrt{E_{s}^{2}-\Delta_{0}(T)^{2}}:\ \psi_{1}=u_{0},\ \psi_{2}=v_{0},\
k_{x+}=\frac{\sqrt{2m_{e}(\mu+\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)})-\hbar^{2}k_{y}^{2}}}{\hbar},$
(30) $\xi=-\sqrt{E_{s}^{2}-\Delta_{0}(T)^{2}}:\ \psi_{1}=v_{0},\
\psi_{2}=u_{0},\
k_{x-}=\frac{\sqrt{2m_{e}(\mu-\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)})-\hbar^{2}k_{y}^{2}}}{\hbar}.$
(31)

On the other hand, for the negative eigenenergy
$E_{s}=-\sqrt{\xi^{2}+\Delta^{2}}$, the corresponding solution of the
eigenstate changes to

$\psi_{1}=-\sqrt{\frac{1}{2}(1+\frac{\xi}{E_{s}})}^{*},\
\psi_{2}=\sqrt{\frac{1}{2}(1-\frac{\xi}{E_{s}})}^{*}.$ (32)

As before, there are two possible choices of $\xi$. We have

$\xi=-\sqrt{E_{s}^{2}-\Delta_{0}(T)^{2}}:\ \psi_{1}=-v_{0}^{*},\
\psi_{2}=u_{0}^{*},\
k_{x-}=\frac{\sqrt{2m_{e}(\mu-\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)})-\hbar^{2}k_{y}^{2}}}{\hbar},$
(33) $\xi=\sqrt{E_{s}^{2}-\Delta_{0}(T)^{2}}:\ \psi_{1}=-u_{0}^{*},\
\psi_{2}=v_{0}^{*},\
k_{x+}=\frac{\sqrt{2m_{e}(\mu+\sqrt{E_{s}^{2}-\Delta_{0}^{2}(T)})-\hbar^{2}k_{y}^{2}}}{\hbar}.$
(34)

Based on the above, the negative energy version of Eq. (3) in the main text is
modified to

$\Psi_{S}=\left[\begin{pmatrix}-v_{0}^{*}\\\ 0\\\ 0\\\
u_{0}^{*}\end{pmatrix}e^{-ik_{x-}x}+a_{1}\begin{pmatrix}-v_{0}^{*}\\\ 0\\\
0\\\ u_{0}^{*}\end{pmatrix}e^{ik_{x-}x}+b_{1}\begin{pmatrix}u_{0}^{*}\\\ 0\\\
0\\\
-v_{0}^{*}\end{pmatrix}e^{-ik_{x+}x}\right]e^{-\frac{iE_{s}t}{\hbar}}e^{ik_{y}y}+\left[c_{1}\begin{pmatrix}0\\\
-v^{{}^{\prime}*}_{0}\\\ -u^{{}^{\prime}*}_{0}\\\
0\end{pmatrix}e^{ik_{x-}^{{}^{\prime}}x}+d_{1}\begin{pmatrix}0\\\
u^{{}^{\prime}*}_{0}\\\ v^{{}^{\prime}*}_{0}\\\
0\end{pmatrix}e^{-ik_{x+}^{{}^{\prime}}x}\right]e^{-\frac{iE_{s}^{{}^{\prime}}t}{\hbar}}e^{ik_{y}y}$
(35)

when both $E_{s}$ and $E_{s}^{{}^{\prime}}$ are negative. If either $E_{s}$ or
$E_{s}^{{}^{\prime}}$ is negative, only the wave vectors and coherence factors
regarding the negative energy should be modified. A similar procedure can be
applied to the $d$-wave SC and the corresponding modifications for the
negative energy wavefunctions are indicated in the main text.

Figure 4: (Color online) Plot of the energy-dependence of the integrands
determining the pumped spin current, demonstrating the qualitative difference
in spectral spin current for the different superconducting order parameter
symmetries. The upper row corresponds to incoming spin$-\uparrow$ electron-
like quasiparticles whereas the lower row corresponds to spin$-\downarrow$
electron-like quasiparticles. The $s$-wave case is shown in the left column,
$d$-wave $\alpha=0$ in the middle column, $d$-wave $\alpha=\pi/4$ in the right
column. Note how the energy-integrands are inverted and shifted by
$\hbar\Omega$ for the two spin species. We have used $T/T_{c}=0.1$ and
$\hbar\Omega/\Delta_{0}=0.5.$

##  Appendix C Total spin pumping current

As discussed in the main text, the spin-pumping current induced by an incoming
electron-up ($e\uparrow$) particle is calculated by Eq. (22), in which $E_{s}$
is integrated from $-W_{l}$ to $W_{h}$ to ensure the wave vectors in the SC
are real and propagating with $E_{s}^{{}^{\prime}}=E_{s}-\hbar\Omega$. The
same integration range applies for spin current induced by $h\downarrow$. As
for the spin current based on $e\downarrow$ and $h\uparrow$ incident
particles, $E_{s}$ should be integrated from $-W_{h}$ to $W_{l}$ since the
energy relation changes to $E_{s}^{{}^{\prime}}=E_{s}+\hbar\Omega$.

In Eq. (22), the quantum mechanical spin-pumping current
$j_{sz}=-\frac{\hbar^{2}}{2m_{e}}{\Im\mathrm{m}}[\Psi_{S}^{\dagger}\partial_{x}\mathbf{\hat{\sigma_{z}}}\Psi_{S}]$
is introduced. Based on the complete wavefunctions in Appendix A and B, it can
be seen that the hole and electron-induced quantum mechanical spin currents
fulfill the relation

$j_{sz,h\uparrow(\downarrow)}(E_{s},\theta)=-j_{sz,e\uparrow(\downarrow)}(-E_{s},\pi-\theta).$
(36)

Let us first consider the resulting spin-pumping current induced by
$e\uparrow$ and $h\uparrow$, i.e.

$I_{s,e\uparrow}=\int^{\theta_{c}}_{-\theta_{c}}\int_{-W_{l}}^{W_{h}}\cos\theta
j_{sz,e\uparrow}(E_{s},\theta)f_{0}(E_{s})N_{s}dE_{s}d\theta,$ (37)
$I_{s,h\uparrow}=\int^{\theta_{c}}_{-\theta_{c}}\int_{-W_{h}}^{W_{l}}\cos\theta
j_{sz,h\uparrow}(E_{s},\theta)f_{0}(E_{s})N_{s}dE_{s}d\theta.$ (38)

Inserting Eq. (36) into Eq. (38) and using $f_{0}(-E_{s})=1-f_{0}(E_{s})$,
$j_{sz}(\theta)=j_{sz}(\theta+\pi)$ for any type of incoming particle, and
$N_{s}(E_{s})=N_{s}(-E_{s})$, we have

$I_{s,e\uparrow}=\frac{1}{2}(I_{s,e\uparrow}+I_{s,h\uparrow}+I_{\text{bg}}),$
(39)

in which the background contribution
$I_{\text{bg}}=\int^{\theta_{c}}_{-\theta_{c}}\int_{-W_{l}}^{Wh}\cos\theta
j_{sz,e\uparrow}(E_{s},\theta)N_{s}dE_{s}d\theta$. Similar relations apply for
incident spin-$\downarrow$ particles. Therefore, there are two equivalent ways
of computing the spin current. The first way considers only electrons and no
hole excitations. The second way considers both electrons and holes, in which
case a background spin current contribution $I_{\text{bg}}$ must also be
included. Consequently, we may calculate the total spin current as the sum of
electron-up and electron-down contributions without involving the hole
incidents, i.e., $I_{s}=I_{s,e\uparrow}+I_{s,e\downarrow}$. The appearance of
the background contribution can be understood most easily by considering a
normal metal with second-quantized Hamiltonian
$\hat{H}=\sum_{{\bm{k}}\sigma}\epsilon_{\bm{k}}c_{{\bm{k}}\sigma}^{\dagger}c_{{\bm{k}}\sigma}$.
Here, $\epsilon_{{\bm{k}}}$ is the energy dispersion relation and
$c,c^{\dagger}$ are fermion annihilation/creation operators. Using that
$c_{{\bm{k}}\sigma}^{\dagger}c_{{\bm{k}}\sigma}=1-c_{{\bm{k}}\sigma}c_{{\bm{k}}\sigma}^{\dagger}$,
we can rewrite the free electron gas form above to a combined electron and
hole gas form plus a background contribution consisting of all electron states
being filled:

$\displaystyle\hat{H}=\frac{1}{2}\sum_{{\bm{k}}\sigma}\epsilon_{\bm{k}}c_{{\bm{k}}\sigma}^{\dagger}c_{{\bm{k}}\sigma}-\frac{1}{2}\sum_{{\bm{k}}\sigma}\epsilon_{\bm{k}}c_{{\bm{k}}\sigma}c_{{\bm{k}}\sigma}^{\dagger}+\frac{1}{2}\sum_{{\bm{k}}\sigma}\epsilon_{\bm{k}}.$
(40)

The last term on the right hand side of Eq. (40) is consistent with the
background spin current term $I_{\text{bg}}$ due to the absence of a Fermi-
Dirac distribution function in its expression.

##  Appendix D Energy-dependence of spin current integrand

It is instructive to consider the energy-resolved quantum mechanical spin
current for the $s$-wave and $d$-wave SCs, respectively, and see how the low-
energy behavior changes depending on $\alpha$. We define the quantity
$i(E_{s})=\frac{2m_{e}}{\hbar^{2}}\int^{\theta_{c}}_{-\theta_{c}}j_{sz}(E_{s},\theta)N_{s}\cos\theta
d\theta$. In Fig. 4, we plot the integrand $i(E_{s})$ determining the pumped
spin current as a function of energy $E_{s}$ for different superconducting
order parameter symmetries. Comparing the integrands corresponding to
$e\uparrow$ and $e\downarrow$ quasiparticles shown in the upper and lower
rows, an energy shift of precisely $\hbar\Omega/\Delta_{0}$ between them can
be observed. This energy shift also appears in the different energy
integration ranges for $e\uparrow$ (i.e., $[-W_{l},W_{h}]$) and $e\downarrow$
(i.e., $[-W_{h},W_{l}]$) incident particles with $W_{h}=W_{l}+\hbar\Omega$. In
addition, due to the specific form of the density of states $N_{s}$ and
coherence factors of the SC, it can be seen that the main contributions of the
integrands well below $T_{c}$ are within a small energy range around the gap
edge. Meanwhile, the sum of integrands of $e\uparrow$ and $e\downarrow$
approaches zero for large $E_{s}$ values. Therefore, it is more efficient to
use smaller energy integration range to cover the integrated spin current and
reduce possible numerical error induced by the extremely sharp peaks appearing
in Fig. 4. On the other hand, a larger energy integration range should be
considered when the temperature approaches $T_{c}$. For $T>T_{c}$, the full
energy range $[-W_{l},W_{h}]$ for $e\uparrow$ and $[-W_{h},W_{l}]$ for
$e\downarrow$ must be considered since the SC then behaves as NM with zero
gap.

## References

  * [1] Y. Tserkovnyak, A. Brataas, and G. E. W. Bauer Phys. Rev. Lett. 88, 117601 (2002)
  * [2] S. Takahashi, Physical Principles of Spin Pumping. In: Xu, Y., Awschalom, D., Nitta, J. (eds) Handbook of Spintronics. Springer, Dordrecht (2016).
  * [3] D. C. Ralph and M. D. Stiles, J. Mag. Mag. Mater. 320, 1190 (2008).
  * [4] I. Zutic, J. Fabian, ans S. Das Sarma, Rev. Mod. Phys. 76, 323 (2004).
  * [5] C. Kitte, Phys. Rev. 73, 115 (1948).
  * [6] Y. Tserkovnyak, A. Brataas, G. E. W. Bauer, and B. I. Halperin, Rev. Mod. Phys. 77, 1375 (2005).
  * [7] C. Bell, S. Milikisyants, M. Huber, and J. Aarts, Phys. Rev. Lett. 100, 047002 (2008).
  * [8] J. P. Morten, A. Brataas, G. E. W. Bauer, W. Belzig, and Y. Tserkovnyak, Europhys. Lett. 84, 57008 (2008).
  * [9] K.-R. Jeon, C. Ciccarelli, A. J. Ferguson, H. Kurebayashi, L. F. Cohen, X. Montiel, M. Eschrig, J. W. A. Robinson and M. G. Blamire, Nature Mater. 17, 499 (2018).
  * [10] K.-R. Jeon, C. Ciccarelli, H. Kurebayashi, J. Wunderlich, L. F. Cohen, S. Komori, J. W. A. Robinson, and M. G. Blamire, Phys. Rev. Applied 10, 014029 (2018).
  * [11] Y. Yao, Q. Song, Y. Takamura, J. P. Cascales, W. Yuan, Y. Ma, Y. Yun, X. C. Xie, J. S. Moodera, and W. Han, Phys. Rev. B 97, 224414 (2018).
  * [12] M. Umeda, Y. Shiomi, T. Kikkawa, T. Niizeki, J. Lustikova, S. Takahashi, and E. Saitoh, App. Phys. Lett. 112, 232601 (2018).
  * [13] M. Inoue, M. Ichioka, and H. Adachi, Phys. Rev. B 96, 024414 (2017).
  * [14] T. Kato, Y. Ohnuma, M. Matsuo, J. Rech, T. Jonckheere, T. Martin, Phys. Rev. B 99, 144411 (2019)
  * [15] R. Cheng, J. Xiao, Q. Niu, and A. Brataas, Phys. Rev. Lett. 113, 057601 (2014)
  * [16] P. Vaidya, S. A. Morley, J. V. Tol, Y. Liu, R. Cheng, A. Brataas, D. Lederman, and E. del Barco, Science 368, 160 (2020)
  * [17] E. H. Fyhn and J. Linder, Phys. Rev. B 103, 134508 (2021).
  * [18] Y. Ominato, A. Yamakage, T. Kato, M. Matsuo, Phys. Rev. B 105, 205406 (2022).
  * [19] Y. Ominato, A. Yamakage, and M. Matsuo, Phys. Rev. B 106, L161406 (2022).
  * [20] S. J. Carreira, D. Sanchez-Manzano, M.-W. Yoo, K. Seurre, V. Rouco, A. Sander, J. Santamaria, A. Anane, and J. E. Villegas, Phys. Rev. B 104, 144428 (2021).
  * [21] Y. Tanaka ans S. Kashiwaya, Phys. Rev. Lett. 74, 3451 (1995).
  * [22] C.-R. Hu, Phys. Rev. Lett. 72, 1526 (1994).
  * [23] S. Kashiwaya, Y. Tanaka, N. Yoshida, and M. R. Beasley, Phys. Rev. B 60, 3572 (1999).
  * [24] I. Eremin, F. S. Nogueira, and R.-J. Tarento Phys. Rev. B 73, 054507 (2006).
  * [25] M. S. Grønsleth, J. Linder, J.-M. Børven, and A. Sudbø Phys. Rev. Lett. 97, 147002 (2006)
  * [26] P. M. R. Brydon Phys. Rev. B 80, 224520 – Published 22 December 2009; Erratum Phys. Rev. B 85, 099905 (2012)
  * [27] Y. Tanaka, T. Yokoyama, A. V. Balatsky, and N. Nagaosa Phys. Rev. B 79, 060505(R) (2009).
  * [28] R. Takashima, S. Fujimoto, and T. Yokoyama, Phys. Rev. B 96, 121203(R) (2017)
  * [29] L. G. Johnsen, H. T. Simensen, A. Brataas, and J. Linder, Phsy. Rev. Lett. 127, 207001 (2021).
  * [30] G. E. Blonder, M. Tinkham, and T. M. Klapwijk, Phys. Rev. B 25, 4515 (1982).