11institutetext: University of Liège, Liège 4000, Belgium
11email<EMAIL_ADDRESS>

# Knowledge-Guided Additive Modeling
For Supervised Regression††thanks: This work was supported by Service Public
de Wallonie Recherche under Grant No. 2010235 - ARIAC by DIGITALWALLONIA4.AI.
Computational resources have been provided by the Consortium des Equipements
de Calcul Intensif (CECI), funded by the Fonds de la Recherche Scientifique de
Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region.

Yann Claes  0009-0005-2551-5152 Vân Anh Huynh-Thu  0000-0001-5492-2498 Pierre
Geurts  0000-0001-8527-5000

###### Abstract

Learning processes by exploiting restricted domain knowledge is an important
task across a plethora of scientific areas, with more and more hybrid methods
combining data-driven and model-based approaches. However, while such hybrid
methods have been tested in various scientific applications, they have been
mostly tested on dynamical systems, with only limited study about the
influence of each model component on global performance and parameter
identification. In this work, we assess the performance of hybrid modeling
against traditional machine learning methods on standard regression problems.
We compare, on both synthetic and real regression problems, several approaches
for training such hybrid models. We focus on hybrid methods that additively
combine a parametric physical term with a machine learning term and
investigate model-agnostic training procedures. We also introduce a new hybrid
approach based on partial dependence functions. Experiments are carried out
with different types of machine learning models, including tree-based models
and artificial neural networks111Our Python implementations of the hybrid
methods are available at https://github.com/yannclaes/kg-regression..

###### Keywords:

Knowledge-guided machine learning Physics-guided machine learning Supervised
regression Tree-based methods Neural networks Partial dependence plot Hybrid
modeling.

##  1 Introduction

For the past decades, machine learning (ML) models have been developed to
tackle a variety of real-life problems, complementing/replacing model-based
(MB) approaches, which mostly remain approximations that make stringent
assumptions about the system under study. Traditional ML approaches are said
to be data-driven, i.e. their prediction model is solely built from some
learning dataset, let it be (deep) neural networks or regression trees. While
their design comes with great expressiveness, they are likely to be subject to
over-fitting without enough training examples and to show a lack of robustness
on unseen samples, with predictions that can be inconsistent w.r.t. domain
knowledge [5, 6, 31]. To overcome this generalization issue, hybrid approaches
have been introduced to incorporate a priori domain knowledge within
statistical models, which can be leveraged in a multitude of ways (see [25,
16, 27] for reviews). The success of these hybrid methods have been shown
empirically on a range of synthetic and real-world problems [31, 1, 18, 7].
However, while these models have been mostly applied to dynamical systems,
they have not been thoroughly studied in the context of standard regression
problems. Furthermore, the majority of ML models that have been considered in
these approaches are neural networks and variants of the latter, leaving aside
other methods. Our contributions are the following:

  * •

We investigate empirically the performance and benefits of hybrid methods
against data-driven methods on static regression problems (in opposition to
dynamical problems). The static context removes a layer of complexity related
to the temporal correlation between observed states, which makes it easier to
assess the impact and interaction between the MB and ML components.
Specifically, we focus on hybrid models that combine in an additive way a
parametric physical term with an ML term.

  * •

We compare different approaches for training such hybrid additive models. We
highlight specific hypotheses under which these approaches are expected to
work well and relate the differences in terms of prediction and parameter
recovery performance. We focus on model-agnostic approaches, where the ML term
can be of any type, and we compare tree-based methods against neural networks.
Tree-based methods have several advantages over neural networks, which
motivate their use on static regression problems: they have much less
hyperparameters to tune, appear robust to the presence of irrelevant features
and have been shown to outperform neural networks on tabular data [12].

  * •

We introduce a new hybrid approach based on partial dependence functions,
which makes it easier to find the right balance between the MB and ML
components, makes less hypotheses than other approaches, and is shown to be
competitive in our experiments.

##  2 Problem statement

Let us define a regression problem, with $y\in\mathbb{R}$ and
$\mathbf{x}\in\mathbb{R}^{d}$, with $d\in\mathbb{N}_{+}$, drawn from a
distribution $p(\mathbf{x},y)$ such that $y=f(\mathbf{x})+\varepsilon$ with
$f:\mathbb{R}^{d}\mapsto\mathbb{R}$ the partially known generating function
and $\varepsilon\sim\mathcal{N}(0,\sigma^{2})$ the noise term. We focus on
problems such that $f(\mathbf{x})$ can be decomposed as:

######  Hypothesis 1 (H1, Additivity)

$y=f_{k}(\mathbf{x}_{k})+f_{a}(\mathbf{x})+\varepsilon,$

where $\mathbf{x}_{k}$ is a subset of $K\leq d$ input variables. We assume
partial knowledge of the generating function through some known algebraic
function $h_{k}^{\theta_{k}}(\mathbf{x}_{k})\in\mathcal{H}_{k}$ with tunable
parameters $\theta_{k}$, such that for the optimal parameters $\theta_{k}^{*}$
we have $h_{k}^{\theta_{k}^{*}}=f_{k}$. The residual term $f_{a}(\mathbf{x})$
is unknown and is approximated in this work through an ML component
$h_{a}^{\theta_{a}}\in\mathcal{H}_{a}$, with parameters $\theta_{a}$222In the
following, $h_{k}^{\theta_{k}}$ and $h_{a}^{\theta_{a}}$ will sometimes be
denoted simply as $h_{k}$ and $h_{a}$ to lighten the notations.. The final
model $h\in{\cal H}$ is denoted
$h(\mathbf{x})=h_{k}^{\theta_{k}}(\mathbf{x}_{k})+h_{a}^{\theta_{a}}(\mathbf{x})$,
with the function space $\cal H$ defined as ${\cal H}_{k}+{\cal H}_{a}$. H1 is
common when MB methods and ML models are combined [22, 30, 7, 26].

Given a learning sample of $N$ input-output pairs
$LS=\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{N}$, drawn from $p(\mathbf{x},y)$, we
seek to identify a function $h=h_{k}^{\theta_{k}}+h_{a}^{\theta_{a}}$, i.e.
parameters $\theta_{k}$ and $\theta_{a}$, that minimizes the following two
distances:

$\displaystyle d(h,y)$ $\displaystyle=\mathbb{E}_{(\mathbf{x},y)\sim
p(\mathbf{x},y)}\\{(h(\mathbf{x})-y)^{2}\\},$ (1) $\displaystyle
d_{k}(h^{\theta_{k}}_{k},f_{k})$ $\displaystyle=\mathbb{E}_{\mathbf{x}_{k}\sim
p(\mathbf{x}_{k})}\\{(h^{\theta_{k}}_{k}(\mathbf{x}_{k})-f_{k}(\mathbf{x}_{k}))^{2}\\}.$
(2)

The first distance measures the standard generalization error of the global
model $h$. The hope is that taking $h_{k}$ into account will help learning a
better global model than fitting directly a pure data-driven model on $y$,
especially in the small sample size regime. The second distance $d_{k}$
measures how well the tuned $h_{k}$ approximates $f_{k}$. The main motivation
for this second objective is interpretability: one expects that the algebraic
form of $h_{k}$ will be derived from first principles by domain experts, who
will be interested in estimating the parameters of this term from data. An
alternative to $d_{k}$ is a loss that would compare the estimated and optimal
parameters $\hat{\theta}_{k}$ and $\theta^{*}_{k}$ (e.g.,
$||\hat{\theta}_{k}-\theta^{*}_{k}||^{2}$). $d_{k}$ however has the advantage
not to require $\theta_{k}^{*}$ to be fully identifiable, i.e. there can exist
several sets of parameters $\theta^{*}_{k}$ such that
$h^{\theta^{*}_{k}}_{k}=f_{k}$. In our experiment, we will report both $d_{k}$
and the relative mean absolute error on the estimated parameters.

The following approximation of (1) can be used as training objective:

$\hat{d}(h,y;LS)=\frac{1}{N}\sum_{i=1}^{N}(h(\mathbf{x}_{i})-y_{i})^{2}.$ (3)

Minimizing the distance in (2) is expected to be challenging and sometimes
even ill-posed. Indeed, if $h_{a}$ is too powerful, it could capture $f$
entirely and leave little room for the estimation of $f_{k}$. Finding the
right balance between $h_{k}$ and $h_{a}$ is thus very challenging, if not
impossible, using only guidance of the learning sample $LS$. Unlike (1), (2)
cannot be estimated from a sample of input-output pairs and hence cannot be
explicitly used to guide model training. There are however several scenarios
that will make the problem easier. In the following, we will discuss the
optimality of the hybrid methods under two additional assumptions:

######  Hypothesis 2 (H2, Disjoint features)

Let $\mathbf{x}_{a}$ be a subset of features disjoint from $\mathbf{x}_{k}$
($\mathbf{x}_{k}\cap\mathbf{x}_{a}=\emptyset$). There exists a function
$f^{r}_{a}(\mathbf{x}_{a})$ such that
$f_{a}(\mathbf{x})=f^{r}_{a}(\mathbf{x}_{a})$ for all $\mathbf{x}$.

######  Hypothesis 3 (H3, Independence)

Features in $\mathbf{x}_{k}$ are independent from features in $\mathbf{x}_{a}$
($\mathbf{x}_{k}\perp\\!\\!\\!\perp\mathbf{x}_{a}$).

H2 makes the problem easier as $f_{k}$ captures all the dependence of $y$ on
$\mathbf{x}_{k}$. In the absence of H3, it might be hard to distinguish real
contributions from $\mathbf{x}_{k}$ to $f$ from those due to correlations with
features not in $\mathbf{x}_{k}$.

##  3 Related work

Hybrid additive modeling methods emerged several decades ago, combining first-
principles models with different ML models. Already in the 1990’s, approaches
in [17, 14, 24] complemented physics-based models with neural networks,
weighting contributions of both components (e.g. through radial basis function
networks), to achieve enhanced physical consistency with better generalization
properties. More recently, other works applied the same principles to model
dynamical systems in various domains, still massively relying on neural
networks [30, 21, 22, 26]. In a more standard regression setting, [3, 32, 33]
combined a linear parametric term with a tree-based ML term.

Previous works have introduced regularization of the ML term to reduce
parameter identification issues in the decomposition [28, 15, 13]. Further
works on this matter introduced physically-motivated constraints in the
learning objective to better control contributions of the MB/ML components
[31, 7]. Elements of discussion about the well-posedness of this additive
decomposition have been introduced in previous works: [31] showed the
existence and uniqueness of an optimal pair ($h_{k}$, $h_{a}$) when the
contributions of $h_{a}$ are constrained to be minimal, and [7] demonstrated
the convergence of an algorithm alternating between the optimization of
$h_{k}$ and the optimization of $h_{a}$, without however any guarantee about
convergence points.

##  4 Methods

We focus on model-agnostic approaches, i.e. that can be applied with any
algebraic function $h_{k}$ and any type of ML model $h_{a}$. For both terms,
we only assume access to training functions, respectively denoted
$\operatorname{fit}^{h_{k}}$, $\operatorname{fit}^{h_{k}+\gamma}$, and
$\operatorname{fit}^{h_{a}}$, that can estimate each model parameters,
respectively $\theta_{k}$, $(\theta_{k},\gamma)$ and $\theta_{a}$, so as to
minimize the mean squared error (MSE) over $LS$ (see below for the meaning of
$\gamma$), where parametric methods rely on gradient descent.

###  4.1 Sequential training of $h_{k}$ and $h_{a}$

This baseline approach first fits $h_{k}$ on the observed output $y$, then
fits $h_{a}$ on the resulting residuals, as done in [33]. More precisely, we
first train $h_{k}^{\theta_{k}}$ on $y$ by introducing a constant term
$\gamma\in\mathbb{R}$, such that

$(\hat{\theta}_{k},\hat{\gamma})=\operatorname{fit}^{h_{k}+\gamma}(LS).$ (4)

Our motivation for introducing the term $\gamma$ will be explained below.
Afterwards, we fit $h_{a}$ on the output residuals:
$\hat{\theta}_{a}=\operatorname{fit}^{h_{a}}\\{(\mathbf{x}_{i},y_{i}-h_{k}^{\hat{\theta}_{k}}(\mathbf{x}_{i})-\hat{\gamma})\\}_{i=1}^{N}$.

