# Pharmacokinetic Measurements in Dose Finding Model Guided by Escalation with
Overdose Control

Arnab Kumar Maity∗, Satrajit Roy Chowdhury, Ray Li, Lada Markovtsova, Roberto
Bugarini
<EMAIL_ADDRESS>
(22 May 2022)

###### Abstract

Oncology drug development starts with a dose escalation phase to find the
maximal tolerable dose (MTD). Dose limiting toxicity (DLT) is the primary
endpoint for dose escalation phase. Traditionally, model-based dose escalation
trial designs recommend a dose for escalation based on an assumed dose-DLT
relationship. Pharmacokinetic (PK) data are often available but are currently
only used by clinical teams in a subjective manner to aid decision making.
Formal incorporation of PK data in dose-escalation models can make the
decision process more efficient and lead to an increase in precision. In this
talk we present a Bayesian joint modeling framework for incorporating PK data
in Oncology dose escalation trials. This framework explores the dose-PK and
PK-DLT relationships jointly for better model informed dose escalation
decisions. Utility of the proposed model is demonstrated through a real-life
case study along with simulation.

Keywords: BLRM, DLT, MTD, Phase I, PK, Oncology.

##  1 Introduction

The early phase of a cancer clinical trial or an oncology clinical trial is
unique by the virtue of phenomenon that the participants who are enrolled are
not the healthy volunteers. The primary reason pertaining with this is severe
toxicities generated as side effects of the investigating molecules which
includes myelosuppression, immune suppression, nausea and vomiting, anorexia,
alopecia and mucositis (Crowley and Hoering,, 2012). While a certain label of
toxicity is acceptable and can not be avoided in order for the drug to show
anti-tumor activities, severe toxicity is not desirable and this can be
controlled by determining the appropriate dose of the investigating drug.
Consequently, the primary objective of an early oncology trial is to find the
highest dose which is well tolerated by the participants. This is often
referred as the maximum tolerated dose (MTD) (Ji et al.,, 2010).

Several designs now exist to guide the dose escalation trial and each design
has large literature associated with it. We refer the readers to Geller,
(1984); Rosenberger and Haines, (2002); Hummel et al., (2009); Zhou et al.,
(2018) for a review of different methods in the literature. For a recent and
comprehensive review see Daimon et al., (2019). These designs can generally be
classified in two categories – conventional methods or rule based designs and
model based designs Rosenberger and Haines, (2002). The conventional methods
estimate MTD from the observed data while the latter approach sets the MTD as
a parameter of the model and thereby tries to estimate from the data.

The widely used rule based design is 3 + 3 design and there exist other
variants of it. These are also referred as the trial and error method or up
and down design due the nature of these designs. On the other hand, continuous
reassessment method (CRM) (O’Quigley et al.,, 1990), escalation with overdose
control (EWOC) (Babb et al.,, 1998), and modified toxicity probability
intervals (mTPI) (Ji et al.,, 2010) are prominent examples of the model based
designs. It can be shown that CRM and EWOC are consistent in estimating the
MTD (O’Quigley and Shen,, 1996; Tighiouart and Rogatko,, 2010). Our
concentration in this paper centers around the model based design approach. In
particular, in this article we develop our method to follow the EWOC
principle.

When applying EWOC criteria, it can be assumed that the dose toxicity
relationship is logistic (Zacks et al.,, 1998). In other words, the
probability that an individual will observe severe toxicity which would
restrict the subject to receive the drug further is a logistic function of the
dose. A dose finding model has generally been agreed to be adaptive such that
it can accommodate all data until the current point of time and can drive the
decision in an informative manner and Bayesian methodologies are particularly
suitable for it. In addition, the non Bayesian parameter estimates often
degenerate toward the boundary regions if there is no heterogeneity observed
in the data which is often the case in the beginning of the trial (O’Quigley
and Shen,, 1996). Toward this end, Neuenschwander et al., (2008) proposed the
Bayesian logistic regression model (BLRM) and provided some guidance on the
prior elicitation. Neuenschwander et al., (2015) further advocated and
extended the mechanism for combination dose finding using the single agent
studies as the prior data.

The logistic regression model combined with EWOC criteria can be shown as
optimal Bayesian feasible design a concept which was introduced by Babb et
al., (1998) and will be illustrated further in Section 2. Furthermore, a
desirable property of a good design is that the design would be coherent.
Tighiouart and Rogatko, (2010) proved that EWOC scheme satisfies the
conditions of a coherent design. In fact, any relationship between dose and
toxicity which is monotone can be shown to follow the coherent design
properties which will be further illustrated in Section 2.

Historically it has been argued that any dose finding study has potential to
become more effective in determining the MTD when the available
pharmacokinetic (PK) data is considered in dose decision strategy Collins et
al., (1986, 1990). Since, the information of PK activity in vivo is an useful
indicator in quantifying the drug amount applied, currently the Phase I
Oncology studies use the PK data in a qualitative manner to aid the dose
finding decisions, for example, see Shaw et al., (2014). However,
incorporating PK data in the dose response model has a potential to better
estimate the MTD (Takeda et al.,, 2018).

Nevertheless, attempts had been made to utilize the PK summary in developing
the model based dose escalation approach. History dates back to Collins et
al., (1986) who suggested that a the MTD finding could be reached by
considering smaller dose levels with a preclinically obtained pharmacokinatic
guided design; since then different proposals of escalation designs involving
PK data have been made in the literature; see Collins et al., (1990); Newell,
(1994); Holford, (1995); Piantadosi and Liu, (1996); Whitehead et al., (2001,
2007); Cotterill et al., (2015); Ursino et al., (2017); Takeda et al., (2018)
and references therein. An R (R Core Team,, 2019) package dfpk (Toumazi et
al.,, 2018) have been developed to fit a PK guided Bayesian dose finding
model. Furthermore, Takeda et al., (2018) argued measuring the exposure data
in terms of pharmacokinetics and taking them into consideration in finding the
MTD is critical because there exist different dose exposure relationship such
as linear and nonlinear, in different dose levels. Moreover, higher area under
concentration (AUC) have been observed in those participants in the trial who
experience dose limiting toxicity – a dose which restricts participants to
receive any further dose, than those who do not experience this event.
Consequently, the thirst of incorporating PK data directly in the dose finding
model is becoming more essential.

In this paper, we propose a shared parameter joint model approach which relies
on EWOC principle and is presented in Figure 1 schematically. The joint
parameter consists of two components – first model measures the dose exposure
relationship and the second model quantifies the probability of experiencing a
DLT as a function of PK exposure data. While we propose a Bayesian Log Normal
regression to model the dose exposure relationship, the second model is fitted
using a Bayesian logistic regression model and hence the inference is based on
the usual EWOC criteria. The primary purpose of this approach is twofold –
first, the Bayesian logistic regression method has been shown to better
estimate the entire dose toxicity curve (Cotterill et al.,, 2015) and one
hopes to borrow the same feature here. Second, use of EWOC principle enables
us to retain the current dose escalation practice which is straightforward and
easy to implement using existing software routines such as East (Corporation,,
1992). Moreover, it is shown that under mild conditions the proposed approach
is coherent and consistent in MTD estimation.

Figure 1: A schematic representation of the proposed BLRM-PK joint model.

Cotterill et al., (2015) proposed a similar two parts models in combining PK
data in dual combination model, however, their model is not a joint model in a
statistical sense and hence the consideration of PK data can only be utilized
in setting the dose escalation rule and in setting the trial stopping rule.
Ursino et al., (2017) discussed various dose escalation models which
incorporates PK data formally, however, their dose allocation rule is based on
probabilistic computation. Conceptually, our proposed method is similar in
spirit of the model formulation and the model fitting with that in Takeda et
al., (2018). Even though they used the summary level PK data such as AUC which
will be used in this article, they adjusted the departure of AUC values for a
particular individual by considering a standardization over AUC measurements;
in the contrary, we tackle this problem imposing a stochastic distribution
assuming a Log Normal random variability. Furthermore, we retain the current
practice of next safe dose declaration using EWOC principle.

The rest of the article is organized as follows. Section 2 formalizes the
model and the method, discusses the parameter estimations, describes some
theoretical properties of our proposed approach. In Section 3 we apply our
proposed methodology in two Phase I trial examples. In Section 4 we provide
some simulation studies and show that incorporating PK data is beneficial in
terms of dose allocation. Essentially, the simulation studies show that more
participants receive the doses which belong to the target toxicity interval
than the case when PK measurements are not considered in the model. Finally in
Section 5 we provide the concluding remarks and outline some future
developments.

##  2 Methodology

###  2.1 Joint Model

Suppose $n$ subjects are dosed $d_{i},i=1,\ldots,n$ at different stages of a
dose finding trial. While the response DLT $y_{i}$ is binary taking values 1
or 0 for an individual is measured, the overall objective is to estimate the
MTD. Suppose the pharmacokinetic exposure data is summarized by $x_{i}$ for
the $i$-th individual. The following models then can be used to describe the
dose exposure relationship.

$\mu_{i}=\log(\gamma_{0})+\gamma_{1}\log(d_{i}/d^{*}),i=1,\ldots,n,$ (1)

where

$x_{i}\sim LN(\mu_{i},\sigma^{2}),i=1,\ldots,n.$ (2)

where, $d^{*}$ is a reference dose which is used for standardization purpose.
A guidance is given how to choose the value of $d^{*}$ in Section 2.2. A
random variable $X$ is said to follow a Log Normal distribution
$LN(\mu,\sigma^{2})$ if and only if $\log(X)$ follows the Normal distribution
$N(\mu,\sigma^{2})$. We assume Log Normal distribution here because from a
practical point of view any summary of PK measurements has support on the
positive part of the real line. In addition, logarithmic transformations are
considered here to bring stability in the computations.

The second part of the joint model borrows the widely popular BLRM method
(Neuenschwander et al.,, 2015) idea to predict the probability of experiencing
a DLT in a given subject given PK originated from the given dose. The usual
logistic regression model suffices to quantify this probability,

$\pi_{i}=\frac{\exp\\{\log\alpha+\beta\log(x_{i})\\}}{1+\exp\\{\log\alpha+\beta\log(x_{i})\\}},i=1,\ldots,n,$
(3)

with $\pi_{i}$ being the probability of experiencing DLT, that is, $y_{i}$
follows a Binary distribution with probability $\pi_{i}$. In essence, together
the models (1), (2), and (3) jointly can be fitted to the accumulated data to
infer dose escalation decision.

Note that, this joint BLRM-PK model will be referred as the BLRM-PK model
hereafter. In addition, for the sake of prior specification convenience, we
work with the logarithmic transformations of the coefficient parameters
$\gamma_{1}$ and $\beta$. Then models (1) and (3) become,

$\mu_{i}=\log(\gamma_{0})+\log(\gamma_{1})\log(d_{i}/d^{*}),i=1,\ldots,n,$ (4)

and

$\pi_{i}=\frac{\exp\\{\log\alpha+\log(\beta)\log(x_{i})\\}}{1+\exp\\{\log\alpha+\log(\beta)\log(x_{i})\\}},i=1,\ldots,n,$
(5)

respectively. Working with logarithm of the coefficients also have the
interpretations such as dose level increases the exposure increases inside of
the body and hence probability of observing a DLT increases. If the trial
operators think that this assumption may not be tenable then it can be
accommodated accordingly.

###  2.2 Prior Elicitation

It has generally been agreed that BLRM-PK model is more useful after a trial
accumulates at least some amount of data (personal communication with clinical
pharmacology scientists). Toward this end, in a Bayesian setting after the
trial starts to observe data then the data should dictate the dose
recommendation by outweighing the prior. In this sense, a noninformative prior
could serve the purpose, however, is not recommended primarily for two
reasons. First, under a noninformative prior the DLT rates take a U-shape when
no data is considered in the analysis, which is not desirable. Second, the
Bayesian MCMC often observes divergence under a noninformative prior with a
high uncertainty, that is, with low precision. To this direction, a weakly
informative prior has been suggested in Neuenschwander et al., (2015). We
follow similar approach here in eliciting an appropriate to carry out the
Bayesian analysis.

We set a Bivariate Normal prior for $(\log\alpha,\log,\beta)$ with mean
(logit(0.33), 0), standard deviation (2, 1) and correlation coefficient 0,
where $\text{logit}(x)=\log(\frac{x}{1-x})$. In addition, another Bivariate
Normal prior is specified for $\log\gamma_{0},\log\gamma_{1}$ with mean (0,
0), standard deviation (2, 1), and correlation coefficient 0. To complete the
prior analysis a Log Normal prior with mean 0.25 and standard deviation 0.35
can be specified on $\sigma^{2}$. We select a low value for the mean for
convergence purpose and the standard deviation 0.35 is obtained by
$\log(2)/1.96$. As it will be shown in Section 3, this prior specification
provides a wide range coverage of the DLT rates on the interval (0, 1). It can
also be shown that this prior elicitation provides the anticipated DLT rate
0.33 at the dose $d^{*}$.

