# Discovering dependencies in complex physical systems using Neural Networks

Sachin Kasture<EMAIL_ADDRESS>OptoAI, 2805DL, Gouda, The
Netherlands

###### Abstract

In todays age of data, discovering relationships between different variables
is an interesting and a challenging problem. This problem becomes even more
critical with regards to complex dynamical systems like weather forecasting
and econometric models, which can show highly non-linear behaviour. A method
based on mutual information and deep neural networks is proposed as a
versatile framework for discovering non-linear relationships ranging from
functional dependencies to causality. We demonstrate the application of this
method to actual multivariable non-linear dynamical systems. We also show that
this method can find relationships even for datasets with small number of
datapoints, as is often the case with empirical data.

††preprint: APS/123-QED

##  I Introduction

Finding relationships between different variables in large datasets [1, 2, 3]
is an important problem that has ramifications in fields ranging from
environmental science to economics and genetic networks. Understanding what
variables affect a certain quantity becomes increasingly challenging when
these relationships are highly non-linear, like those occurring in dynamical
systems with several variables. Quite often in a large dataset with several
variables, only a few variables maybe significantly affecting the target
variable and identifying these variables is first vital step in exploring
these dependencies in more detail.

Several methods exist which can help find dependencies and correlations
between variables. However most of these methods are good at detecting a
certain class of functions while they fail for others. There are some methods
which are quite good at detecting functional dependencies between 2 variables
[1, 4], they have however not been demonstrated in a multi-variable scenario
where a target variable depends on several input variables. Finding functional
dependencies has been a topic explored extensively in context of relational
databases[5, 6]. However these methods rely on finding exact functional
relationships by finding all attributes which have a one to one or one to many
relationship with a certain column Y. But this approach does not work well for
small databases which are just a sample of the true distribution as in these
cases one to one relations are more likely to occur. Also in such cases, it is
difficult to reliably find the smallest subset of variables which are
sufficient to describe Y. These methods do not offer any control over what
kind of functional relationships maybe considered intuitively as good or
interesting candidates. Also, these methods do not provide any kind of score
to evaluate functional dependencies.

In this paper, we use Neural networks as devices to model nonlinear behavior
and find complex non-linear relationships. Especially deep neural networks
(DNN) which consist of more than 1 hidden layer are excellent candidates for
efficiently modelling multi-variable non-linear polynomial functions with
small number of neurons [7, 8]. Additionally a regularization mechanism allows
us to control the complexity of the model we wish to consider [9]. Neural
networks have been used recently to discover physical concepts, identify phase
transitions and design quantum experiments[10, 11, 12]. To help find
dependencies, we use an DNN based autoencoder architecture which consists of
an encoder-decoder pair. The encoder maps the input space to a latent space,
while the decoder maps the latent space to the output space. This architecture
has been used, amongst other applications, for non-linear Principle Component
analysis (PCA) where the goal is to find a compressed representation of data
[13]. As such the input and the output of the autoencoder is conventionally
the same. In our method the input will be $X$, which is the set of input
features and $Y$ is the target feature or the set of features. We then use
compression of mutual information in the latent space to derive a loss
function which can be minimized to find the smallest set of features in $X$
which can be used to reliably reconstruct $Y$. The loss function can be used
to assign a score to compare the functional dependencies on different set of
input parameters.We then demonstrate this method to find dependencies in
chaotic dynamical systems. Also we show that this method can be used to find
non-linear causal connections in the Grangier sense for chaotic systems [14,
15, 16], even for a small dataset of 100 samples.

Figure 1:  Plot shows comparison between $x_{i}$ and the corresponding scaled
version of $l_{i}$ for (a)-(d) different values of $y_{i}=dx_{i}/dt$ for
equation 17. In the plots where $l_{i}$ is essentially noise, information from
the corresponding $x_{i}$ is not used to reconstruct $y_{i}$ using the
decoder. $fac$ is a scaling factor chosen so that $x_{i}$ and $l_{i}/fac$ are
comparable

##  II Theory

We now derive a loss function using the information bottleneck method [17]
based on the fact that the latent intermediate layer can be used to extract
only relevant information from $X$ and used to reconstruct $Y$. We represent
this latent representation by $L$. We also now assume a Markov chain
$Y\rightarrow X\rightarrow\ L$. This means $P(Y|X,L)=P(Y|X)$. This is because
$X,Y$ correspond to observed ground truth data.We now use the fact that we
want to extract only relevant information from $X$ which can reconstruct $Y$.
We use Shannon mutual information to quantify this information [17, 18].
Therefore want to maximize the quantity $I(L,Y)-\lambda_{enc}I(L,X)$. The
first term and the second term describe the capacity of the encoder and the
decoder respectively with $\lambda_{enc}$ determining the relative weight
between the two terms. We can write $I(L,Y)$ as:

$\begin{split}&I(L,Y)=\int dydlp(y,l)log\frac{p(y|l)}{p(y)}\\\ &=\int
dlp(l)\int dyp(y|l)log(p(y|l)+H(Y)\\\ \end{split}$ (1)

where $H(Y)$ is the Shannon entropy. We neglect $H(Y)$ since it is fixed by
the data. Since it is very difficult to calculate $p(y|l)$, we can approximate
it by another analytic function $\phi(y|l)$. Using the fact that the KL
divergence which measures the ‘distance’ between 2 probability distributions
is always non-negative:

$\begin{split}&KL(p(y|l),\phi(y|l))\geq 0\\\ &\implies\int
dyp(y|l)logp(y|l)\geq\int dyp(y|l)log\phi(y|l)\\\ \end{split}$ (2)

we can write

$I(L,Y)\geq\int dydlp(y,l)log\phi(y|l)$ (3)

We can now choose an appropriate function for $\phi(y|l)$ which allows us to
derive a suitable loss function as well as allows us to tune the complexity of
the decoder. The output of the decoder is given by $\theta_{dec}(l)$ which
describes the composite function of the decoder neural network which acts on
the latent variable $l$. To also include an additional L1 [9]regulation
parameter which helps restrict the magnitude of the weights in the decoder
neural network, we use the following function for $\phi(y|l)$

$\phi(y|l)=e^{-(\theta_{dec}(l)-y)^{2}/\sigma_{dec}^{2}-\lambda_{dec}(|\theta_{d1}|+|\theta_{d2}|+..)}$
(4)

where $\theta_{d1},\theta_{d2}..$ etc. are weights of different neurons in the
decoder network. Therefore we can write

$\begin{split}I(L,Y)&\geq-\int
dydlp(y,l)[\frac{(\theta_{dec}(l)-y)^{2}}{\sigma_{dec}^{2}}\\\
&+\lambda_{dec}(|\theta_{d1}|+|\theta_{d2}|+..)]\\\ \end{split}$ (5)

Now we use the fact that $p(y,l)=\int dxp(x,y,l)=\int dxp(l|x,y)p(x,y)$. Using
the Markov chain condition, this can be written as $p(y,l)=\int
dxp(l|x)p(x,y)$. Approximating $\int
dxdyp(x,y)A(x,y)=(1/M)\sum_{k=1}^{M}A(x^{k},y^{k})$ where $M$ is the number of
distinct data points, we can write

$\begin{split}I(L,Y)&\geq-(1/M)\sum_{k=1}^{M}\int
dlp(l|x)[\frac{(\theta_{dec}(l)-y^{k})^{2}}{\sigma_{dec}^{2}}\\\
&+\lambda_{dec}(|\theta_{d1}|+|\theta_{d2}|+..)]\\\ \end{split}$ (6)

Similarly we can define $I(L,X)$ as:

$\begin{split}I(L,X)&=\int dldxp(x,l)log\frac{p(l|x)}{p(l)}\\\ &=\int
dxdlp(x,l)logp(l|x)-\int dlp(l)logp(l)\\\ \end{split}$ (7)

We now again use another analytical function $g(l)$ in place of $p(l)$ and use
the result on positivity of KL divergence and get:

$\begin{split}I(L,X)&=\int dldxp(x,l)logp(l|x)-\int p(l)logp(l)\\\ &\leq\int
dxdlp(x,l)log\frac{p(l|x)}{g(l)}\\\ \end{split}$ (8)

For convenience we use a Gaussian function centred at 0.

$g(l)=e^{-\sum_{i}l_{i}^{2}/\sigma_{enc}^{2}}$ (9)

where $l=(l_{1},l_{2}..)$ are different components of $l$ and $\sigma_{enc}$
is an adjustable parameter. For $p(l|x)$ we can use:

$p(l|x)=\prod_{i}e^{-(l_{i}-W_{i}x_{i})^{2}/\sigma_{enc}^{2}}$ (10)

where $x=(x_{1},x_{2},..)$ This means we use a linear transformation from $X$
and add a independent Gaussian noise with variance $\sigma_{enc}^{2}$ and mean
0 to each component. We now plug in definitions 9,10 into equation 8 and
obtain:

$I(L,X)\leq\int
dxdlp(x,l)loge^{-\sum_{i}W_{i}x_{i}(W_{i}x_{i}-2l_{i})/\sigma_{enc}^{2}}$ (11)

Writing $p(x,l)=p(x)p(l|x)$ we can write the above equation as

$\begin{split}I(L,X)&\leq-\int
dxdlp(x)\prod_{i}e^{-(l_{i}-W_{i}x_{i})^{2}/\sigma_{enc}^{2}}\\\
&[\frac{\sum_{i}W_{i}x_{i}(W_{i}x_{i}-2l_{i})}{\sigma_{enc}^{2}}]\\\
\end{split}$ (12)

Using the approximation $\int dxp(x)A(x)=(1/M)\sum_{k=1}^{M}A(x^{k})$, we can
write

$\begin{split}I(L,X)&\leq-(1/M)\sum_{k=1}^{M}\int
dl\prod_{i}e^{-(l_{i}-W_{i}x_{i}^{k})^{2}/\sigma_{enc}^{2}}\\\
&[\frac{\sum_{i}W_{i}x_{i}^{k}(W_{i}x_{i}^{k}-2l_{i})}{\sigma_{enc}^{2}}]\\\
\end{split}$ (13)

Similarly substituting equation 10 into equation 6 and assuming
$\sigma_{enc}^{2}$ to be small enough so that
$e^{-(l_{i}-W_{i}x_{i})^{2}/\sigma_{enc}^{2}}\approx\delta(l_{i}-W_{i}x_{i})$we
obtain:

$\begin{split}I(L,Y)&-\lambda_{enc}I(L,X)\geq-(1/M)\sum_{k=1}^{M}[\frac{(\theta_{dec}(l)-y^{k})^{2}}{\sigma_{dec}^{2}}+\\\
&\lambda_{dec}(|\theta_{d1}|+|\theta_{d2}|+..)+\lambda_{enc}\sum_{i}\frac{(W_{i}x_{i}^{k})^{2}}{\sigma_{enc}^{2}}]\\\
\end{split}$ (14)

Figure 2:  Plots shows the case of fan-in causality pattern for set of delay
equations in equation 18 for set of $\xi_{ij}$ values used to obtain results
in Figure 3 Figure 3:  Plot shows comparison between $Y_{i}$ and the
corresponding scaled version of $l_{i}$ for (a)-(c) different values of
$y_{i}=Y_{i}$ for the set of delay equations 18. In the plots where $l_{i}$ is
noise, information from the corresponding $x_{i}$ is not used to reconstruct
$y_{i}$ using the decoder

Therefore we can define a loss function to be minimized as

$\begin{split}\mathcal{L}&=(1/M)\sum_{k=1}^{M}[\frac{(\theta_{dec}(l)-y^{k})^{2}}{\sigma_{dec}^{2}}+\\\
&\lambda_{dec}(|\theta_{d1}|+|\theta_{d2}|+..)+\lambda_{enc}\sum_{i}\frac{(W_{i}x_{i}^{k})^{2}}{\sigma_{enc}^{2}}]\\\
\end{split}$ (15)

We observe that the first term tries to minimize the least squares difference
between $\theta_{dec}(l)$ and $y$ and the second term controls the size of the
weights of the decoder which in turn controls the maximum degree polynomials
the decoder NN can approximate. For the third term we see that as we increase
the $\lambda_{inc}$, the NN will try to keep $(W_{i}x_{i}^{k})^{2}$ small to
keep the total loss function small. Assuming now that we standardize our data
so that $x_{i}^{\prime}s$ on an average have similar magnitudes, we absorb it
into $\lambda_{enc}$. The third term will now be smallest when only
$W_{i}^{\prime}s$ corresponding to those $x_{i}^{\prime}s$ are non-zero, which
are required to reproduce $Y$. Using this intution and the fact that term
inside the summation over $i$ in equation 17 is always $\geq 0$, we can
further simplify the loss function as

$\begin{split}\mathcal{L}&=(1/M)\sum_{k=1}^{M}[\frac{(\theta_{dec}(l)-y^{k})^{2}}{\sigma_{dec}^{2}}]+\\\
&\lambda_{dec}(|\theta_{d1}|+|\theta_{d2}|+..)+\lambda_{enc}\sum_{i}(|W_{i}|)\\\
\end{split}$ (16)

where we have merged $\sigma_{enc}^{2}$ with $\lambda_{enc}$. This way we
treat both the encoder and decoder weights on equal terms using L1
regularization. From a practical standpoint L1 is advantageous since it can
shrink weights faster.

##  III Application

For further study we use a NN in which the encoder has 2 linear layers. This
gives us a mapping $X\rightarrow L$. We then add Gaussian noise to the latent
variables $l_{i}=l_{i}+N(0,\sigma^{2}_{enc})$. The latent code is then sent
through a multilayer decoder network with non-linear activation functions to
give the output $\theta_{dec}(l)$. We perform batch-normalization in between
intermediate neural network layers [19]. This layers prevents change in data
distributions between adjacent layers and allows neural network learning at a
higher learning rate. We then minimize the loss function in equation 16 using
Stochastic gradient descent with different batch sizes. We can tune the values
of $\lambda_{enc},\lambda_{dec}$ (regularization parameters) to obtain as low
values of loss function as possible. This choice of regularization parameters
may also depend on our prior knowledge about the complexity of the system. The
data is split into the training and validation set. The training data is used
to build the model and validation set checks how well the model generalizes.
The basic heuristic for tuning these parameters is as follows: after fixing
the learning rate for the gradient descent, we first increase the value of
$\lambda_{dec}$ which basically fixes the complexity of functions the decoder
can simulate. We then increase the value of $\lambda_{enc}$ and look at the
value of the mean square error and stop when the mean square error is as small
as possible for both the training and the validation set. We now use this
method to infer relationships in well known non-linear systems. We first
consider a Lorenz96 non-linear system which is defined as:

$\frac{dx_{i}}{dt}=(x_{i+1}-x_{i-2}x_{i-1}-x_{i}+F)$ (17)

where $i$ goes from $1$ to $N$ where $N$ is the number of oscillators and
$x_{N+1}=x_{1}$,$x_{-1}=x_{N-1}$, $x_{0}=x_{N}$. $F$ is the driving term and
we choose $F=8$ where the system behaves in the chaotic regime. Figure 1 shows
the results for N=5. We run N=5 times with each time $y=\frac{dx_{i}}{dt}$ for
i from 1 to 5. We see that the latent representation $l_{i}$ is basically just
the added Gaussian noise when the corresponding $y$ has no dependency on
$l_{i}$. The number of data points was 3000 and learning rate was 0.0001 and
values of $\lambda_{dec},\lambda_{enc}$ where 0 and 0.1 respectively. The
training was run for 1000 epochs with a batch size of 300.
Next we apply NN to infer causal relationship in a set of non-linear delay
equations. For this we look at the following set of equations:

$Y_{i}(t+1)=(\xi_{ii}-\sum_{j=1,2,3}(\xi_{jj}-\xi_{ij}Y_{j}(t)))Y_{i}(t)$ (18)

for i=1,2,3. We choose to choose parameters $\xi_{ij}$ which correspond to a
fan-in pattern shown in Figure 2. The values of $\xi$ are as follows
$\xi_{11}=4,\xi_{22}=3,\xi_{33}=2,\xi_{31}=0.6,\xi_{32}=-0.6$. These
parameters corresponds to a chaotic regime. In this case both $Y_{2}$ and
$Y_{3}$ are causally driven by $Y_{1}$. A fan-in pattern is a good test
because correlation based tests would falsely infer a causal relationship
between $Y_{2}$ and $Y_{3}$ [2]. To infer the causal relationships, we run the
NN with $y=Y_{i}(t+1)$ and input $X=[Y_{1}(t),Y_{2}(t),Y_{3}(t)]$. From Figure
3 we can see that we are able to correctly infer the dependencies, even for a
very small data-set of 50 points. The plots were obtained for a learning rate
of 0.001 and $\lambda_{enc},\lambda_{dec}$ values of 0.1 and 0.005
respectively.The number of epochs was 1500 with a batch size of 32.

Figure 4:  Plot shows the plot for FD vs MR for different values of
$\lambda_{enc}$. The legend also mentions the non-linear system for the
plotted data. ‘dde’ stands for the delay difference equations in equation 18

We also summarize the performance of this method using 2 metrics False
discovery (FD) and Miss rate (MR) which are defined as:

$\begin{split}&FD=\frac{FP}{FP+TP}\\\ &MR=\frac{FN}{FN+TP}\\\ \end{split}$
(19)

where FN, FP, TP are False negatives, false positives and true positives
respectively. Here a positive means a certain variable has been discovered to
be independent of the output. The negative means a variable has been
discovered to be related to the output.This data is obtained by obtaining
results over 20 independent runs of the model. For the Lorenz96 model, the
best result is obtained with $\lambda_{enc}=0.2$ while for the set of
equations 18, best results are obtained for $\lambda_{enc}=0.1$

##  IV Conclusion

The proposed approach using NN is a versatile platform for inferring
relationships, especially in complex non-linear systems. This is because NN
are a powerful tool to model such non-linear functions. Even though it is
difficult to infer the exact functional form using a NN, this method can help
locate functional dependencies between variables in a multivariable system.
These variables can then be probed more extensively to find the functional (or
approximate functional) form of the relationships. Methods based on sparse
regression have been used in the past to find functional relationships.
However they rely on pre-knowledge of the set of basis functions to use for
the regression. The proposed method has no such requirement and with a large
enough NN, can simulate any complex non-linear function. Besides locating
functional relationships, it can also help infer causal relationships in non-
linear data as seen in the discussed example, where it correctly inferred
causal relationship even for a small dataset of 50 samples.

##  V Acknowledgements

The author would like to thank Akshatha Mohan for helpful comments and
critical assessment of the manuscript.

## References

  * Reshef  _et al._ [2011] D. N. Reshef, Y. A. Reshef, H. K. Finucane, S. R. Grossman, G. McVean, P. J. Turnbaugh, E. S. Lander, M. Mitzenmacher, and P. C. Sabeti, Science 334, 1518 (2011).
  * Marbach  _et al._ [2010] D. Marbach, R. J. Prill, T. Schaffter, C. Mattiussi, D. Floreano, and G. Stolovitzky, Proceedings of the National Academy of Sciences 107, 6286 (2010).
  * Brunton  _et al._ [2016] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Proceedings of the National Academy of Sciences 113, 3932 (2016).
  * Dembo  _et al._ [2001] A. Dembo, A. Kagan, and L. A. Shepp, Bernoulli 7, 343 (2001).
  * Liu  _et al._ [2012] J. Liu, J. Li, C. Liu, and Y. Chen, IEEE Transactions on Knowledge and Data Engineering 24, 251 (2012).
  * Huhtala [1999] Y. Huhtala, The Computer Journal 42, 100 (1999).
  * Lin  _et al._ [2017] H. W. Lin, M. Tegmark, and D. Rolnick, Journal of Statistical Physics 168, 1223 (2017).
  * Rolnick and Tegmark [2018] D. Rolnick and M. Tegmark, arXiv:1705.05502 [cs, stat] (2018), arXiv: 1705.05502.
  * Tibshirani [1996] R. Tibshirani, Journal of the Royal Statistical Society: Series B (Methodological) 58, 267 (1996).
  * Iten  _et al._ [2020] R. Iten, T. Metger, H. Wilming, L. del Rio, and R. Renner, Physical Review Letters 124, 010508 (2020).
  * Rem  _et al._ [2019] B. S. Rem, N. Käming, M. Tarnowski, L. Asteria, N. Fläschner, C. Becker, K. Sengstock, and C. Weitenberg, Nature Physics 15, 917 (2019).
  * Melnikov  _et al._ [2018] A. A. Melnikov, H. Poulsen Nautrup, M. Krenn, V. Dunjko, M. Tiersch, A. Zeilinger, and H. J. Briegel, Proceedings of the National Academy of Sciences 115, 1221 (2018).
  * Hinton [2006] G. E. Hinton, Science 313, 504 (2006).
  * Detto  _et al._ [2012] M. Detto, A. Molini, G. Katul, P. Stoy, S. Palmroth, and D. Baldocchi, The American Naturalist 179, 524 (2012).
  * Runge  _et al._ [2012] J. Runge, J. Heitzig, V. Petoukhov, and J. Kurths, Physical Review Letters 108, 258701 (2012).
  * Ma  _et al._ [2015] H. Ma, K. Aihara, and L. Chen, Scientific Reports 4, 7464 (2015).
  * Tishby  _et al._ [2000] N. Tishby, F. C. Pereira, and W. Bialek, arXiv:physics/0004057 (2000), arXiv: physics/0004057.
  * Giannella and Robertson [2004] C. Giannella and E. Robertson, Information Systems 29, 483 (2004).
  * Ioffe and Szegedy [2015] S. Ioffe and C. Szegedy, arXiv:1502.03167 [cs] (2015), arXiv: 1502.03167.


# Towards preserving word order importance through
Forced Invalidation

Hadeel Al-Negheimish1, Pranava Madhyastha2,1 Alessandra Russo1
1 Imperial College London 2 City, University of London
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

Large pre-trained language models such as BERT have been widely used as a
framework for natural language understanding (NLU) tasks. However, recent
findings have revealed that pre-trained language models are insensitive to
word order. The performance on NLU tasks remains unchanged even after randomly
permuting the word of a sentence, where crucial syntactic information is
destroyed. To help preserve the importance of word order, we propose a simple
approach called Forced Invalidation (FI): forcing the model to identify
permuted sequences as invalid samples. We perform an extensive evaluation of
our approach on various English NLU and QA based tasks over BERT-based and
attention-based models over word embeddings. Our experiments demonstrate that
FI significantly improves the sensitivity of the models to word order.111Our
code and data for replication are available at
https://github.com/halnegheimish/ForcedInvalidation

##  1 Introduction

Ordering of words in a sentence is an important structural attribute for
natural languages such as English, where subject-verb-object (SVO) structure
is common and important to convey the meaning of the sentence. Understanding
and comprehending natural language without strict adherence to a systematic
ordering of words would make it an extremely challenging task. Recent work has
investigated a surprising lack of sensitivity to word order information in
state-of-the-art masked language models.

Recent research has focused on the impact of word order perturbation during
the evaluation of models trained on well-ordered data. The results show that
masked language models exhibit a catastrophic lack of sensitivity to word
order permutations or shuffles, even for complex tasks in which task-relevant
syntactic properties are completely destroyed pham-etal-2020-order; al-
negheimish-etal-2021-numerical; sinha-etal-2021-unnatural; Gupta-
etal-2021-bert. These studies show that models are still predicting the gold
label for examples even after sequences have been permuted, and they also do
so with high confidence (sinha-etal-2021-unnatural; Gupta-etal-2021-bert).
This anomalous behaviour can potentially result in undesirable shortcuts or
can cause models to fail catastrophically in simple adversarial settings.
Furthermore, sinha-etal-2021-masked study the effect of pre-training masked
language models on shuffled data, and suggest that the model might simply be
capturing higher-order word co-occurrence statistics, rather than uncovering
sophisticated semantic and syntactic structures necessary for language
understanding.

In this paper, we present a simple, yet general approach called Forced
Invalidation (FI), where we force models to explicitly identify sequence
permutations (§3). While our proposal is extensible to multiple types of
models, here we present a controlled study over masked language model-based
BERT models, and attention-based models over word-embeddings (§4). We present
a large battery of experiments over a variety of natural language
understanding tasks in the English language, including complex question
answering based tasks, natural language inference based tasks and commonsense
reasoning based tasks. Results show that our proposal significantly improves
the sensitivity of the models to word-order information (§5).

##  2 Related Work

Recent work has proposed a few mitigation strategies for classification-type
tasks: pham-etal-2020-order proposes improving word-order sensitivity by
including a precursor fine-tuning step on synthetic CoLA-like tasks before
finetuning for downstream tasks. While this improves their defined word order
sensitivity score, accuracy on permuted samples remains significantly above
chance, making this approach unreliable. Gupta-etal-2021-bert present three
approaches: one based on entropy regularisation, another on model
probabilities thresholding, and finally an approach based on augmenting
additional data consisting of destructive transformations, which include
1-gram permutations. The model is modified with an additional class to
identify these destructive transformations. While the first two approaches
require changes to model setup, the last one is based on a set of manual
heuristics. Our approach has some similarity with the latter, however, our
permutations are based on the principals borrowed from the $n$-gram language
modelling literature roark2007discriminative, where the $n$-grams capture
sufficient first-order statistics of the language. Our approach significantly
reduces the models’ reliance only on simple first-order statistics of
language, and our empirical observations demonstrate the generalisability of
our approach to various settings.

Figure 1: Accuracy using exact match (EM) on the DROP QA dataset, (a) NAQANet,
(b) MTMSN Base, and (c) MTMSN Large. dev, numset and contrast sets are
unperturbed well-ordered datasets, while the other bars show {1,2,3}-gram
permutations on numset. FI models exhibit clear sensitivity to word order, as
they no longer predict the original answer in most perturbed cases. FI
accuracy is averaged over training with 3-random seeds, and the standard
deviation is shown over the yellow bars. FI-models are confident in
identifying invalid samples.

##  3 Approach: Forced Invalidation

Our method is grounded on recent observations which show that masked language
models and other similar models tend to exploit shortcuts based on information
about distributional word vicinity information for a diverse set of natural
language processing tasks al-negheimish-etal-2021-numerical; sinha-
etal-2021-masked. These shortcuts tend to make the models less sensitive to
word-order information even for tasks that require the preservation of word
ordering for accurate recovery of meaning sinha-etal-2021-unnatural. al-
negheimish-etal-2021-numerical, in particular, show the insensitivity on a
variety of $n$-gram based permutations of samples, where the authors
specifically experiment with {1,2,3}-gram permutations. Further, {1,2,3}-gram
based permutations capture a variety of word vicinity based patterns and also
capture some of the most frequently occurring unigram, bigram and trigram
patterns in a variety of benchmark datasets. Based on these salient
observations, for each given dataset, the _Forced Invalidation_ (FI)
methodology consists of the following two steps:

  1. 1.

Augmenting training data with {1,2,3}-gram permutation samples (sampled from
trainset) labelled with _invalid_ as the additional label.

  2. 2.

Modifying models to account for the new label and training them in the
standard setting with a combination of standard training examples and the
augmented _invalid_ samples.

We observe that this simple FI approach improves the sensitivity of the model
to word order and also improves the robustness of the models across a variety
of tasks to first-order shortcuts. In the following sections, we present a
rigorous experimental study that showcases the utility of FI.

##  4 Experimental settings

#### Data

To generate $n$-gram permutations, we simply subsample from the training
dataset, such that the ratio to the valid samples (samples with correct word
order and the task-specific label) and invalid samples (samples with permuted
n-grams and the _invalid_ label) is 1-1. The invalid samples are generated
such that they contain a uniform distribution of {1,2,3}-gram
permutations.222Input string is divided to n-grams (based on white-space), and
permuted, preserving the final punctuation. The only condition is that it
varies from the original string, which is a weaker constraint than previous
studies that require that no n-gram stays in its original position.
Furthermore, we split the training set such that we use 90% of the samples for
training the models and the remaining 10% is used as a development set. The
development set is used to monitor training and perform early-stopping.
Evaluations are done on a separate unseen task-specific validation set
provided by the dataset creators.

We perform two evaluations: well-ordered and permuted. The first one is the
standard evaluation of the model over the original unperturbed task-specific
unseen validation set. For our experimental evaluations over permutations, we
retain the original label of the same unseen validation set, but we permute
the specific components of samples (e.g., premise or the hypothesis, question
or the passage, etc.); we expect the models to reject the permuted sample
instead of predicting the same ground-truth label. This setup allows us to
evaluate the sensitivity of the model to word-order permutations of various
degrees over the different components of the data.

#### Models

We predominantly experiment with BERT-based models that either use BERT
representations as the contextualised embeddings or classifiers that are
directly trained with BERT. We also experiment with an additional simpler
model that largely exploits attention over word-embeddings.333Details about
training parameters can be found in the appendix LABEL:app:training We will
expand on the specific models for each of the tasks in the following section.
Results for FI are reported over the averages and standard deviation of models
trained with three random seeds.

##  5 Results and Observations

###  5.1 Unconstrained Question Answering

DROP Dua2019DROP is a reading-comprehension dataset with unconstrained
answers. It comprises questions that require reasoning over the content of
different paragraphs. Even though this is designed to be a challenging task,
al-negheimish-etal-2021-numerical show that for most models, permuting
questions had little impact on the model’s ability to predict the correct
answer for numerical reasoning questions. This was found specifically
problematic as the task consists of complex questions, and permutations render
the questions syntactically and semantically redundant. We apply FI for two
module-based models designed specifically for this task, NAQANet
(Dua2019DROP), an attention-based model with GloVe-based embeddings
(pennington-etal-2014-glove), which is the original model proposed by the
authors of the dataset; and MTMSN (hu-etal-2019-multi), which is based on
BERT-based (devlin2018bert) contextual representations. These models have
separate modules targeting different kinds of reasoning, e.g. a module for
counting and a module for arithmetic expressions. We augment these with an
additional module to force invalidation over invalid permuted samples, a two-
way classification module that learns to distinguish between permuting either
the question or passage. FI models should now choose the invalid type for
permuted samples, instead of giving the same answer as well-ordered samples.
Fig 1 shows a comparison of the Exact Match accuracy of NAQANet and MTMSN,
between original and FI models. We observe that question permutations have
little effect on the original model, as previously noted in (al-negheimish-
etal-2021-numerical). Interestingly, passage permutations can drastically
reduce performance by a third. While performance degrades for passage
permutations, it remains unacceptably high, as we note that DROP is an
unconstrained-QA task, so the space of possible answers is large. FI, on the
other hand, succeeds in making the model sensitive to almost all permutations.
We see that this generalises across BERT${}_{{\textsc{base}}}$ and
BERT${}_{{\textsc{large}}}$. More importantly, we observe that while the model
with FI is sensitive to word order (no longer predicting original answers for
permuted examples), the model’s performance on well-ordered data is largely
retained. To demonstrate that FI are correctly predicting invalid, table
LABEL:table:percentage-invalid shows the percentage of the data predicted as
invalid, which is near-perfect for permuted samples.


# Anisotropic resistivity tensor from disk geometry magneto-transport

Oskar Vafek<EMAIL_ADDRESS>National High Magnetic Field Laboratory,
Tallahassee, Florida, 32310, USA  Department of Physics, Florida State
University, Tallahassee, Florida 32306, USA

###### Abstract

Magneto-transport measurements on two dimensional van der Waals
heterostructures have recently shown signatures of uniaxial anisotropy. Such
measurements are almost exclusively performed in a Hall bar geometry which
makes it difficult to extract the full resistivity tensor. The goal of this
paper is to theoretically analyze anisotropic magneto-transport in a
homogeneous disk geometry and to provide a closed form expression for the
electrical potential anywhere on the disk if the current source and drain are
located somewhere on the circumference. This expression can then be used to
experimentally extract the full resistivity tensor.

##  I Introduction

Two dimensional van der Waals (vdW) heterostructures host a broad range of
interesting physical phenomena[1], including anisotropic magnetotransport.
With rare exceptions[2, 3], the transport measurements are performed in a Hall
bar geometry, making it difficult to extract the full resistivity tensor
particularly if the transport principal axis is misaligned with the current
flow. For example, the heterostructures can be subject to an unintentional
strain, in which case the misalignment is not directly controlled in an
experiment. Moreover, the orientation of the transport principal axis can be
carrier concentration (filling) dependent as was recently shown [4] in
numerical solutions of the Boltzman equation for twisted bilayer graphene
subject to heterostrain, even if the strain tensor and the transport
relaxation time are momentum and filling independent. For open Fermi surfaces,
the magneto-resistance is expected to grow with the magnetic field $B$ without
saturation along one of the principal axis, but to saturate with increasing
$B$ along the perpendicular principal axis[5]. Direct measurement of the full
anisotropic resistivity tensor in the vdW heterostructures as a function of
filling and $B$ would therefore help in understanding the complex transport
phenomena in these materials.

In this paper we analyze the solution to the magneto-transport equations in a
uniform disk of radius $a$. Inside the disk the anisotropic conductivity
tensor is assumed to be homogeneous, while outside the disk there is no
conduction. Thus,

$\displaystyle\sigma=D(x,y)\left(\sigma_{+}\hat{\mathbf{x}}\hat{\mathbf{x}}+\sigma_{-}\hat{\mathbf{y}}\hat{\mathbf{y}}+\sigma_{H}\left(\hat{\mathbf{x}}\hat{\mathbf{y}}-\hat{\mathbf{y}}\hat{\mathbf{x}}\right)\right),$
(1)

where $D(x,y)=\Theta(a^{2}-x^{2}-y^{2})$ and $\Theta$ is the Heaviside step
function, restricting the conduction to the interior of the circle. Without
loss of generality, we also assume that the principal axes are aligned with
the $x$ and $y$ axes of the coordinate system and adopt the dyadic product to
represent the conductivity tensor. Here $\sigma_{\pm}$ are the two components
of the longitudinal conductivity along the principal axes and $\sigma_{H}$ is
the Hall conductivity.

We express the longitudinal conductivities as
$\sigma_{\pm}=\bar{\sigma}\pm\Delta\sigma$ and without loss of generality take
the $x$-axis to be along the axis with larger resistivity i.e.
$\Delta\sigma/\bar{\sigma}<0$. The analysis below then provides an expression
in the form of a rapidly convergent series which can be used to extract the
resistivity tensor for a point current source/drain at
$\mathbf{r}_{S,D}=a(\cos\theta_{A,B},\sin\theta_{A,B})$. The expression for
the electrical potential at $x,y$ reads

$\displaystyle V(x,y;\mathbf{r}_{S},\mathbf{r}_{D})=$ (2)
$\displaystyle\frac{I}{\pi}\frac{\sqrt{\sigma_{+}\sigma_{-}}}{\sigma_{+}\sigma_{-}+\sigma^{2}_{H}}\left(\sum_{n=0,2,4,\ldots}^{\infty}\ln\frac{\left|1+e^{-2i\theta_{B}}\Omega^{2+4n}-e^{-i\theta_{B}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right|}{\left|1+e^{-2i\theta_{A}}\Omega^{2+4n}-e^{-i\theta_{A}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right|}+\sum_{n=1,3,5,\ldots}^{\infty}\ln\frac{\left|1+e^{2i\theta_{B}}\Omega^{2+4n}-e^{i\theta_{B}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right|}{\left|1+e^{2i\theta_{A}}\Omega^{2+4n}-e^{i\theta_{A}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right|}\right)$
$\displaystyle+$
$\displaystyle\frac{I}{\pi}\frac{\sigma_{H}}{\sigma_{+}\sigma_{-}+\sigma^{2}_{H}}\left(\sum_{n=0,2,4,\ldots}^{\infty}\arg\left(1+e^{-2i\theta_{B}}\Omega^{2+4n}-e^{-i\theta_{B}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right)-\arg\left(1+e^{-2i\theta_{A}}\Omega^{2+4n}-e^{-i\theta_{A}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right)\right.$
$\displaystyle+$
$\displaystyle\left.\sum_{n=1,3,5,\ldots}^{\infty}\arg\left(1+e^{2i\theta_{B}}\Omega^{2+4n}-e^{i\theta_{B}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right)-\arg\left(1+e^{2i\theta_{A}}\Omega^{2+4n}-e^{i\theta_{A}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right)\right).$

where the $x,y$ position enters via the complex variable
$Z=X+iY=\frac{x}{\sqrt{1+\frac{\Delta\sigma}{\bar{\sigma}}}}+i\frac{y}{\sqrt{1-\frac{\Delta\sigma}{\bar{\sigma}}}}$,
and the parameters
$\alpha_{+}=\frac{a}{2}\left(\frac{1}{\sqrt{1+\frac{\Delta\sigma}{\bar{\sigma}}}}+\frac{1}{\sqrt{1-\frac{\Delta\sigma}{\bar{\sigma}}}}\right)$
and
$\Omega=\sqrt{\frac{\sqrt{1-\frac{\Delta\sigma}{\bar{\sigma}}}-\sqrt{1+\frac{\Delta\sigma}{\bar{\sigma}}}}{\sqrt{1-\frac{\Delta\sigma}{\bar{\sigma}}}+\sqrt{1+\frac{\Delta\sigma}{\bar{\sigma}}}}}$.
The function $\arg$ is the argument of a complex number. Note that because
$|\Delta\sigma|<\bar{\sigma}$, the parameter $0\leq\Omega<1$ and therefore the
above sum converges (the convergence is rapid unless $\Omega$ is very close to
$1$). Illustrative contour plots of $V(x,y;\mathbf{r}_{S},\mathbf{r}_{D})$ for
several parameters are shown in the Fig.1.

Although the above expression is obtained for a point current source/drain,
the linearity of the differential equation whose solution it is allows direct
determination of the formula for multiple point, as well as spatially
extended, current sources/drains. Such formula is presented in the discussion
section.

(a) (b) (c) (d)

Figure 1: Equipotential contours computed using Eq.(2) for isotropic
conductivity tensor (a) for $\sigma_{H}=0$ and (b) for
$\sigma_{H}=0.4\bar{\sigma}$; in each case the point source is at
$\theta_{A}=\pi/2$ and the point drain at $\theta_{B}=-\pi/2$. Equipotential
contours for anisotropic conductivity tensor for $\sigma_{H}=0$ and
$\Delta\sigma=-0.7\bar{\sigma}$ for point drain at $\theta_{B}=-\pi/2$ and (c)
point source at $\theta_{A}=\pi/2$ and (d) $\theta_{A}=\pi/4$; ten terms in
the sum were kept.

##  II Analysis

The starting assumption is that Ohm’s law holds, i.e.

$\displaystyle\mathbf{j}$ $\displaystyle=$
$\displaystyle\sigma\cdot\mathbf{E}=-\sigma\cdot\nabla V,$ (3)

where $\mathbf{j}$ is the current density, $\mathbf{E}$ is the electric field
and $V$ is the electrical potential, all of which are assumed to be position
dependent. For an idealized point current source and drain, the continuity
equation gives

$\displaystyle\nabla\cdot\mathbf{j}$ $\displaystyle=$ $\displaystyle
I\left(\delta(\mathbf{r}-\mathbf{r}_{A})-\delta(\mathbf{r}-\mathbf{r}_{B})\right),$
(4)

where $I$ is the current, its source is at $\mathbf{r}_{A}$, its drain at
$\mathbf{r}_{B}$, and $\delta(\mathbf{r})$ is the Dirac delta function.
Combining Eqs.(1-4) gives

$\displaystyle-\frac{\partial}{\partial x}\left(D\sigma_{+}\frac{\partial
V}{\partial x}\right)-\frac{\partial}{\partial
y}\left(D\sigma_{-}\frac{\partial V}{\partial
y}\right)-\frac{\partial}{\partial x}\left(D\sigma_{H}\frac{\partial
V}{\partial y}\right)+\frac{\partial}{\partial
y}\left(D\sigma_{H}\frac{\partial V}{\partial
x}\right)=I\left(\delta(\mathbf{r}-\mathbf{r}_{A})-\delta(\mathbf{r}-\mathbf{r}_{B})\right).$
(5)

The solution to the above inhomogeneous linear partial differential equation
gives $V$ as a function of $\mathbf{r}$.

Expressing the longitudinal conductivities as

$\displaystyle\sigma_{\pm}$ $\displaystyle=$
$\displaystyle\bar{\sigma}\pm\Delta\sigma,$ (6)

it will be convenient to rescale the coordinate axes according to

$\displaystyle X$ $\displaystyle=$
$\displaystyle\frac{x}{\sqrt{1+\frac{\Delta\sigma}{\bar{\sigma}}}},$ (7)
$\displaystyle Y$ $\displaystyle=$
$\displaystyle\frac{y}{\sqrt{1-\frac{\Delta\sigma}{\bar{\sigma}}}},$ (8)

so that Eq.(5) becomes

$\displaystyle-\left(\frac{\partial}{\partial X}\left(D\frac{\partial
V}{\partial X}\right)+\frac{\partial}{\partial Y}\left(D\frac{\partial
V}{\partial
Y}\right)\right)-\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\left(\frac{\partial}{\partial
X}\left(D\frac{\partial V}{\partial Y}\right)-\frac{\partial}{\partial
Y}\left(D\frac{\partial V}{\partial X}\right)\right)=$
$\displaystyle\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\left(\delta(X-X_{A})\delta(Y-Y_{A})-\delta(X-X_{B})\delta(Y-Y_{B})\right).$
(9)

The new domain, specified by
$D\left(\sqrt{1+\frac{\Delta\sigma}{\bar{\sigma}}}X,\sqrt{1-\frac{\Delta\sigma}{\bar{\sigma}}}Y\right)$,
is given by
$\Theta(a^{2}-\left(1+\frac{\Delta\sigma}{\bar{\sigma}}\right)X^{2}-\left(1-\frac{\Delta\sigma}{\bar{\sigma}}\right)Y^{2})$,
i.e. it is an ellipse. If $\Delta\sigma/\bar{\sigma}>0$, the ellipse is
elongated along the $Y$-direction, if $\Delta\sigma/\bar{\sigma}<0$, then the
ellipse is elongated along the $X$-direction. Without loss of generality we
can choose the $x$-axis to be along the axis with larger resistivity, i.e. it
will be assumed from now on that

$\Delta\sigma/\bar{\sigma}<0.$ (10)

The equation (9) can be expressed using complex coordinates

$\displaystyle Z=X+iY,$ (11)

when, after some simplification, it becomes

$\displaystyle-2\left(\frac{\partial}{\partial Z}\left(D\frac{\partial
V}{\partial\bar{Z}}\right)+\frac{\partial}{\partial\bar{Z}}\left(D\frac{\partial
V}{\partial
Z}\right)\right)-\frac{2i\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\left(\frac{\partial}{\partial\bar{Z}}\left(D\frac{\partial
V}{\partial Z}\right)-\frac{\partial}{\partial Z}\left(D\frac{\partial
V}{\partial\bar{Z}}\right)\right)=$
$\displaystyle\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\left(\delta(X-X_{A})\delta(Y-Y_{A})-\delta(X-X_{B})\delta(Y-Y_{B})\right).$
(12)

To avoid confusion, the right-hand-side is kept in terms of the real and
imaginary parts of $Z$. This form makes it clear that inside the ellipse where
$D=1$, the solution can be written in terms of a sum of a function of $Z$ and
a function of $\bar{Z}$. The boundary conditions are determined from the right
hand side and the derivatives of the boundary function $D$.

###  II.1 Zhukovsky conformal mapping of the ellipse to annulus

It will be convenient to perform a conformal map transforming the boundary of
the ellipse to the boundary of the circle. This can be done using the
Zhukovsky transformation

$\displaystyle Z$ $\displaystyle=$
$\displaystyle\alpha_{+}w+\frac{\alpha_{-}}{w},$ (13) $\displaystyle w$
$\displaystyle=$ $\displaystyle u+iv,$ (14)

where $u(X,Y)$ and $v(X,Y)$ are purely real. To determine the coefficients
$\alpha_{+}$ and $\alpha_{-}$ we demand that

$\displaystyle\left(1+\frac{\Delta\sigma}{\bar{\sigma}}\right)X_{0}^{2}+\left(1-\frac{\Delta\sigma}{\bar{\sigma}}\right)Y_{0}^{2}$
$\displaystyle=$ $\displaystyle a^{2},$ (15)

implies

$\displaystyle u_{0}^{2}+v_{0}^{2}=1,$ (16)

i.e. if $X_{0}$ and $Y_{0}$ lie on the ellipse, then $u_{0}$ and $v_{0}$ are
forced to lie on the unit circle. From Eq.(13), we have

$\displaystyle X_{0}+iY_{0}$ $\displaystyle=$
$\displaystyle\alpha_{+}(u_{0}+iv_{0})+\alpha_{-}(u_{0}-iv_{0}),$ (17)

because, being on unit circle, $1/(u_{0}+iv_{0})=u_{0}-iv_{0}$. Therefore,

$\displaystyle X_{0}$ $\displaystyle=$
$\displaystyle(\alpha_{+}+\alpha_{-})u_{0},$ (18) $\displaystyle Y_{0}$
$\displaystyle=$ $\displaystyle(\alpha_{+}-\alpha_{-})v_{0}.$ (19)

So, from (15)

$\displaystyle\left(1+\frac{\Delta\sigma}{\bar{\sigma}}\right)(\alpha_{+}+\alpha_{-})^{2}u^{2}_{0}+\left(1-\frac{\Delta\sigma}{\bar{\sigma}}\right)(\alpha_{+}-\alpha_{-})^{2}v^{2}_{0}=a^{2},$

which implies

$\displaystyle\alpha_{\pm}=\frac{a}{2}\left(\frac{1}{\sqrt{1+\frac{\Delta\sigma}{\bar{\sigma}}}}\pm\frac{1}{\sqrt{1-\frac{\Delta\sigma}{\bar{\sigma}}}}\right).$
(21)

This fixes the conformal map. Having established that the ellipse in the
$(X,Y)$-plane maps onto the unit circle in the $(u,v)$-plane, we wish to know
where does the interior of the ellipse map. To this end, seek such $w=\Omega$
that would give

$\displaystyle\alpha_{+}\Omega$ $\displaystyle=$
$\displaystyle\frac{\alpha_{-}}{\Omega}\in\Re e,$ (22)

for $\Delta\sigma/\bar{\sigma}<0$. This gives

$\displaystyle\Omega$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\sqrt{1-\frac{\Delta\sigma}{\bar{\sigma}}}-\sqrt{1+\frac{\Delta\sigma}{\bar{\sigma}}}}{\sqrt{1-\frac{\Delta\sigma}{\bar{\sigma}}}+\sqrt{1+\frac{\Delta\sigma}{\bar{\sigma}}}}}.$
(23)

So, letting $w=\Omega e^{i\phi}$ where $\phi$ is the polar angle in the
$u,v$-plane and using (22) results in

$\alpha_{+}\Omega e^{i\phi}+\frac{\alpha_{-}}{\Omega
e^{i\phi}}=a\frac{\sqrt{-2\frac{\Delta\sigma}{\bar{\sigma}}}}{\sqrt{1-\left(\frac{\Delta\sigma}{\bar{\sigma}}\right)^{2}}}\cos\phi.$
(24)

This means that the circle of radius $\Omega$ in $u,v$-plane maps onto the
line segment connecting the foci of the ellipse in the $X,Y$-plane. For
$\Delta\sigma/\bar{\sigma}<0$, the foci lie on the x-axis. Therefore, the
ellipse in $X,Y$-plane, including its interior, maps onto an annulus in the
$u,v$-plane with the outer radius $1$ and the inner radius $\Omega$ as
illustrated in the Figure 2.

(a) (b)

Figure 2: For $\Delta\sigma/\bar{\sigma}<0$, (a) the interior of the circular
device with radius a maps onto the interior of the ellipse after the
coordinate rescaling. (b) The interior of the ellipse in the $X,Y$-plane
conformally maps onto the annulus in the $u,v$-plane with unit outer radius
and inner radius set by $\Omega$.

Because $Z=f(w)$ i.e. $Z$ is a function of $w$, $w$ is in turn a function of
$Z$, i.e. $w=g(Z)$, the left hand side of the differential equation can be
written as

$\displaystyle-2\left(\frac{\partial}{\partial Z}\left(D\frac{\partial
V}{\partial\bar{Z}}\right)+\frac{\partial}{\partial\bar{Z}}\left(D\frac{\partial
V}{\partial
Z}\right)\right)-\frac{2i\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\left(\frac{\partial}{\partial\bar{Z}}\left(D\frac{\partial
V}{\partial Z}\right)-\frac{\partial}{\partial Z}\left(D\frac{\partial
V}{\partial\bar{Z}}\right)\right)=$ $\displaystyle\frac{\partial w}{\partial
Z}\frac{\partial\bar{w}}{\partial\bar{Z}}\left(-\left(\frac{\partial}{\partial
u}D\frac{\partial V}{\partial u}+\frac{\partial}{\partial v}D\frac{\partial
V}{\partial
v}\right)+\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\left(\frac{\partial
D}{\partial v}\frac{\partial V}{\partial u}-\frac{\partial D}{\partial
u}\frac{\partial V}{\partial v}\right)\right).$ (25)

Using Cauchy-Riemann conditions, it can be readily shown that

$\displaystyle\frac{\partial w}{\partial
Z}\frac{\partial\bar{w}}{\partial\bar{Z}}=J\left(\frac{u,v}{X,Y}\right),$ (26)

where $J\left(\frac{u,v}{X,Y}\right)$ is the Jacobian determinant. Eq.(12)
therefore gives

$-\left(\frac{\partial}{\partial u}D\frac{\partial V}{\partial
u}+\frac{\partial}{\partial v}D\frac{\partial V}{\partial
v}\right)+\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\left(\frac{\partial
D}{\partial v}\frac{\partial V}{\partial u}-\frac{\partial D}{\partial
u}\frac{\partial V}{\partial
v}\right)=\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{\delta(X-X_{A})\delta(Y-Y_{A})-\delta(X-X_{B})\delta(Y-Y_{B})}{J\left(\frac{u,v}{X,Y}\right)}.$
(27)

But, by the properties of the Dirac $\delta$ function under coordinate
transformation, it follows that

$\displaystyle-\left(\frac{\partial}{\partial u}D\frac{\partial V}{\partial
u}+\frac{\partial}{\partial v}D\frac{\partial V}{\partial
v}\right)+\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\left(\frac{\partial
D}{\partial v}\frac{\partial V}{\partial u}-\frac{\partial D}{\partial
u}\frac{\partial V}{\partial v}\right)=$
$\displaystyle\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\left(\delta(u-u_{A})\delta(v-v_{A})-\delta(u-u_{B})\delta(v-v_{B})\right),$
(28)

where $D=\Theta\left(1-u^{2}-v^{2}\right)$. Now, because $X_{A,B},Y_{A,B}$ lie
on the ellipse, $u_{A,B},v_{A,B}$ must lie on the unit circle. From Eq.17

$\displaystyle\frac{x_{A,B}}{\sqrt{1+\frac{\Delta\sigma}{\bar{\sigma}}}}$
$\displaystyle=$
$\displaystyle\left(\alpha_{+}+\alpha_{-}\right)u_{A,B}\Rightarrow
x_{A,B}=au_{A,B},$ (29)
$\displaystyle\frac{y_{A,B}}{\sqrt{1-\frac{\Delta\sigma}{\bar{\sigma}}}}$
$\displaystyle=$
$\displaystyle\left(\alpha_{+}-\alpha_{-}\right)v_{A,B}\Rightarrow
y_{A,B}=av_{A,B}.$ (30)

###  II.2 Polar coordinates in the $u,v$-plane

Switching to polar coordinates in the $u,v$-plane

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle\sqrt{u^{2}+v^{2}},$ (31)
$\displaystyle\phi$ $\displaystyle=$ $\displaystyle\tan^{-1}\frac{v}{u},$ (32)

gives

$\displaystyle\frac{\partial}{\partial u}$ $\displaystyle=$
$\displaystyle\cos\phi\frac{\partial}{\partial\rho}-\frac{\sin\phi}{\rho}\frac{\partial}{\partial\phi},$
(33) $\displaystyle\frac{\partial}{\partial v}$ $\displaystyle=$
$\displaystyle\sin\phi\frac{\partial}{\partial\rho}+\frac{\cos\phi}{\rho}\frac{\partial}{\partial\phi}.$
(34)

Therefore, the derivatives of the boundary function are

$\displaystyle\frac{\partial D}{\partial u}$ $\displaystyle=$
$\displaystyle-\cos\phi\delta(\rho-1),$ (35) $\displaystyle\frac{\partial
D}{\partial v}$ $\displaystyle=$ $\displaystyle-\sin\phi\delta(\rho-1),$ (36)

and the differential equation (II.1) becomes

$\displaystyle-\left(\frac{\partial}{\partial\rho}\left(D\rho\frac{\partial
V}{\partial\rho}\right)+\frac{D}{\rho}\frac{\partial^{2}V}{\partial\phi}\right)+\frac{\sigma_{H}\delta(\rho-1)}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{\partial
V}{\partial\phi}=\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\delta(\rho-1)\left(\delta(\phi-\theta_{A})-\delta(\phi-\theta_{B})\right).$
(37)

###  II.3 Homogeneous solution and the boundary conditions

A general solution of Eq.(37) for $\rho<1$ where the terms containing
$\delta(\rho-1)$ vanish can be written as

$\displaystyle V(\rho,\phi)$ $\displaystyle=$
$\displaystyle\sum_{m=1}^{\infty}\left(A_{|m|}\left(\frac{\rho^{m}}{\Omega^{m}}+\frac{\Omega^{m}}{\rho^{m}}\right)\cos
m\phi+B_{|m|}\left(\frac{\rho^{m}}{\Omega^{m}}-\frac{\Omega^{m}}{\rho^{m}}\right)\sin
m\phi\right).$ (38)

This form satisfies the homogeneous differential equation and is continuous
and differentiable across the line cut joining the foci. To see this, notice
that the points on the circle of radius $\Omega$ in the $u,v$-plane map onto
the line segment joining the foci
$X\in(-2\alpha_{+}\Omega,2\alpha_{+}\Omega)$, $Y=0$, as we saw in the Eq.24.
Therefore, the points on the inner circle in the $u,v$ plane which are related
by the mirror reflection about the $v=0$ axis should be identified as the same
points. In other words, $\rho=\Omega$ and $\phi$, and $\rho=\Omega$ and
$-\phi$ map onto the same physical point in the $X,Y$ and therefore $x,y$
plane. We therefore want the potential at $\Omega^{+}$ and $\phi$ to either be
the same at $-\phi$ which is accomplished by
$\left(\frac{\rho^{m}}{\Omega^{m}}+\frac{\Omega^{m}}{\rho^{m}}\right)\cos
m\phi$, or we want it to vanish at $\Omega^{+}$ with a continuous slope.
Vanishing at $\Omega$ is accomplished by
$\frac{\rho^{m}}{\Omega^{m}}-\frac{\Omega^{m}}{\rho^{m}}$, and the reason why
only $\sin m\phi$ can multiply it is that multiplying it by $\cos m\phi$ would
introduce a cusp across the line segment.

Integrating both sides of Eq.37 over an infinitesimal interval straddling
$\rho=1$ gives the boundary condition

$\displaystyle\frac{\partial
V}{\partial\rho}|_{\rho=1}+\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{\partial
V}{\partial\phi}|_{\rho=1}=\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\left(\delta(\phi-\theta_{A})-\delta(\phi-\theta_{B})\right).$
(39)

Substituting Eq.38 into the above results in

$\displaystyle\frac{\partial V(\rho,\phi)}{\partial\rho}|_{\rho=1}$
$\displaystyle=$
$\displaystyle\sum_{m=1}^{\infty}\left(mA_{|m|}\left(\frac{1}{\Omega^{m}}-\Omega^{m}\right)\cos
m\phi+mB_{|m|}\left(\frac{1}{\Omega^{m}}+\Omega^{m}\right)\sin m\phi\right),$
(40) $\displaystyle\frac{\partial V(\rho,\phi)}{\partial\phi}|_{\rho=1}$
$\displaystyle=$
$\displaystyle\sum_{m=1}^{\infty}\left(-mA_{|m|}\left(\frac{1}{\Omega^{m}}+\Omega^{m}\right)\sin
m\phi+mB_{|m|}\left(\frac{1}{\Omega^{m}}-\Omega^{m}\right)\cos m\phi\right),$
(41)

and the differential equation (37) becomes

$\displaystyle\sum_{m=1}^{\infty}\left(mA_{|m|}\left(\frac{1}{\Omega^{m}}-\Omega^{m}\right)\cos
m\phi+mB_{|m|}\left(\frac{1}{\Omega^{m}}+\Omega^{m}\right)\sin m\phi\right)+$
(42) $\displaystyle+$
$\displaystyle\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\sum_{m=1}^{\infty}\left(-mA_{|m|}\left(\frac{1}{\Omega^{m}}+\Omega^{m}\right)\sin
m\phi+mB_{|m|}\left(\frac{1}{\Omega^{m}}-\Omega^{m}\right)\cos m\phi\right)=$
$\displaystyle\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{1}{\pi}\sum_{m=1}^{\infty}\left(\cos
m\phi\left(\cos m\theta_{A}-\cos m\theta_{B}\right)+\sin m\phi\left(\sin
m\theta_{A}-\sin m\theta_{B}\right)\right),$

where the following identity was used for the right hand side

$\displaystyle\delta\left(\phi-\theta_{A,B}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\sum_{m=-\infty}^{\infty}e^{im\phi}e^{-im\theta_{A,B}}$
(43) $\displaystyle=$
$\displaystyle\frac{1}{2\pi}+\frac{1}{\pi}\sum_{m=1}^{\infty}\left(\cos
m\phi\cos m\theta_{A,B}+\sin m\phi\sin m\theta_{A,B}\right).$ (44)

Matching the coefficients of $\cos m\phi$ and $\sin m\phi$ and solving for
$A_{|m|}$ and $B_{|m|}$ gives

$\displaystyle
A_{|m|}=\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{1}{\pi}\frac{1}{1+\frac{\sigma^{2}_{H}}{\sigma_{+}\sigma_{-}}}\frac{1}{m}\left(\frac{\cos
m\theta_{A}-\cos
m\theta_{B}}{\Omega^{-m}-\Omega^{m}}-\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{\sin
m\theta_{A}-\sin m\theta_{B}}{\Omega^{-m}+\Omega^{m}}\right)$ (45)
$\displaystyle
B_{|m|}=\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{1}{\pi}\frac{1}{1+\frac{\sigma^{2}_{H}}{\sigma_{+}\sigma_{-}}}\frac{1}{m}\left(\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{\cos
m\theta_{A}-\cos m\theta_{B}}{\Omega^{-m}-\Omega^{m}}+\frac{\sin
m\theta_{A}-\sin m\theta_{B}}{\Omega^{-m}+\Omega^{m}}\right).$ (46)

Thus,

$\displaystyle
V(\rho,\phi)=\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{1}{\pi}\frac{1}{1+\frac{\sigma^{2}_{H}}{\sigma_{+}\sigma_{-}}}\sum_{m=1}^{\infty}\frac{1}{m}\left(\frac{\cos
m\theta_{A}-\cos
m\theta_{B}}{\Omega^{-m}-\Omega^{m}}-\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{\sin
m\theta_{A}-\sin
m\theta_{B}}{\Omega^{-m}+\Omega^{m}}\right)\left(\frac{\rho^{m}}{\Omega^{m}}+\frac{\Omega^{m}}{\rho^{m}}\right)\cos
m\phi$ (47) $\displaystyle+$
$\displaystyle\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{1}{\pi}\frac{1}{1+\frac{\sigma^{2}_{H}}{\sigma_{+}\sigma_{-}}}\sum_{m=1}^{\infty}\frac{1}{m}\left(\frac{\sin
m\theta_{A}-\sin
m\theta_{B}}{\Omega^{-m}+\Omega^{m}}+\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{\cos
m\theta_{A}-\cos
m\theta_{B}}{\Omega^{-m}-\Omega^{m}}\right)\left(\frac{\rho^{m}}{\Omega^{m}}-\frac{\Omega^{m}}{\rho^{m}}\right)\sin
m\phi$ $\displaystyle-$ $\displaystyle(A\rightarrow B).$

###  II.4 Summing over the angular momenta

The sum over $m$ converges slowly. In order to convert it into a rapidly
convergent sum, we first Taylor expand the denominators involving $\Omega^{m}$
and $\Omega^{-m}$, in powers of $\Omega$ as

$\displaystyle
V(\rho,\phi)=\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{1}{\pi}\frac{1}{1+\frac{\sigma^{2}_{H}}{\sigma_{+}\sigma_{-}}}\sum_{n=0}^{\infty}\sum_{m=1}^{\infty}\frac{\left(\Omega^{1+2n}\right)^{m}}{m}\left(\cos
m\theta_{A}-(-1)^{n}\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\sin
m\theta_{A}\right)\left(\frac{\rho^{m}}{\Omega^{m}}+\frac{\Omega^{m}}{\rho^{m}}\right)\cos
m\phi$ (48) $\displaystyle+$
$\displaystyle\frac{I}{\sqrt{\sigma_{+}\sigma_{-}}}\frac{1}{\pi}\frac{1}{1+\frac{\sigma^{2}_{H}}{\sigma_{+}\sigma_{-}}}\sum_{n=0}^{\infty}\sum_{m=1}^{\infty}\frac{\left(\Omega^{1+2n}\right)^{m}}{m}\left((-1)^{n}\sin
m\theta_{A}+\frac{\sigma_{H}}{\sqrt{\sigma_{+}\sigma_{-}}}\cos
m\theta_{A}\right)\left(\frac{\rho^{m}}{\Omega^{m}}-\frac{\Omega^{m}}{\rho^{m}}\right)\sin
m\phi$ $\displaystyle-$ $\displaystyle(A\rightarrow B).$

Then the resulting sum over $m$ is related to the geometric series by
integration, and since $\Omega<1$, the sum over $n$ will be rapidly
convergent. Therefore, for $C=A,B$, we have

$\displaystyle\sum_{m=1}^{\infty}e^{im\theta_{C}}e^{im\phi}\frac{\Omega^{m(1+2n)}}{m}\left(\frac{\rho^{m}}{\Omega^{m}}\pm\frac{\Omega^{m}}{\rho^{m}}\right)$
$\displaystyle=$
$\displaystyle-\ln\left(1-e^{i\theta_{C}}\frac{w}{\Omega}\Omega^{1+2n}\right)\mp\ln\left(1-e^{i\theta_{C}}\frac{\Omega}{\bar{w}}\Omega^{1+2n}\right).$
(49)

Adding longitudinal and Hall contributions finally gives Eq.(2).

##  III Discussion

Because the differential equation is linear, it is straightforward to
generalize the expression derived above to the case with multiple sources and
drains. In such a case the continuity equation reads

$\nabla\cdot\mathbf{j}=\sum_{j=1}^{n_{S}}I^{S}_{j}\delta(\mathbf{r}-\mathbf{r}_{A,j})-\sum_{j=1}^{n_{D}}I^{D}_{j}\delta(\mathbf{r}-\mathbf{r}_{B,j}),$
(50)

where $n_{S}$ is the number of point sources and $n_{D}$ is the number of
point drains, and $\sum_{j=1}^{n_{S}}I^{S}_{j}=\sum_{j=1}^{n_{D}}I^{D}_{j}=I$.
The resulting expression is

$\displaystyle V(x,y;\\{\mathbf{r}_{A,j}\\},\\{\mathbf{r}_{B,j}\\})=$ (51)
$\displaystyle\frac{1}{\pi}\frac{\sqrt{\sigma_{+}\sigma_{-}}}{\sigma_{+}\sigma_{-}+\sigma^{2}_{H}}\sum_{j=1}^{n_{D}}I^{D}_{j}\left(\sum_{n=0,2,4,\ldots}^{\infty}\ln\left|1+e^{-2i\theta_{B,j}}\Omega^{2+4n}-e^{-i\theta_{B,j}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right|+\sum_{n=1,3,5,\ldots}^{\infty}\ln\left|1+e^{2i\theta_{B,j}}\Omega^{2+4n}-e^{i\theta_{B,j}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right|\right)$
$\displaystyle-$
$\displaystyle\frac{1}{\pi}\frac{\sqrt{\sigma_{+}\sigma_{-}}}{\sigma_{+}\sigma_{-}+\sigma^{2}_{H}}\sum_{j=1}^{n_{S}}I^{S}_{j}\left(\sum_{n=0,2,4,\ldots}^{\infty}\ln\left|1+e^{-2i\theta_{A,j}}\Omega^{2+4n}-e^{-i\theta_{A,j}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right|+\sum_{n=1,3,5,\ldots}^{\infty}\ln\left|1+e^{2i\theta_{A,j}}\Omega^{2+4n}-e^{i\theta_{A,j}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right|\right)$
$\displaystyle+$
$\displaystyle\frac{1}{\pi}\frac{\sigma_{H}}{\sigma_{+}\sigma_{-}+\sigma^{2}_{H}}\left(\sum_{n=0,2,4,\ldots}^{\infty}\sum_{j=1}^{n_{D}}I^{D}_{j}\arg\left(1+e^{-2i\theta_{B,j}}\Omega^{2+4n}-e^{-i\theta_{B,j}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right)-\sum_{j=1}^{n_{S}}I^{S}_{j}\arg\left(1+e^{-2i\theta_{A,j}}\Omega^{2+4n}-e^{-i\theta_{A,j}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right)\right.$
$\displaystyle+$
$\displaystyle\left.\sum_{n=1,3,5,\ldots}^{\infty}\sum_{j=1}^{n_{D}}I^{D}_{j}\arg\left(1+e^{2i\theta_{B,j}}\Omega^{2+4n}-e^{i\theta_{B,j}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right)-\sum_{j=1}^{n_{S}}I^{S}_{j}\arg\left(1+e^{2i\theta_{A,j}}\Omega^{2+4n}-e^{i\theta_{A,j}}\frac{Z}{\alpha_{+}}\Omega^{2n}\right)\right).$

The expression for extended source/drain can be found by treating
$I^{S/D}_{j}$ as infinitesimal and then converting the Riemann sum into an
integral. The obtained expression can now be used to fit measurements with
multiple voltage probes for arbitrary current source and drain placed on the
perimeter of the disk.

###### Acknowledgements.

I wish to express sincere gratitude to Prof. J.I.A. Li for sharing their
unpublished results and for his encouragement to publish this article. I would
also like to thank Prof. Jian Kang for going over the calculations in the
manuscript. O. V. is supported by NSF DMR-1916958 and is partially funded by
the Gordon and Betty Moore Foundation’s EPiQS Initiative Grant GBMF11070,
National High Magnetic Field Laboratory through NSF Grant No. DMR-1157490 and
the State of Florida.

## References

  * [1] L. Balents, C.R. Dean, D.K. Efetov, and A.F. Young “Superconductivity and strong correlations in moiré flat bands,” Nature Physics 16, 725 (2020).
  * [2] K. Kang, T. Li, E. Sohn, J. Shan, and Kin Fai Mak, “Nonlinear anomalous Hall effect in few-layer WTe2” Nature Materials 18, 324 (2019).
  * [3] J.I.A. Li et.al. (unpublished).
  * [4] X. Wang et.al. arXiv:2209.xxxxx.
  * [5] I. Lifshits, M. Azbel, and M. Kaganov, “Electron Theory of Metals” (Springer, 1973).


# Power efficient ReLU design for neuromorphic computing using spin Hall
effect

Venkatesh Vadde1,1, Bhaskaran Muralidharan1,2, Abhishek Sharma2,3
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>1Department of Electrical Engineering, Indian Institute of Technology Bombay,
Powai, Mumbai-400076, India  2Department of Electrical Engineering, Indian
Institute of Technology Ropar, Rupnagar, Punjab-140001, India

###### Abstract

We demonstrate a magnetic tunnel junction injected with spin Hall current to
exhibit linear rotation of magnetization of the free-ferromagnet using only
the spin current. Using the linear resistance change of the MTJ, we devise a
circuit for the rectified linear activation (ReLU) function of the artificial
neuron. We explore the role of different spin Hall effect (SHE) heavy metal
layers on the power consumption of the ReLU circuit. We benchmark the power
consumption of the ReLU circuit with different SHE layers by defining a new
parameter called the spin Hall power factor. It combines the spin Hall angle,
resistivity, and thickness of the heavy metal layer, which translates to the
power consumption of the different SHE layers during spin-orbit
switching/rotation of the free FM. We employ a hybrid spintronics-CMOS
simulation framework that couples Keldysh non-equilibrium Green’s function
formalism with Landau-Lifshitz-Gilbert-Slonzewski equations and the HSPICE
circuit simulator to account for diverse physics of spin-transport and the
CMOS elements in our proposed ReLU design. We also demonstrate the robustness
of the proposed ReLU circuit against thermal noise and non-trivial power-error
trade-off that enables the use of an unstable free-ferromagnet for energy-
efficient design. Using the proposed circuit, we evaluate the performance of
the convolutional neural network for MNIST datasets and demonstrate comparable
classification accuracies to the ideal ReLU with an energy consumption of 75
$pJ$ per sample.

##  I Introduction

Artificial neural networks (ANNs) are widely used by machine learning and data
science communities to solve complex problems. The ANNs are inspired by the
biological brains, which have memory and computing intertwined to solve
diverse problems while consuming low energy [1, 2]. The von Neumann-based
modern computers that separate memory and computing are not suitable for
hardware implementation of neural networks [3]. Neural networks contain highly
interconnected perceptrons, which define the mathematical model of biological
neurons as the sum of weighted inputs passed through an activation function
[4].

Learning in neural networks can be achieved by using activation functions [5].
The activation function introduces non-linearity to the network and enables
the network to learn complex data structures and differentiate between
outputs. Traditionally, sigmoid and tanh activation functions have been widely
utilized. But these standard functions limit the network’s ability to learn
since they saturate when the input is very high or very low[6, 7]. The sigmoid
and tanh functions also face the vanishing gradient problem, where the
gradient information used to learn networks becomes almost zero for deep
networks, thus affecting the deep network’s learning capacity[6]. Glorot et
al.[7] showed that the rectified linear unit (ReLU) activation function can
improve the learning speed of various neural networks. The ReLU function also
overcomes the vanishing gradient and saturation problems that tanh and sigmoid
functions face. It often produces better results than traditional functions
such as tanh and sigmoid in neural networks [7, 6]. Thus, the ReLU has become
a default activation function for various neural networks [8, 9, 5].

The ReLU function is described as

$f(x)=\begin{cases}0&\text{if $x\leq 0$}\\\ x&\text{if $x>0$}\end{cases}$

The activation function using CMOS technology has been explored by a few works
[10, 11, 12], but these realizations of the activation function require
additional interconnect circuits to interface with synaptic and max-pooling
layers[13] of the neural network. The CMOS implementation is also limited by
the area and energy requirements [3, 14]. On the other hand, spintronics
provides a wide range of devices that can be engineered to have non-
volatility, plasticity, oscillatory and stochastic behavior[15, 16, 17, 18,
19, 20]. These properties suit well for in-memory computing, enabling
neuromorphic computing and taking advantage of the paradigm “let physics do
the computing”[21]. In this paper, we demonstrate an MTJ-based design to
emulate the ReLU function, which can be easily connected to cross-bar-based
synaptic layers[22, 23] and the max-pooling layers[13].

Current induced spin-orbit torques (SOTs) originating from spin Hall
effect[24, 25, 26, 27] in heavy metal(HM)/ferromagnet(FM) hetero-structure,
have recently emerged as an energy-efficient manipulation of magnetization at
the nanoscale. The efficient conversion of the charge current to the spin
current via SHE is quantified by the spin Hall angle($\theta_{SH}$). There
have been consistent efforts [28, 29, 30, 31, 32, 33, 34] to increase
$\theta_{SH}$, but it came at the cost of increased resistivity, which
increases the power consumption of the HM. We define the spin Hall power
parameter($\epsilon_{SHE}$) that directly relates to the power consumption of
the free FM switching process, accounting for the $\theta_{SH}$, $\rho$, and
the thickness(t) of the HM. The $\epsilon_{SHE}$ defined in this work can be
used to compare the power consumption of different SHE layers for the free FM
switching/rotation via spin-orbit torque.

This paper is organized as follows. In Sec. II, we give the details on the
design of the ReLU circuit and introduce our spin Hall power factor. In Sec.
III, we describe our simulation platform where we couple Keldysh non-
equilibrium Green’s function formalism with Landau-Lifshitz-Gilbert-Slonzewski
equations and the HSPICE circuit simulator to account for diverse physics of
spin-transport and the CMOS elements in our proposed ReLU design. In Sec. IV,
we present the results of our ReLU design and the performance of the proposed
circuit against the thermal stability factor. Here we also show the result of
our investigation into various heavy metals. We show that our design is
resistant to thermal noise and that there exists a non-trivial power-error
trade-off that leads to the energy-efficient circuit design using unstable
free-ferromagnets. In Sec. V we explore the potential application of our ReLU
design using convolutional neural networks and show that our network achieves
practically the same classification accuracy as ideal ReLU implementation. We
conclude in Sec. VI.

##  II Design

(a)

(b)

Figure 1: Design schematics. (a) The MTJ device is stacked on top of the heavy
metal layer. Charge current is injected into the HM layer along $\hat{x}$,
which injects $\hat{y}$-polarized spin current into the free-FM of the MTJ.
(b). Circuit design for ReLU functionality. The current source $I_{bias}$
converts the change in resistance to a change in voltage that is connected to
the CMOS inverter.

###  II-A ReLU circuit

MTJs are traditionally used as binary memories, here we show that an MTJ can
be designed to have linear functionality. This is achieved when a spin current
whose polarization is orthogonal to the anisotropy direction is applied to the
free-FM. In this work, we inject a $\hat{y}-$ polarized spin current to the
perpendicular magnetic anisotropy (PMA)-FM(CoFeB) in order to produce a linear
rotation in $\hat{x}$\- component of the magnetization as shown in Fig. 1a.

The linear rotation in magnetization is translated to the change in resistance
via the TMR effect of the MTJ. The injected current($I_{bias}$) renders the
resistance change to the voltage change across the MTJ, which drives the CMOS
inverter to obtain the ReLU functionality, as shown in Fig. 1b. The CMOS
inverter operates in the linear region to invert and amplifies the voltage
change across the MTJ, which results in ReLU function output as shown in Fig.
3b. The current source $I_{bias}$ can be replaced with a resistor to obtain
the ReLU functionality at the cost of decreased output ($V_{out}$) swing.

###  II-B Spin Hall power factor

The SHE-driven MTJs are being explored extensively for low-energy applications
using spintronic devices. Charge-to-spin conversion via spin Hall effect
exhibited by the heavy metal is utilized for switching the free-FM layer in
SHE-driven MTJs[35, 29]. The power consumption of the switching can be
minimized using a large charge-to-spin conversion factor (Eq. 2). There have
been consistent efforts [28, 29, 30, 31, 32, 33, 34] to increase the spin Hall
angle ($\theta$) through heavy metal engineering for enhanced charge-to-spin
conversion factor. But an increase in $\theta$ usually comes with an increase
in the resistivity of the heavy metal resulting in large power consumption.
Traditionally, the spin Hall conductivity[34] $\sigma_{SH}=\theta/\rho$ has
been used to characterize the SHE, although this includes the effect of
resistivity and spin Hall angle, it lacks analytical reasoning to represent
power consumption of the SHE layer. The spin Hall angle and resistivity also
depend on the thickness of heavy metal [32, 31], compelling us to define a
parameter that can unequivocally benchmark the various heavy metals for SHE
switching power consumption.

The charge to spin conversion[24, 25, 36] of the SHE layer and the
polarization and direction[24, 25, 26, 27] of the generated spin current are
given by

$\theta_{SH}=\frac{J_{s}}{J_{c}}$ (1)

$I_{s}=\theta_{SH}\frac{L}{t}I_{c}$ (2)

$\hat{I_{s}}=\hat{I_{c}}\times\sigma$ (3)

Here, $I_{s}$ is the spin current generated, $\theta_{SH}$ is the spin Hall
angle of the heavy metal, L is the length of the heavy metal, t is the
thickness of the heavy metal, and $I_{c}$ is the charge current injected.
$\hat{I_{s}}$ is the direction of generated spin current flow, $\hat{I_{c}}$
is the direction of input charge current, and $\sigma$ is the polarization of
the generated spin current. From Eq. 3, injection of charge current to heavy
metal in $\hat{x}-$ direction results in y-polarized spin current injection to
the free-FM (z-direction) on top of the HM layer.

The resistance (R) of the heavy metal is given by

$R=\rho\frac{L}{Wt}$ (4)

Here, $\rho$ and W are the resistivity and width of the heavy metal
respectively.

The power consumed by the heavy metal is given by $P_{HM}=I_{c}^{2}R$. Here,
$I_{c}$ can be written as

$I_{c}=\frac{I_{s}}{\sqrt{RV/d}}\frac{\sqrt{\rho t}}{\theta_{SH}}$ (5)

where V, d are the volume and thickness of the free-FM layer. The power
consumed by the heavy metal is given by

$P_{HM}=\frac{I_{s}^{2}}{V/d}\frac{\rho t}{\theta_{SH}^{2}}$ (6)

Here the $I_{s}$ represents the spin current needed for switching the
ferromagnet. For a given free-FM layer $I_{s}$, V and d are constants, so from
the above derivation we define a spin Hall power factor
$\epsilon_{SHE}(=\frac{\sqrt{\rho t}}{\theta_{SH}})$ that can be used to
compare different heavy metal. The material with the lowest $\epsilon_{SHE}$
will consume less power.

In our proposed circuit, $440\mu A$ is required to achieve an output voltage
of $0.35V$ for free-FM with a thermal stability factor of 45. The proposed
factor is not limited to this work and can also be used in SHE-driven FM
switching mechanisms. Some HMs such as Pt [28] affect the damping factor
($\alpha$) of the free layer as well, leading to an increase in the spin
current ($I_{s}$) required. In such cases, the $\epsilon_{SHE}$ needs to be
multiplied by the change in the damping factor
($\frac{\alpha_{old}}{\alpha_{new}}$). In our proposed ReLU circuit, the
increased $\alpha$ has a negligible effect on the spin current requirement for
the linear rotation.

TABLE I:  Simulation Parameters. Symbol |  Quantity | Value
---|---|---
$M_{s}$ |  saturation magnetization | 1150 emu/$cm^{3}$
$H_{k}$ |  anisotropy field | 1070 - 3670 Oe[37]
V |  volume of ferromagnet | 1000 $nm^{3}$
d |  thickness of the ferromagnet | 1 nm
$\Delta$ |  thermal stability factor | 14.85 - 50.95
$\alpha$ |  Gilbert damping | 0.01
$C_{MTJ}$ |  MTJ capacitance | 26.56 aF
$I_{bias}$ |  bias current | 3.13 $\mu A$
$V_{DD}$ |  voltage source | 0.5 V
$C_{g}$ |  CMOS inverter input capacitance | 0.175 fF
$C_{o}$ |  CMOS inverter output capacitance | 0.305 fF
$\Delta t$ |  simulation time step | 0.5 pS
$\hbar$ |  reduced Plank’s constant |  $1.055\times 10^{-34}$ Js
$k_{B}$ |  Boltzmann constant | $1.38\times 10^{-16}ergK^{-1}$
T |  temperature | 300 K

##  III Simulation Methods

Figure 2 shows the schematic overview of the hybrid spintronics-CMOS
simulation framework. The MTJ and current source ($I_{bias}$) parameters are
given to the NEGF simulator as shown in Fig. 2. The NEGF simulator is self-
consistently coupled with $I_{bias}$, since the device resistance depends on
the MTJ angle and the voltage across the MTJ. The NEGF produces the device
resistance vs MTJ angle plot. This result is coupled to the HSPICE circuit
simulator via VerilogA. The circuit simulator simulates the entire circuit
including the magnetization dynamics[38, 39] and it also simulates the
approximation of the CMOS inverter pair based on a 16nm predictive technology
model (PTM)[40].

(a)

(b)

Figure 2: Hybrid NEGF-CMOS simulation platform setup. (a) The NEGF is self-
consistently coupled with $I_{bias}$ to calculate the resistance of the MTJ.
(b) The MTJ resistance is coupled to the HSPICE circuit simulator using
VerilogA. The LLGS equation is coupled with other parts of the circuit to
calculate free-FM magnetization.

###  III-A Quantum Transport

We use the Keldysh NEGF technique [41, 17, 42, 43] to describe the transport
through MTJ that has MgO sandwiched between free and fixed CoFeB FM layers.
The NEGF formalism is given by

$\displaystyle G(E)=[EI-H-\Sigma]^{-1},$ (7) $\displaystyle
A(E)=i[G-G^{\dagger}],$ (8)
$\displaystyle\Gamma_{T,B}(E)=i([\Sigma_{T,B}(E)]-[\Sigma_{T,B}(E)]^{\dagger}),$
(9)
$\displaystyle\Sigma^{in}(E)=[\Gamma_{T}(E)]f_{T}(E)+[\Gamma_{B}(E)]f_{B}(E),$
(10) $\displaystyle G^{n}=\int dE[G(E)][\Sigma^{in}(E)][G(E)]^{\dagger},$ (11)
$\displaystyle\Sigma=\Sigma_{T}+\Sigma_{B}$ (12)

Here $[H]$ is the device Hamiltonian, $[H]=[H_{0}]+[U]$, comprising device
tight-binding matrix $[H_{0}]$ and the Coulomb charging matrix $[U]$, and
$[I]$ is the identity matrix, $E$ is the energy variable. The charging matrix
$[U]$ is calculated self-consistently using Poisson’s equation. $G(E)$ is the
Green’s function matrix, $\Gamma_{T,B},f_{T,B},\Sigma_{T,B}$ are the
broadening matrix, the Fermi function, and the self-energy matrices for the
top (fixed) and bottom (free) FM layers respectively. $A$ is the spectral
function, $\Sigma^{in}$ is the in-scattering function, and $G^{n}$ is the
correlation matrix.

The quantum transport segment culminates with the calculation of the current
operator ($I_{op}$) that represents the charge current between two lattice
points i and i+1 is given by

$I_{op}=\frac{i}{\hbar}(H_{i,i+1}G^{n}_{i+1,i}-H_{i+1,i}G^{n}_{i,i+1})$ (13)

The current operator $I_{op}$ is $2\times 2$ matrix in the spin space of the
lattice point. Using this the charge current can be evaluated as

$I=q\int Real[Trace(\hat{I}_{op})]dE,$ (14)

where q is the quantum of electronic charge.

(a)

(b)

Figure 3: (a) The free-FM magnetization($m_{x}$) and the MTJ resistance as the
normalized input ($I_{in}/I_{0}$) is varied with an initial magnetization
oriented in the +z direction. $I_{in}$ is applied along $\hat{x}$. (b) The
output of the ReLU circuit with varied normalized input. $I_{0}=380\mu A$ and
$\Delta=45$ with $Au_{0.25}Pt_{0.75}$ heavy metal.

###  III-B Current injected MTJ

In the ReLU circuit, we have employed the current source $I_{bias}$ to
translate the change in MTJ resistance to voltage variation. The MTJ
resistance depends not just on the free-FM magnetization but also on the
voltage across MTJ [13]. The MTJ voltage itself depends on the resistance and
the current source $I_{bias}$, so the device resistance and voltage need to be
found self consistently. Figure 2(a) shows the algorithm for this self-
consistent calculation of the resistance while accounting for its voltage and
current source $I_{bias}$ dependence. The self-consistent loop is run until
the MTJ current is equal to the applied $I_{bias}$.

###  III-C Magnetization dynamics

The LLGS equation [44, 45] is used to describe the magnetization dynamics of
the free-FM. The LLGS equation is given by

$(\frac{1+\alpha^{2}}{\gamma
H_{k}})\frac{d\hat{m}}{dt}=-\hat{m}\times\vec{h}_{eff}-\alpha\hat{m}\times\hat{m}\times\vec{h}_{eff}\\\
-\hat{m}\times\hat{m}\times\vec{i}_{s}+\alpha\hat{m}\times\vec{i}_{s},$ (15)

where $\hat{m}$ is the unit vector along the direction of magnetization of the
free magnet, $\gamma$ is the gyromagnetic ratio, $\alpha$ is the Gilbert
damping parameter, $\vec{h}_{eff}=\frac{\vec{H}_{eff}}{H_{k}}$ is the reduced
effective field and $\vec{i}_{s}=\frac{\hbar\vec{I}_{s}}{2qM_{s}VH_{k}}$ is
the normalized spin current. The term $\vec{H}_{eff}$ includes the
contribution of the anisotropy field ($H_{k}$) and the thermal noise
($H_{th}$). The thermal noise is given by $\langle
H_{th}^{2}\rangle=\frac{2\alpha k_{B}T}{\gamma M_{s}V}$ and $\langle\rangle$
represents the ensemble average [46].

##  IV Results

TABLE II: Parameters for different HMs. Heavy Metal |  Spin Hall angle $\theta$ |  Resistivity $\rho$ ($\mu\Omega cm$) |  HM thickness t (nm) |  HM Length L (nm) |  HM Width W (nm) |  Spin Hall power factor $\epsilon_{SHE}=\frac{\sqrt{\rho t}}{\theta}$ ($\sqrt{\Omega}nm$) |  $I_{0}$ ($\mu A$)
---|---|---|---|---|---|---|---
$Pt$ [28] |  0.07 |  20 |  6 |  38.73 |  25.82 |  494.87 |  980
$\beta-Ta$ [29] |  0.12 |  190 |  4 |  10.26 |  97.46 |  726.48 |  1442
$\alpha+\beta-W$ [30] |  0.18 |  80 |  6.2 |  19.68 |  50.80 |  391.26 |  775
$Au_{0.25}Pt_{0.75}$ [31] |  0.3 |  83 |  4 |  15.52 |  64.42 |  192.06 |  380
$\beta-W$ [32] |  0.35 |  210 |  9 |  14.64 |  68.31 |  392.79 |  778
$W(O)$ [33] |  0.5 |  200 |  6 |  12.25 |  81.65 |  219.09 |  434
$W_{0.88}Ta_{0.12}$ [34] |  0.58 |  260 |  5 |  9.81 |  101.98 |  196.58 |  389

We show in Fig. 3a the linear rotation in magnetization of the free FM layer
of the MTJ. The linear rotation is achieved by injecting $\hat{y}-$polarized
spin current into the free FM layer. The spin current is generated by applying
charge current to the HM layer along $\hat{x}$. The TMR effect of the MTJ
translates the linear magnetization change into a linear change in the
resistance, this is shown in Fig. 3a. The resulting linear variation of MTJ
resistance is employed in the circuit design (Fig. 1) to realize the ReLU
output as shown in Fig. 3b. The output closely emulates the ReLU activation
function for normalized inputs of less than 1. The parameters used in this
design are given in Tab. I.

We evaluate the role of different HMs in our proposed ReLU design. The ReLU
circuit’s performance is assessed against the thermal stability factor
($\Delta=\frac{H_{k}M_{s}V}{2k_{B}T}$) of the free FM layer. The $\Delta$
factor not only captures the stability of the free FM against thermal noise
but also determines the spin current required for MTJ switching. We vary the
$\Delta$ factor of the free FM by changing the anisotropy field while keeping
the dimensions of the free FM fixed. A decrease in the $\Delta$ factor reduces
the spin current needed for linear rotation, diminishing the HM’s input charge
current and power consumption.

(a)

(b)

Figure 4: (a) Power consumption of the ReLU circuit with different HMs with
varying $\Delta$. With y-axis in log scale (b) Power consumption of the entire
ReLU circuit and power consumption of only the HM layer using
$Au_{0.25}Pt_{0.75}$

The heavy metals used in this paper for the ReLU design are given in Tab. II
along with respective parameters such as spin Hall angle $\theta$, resistivity
$\rho$, and thickness $t$ taken from experimental works [28, 29, 30, 31, 32,
33, 34]. Using these parameters, the length $L$ and width $W$ are calculated
such that the HM offers a resistance of 50 $\Omega$. The spin Hall power
factor ($\epsilon_{SHE}$) and the normalizing current $I_{0}$ (current
required to achieve maximum magnetization rotation with the ReLU output of
0.35V) are also shown in the table. It is inferred from Tab. II that the
current $I_{0}$ accompanies $\epsilon_{SHE}$. We can see from Fig. 4a and Tab.
II, that the average static power consumption of the ReLU circuit for
different HMs also trails the $\epsilon_{SHE}$. The static power consumption
decreases with $\Delta$ as it decimates the spin current required for the ReLU
functionality. We show in Fig. 4b the contribution of HM power consumption to
that of the entire ReLU circuit for $Au_{25}Pt_{75}$ as the SHE layer. Along
with the HM, MTJ consumes a fixed amount of power for translating the
magnetization changes to voltage changes hence not affected by the change in
$\Delta$. The MTJ sensing power consumption dominates the total ReLU power at
lower $\Delta$s.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 5: Static power consumption of the ReLU circuit averaged over the
entire input range, and the average absolute percentage error for a normalized
input of 0.5 with (a) Pt, (b) $\beta-$Ta, (c) $\alpha+\beta-$W, (d)
$Au_{0.25}Pt_{0.75}$, (e) $\beta-$W, (f) W(O), and (g) $W_{0.88}Ta_{0.12}$ as
SHE layer. Figure 6: Schematic of CNN architecture used for training in
TensorFlow for the image classification task. The architecture has a $3\times
3$ convolution layer followed by the ReLU activation function, max-pooling
layers, and finally the fully connected layer.

We show in Fig. 5 the static power consumption and the average absolute
percentage error of the ReLU circuit for different HMs. The average of
absolute error increases as $\Delta$ is decreased since the contribution of
the thermal noise ($\langle H_{th}^{2}\rangle=\frac{2\alpha k_{B}T}{\gamma
M_{s}V}$) increases in the effective magnetic field ($H_{eff}$). The $H_{eff}$
includes the effect from both thermal noise and the anisotropy field. As the
$\Delta$ is reduced, the $\langle H_{th}^{2}\rangle$ stays constant, but the
$H_{k}$ decreases. The effect of the thermal noise on the output is estimated
using 100 Monte Carlo [47] simulations, with a normalized input
($I_{in}/I_{0}$) of 0.5. All the HMs show a decrease in power consumption and
an increase in the error while the $\Delta$ is reduced, presenting an
opportunity to optimize the circuit to consume less power while obtaining
reliable results. For $Au_{0.25}Pt_{0.75}$ the static power consumption is
1.37 $\mu W$ while the absolute error percentage is $2.98\%$ for $\Delta=20$
as shown in Fig. 5d. These results suggest that an unstable ($\Delta<40$) free
FM-based MTJ can be used to obtain reliable results for the ReLU circuit.

##  V Application: Convolutional Neural Networks

Convolutional Neural Networks (CNN) are a class of artificial neural networks
that produces excellent results in problems involving image data. Figure 6
shows the architecture of the CNN used for classifying the MNIST and fashion
MNIST data sets. Here we use our developed ReLU circuit instead of the ideal
ReLU in the feature extraction stage of the CNN to train the network.

(a)

(b)

Figure 7: Classification accuracy(%) and loss of the CNN as a function of
epochs for different datasets: (a) MNIST (b) Fashion MNIST.

We show in Fig 7a $\&$ 7b the accuracy and loss of the network during training
and testing for MNIST and fashion MNIST datasets respectively. We notice an
accuracy of $98.76\%$ on the test data for the MNIST dataset and an accuracy
of $91.3\%$ for the fashion MNIST dataset. The accuracies for full software
implementation(with ideal ReLU) are $98.86\%$ and $90.41\%$ for MNIST and
fashion MNIST datasets. The accuracies using the non-ideal ReLU (our developed
ReLU circuit) are closer to the ideal ReLU, this demonstrates the robustness
of our ReLU circuit for neuromorphic applications.

Figure 8: Energy consumption of the ReLU implementation in CNN architecture
for testing a single sample for different heavy metals and the spin Hall power
factors of the heavy metals.

Figure 8 shows the energy consumption for ReLU implementation in our CNN
architecture for different SHE layers, along with their spin Hall power
factors. Here the free-FMs have a thermal stability factor of $20$. Here we
see that the energy consumption of the ReLU implementation follows the spin
Hall power factor. Our results suggest that $75pJ$ of energy is consumed by
the ReLU implementation in testing a single sample in our CNN architecture
with $Au_{0.25}Pt_{0.75}$ as heavy metal and a free-FM with $\Delta=20$.

##  VI Conclusion

In this paper, we showcased the linear rotation in the magnetization of free-
FM and proposed a circuit design that effectively emulates the ReLU function,
a fundamental component of deep learning neural networks. We introduced a new
metric, the spin Hall power factor, to unequivocally quantify the SHE layers’
power consumption. Our simulation results not only confirm the validity of
this factor but also demonstrate its potential to significantly impact the
design of SHE driven devices and circuits. We deploy our developed simulation
framework that combines current injected MTJ with NEGF, LLG, and HSPICE
circuit simulator, enabling us to design and analyze the proposed ReLU circuit
with varying HMs.

We demonstrated the existence of a non-trivial power error trade-off that
enables the use of unstable free-FM for energy-efficient ReLU design. We
demonstrated that the most energy-efficient realization of the proposed ReLU
circuit consumes 1.37 $\mu W$ of static power with a low error rate of 2.98%
using the HM $Au_{0.25}Pt_{0.75}$ at $\Delta=20$. Furthermore, we showed the
potential of our ReLU design in CNNs, producing classification accuracies
close to the software ReLU implementation with an energy consumption of 75
$pJ$ per sample.

## Acknowledgements

The author AS acknowledges the support of ISIRD phase-1 project of IIT Ropar.
The author BM wishes to acknowledge the Science and Engineering Board (SERB),
Government of India, for funding under the MATRICS grant (Grant No.
MTR/2021/000388).

## Conflict of Interest

The authors have no conflicts to disclose.

## Data availability

Data is available on request from the authors.

## References

  * [1] D. Marković, A. Mizrahi, D. Querlioz, and J. Grollier, “Physics for neuromorphic computing,” _Nature Reviews Physics_ , vol. 2, no. 9, pp. 499–510, 2020.
  * [2] C. D. Schuman, T. E. Potok, R. M. Patton, J. D. Birdwell, M. E. Dean, G. S. Rose, and J. S. Plank, “A survey of neuromorphic computing and neural networks in hardware,” _arXiv preprint arXiv:1705.06963_ , 2017.
  * [3] J. Grollier, D. Querlioz, K. Camsari, K. Everschor-Sitte, S. Fukami, and M. D. Stiles, “Neuromorphic spintronics,” _Nature electronics_ , vol. 3, no. 7, pp. 360–370, 2020.
  * [4] H. Taud and J. Mas, “Multilayer perceptron (mlp),” in _Geomatic approaches for modeling land change scenarios_. Springer, 2018, pp. 451–455.
  * [5] H. Ide and T. Kurita, “Improvement of learning for cnn with relu activation by sparse regularization,” in _2017 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2017, pp. 2684–2691.
  * [6] I. Goodfellow, Y. Bengio, and A. Courville, _Deep learning_. MIT press, 2016.
  * [7] X. Glorot, A. Bordes, and Y. Bengio, “Deep sparse rectifier neural networks,” in _Proceedings of the fourteenth international conference on artificial intelligence and statistics_. JMLR Workshop and Conference Proceedings, 2011, pp. 315–323.
  * [8] V. Nair and G. E. Hinton, “Rectified linear units improve restricted boltzmann machines,” in _Icml_ , 2010.
  * [9] K. Jarrett, K. Kavukcuoglu, M. Ranzato, and Y. LeCun, “What is the best multi-stage architecture for object recognition?” in _2009 IEEE 12th international conference on computer vision_. IEEE, 2009, pp. 2146–2153.
  * [10] C.-H. Chang, H.-Y. Kao, and S.-H. Huang, “Hardware implementation for multiple activation functions,” in _2019 IEEE International Conference on Consumer Electronics-Taiwan (ICCE-TW)_. IEEE, 2019, pp. 1–2.
  * [11] C. Geng, Q. Sun, and S. Nakatake, “An analog cmos implementation for multi-layer perceptron with relu activation,” in _2020 9th International conference on modern circuits and systems technologies (MOCAST)_. IEEE, 2020, pp. 1–6.
  * [12] P. Priyanka, G. Nisarga, and S. Raghuram, “Cmos implementations of rectified linear activation function,” in _International Symposium on VLSI Design and Test_. Springer, 2018, pp. 121–129.
  * [13] V. Vadde, B. Muralidharan, and A. Sharma, “Orthogonal spin current injected magnetic tunnel junction for convolutional neural networks,” _arXiv preprint arXiv:2207.14603_ , 2022.
  * [14] G. Indiveri, B. Linares-Barranco, T. J. Hamilton, A. v. Schaik, R. Etienne-Cummings, T. Delbruck, S.-C. Liu, P. Dudek, P. Häfliger, S. Renaud _et al._ , “Neuromorphic silicon neuron circuits,” _Frontiers in neuroscience_ , vol. 5, p. 73, 2011.
  * [15] S. Datta, _Lessons from Nanoelectronics: A New Perspective on Transport—Part B: Quantum Transport_. World Scientific, 2018.
  * [16] D. Apalkov, B. Dieny, and J. M. Slaughter, “Magnetoresistive random access memory,” _Proceedings of the IEEE_ , vol. 104, no. 10, pp. 1796–1830, 2016\.
  * [17] A. Sharma, A. A. Tulapurkar, and B. Muralidharan, “Resonant spin-transfer-torque nano-oscillators,” _Physical Review Applied_ , vol. 8, no. 6, p. 064014, 2017.
  * [18] K. Y. Camsari, S. Salahuddin, and S. Datta, “Implementing p-bits with embedded mtj,” _IEEE Electron Device Letters_ , vol. 38, no. 12, pp. 1767–1770, 2017\.
  * [19] A. Hirohata, K. Yamada, Y. Nakatani, I.-L. Prejbeanu, B. Diény, P. Pirro, and B. Hillebrands, “Review on spintronics: Principles and device applications,” _Journal of Magnetism and Magnetic Materials_ , vol. 509, p. 166711, 2020.
  * [20] A. Sengupta and K. Roy, “Neuromorphic computing enabled by physics of electron spins: Prospects and perspectives,” _Applied Physics Express_ , vol. 11, no. 3, p. 030101, 2018.
  * [21] A. Parihar, N. Shukla, M. Jerry, S. Datta, and A. Raychowdhury, “Computing with dynamical systems based on insulator-metal-transition oscillators,” _Nanophotonics_ , vol. 6, no. 3, pp. 601–611, 2017.
  * [22] S. H. Jo, T. Chang, I. Ebong, B. B. Bhadviya, P. Mazumder, and W. Lu, “Nanoscale memristor device as synapse in neuromorphic systems,” _Nano letters_ , vol. 10, no. 4, pp. 1297–1301, 2010.
  * [23] D. Fan, M. Sharad, and K. Roy, “Design and synthesis of ultralow energy spin-memristor threshold logic,” _IEEE Transactions on Nanotechnology_ , vol. 13, no. 3, pp. 574–583, 2014.
  * [24] C. Song, R. Zhang, L. Liao, Y. Zhou, X. Zhou, R. Chen, Y. You, X. Chen, and F. Pan, “Spin-orbit torques: Materials, mechanisms, performances, and potential applications,” _Progress in Materials Science_ , vol. 118, p. 100761, 2021.
  * [25] L. Liu, O. Lee, T. Gudmundsen, D. Ralph, and R. Buhrman, “Current-induced switching of perpendicularly magnetized magnetic layers using spin torque from the spin hall effect,” _Physical review letters_ , vol. 109, no. 9, p. 096602, 2012.
  * [26] S. Takahashi and S. Maekawa, “Spin current, spin accumulation and spin hall effect,” _Science and Technology of Advanced Materials_ , vol. 9, no. 1, p. 014105, 2008.
  * [27] J. Hirsch, “Spin hall effect,” _Physical review letters_ , vol. 83, no. 9, p. 1834, 1999.
  * [28] L. Liu, T. Moriyama, D. Ralph, and R. Buhrman, “Spin-torque ferromagnetic resonance induced by the spin hall effect,” _Physical review letters_ , vol. 106, no. 3, p. 036601, 2011.
  * [29] 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_ , vol. 336, no. 6081, pp. 555–558, 2012.
  * [30] C.-F. Pai, L. Liu, Y. Li, H. Tseng, D. Ralph, and R. Buhrman, “Spin transfer torque devices utilizing the giant spin hall effect of tungsten,” _Applied Physics Letters_ , vol. 101, no. 12, p. 122404, 2012.
  * [31] L. Zhu, D. C. Ralph, and R. A. Buhrman, “Highly efficient spin-current generation by the spin hall effect in au 1- x pt x,” _Physical Review Applied_ , vol. 10, no. 3, p. 031001, 2018.
  * [32] Q. Hao and G. Xiao, “Giant spin hall effect and switching induced by spin-transfer torque in a w/co 40 fe 40 b 20/mgo structure with perpendicular magnetic anisotropy,” _Physical Review Applied_ , vol. 3, no. 3, p. 034009, 2015.
  * [33] K.-U. Demasius, T. Phung, W. Zhang, B. P. Hughes, S.-H. Yang, A. Kellock, W. Han, A. Pushp, and S. S. Parkin, “Enhanced spin–orbit torques by oxygen incorporation in tungsten films,” _Nature communications_ , vol. 7, no. 1, pp. 1–7, 2016.
  * [34] N. Behera, H. Fulara, L. Bainsla, A. Kumar, M. Zahedinejad, A. Houshang, and J. Åkerman, “Energy-efficient w 100- x ta x/co-fe-b/mgo spin hall nano-oscillators,” _Physical Review Applied_ , vol. 18, no. 2, p. 024017, 2022.
  * [35] 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_ , vol. 476, no. 7359, pp. 189–193, 2011.
  * [36] A. Sengupta, C. M. Liyanagedera, B. Jung, and K. Roy, “Magnetic tunnel junction as an on-chip temperature sensor,” _Scientific reports_ , vol. 7, no. 1, pp. 1–8, 2017.
  * [37] M. Gajek, J. Nowak, J. Sun, P. Trouilloud, E. O’sullivan, D. Abraham, M. Gaidis, G. Hu, S. Brown, Y. Zhu _et al._ , “Spin torque switching of 20 nm magnetic tunnel junctions with perpendicular anisotropy,” _Applied Physics Letters_ , vol. 100, no. 13, p. 132408, 2012.
  * [38] G. D. Panagopoulos, C. Augustine, and K. Roy, “Physics-based spice-compatible compact model for simulating hybrid mtj/cmos circuits,” _IEEE Transactions on Electron Devices_ , vol. 60, no. 9, pp. 2808–2814, 2013.
  * [39] J. Z. Sun, “Spin-current interaction with a monodomain magnetic body: A model study,” _Physical Review B_ , vol. 62, no. 1, p. 570, 2000.
  * [40] “Predictive technology model (ptm),” https://ptm.asu.edu/, (Accessed on 09/12/2022).
  * [41] K. Y. Camsari, S. Chowdhury, and S. Datta, “The non-equilibrium green function (negf) method,” _arXiv preprint arXiv:2008.01275_ , 2020.
  * [42] D. Datta, “Modeling of spin transport in mtj devices,” Ph.D. dissertation, Purdue University, 2012.
  * [43] D. Datta, B. Behin-Aein, S. Datta, and S. Salahuddin, “Voltage asymmetry of spin-transfer torques,” _IEEE Transactions on Nanotechnology_ , vol. 11, no. 2, pp. 261–272, 2011.
  * [44] J. C. Slonczewski, “Current-driven excitation of magnetic multilayers,” _Journal of Magnetism and Magnetic Materials_ , vol. 159, no. 1-2, pp. L1–L7, 1996.
  * [45] A. Brataas, A. D. Kent, and H. Ohno, “Current-induced torques in magnetic materials,” _Nature materials_ , vol. 11, no. 5, pp. 372–381, 2012.
  * [46] J. Z. Sun, T. Kuan, J. Katine, and R. H. Koch, “Spin angular momentum transfer in a current-perpendicular spin-valve nanomagnet,” in _Quantum Sensing and Nanophotonic Devices_ , vol. 5359. SPIE, 2004, pp. 445–455.
  * [47] J. Hammersley, _Monte carlo methods_. Springer Science & Business Media, 2013.