Let $\hat{\mathcal{F}}_{k}$ be the set of all functions $\hat{f}_{k}$ mapping
$\mathbf{x}_{k}\in\mathcal{X}_{k}$ to some value $y\in\mathbb{R}$, i.e.
$\hat{\mathcal{F}}_{k}=\\{\hat{f}_{k}:\mathcal{X}_{k}\mapsto\mathbb{R}\\}$.
Under H2 and H3, it can be shown that
$\hat{f}_{k}^{*}=\arg\min_{\hat{f}_{k}\in\hat{\mathcal{F}}_{k}}d(\hat{f}_{k},y)$
is such that $\hat{f}_{k}^{*}(\mathbf{x}_{k})=f_{k}(\mathbf{x}_{k})+C$, for
every $\mathbf{x}_{k}\in\mathcal{X}_{k}$, with
$C=\mathbb{E}_{\mathbf{x}_{a}}\left\\{f^{r}_{a}(\mathbf{x}_{a})\right\\}$ (see
Appendix 0.A). Hence, this approach is sound at least asymptotically and
justifies the introduction of $\gamma$. Note however that even under H2 and
H3, we have no guarantee that this approach produces the best estimator for a
finite sample size, as $f^{r}_{a}(\mathbf{x}_{a})+\epsilon$ acts as a pure
additive noise term that needs to be averaged out during training. The
approaches described in Sections 4.2 and 4.3 try to overcome this issue by
fitting $h^{\theta_{k}}_{k}$ on corrected outputs that are expected to be
closer to $f_{k}(\mathbf{x}_{k})$. Without H2 and H3, the quality of the
estimation of $f_{k}$ by $h_{k}^{\hat{\theta}_{k}}$, according to (2), is not
guaranteed as there are regression problems satisfying H1 such that:

$\nexists\gamma\in\mathbb{R}:\arg\min_{\hat{f}_{k}\in\hat{\mathcal{F}}_{k}}d(\hat{f}_{k},y)=f_{k}+\gamma.$
(5)

An example will be given in Section 5.2.

###  4.2 Alternate training of $h_{k}$ and $h_{a}$

A hybrid additive approach was proposed in [7] that alternates between
updating $h_{k}$ and updating $h_{a}$, using neural networks for $h_{a}$. Such
alternate training was also proposed in [3, 32] with a single decision tree as
$h_{a}$ and a linear $h_{k}$. We include this approach in our comparison, but
also investigate it with random forests [2] and tree gradient boosting [9].
$\hat{\theta}_{k}$ is initialized by (fully) fitting
$h_{k}^{{\theta}_{k}}+\gamma$ on $y$. Then, we alternate between: (1) a single
epoch of gradient descent on $h_{k}^{{\theta}_{k}}+\gamma$ and (2) either a
single epoch for $h_{a}$ (in the case of neural networks, as in [7]) or a
complete fit of $h_{a}$ (in the case of tree-based models).

While some theoretical results are provided in [7], convergence of the
alternate method towards the optimal solution is not guaranteed in general.
Despite an initialization favoring $h_{k}$, it is unclear whether a too
expressive $h_{a}$ will not dominate $h_{k}$ and finding the right balance
between these two terms, e.g. by regularizing further $h_{a}$, is challenging.
Under H2 and H3 however, the population version of the algorithm produces an
optimal solution. Indeed, $h_{k}$ will be initialized as the true $f_{k}$, as
shown previously, making the residuals $y-h_{k}$ at the first iteration, as
well as $h_{a}$, independent of $\mathbf{x}_{k}$. $h_{k}$ will thus remain
unchanged (and optimal) at subsequent iterations.

###  4.3 Partial Dependence-based training of $h_{k}$ and $h_{a}$

We propose a novel approach relying on partial dependence (PD) functions [8]
to produce a proxy dataset depending only on $\mathbf{x}_{k}$ to fit $h_{k}$.
PD measures how a given subset of features impact the prediction of a model,
on average. Let $\mathbf{x}_{k}$ be the subset of interest and
$\mathbf{x}_{-k}$ its complement, with
$\mathbf{x}_{k}\cup\mathbf{x}_{-k}=\mathbf{x}$, then the PD of a function
$f(\mathbf{x})$ on $\mathbf{x}_{k}$ is:

$PD(f,\mathbf{x}_{k})=\mathbb{E}_{\mathbf{x}_{-k}}\left[f(\mathbf{x}_{k},\mathbf{x}_{-k})\right]=\int
f(\mathbf{x}_{k},\mathbf{x}_{-k})p(\mathbf{x}_{-k})d\mathbf{x}_{-k},$ (6)

where $p(\mathbf{x}_{-k})$ is the marginal distribution of $\mathbf{x}_{-k}$.
Under H1 and H2, the PD of
$f(\mathbf{x})=f_{k}(\mathbf{x}_{k})+f^{r}_{a}(\mathbf{x}_{a})$ is [8]:

$PD(f,\mathbf{x}_{k})=f_{k}(\mathbf{x}_{k})+C,\mbox{ with
}C=E_{\mathbf{x}_{a}}\\{f^{r}_{a}(\mathbf{x}_{a})\\}.$ (7)

The idea of our method is to first fit any sufficiently expressive ML model
$h_{a}(\mathbf{x})$ on $LS$ and to compute its PD w.r.t. $\mathbf{x}_{k}$ to
obtain a first approximation of $f_{k}(\mathbf{x}_{k})$ (up to a constant).
Although computing the actual PD of a function using (6) requires in principle
access to the input distribution, an approximation can be estimated from $LS$
as follows:

$\widehat{PD}(h_{a},\mathbf{x}_{k};LS)=\frac{1}{N}\sum_{i=1}^{N}h_{a}(\mathbf{x}_{k},\mathbf{x}_{i,-k}),$
(8)

where $\mathbf{x}_{i,-k}$ denotes the values of $\mathbf{x}_{-k}$ in the
$i$-th sample of $LS$. A new dataset of pairs
$(\mathbf{x}_{k},\widehat{PD}(h_{a},\mathbf{x}_{k};LS))$ can then be built to
fit $h_{k}$. In our experiments, we consider only the $\mathbf{x}_{k}$ values
observed in the learning sample but $\widehat{PD}(h_{a},\mathbf{x}_{k};LS)$
could also be estimated at other points $\mathbf{x}_{k}$ to artificially
increase the size of the proxy dataset.

In practice, optimizing $\theta_{k}$ only once on the PD of $h_{a}$ could
leave residual dependence of $\mathbf{x}_{k}$ on the resulting
$y-h_{k}^{\hat{\theta}_{k}}(\mathbf{x}_{k})-\hat{\gamma}$. We thus repeat the
sequence of fitting $h_{a}$ on the latter residuals, then fitting $h_{k}$ on
the obtained
$\widehat{PD}(h_{a}^{\hat{\theta}_{a}},\mathbf{x}_{k};LS)+h_{k}^{\hat{\theta}_{k}}(\mathbf{x}_{k})+\hat{\gamma}$,
with $\hat{\theta}_{k}$ and $\hat{\theta}_{a}$ the current optimized parameter
vectors (see Algorithm 1).

Algorithm 1 Partial Dependence Optimization

Input: $LS={(\mathbf{x}_{i},y_{i})}_{i=1}^{N}$

$\hat{\theta}_{a}\leftarrow\operatorname{fit}^{h_{a}}(LS)$

$(\hat{\theta}_{k},\hat{\gamma})\leftarrow\operatorname{fit}^{h_{k}+\gamma}(\\{(\mathbf{x}_{k,i},\widehat{PD}(h_{a}^{\hat{\theta}_{a}},\mathbf{x}_{k,i};LS))\\}_{i=1}^{N})$

for $n=1$ to $N_{repeats}$ do

$\hat{\theta}_{a}\leftarrow\operatorname{fit}^{h_{a}}(\\{(\mathbf{x}_{i},y_{i}-h_{k}^{\hat{\theta}_{k}}(\mathbf{x}_{k,i})-\hat{\gamma})\\}_{i=1}^{N})$

$(\hat{\theta}_{k},\hat{\gamma})\leftarrow\operatorname{fit}^{h_{k}+\gamma}(\\{(\mathbf{x}_{k,i},h_{k}^{\hat{\theta}_{k}}(\mathbf{x}_{k,i})+\hat{\gamma}+\widehat{PD}(h_{a}^{\hat{\theta}_{a}},\mathbf{x}_{k,i};LS))\\}_{i=1}^{N})$

end for

$\hat{\theta}_{a}\leftarrow\operatorname{fit}^{h_{a}}(\\{(\mathbf{x}_{i},y_{i}-h_{k}^{\hat{\theta}_{k}}(\mathbf{x}_{k,i})-\hat{\gamma})\\}_{i=1}^{N})$

The main advantage of this approach over the alternate one is to avoid
domination of $h_{a}$ over $h_{k}$. Unlike the two previous approaches, this
one is also sound even if H3 is not satisfied as it is not a requirement for
(7) to hold. One drawback is that it requires $h_{a}$ to capture well the
dependence of $f$ on $\mathbf{x}_{k}$ so that its PD is a good approximation
of $f_{k}$. The hope is that even if it is not the case at the first
iteration, fitting $h_{k}$, that contains the right inductive bias, will make
the estimates better and better over the iterations.

##  5 Experiments

We compare the different methods on several regression datasets, both
simulated and real. Performance is measured through estimates of (1) and (2)
(the latter only on simulated problems) on a test set $TS$, respectively
denoted $\hat{d}(h,y;TS)$ and $\hat{d}_{k}(h_{k}^{\theta_{k}},f_{k};TS)$. In
some cases, we also report $\operatorname{rMAE}(\theta_{k}^{*},\theta_{k})$,
the relative mean absolute error between $\theta_{k}^{*}$ and $\theta_{k}$
(lower is better for all measures). For the hybrid approaches, we use as
$h_{a}$ either a multilayer perceptron (MLP), gradient boosting with decision
trees (GB) or random forests (RF). We compare these hybrid models to a
standard data-driven model that uses only $h_{a}$. We also compare fitting
$h_{a}$ with and without input filtering, i.e. respectively removing or
keeping $\mathbf{x}_{k}$ from its inputs, to verify convergence claims about
$h_{k}$ in Section 4.2. Architectures (e.g. for MLP, the number of layers and
neurons) are kept fixed across training methods to allow a fair comparison
between them, and are given in Appendix 0.B. We use early stopping of gradient
descent training by monitoring the loss on a validation set (except for pure
tree-based models, which are trained in the standard way, hence not using
gradient descent).

###  5.1 Friedman problem (H2 and H3 satisfied)

We consider the following synthetic regression problem:

$y=\theta_{0}\sin(\theta_{1}x_{0}x_{1})+\theta_{2}(x_{2}-\theta_{3})^{2}+\theta_{4}x_{3}+\theta_{5}x_{4}+\varepsilon,$

where $x_{j}\sim\mathcal{U}(0,1),j=0,\dots 9$, and
$\varepsilon\sim\mathcal{N}(0,1)$ [10]. We generate 10 different datasets
using 10 different sets of values for $\theta_{0},\dots,\theta_{5}$, each with
300, 300 and 600 samples for respectively the training, validation and test
sets. For the hybrid approaches, we use the first term as prior knowledge,
i.e. $f_{k}(\mathbf{x}_{k})=\theta_{0}\sin(\theta_{1}x_{0}x_{1})$.

Table 1: Results on the Friedman problem. We report the mean and standard deviation of $\hat{d}$ and $\hat{d}_{k}$ over the 10 test sets (TS). “$f_{k}\rightarrow h_{a}$” fits $h_{a}$ on $y-f_{k}(\mathbf{x}_{k})$. “Unfiltered” indicates that all the features are used as inputs of $h_{a}$, while “Filtered” indicates that the features $\mathbf{x}_{k}$ are removed from the inputs of $h_{a}$. |  | $\hat{d}(h,y;TS)$ | $\hat{d}_{k}(h_{k}^{\theta_{k}},f_{k};TS)$
---|---|---|---
| Method | Unfiltered | Filtered | Unfiltered | Filtered
| $f_{k}\rightarrow h_{a}$ | $1.58\pm 0.33$ | $1.23\pm 0.10$ | -
| Sequential | $1.54\pm 0.31$ | $1.43\pm 0.13$ | $0.18\pm 0.16$
MLP | Alternate | $1.43\pm 0.09$ | $1.32\pm 0.09$ | $0.10\pm 0.09$ | $0.02\pm 0.02$
| PD-based | $1.54\pm 0.12$ | $1.38\pm 0.09$ | $0.06\pm 0.07$
|  $h_{a}$ only | $2.62\pm 0.75$ | -
| $f_{k}\rightarrow h_{a}$ | $1.73\pm 0.09$ | $1.75\pm 0.12$ | -
| Sequential | $1.74\pm 0.11$ | $1.81\pm 0.14$ | $0.18\pm 0.16$
GB | Alternate | $1.79\pm 0.11$ | $1.78\pm 0.15$ | $0.91\pm 1.45$ | $0.06\pm 0.06$
| PD-based | $1.77\pm 0.13$ | $1.78\pm 0.12$ | $0.03\pm 0.02$
|  $h_{a}$ only | $3.43\pm 0.94$ | -
| $f_{k}\rightarrow h_{a}$ | $2.03\pm 0.18$ | $1.96\pm 0.17$ | -
| Sequential | $2.11\pm 0.23$ | $2.05\pm 0.24$ | $0.18\pm 0.16$
RF | Alternate | $2.03\pm 0.19$ | $1.98\pm 0.17$ | $0.04\pm 0.03$ | $0.04\pm 0.04$
| PD-based | $2.16\pm 0.27$ | $2.09\pm 0.26$ | $0.16\pm 0.15$
|  $h_{a}$ only | $5.58\pm 1.91$ | -