###  2.3 Inference and Decisions

For a dose escalation (or de-escalation) decision at any given point of time
during a clinical we follow the process described by Neuenschwander et al.,
(2015). Hence the DLT rates at any dose can be categorized into three
subgroups implying two extreme intervals are not intended as they either
reflect an ineffective dose or they reflect a toxic dose. Formally, it is said
that $i$-th participant is under dosed if $\pi_{i}<0.16$ and the participant
is over dosed if $\pi_{i}>0.33$. As a consequence one wants the more
participants receive doses having DLT rate in between 0.16 and 0.33. In other
words, it is a good practice of any Phase I trial that the dose finding method
would maximize the probability $\pi_{i}\in(0.16,0.33)$. As a result, a good
dose escalation design would be tailored by a good matrix of target toxicity
(TT) probability which often ranges between 0.60 and 0.80 (Neuenschwander et
al.,, 2015). As a byproduct, we can also compute the probability of under
dosing (UD) and probability of overdosing (OD) which are targeted to be
minimized in a good design.

The posterior probability of target toxicity can be computed by integrating
the posterior DLT rate over the target toxicity region with respect to the
posterior distributions, that is,

$\Pr(\pi_{d}\in(0.16,0.33)|\bm{x},\bm{y},\alpha,\beta,\gamma_{0},\gamma_{1})=\int_{\pi_{d}\in(0.16,0.33)}\Pr(\bm{\pi}|\bm{x},\bm{y},\alpha,\beta,\gamma_{0},\gamma_{1})d\bm{\pi}d{\alpha}d{\beta}d{\gamma_{0}}d{\gamma_{1}},$

where $\Pr(\pi_{d}|.)$ is the posterior DLT rate at dose $d$.

In a Bayesian design these probabilities can be obtained from the Markov Chain
Monte Carlo (MCMC) as the exact analytical form of the posterior distribution
of $\pi_{i}$ is not available and we use JAGS (Plummer,, 2003) and rjags
(Plummer,, 2019). In particular, the following approximation can be used to
estimate the above integral:

$\Pr(\pi_{d}\in(0.16,0.33)|\bm{x},\bm{y},\alpha,\beta,\gamma_{0},\gamma_{1})=\frac{1}{M}\sum_{m=1}^{M}\Pr(\bm{\pi}^{(m)}|\bm{x},\bm{y},\alpha^{(m)},\beta^{(m)},\gamma_{0}^{(m)},\gamma_{1}^{(m)})I(\pi_{d}^{(m)}\in(0.16,0.33)),$

where $M$ is the total number of MCMC samples and
$\alpha^{(m)},\beta^{(m)},\gamma_{0}^{(m)},\gamma_{1}^{(m)}$ are the posterior
estimates of $\alpha,\beta,\gamma_{0},\gamma_{1}$ at $m$-th iteration
respectively. In a very similar fashion, one can compute the probability of
underdosing and probability overdosing.

Recall that the primary objective of the trial is to estimate the MTD, denoted
by $\xi$ hereafter, and we define the proportion of subjects receiving MTD as,

$\theta=\Pr(y_{i}=1|d=\xi),i=1,\ldots,n.$

In our formulation, the sequence of dose levels can be shown to be Bayesian-
feasible at level $(1-\alpha)$ (see Section 2.4). Following Babb et al.,
(1998) we set $\alpha$ at 0.25 meaning that at any stage of the dose finding
trial a next dose will be statistically recommended if the overdose
probability at the next dose is less than 0.25. This fully defines the
Bayesian adaptive and sequential dose escalation (or de-escalation) process
which will be followed in this article.

###  2.4 Statistical Properties

####  2.4.1 Consistency

In what follows we use similar notations as in Tighiouart and Rogatko, (2010).
Specifically, we assume that there exist $d^{1}$ and $d^{2}$, $d^{1}<d^{2}$,
such that

$\Pr(y_{i}=1|x_{i},d_{i}=d^{1})=0,i=1,\ldots,n,$
$\Pr(y_{i}=1|x_{i},d_{i}=d^{2})=1-\epsilon,i=1,\ldots,n,$

where $0<\epsilon<1$ and $\theta<1-\epsilon$.

Moreover, we consider the dose toxicity relationship of the form

$\Pr(y_{i}=1|x_{i},d_{i})=F\biggl{(}F^{-1}(1-\epsilon)+\psi\bigl{(}\frac{d_{i}-d^{1}}{d^{2}-d^{1}}\bigr{)}\biggr{)},$
(6)

where $0<\psi^{1}\leq\psi\leq\psi^{2}$ for some $\psi^{1}$ and $\psi^{2}$.
This model assumes that the quantiles of $F$ are linear in the log-standarized
dose $z=\log[(d-d^{1})/(d^{2}-d)]$. In our setup $F=F_{1}*F_{2}$, where
$F_{1}(x)=e^{x}/(1+e^{x})$, $F_{2}(.)$ is an identity map, and $F$ is a
convolution of $F_{1}$ and $F_{2}$.

Now letting $\phi=F^{-1}(1-\epsilon)-F^{-1}(\theta)$ the MTD can be shown to
be equal to $-\phi/\psi$ on the log standardized scale. Moreover, let
$G(x)=F(F^{-1}(\theta)+\phi+x),g(x)G^{\prime}(x),$ and $\log d_{1}=-\phi/\psi$
be the log dose administered to the first participant. Then for subject $i$,
log dose $\log d_{i}\in
L^{*}=[-\frac{\phi}{\psi^{1}},-\frac{\phi}{\psi^{2}}].$ Also let $h(\psi)$ be
a prior density function for the parameter $\psi$ on $[\psi^{1},\psi^{2}]$.
Then the sequence of dose levels ${d_{n}}$ such that $\Pr(\log
d_{n}\leq-\phi/\psi|\bm{x},\bm{y})\geq 1-\alpha$ for all $n\geq 2$ is called
Bayesian feasible at level $(1-\alpha)$. Now, using Theorem 2.1 in Tighiouart
and Rogatko, (2010) it can be shown that the sequence of dose levels
considering our formulation will be Bayesian-feasible at level $(1-\alpha)$.

####  2.4.2 Coherence

A desirable property of good design, primarily because of ethical concern, is
that to allow a dose escalation if there is no DLT observed in the current
dose level and that to allow a dose de-escalation if there is a DLT observed
in the current dose level. This property is known as coherent property and
Cheung, (2005) argued how a practical adaptive model based approach such as
CRM is designed to satisfy this property.

Suppose after reparameterization in terms of the MTD $\xi$,
$F(d_{i},\xi)=\Pr(y_{i}=1|x_{i},d_{i})$ is the model given in (6). The BLRM-PK
method will be said to be coherent in escalation if $d_{n}\geq d_{n-1}$
whenever $y_{n}-1=0$, and it will be coherent in de-escalation if $d_{n}\leq
d_{n-1}$ whenever $y_{n-1}=1$ for all $n=2,3,\ldots$. In essence, the BLRM-PK
design will be coherent if it is coherent in both dose escalation and dose de-
escalation. The following corollary then establishes that the joint BLRM-PK
model preserves this coherence principle.

######  Corollary 1.

If $F(d,\xi)$ is nonincreasing in $\xi$ for fixed dose $d$, then the BLRM-PK
method described in Sections 2.1 and 2.4 is coherent.

By virtue of constructions in (4) and (5), $F(d,\xi)$ is in fact nonincreasing
in $\xi$ for fixed dose $d$ and hence the proof follows by similar argument in
Tighiouart and Rogatko, (2010).

##  3 Real Data Example

In this Section we apply our proposed methodology in two real world data
originated from two clinical trials and carry out the data analysis
respectively.

###  3.1 Application 1: Phase I Trial (Exposure is an Increasing Function in
Dose)

This study enrolls 20 volunteers in the dose escalation part where the drug is
given intravenously to the participants. The full data along with the DLT
information and Cmax summaries are provided in Table 1. Cmax is computed by
the clinical pharmacology scientists as a summary measure of the drug
interactions with the body and is defined by the maximum serum concentration
that a drug achieves in a specified compartment or test area of the body after
the drug has been administered and before the administration of a second dose.
Figure 2 provides the Camx measurements of the participants given the drug
doses received by them. It can be noted that Cmax (exposure) is an
approximately increasing function in the doses received by the subjects and
the two DLTs (in red) have occurred in the region that belongs to higher
values of doses and higher values of Cmax.

Table 1: Dose DLT Data of Application 1 – The columns describe the dose given to the participants, DLT: 1 or 0 according to DLT is experienced or not, and Cmax summary. Participant ID | Dose (unit) | DLT | Cmax | Participant ID | Dose (unit) | DLT | Cmax
---|---|---|---|---|---|---|---
1 | 0.1 | 0 | 2.82 | 11 | 10.0 | 0 | 297.00
2 | 0.1 | 0 | 1.68 | 12 | 10.0 | 0 | 164.00
3 | 0.3 | 0 | 2.87 | 13 | 30.0 | 0 | 401.00
4 | 0.3 | 0 | 8.14 | 14 | 30.0 | 0 | 766.00
5 | 0.3 | 0 | 4.39 | 15 | 30.0 | 0 | 592.00
6 | 1.0 | 0 | 13.30 | 16 | 30.0 | 1 | 769.00
7 | 1.0 | 0 | 12.20 | 17 | 30.0 | 0 | 577.00
8 | 3.0 | 0 | 79.70 | 18 | 50.0 | 0 | 834.00
9 | 3.0 | 0 | 47.20 | 19 | 50.0 | 0 | 1160.00
10 | 3.0 | 0 | 41.00 | 20 | 50.0 | 1 | 1160.00
Figure 2: Cmax of the the participants in Example 1 have been plotted in the
log scale. Participants who have experienced DLT (in red) have higher Cmax
(exposure) values than those who have not experienced DLT (in green) for
different administered doses.

We fit both BLRM and BLRM-PK models to the data. To infer the safe dose
recommendation we compute three probabilities – probability of underdosing,
probability of target toxicity, and probability of overdosing – as discussed
in Section 2.3. To carry out the safe dose recommendation in optimize manner
we follow the EWOC principle with $\alpha=0.25$, that is, a dose will be said
to be safe to receive if the overdose probability is less than 0.25. The
resulting probabilities for each doses is plotted in Figure 3 after fitting
BLRM (left hand side plot) and BLRM-PK model (right hand side plot). From the
figures it is evident that dose 50 units can not be recommended as the next
safe dose (in red) as the corresponding overdosing probability is beyond 0.25.
This is confirmed by both the methods – BLRM and the proposed BLRM-PK method.
Intuitively, it means that higher the dose levels, the drug triggers higher
Cmax in the body and thereby it increases the risk of toxicity observed by the
volunteers.

Figure 3: Posterior DLT rates after fitting BLRM (left hand side) and BLRM-PK
(right hand side) in Example 1. The posterior probabilities are probability of
underdosing (bottom panel), probability of target toxicity (TT), and
probability of overdosing (OD). Dose that has more than 0.25 probability of
overdosing is shown in red (EWOC criteria).

###  3.2 Application 2: Phase I Trial (Similar Exposure for Different
Toxicity)

The second application is a phase I clinical trial. Even though the original
study had different parts depending on different regimanes, we concentrate on
the Phase I dose finding study which consists of the experiment where 39 study
participants receives the investigational drug once in three weeks
intravenously. The full dose information and dose limiting toxicity
information is provided in Table 2. For the analysis we consider Cmax, a
standard summary measurement in pharmaceutics, as in the previous example.