We see in Table 1 that all hybrid training schemes outperform their data-
driven counterpart. They come very close to the ideal $f_{k}\rightarrow h_{a}$
method, and sometimes even slightly better, probably due to small overfitting
issues. Sequential fitting of $h_{k}$ and $h_{a}$ performs as well as the
alternate or PD-based approaches, as H2 and H3 are satisfied for this problem
(see Section 4.1). Filtering generally improves the performance of hybrid
schemes as H2 is verified. PD-based optimization yields good approximations of
$f_{k}$ (as shown by a low $\hat{d}_{k}$). The alternate approach follows
closely whereas the sequential one ends up last, which can be expected as
fitting $h_{k}$ only on $y$ induces a higher noise level centered around
$\mathbb{E}_{\mathbf{x}_{a}}\left\\{f_{a}(\mathbf{x}_{a})\right\\}$, while the
other approaches benefit from reduced perturbations through $h_{a}$
estimation, as explained in Section 4.1. Filtering vastly decreases
$\hat{d}_{k}$ for alternate approaches, supporting claims introduced in
Section 4.2, while this measure remains unimpaired for sequential and PD-based
training by construction.

###  5.2 Correlated input features (H3 not satisfied)

#### Correlated linear model.

Let $y=\beta_{0}x_{0}+\beta_{1}x_{1}+\varepsilon$, with
$\beta_{0}=-0.5,\beta_{1}=1,\mathbf{x}\sim\mathcal{N}(\mathbf{0},\Sigma),$ and
$\varepsilon\sim\mathcal{N}(0.5^{2},1)$. We generate 50, 50 and 600 samples
respectively for the training, validation and test sets. We use as known term
$f_{k}(\mathbf{x}_{k})=\beta_{0}x_{0}$. Regressing $y$ on $x_{0}$ yields the
least-squares solution [11]:

$\mathbb{E}\left[\hat{\beta}_{0}\right]=\beta_{0}+\frac{\operatorname{cov}(x_{0},x_{1})}{\operatorname{var}(x_{0})}\beta_{1}.$
(9)

We set $\operatorname{cov}(x_{0},x_{1})=2.25$ and
$\operatorname{var}(x_{0})=2$ so that (9) reverses the sign of $\beta_{0}$ and
(5) is satisfied. The sequential approach should hence yield parameter
estimates of $\beta_{0}$ close to (9) while we expect the others to correct
for this bias.

Table 2: Results for the correlated linear problem. We report $\hat{d}$ and $\operatorname{rMAE}(\beta_{0}^{*},\hat{\beta}_{0})$, over 10 different datasets. | $\hat{d}(h,y;TS)$ | $\operatorname{rMAE}(\beta_{0}^{*},\hat{\beta}_{0})$
---|---|---
Method | Unfiltered | Filtered | Unfiltered | Filtered
Seq. (MLP) | $0.30\pm 0.03$ | $0.74\pm 0.09$ | $224.14\pm 13.48$
Alt. (MLP) | $0.30\pm 0.02$ | $0.31\pm 0.04$ | $186.65\pm 21.31$ | $15.53\pm 13.57$
PD (MLP) | $0.30\pm 0.03$ | $0.29\pm 0.02$ | $26.47\pm 17.32$
Seq. (GB) | $0.59\pm 0.06$ | $1.38\pm 0.11$ | $224.14\pm 13.48$
Alt. (GB) | $0.57\pm 0.06$ | $0.60\pm 0.09$ | $148.75\pm 67.35$ | $24.58\pm 12.20$
PD. (GB) | $0.56\pm 0.05$ | $0.64\pm 0.13$ | $36.05\pm 17.50$
Seq. (RF) | $0.53\pm 0.05$ | $0.90\pm 0.07$ | $224.14\pm 13.48$
Alt. (RF) | $0.43\pm 0.04$ | $0.42\pm 0.04$ | $111.04\pm 52.78$ | $45.38\pm 22.39$
PD (RF) | $0.41\pm 0.03$ | $0.43\pm 0.04$ | $57.47\pm 15.55$

From Table 2, we observe that, contrary to the PD-based approach, the
sequential and alternate methods return very bad estimations of $\beta_{0}$,
as H3 is no longer verified. Filtering corrects the bias for the alternate
approach but degrades the MSE performance for the sequential method as it
removes the ability to compensate for the $h_{k}$ misfit.

#### Correlated Friedman problem.

The structure is identical to the one in Section 5.1 but with correlated
inputs drawn from a multivariate normal distribution where $\mu_{i}=0.5$ and
$\operatorname{var}(x_{i})=0.75,\forall i$, and
$\operatorname{cov}(x_{i},x_{j})=\pm 0.3,\forall i\neq j$ (the covariance sign
being chosen randomly). Sizes of the training, validation and test sets are
identical to those of Section 5.1. Inputs are then scaled to be roughly in
$[-1,1]$. Here again, we use
$f_{k}(\mathbf{x}_{k})=\theta_{0}\sin(\theta_{1}x_{0}x_{1})$.

Table 3: Results for the correlated Friedman problem. |  | $\hat{d}(h,y;TS)$ | $\hat{d}_{k}(h_{k}^{\theta_{k}},f_{k};TS)$
---|---|---|---
| Method | Unfiltered | Filtered | Unfiltered | Filtered
| $f_{k}\rightarrow h_{a}$ | $1.64\pm 0.23$ | $1.51\pm 0.17$ | -
| Sequential | $2.07\pm 0.40$ | $2.68\pm 1.38$ | $1.35\pm 1.42$
MLP | Alternate | $1.95\pm 0.33$ | $1.62\pm 0.24$ | $0.49\pm 0.44$ | $0.14\pm 0.19$
| PD-based | $2.24\pm 0.31$ | $1.78\pm 0.30$ | $0.17\pm 0.23$
|  $h_{a}$ only | $2.77\pm 0.73$ | -
| $f_{k}\rightarrow h_{a}$ | $2.58\pm 0.45$ | $2.53\pm 0.44$ | -
| Sequential | $2.90\pm 0.39$ | $3.91\pm 1.49$ | $1.35\pm 1.42$
GB | Alternating | $2.67\pm 0.38$ | $2.62\pm 0.43$ | $0.51\pm 0.53$ | $0.22\pm 0.25$
| PD-based | $2.54\pm 0.35$ | $2.47\pm 0.36$ | $0.03\pm 0.02$
|  $h_{a}$ only | $4.49\pm 0.66$ | -
| $f_{k}\rightarrow h_{a}$ | $3.02\pm 0.45$ | $2.93\pm 0.45$ | -
| Sequential | $3.78\pm 0.78$ | $4.04\pm 1.30$ | $1.35\pm 1.42$
RF | Alternating | $3.06\pm 0.39$ | $2.99\pm 0.38$ | $0.14\pm 0.16$ | $0.15\pm 0.18$
| PD-based | $3.24\pm 0.38$ | $3.16\pm 0.37$ | $0.27\pm 0.20$
|  $h_{a}$ only | $6.70\pm 1.47$ | -

As in Section 5.1, Table 3 shows that hybrid models outperform their data-
driven equivalents. PD-based methods usually yield more robust $h_{k}$
estimations in the general unfiltered case, but struggle to line up with the
alternate scheme in terms of predictive performance, except for GB-related
models. For RF, this can be explained by a worse $h_{k}$ estimation while for
MLP we assume that it is due to $h_{a}$ overfitting: in the alternate
approach, it is optimized one epoch at a time, interleaved with one step on
$h_{k}$, whereas that of PD-based methods is fully optimized (with identical
complexities). Sequential and alternate approaches undergo stronger $h_{k}$
misparameterization without filtering since H3 is not met, but the latter
mitigates this w.r.t. the former, as was already observed in Section 5.1.
Input filtering degrades predictive performance for the sequential methods as
they cannot counterbalance a poor $h_{k}$.

###  5.3 Overlapping additive structure (H2 and H3 not satisfied)

Let $y=\beta x_{0}^{2}+\sin(\gamma x_{0})+\delta x_{1}+\varepsilon$ with
$\varepsilon\sim\mathcal{N}(0,0.5^{2}),\beta=0.2,\gamma=1.5,\delta=1$ and
$\mathbf{x}$ sampled as in the correlated linear problem. We generate 50, 50
and 600 samples respectively for the training, validation and test sets. We
define $f_{k}(\mathbf{x}_{k})=\beta x_{0}^{2}$ and
$f_{a}(\mathbf{x})=\sin(\gamma x_{0})+\delta x_{1}+\varepsilon$. Hence, H2 and
H3 do not hold. Even with $\hat{\beta}=\beta^{*}$, $h_{a}$ still needs to
compensate for $\sin(\gamma x_{0})$. Filtering is thus expected to degrade
performance for all hybrid approaches as $h_{a}(x_{1})$ will never compensate
this gap, which is observed in Table 4. Results for RF are not shown for the
sake of space, but are similar to GB.

Table 4: Results for the overlapping problem. | $\hat{d}(h,y;TS)$ | $\hat{d}(h,y;TS)$
---|---|---
| Unfiltered | Filtered | Unfiltered | Filtered
Method | MLP | GB
$f_{k}\rightarrow h_{a}$ | $0.35\pm 0.02$ | $0.54\pm 0.04$ | $0.51\pm 0.04$ | $1.00\pm 0.12$
Sequential | $0.35\pm 0.01$ | $0.59\pm 0.05$ | $0.55\pm 0.07$ | $1.07\pm 0.11$
Alternate | $0.35\pm 0.02$ | $0.56\pm 0.05$ | $0.54\pm 0.09$ | $1.01\pm 0.11$
PD-based | $0.34\pm 0.02$ | $0.56\pm 0.05$ | $0.53\pm 0.05$ | $0.99\pm 0.12$
$h_{a}$ only | $0.37\pm 0.02$ | $0.55\pm 0.07$

###  5.4 Real regression problems

We now apply all methods on two real-world static datasets. As algebraic term
$h_{k}$, we chose to use a linear prior on $x_{k}$, where $x_{k}$ is the
feature with the highest importance score in a RF model trained on the full
dataset (assumed to be inaccessible at training time). As there is no
guarantee that any of our hypotheses are met (and in particular H1), we do not
measure the distance $d_{k}$ in (2) as it is deemed irrelevant.

We consider two settings for each dataset, inspired by [33]. In the first
setting (INT), the training and test sets are sampled from the same
distribution $p(\mathbf{x},y)$ whereas the second one (EXT) evaluates
extrapolation performance for samples with unseen target values. If the linear
prior is a reasonable assumption or, at the very least, the target increases
or decreases monotonically with $x_{k}$, then we can expect hybrid methods to
yield better results in the latter setting. In the INT setting for each
dataset, we randomly select 100 samples for the learning set, 100 samples for
the validation set and keep the rest as test set. For the EXT setting, we
select the samples (one fourth of the dataset) with the lowest output values
as test set. From the remaining samples, we randomly select 100 samples for
the learning set and 100 samples for the validation set. For both INT and EXT
settings, performance metrics are averaged over 10 different splits. We
standardize both input and output variables.

The features for both datasets are described in Table 5. The Combined Cycle
Power Plant (CCPP) dataset [23] collects 9,568 measurements of net hourly
electrical energy output for a combined cycle power plant, along with four
hourly average input variables. The Concrete Compressive Strength (CCS)
dataset [29] is composed of 1,030 samples relating amounts of concrete
components with the resulting compressive strength. As done in [33], we
introduce a new feature corresponding to the cement-to-water ratio.