Table 2: Dose DLT Data of Application 2 – The columns describe the dose given to the participants, DLT: 1 or 0 according to DLT is experienced or not, and Cmax summary. Participant ID | Dose (unit) | DLT | Cmax | Participant ID | Dose (unit) | DLT | Cmax
---|---|---|---|---|---|---|---
1 | 0.13 | 0 | 5170 | 21 | 3.20 | 0 | 120000
2 | 0.13 | 0 | 4220 | 22 | 2.80 | 0 | 57800
3 | 0.33 | 0 | 8500 | 23 | 2.80 | 0 | 118000
4 | 0.33 | 0 | 14500 | 24 | 2.80 | 1 | 73900
5 | 0.83 | 0 | 29600 | 25 | 2.80 | 0 | 113000
6 | 1.40 | 0 | 71600 | 26 | 2.80 | 0 | 96000
7 | 1.40 | 0 | 60200 | 27 | 2.80 | 0 | 85200
8 | 1.87 | 0 | 78400 | 28 | 3.20 | 0 | 104000
9 | 2.47 | 0 | 102000 | 29 | 3.20 | 0 | 89600
10 | 2.47 | 0 | 115000 | 30 | 2.10 | 0 | 109000
11 | 0.83 | 1 | 20100 | 31 | 2.80 | 0 | 92100
12 | 1.87 | 0 | 50700 | 32 | 2.10 | 0 | 84100
13 | 2.47 | 0 | 110000 | 33 | 2.10 | 0 | 102000
14 | 1.40 | 1 | 49400 | 34 | 2.80 | 0 | 179000
15 | 1.87 | 0 | 31200 | 35 | 2.80 | 1 | 98200
16 | 2.47 | 0 | 40600 | 36 | 2.80 | 0 | 129000
17 | 1.40 | 0 | 87900 | 37 | 3.20 | 1 | 93300
18 | 1.87 | 0 | 81400 | 38 | 3.20 | 1 | 123000
19 | 2.47 | 0 | 101000 | 39 | 3.20 | 0 | 90300
20 | 2.47 | 0 | 149000 |  |  |  |
Figure 4: Cmax of the the participants in Example 2 have been plotted in the
log scale. Participants who have experienced DLT (in red) and those who have
not experienced DLT (in green) have similar Cmax (exposure) profile for
different administered doses.

In Figure 4 we depict the observed Cmax of the participants across different
doses administered on the subjects. It is interesting to note that the
participants who have experienced DLT (in red) and those who have not
experienced DLT (in green) have similar Cmax (exposure) profile for different
administered doses. Hence, it would be worth to carry out the BLRM-PK analysis
if the phramacokinatic data has any impact on the inference of the dose
escalation decision when incorporating the pharmacokinetic data.

Toward this end, first the usual Bayesian logistic regression model (BLRM) is
fitted to the data. The resulting posterior DLT rates can be summarized by
left hand side plot in Figure 5 classified into probability of underdosing,
probability of target toxicity, and probability of overdosing. According to
the EWOC criteria described in Section 2 a dose is said to be safe to
administer at level $\alpha$ if probability of overdosing is less than
$\alpha$. As in the previous example $\alpha$ is set to 0.25. We note that the
dose 3.2 units is not a safe dose according to this criteria. On the other
hand, when the proposed BLRM-PK model is fitted in this data, the same EWOC
criteria implies that dose 3.2 mg is a safe dose (Right panel of Figure 5).
Since the actual true scenario is not known in the real data, further study,
for instance, a rigorous simulation study is needed to better understand the
utility of the BLRM-PK method which is conducted in the next Section.

Figure 5: Posterior DLT rates after fitting BLRM (left hand side) and BLRM-PK
(right hand side) in Example 2. The posterior probabilities are probability of
underdosing (bottom panel), probability of target toxicity (TT), and
probability of overdosing (OD). Dose that has more than 0.25 probability of
overdosing is shown in red (EWOC criteria).

##  4 Simulation Study

In this section we conduct a simulation study to illustrate the effectiveness
of incorporating the pharmacokinetic information in the formal Statistical
modeling via the proposed BLRM-PK model. We start with known exposure profile,
see Table 3 for each dose. In essence, at each dose we assume that it
generates a given exposure according to Table 3. To be consistent with the
model formulation in Section 2.1 we assume the exposure data is in log scale
in the sense that the actual exposure can be obtained by raising them in the
power of exponential function. By setting
$(\log(\alpha),\log(\beta))=(-2.5,2)$ in (4) the data generating toxicity
profile is obtained as in Table 3.

Table 3: Exposure and Toxicity Profile for Data Generation in Simulation. Dose (unit) | Exposure | Toxicity | Dose (unit) | Exposure | Toxicity
---|---|---|---|---|---
First Scenario | Second Scenario
0.1 | 0.40 | 0.15 | 0.13 | 0.05 | 0.125
0.3 | 0.47 | 0.17 | 0.33 | 0.13 | 0.133
1.0 | 0.53 | 0.19 | 0.83 | 0.28 | 0.151
3.0 | 0.60 | 0.21 | 1.40 | 0.67 | 0.210
10.0 | 0.67 | 0.24 | 1.87 | 0.75 | 0.223
30.0 | 0.73 | 0.26 | 2.10 | 1.10 | 0.289
50.0 | 0.80 | 0.29 | 2.47 | 1.15 | 0.299
|  |  | 2.80 | 1.15 | 0.300
|  |  | 3.2 | 1.16 | 0.302

The data generation and model fitting scheme proceeds in the following way. At
dose level 0.1 we generate three DLT data points from a Binomial distribution
with size 3 and probability 0.15. Moreover, we generate the pharmacokinetic
summary from a log normal distribution with mean 0.40 and standard deviation
0.5. Then a BLRM is fitted using the dose DLT data and proposed BLRM-PK method
is fitted using dose-DLT-pharmacokinatic data. Then the next safe dose is
selected according to the EWOC criteria. The next dose recommendation is the
maximum dose which satisfies the EWOC criteria but does not exceed maximum
allowable increments. In this example the maximal allowable increment is set
to 200%, that is, the experiment can jump from the dose level 0.1 units to 0.3
units if EWOC criteria is satisfied. When the next dose is suggested by BLRM
or BLRM-PK is selected, then again 3 data points are generated in a very
similar manner as in first dose level. In this way the trial is carried out
until the stopping criteria is met:

  1. 1.

No further escalation is possible.

  2. 2.

At least 6 patients have been treated at the recommended MTD $\overline{d}$.

  3. 3.

The dose $\overline{d}$ satisfies one of the following conditions:

     * •

The probability of target toxicity at dose $\overline{d}$ exceeds 0.5, i.e.
$\Pr(0.16\leq\pi_{\overline{d}}<0.33)\geq 0.5$;

     * •

A minimum of 15 patients have been treated in the trial.

The maximum number of patients in dose escalation part of the trial is set to
50. We generate 1000 trials in this fashion and assess the performance of the
proposed BLRM-PK method with that of BLRM. When assessing the performance of
the methods, following matrices are populated and compared:

  * •

Percentage of participants receiving dose in the target toxicity interval

  * •

Percentage of participants receiving an overdose

  * •

Percentage of participants receiving an under dose

  * •

Probability that recommended MTD at the end of the trial is in the target
toxicity interval

  * •

Probability that recommended MTD is an overdose

  * •

Probability that recommended MTD is an under dose

  * •

Percentage of trials with MTD below starting dose

  * •

Average sample size.

Table 4: Simulation Result. TT: Target Toxicity, OD: Overdose, UD: Underdose, SD: Starting Dose. Method | Participant allocation |  $\Pr$(MTD |  $\Pr$(MTD |  $\Pr$(MTD |  $\Pr$(MTD | Average
---|---|---|---|---|---|---
| in TT (%) | in OD (%) | in UD (%) | in TT) | in OD) | in UD) | is below | sample
|  |  |  |  |  |  | SD) | size
Simulation Example 1
BRLM | 51.1 | 0.0 | 48.9 | 0.69 | 0.0 | 0.02 | 0.29 | 10.0
BLRM-PK | 71.7 | 0.0 | 28.3 | 0.90 | 0.0 | 0.03 | 0.07 | 14.0
Simulation Example 2
BRLM | 51.1 | 0.0 | 48.9 | 0.69 | 0.0 | 0.02 | 0.29 | 10.0
BLRM-PK | 69.2 | 0.0 | 30.8 | 0.86 | 0.0 | 0.05 | 0.09 | 13.8
Simulation Example 3
BRLM | 25.3 | 0.0 | 74.7 | 0.41 | 0.0 | 0.28 | 0.31 | 15.9
BLRM-PK | 48.9 | 0.0 | 51.1 | 0.80 | 0.0 | 0.16 | 0.04 | 22.2
Simulation Example 4
BRLM | 25.3 | 0.0 | 74.7 | 0.41 | 0.0 | 0.28 | 0.31 | 15.9
BLRM-PK | 45.3 | 0.0 | 54.6 | 0.74 | 0.0 | 0.21 | 0.05 | 21.5

We carry out four simulation examples and present the frequentist operating
characteristics in Table 4.

  * •

Simulation Example 1: PK measurements are generated from scenario 1 in Table 3
with standard deviation 0.5.

  * •

Simulation Example 2: PK measurements are generated from scenario 1 in Table 3
with standard deviation 1.

  * •

Simulation Example 3: PK measurements are generated from scenario 2 in Table 3
with standard deviation 0.5.

  * •

Simulation Example 4: PK measurements are generated from scenario 2 in Table 3
with standard deviation 1.

From the numbers presented in Table 4 it is evident that BLRM-PK method
recovers the MTD which falls in the target toxicity region more than that of
BLRM method. For instance, in case of Simulation Example 1, the estimated
probability from the 1000 trials that the the dose identified by the MTD by
BLRM-PK model is a dose from the target toxicity region is 0.90 while the same
due to the BLRM method is 0.69. Furthermore, when BLRM-PK method is fitted
then 71.1% of the total participants receive a dose belonging to the target
toxicity region while the BLRM method allocates 51.1% of the total
participants in the target toxicity region. Since it is the intend to maximize
number of participants in the target toxicity region in an Oncology Phase I
trial, obviuosly BLRM-PK method outperfoms the traditional BLRM method in this
regard. We note that, the conclusion remains same across all the Simulation
Examples.

##  5 Conclusion

In this analysis we have shown that the joint modeling of dose-exposure and
exposure-DLT approach has the potential to provide better understanding about
the dose escalation curve. Furthermore it is more consistent design in
estimating the correct MTD which is evident from the simulation. One
assumption in the proposed model is that the dose exposure relationship is
monotone which is a natural assumption according to the clinical pharmacology
scientists (personal communication). Moreover, it allows to satisfy some
sufficient conditions to consistently estimate the MTD. Additionally, by
virtue of these assumption, it can be theoretically shown that the proposed
design is coherent which is desirable property in dose escalation trials. In
case, this assumption is tenable, then this can be relaxed and the model can
be easily written to accommodate that. In our experience, we have seen that,
in the long run, such as in simulation study this does not affect the
conclusion about the superiority of the BLRM-PK approach.

Our research is thus an important contribution to the overreaching goal of the
ongoing efforts to model other important clinical measures and thereby
informing the researchers quantitatively when a dose finding study is run.
Nonetheless, it opens up the other avenues for other future research such as
eliciting an appropriate prior for the dose exposure model parameters when and
if any knowledge is available about the dose exposure relationship.