Table 5: Variables used for the real-world datasets. For each dataset, the variable indicated in bold type is the one used in the linear prior ($x_{k}$). Dataset | Name | Description
---|---|---
CCPP | T | Ambient temperature [°C]
| AP | Ambient pressure [mbar]
| RH | Relative humidity [-]
| V | Exhaust vacuum [cmHg]
CCS | Cement | Amount of cement in the mixture $\operatorname{[kg/m^{3}]}$
| Blast Furnace Slag | Amount of blast furnace slag in the mixture $\operatorname{[kg/m^{3}]}$
| Fly Ash | Amount of fly ash in the mixture $\operatorname{[kg/m^{3}]}$
| Water | Amount of water in the mixture $\operatorname{[kg/m^{3}]}$
| Superplasticizer | Amount of superplasticizer in the mixture $\operatorname{[kg/m^{3}]}$
| Coarse Aggregate | Amount of coarse aggregate in the mixture $\operatorname{[kg/m^{3}]}$
| Fine Aggregate | Amount of fine aggregate in the mixture $\operatorname{[kg/m^{3}]}$
| Age | Day (1-365)
| Cement/Water | Cement to water ratio [-]
Table 6: Results for the real-world datasets. We report the mean and standard deviation of $\hat{d}$ over the test sets. |  | CCPP | CCS
---|---|---|---
| Method | INT | EXT | INT | EXT
| Sequential | $0.07\pm 0.01$ | $0.23\pm 0.13$ | $0.25\pm 0.04$ | $0.74\pm 0.24$
MLP | Alternating | $0.07\pm 0.01$ | $0.24\pm 0.10$ | $0.24\pm 0.03$ | $0.89\pm 0.36$
| PD-based | $0.07\pm 0.01$ | $0.07\pm 0.01$ | $0.25\pm 0.03$ | $0.38\pm 0.04$
|  $h_{a}$ only | $0.07\pm 0.01$ | $0.36\pm 0.21$ | $0.25\pm 0.05$ | $0.83\pm 0.30$
| Sequential | $0.08\pm 0.01$ | $0.35\pm 0.07$ | $0.21\pm 0.02$ | $0.91\pm 0.24$
GB | Alternating | $0.08\pm 0.01$ | $0.41\pm 0.17$ | $0.20\pm 0.03$ | $0.85\pm 0.28$
| PD-based | $0.08\pm 0.02$ | $0.10\pm 0.01$ | $0.22\pm 0.02$ | $0.31\pm 0.06$
|  $h_{a}$ only | $0.10\pm 0.01$ | $0.95\pm 0.18$ | $0.25\pm 0.02$ | $1.54\pm 0.28$
| Sequential | $0.07\pm 0.01$ | $0.28\pm 0.07$ | $0.20\pm 0.02$ | $1.06\pm 0.31$
RF | Alternating | $0.07\pm 0.01$ | $0.32\pm 0.12$ | $0.20\pm 0.02$ | $1.00\pm 0.32$
| PD-based | $0.07\pm 0.01$ | $0.08\pm 0.01$ | $0.20\pm 0.02$ | $0.30\pm 0.08$
|  $h_{a}$ only | $0.08\pm 0.01$ | $0.99\pm 0.19$ | $0.25\pm 0.02$ | $1.89\pm 0.39$

From Table 6, it seems that introducing a linear prior does not yield any
benefit in the interpolation setting, as all models perform equivalently good,
which suggests that either the prior is not adequate or that H1 is not
verified. In the extrapolation scenario, we can however observe that the
linear prior allows to mitigate the impact of moving out of the distribution
compared to data-driven models. Indeed, compared to the INT setting, the
performance of all purely data-driven methods degrades in the EXT scenario,
especially for GB and RF as their output predictions are bounded by the
minimum target value observed in the training set. PD-based hybrid methods
consistently outperform other hybrid approaches and are only slightly
impacted, while sequential and alternating methods attain similar results.

##  6 Conclusion

We study several hybrid methods on supervised regression problems modeled in
an additive way, using neural network models and tree-based approaches. We
empirically show that trends observed for neural networks also apply for the
non-parametric tree-based approaches, in terms of predictive performance as in
the estimation of the algebraic known function. We introduce claims related to
the convergence of these hybrid approaches, under mild hypotheses, and verify
their soundness on illustrative experiments. We present a new hybrid approach
leveraging partial dependence and show its competitivity against sequential
and alternate optimization schemes on both synthetic and real-world problems.
We highlight its benefits in estimating the parametric prior and show that it
alleviates both the risk of the ML term to dominate the known term and the
need for assuming independent input features sets. As future work, we plan to
investigate further the theoretical properties of the PD-based approach and
extend it to dynamical problems.

## References

  * [1] Ayed, I., de Bézenac, E., Pajot, A., Brajard, J., Gallinari, P.: Learning dynamical systems from partial observations. Machine Learning and the Physical Sciences: Workshop at the 33rd Conference on Neural Information Processing Systems (NeurIPS) (2019)
  * [2] Breiman, L.: Random forests. Machine learning 45(1), 5–32 (2001)
  * [3] Chen, J., Yu, K., Hsing, A., Therneau, T.M.: A partially linear tree-based regression model for assessing complex joint gene–gene and gene–environment effects. Genetic Epidemiology: The Official Publication of the International Genetic Epidemiology Society 31(3), 238–251 (2007)
  * [4] Chen, T., Guestrin, C.: XGBoost: A scalable tree boosting system. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 785–794. KDD ’16, ACM, New York, NY, USA (2016)
  * [5] Daw, A., Karpatne, A., Watkins, W., Read, J., Kumar, V.: Physics-guided neural networks (pgnn): An application in lake temperature modeling. arXiv preprint arXiv:1710.11431 (2017)
  * [6] De Bézenac, E., Pajot, A., Gallinari, P.: Deep learning for physical processes: Incorporating prior scientific knowledge. Journal of Statistical Mechanics: Theory and Experiment 2019(12), 124009 (2019)
  * [7] Donà, J., Déchelle, M., Levy, M., Gallinari, P.: Constrained physical-statistics models for dynamical system identification and prediction. In: ICLR 2022-The Tenth International Conference on Learning Representations (2022)
  * [8] Friedman, J.H.: Greedy function approximation: a gradient boosting machine. Annals of statistics pp. 1189–1232 (2001)
  * [9] Friedman, J.H.: Stochastic gradient boosting. Computational statistics & data analysis 38(4), 367–378 (2002)
  * [10] Friedman, J.H., Grosse, E., Stuetzle, W.: Multidimensional additive spline approximation. SIAM Journal on Scientific and Statistical Computing 4(2), 291–301 (1983)
  * [11] Greene, W.H.: Econometric analysis. Pearson Education India (2003)
  * [12] Grinsztajn, L., Oyallon, E., Varoquaux, G.: Why do tree-based models still outperform deep learning on typical tabular data? In: Thirty-sixth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (2022)
  * [13] Hu, G., Mao, Z., He, D., Yang, F.: Hybrid modeling for the prediction of leaching rate in leaching process based on negative correlation learning bagging ensemble algorithm. Computers & chemical engineering 35(12), 2611–2617 (2011)
  * [14] Johansen, T.A., Foss, B.A.: Representing and learning unmodeled dynamics with neural network memories. In: 1992 American Control Conference. pp. 3037–3043. IEEE (1992)
  * [15] Kahrs, O., Marquardt, W.: Incremental identification of hybrid process models. Computers & Chemical Engineering 32(4-5), 694–705 (2008)
  * [16] Karniadakis, G.E., Kevrekidis, I.G., Lu, L., Perdikaris, P., Wang, S., Yang, L.: Physics-informed machine learning. Nature Reviews Physics 3(6), 422–440 (2021)
  * [17] Kramer, M.A., Thompson, M.L., Bhagat, P.M.: Embedding theoretical models in neural networks. In: 1992 American Control Conference. pp. 475–479. IEEE (1992)
  * [18] Mehta, V., Char, I., Neiswanger, W., Chung, Y., Nelson, A., Boyer, M., Kolemen, E., Schneider, J.: Neural dynamical systems: Balancing structure and flexibility in physical prediction. In: 2021 60th IEEE Conference on Decision and Control (CDC). pp. 3735–3742. IEEE (2021)
  * [19] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., Chintala, S.: Pytorch: An imperative style, high-performance deep learning library. In: Advances in Neural Information Processing Systems 32, pp. 8024–8035. Curran Associates, Inc. (2019)
  * [20] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E.: Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825–2830 (2011)
  * [21] Qian, Z., Zame, W., Fleuren, L., Elbers, P., van der Schaar, M.: Integrating expert odes into neural odes: Pharmacology and disease progression. Advances in Neural Information Processing Systems 34, 11364–11383 (2021)
  * [22] Takeishi, N., Kalousis, A.: Physics-integrated variational autoencoders for robust and interpretable generative modeling. Advances in Neural Information Processing Systems 34, 14809–14821 (2021)
  * [23] Tfekci, P., Kaya, H.: Combined Cycle Power Plant. UCI Machine Learning Repository (2014), DOI: https://doi.org/10.24432/C5002N
  * [24] Thompson, M.L., Kramer, M.A.: Modeling chemical processes using prior knowledge and neural networks. AIChE Journal 40(8), 1328–1340 (1994)
  * [25] Von Rueden, L., Mayer, S., Beckh, K., Georgiev, B., Giesselbach, S., Heese, R., Kirsch, B., Pfrommer, J., Pick, A., Ramamurthy, R., et al.: Informed machine learning–a taxonomy and survey of integrating prior knowledge into learning systems. IEEE Transactions on Knowledge and Data Engineering 35(1), 614–633 (2021)
  * [26] Wehenkel, A., Behrmann, J., Hsu, H., Sapiro, G., Louppe, G., Jacobsen, J.H.: Improving generalization with physical equations. Machine Learning and the Physical Sciences: Workshop at the 36th Conference on Neural Information Processing Systems (NeurIPS) (2022)
  * [27] Willard, J., Jia, X., Xu, S., Steinbach, M., Kumar, V.: Integrating scientific knowledge with machine learning for engineering and environmental systems. ACM Computing Surveys 55(4), 1–37 (2022)
  * [28] Wouwer, A.V., Renotte, C., Bogaerts, P.: Biological reaction modeling using radial basis function networks. Computers & chemical engineering 28(11), 2157–2164 (2004)
  * [29] Yeh, I.C.: Concrete Compressive Strength. UCI Machine Learning Repository (2007), DOI: https://doi.org/10.24432/C5PK67
  * [30] Yin, Y., Ayed, I., de Bézenac, E., Baskiotis, N., Gallinari, P.: Leads: Learning dynamical systems that generalize across environments. Advances in Neural Information Processing Systems 34, 7561–7573 (2021)
  * [31] Yin, Y., Le Guen, V., Dona, J., de Bézenac, E., Ayed, I., Thome, N., Gallinari, P.: Augmenting physical models with deep networks for complex dynamics forecasting. Journal of Statistical Mechanics: Theory and Experiment 2021(12), 124012 (2021)
  * [32] Yu, K., Wheeler, W., Li, Q., Bergen, A.W., Caporaso, N., Chatterjee, N., Chen, J.: A partially linear tree-based regression model for multivariate outcomes. Biometrics 66(1), 89–96 (2010)
  * [33] Zhang, H., Nettleton, D., Zhu, Z.: Regression-enhanced random forests. arXiv preprint arXiv:1904.10416 (2019)

##  Appendix 0.A Optimal model under H2 and H3

Let us recall the regression problem, where $y\in\mathbb{R}$ can be decomposed
into the addition of two independent terms:

$y=f_{k}(\mathbf{x}_{k})+f_{a}^{r}(\mathbf{x}_{a})+\varepsilon,\quad\varepsilon\sim\mathcal{N}(0,\sigma^{2}),\quad\mathbf{x}_{k}\cup\mathbf{x}_{a}=\mathbf{x},\mathbf{x}_{k}\cap\mathbf{x}_{a}=\emptyset,\mathbf{x}_{k}\perp\\!\\!\\!\perp\mathbf{x}_{a}.$

For clarity, let us denote by $\mathbb{E}_{\mathbf{x}}$ the subsequent
expectations over the input space $\mathbb{E}_{\mathbf{x}\sim
p(\mathbf{x})}\left\\{\cdot\right\\}$. We have:

$\displaystyle\hat{f}_{k}^{*}$
$\displaystyle=\arg\min_{\hat{f}_{k}\in\hat{\mathcal{F}}_{k}}d(\hat{f}_{k},y)=\arg\min_{\hat{f}_{k}\in\hat{\mathcal{F}}_{k}}\mathbb{E}_{\mathbf{x},\varepsilon}\left\\{\left(\hat{f}_{k}(\mathbf{x}_{k})-f_{k}(\mathbf{x}_{k})-f_{a}^{r}(\mathbf{x}_{a})-\varepsilon\right)^{2}\right\\}$
$\displaystyle=\arg\min_{\hat{f}_{k}\in\hat{\mathcal{F}}_{k}}\mathbb{E}_{\mathbf{x}_{k}}\left\\{(\hat{f}_{k}(\mathbf{x}_{k})-f_{k}(\mathbf{x}_{k}))^{2}\right\\}+\mathbb{E}_{\mathbf{x}_{a},\varepsilon}\left\\{(f_{a}^{r}(\mathbf{x}_{a})+\varepsilon)^{2}\right\\}$
$\displaystyle-\mathbb{E}_{\mathbf{x},\varepsilon}\left\\{2(\hat{f}_{k}(\mathbf{x}_{k})-f_{k}(\mathbf{x}_{k}))(f_{a}^{r}(\mathbf{x}_{a})+\varepsilon)\right\\}.$

The second term is independent w.r.t. $\hat{f}_{k}$ and thus has no impact on
the minimization. Moreover, since
$\mathbf{x}_{k}\perp\\!\\!\\!\perp\mathbf{x}_{a}$, the last term writes as the
product of two expectations, one of which is constant w.r.t. $\hat{f}_{k}$. We
thus have:

$\hat{f}_{k}^{*}=\arg\min_{\hat{f}_{k}\in\hat{\mathcal{F}}_{k}}\mathbb{E}_{\mathbf{x}_{k}}\left\\{(\hat{f}_{k}(\mathbf{x}_{k})-f_{k}(\mathbf{x}_{k}))^{2}\right\\}-2C\mathbb{E}_{\mathbf{x}_{k}}\left\\{(\hat{f}_{k}(\mathbf{x}_{k})-f_{k}(\mathbf{x}_{k}))\right\\},$
(10)

with
$C=\mathbb{E}_{\mathbf{x}_{a},\varepsilon}\left\\{(f_{a}^{r}(\mathbf{x}_{a})+\varepsilon)\right\\}=\mathbb{E}_{\mathbf{x}_{a}}\left\\{f_{a}^{r}(\mathbf{x}_{a})\right\\}$.
Cancelling the derivative of (10) w.r.t. $\hat{f}_{k}$, we obtain the optimal
model $\hat{f}_{k}^{*}(\mathbf{x}_{k})=f_{k}(\mathbf{x}_{k})+C$, for every
$\mathbf{x}_{k}\in\mathcal{X}_{k}$.

##  Appendix 0.B Model architectures

Model hyperparameters are reported in Table 7. We used PyTorch [19] for MLP,
scikit-learn [20] for RF and xgboost [4] for GB. Unspecified parameters keep
their default values. Learning rates for training $h_{k}$ and MLP are set to
$0.005$. $\operatorname{H}$ is the number of hidden layers in MLP and
$\operatorname{W}$ the number of neurons per hidden layer. $\operatorname{T}$
is the number of trees in GB and RF, $\operatorname{d}$ the maximum tree depth
and $\operatorname{mss}$ the minimum number of samples required to split an
internal tree node.

Table 7: Model hyperparameters, for each experiment. Problem | MLP | GB | RF
---|---|---|---
(Correlated) Friedman | $\operatorname{H}=2,\operatorname{W}=15$ | $\operatorname{T}=700,\operatorname{d}=2$ | $\operatorname{T}=500,\operatorname{mss}=5$
Linear & Overlap | $\operatorname{H}=2,\operatorname{W}=10$ | $\operatorname{T}=400,\operatorname{d}=2$ | $\operatorname{T}=500,\operatorname{mss}=5$
Real-world problems | $\operatorname{H}=2,\operatorname{W}=30$ | $\operatorname{T}=300,\operatorname{d}=2$ | $\operatorname{T}=200,\operatorname{mss}=5$


# Spin-orbit torques and magnetization switching in Gd/Fe multilayers
generated by current injection in NiCu alloys

Federica Nasr  Federico Binda  Charles-Henri Lambert  Giacomo Sala  Paul Noël
Pietro Gambardella  Department of Materials, ETH Zurich, CH-8093 Zurich,
Switzerland  Authors to whom correspondence should be addressed:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
(29 December 2023)

###### Abstract

Light transition metals have recently emerged as a sustainable material class
for efficient spin-charge interconversion. We report measurements of current-
induced spin-orbit torques generated by Ni1-xCux alloys in perpendicularly-
magnetized ferrimagnetic Gd/Fe multilayers. We show that the spin-orbit torque
efficiency of Ni1-xCux increases with the Ni/Cu atomic ratio, reaching values
comparable to those of Pt for Ni55Cu45. Furthermore, we demonstrate
magnetization switching of a 20-nm-thick Gd/Fe multilayer with a threshold
current that decreases with increasing Ni concentration, similar to the spin-
orbit torque efficiency. Our findings show that Ni1-xCu${}_{x}-$based magnetic
heterostructures allow for efficient control of the magnetization by electric
currents.

††preprint: AIP/123-QED

Materials that allow for converting charge currents into spin currents are
used to generate spin-orbit torques (SOTs) in different classes of magnetic
systems and devices.Manchon  _et al._ (2019) SOTs enable the energy-efficient
manipulation of the magnetization in ferromagnets,Miron  _et al._ (2011); Liu
_et al._ (2012); Grimaldi  _et al._ (2020) ferrimagnets,Finley and Liu (2016);
Caretta  _et al._ (2018); Pham  _et al._ (2018); Sala and Gambardella (2022a)
and antiferromagnets,Olejník  _et al._ (2018); Tsai  _et al._ (2020);
Krishnaswamy  _et al._ (2022) with applications in magnetic random access
memories,Krizakova  _et al._ (2022) racetrack memories,Bläsing  _et al._
(2020) logic devices,Luo  _et al._ (2020) and spin-torque nano-
oscillators.Demidov  _et al._ (2017); Zahedinejad  _et al._ (2022) SOTs have
been extensively studied in heavy metal/ferromagnet bilayers,Miron  _et al._
(2011); Liu  _et al._ (2012); Pai  _et al._ (2012); Garello  _et al._ (2013);
Kim  _et al._ (2013); Manchon  _et al._ (2019); Ishikuro  _et al._ (2019);
Miron  _et al._ (2010) where the charge-to-spin conversion originates from the
interplay of bulk and interfacial spin-orbit coupling, as exemplified by the
spin Hall Sinova  _et al._ (2015) and Rashba-Edelstein effects.Bihlmayer  _et
al._ (2022) The two figures of merit describing the charge-to-spin conversion
efficiency are the spin Hall angle and the spin Hall conductivity. The former
is defined as the ratio of the transverse spin current generated in the spin
injector to the applied charge current, whereas the latter is the ratio of the
spin current to the external electric field driving the charge current.Manchon
_et al._ (2019); Nguyen, Ralph, and Buhrman (2016)

A large spin Hall angle ($\geq 0.1$) is usually observed in the _5d_
transition metals, such as Ta, W, and Pt, owing to their strong atomic spin-
orbit coupling.Sinova  _et al._ (2015) Whether light metals can also show
large charge-to-spin conversion is a question of both fundamental and
practical impact on developing new and low-cost spin injector materials. In
this context, _3d_ transition metals and their alloys have recently attracted
interest because they present appreciable spin Hall angles despite their weak
spin-orbit coupling.Miao  _et al._ (2013); Du  _et al._ (2014); Qu, Huang, and
Chien (2015); An  _et al._ (2016); Wang  _et al._ (2017); h. C. Baek  _et al._
(2018); Hibino  _et al._ (2020) The interest in light metals has further
surged following the theoretical prediction of a large orbital Hall
conductivity in the nonmagnetic _3d_ elements,Kontani  _et al._ (2009); Jo,
Go, and Lee (2018) which can lead to the generation of strong SOTs in adjacent
ferromagnetic layers.Go  _et al._ (2020); Lee  _et al._ (2021a, b); Sala and
Gambardella (2022b)

For transition metals, besides the spin-orbit coupling, the _d_ -orbital
occupation plays an important role in determining the magnitude and sign of
the spin Hall angle.Sinova  _et al._ (2015); Kontani  _et al._ (2009) Ni is
particularly interesting in this respect because it has the same valency and a
similar band structure as Pt. The calculated spin Hall conductivity of Ni is
$\sim 1.2\times
10^{5}~{}$\text{\,}\mathrm{\SIUnitSymbolOhm}^{-1}\mathrm{m}^{-1}$$, about half
of Pt, whereas the orbital Hall conductivities of Ni and Pt are very
similar.Jo, Go, and Lee (2018) Extrinsic spin scattering effects due to
alloying of different metallic elements can further enhance the charge-spin
interconversion efficiency,Zhu, Ralph, and Buhrman (2021); Hu and Pai (2020);
Chi  _et al._ (2021); Wang  _et al._ (2022) as shown, e.g., for Cu1-xBix and
Au1-xPtx.Niimi  _et al._ (2012); Zhu, Ralph, and Buhrman (2018) In light of
these observations, experimental work has focused on NixCu1-x alloys, in which
the introduction of Cu reduces the Curie temperature relative to pure Ni,
making the alloy paramagnetic at room temperature for $x\leq 70\%$.Bozorth
(1993); Ahern, Martin, and Sucksmith (1958) Despite the alloying, the band
structure of NixCu1-x remains similar to that of paramagnetic Ni, suggesting a
significant intrinsic spin Hall effect regardless of chemical disorder.Keller
_et al._ (2019) Recent experiments have reported large effective spin Hall
angles comparable to 5$d$ metals in Permalloy/Ni60Cu40, CoFeB/Cu/NixCu1-x, and
YIG/Ni80Cu20 bilayers.Keller  _et al._ (2019); Varotto  _et al._ (2020); Wu
_et al._ (2022) Additionally, spin-to-charge conversion in CoFeB/Ni70Cu30 was
shown to induce $\text{\,}\mathrm{THz}$ emission upon fs laser irradiation
with intensity up to half that of CoFeB/Pt.Cheng  _et al._ (2022)

These studies show that NixCu1-x alloys constitute an interesting material
system for spin-charge interconversion and the realization of spin current
detectors. Here, we investigate NixCu1-x alloys as spin injectors for
generating SOTs in perpendicularly magnetized NixCu1-x/Ti/[Gd/Fe]30
heterostructures. We select three different compositions ($x=40,\,50,\,55\%$)
for which NixCu1-x is paramagnetic at room temperature and report composition-
dependent SOT efficiencies, reaching values of up to 0.06 for $x=55\%$, which
corresponds to a spin Hall conductivity of $(1.04\pm 0.07)\times
10^{5}~{}$\text{\,}\mathrm{\SIUnitSymbolOhm}^{-1}\mathrm{m}^{-1}$$ without
taking into account corrections due to spin memory loss and interface
transparency. In addition, we demonstrate the efficient magnetization
switching capability of these systems, including a fourth alloy’s composition
$x=70\%$ at the ferromagnetic phase transition. Our results show that NixCu1-x
alloys are efficient spin current generators that can be used to achieve
electrical control of the magnetization in spintronic devices.

We fabricated three different sets of samples, all grown on SiO2 substrates by
magnetron sputtering: (i) the full heterostructures
NixCu1-x($6\text{\,}\mathrm{nm}$)/Ti($1\text{\,}\mathrm{nm}$)/FiM/Si3N4
($x=40,\,50,\,55,\,70\%$), referred to as Ni$x$, where FiM stands for the
ferrimagnetic multilayer
[Gd($0.4\text{\,}\mathrm{nm}$)/Fe($0.27\text{\,}\mathrm{nm}$)]30, (ii) the
reference layers NixCu1-x($6\text{\,}\mathrm{nm}$)/Si3N4, referred to as
Ni$x-r$, and (iii) the single FiM layer Ti($1\text{\,}\mathrm{nm}$)/FiM/Si3N4.
The NixCu1-x layers are amorphous, as revealed by X-ray diffraction. The
NixCu1-x layers in Ni40, Ni50, and Ni55 are paramagnetic at room temperature,
which avoids spurious magnetotransport signals in the electrical measurements
of SOTs. The Ti layer in all the heterostructures serves as a spacer to
suppress possible magnetic proximity effects between NixCu1-x and the FiM. A
sketch of the Ni$x$ samples is shown in Fig. 1(a). The samples were patterned
in $7.5\text{\,}$-$\text{\,}\mathrm{\SIUnitSymbolMicro m}$-wide and
$11\text{\,}$-$\text{\,}\mathrm{\SIUnitSymbolMicro m}$-long Hall bar devices
using photolithography and lift-off, as shown in Fig. 1(b).

Figure 1:  (a) Schematic of the Ni$x$ heterostructure. (b) Optical image of a
Hall bar device with superimposed electrical connections and coordinate
system. (c) Magnetic moment per unit area as a function of out-of-plane
magnetic field for Ni55 (green dots), its reference Ni55$-r$ (open circles),
and the single FiM sample (purple hexagons). (d) Anomalous Hall resistance as
a function of out-of-plane magnetic field for Ni55 (green dots) and the single
FiM (purple hexagons). The current density in the FiM layer was
$j_{\textrm{FiM}}\sim$3.5\text{\times}{10}^{10}\text{\,}\mathrm{A}\mathrm{/}\mathrm{m}^{2}$$
in both measurements.

Figure 1(c) shows $M_{\textrm{s}}\,t$, the magnetic moment per unit area, as a
function of out-of-plane magnetic field of Ni$55$, its reference Ni$55-r$, and
the single FiM sample. Here, $M_{\textrm{s}}$ represents the saturation
magnetization, $t$ the thickness of the magnetic layer, and
$M_{\textrm{s}}\,t$ is the total magnetic moment of the as-grown layers
measured at room temperature by a superconducting quantum interference device
divided by the sample area. The comparison between the different curves shows
that Ni$55-r$ develops a sizable magnetization in a magnetic field, which adds
to that of the FiM in the Ni$55$ sample. Because we consider the effect of the
SOTs on the sole FiM, we estimate an effective magnetic moment per unit area
for each Ni$x$ sample as the difference between the moments of the full
heterostructures Ni$x$ and their references Ni$x-r$, namely
$M_{\textrm{s}}\,t=(M_{\textrm{s}}\,t)_{\textrm{Ni}x}-(M_{\textrm{s}}\,t)_{{\textrm{Ni}x-r}}$.
The magnetotransport measurements were performed by injecting an alternate
current of frequency $\omega/(2\pi)=$10\text{\,}\mathrm{Hz}$$ and density $j$
and computing the respective Hall ($R^{\textrm{H}}_{1,2\omega}$) and
longitudinal ($R^{\textrm{L}}_{1,2\omega}$) components of the first and second
harmonic resistances from the measured harmonic voltages.Garello  _et al._
(2013); Avci  _et al._ (2014) The first harmonic Hall resistance accounts for
the contributions from the anomalous Hall effect (AHE) and planar Hall effect
(PHE) and is given by

$R^{\textrm{H}}_{1\omega}=R_{\textrm{AHE}}\cos\vartheta+R_{\textrm{PHE}}\sin^{2}\vartheta\sin
2\varphi,$ (1)

where $R_{\textrm{AHE}}$ and $R_{\textrm{PHE}}$ are the anomalous and planar
Hall resistance coefficients, and $\vartheta$ and $\varphi$ are the polar and
azimuthal angles between the magnetization and the $z$\- and $x$-axis,
respectively. Figure 1(d) shows $R^{\textrm{H}}_{1\omega}$ as a function of
out-of-plane field for Ni55 and the single FiM sample. The signal is dominated
by the AHE of the FiM in both samples; the remanence and polarity of the
curves indicate that the FiM has perpendicular magnetic anisotropy and its
magnetization is Fe-like. The amplitude of the AHE in Ni55 is reduced compared
to the single FiM sample because part of the current flows in the paramagnetic
NixCu1-x layer and Ni has opposite sign of $R_{\textrm{AHE}}$ compared to Fe.
Omori  _et al._ (2019) Moreover, the magnetic properties and AHE of FiM layers
are known to depend on the adjacent layers, which are different in the two
samples.Haltz  _et al._ (2018); Zhang  _et al._ (2020); Dang  _et al._ (2020);
Bello  _et al._ (2022)

The SOTs were measured using the harmonic Hall voltage method, which relies on
$R^{\textrm{H}}_{2\omega}$ as a measure of the current-induced oscillations of
the magnetization.Garello  _et al._ (2013); Hayashi  _et al._ (2014) The SOTs
can be decomposed into two components orthogonal to the unitary magnetization
m, to define a damping-like (DL) torque,
$\textbf{T}_{\text{DL}}=\tau_{\text{DL}}\,\textbf{m}\times(\bm{\zeta}\times\textbf{m})$,
and a field-like (FL) torque,
$\textbf{T}_{\text{FL}}=\tau_{\text{FL}}\,\textbf{m}\times\bm{\zeta}$. Here,
the torques are normalized by the magnetization, $\bm{\zeta}$ is the unit
vector representing the spin polarization of the injected spin current, and
$\tau_{\text{DL}}$ and $\tau_{\text{FL}}$ are the amplitudes of the effective
magnetic fields describing the action of the SOTs on the magnetization. We
estimated $\tau_{\text{DL}}$ and $\tau_{\text{FL+Oe}}$, with Oe denoting the
Oersted field, using the following equations valid under the assumption
$\vartheta\to 0$ Hayashi  _et al._ (2014)

$\displaystyle\tau_{\textrm{DL}}=2\frac{B_{x}\pm 2\chi B_{y}}{1-4\chi^{2}},$
(2) $\displaystyle\tau_{\textrm{FL+Oe}}=2\frac{B_{y}\pm 2\chi
B_{x}}{1-4\chi^{2}},$ (3)

where $\chi=R_{\textrm{PHE}}/R_{\textrm{AHE}}$ ($\sim 10^{-2}$ in our samples)
and $B_{x,y}\equiv\bigl{(}\frac{\partial R^{\textrm{H}}_{2\omega}}{\partial
H}/\frac{\partial^{2}R^{\textrm{H}}_{1\omega}}{\partial^{2}H}\bigr{)}\big{|}_{\hat{m}\parallel{x,y\,}}$.
As an illustrative case, Figs. 2(a) and (b) show $R^{\textrm{H}}_{1,2\omega}$
of Ni40 as a function of a nearly in-plane external magnetic field oriented
along the $x$\- and $y$-axes measured at
$j=$7.4\text{\times}{10}^{10}\text{\,}\mathrm{A}\mathrm{/}\mathrm{m}^{2}$$
($j$NiCu
$=$1.9\text{\times}{10}^{11}\text{\,}\mathrm{A}\mathrm{/}\mathrm{m}^{2}$$) and
the fits (solid lines) used to estimate the SOT effective fields. Figure 2(c)
shows $\tau_{\textrm{DL}}$ and $\tau_{\textrm{FL+Oe}}$ vs $j_{\text{NiCu}}$
for the different Ni$x$ samples. Note that the total FL+Oe field is similar to
the expected Oersted field from Biot-Savart’s law, evidencing the small
contribution of the FL torque in our samples.

Figure 2:  (a) and (b) $R^{\textrm{H}}_{1\omega}$ (left axis) and
$R^{\textrm{H}}_{2\omega}$ (right axis) of Ni40 as a function of external
magnetic field applied at an angle $\vartheta_{\textrm{B}}=87^{\circ}$ along
(a) the $x$-axis (left plot) and (b) the $y$-axis (right plot). The solid
lines are fits used to estimate
$\partial^{2}R^{\textrm{H}}_{1\omega}/\partial^{2}H$ and $\partial
R^{\textrm{H}}_{2\omega}/\partial H$ in Eq. 3. (c) DL (black) and FL (red)
effective fields of Ni40, Ni50, and Ni55 as a function of the current density
in the NixCu1-x layer. The solid lines are linear fits with intercept forced
to zero. The error bars, which include the fit errors in (a) and (b) and the
uncertainties in $R_{\textrm{AHE}}$ and $R_{\textrm{PHE}}$, are smaller than
the symbols. (d) DL (black) and FL (red) SOT efficiencies of Ni$x$. The left
panel shows the SOT efficiencies with respect to the total current density
$\upxi^{j}$ (empty columns) and to the current density in the sole NixCu1-x
layer, $\upxi^{j_{\textrm{{NiCu}}}}$ (full columns); the right panel shows the
SOT efficiencies with respect to the electric field $\upxi^{E}$. The errors
are calculated from the uncertainty of the fits in (c) and
$M_{\textrm{s}}\,t$.

The SOT efficiencies with respect to the applied current density ($j$) and
electric field driving the current ($E$) are given byManchon  _et al._ (2019);
Nguyen, Ralph, and Buhrman (2016)

$\displaystyle\upxi^{j}_{\textrm{DL,FL}}$
$\displaystyle=\frac{2e}{\hbar}M_{\textrm{s}}t\frac{\tau_{\textrm{DL,FL}}}{j},$
(4) $\displaystyle\upxi^{E}_{\textrm{DL,FL}}$
$\displaystyle=\frac{2e}{\hbar}M_{\textrm{s}}t\frac{\tau_{\textrm{DL,FL}}}{E},$
(5)

where $e$ is the elementary charge and $\hbar$ is the reduced Planck constant.
$\upxi^{j}$ ($\upxi^{E}$) represents the ratio of the effective spin current
absorbed by the magnetic layer relative to the applied charge current
(electric field) and can be thus regarded as an effective spin Hall angle
(effective spin Hall conductivity) for a specific combination of a spin
injector and magnetic layer. Because distinguishing between SOTs generated by
spin and orbital currents is not straightforward,Sala and Gambardella (2022b)
$\upxi^{j}$ and $\upxi^{E}$ include also the effect of orbital currents. The
efficiencies calculated by Eqs. 4 and 5 constitute a lower bound for
quantifying the charge-to-spin conversion of a spin injector since they do not
account for the spin transparency of the interface between the spin injector
and the magnetic material.

Figure 2(d) presents the SOT efficiencies calculated by averaging over two
devices for each Ni$x$ sample. The sign of the DL torque is consistent with
that of Pt/Co, considering the same stacking order between the spin injector
and the magnetic layer. Both $\upxi_{\textrm{DL}}$ and
$\upxi_{\textrm{FL+Oe}}$ increase monotonically with Ni content and reach a
maximum at $x=55\%$ in the measured composition range:
$\upxi^{j}_{\textrm{DL}}=0.16$ ($\upxi^{E}_{\textrm{DL}}=1.04\times
10^{5}~{}$\text{\,}\mathrm{\SIUnitSymbolOhm}^{-1}\mathrm{m}^{-1}$$) and
$\upxi^{j}_{\textrm{FL+Oe}}=0.08$ ($\upxi^{E}_{\textrm{FL+Oe}}=0.48\times
10^{5}~{}$\text{\,}\mathrm{\SIUnitSymbolOhm}^{-1}\mathrm{m}^{-1}$$) at
$x=55\%$ are comparable with the values reported for 5$d$ metal system such as
Pt/Co.Manchon  _et al._ (2019) In making such a comparison, one should note
that the effective spin Hall angle changes if it is calculated with respect to
the average current density ($\upxi^{j}$) or to the current density flowing in
the spin injector ($\upxi^{j_{\textrm{NiCu}}}$) estimated using a parallel
resistor model. In the latter case, we obtain
$\upxi^{j_{\textrm{NiCu}}}_{\textrm{DL}}=0.06$ and
$\upxi^{j_{\textrm{NiCu}}}_{\textrm{FL+Oe}}=0.03$, which are about half of the
values reported for Pt/Co. The resistivity and magnetization parameters used
to calculate the different $\upxi$ are reported in Table 1.

We now compare the SOT efficiencies of the Ni$x$ samples with the measurements
of spin-charge interconversion reported for NixCu1-x alloys in the literature.
Keller _et al._ Keller  _et al._ (2019) reported $\upxi^{E}_{\textrm{DL}}\sim
1\times 10^{5}~{}$\text{\,}\mathrm{\SIUnitSymbolOhm}^{-1}\mathrm{m}^{-1}$$ at
$x=60\%$, before considering any correction due to interfacial spin memory
loss. This efficiency, normalized by the reported electrical conductivity of
the NiCu layer of
$1.93\text{\times}{10}^{6}\text{\,}\mathrm{\SIUnitSymbolOhm}^{-1}\mathrm{m}^{-1}$,
corresponds to a spin Hall angle of 0.05. Varotto _et al._ Varotto  _et al._
(2020) reported a spin Hall angle of $0.04$ at $x=60\%$, in line with our
results. Because the spin diffusion length $\lambda_{s}$ decreases with
increasing Ni content Hsu  _et al._ (1996) and the measured values of
$\tau_{\text{DL}}$ and $\tau_{\text{FL}}$ depend on $\lambda_{s}$ if $t_{\rm
NiCu}\lesssim 2\lambda_{s}$, an alternative figure of merit is the product
$\upxi^{j_{\textrm{NiCu}}}\lambda_{s}$. Assuming
$\lambda_{s}\sim$4\text{\,}\mathrm{nm}$$ at $x=40\%$ and
$3\text{\,}\mathrm{nm}$ at $x=55\%$,Hsu  _et al._ (1996) we obtain
$\upxi_{\textrm{DL}}^{j_{\textrm{NiCu}}}\lambda_{s}=$0.12\text{\,}\mathrm{nm}$$
for Ni40 and $0.18\text{\,}\mathrm{nm}$ for Ni55. These values are similar to
$\upxi_{\textrm{DL}}\lambda_{s}=$0.1\text{\,}\mathrm{nm}$$ obtained in Ref.
Varotto  _et al._ , 2020. Overall, this comparison shows that NixCu1-x alloys
present large and comparable figures of merit for charge-to-spin and spin-to-
charge conversion.