## References

  * Babb et al., (1998) Babb, J., Rogatko, A., and Zacks, S. (1998).  Cancer phase I clinical trials: efficient dose escalation with overdose control.  Statistics in Medicine, 17(10):1103–1120.
  * Cheung, (2005) Cheung, Y. K. (2005).  Coherence principles in dose-finding studies.  Biometrika, 92(4):863–873.
  * Collins et al., (1990) Collins, J. M., Grieshaber, C. K., and Chabner, B. A. (1990).  Pharmacologically guided phase I clinical trials based upon preclinical drug development.  JNCI: Journal of the National Cancer Institute, 82(16):1321–1326.
  * Collins et al., (1986) Collins, J. M., Zaharko, D. S., Dedrick, R. L., and Chabner, B. A. (1986).  Potential roles for preclinical pharmacology in phase I clinical trials.  Cancer Treat Rep, 70(1):73–80.
  * Corporation, (1992) Corporation, C. S. (1992).  EaSt: A software package for the design and interim monitoring of group sequential clinical trials.
  * Cotterill et al., (2015) Cotterill, A., Lorand, D., Wang, J., and Jaki, T. (2015).  A practical design for a dual-agent dose-escalation trial that incorporates pharmacokinetic data.  Statistics in Medicine, 34(13):2138–2164.
  * Crowley and Hoering, (2012) Crowley, J. and Hoering, A. (2012).  Handbook of Statistics in Clinical Oncology.  CRC Press.
  * Daimon et al., (2019) Daimon, T., Hirakawa, A., and Matsui, S. (2019).  Dose-Finding Designs for Early-Phase Cancer Clinical Trials: A Brief Guidebook to Theory and Practice.  Springer.
  * Ferlay et al., (2015) Ferlay, J., Soerjomataram, I., Dikshit, R., Eser, S., Mathers, C., Rebelo, M., Parkin, D. M., Forman, D., and Bray, F. (2015).  Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012.  International Journal of Cancer, 136(5):E359–E386.
  * Geller, (1984) Geller, N. L. (1984).  Design of phase I and II clinical trials in cancer: a statistician’s view.  Cancer Investigation, 2(6):483–491.
  * Holford, (1995) Holford, N. H. (1995).  The target concentration approach to clinical drug development.  Clinical Pharmacokinetics, 29(5):287–291.
  * Hummel et al., (2009) Hummel, J., McKendrick, S., Brindley, C., and French, R. (2009).  Exploratory assessment of dose proportionality: review of current approaches and proposal for a practical criterion.  Pharmaceutical Statistics: The Journal of Applied Statistics in the Pharmaceutical Industry, 8(1):38–49.
  * Ji et al., (2010) Ji, Y., Liu, P., Li, Y., and Nebiyou Bekele, B. (2010).  A modified toxicity probability interval method for dose-finding trials.  Clinical Trials, 7(6):653–663.
  * Neuenschwander et al., (2008) Neuenschwander, B., Branson, M., and Gsponer, T. (2008).  Critical aspects of the Bayesian approach to phase I cancer trials.  Statistics in Medicine, 27(13):2420–2439.
  * Neuenschwander et al., (2015) Neuenschwander, B., Matano, A., Tang, Z., Roychoudhury, S., Wandel, S., and Bailey, S. (2015).  A Bayesian industry approach to phase I combination trials in oncology.  In Zhao, W. and Yang, H., editors, Statistical Methods in Drug Combination Studies, pages 95–135. Chapman & Hall/CRC Press: Boca Raton, FL.
  * Newell, (1994) Newell, D. R. (1994).  Pharmacologically based phase I trials in cancer chemotherapy.  Hematology/Oncology Clinics of North America, 8(2):257–275.
  * O’Quigley et al., (1990) O’Quigley, J., Pepe, M., and Fisher, L. (1990).  Continual reassessment method: a practical design for phase 1 clinical trials in cancer.  Biometrics, pages 33–48.
  * O’Quigley and Shen, (1996) O’Quigley, J. and Shen, L. Z. (1996).  Continual reassessment method: a likelihood approach.  Biometrics, pages 673–684.
  * Piantadosi and Liu, (1996) Piantadosi, S. and Liu, G. (1996).  Improved designs for dose escalation studies using pharmacokinetic measurements.  Statistics in Medicine, 15(15):1605–1618.
  * Plummer, (2003) Plummer, M. (2003).  JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling.  In Proceedings of the 3rd International Workshop on Distributed Statistical Computing, volume 124, pages 1–10. Vienna, Austria.
  * Plummer, (2019) Plummer, M. (2019).  rjags: Bayesian Graphical Models using MCMC.  R package version 4-10.
  * R Core Team, (2019) R Core Team (2019).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rosenberger and Haines, (2002) Rosenberger, W. F. and Haines, L. M. (2002).  Competing designs for phase I clinical trials: a review.  Statistics in Medicine, 21(18):2757–2770.
  * Shaw et al., (2014) Shaw, A. T., Kim, D.-W., Mehra, R., Tan, D. S., Felip, E., Chow, L. Q., Camidge, D. R., Vansteenkiste, J., Sharma, S., De Pas, T., et al. (2014).  Ceritinib in ALK-rearranged non–small-cell lung cancer.  New England Journal of Medicine, 370:1189–1197.
  * Takeda et al., (2018) Takeda, K., Komatsu, K., and Morita, S. (2018).  Bayesian dose-finding phase I trial design incorporating pharmacokinetic assessment in the field of oncology.  Pharmaceutical Statistics, 17(6):725–733.
  * Tighiouart and Rogatko, (2010) Tighiouart, M. and Rogatko, A. (2010).  Dose finding with escalation with overdose control (EWOC) in cancer clinical trials.  Statistical Science, pages 217–226.
  * Toumazi et al., (2018) Toumazi, A., Comets, E., Alberti, C., Friede, T., Lentz, F., Stallard, N., Zohar, S., and Ursino, M. (2018).  dfpk: an R-package for Bayesian dose-finding designs using pharmacokinetics (PK) for phase I clinical trials.  Computer Methods and Programs in Biomedicine, 157:163–177.
  * Ursino et al., (2017) Ursino, M., Zohar, S., Lentz, F., Alberti, C., Friede, T., Stallard, N., and Comets, E. (2017).  Dose-finding methods for phase I clinical trials using pharmacokinetics in small populations.  Biometrical Journal, 59(4):804–825.
  * Whitehead et al., (2001) Whitehead, J., Patterson, S., Webber, D., Francis, S., and Zhou, Y. (2001).  Easy-to-implement Bayesian methods for dose-escalation studies in healthy volunteers.  Biostatistics, 2(1):47–61.
  * Whitehead et al., (2007) Whitehead, J., Zhou, Y., Hampson, L., Ledent, E., and Pereira, A. (2007).  A Bayesian approach for dose-escalation in a phase I clinical trial incorporating pharmacodynamic endpoints.  Journal of Biopharmaceutical Statistics, 17(6):1117–1129.
  * Zacks et al., (1998) Zacks, S., Rogatko, A., and Babb, J. (1998).  Optimal Bayesian-feasible dose escalation for cancer phase I trials.  Statistics and Probability Letters, 38(3):215–220.
  * Zhou et al., (2018) Zhou, H., Murray, T. A., Pan, H., and Yuan, Y. (2018).  Comparative eview of novel model-assisted designs for phase I clinical trials.  Statistics in Medicine, 37(14):2208–2222.


# PRESERVE PRE-TRAINED KNOWLEDGE: TRANSFER LEARNING WITH SELF-DISTILLATION FOR
ACTION RECOGNITION

###### Abstract

Video-based action recognition is one of the most popular topics in computer
vision. With recent advances of self-supervised video representation learning
approaches, action recognition usually follows a two-stage training framework,
i.e., self-supervised pre-training on large-scale unlabeled sets and transfer
learning on a downstream labeled set. However, catastrophic forgetting of the
pre-trained knowledge becomes the main issue in the downstream transfer
learning of action recognition, resulting in a sub-optimal solution. In this
paper, to alleviate the above issue, we propose a novel transfer learning
approach that combines self-distillation in fine-tuning to preserve knowledge
from the pre-trained model learned from the large-scale dataset. Specifically,
we fix the encoder from the last epoch as the teacher model to guide the
training of the encoder from the current epoch in the transfer learning. With
such a simple yet effective learning strategy, we outperform state-of-the-art
methods on widely used UCF101 and HMDB51 datasets in action recognition task.

11footnotetext: Corresponding author.

Index Terms—  Video action recognition, self-distillation, transfer learning

##  1 Introduction

Action recognition is a fundamental topic of computer vision. It has gained
significant attention due to its wide range of applications in video
surveillance, video recommendation, and human-computer interaction [1]. Early
works [2, 3] rely on a large amount of labeled data for learning spatial-
temporal representations, requiring huge annotation costs. To alleviate the
dependence on the labeled data, pre-training strategies for self-supervised
video representation learning have attracted extensive attention. These
approaches, such as IIC[4] and CVRL[5], require only a large-scale set with
few or no annotations.

After self-supervised pre-training, usually, transfer learning is employed for
downstream tasks on the labeled set. However, catastrophic forgetting of the
pre-trained knowledge would become the main issue in the transfer learning
stage. When adapting the model to the labeled set, the knowledge of large-
scale unlabeled data will be gradually forgotten.

To address the above issue, we propose a novel transfer learning method that
combines self-distillation in the fine-tuning stage. Specifically, at the
beginning of each training epoch, we make a copy of the base encoder network
from the last epoch as the teach model and freeze the model parameters. Except
for the conventional cross-entropy loss used for classification in training,
we constrain the difference between representations from the teacher model and
the base encoder with the Euclidean distance as the regularization term. As a
result, the previous knowledge of the model can be largely preserved. With
such a simple yet effective transfer learning approach, we achieve 82.3$\%$
and 47.3$\%$ top-1 accuracy on UCF101[6] and HMDB51[7] datasets, respectively,
outperforming state-of-the-art methods.

In summary, the main contributions of this paper are as follows:

  * •

We propose a novel self-distillation guided transfer knowledge approach for
action recognition that can alleviate the catastrophic forgetting phenomenon
and preserve the pre-trained knowledge learned from large-scale datasets.

  * •

Our method can be fitted for different architectures and models, and easily
applied to any existing transfer learning frameworks.

  * •

We achieve state-of-the-art performance on commonly used action recognition
datasets to demonstrate the effectiveness of our proposed method.

Fig. 1: Overview of our proposed self-distillation guided transfer learning
framework. At the beginning of each epoch, we make a copy of the spatial-
temporal encoder from the last epoch as the teacher network
$f^{\prime}(\cdot)$. The Euclidean distance between representations of both
encoders is treated as a regularization in transfer learning. Besides that,
the conventional cross-entropy loss is employed on the head network $g(\cdot)$
for action recognition.

##  2 Related Work

###  2.1 Encoder Models for Action Recognition

Encoder models play an essential role in representation learning in video-
based action recognition. To learn spatial-temporal information of videos,
different types of encoding architecture are widely adopted, including two-
stream networks, CNN-LSTM networks, 3D CNN, and vision transformers. For
example, the two-stream architecture, which was first proposed in [2],
incorporates spatial and temporal streams based on CNN. [8] presents LRCN,
which combines CNN and LSTM to learn the spatial features and temporal
features more thoroughly. [3] proposes 3D CNN to learn spatial-temporal
features directly. It performs better compared with 2D CNN. Moreover,
transformer-based encoder has been applied to action recognition and achieves
brilliant results. For example, [9] handles the large number of spatial-
temporal tokens extracted from multiple modalities and promotes learning the
inter-modal interactions. Additionally, in recent years, light-weighted models
with less computational cost have been explored in action recognition.
Specifically, [10] replaces computational expensive 3D convolutions with
lower-cost 2D convolutions. [11] proposes SlowFast model, which applies a Slow
pathway to capture spatial semantics and applies a Fast pathway to capture
temporal information, like a two-stream model. [12] suggests factorizing 3D
convolutions by separating channel interactions and spatial-temporal
interactions. [13] expands a tiny 2D image classification architecture. The
current study shows that it is still a challenge to learn good representations
with long periods of video actions. Without guidance on action labels, the
encoders of action recognition learn background features rather than
foreground actions. They analyze the background to assist action recognition,
which deviates from the purpose of the representation learning task.

###  2.2 Self-Supervised Learning for Action Recognition Pre-Training

With a large-scale unlabeled dataset, self-supervised video representation
learning has gained more attention from the community. For example, [5]
presents a self-supervised video representation learning method with
contrastive loss, where two augmented clips from the same short video are
pulled together, and clips from different videos are pushed away. [14]
presents a new way to detect the playback speed and exploits the relative
speed between two video clips as labels. [14] also proposes to enforce the
model to capture the appearance difference between two video clips. [15]
decomposes the representation learning objective into two sub-tasks,
respectively, emphasizing spatial and temporal features. [16] proposes to
track cycle-consistency as free supervisory signal for learning visual
representations. [17] proposes a new architecture and learning framework named
Memory-augmented Dense Predictive Coding for self-supervised video
representation learning task. It trains the network with a predictive
attention mechanism over the set of compressed memories and represents the
future states with the combination of condensed representations in memory.
Current pre-training-based action recognition methods usually adopt simple
transfer learning on the downstream labeled dataset. Designing more effective
transfer learning approaches to address the catastrophic forgetting phenomenon
of pre-trained knowledge is still under-investigated.

###  2.3 Self-Distillation

Knowledge distillation was firstly proposed in [18] to transfer the knowledge
from the cumbersome model to a small model. Self-distillation can be treated
as a regularization strategy when the teacher model and the student model
share the same architecture. For example, [19] first uses self-distillation,
which transfers the knowledge in the deeper portion of the networks to the
shallow sections. [20] improves representation quality of the smaller models
by single-stage online knowledge distillation. [21] uses self-distillation
between different epochs with soft targets. However, combining self-
distillation as the regularization for transfer learning is not well-
exploited.

##  3 Method

In this section, we demonstrate our proposed self-distillation guided transfer
learning framework for action recognition in detail. We first show the
conventional transfer learning framework for downstream tasks with pre-trained
encoder network. Then, we illustrate how to incorporate the self-distillation
to adapt the network on the downstream dataset and memorize the pre-trained
knowledge simultaneously.

###  3.1 Transfer Learning with Pre-Trained Model

Denote $f_{\boldsymbol{\theta}}(\cdot)$ as the base encoder network with
parameters $\boldsymbol{\theta}$ that are pre-trained on the large-scale
unlabeled dataset via self-supervised learning [14]. To achieve the goal of
action recognition as the downstream task, we stack the base encoder with an
extra classification head network $g_{\boldsymbol{\phi}}(\cdot)$, where
$\boldsymbol{\phi}$ represents the parameters of the head network. Then we
apply transfer learning to the entire network with the following cross-entropy
loss,

$\mathcal{L}_{CE}(\mathcal{X}|\boldsymbol{\phi},\boldsymbol{\theta})=-\sum_{i=1}^{N}\log(g^{y_{i}}_{\boldsymbol{\phi}}(f_{\boldsymbol{\theta}}(\boldsymbol{x}_{i}))),$
(1)

where $\mathcal{X}$ represents the training data, $(\boldsymbol{x}_{i},y_{i})$
is the individual data-label pair,
$g^{y_{i}}_{\boldsymbol{\phi}}(f_{\boldsymbol{\theta}}(\boldsymbol{x}_{i}))$
represents the prediction of the ground truth class $y_{i}$ after softmax
operation, and $N$ is the number of training samples. We can obtain the
learned parameters $\hat{\boldsymbol{\phi}}$ and $\hat{\boldsymbol{\theta}}$
by minimizing the above loss.

###  3.2 Self-Distillation Guided Transfer Learning

Eq. (1) is the conventional transfer learning for the downstream task.
However, it has the issue of catastrophic forgetting of the pre-trained
knowledge. To address the above issue, we propose a self-distillation guided
transfer learning, aiming to adapt the network on the downstream dataset and
memorize the pre-trained knowledge simultaneously. Specifically, at the
beginning of each epoch $t$, we make a copy of the encoder network from the
last epoch $t-1$ and treat the encoder
$f^{\prime}_{\boldsymbol{\theta}^{t-1}}(\cdot)$ as the teacher model. We fix
the parameters of the teacher model at the current epoch in training. Then the
self-distillation is defined between the encoder model and the teacher model
as follows,

$\mathcal{L}_{SD}(\mathcal{X}|\boldsymbol{\theta}^{t})=\sum_{i=1}^{N}||f_{\boldsymbol{\theta}^{t}}(\boldsymbol{x}_{i})-f^{\prime}_{\boldsymbol{\theta}^{t-1}}(\boldsymbol{x}_{i})||^{2}.$
(2)

The base encoder model plays a role as the student model in the self-
distillation. The self-distillation is treated as a regularization term in
transfer learning and used for memorizing the knowledge from the last epoch.
At the first epoch, i.e., $t=1$, we use the pre-trained base encoder as the
teacher model $f^{\prime}_{\boldsymbol{\theta}^{0}}(\cdot)$.

Combined with cross-entropy loss in Eq. (1) and self-distillation loss in Eq.
(2), the total loss $L$ for the $t$-th epoch can be written as follows,

$\mathcal{L}(\mathcal{X}|\boldsymbol{\phi}^{t},\boldsymbol{\theta}^{t})=\mathcal{L}_{CE}(\mathcal{X}|\boldsymbol{\phi}^{t},\boldsymbol{\theta}^{t})+\lambda\mathcal{L}_{SD}(\mathcal{X}|\boldsymbol{\theta}^{t}),$
(3)

where $\lambda$ is the weight for self-distillation loss. The training details
are summarized in Algorithm 1.

1Input: Labeled set $\mathcal{X}$ and pre-trained model
$f_{\boldsymbol{\theta}}$.

2

3Initialization: Initialize head classifier $g_{\boldsymbol{\phi}}$.

4

5Set $\boldsymbol{\theta}^{0}\leftarrow\boldsymbol{\theta}$.

6 for  _$t=1,2,...,t_{max}$_ do

7 $f^{\prime}_{\boldsymbol{\theta}^{t-1}}\leftarrow
f_{\boldsymbol{\theta}^{t-1}}$ and freeze
$f^{\prime}_{\boldsymbol{\theta}^{t-1}}$.

8 for  _Batch data $\mathcal{B}$ in $\mathcal{X}$_ do

9
$\boldsymbol{\theta}^{t}\leftarrow\boldsymbol{\theta}^{t-1}-\eta\boldsymbol{\bigtriangledown}_{\boldsymbol{\theta}^{t-1}}\mathcal{L}(\mathcal{B}|\boldsymbol{\phi}^{t},\boldsymbol{\theta}^{t})$.

10

11
$\boldsymbol{\phi}^{t}\leftarrow\boldsymbol{\phi}^{t-1}-\eta\boldsymbol{\bigtriangledown}_{\boldsymbol{\phi}^{t-1}}\mathcal{L}(\mathcal{B}|\boldsymbol{\phi}^{t},\boldsymbol{\theta}^{t})$.

12

13

14

15Set $\hat{\boldsymbol{\theta}}\leftarrow{\boldsymbol{\theta}^{t_{max}}}$,
$\hat{\boldsymbol{\phi}}\leftarrow{\boldsymbol{\phi}^{t_{max}}}$.

16

17Output: $\hat{\boldsymbol{\theta}},\hat{\boldsymbol{\phi}}$.

18

Algorithm 1 Training

##  4 Experiments and Results

###  4.1 Datasets

We adopt top-1 and top-5 accuracy (denoted as Acc@1 and Acc@5, respectively)
as evaluation metrics. To verify the effectiveness of our proposed method, we
use UCF101[6] and HMDB51[7] datasets in the evaluation. The statistics of the
two datasets are summarized as follows.

UCF101 [6] is a widely used dataset which contains 13,320 clips with 101
action categories. Following previous work [22], we use the training split-1
as the finetuning dataset and testing split-1 for evaluation.

HMDB51 [7] is a dataset which consists of 6,849 clips with 51 human action
classes from YouTube. It is divided into three training/testing splits.
Following prior work [22], we use training/testing split-1 in the experiment.

###  4.2 Implementation Details

We adopt the pre-trained model based on RSPNet [14]. The feature generated by
spatial-temporal encoder $f_{\boldsymbol{\theta}}(\cdot)$ is of size $24\times
512$ and $32\times 512$ for C3D[3] and R(2+1)D[23], respectively.
Classification head $g_{\boldsymbol{\phi}}(\cdot)$ is composed of one fully
connected layer. In the transfer learning stage, we use stochastic gradient
descent (SGD) optimizer with learning rate 0.005 for C3D model and 0.1 for
R(2+1)D model. All the models are fine-tuned for 100 epochs. The setting of
the hyper-parameter $\lambda$ is illustrated in the ablation study.

Table 1: Performance on UCF101 and HMDB51 compared with RSP. * symbol indicates that we fine-tune 100 epochs on original RSP basis. | UCF101 | HMDB51
---|---|---
Method |  Acc@1 |  Acc@5 |  Acc@1 |  Acc@5
C3D Architecture
3D ST [24] |  60.6 |  - |  28.3 |  -
MAS [25] |  61.2 |  - |  33.4 |  -
RSPNet [14] |  76.7 |  - |  44.6 |  -
RSPNet*[14] |  77.1 |  93.5 |  43.9 |  76.5
Ours |  78.6 |  94.6 |  45.8 |  76.5
R(2+1)D Architecture
Pace [26] |  77.1 |  - |  36.6 |  -
RSPNet [14] |  81.1 |  - |  44.6 |  -
RSPNet* [14] |  80.2 |  95.0 |  44.1 |  75.6
Ours |  82.3 |  96.3 |  47.3 |  77.8

###  4.3 Comparison with State-of-the-Art Methods

The comparison with other state-of-the-art methods is shown in Table 1,
including Acc@1 and Acc@5 results with C3D [3] and R(2+1)D [23] architectures
on UCF101 and HMDB51 datasets, respectively. The compared methods include 3D
ST [24], MAS [25], Pace [26] and RSPNet [14]. Specifically, 3D ST [24] solves
3D video cubic puzzles by using pretext task in self-supervised video
represent learning. Instead of designing new self-supervised task, MAS [25]
regresses both motion and appearance statistics along spatial and temporal
dimensions to learn visual features. Pace[26] novelly uses pace prediction to
learn spatio-temporal features, and RSPNet [14] further proposes relative
speed perception as a more effective supervision method for video presentation
learning. For a fair comparison, we re-implement RSPNet with 100-epoch fine-
tuning, denoted as RSPNet*, where RSPNet represents originally reported
results in [14] that are based on 30-epoch fine-tuning. From the results, our
method consistently improves the Acc@1 and Acc@5 performance on two datasets,
verifying the effectiveness of our method. Specifically, when using R(2+1)
model, we achieve 82.3$\%$ Acc@1 on UCF101 dataset and 47.3$\%$ on HMDB51
dataset, which is significantly higher than other state-of-the-art methods.

###  4.4 Ablation Study

To analyze the influence of the self-distillation loss, we conduct ablation
study with different settings of the loss weight $\lambda$. We vary $\lambda$
from 0.1 to 1000 and report Acc@1 and Acc@5 results in Table 2 with C3D and
R(2+1)D architectures on UCF101 and HMDB51 datasets, respectively. From the
results, even with an extensive range of choices of $\lambda$, our proposed
method outperforms the baseline RSPNet in most settings, verifying the
effectiveness of the proposed self-distillation guided transfer learning
framework.

Table 2: Top-1 and Top-5 accuracy on UCF101 and HMDB51 datasets. $\lambda$ is varying from 0.1 to 1000. | UCF101 | HMDB51
---|---|---
Method |  Acc@1 |  Acc@5 |  Acc@1 |  Acc@5
C3D Architecture
RSPNet baseline |  77.1 |  93.5 |  43.9 |  76.5
Ours ($\lambda$ = 0.1) |  76.7 |  93.7 |  45.8 |  76.5
Ours ($\lambda$ = 1) |  77.1 |  94.1 |  43.9 |  76.7
Ours ($\lambda$ = 10) |  77.9 |  94.1 |  44.2 |  76.7
Ours ($\lambda$ = 100) |  78.6 |  94.6 |  42.9 |  77.4
Ours ($\lambda$ = 1000) |  77.5 |  94.7 |  41.4 |  75.9
R(2+1)D Architecture
RSPNet baseline |  80.2 |  95.0 |  44.1 |  75.6
Ours ($\lambda$ = 0.1) |  80.8 |  94.4 |  44.8 |  74.3
Ours ($\lambda$ = 1) |  80.3 |  94.8 |  43.3 |  75.1
Ours ($\lambda$ = 10) |  81.4 |  95.4 |  45.6 |  76.1
Ours ($\lambda$ = 100) |  82.3 |  96.3 |  46.3 |  76.3
Ours ($\lambda$ = 1000) |  81.5 |  96.1 |  47.3 |  77.8

##  5 Conclusion

In this paper, we propose a novel, simple yet effective transfer learning
method with self-distillation to preserve pre-training knowledge for action
recognition. The self-distillation is conducted as a regularization term in
the fine-tuning to memorize the knowledge learned from the previous model. It
alleviates the phenomenon of catastrophic forgetting on the labeled set.
Experiments on UCF101 and HMDB51 datasets show that our proposed transfer
learning method consistently improves the top-1 and top-5 accuracy for action
recognition task.