Figure 3:  SOT-induced magnetization switching. (a) and (b) Anomalous Hall resistance as a function of current density in NixCu1-x. The switching loops are acquired using $10\text{\,}$-$\text{\,}\mathrm{ms}$-long current pulses in the presence of an in-plane magnetic field of $10\text{\,}\mathrm{mT}$ oriented along (a) $-x$ and (b) $+x$. (c) Anomalous Hall resistance as a function of current density for the single FiM sample. The loop is acquired using $10\text{\,}$-$\text{\,}\mathrm{ms}$-long current pulses and an in-plane magnetic field of $170\text{\,}\mathrm{mT}$ oriented along $+x$. (d) Ratio of the magnetic moment per unit area to the threshold switching current as a function of Ni content. Sample |  $M_{\textrm{s}}\,t$ ($$\text{\,}\mathrm{mA}$)$ |  $R_{\textrm{AHE}}$ $($\text{\,}\mathrm{\SIUnitSymbolOhm}$)$ |  $\rho$ $($\text{\,}\mathrm{\SIUnitSymbolMicro\SIUnitSymbolOhm}\text{\,}\mathrm{cm}$)$ |  $\rho$NiCu $($\text{\,}\mathrm{\SIUnitSymbolMicro\SIUnitSymbolOhm}\text{\,}\mathrm{cm}$)$ |  $\rho$FiM $($\text{\,}\mathrm{\SIUnitSymbolMicro\SIUnitSymbolOhm}\text{\,}\mathrm{cm}$)$ |  $R/R\textsubscript{NiCu}$ $(\%)$
---|---|---|---|---|---|---
Ni40 | 1.0 | 0.8 | 170 | 64 | 320 | 58
Ni50 | 1.1 | 0.7 | 181 | 60 | 423 | 67
Ni55 | 1.7 | 1.0 | 157 | 58 | 308 | 60
Ni70 | 2.9 | 0.5 | 129 | 49 | 242 | 58
FiM | 1.3 | 7.5 | 335 | - | 335 | -
Table 1: Summary of the magnetic and electrical properties of Ni$x$ and the
single FiM sample at room temperature. From left to right: $M_{\textrm{s}}\,t$
of the FiM, $R_{\rm AHE}$, resistivity of Ni$x$, Ni$x-r$, and FiM layer, and
current branching ratio in the NixCu1-x layer. For the Ni$x$ samples,
$\rho_{\textrm{FiM}}$ was estimated using a parallel resistor model.

Based on the large charge-to-spin conversion efficiency, we further
demonstrate SOT-induced magnetization switching of Ni$x$ for four different
alloy’s compositions ($x=40,\,50,\,55,\,70\%$). To do so, we injected
$10\text{\,}$-$\text{\,}\mathrm{ms}$-long current pulses to induce
magnetization switching and inspected the magnetization state by reading
$R^{H}_{1\omega}$ using an alternate current of frequency
$\omega/(2\pi)=$115\text{\,}\mathrm{Hz}$$ after each pulse. As required for
samples with perpendicular magnetic anisotropy, we apply an in-plane magnetic
field to break the symmetry of the DL torque and define the switching
polarity.Miron  _et al._ (2011) Figures 3(a) and (b) show the switching loops
of Ni$x$ acquired for negative and positive external magnetic fields of
amplitude $10\text{\,}\mathrm{mT}$ applied along the $x$-axis. For all the
samples, we observed full magnetization switching, as demonstrated by the
semi-amplitude of the switching loops, matching the saturation
$R_{\text{AHE}}$ values. As expected, the switching polarity changes upon
inversion of the in-plane magnetic field. In addition, the switching polarity
is the same as observed in Pt/Co,Miron  _et al._ (2011) consistently with the
sign of the SOTs. We further verified that switching occurs reliably using
$100\text{\,}$-$\text{\,}\mathrm{ns}$-long current pulses and that the
threshold switching current $I_{\textrm{th}}$ decreases monotonically with the
applied in-plane field. Measurements of the single FiM sample gave null
results, as shown in Fig. 3(c), excluding the possibility of switching induced
by self-torques in the FiM material.Céspedes-Berrocal  _et al._ (2021)

$I_{\textrm{th}}$ is generally influenced by the coercivity and magnetic
anisotropy of the FiM layer, which can vary from sample to sample. If such
variations are neglected, the ratio $M_{\text{s}}\,t/I_{\textrm{th}}$ can be
taken as a rough measure of the SOT switching efficiency. Figure 3(d) shows
that $M_{\text{s}}\,t/I_{\textrm{th}}$ measured in the different Ni$x$ samples
improves with increasing Ni content, in agreement with the trend of the SOT
efficiency reported in Fig. 2(d). The same trend is observed if
$I_{\textrm{th}}$ is taken as the total current (squares) or only the current
flowing in the NixCu1-x layer (dots). The improved switching efficiency at
$x=70\%$ is also in agreement with the higher charge-to-spin conversion
efficiency measured by ferromagnetic resonance.Keller  _et al._ (2019)

In conclusion, our results show that NixCu1-x alloys are very efficient
generators of SOTs. Current injection in NixCu1-x is further shown to induce
reliable switching of the magnetization of a 20-nm-thick [Gd/Fe] multilayer at
a current density of about $5\times 10^{11}$ A/m2. The SOT efficiency
increases with Ni content, reaching a maximum spin Hall conductivity of
$\upxi^{E}_{\textrm{DL}}=1.04\times
10^{5}~{}$\text{\,}\mathrm{\SIUnitSymbolOhm}^{-1}\mathrm{m}^{-1}$$ in Ni55,
which is comparable to that of heavy metal systems.

###### Acknowledgements.

This work was supported by the Swiss National Science Foundation (Grant No.
200020-200465). P.N. acknowledges the support of the ETH Zurich Postdoctoral
Fellowship Program 19-2 FEL-61. We thank Laura van Schie for her assistance in
performing switching measurements using
$100\text{\,}$-$\text{\,}\mathrm{ns}$-long current pulses.

The authors have no conflicts to disclose.

The data that support the findings of this study are openly available from the
ETH Research Collection at https://doi.org/10.3929/ethz-b-000633473.