## References

  * [1] De-An Huang, Vignesh Ramanathan, Dhruv Mahajan, Lorenzo Torresani, Manohar Paluri, Li Fei-Fei, and Juan Carlos Niebles,  “What makes a video a video: Analyzing temporal information in video understanding models and datasets,”  in CVPR, 2018, pp. 7366–7375.
  * [2] Karen Simonyan and Andrew Zisserman,  “Two-stream convolutional networks for action recognition in videos,”  NeurIPS, vol. 27, 2014.
  * [3] Du Tran, Lubomir Bourdev, Rob Fergus, Lorenzo Torresani, and Manohar Paluri,  “Learning spatiotemporal features with 3d convolutional networks,”  in ICCV, 2015, pp. 4489–4497.
  * [4] Li Tao, Xueting Wang, and Toshihiko Yamasaki,  “Self-supervised video representation learning using inter-intra contrastive framework,”  in Proceedings of the 28th ACM International Conference on Multimedia, 2020, pp. 2193–2201.
  * [5] Rui Qian, Tianjian Meng, Boqing Gong, Ming-Hsuan Yang, Huisheng Wang, Serge Belongie, and Yin Cui,  “Spatiotemporal contrastive video representation learning,”  in CVPR, 2021, pp. 6964–6974.
  * [6] Khurram Soomro, Amir Roshan Zamir, and Mubarak Shah,  “Ucf101: A dataset of 101 human actions classes from videos in the wild,”  arXiv preprint arXiv:1212.0402, 2012.
  * [7] Hildegard Kuehne, Hueihan Jhuang, Estíbaliz Garrote, Tomaso Poggio, and Thomas Serre,  “Hmdb: a large video database for human motion recognition,”  in ICCV. IEEE, 2011, pp. 2556–2563.
  * [8] Jeffrey Donahue, Lisa Anne Hendricks, Sergio Guadarrama, Marcus Rohrbach, Subhashini Venugopalan, Kate Saenko, and Trevor Darrell,  “Long-term recurrent convolutional networks for visual recognition and description,”  in CVPR, 2015, pp. 2625–2634.
  * [9] Jiawei Chen and Chiu Man Ho,  “Mm-vit: Multi-modal video transformer for compressed video action recognition,”  in WACV, 2022, pp. 1910–1921.
  * [10] Saining Xie, Chen Sun, Jonathan Huang, Zhuowen Tu, and Kevin Murphy,  “Rethinking spatiotemporal feature learning: Speed-accuracy trade-offs in video classification,”  in ECCV, 2018, pp. 305–321.
  * [11] Christoph Feichtenhofer, Haoqi Fan, Jitendra Malik, and Kaiming He,  “Slowfast networks for video recognition,”  in ICCV, 2019, pp. 6202–6211.
  * [12] Du Tran, Heng Wang, Lorenzo Torresani, and Matt Feiszli,  “Video classification with channel-separated convolutional networks,”  in ICCV, 2019, pp. 5552–5561.
  * [13] Christoph Feichtenhofer,  “X3d: Expanding architectures for efficient video recognition,”  in CVPR, 2020, pp. 203–213.
  * [14] Peihao Chen, Deng Huang, Dongliang He, Xiang Long, Runhao Zeng, Shilei Wen, Mingkui Tan, and Chuang Gan,  “Rspnet: Relative speed perception for unsupervised video representation learning,”  in AAAI, 2021, vol. 1.
  * [15] Zehua Zhang and David Crandall,  “Hierarchically decoupled spatial-temporal contrast for self-supervised video representation learning,”  in WACV, 2022, pp. 3235–3245.
  * [16] Xiaolong Wang, Allan Jabri, and Alexei A. Efros,  “Learning correspondence from the cycle-consistency of time,”  in CVPR, 2019.
  * [17] Tengda Han, Weidi Xie, and Andrew Zisserman,  “Memory-augmented dense predictive coding for video representation learning,”  in ECCV. Springer, 2020, pp. 312–329.
  * [18] Geoffrey Hinton, Oriol Vinyals, Jeff Dean, et al.,  “Distilling the knowledge in a neural network,”  arXiv preprint arXiv:1503.02531, vol. 2, no. 7, 2015.
  * [19] Linfeng Zhang, Jiebo Song, Anni Gao, Jingwei Chen, Chenglong Bao, and Kaisheng Ma,  “Be your own teacher: Improve the performance of convolutional neural networks via self distillation,”  in ICCV, 2019, pp. 3713–3722.
  * [20] Prashant Bhat, Elahe Arani, and Bahram Zonooz,  “Distill on the go: Online knowledge distillation in self-supervised learning,”  in CVPR, 2021, pp. 2678–2687.
  * [21] Kyungyul Kim, ByeongMoon Ji, Doyoung Yoon, and Sangheum Hwang,  “Self-knowledge distillation with progressive refinement of targets,”  in ICCV, 2021, pp. 6567–6576.
  * [22] Dejing Xu, Jun Xiao, Zhou Zhao, Jian Shao, Di Xie, and Yueting Zhuang,  “Self-supervised spatiotemporal learning via video clip order prediction,”  in CVPR, 2019, pp. 10334–10343.
  * [23] Du Tran, Heng Wang, Lorenzo Torresani, Jamie Ray, Yann LeCun, and Manohar Paluri,  “A closer look at spatiotemporal convolutions for action recognition,”  in CVPR, 2018, pp. 6450–6459.
  * [24] Dahun Kim, Donghyeon Cho, and In So Kweon,  “Self-supervised video representation learning with space-time cubic puzzles,”  in AAAI, 2019, vol. 33, pp. 8545–8552.
  * [25] Jiangliu Wang, Jianbo Jiao, Linchao Bao, Shengfeng He, Yunhui Liu, and Wei Liu,  “Self-supervised spatio-temporal representation learning for videos by predicting motion and appearance statistics,”  in CVPR, 2019, pp. 4006–4015.
  * [26] Jiangliu Wang, Jianbo Jiao, and Yun-Hui Liu,  “Self-supervised video representation learning by pace prediction,”  in ECCV. Springer, 2020, pp. 504–521.


# Spin Torque Generated by Valley Hall Effect in WSe2

D. J. P. de Sousa<EMAIL_ADDRESS>M. J. Sammon  Department of Electrical and
Computer Engineering, University of Minnesota, Minneapolis, Minnesota 55455,
USA  Raseong Kim  Hai Li  Ian A. Young  Components Research, Intel
Corporation, Hillsboro, OR 97124 USA  Tony Low<EMAIL_ADDRESS>Department of
Electrical and Computer Engineering, University of Minnesota, Minneapolis,
Minnesota 55455, USA

###### Abstract

Monolayer transition metal dichalcogenides are promising materials for
spintronics due to their robust spin-valley locked valence states, enabling
efficient charge-to-spin conversion via valley Hall effect with non-
equilibrium spins possessing long spin diffusion lengths of hundreds of
nanometers. In this work, we show that the injection of a pure valley current,
induced by valley Hall effect in a WSe2 monolayer, imparts a spin torque on
the magnetization of an overlaid Fe or CoFe in a tunneling structure. The
torque efficiency is found to be comparable to that in conventional
perpendicular magnetic tunnel junctions and can be further optimized with
valley Hall angle in WSe2. The valley nature of the spin torque gives rise to
out-of-plane damping-like torques in a current-in-plane configuration,
vanishing charge transport perpendicular-to-the-plane as well as torque
efficiency tunable through gating.

###### pacs:

71.10.Pm, 73.22.-f, 73.63.-b

_Introduction._ The ability to electrically manipulate the magnetization of a
ferromagnetic thin film with perpendicular magnetic anisotropy in an efficient
manner is envisioned to enable unprecedented technological advances through
the implementation of in-memory computing technologies [1, 2, 3, 4, 5, 6].
Deterministic field-free switching of a perpendicular magnetization requires
incident spin currents with a non-vanishing spin component aligned with the
perpendicular anisotropy axis. State-of-the-art approaches to the generation
of such out-of-plane spins include the charge-to-spin conversion through
unconventional spin Hall effects in low-symmetry non-magnets [7, 8, 9] and the
injection of spins by additional ferromagnet layers in fully perpendicular
magnetic tunnel junctions [1, 10, 11, 12, 13]. Recently, a new subfield of
spintronics has emerged [14, 15, 16, 17] which exploits the valley degree of
freedom of electrons in two-dimensional materials. Here, we examine the spin
torque generated through a valley polarization on an overlaid ferromagnet, and
its efficiency in switching a perpendicular magnetization. Addressing this
question is key in evaluating the potential of valley physics for spintronics.

In monolayer transition-metal dichalcogenides (TMDs), the spin-orbit
interaction induces a spin-valley locking, with $\rm{K}$ and $\rm{K}^{\prime}$
valleys supporting opposite spin states as required by time-reversal symmetry.
Spin-resolved photoemission spectroscopy measurements have revealed valley-
dependent out-of-plane spin polarized valence states in TMDs [18, 19, 20, 21,
22]. Subsequently, transport experiments unambiguously demonstrated that a
flow of out-of-plane spins in $\rm{WSe}_{2}$ monolayers can be electrically
generated through the valley Hall effect (VHE) [23], a topological Hall
response driven by the finite and opposite Berry curvatures in the two valleys
[23, 24, 25, 26, 27]. In these reports, an out-of-plane spin/valley
polarization of 70% was observed at the edges of p-doped monolayer
$\rm{WSe}_{2}$ due to VHE [27]. A similar electrically generated interfacial
spin polarization of 38% was also observed in $\rm{WSe}_{2}$/graphene
heterostructure [23]. Spin-valley locked states in TMDs are also long-lived
[28, 29] with large out-of-plane spin diffusion lengths of hundreds of
nanometers [25]. These features indicate the potential of spin-valley locking
physics to enable unprecedented applications in spintronics.

In this letter, we address the feasibility of utilizing the spin-valley
locking of TMDs to induce reversal of an adjacent perpendicular magnetization.
By performing transport calculations on a TMD/insulator/ferromagnetic tunnel
junction, we show the existence of a sizable spin torque (ST) acting on the
magnetization of a $\rm{Fe}$ or $\rm{CoFe}$ slab originating solely from a
non-equilibrium valley polarization in a $\rm{WSe}_{2}$ monolayer. We found
that the torque efficiency due to spin-valley locking depends sensitively on
the $\rm{WSe}_{2}$ doping levels and on the lattice misalignment with the
ferromagnet and is comparable to that in perpendicular magnetic tunnel
junction. The fact that the ST arises from the flow of out-of-plane spins
constitutes a novel mechanism suitable for switching a perpendicular
magnetization. Our findings point toward the utilization of spin-valley
locking physics of TMDs in novel efficient quantum spintronic devices.

_Theoretical method._ Figure 1(a) displays a
$\rm{WSe}_{2}$/insulator/$\rm{CoFe}$ stack as a prototypical example of the
tunneling junction studied in this work. In this system, a non-equilibrium
valley population induced in the WSe2 portion underneath the CoFe slab gives
rise to vertical tunneling valley currents. This situation can be realized by
positioning the CoFe slab on one of the side arms of the WSe2 crossbar and
taking advantage of the VHE to produce a local non-equilibrium valley density
[See Fig. 1(b)]. The associated non-equilibrium valley chemical potential
difference, $\delta\mu=\mu_{\rm{K}}-\mu_{\rm{K^{\prime}}}$, is well described
by the drift-diffusion approach [30], which renders

$\displaystyle\displaystyle\delta\mu=\frac{eI}{\sigma_{xx}}\frac{\theta_{H}}{1+\theta_{H}^{2}},$
(1)

where $I$ is the electric field-accompanying in-plane charge current and
$\theta_{H}$ the valley Hall angle of $\rm{WSe}_{2}$, defined as
$\theta_{H}=\sigma_{H}/\sigma_{xx}$ where $\sigma_{H}$ and $\sigma_{xx}$ are
the valley Hall and longitudinal charge conductivities, respectively. Our
drift-diffusion analysis indicates that Eq. (1) is the optimized valley
potential induced by VHE in a crossbar geometry. The optimization takes place
in the limit where the crossbar arm dimensions, $W_{1}$, $W_{2}$ and $L$,
become smaller than the spin diffusion length in WSe2, where $\delta\mu$
becomes independent of the crossbar dimensions. In the supplementary material,
we present numerical results showing the validity of Eq. (1) and how the
optimization comes about in the large spin/valley diffusion length limit [30].
The electric field is assumed to be applied along the $x$ direction and the
longitudinal charge conductivity is
$\sigma_{xx}=(4e^{2}/h)(\eta\beta)^{-1}\ln[1+\exp(\beta\mu_{\textrm{2D}})]$,
where $\mu_{\textrm{2D}}$ characterizes the doping levels of WSe2 and
$\beta=1/k_{B}T$, with $k_{B}$ being the Boltzmann constant and $T$ the
temperature [31]. The valley polarization is defined as
$(n_{K}-n_{K^{\prime}})/(n_{K}+n_{K^{\prime}})\times 100\%$ with valley-
dependent carrier density
$n_{K(K^{\prime})}=(m^{*}/\beta\pi\hbar^{2})\ln[1+\exp(\beta[\mu_{\textrm{2D}}\pm\delta\mu/2])]$
where $m^{*}$ is the effective mass of the WSe2 valence states.

Figure 1: (a) Representative $\rm{WSe}_{2}$/insulator/$\rm{CoFe}$ tunnel
junction with vertical tunneling valley/spin currents induced by an in-plane
electric field in the WSe2 monolayer. (b) The valley Hall effect causes a non-
equilibrium valley imbalance in the transverse side arms of a WSe2 crossbar in
response to a longitudinal in-plane electric field. The induced valley voltage
is optimized when the dimensions $W_{1}$, $W_{2}$ and $L$ are smaller than the
spin/valley diffusion length in WSe2. (c) Field-induced non-equilibrium valley
polarization in WSe2 as a function of the valley Hall angle at several
equilibrium WSe2 doping levels $\mu_{2D}$. Figure 2: The first principles
spin-resolved band structure (top panel) and spin density-of-states (bottom
panel) of a $\rm{WSe}_{2}$ monolayer and a 0.7 nm thick $\rm{CoFe}$ slab are
shown in panels (a) and (b), respectively. The horizontal dashed line in the
band structure plots indicate the energy used in the corresponding spin
density-of-states calculations. The valence states of $\rm{WSe}_{2}$ display a
valley-dependent out-of-plane spin polarization ($S_{z}$) that can be
extracted through valley Hall effect [23].

Figure 1(c) reveals that the total electrically generated valley polarization
induced at the side arms of the crossbar is a non-monotonic function of
$\theta_{H}$ with a maximum at $\theta_{H}\approx 1$ for all doping levels,
$\mu_{\textrm{2D}}$. Similar non-monotonic behavior was reported in Ref. [32].
Substantial valley polarizations are easily achieved at smaller
$\mu_{\textrm{2D}}$ due to the fact that valley depopulation takes place with
greater ease at lower $\mu_{\textrm{2D}}$ owing to the larger resistivity. The
vertical flow of such non-equilibrium valley-polarized electrons can be
addressed through a quantum mechanical tunneling approach, as described in the
following.

We employ the Bardeen transfer Hamiltonian formalism [33, 34] to describe the
tunneling process between the $\rm{WSe}_{2}$ monolayer and ferromagnetic thin
film. In Bardeen’s approach, quantum tunneling is treated perturbatively with
transition rates being fully described by the electronic ground state of
isolated contacts [35, 36]. The electronic states were obtained with the
pseudopotential/plane-wave method employed in QUANTUM ESPRESSO [37] and
subsequently converted to maximally-localized Wannier function basis with the
WANNIER90 package [38]. Figures 2(a) and (b) show the spin-resolved band
structure and spin density-of-states of a $\rm{WSe}_{2}$ monolayer and a 0.7
nm thick $\rm{CoFe}$ slab, respectively. The valence states of the
$\rm{WSe}_{2}$ monolayer display giant spin splittings of $\approx 468$ meV,
in agreement with previous results [18, 21], with nonequivalent valleys
$\rm{K}$ and $\rm{K}^{\prime}$ hosting opposite out-of-plane spin
polarizations. This strong spin-valley locking is better appreciated in the
bottom panel of Fig. 2(a), where we show the momentum-resolved spin density-
of-states for the energy indicated by the horizontal dashed line in the band
structure plot. The Fermi level momentum-resolved spin density-of-states of
the $\rm{CoFe}$ slab, shown at bottom panel of Fig. 2(b), indicates the
presence of electron states throughout the whole Brillouin zone including
momenta coinciding the $\rm{K}$ and $\rm{K}^{\prime}$ valleys. These finite
vertical tunneling currents require momentum matched states in the two layers.

Figure 3: (a) Energy-resolved spin torque acting on the magnetization of an Fe
(red) and CoFe (blue) slab. The energies are measured in relation to the
equilibrium chemical potential of the system $\mu_{0}$. (b) WSe2 chemical
potential dependence of the spin torque. Here, $\mu_{\textrm{2D}}$
characterizes the doping levels of WSe2 only, i.e., the chemical potential of
the ferromagnet is kept constant at $\mu_{0}$ in these calculations. The
energy-resolved tunneling charge-current passing through the insulating
barrier is shown in panel (c) for both systems. (d) Momentum-resolved non-
equilibrium chemical potential of the $\rm{WSe}_{2}$ monolayer.

The spin-space components of the tunneling charge-current density reads

$\displaystyle J^{ss^{\prime}}=$
$\displaystyle\displaystyle\frac{e}{\hbar}\int
d\epsilon\int\frac{d\textbf{k}_{\textrm{T}}}{(2\pi)^{2}}\int\frac{d\textbf{k}_{\textrm{B}}}{(2\pi)^{2}}[f(\textbf{k}_{\textrm{T}})-f(\textbf{k}_{\textrm{B}})]T^{ss^{\prime}}(\epsilon,\textbf{k}_{\textrm{T}},\textbf{k}_{\textrm{B}}),$
(2)

where
$f(\textbf{k}_{\textrm{T(B)}})=[1+\exp(\beta[\epsilon-\mu(\textbf{k}_{\textrm{T(B)}})])]^{-1}$
is the Fermi-Dirac distribution of the top (bottom) electrode and
$T^{ss^{\prime}}(\epsilon,\textbf{k}_{\textrm{T}},\textbf{k}_{\textrm{B}})$
are the spin space components of the momentum-dependent tunneling rates. The
momentum-dependent chemical potentials are
$\mu(\textbf{k}_{\textrm{B}})=\mu_{0}+\delta\mu(\textbf{k}_{\textrm{B}})$ and
$\mu(\textbf{k}_{\textrm{T}})=\mu_{0}$, where $\mu_{0}$ is the equilibrium
chemical potential of the whole structure and
$\delta\mu(\textbf{k}_{\textrm{B}})$ equals $+\delta\mu/2$ ($-\delta\mu/2$) at
the $\rm{K}$ ($\rm{K}^{\prime}$) valley due to the VHE. The function
$\delta\mu(\textbf{k}_{\textrm{B}})$ is related to the electronic properties
of $\rm{WSe}_{2}$ and applied in-plane current through Eq. (1).

The tunneling rate components appearing in Eq. (2) are

$\displaystyle
T^{ss^{\prime}}(\epsilon,\textbf{k}_{\textrm{T}},\textbf{k}_{\textrm{B}})=$
$\displaystyle\displaystyle
2\pi\sum_{n,m}|M_{nm}^{ss^{\prime}}(\textbf{k}_{\textrm{T}},\textbf{k}_{\textrm{B}})|^{2}\delta(\epsilon-\epsilon_{ns\textbf{k}_{\textrm{T}}})\delta(\epsilon-\epsilon_{ms^{\prime}\textbf{k}_{\textrm{B}}}),$
(3)

with matrix elements expressed in terms of the electronic eigenstates of
isolated top and bottom electrodes through

$\displaystyle
M_{nm}^{ss^{\prime}}(\textbf{k}_{\textrm{T}},\textbf{k}_{\textrm{B}})=$
$\displaystyle\displaystyle\frac{\hbar^{2}}{2m}\int_{z_{0}}dS\left(\frac{\partial\psi_{ns\textbf{k}_{\textrm{T}}}}{\partial
z}\psi_{ms^{\prime}\textbf{k}_{\textrm{B}}}^{\dagger}-\psi_{ns\textbf{k}_{\textrm{T}}}\frac{\partial\psi_{ms^{\prime}\textbf{k}_{\textrm{B}}}^{\dagger}}{\partial
z}\right),$ (4)

where the surface integral is performed halfway inside the barrier at $z_{0}$
and $m$ is the effective mass of electrons in the insulating region [30]. We
use the embedded Green’s function formalism for computing the matrix element
[34]. In Ref. [30] we show how this approach can be used to model tunneling
through an amorphous tunnel barrier with incommensurate contacts.

The total tunneling charge-current and out-of-plane polarized spin currents
are $J_{c}=\bar{J}+\Delta J\cos(\theta)$ and $Q^{z}=\bar{Q}^{z}+\Delta
Q^{z}\cos(\theta)$, respectively, where $\theta$ is the angle between the
magnetization of the ferromagnetic slab and the $z$ axis. We have defined
$\bar{J}=(J^{\uparrow\uparrow}+J^{\uparrow\downarrow}+J^{\downarrow\uparrow}+J^{\downarrow\downarrow})/2$,
$\Delta
J=(J^{\uparrow\uparrow}-J^{\uparrow\downarrow}-J^{\downarrow\uparrow}+J^{\downarrow\downarrow})/2$,
$\bar{Q}^{z}=(\hbar/4e)[J^{\uparrow\uparrow}+J^{\uparrow\downarrow}-J^{\downarrow\uparrow}-J^{\downarrow\downarrow}]$
and $\Delta
Q^{z}=(\hbar/4e)[J^{\uparrow\uparrow}-J^{\uparrow\downarrow}+J^{\downarrow\uparrow}-J^{\downarrow\downarrow}]$.
The damping-like torque acting on the magnetization of the ferromagnetic slab
is
$\boldsymbol{\tau}=\tau\hat{\textbf{m}}\times(\hat{\textbf{m}}\times\hat{\textbf{z}})$
where $\hat{\textbf{m}}$ is the unit vector along the magnetization direction
and $\tau=(Q^{z}(0)-Q^{z}(\pi))/2\rightarrow\tau=\Delta Q^{z}$ [39]. In the
following, we discuss the main features of the ST as well as the tunneling
charge and spin current originating from a valley polarization in the
$\textrm{WSe}_{2}$.

_Results and discussions._ We initially take the zigzag direction of the
$\rm{WSe}_{2}$ monolayer to be aligned with the [100] direction of the
$\rm{Fe}$ and $\rm{CoFe}$ slabs, an assumption that will be relaxed later. We
assume $\theta_{H}=0.4$ throughout [25]. Figure 3(a) displays the energy
dependence of the damping-like torque acting on the magnetization of a
$\rm{Fe}$ (red curve) or $\rm{CoFe}$ (blue curve) slab. We have taken the in-
plane charge current to be $I=1\mu$A and the ST was evaluated at
$\theta=\pi/2$, since its angular dependence is fully predetermined. As shown,
all contributions to the ST take place within a small energy window around the
equilibrium chemical potential, where states with energy above and below
$\mu_{0}$ contribute positively. Hence, a finite integrated ST acts on the
magnetization of the ferromagnetic slab originating solely from the non-
equilibrium valley-polarized electrons in the $\rm{WSe}_{2}$ monolayer. Figure
3(a) also indicates a relatively stronger (weaker) ST in the
$\rm{WSe}_{2}/insulator/\rm{Fe}(\rm{CoFe})$ tunnel junction, which we will
revisit later.

Figure 4: Dependence of the spin torque efficiency due to a lattice
misalignment. The lattice misalignment is quantified by the angle $\phi$
between the [100] direction of the ferromagnetic layer relative to the zigzag
direction of the $\rm{WSe}_{2}$ monolayer, which is held fixed along the $x$
direction. Panel (a) shows the angular dependence of the torque efficiency in
the $\rm{WSe}_{2}$/insulator/$\rm{CoFe}$ system while panels (b) and (c) show
the momentum-resolved density-of-states map of the $\rm{CoFe}$ at
$\phi=0^{\circ}$ and $\phi=7.5^{\circ}$, respectively. Panels (d), (e) and (f)
display similar results but for the $\rm{WSe}_{2}$/insulator/$\rm{Fe}$ system,
where the misalignment angles in the two maps are taken to be $\phi=0^{\circ}$
and $\phi=15^{\circ}$. For visualization purposes, we superposed the hexagonal
Brillouin zone of $\rm{WSe}_{2}$, highlighting the valleys.

The ST efficiency in a ferromagnetic element of area
$\textrm{A}_{\textrm{FM}}\approx 196\rm{nm}\times 196\rm{nm}$, defined as
$\xi=\textrm{A}_{\textrm{FM}}(e/\hbar)\tau/I$, due to a valley polarization of
70% is found to be $\xi_{\rm{CoFe}}\approx 0.07$ and $\xi_{\rm{Fe}}\approx
0.2$ for CoFe and Fe, respectively, for typical $\rm{WSe}_{2}$ hole
concentrations of $p\approx 2.75\times 10^{11}$ cm-2. For the sake of
comparison, the spin transfer torque efficiency in perpendicular magnetic
tunnel junctions range in the interval $0.01-0.1$ with a slight voltage
dependence, as estimated from the typical critical switching current density
$10^{6}$ A/cm2 and thermal stability factor of $\approx 60$ [13, 41, 40].
Therefore, such valley-induced torque efficiency is large enough to excite
magnetization dynamics or reversal of a perpendicular ferromagnetic thin film.
We also emphasize that the valley-induced torque efficiency can be further
optimized with the valley Hall angle and doping concentrations, as discussed
next.