## References

  * Manchon  _et al._ (2019) A. Manchon, J. Železný, I. M. Miron, T. Jungwirth, J. Sinova, A. Thiaville, K. Garello, and P. Gambardella, “Current-induced spin-orbit torques in ferromagnetic and antiferromagnetic systems,” Rev. Mod. Phys. 91(3), 035004 (2019).
  * Miron  _et al._ (2011) I. M. Miron, K. Garello, G. Gaudin, P.-J. Zermatten, M. V. Costache, S. Auffret, S. Bandiera, B. Rodmacq, A. Schuhl, and P. Gambardella, “Perpendicular switching of a single ferromagnetic layer induced by in-plane current injection,” Nature 476, 189–193 (2011).
  * Liu  _et al._ (2012) L. Liu, C.-F. Pai, Y. Li, H. Tseng, D. Ralph, and R. Buhrman, “Spin-torque switching with the giant spin Hall effect of tantalum,” Science 336, 555 (2012).
  * Grimaldi  _et al._ (2020) E. Grimaldi, V. Krizakova, G. Sala, F. Yasin, S. Couet, G. S. Kar, K. Garello, and P. Gambardella, “Single-shot dynamics of spin–orbit torque and spin transfer torque switching in three-terminal magnetic tunnel junctions,” Nat. Nanotechnol. 15, 111–117 (2020).
  * Finley and Liu (2016) J. Finley and L. Liu, “Spin-Orbit-Torque Efficiency in Compensated Ferrimagnetic Cobalt-Terbium Alloys,” Phys. Rev. Appl. 6, 054001 (2016).
  * Caretta  _et al._ (2018) L. Caretta, M. Mann, F. Büttner, K. Ueda, B. Pfau, C. M. Günther, P. Hessing, A. Churikova, C. Klose, M. Schneider, D. Engel, C. Marcus, D. Bono, K. Bagschik, S. Eisebitt, and G. S. D. Beach, “Fast current-driven domain walls and small skyrmions in a compensated ferrimagnet,” Nat. Nanotechnol. 13, 1154–1160 (2018).
  * Pham  _et al._ (2018) T. H. Pham, S.-G. Je, P. Vallobra, T. Fache, D. Lacour, G. Malinowski, M. C. Cyrille, G. Gaudin, O. Boulle, M. Hehn, J.-C. Rojas-Sánchez, and S. Mangin, “Thermal Contribution to the Spin-Orbit Torque in Metallic-Ferrimagnetic Systems,” Phys. Rev. Appl. 9, 064032 (2018).
  * Sala and Gambardella (2022a) G. Sala and P. Gambardella, “Ferrimagnetic Dynamics Induced by Spin-Orbit Torques,” Adv. Mater. Interfaces 9, 2201622 (2022a).
  * Olejník  _et al._ (2018) K. Olejník, T. Seifert, Z. Kašpar, V. Novák, P. Wadley, R. P. Campion, M. Baumgartner, P. Gambardella, P. Němec, J. Wunderlich, J. Sinova, P. Kužel, M. Müller, T. Kampfrath, and T. Jungwirth, “Terahertz electrical writing speed in an antiferromagnetic memory,” Sci. Adv. 4, eaar3566 (2018).
  * Tsai  _et al._ (2020) H. Tsai, T. Higo, K. Kondou, T. Nomoto, A. Sakai, A. Kobayashi, T. Nakano, K. Yakushiji, R. Arita, S. Miwa, Y. Otani, and S. Nakatsuji, “Electrical manipulation of a topological antiferromagnetic state,” Nature 580, 608–613 (2020).
  * Krishnaswamy  _et al._ (2022) G. K. Krishnaswamy, G. Sala, B. Jacot, C.-H. Lambert, R. Schlitz, M. D. Rossell, P. Nöel, and P. Gambardella, “Time-Dependent Multistate Switching of Topological Antiferromagnetic Order in Mn3Sn,” Phys. Rev. Appl. 18, 024064 (2022).
  * Krizakova  _et al._ (2022) V. Krizakova, M. Perumkunnil, S. Couet, P. Gambardella, and K. Garello, “Spin-orbit torque switching of magnetic tunnel junctions for memory applications,” J. Magn. Magn. Mater. 562, 169692 (2022).
  * Bläsing  _et al._ (2020) R. Bläsing, A. A. Khan, P. C. Filippou, C. Garg, F. Hameed, J. Castrillon, and S. S. P. Parkin, “Magnetic Racetrack Memory: From Physics to the Cusp of Applications Within a Decade,” Proc. IEEE 108, 1303–1321 (2020).
  * Luo  _et al._ (2020) Z. Luo, A. Hrabec, T. P. Dao, G. Sala, S. Finizio, J. Feng, S. Mayr, J. Raabe, P. Gambardella, and L. J. Heyderman, “Current-driven magnetic domain-wall logic,” Nature 579, 214–218 (2020).
  * Demidov  _et al._ (2017) V. E. Demidov, S. Urazhdin, G. de Loubens, O. Klein, V. Cros, A. Anane, and S. O. Demokritov, “Magnetization oscillations and waves driven by pure spin currents,” Phys. Rep. 673, 1–31 (2017).
  * Zahedinejad  _et al._ (2022) M. Zahedinejad, H. Fulara, R. Khymyn, A. Houshang, M. Dvornik, S. Fukami, S. Kanai, H. Ohno, and J. Åkerman, “Memristive control of mutual spin Hall nano-oscillator synchronization for neuromorphic computing,” Nat. Mater. 21, 81–87 (2022).
  * Pai  _et al._ (2012) C.-F. Pai, L. Liu, Y. Li, H. W. Tseng, D. C. Ralph, and R. A. Buhrman, “Spin transfer torque devices utilizing the giant spin Hall effect of tungsten,” Appl. Phys. Lett. 101 (2012).
  * Garello  _et al._ (2013) K. Garello, I. M. Miron, C. O. Avci, F. Freimuth, Y. Mokrousov, S. Blügel, S. Auffret, O. Boulle, G. Gaudin, and P. Gambardella, “Symmetry and magnitude of spin–orbit torques in ferromagnetic heterostructures,” Nat. Nanotechnol. 8, 587–593 (2013).
  * Kim  _et al._ (2013) J. Kim, J. Sinha, M. Hayashi, M. Yamanouchi, S. Fukami, T. Suzuki, S. Mitani, and H. Ohno, “Layer thickness dependence of the current-induced effective field vector in Ta/CoFeB/MgO,” Nat. Mater. 12, 240 (2013).
  * Ishikuro  _et al._ (2019) Y. Ishikuro, M. Kawaguchi, N. Kato, Y.-C. Lau, and M. Hayashi, “Dzyaloshinskii-Moriya interaction and spin-orbit torque at the Ir/Co interface,” Phys. Rev. B 99, 134421 (2019).
  * Miron  _et al._ (2010) I. M. Miron, G. Gaudin, S. Auffret, B. Rodmacq, A. Schuhl, S. Pizzini, J. Vogel, and P. Gambardella, “Current-driven spin torque induced by the Rashba effect in a ferromagnetic metal layer,” Nat. Mater. 9, 230–234 (2010).
  * Sinova  _et al._ (2015) J. Sinova, S. O. Valenzuela, J. Wunderlich, C.-H. Back, and T. Jungwirth, “Spin Hall effects,” Rev. Mod. Phys. 87, 1213 (2015).
  * Bihlmayer  _et al._ (2022) G. Bihlmayer, P. Noël, D. V. Vyalikh, E. V. Chulkov, and A. Manchon, “Rashba-like physics in condensed matter,” Nat. Rev. Phys. 4, 642–659 (2022).
  * Nguyen, Ralph, and Buhrman (2016) M.-H. Nguyen, D. C. Ralph, and R. A. Buhrman, “Spin Torque Study of the Spin Hall Conductivity and Spin Diffusion Length in Platinum Thin Films with Varying Resistivity,” Phys. Rev. Lett. 116, 126601 (2016).
  * Miao  _et al._ (2013) B. F. Miao, Y. Huang, D. Qu, and C. L. Chien, “Inverse Spin Hall Effect in a Ferromagnetic Metal,” Phys. Rev. Lett. 111, 066602 (2013).
  * Du  _et al._ (2014) C. Du, H. Wang, F. Yang, and P. C. Hammel, “Systematic variation of spin-orbit coupling with d-orbital filling: Large inverse spin Hall effect in 3d transition metals,” Phys. Rev. B 90, 140407(R) (2014).
  * Qu, Huang, and Chien (2015) D. Qu, S. Y. Huang, and C. L. Chien, “Inverse spin Hall effect in Cr: Independence of antiferromagnetic ordering,” Phys. Rev. B 92, 020418(R) (2015).
  * An  _et al._ (2016) H. An, Y. Kageyama, Y. Kanno, N. Enishi, and K. Ando, “Spin–torque generator engineered by natural oxidation of Cu,” Nat. Commun. 7, 13069 (2016).
  * Wang  _et al._ (2017) T. Wang, W. Wang, Y. Xie, M. A. Warsi, J. Wu, Y. Chen, V. O. Lorenz, X. Fan, and J. Q. Xiao, “Large spin Hall angle in vanadium film,” Sci. Rep. 7, 1306 (2017).
  * h. C. Baek  _et al._ (2018) S. h. C. Baek, V. P. Amin, Y.-W. Oh, G. Go, S.-J. Lee, G.-H. Lee, K.-J. Kim, M. D. Stiles, B.-G. Park, and K.-J. Lee, “Spin currents and spin–orbit torques in ferromagnetic trilayers,” Nat. Mater. 17, 509–513 (2018).
  * Hibino  _et al._ (2020) Y. Hibino, T. Taniguchi, K. Yakushiji, A. Fukushima, H. Kubota, , and S. Yuasa, “Large Spin-Orbit-Torque Efficiency Generated by Spin Hall Effect in Paramagnetic Co-Ni-B Alloys,” Phys. Rev. Appl. 14, 064056 (2020).
  * Kontani  _et al._ (2009) H. Kontani, T. Tanaka, D. S. Hirashima, K. Yamada, and J. Inoue, “Giant Orbital Hall Effect in Transition Metals: Origin of Large Spin and Anomalous Hall Effects,” Phys. Rev. Lett. 102, 016601 (2009).
  * Jo, Go, and Lee (2018) D. Jo, D. Go, and H.-W. Lee, “Gigantic intrinsic orbital hall effects in weakly spin-orbit coupled metals,” Phys. Rev. B 98, 214405 (2018).
  * Go  _et al._ (2020) D. Go, F. Freimuth, J.-P. Hanke, F. Xue, O. Gomonay, K.-J. Lee, S. Blügel, P. M. Haney, H.-W. Lee, and Y. Mokrousov, “Theory of current-induced angular momentum transfer dynamics in spin-orbit coupled systems,” Phys. Rev. Res. 2, 033401 (2020).
  * Lee  _et al._ (2021a) D. Lee, D. Go, H.-J. Park, W. Jeong, H.-W. Ko, D. Yun, D. Jo, S. Lee, G. Go, J. H. Oh, K.-J. Kim, B.-G. Park, B.-C. Min, H. C. Koo, H.-W. Lee, O. Lee, and K.-J. Lee, “Orbital torque in magnetic bilayers,” Nat. Commun. 12, 6710 (2021a).
  * Lee  _et al._ (2021b) S. Lee, M.-G. Kang, D. Go, D. Kim, J.-H. Kang, T. Lee, G.-H. Lee, J. Kang, N. J. Lee, Y. Mokrousov, S. Kim, K.-J. Kim, K.-J. Lee, and B.-G. Park, “Efficient conversion of orbital Hall current to spin current for spin-orbit torque switching,” Commun. Phys. 4, 234 (2021b).
  * Sala and Gambardella (2022b) G. Sala and P. Gambardella, “Giant orbital Hall effect and orbital-to-spin conversion in $3d$, $5d$, and $4f$ metallic heterostructures,” Phys. Rev. Res. 4, 033037 (2022b).
  * Zhu, Ralph, and Buhrman (2021) L. Zhu, D. C. Ralph, and R. A. Buhrman, “Maximizing spin-orbit torque generated by the spin Hall effect of Pt,” Appl. Phys. Rev. 8 (2021).
  * Hu and Pai (2020) C.-Y. Hu and C.-F. Pai, “Benchmarking of Spin–Orbit Torque Switching Efficiency in Pt Alloys,” Adv. Quantum Technol. 3, 2000024 (2020).
  * Chi  _et al._ (2021) Z. Chi, Y.-C. Lau, M. Kawaguchi, and M. Hayashi, “Charge–spin conversion in Pt1-xBix alloys for spin–orbit torque switching,” APL Mater. 9 (2021).
  * Wang  _et al._ (2022) P. Wang, A. Migliorini, S.-H. Yang, J.-C. Jeon, I. Kostanovskiy, H. Meyerheim, H. Han, H. Deniz, and S. S. P. Parkin, “Giant Spin Hall Effect and Spin–Orbit Torques in 5d Transition Metal–Aluminum Alloys from Extrinsic Scattering,” Adv. Mater. 34, 2109406 (2022).
  * Niimi  _et al._ (2012) Y. Niimi, Y. Kawanishi, D. H. Wei, C. Deranlot, H. X. Yang, M. Chshiev, T. Valet, A. Fert, and Y. Otani, “Giant Spin Hall Effect Induced by Skew Scattering from Bismuth Impurities inside Thin Film CuBi Alloys,” Phys. Rev. Lett. 109, 156602 (2012).
  * Zhu, Ralph, and Buhrman (2018) L. Zhu, D. C. Ralph, and R. A. Buhrman, “Highly Efficient Spin-Current Generation by the Spin Hall Effect in ${\mathrm{Au}}_{1-x}{\mathrm{Pt}}_{x}$,” Phys. Rev. Appl. 10, 031001 (2018).
  * Bozorth (1993) R. M. Bozorth,  _Ferromagnetism_ (Wilei, New York, 1993).
  * Ahern, Martin, and Sucksmith (1958) S. A. Ahern, M. J. C. Martin, and W. Sucksmith, “The spontaneous magnetization of nickel + copper alloys,” Proc. R. Soc. Lond. A 248, 145–152 (1958).
  * Keller  _et al._ (2019) M. W. Keller, K. S. Gerace, M. Arora, E. K. Delczeg-Czirjak, J. M. Shaw, and T. J. Silva, “Near-unity spin Hall ratio in NixCu1-x alloys,” Phys. Rev. B 99, 214411 (2019).
  * Varotto  _et al._ (2020) S. Varotto, M. Cosset-Chéneau, C. Grèzes, Y. Fu, P. Warin, A. Brenac, J.-F. Jacquot, S. Gambarelli, C. Rinaldi, V. Baltz, J.-P. Attané, L. Vila, and P. Noël, “Independence of the Inverse Spin Hall Effect with the Magnetic Phase in Thin NiCu Films,” Phys. Rev. Lett. 125, 267204 (2020).
  * Wu  _et al._ (2022) P.-H. Wu, D. Qu, Y.-C. Tu, Y.-Z. Lin, C. Chien, and S.-Y. Huang, “Exploiting Spin Fluctuations for Enhanced Pure Spin Current,” Phys. Rev. Lett. 128, 227203 (2022).
  * Cheng  _et al._ (2022) H. Cheng, Y. Wang, H. He, Q. Huang, and Y. Lu, “Efficient and temperature-independent terahertz emission from CoFeB/NiCu heterostructures,” Phys. Rev. B 105, 155141 (2022).
  * Avci  _et al._ (2014) C. O. Avci, K. Garello, M. Gabureac, A. Ghosh, A. Fuhrer, S. F. Alvarado, and P. Gambardella, “Interplay of spin-orbit torque and thermoelectric effects in ferromagnet/normal-metal bilayers,” Phys. Rev. B 90, 224427 (2014).
  * Omori  _et al._ (2019) Y. Omori, E. Sagasta, Y. Niimi, M. Gradhand, L. E. Hueso, F. Casanova, and Y. Otani, “Relation between spin Hall effect and anomalous Hall effect in 3d ferromagnetic metals,” Phys. Rev. B 99, 014403 (2019).
  * Haltz  _et al._ (2018) E. Haltz, R. Weil, J. Sampaio, A. Pointillon, O. Rousseau, K. March, N. Brun, Z. Li, E. Briand, C. Bachelet, Y. Dumont, and A. Mougin, “Deviations from bulk behavior in TbFe(Co) thin films: Interfaces contribution in the biased composition,” Phys. Rev. Mater. 2, 104410 (2018).
  * Zhang  _et al._ (2020) Q. Zhang, D. Zheng, Y. Wen, Y. Zhao, W. Mi, A. Manchon, O. Boulle, and X. Zhang, “Effect of surface roughness on the anomalous Hall effect in Fe thin films,” Phys. Rev. B 101, 134412 (2020).
  * Dang  _et al._ (2020) T. H. Dang, Q. Barbedienne, D. Q. To, E. Rongione, N. Reyren, F. Godel, S. Collin, J. M. George, and H. Jaffrès, “Anomalous Hall effect in $3d/5d$ multilayers mediated by interface scattering and nonlocal spin conductivity,” Phys. Rev. B 102, 144405 (2020).
  * Bello  _et al._ (2022) J. Bello, D. Lacour, S. Migot, J. Ghanbaja, S. Mangin, and M. Hehn, “Impact of interfaces on magnetic properties of Gdx(Fe90Co10)1-x alloys,” Appl. Phys. Lett. 121, 212402 (2022).
  * Hayashi  _et al._ (2014) M. Hayashi, J. Kim, M. Yamanouchi, and H. Ohno, “Quantitative characterization of the spin-orbit torque using harmonic Hall voltage measurements,” Phys. Rev. B 89, 144425 (2014).
  * Hsu  _et al._ (1996) S.-Y. Hsu, P. Holody, R. Loloee, J. M. Rittner, W. P. Pratt, and P. A. Schroeder, “Spin-diffusion lengths of Cu1-xNix using current perpendicular to plane magnetoresistance measurements of magnetic multilayers,” Phys. Rev. B 54, 9027 (1996).
  * Céspedes-Berrocal  _et al._ (2021) D. Céspedes-Berrocal, H. Damas, S. Petit-Watelot, D. Maccariello, P. Tang, A. Arriola-Córdova, P. Vallobra, Y. Xu, J.-L. Bello, E. Martin, S. Migot, J. Ghanbaja, S. Zhang, M. Hehn, S. Mangin, C. Panagopoulos, V. Cros, A. Fert, and J.-C. Rojas-Sánchez, “Current-Induced Spin Torques on Single GdFeCo Magnetic Layers,” Adv. Mater. 33, 2007047 (2021).