Our results also indicate that the ST efficiency depends sensitively on the
doping levels of $\rm{WSe}_{2}$, $\mu_{\textrm{2D}}$, as measured from the
valence band top. This is shown in Fig. 3(b) for the $\rm{Fe}$ (red symbols)
and $\rm{CoFe}$ (blue symbols) cases. Here, the ST is a monotonically
decreasing function of $\mu_{\textrm{2D}}$ with a stronger dependence at lower
doping levels for both ferromagnetic slabs. This behavior is a signature of
the valley polarization dependence of the ST and can be understood in the low
temperature limit as follows: the Fermi-Dirac distribution can approximated as
$f(\epsilon-\mu_{\textrm{T(B)}})\approx\Theta(\mu_{\textrm{T(B)}}-\epsilon)$,
where $\Theta$ is the Heaviside function, leading to $\tau\propto\delta\mu$.
By explicitly writing $\delta\mu$ in terms of $\mu_{\textrm{2D}}$ through Eq.
(1) we find $\tau/I\propto 1/\mu_{\textrm{2D}}$, in agreement with the
behavior of Fig. 3(b). This result derives exclusively from the valley physics
of $\rm{WSe}_{2}$ monolayer and suggests that the ST efficiency can be
modulated through gating, a feature that has no analogue in conventional
spintronics.

Figure 3(c) shows that for both ferromagnetic slabs, states with energies
above and below $\mu_{0}$ contribute oppositely to the vertical total charge-
current density $J_{c}$. This implies in the scenario of Fig. 1(a); The ST
acting on the ferromagnet originates from an incoming flow of spin up
electrons and an outflow of spin down electrons such that the net tunneling
charge current vanishes while a net out-of-plane polarized spin current
penetrates the ferromagnet. This is possible by the non-equilibrium valley-
dependent chemical potential established in the $\rm{WSe}_{2}$ monolayer as
shown in Fig. 3(d), where the chemical potential of electrons with momenta at
the vicinity of the $\rm{K}$ ($\rm{K}^{\prime}$) is slightly lower (higher)
than that of the ferromagnet, leading to a vertical outflow (inflow) of
$\rm{K}$ ($\rm{K}^{\prime}$) valley-polarized electrons.

To further explore the nature and behavior of the ST, we study how its
magnitude is affected by a lattice misorientation. This is done by rotating
the [100] orientation of the ferromagnetic layer by an angle $\phi$ with
respect to the zigzag direction of the $\rm{WSe}_{2}$ monolayer, which is
maintained fixed along the $x$ direction. The angular dependence of the torque
efficiencies acting on the $\rm{CoFe}$ and $\rm{Fe}$ slabs are shown in Fig.
4(a) and (d), respectively. The results reveal an oscillatory behavior with a
30∘ period for both cases, where the maximum-to-minimum torque ratio,
$\tau_{\textrm{max}}/\tau_{\textrm{min}}$, is approximately $1.5$ and $3$ for
the $\rm{CoFe}$ and $\rm{Fe}$ cases, respectively. The stronger angular
dependence in the $\rm{Fe}$ case is due to more Fermi surface states at
momenta coinciding with the $\rm{K}$ and $\rm{K}^{\prime}$ valleys as compared
to $\rm{CoFe}$ system. This is shown in panels (b) and (e) of Fig. 4, where we
show the Fermi level momentum-resolved density-of-states of $\rm{CoFe}$ and
$\rm{Fe}$, respectively, at $\phi=0^{\circ}$. The superposed hexagonal
Brillouin zone of the $\rm{WSe}_{2}$ monolayer highlights the location of
$\rm{K}$ and $\rm{K}^{\prime}$ valleys. As seen in Fig. 4(e), the two valleys
located at $k_{y}=0$ Å -1 coincide with momentum space density-of-states hot
spots of $\rm{Fe}$, giving rise to a larger tunneling spin current. Such Fermi
surface matching is less ideal for the $\rm{CoFe}$ system. Figures 4(b) and
(e) show the hotspot misalignments giving rise to the maximum ST efficiency in
CoFe and Fe systems, respectively, while Figs. 4 (c) and (f) display the
configurations giving rise to the minimum ST efficiency [See Supplementary
material [30] for more discussion].

_Conclusions._ We have shown that a valley Hall effect-induced non-equilibrium
valley polarization in a $\rm{WSe}_{2}$ monolayer results in a spin torque
(ST) acting on the magnetization of a ferromagnetic thin film in
$\rm{WSe}_{2}$/insulator/ferromagnetic tunnel junctions. The valley-induced
damping-like torque arises from the flow of out-of-plane spin polarized
electrons and, therefore, is suitable for exciting magnetization dynamics of
thin films with perpendicular magnetic anisotropy, a novel approach that does
not require the presence of additional ferromagnetic elements or mirror-broken
spin Hall materials. The valley-induced ST was shown to display an efficiency
comparable to that in perpendicular magnetic tunnel junctions for typical
valley polarizations reported from experiments, with further room for
improvement through the valley Hall angle and doping concentrations in the
$\rm{WSe}_{2}$ monolayer, a feature that enables one to modulate the ST
efficiency through gating. Finally, we demonstrate how the lattice
misalignment and different magnetic layers affect the spin torque efficiency.

Acknowledgments. This material is based upon work supported by Intel
Corporation from the University Center (Valleytronics) program. We acknowledge
the Minnesota supercomputing institute (MSI) for providing the computational
resources. We thank Punyashloka Debashis at Intel Corporation for helpful
discussions.

## References

  * [1] H. Meng and J.-P. Wang, Appl. Phys. Lett. 88, 172506 (2006).
  * [2] S. Mangin, D. Ravelosona, J. A. Katine, M. J. Carey, B. D. Terris, and E. E. Fullerton, Nat. Mater. 5, 210 (2006).
  * [3] M. Gajek, J. J. Nowak, J. Z. Sun, P. L. Trouilloud, E. J. O’Sullivan, D. W. Abraham, M. C. Gaidis, G. Hu, S. Brown, Y. Zhu, et al., Appl. Phys. Lett 100, 132408 (2012).
  * [4] D. Zhang, M. Bapna, W. Jiang, D. Sousa, Y.-C. Liao, Z. Zhao, Y. Lv, P. Sahu, D. Lyu, A. Naeemi, Tony Low, S. A. Majetich, and J.-P. Wang, Nano Lett. 22, 2, 622 (2022).
  * [5] A. Sebastian, M. L. Gallo, R. K.-Aljameh and E. Eleftheriou, Nature Nanotechnology 15, 529 (2020).
  * [6] W. Zhang, B. Gao, J. Tang, P. Yao, S. Yu, M.-F. Chang, H.-J. Yoo, H. Qian and H. Wu, Nature Electronics 3, 371 (2020).
  * [7] D. MacNeill, G. M. Stiehl, M. H. D. Guimaraes, R. A. Buhrman, J. Park and D. C. Ralph, Nat. Phys. 13, 300 (2017).
  * [8] Fei Xue , C.Rohmann, J. Li, V. Amin, and P. M. Haney, Phys. Rev. B 102, 014401 (2020).
  * [9] S. Shi, J. Li, C.-H. Hsu, K. Lee, Y. Wang, L. Yang, J. Wang, Q.Wang, H. Wu, W. Zhang, G. Eda, G. Liang, H. Chang, H. Yang, Advanced Quantum Technologies 4, 2100038 (2021).
  * [10] H. Suto, T. Nagasawa, T. Kanao, K. Yamada and K. Mizushima, Sci. Rep. 9, 19543 (2019).
  * [11] D. J. P. de Sousa, P. M. Haney, D. L. Zhang, J. P. Wang, and Tony Low, Phys. Rev. B 101, 081404(R) (2020).
  * [12] H. Liu, D. Bedau, J. Z. Sun, S. Mangin, E. E. Fullerton, J. A. Katine, A. D. Kent, Journal of Magnetism and Magnetic Materials 358, 233 (2014).
  * [13] J. Z. Sun, IBM Journal of Research and Development, 50, 81 (2006).
  * [14] H. Yuan, M. S. Bahramy, K. Morimoto, S. Wu, K. Nomura, B.-J. Yang, H. Shimotani, R. Suzuki, M. Toh, C. Kloc, X. Xu, R. Arita, N. Nagaosa and Y. Iwasa, Nature Physics 9, 563 (2013).
  * [15] V. T. Renard, B. A. Piot, X. Waintal, G. Fleury, D. Cooper, Y. Niida, D. Tregurtha, A. Fujiwara, Y. Hirayama and K. Takashina, Nature Communications 6, 7230 (2015).
  * [16] L. A. Benítez, J. F. Sierra, W. S.Torres, A. Arrighi, F. Bonell, M. V. Costache and S. O. Valenzuela, Nature Physics 14, 303 (2018).
  * [17] R. Ahammed and A. De Sarkar, Phys. Rev. B 105, 045426 (2022).
  * [18] Z. Y. Zhu, Y. C. Cheng, and U. Schwingenschlögl Phys. Rev. B 84, 153402 (2011).
  * [19] R. Suzuki, M. Sakano, Y. J. Zhang, R. Akashi, D. Morikawa, A. Harasawa, K. Yaji, K. Kuroda, K. Miyamoto, T. Okuda, et al., Nat. Nanotech.y 9, 611 (2014).
  * [20] J. M. Riley, F. Mazzola, M. Dendzik, M. Michiardi, T. Takayama, L. Bawden, C. Granerød, M. Leandersson, T. Balasubramanian, M. Hoesch, et al., Nat. Phys., 10, 835 (2014).
  * [21] Di Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao Phys. Rev. Lett. 108, 196802 (2012).
  * [22] C. J. Ciccarino, T. Christensen, R. Sundararaman and P. Narang, Nano Lett. 18, 5709 (2018).
  * [23] T. Y. T. Hung, A. Rustagi, S. Zhang, P. Upadhyaya, Z. Chen, InfoMat. 2, 968 (2020).
  * [24] K. F. Mak, K. L. McGill, J. Park, P. L. McEuen, Science 344, 1489 (2014).
  * [25] Z. Wu, B. T. Zhou, X. Cai, .k Cheung, G.-B. Liu, M. Huang, J. Lin, T. Han, L. An, Y. Wang, et al., Nat. Commu. 10, 611 (2019).
  * [26] J. Lee, K. F. Mak and J. Shan, Nat. Nanotech. 11, 421 (2016).
  * [27] E. Barré, J. A. C. Incorvia, S. H. Kim, C. J. McClellan, E. Pop, H.-S. P. Wong, and T. F. Heinz Nano Letters 19 (2), 770 (2019).
  * [28] L. Yang, N. A. Sinitsyn, W. Chen, J. Yuan, J. Zhang, J. Lou and S. A. Crooker Nature Physics 11, 830 (2015).
  * [29] J. Kim, C. Jin, B. Chen, H. Cai, T. Zhao, P. Lee, S. Kahn, K. Watanabe, T. Taniguchi, S. Tongay, M. F. Crommie and F. Wang, Science Advances 3, e170051 (2017)
  * [30] See Supplemental Material at http://xxxx for more details on the theoretical model.
  * [31] In our calculations, we assume $T=1$ K for the temperature, $\eta$ = 0.1 eV for the spectral broadening and $I=1\mu$A for the in-plane electric field-accompanying charge current.
  * [32] H. Li, D. E. Nikonov, C.C. Lin, K. Camsari, Y.-C. Liao, C.-S. Hsu, A. Naeemi and I. A. Young, IEEE Journal on Exploratory Solid-State Computational Devices and Circuits,8, 10 (2022).
  * [33] J. Bardeen, Phys. Rev. Lett. 6, 57 (1961).
  * [34] D. Wortmann, H. Ishida, and S. Blügel, Phys. Rev. B 72, 235113 (2005).
  * [35] M. T. Greenaway, E. E. Vdovin, A. Mishchenko, O. Makarovsky, A. Patanè, J. R. Wallbank, Y. Cao, A. V. Kretinin, M. J. Zhu, S. V. Morozov, et al., Nat. Phys. 11, 1057 (2015).
  * [36] R. M. Feenstra, D. Jena and G. Gu, Journ. Appl. Phys. 111, 043711 (2012).
  * [37] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni et al., J. Phys.: Condens. Matter 29, 465901 (2017).
  * [38] A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 185, 2309 (2014).
  * [39] I. Theodonis, N. Kioussis, A. Kalitsov, M. Chshiev, and W. H. Butler, Phys. Rev. Lett. 97, 237205 (2006).
  * [40] In these estimations, we assume $\alpha=0.001$ is the intrinsic damping of the ferromagnetic slab.
  * [41] A. Kalitsov, W. Silvestre, M. Chshiev and J. P. Velev, Phys. rev. B 88, 104430 (2013).