# On the Challenges of using Reinforcement Learning in Precision Drug Dosing:
Delay and Prolongedness of Action Effects

Sumana Basu1,2, Marc-André Legault1,2, Adriana Romero-Soriano1,2,3, Doina
Precup1,2

###### Abstract

Drug dosing is an important application of AI, which can be formulated as a
Reinforcement Learning (RL) problem. In this paper, we identify two major
challenges of using RL for drug dosing: delayed and prolonged effects of
administering medications, which break the Markov assumption of the RL
framework. We focus on prolongedness and define PAE-POMDP (Prolonged Action
Effect-Partially Observable Markov Decision Process), a subclass of POMDPs in
which the Markov assumption does not hold specifically due to prolonged
effects of actions. Motivated by the pharmacology literature, we propose a
simple and effective approach to converting drug dosing PAE-POMDPs into MDPs,
enabling the use of the existing RL algorithms to solve such problems.

We validate the proposed approach on a toy task, and a challenging glucose
control task, for which we devise a clinically-inspired reward function. Our
results demonstrate that: (1) the proposed method to restore the Markov
assumption leads to significant improvements over a vanilla baseline; (2) the
approach is competitive with recurrent policies which may inherently capture
the prolonged effect of actions; (3) it is remarkably more time and memory
efficient than the recurrent baseline and hence more suitable for real-time
dosing control systems; and (4) it exhibits favourable qualitative behavior in
our policy analysis.

##  1 Introduction

Drug dosing plays an important role in human health– e.g.  individuals with
type 1 diabetes require regular insulin injections to manage their blood
glucose levels, intensive care patients require continuous monitoring and
administration of drugs, optimal doses of anaesthesia are required during
operative procedures, etc. Optimal drug dosing is most important in cases
where the therapeutic window is narrow, meaning that small deviations from the
therapeutic range of drug concentration may lead to serious clinical
complications (FDA 2015; Maxfield and Zineh 2021). These problems are
compounded by idiosyncratic differences in the dynamics of drug absorption,
distribution, metabolism or excretion (collectively referred to as
pharmacokinetics) and drug sensitivity. Therefore, one important goal of
precision medicine is to tailor patient care while accounting for individual
characteristics, and developing algorithmic solutions for drug dosing is a
contribution towards that broader goal.

Reinforcement learning (RL) offers a framework to account for individual
characteristics and automatically derive personalized treatment policies in
line with the objective of precision medicine (Ribba et al. 2020). However,
RL-based algorithms cannot be applied off-the-shelf to tackle precision dosing
since all drugs are known to have a _delayed_ and _prolonged_ effect from the
point of medication (Holford 2018) (see Figure 1). The delay is attributed to
the time it takes for the drug to distribute to the target site, bind to the
receptor and finally to change physiological substances _before_ its response
can be observed. This can vary between minutes to hours or even longer
(Holford 2018). The prolongedness is due to individual variation in
pharmacokinetics (Vogenberg, Barash, and Pursel 2010). In this delayed and
prolonged action effects scenario, the future depends on the previous drug
dosages and their effect, and therefore the Markov assumption usually made by
RL algorithms no longer holds. Although RL has been applied to address drug
dosing problems such as controlling glucose levels for closed loop artificial
pancreas (Tejedor, Woldaregay, and Godtliebsen 2020), the violation of the
Markov assumption is in this case not only problematic from an RL research
perspective, but also from a safety perspective, as ignoring the delayed and
prolonged effects of a drug can lead to drug overdosing related toxicity
(Guengerich 2011).

Figure 1: Delayed and prolonged effects of drugs. Drug is only administered at
the first step, but no immediate effect is observed. After the initial delay,
the initial dose keeps altering the patient status for a while, even in the
absence of additional dosage. The duration of delay and prolongedness is
individual specific.

##### Contributions.

In this paper, we identify prolongedness and delay as fundamental roadblocks
to using RL in precision drug dosing, and focus on addressing the former. To
that end, we introduce the _prolonged action effect partially observable
Markov decision process_ (PAE-POMDP), a framework for modeling delayed action
effects in decision making. We then present assumptions inspired by the
pharmacology literature to convert the _non-Markovian prolonged action effect_
problem into a Markov decision process (MDP), therefore enabling the use of
recent advances in model-free RL. To the best of our knowledge, our work is
the _first_ to explore the prolonged effect of actions. We validate the
proposed approach on a toy task, where the only violation of the Markov
property comes from the prolonged effect of actions, and show that restoring
the Markov assumption allows the RL agent we develop to significantly
outperform previous baselines. We also address the challenging task of glucose
control by optimal insulin dosing for type 1 diabetes patients by leveraging
the open-sourced version of the FDA-approved UDA/Padova simulator (Dalla, MD,
and C. 2009; Xie 2018). Although the glucose control task has garnered
interest in the RL literature in the past (Tejedor, Woldaregay, and
Godtliebsen 2020), there does not appear to be a widely adopted reward
function. Therefore, we design a clinically motivated reward function that
explicitly avoids overdosing while effectively controlling glucose levels. Our
results show that our approach of converting the PAE-POMDP into a MDP is not
only competitive in terms of performance, but is also remarkably more time and
memory efficient than the baselines, making it more suitable for real-time
dosing control systems. With these contributions, we aim to raise awareness of
both the delayed and prolonged effects of drugs while tackling personalized
drug dosing with RL, and hope that the proposed reward function will help
break the entry barrier to control glucose levels for a closed-loop artificial
pancreas, while fostering future research in this direction.

##### What this paper does not do.

This paper does not aim to provide a general solution to tackle the prolonged
effect of actions, but one crafted specifically for precision drug dosing that
allows to quickly bring the problem back to the MDP framework in a compute-
and memory-efficient way. We do not provide a new algorithm to tackle
prolongedness, but instead propose a pharmacologically motivated effective way
to enable the use of already existing RL algorithms for precision drug dosing.
We are assuming drug action in isolation, instead of in combination with other
drugs. Learning and predicting drug synergies is a growing body of research
that can provide interesting future work.

##  2 Related Work

Blood glucose control in individuals with Type 1 diabetes is a longstanding
problem in medicine that has stimulated the interest of RL researchers for a
long time. In a systematic review of reinforcement learning approaches for
blood glucose management, Tejedor, Woldaregay, and Godtliebsen (2020) reported
347 papers between 1990 and 2019 on the topic, of which 11 used the UVA/Padova
simulator (Dalla, MD, and C. 2009; Dalla et al. 2014; Xie 2018). While most of
these work explore RL algorithms – e.g.  Actor Critic, Q-learning, SARSA, and
DQN among others – to address the glucose control problem, none of them
acknowledges the fundamental challenge of prolongedness and most consider the
current blood glucose level as a sufficient statistic for state information.
Even in the cases where actions are considered as part of the state
information, only the most recent action is taken into account, which is not
enough to restore the Markov assumption violated due to prolongedness (see
appendix). Fox et al. (2020) formalized the problem as a POMDP by augmenting 4
hours of blood glucose and insulin information, but they did not recognize
prolongeness of drug effects as the reason for partial observability. In fact,
to handle prolongedness, the most recent blood glucose measurement and
residual active insulin are sufficient, as we show in the experiments section.

There is a similar theme of exploring various state of the art deep RL
algorithms for other drug dosing problems. However, these works also miss the
crucial point of recognizing and handling the prolonged effect of drugs. In a
non-exhaustive list, Nemati, Ghassemi, and Clifford (2016) use a Q-Network to
learn a personalized heparin dosing policy, Weng et al. (2017) use policy
iteration to manage glycemia for septic patients in critical care, Lin et al.
(2018) propose the use of Deep Deterministic Policy Gradient (DDPG, (Lillicrap
et al. 2016)) for heparin dosing in critical care patients, and Lopez-Martinez
et al. (2019) use Double Deep Q-Networks (DDQN) to administer opioids for pain
management in critical care patients. But none of these papers discusses the
prolonged effect of drug doses. Zadeh, Street, and Thomas (2022) use Deep
Q-Networks (DQN) to administer the anticoagulant warfarin, and they consider a
pre-defined duration along with a dose. They use a
Pharmacokinetic/Pharmacodynamic (PK/PD) model to determine the duration and
add it to the state information, along with the patient information, blood
coagulability measure and dosing history. However, the authors do not
explicitly recognize prolonged drug effect as one of the reasons for adding
the duration information. In addition, since dose response is individual
specific, finding the right duration is a challenge in itself.

##  3 Background

##### Markov Decision Process (MDP).

Reinforcement Learning (RL) is a framework for solving sequential decision
making problems where an agent interacts with the environment and receives
feedback in the form of reward. The typical formal framework for RL is the
Markov Decision Process (MDP). A MDP $\mathcal{M}$ is a 5-tuple
$(\mathcal{S},\mathcal{A},r,\mathcal{P},\gamma)$, where $\mathcal{S}$ is a
(finite) set of states, $\mathcal{A}$ is a (finite) set of actions,
$\mathcal{P}$ is the state transition probability
$\mathcal{P}(s_{t+1}=s^{\prime}|s_{t}=s,a_{t}=a)$,
$r:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow\mathbb{R}$ is the
reward function and $\gamma\in[0,1)$ is the discount factor. As the name
suggests, a MDP obeys the Markov assumption that the future is in independent
of the past given the present, which means that transitions and rewards depend
only on the current state and action and not on the past history. The goal of
an RL agent is to find a policy
$\pi:\mathcal{S}\times\mathcal{A}\rightarrow[0,1]$ that maximizes the
cumulative discounted return $\sum_{t=0}^{\infty}\gamma^{t}r_{t}$.

##### POMDP and Recurrent Policy.

A Partially Observable MDP (POMDP) is a 7-tuple
$(\mathcal{S},\mathcal{A},r,\mathcal{P},\Omega,\mathcal{O},\gamma)$, where
$\Omega$ is the set of observations, $\mathcal{O}$ is the set of conditional
observation probabilities, $O(\omega|s)$, also known as the emission function,
and the rest of the elements are the same as in an MDP. In a POMDP, the agent
does not have direct access to the identity of the states, instead needing to
infer them through the observations. Note that while in a MDP, an agent which
aims to act optimally with respect to the expected long-term return only needs
to consider Markovian policies,
$\pi:\mathcal{S}\times\mathcal{A}\rightarrow[0,1]$, in a POMDP, policies need
to either rely on the entire history of action and observations, or to infer
the distribution of hidden states from this history (as in belief-based POMDP
dolution methods). In recent work, recurrent networks have become the standard
for implementing POMDP policies. We call recurrent policy a mapping
$\pi:\mathcal{T}\times\mathcal{A}\rightarrow[0,1]$ where $\mathcal{T}$ is the
space of trajectories $\tau=\\{(o_{t},a_{t},r_{t})\\}_{t=0}^{T}$ of up to $T$
time steps.

##### Q-Learning and Deep Q-Networks (DQN).

Q-Learning (Watkins and Dayan 1989) is a model-free RL algorithm that
estimates $\mathcal{Q}:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$, a
value function assessing the quality of action $a_{t}\in\mathcal{A}$ at state
$s_{t}\in\mathcal{S}$. The value function can be estimated as:
$Q(s_{t},a_{t})=Q(s_{t},a_{t})+\alpha\big{(}r_{t}+\gamma\max_{a}Q(s_{t+1},a)-Q(s_{t},a_{t})\big{)}$,
where $\alpha\in(0,1)$ is the learning rate. In DQN (Mnih, Kavukcuoglu, and et
al. 2015), the $Q$-values are estimated by a neural network minimizing the
Mean Square Bellman Error:
$L_{i}(\theta_{i})=\mathbb{E}_{s_{t},a_{t}\sim\pi_{b}}\big{[}(y_{i}-Q(s_{t},a_{t};\theta_{i}))^{2}\big{]}\\\
y_{i}=(r_{t}+\gamma\max_{a_{t+1}}Q(s_{t+1},a_{t+1};\theta_{i-1})$, where
$\pi_{b}$ is the behavior policy, $y_{i}$ is the target and $\theta_{i}$ are
the parameters of the $Q$-network at iteration $i$.

##  4 Method

###  4.1 Problem Formulation

Drug dosing can be formalized as a POMDP, because the effect of a medication
is felt over a period of time after its administration. In this section, we
introduce _prolonged action effect_ POMDPs, referred to as PAE-POMDPs, a
subclass of POMDPs in which an action’s effect lasts more than one time step.
More precisely, the action $a_{t}$ taken at time step $t$ continues to affect
the future states of the environment for $\kappa\in\mathbb{Z^{+}}$ time steps,
where $\kappa$ is the time interval necessary for $a_{t}$’s effect to fall
below a given threshold. The value of $\kappa$ is environment specific and
also depends on the amplitude of the action $a_{t}$. The above-defined
prolongedness offers a _forward view_ in time. From a _backward view_
perspective, the state at a time step $t$ is the result of super-imposing the
effects of the actions from several preceding time steps, which are still felt
at time $t$ i.e. $s_{t+1}=s_{t}+\sum_{k=0}^{\kappa}\Delta_{a_{t-k}}$, where
$\Delta_{a_{t}}$ is the effect of action $a_{t}$ on the state. In other words,
the Markov transition function $P(s_{t+1}|s_{t},a_{t})$ is no longer valid,
and instead the transition function becomes
$P(s_{t+1}|s_{t},(a_{t-k})_{k=0}^{\kappa})$. In this case, the reward function
can still be assumed to be Markov $R(s_{t},a_{t},s_{t+1})$ (see Equation 2).
From the perspective of the agent, this means that it needs to keep track of
the history of actions over a preceding period of time. In general, this
problem is no simpler than a regular POMDP, as an agent that keeps track of
its history may need to remember all the actions taken since the beginning of
time. However, in this work, we consider a more circumscribed problem
formulation which is relevant for therapeutic dosing. Specifically, we use
knowledge from pharmacology to enable us to restore the Markov assumptions for
drug-dosing PAE-POMDP.

###  4.2 Converting PAE-POMDP into MDP

In this section, we present the specific knowledge used to restore the Markov
assumptions in drug dosing problems and show how by leveraging them we can
effectively convert a PAE-POMDP into an MDP.

##### Action and effect equivalence.

We start by noting that as per the pharmacodynamics axioms, drug effects are
determined by drug concentration at the site of action (Holford 1984, 2018).
Therefore, action and effect are considered to be equivalent and used
interchangeably hereinafter.

##### Exponential rate of decay.

Motivated by the pharmacology (Benet and Zia-Amirhosseini 1995; Dasgupta and
Krasowski 2020), biotechnology (Hobbie and Roth 2007), as well as chemical-
kinetics (Peter and de Paula Julio 2006) literature, we adopt an exponential
decay model of drug concentration over time (Annamalai 2010). In particular,
we assume that the initial action effect will decay at a constant rate
$\lambda$, which is specific to the environment (drug as well as individual).
Formally, we assume that $a_{>t\leq\kappa}=\lambda a_{t}$. Note that here we
implicitly assume that the action’s effects can safely be ignored after
$\kappa$ time-steps. This assumption is valid in the context of drug dosing,
since drugs are ineffective below a certain concentration. The observation of
the agent at a particular time step can therefore be defined conditional on
the current state and on all the past actions:
$o_{t+1}\sim\mathcal{O}(.|s_{t},a_{<t})$. Note that although this is a
generalized notation, not all the actions $a_{<t}$ that occur before $t$ have
an effect on $o_{t+1}$ (only those within a $\kappa$ window). If the amplitude
of an action decays by a factor of $\lambda$ for each of the $\kappa$ time
steps during which the action’s effect is felt, then at any time-step $t$, the
environment state is a function of a sub-sequence of the previous actions
$\\{a_{\leq t}\\}$, which we call _effective action_ , $a_{t_{\mathcal{E}}}$,
and define as:

$a_{t_{\mathcal{E}}}=\begin{cases}0,&\text{if $t=0$}\\\ \lambda
a_{(t-1)_{\mathcal{E}}}+a_{t-1},&\text{otherwise}\\\ \end{cases}$ (1)

where $\lambda\in(0,1)$ enforces the exponential decay assumption on the
actions. Note that setting $\lambda=0$ falls back to the classical RL
scenario, whereas setting $\lambda=1$ entails a infinitely long action
effects.

##### Additive composition of action effects.

We note that Equation 1 leverages the additive property of drug concentration,
and therefore assumes that action effects are additive and independent of each
other. With the additive action effects assumption, the conditional
distribution $\mathcal{O}(.|s_{t},a_{<t})$ could be modelled as a sum of
effects due to the actions.

##### From PAE-POMDP to MDP.

Finally, we convert the PAE-POMDP into a MDP, by augmenting the states with
effective actions. Formally, we define a revised MDP,
$\mathcal{M_{E}}=(\mathcal{S_{E}},\mathcal{A},r,\mathcal{P},\gamma)$, where
the states in $\mathcal{S_{E}}$ are defined as
$s_{t_{\mathcal{E}}}=(s_{t},a_{t_{\mathcal{E}}})$. The updated state
definition restores the Markov assumption, and therefore $\mathcal{M_{E}}$ can
be solved as a proxy of the original MDP $\mathcal{M}$. Note that
$\mathcal{M_{E}}$ and $\mathcal{M}$ only differ in their state space. Since we
have access to $\mathcal{M_{E}}$, we can use traditional RL algorithms to
solve it.

##  5 Experiments

Figure 2: Performance on MoveBlock. Cumulative discounted return of tabular
Q-learning and Effective Q-learning. Mean and std of performance over 10
seeds. Performance is smoothed with a window of size 1000.

We validate the proposed approach on a toy environment, as a sanity check to
prototype our method, and on a challenging glucose control environment. This
section describes both environments and presents the obtained results.

###  5.1 Toy Environment: _MoveBlock_

##### Environment.

We created the MoveBlock environment to assess our proposed approach in a
controlled manner. In particular, we built the environment in such a way that
the only violation of the Markov property comes from the prolonged effect of
actions. In MoveBlock, the task is to move a block from its random initial
position to a fixed final position on a slippery surface with minimum possible
force applied along the horizontal axis in the direction of the goal. Due to
the slippery nature of the floor, an exerted force will cause the object to
move until it is stopped by friction, prolonging the effect of actions over
multiple time-steps. The reward system is the same as the one in the Mountain
Car environment (Brockman et al. 2016) – i.e.  there is a high reward for
reaching the goal position, but action magnitude proportional penalties to
encourage policies where minimum effort is used. The observation is the
position (continuous) while the velocity remains unobserved to the agent.

##### Experimental setup and results.

To keep the states tractable, we discretized the continuous position values
into 10k states. The action space is discrete, but to keep the effective
actions tractable, we clipped the effective action value to the maximum action
value. We then trained a tabular Q-learning agent, which serves as baseline,
with $\epsilon$-greedy exploration. Next, we modified this tabular Q-learning
agent with the method proposed in Section 4.2, and trained an _Effective
Q-learning_ agent. The results are shown in Figure 2. We observe that the
Effective Q-learning agent, which restores the Markov assumptions broken by
the prolonged action effects, consistently outperforms the Q-learning agent
despite using the same Q-learning algorithm, which highlights the importance
of considering the prolongedness of actions when solving PAE-POMDPs.

###  5.2 Glucose Control Environment

##### Environment.

We chose blood glucose control as a real-world example of PAE-POMDP given the
low therapeutic index of insulin, which makes it a good candidate for
precision dosing. We used Simglucose (Xie 2018), the OpenAI Gym (Brockman et
al. 2016) implementation of the FDA approved UVA/Padova Type1 Diabetes
Mellitus Simulator (T1DMS) (Dalla, MD, and C. 2009; Dalla et al. 2014) built
on an in-silico population of 30 patients. In this case, the task is to manage
blood glucose by administering insulin to type 1 diabetes patients, who lose
their ability to produce insulin for the rest of their lives. In particular,
the goal is to maintain blood glucose as close as possible to the non-diabetic
range for a given patient. However, the environment does not come with a
recognized pre-defined reward function.

##### Experimental setup.

We are interested in training a personalized glucose control policy instead of
a population-level policy. To this end, we randomly selected a virtual patient
from the simulator and modified the simulator for the experiments presented in
this paper (See Appendix for more details). In this PAE-POMDP setting, the
states are hidden patient parameters (Kovatchev et al. 2009), observations are
continuous blood glucose readings, actions are discrete insulin doses in the
range of $[0,5]$, transition dynamics are specific to the chosen patient but
unknown to the agent, and the reward function is described in the next
paragraph. We appended the current observations with the effective actions as
described in section 4.2.

##### Reward.

Reward design in RL is a challenging task in itself. The task of maintaining
blood glucose as close as possible to the non-diabetic range, can be addressed
by maximizing the time spent in the target zone with as little insulin as
possible. The minimal insulin dose requirement is due to the need to avoid
building insulin resistance, a condition where a person no longer responds to
small doses of insulin. Note that the tendency to develop resistance is
applicable to any drug. Unfortunately, recently introduced reward functions do
not fulfill the above-discussed criteria. Therefore, we designed a
biologically inspired custom zone reward that incentivizes the time spent in
the target zone and penalizes hyperglycemia, hypoglycemia and high insulin
doses. Formally, the reward is defined as:

$r_{t}(s_{t-1},a_{t},s_{t})=r_{t,state}(s_{t-1},s_{t})-r_{t,action}(a_{t})$
(2)

where $r_{t,state}$ is a state reward which encourages the agent to maintain
the target vital statistics within a healthy range as follows:

$\begin{split}r&{}_{t,state}(s_{t-1},s_{t})=\\\
&\begin{cases}-100&s_{t}<70\text{ or }s_{t}>200\text{ (episode
termination)}\\\ -1&s_{t}<100\text{ and }s_{t}-s_{t-1}<0.5\text{
(hypoglycemia)}\\\ -1&s_{t}>150\text{ and }s_{t}-s_{t-1}>0.5\text{
(hyperglycemia)}\\\ 10&100<=s_{t}<=150\text{ (target blood
glucose)}\end{cases}\end{split}$ (3)

and $r_{t,action}$ is a penalty term which encourages the agent to avoid
overdosing:

$r_{t,action}(a_{t})=0.1*a_{t}^{2}\quad\text{(penalize overdosing)}$ (4)

Note that this kind of reward function is well-suited for any homeostatic
control problem in drug dosing.

Figure 3: Performance (Cumulative discounted sum of rewards) on the glucose
control task. Since DQN assumes the environment to be Markov, its performance
is almost as poor as the fixed policy. The performance of Effective-DQN is at
par with ADRQN (history length 15); both consider prolonged effects of
actions. Figure 4: Glucose control performance breakdown. Average number of
steps in each blood glucose zone over training episodes (std over 5 seeds).
Due to the prolonged action effect, Effective-DQN spends more time in the
hyperglycemic zone as well as in the target zone than ADRQN. DQN learns a sub-
optimal policy that tends to overdose in order to to quickly get out of
hyperglycemia and eventually leads to episode termination due to hypoglycemia.

Model | Runtime | Max
---|---|---
|  | memory
ADRQN | $359.64\pm 19.25$ | 1623
Effective-DQN | $\textbf{129.04}\pm\textbf{42.73}$ | 1513

Table 1: Runtime and memory requirements. 5 agents trained on a Nvidia RTX
8000 GPU with the same 5 seeds (mean $\pm$ std). Runtime is in mins/10K
episode and memory in MiB.

Figure 5: Effect of history length on compute (left) and memory requirement
(right). Compute time increases exponentially and memory requirement linearly
with the history length for ADRQN agents, but Effective-DQN compute and memory
are independent of the history length and hence remain constant.

#### Policies under comparison.

All agents are trained using the zone reward function introduced in the
previous section.

  1. 1.

Fixed Dose. This is a handcrafted step-wise policy with the underlying
assumption that higher blood glucose levels call for higher doses.
Specifically, we inject the maximum insulin (5 units) for blood glucose in the
range of 190 to 200. For every 10 units drop in blood glucose level from 190,
we reduce the insulin dose by 1 unit until it is cut off at blood glucose 150.

  2. 2.

DQN. The DQN agent assumes the environment to be Markov, and hence is a good
baseline to capture the effect of ignoring the inherent prolongedness of the
actions.

  3. 3.

Action-specific Deep Recurrent Q-Network (ADRQN) (Zhu, Li, and Poupart 2017).
This baseline leverages a recurrent neural network that encodes the action
history and appends it with the current observation to form a state. The
action history length is set to 15, as per hyper-parameter tuning. Note that
since we are the first to acknowledge the prolonged effect of drug dosing, we
are also the first to propose an action history based baseline in this
context. In a scenario where compute and memory is not a limiting factor,
ADRQN (or any other action history based policy) can serve as the state of the
art.

  4. 4.

Effective-DQN. This approach enhances the DQN baseline, but with the decay
assumption introduced in the previous section i.e. the blood glucose values
are augmented with the effective insulin dose to represent the state.

##### Quantitative Results.

We compare the cumulative discounted sum of rewards for all the agents under
consideration in Figure 3. The fixed dose policy, which undergoes no learning,
achieves the lowest performance. Unsurprisingly, the DQN agent only slightly
improves upon this naive hand-crafted policy, since it assumes the environment
is Markov and thus ignores the prolonged effect of actions. The ADRQN agent
exhibits a significant performance boost as compared to the DQN agent, since
the learned policy depends not only on the current observation, but also on
the action history responsible for the prolonged effect. The Effective-DQN
agent’s performance is on par with the ADRQN agent, suggesting that our decay
assumption is a good approximation for the prolonged action effect observed in
the diabetes simulator (see Appendix for details). To gain further insights,
we stratify the performance of each RL agent into the hyperglycemic,
hypoglycemic, and target zones, and assess the number of steps each agent
spends in each blood glucose zone over training episodes (see Figure 4). We
observe that the DQN agent minimizes the negative reward for hyperglycemia and
spends minimum time in the hyperglycemic zone. However, due to the prolonged
effect of the actions taken in hyperglycemia, the agent spends very little
time in the target zone as well as hypoglycemic zone, suggesting that the
agent learns to overdose to quickly get out of hyperglycemia, and eventually
leading to episode termination due to hypoglycemia. Both ADRQN and Effective-
DQN spend a similar amount of time in the hypoglycemic zone, whereas
Effective-DQN spends slightly more time in hyperglycemia than ADRQN,
increasing the time it spends in the target zone. This explains the
performance gain of the Effective-DQN agent over ADRQN.

##### Efficiency analysis.

Since time sensitivity and limited memory are crucial aspects of real-time
control systems, we compare the best performing agents (ADRQN and Effective-
DQN) in terms of compute and memory efficiency. In Table 1, we present the
compute and maximum memory requirement to train the models on the glucose
control task. The _Runtime_ column presents the time required to train the
models for 10K episodes on a single Nvidia RTX 8000 GPU. We used the readily
available ADRQN PyTorch
implementation111https://mlpeschl.com/post/tiny_adrqn/, and ensured both ADRQN
and Effective-DQN are comparable in terms of capacity. Although some of the
differences might still be attributed to the implementation, many are due to
the individual properties of each model. In particular, ADRQN has the longest
runtime as it needs to learn the action history through the LSTM (Hochreiter
and Schmidhuber 1997) unit, and as shown in Figure 5 (left), the compute time
grows exponentially with the history length. History length depends on how
prolonged the action effects are, which in turn are individual and drug
specific factors. However, the compute time for Effective-DQN remains constant
no matter the history length, making the compute efficiency gap between ADRQN
and Effective-DQN increase for longer history lengths. In the _max memory_
column of Table 1, we present the maximum memory required by the models during
training. We keep a replay buffer of the same size to train all the models. As
shown in the table, ADRQN requires slightly more memory than Effective-DQN.
While this difference does not appear significant, it is worth noting that
ADRQN requires memory proportional to the history length (see Figure 5
(right)). Moreover, in order to benefit from longer histories, ADRQN might
also require larger models, further increasing its memory requirement. In
contrast, the Effective-DQN’s memory requirement remains constant as we
increase the history length.

##### Qualitative results.

(a) Fixed Dose

(b) DQN

(c) ADRQN

(d) Effective-DQN

Figure 6: Policy visualization. Trained greedy agents evaluated on 10 episodes
with the same starting blood glucose. Red crosses mark the end of episodes.
The insulin administered is as prescribed by the respective policy. In the
blood glucose plot, the space between the two red dotted lines is the target
zone. Effective-DQN appears to be conservative and thus does not allow the
patient to get into hypoglycemia. The ADRQN agent gets into both hypo- and
hyperglycemia.

We present a qualitative assessment of the policies under comparison to better
understand their behavior. Figure 6 displays the effect of insulin intake
dictated by each policy, and food intake, on the blood glucose level over 10
evaluation episodes. Note that the food intake schedule is fixed across all
agents and episodes, and that food occurrence might look different across
different plots due to episodes terminating, in some cases, before mealtime.
The starting blood glucose level is also the same in all cases, but from then
on, the trajectories are generated based on each policy. The fixed policy
(Figure 6(a)) administers the same dose given the same blood glucose level,
and hence it is not surprising to see all the episodes terminate in
hypoglycemia due to the prolonged effect of the doses injected at higher blood
glucose levels. The DQN agent (Figure 6(b)) learns to administer a constant
low dose of insulin due to the various penalties in the reward function. A
probable interpretation is that: (1) The agent learns that not doing anything
– i.e.  setting the dose to $0$ – makes the agent spend more time in
hyperglycemia, and incurs penalties (2) Due to the inherent prolongedness of
actions, higher doses quickly drive the agent into hypoglycemia, which also
leads to penalties. Since the agent has no prior information about the
prolonged action effects, it may interpret this penalty as indicating that
choosing high doses is harmful. (3) There is also an additional penalty that
comes from overdosing. Hence, the DQN agent chooses a middle ground of
constantly administering the minimum non-zero dose available (dose = 1) to
maximize its duration into the target zone. The ADRQN agent (Figure 6(c)) is
able to leverage insulin history through its state description and hence
mostly administers new doses sparingly. As shown in the figure, it chooses all
possible doses between 0 and 5. We observe frequent doses when blood glucose
spikes around food intake, but in an attempt to quickly regulate the spike in
blood glucose, the agent often ends up over-dosing. This leads to quick
episode termination due to hypoglycemia right after food consumption. This is
an undesirable behavior, evoking safety concerns, as in the clinical context,
hypoglycemia is more fatal than hyperglycemia (McCrimmon and Sherwin 2010),
and requires immediate medical attention. Finally, the Effective-DQN agent
(Figure 6(d)) appears to be conservative; like the DQN agent, it only injects
the smallest possible dose when it decides to administer insulin. However,
unlike the DQN agent, the Effective-DQN agent has information about the
residual insulin in the blood from previous doses and hence learns when not to
administer further doses. This is an important improvement over DQN, which
closes the performance gap with ADRQN. It is interesting to note that the
Effective-DQN agent also chooses to administer insulin around food intake when
blood glucose spikes, but the dose magnitude remains small. A noteworthy
behavior is that the Effective-DQN agent never ended an episode due to
hypoglycemia. A potential explanation is that either administering a slightly
higher does of 2 units or more frequently administering 1-unit doses might
push the agent into hypoglycemia and eventually death. This hypothesis could
be verified in the future by training an agent with access to smaller insulin
doses. Note however that this behavior is an observation, not a constraint,
and hence not a guarantee on the policy. To have such guarantees one must add
desired conservative properties as constraints, which we leave for future
work.

##  6 Discussion

##### Conclusion.

In this paper, we identified delay and prolongedness as a common roadblock for
using RL in drug dosing. We defined PAE-POMDP, a sub-class of POMDPs in which
action effects persist for multiple time steps, and introduced a simple and
effective framework to convert PAE-POMDPs into MDPs which can be subsequently
solved with traditional RL algorithms. We evaluated our proposed approach on a
toy task and on a glucose control task, for which we proposed a clinically-
inspired reward function for the glucose control simulator, which we hope will
facilitate further research. Our quantitative results have shown that our
proposed approach to convert PAE-POMDPs into MDPs is competitive and offers
important advantages over the baselines. In particular, we have shown that by
introducing domain knowledge on the prolongedness of action effects, we could
build a solution capable of matching state-of-the-art recurrent policies such
as ADRQN while being remarkably more time and memory efficient. Our
qualitative results have further emphasized the benefits of the introduced
Effective-DQN policy, which appeared to be more conservative than the
recurrent policy, despite using the same reward function.

##### Limitations.

Our glucose control experiments are limited to a single virtual patient
environment. Since humans exhibit similar profiles for prolonged action
effects, similar performance on patients from other demographics is expected
but should be evaluated. Although based on pharmacological literature we
acknowledge that delay and prolongedness are present in any drug dosing
scenario, we did not find other simulated environments to test our hypothesis.
Hence, our experiments at the moment are limited to a single drug effect
simulation. Although this work is only in the context of drug dosing,
prolonged action effects seem to be present in other application areas such as
robotics, where the decay assumption may not be justified. The recurrent
solution is still the most general solution in such cases. This work does not
explore those applications.

##### Future work.

In this work, we have only considered discrete action spaces. Exploring
continuous actions might lead to more interesting insights in the future.
Delay and prolongedness appear across all drug dosing scenarios, so we hope
future RL applications of autonomous drug dosing will consciously account for
them. More broadly, future work includes coming up with more generalized
solutions for action superposition that accommodate qualitative actions and
are suitable for application areas beyond drug dosing.

## Acknowledgements

The first author is supported by IVADO PhD fellowship. We are grateful to Mila
IDT (specially Olexa Bilaniuk) for their extraordinary support, and Emmanuel
Bengio for the productive discussions.

## References

  * Annamalai (2010) Annamalai, C. 2010.  Applications of exponential decay and geometric series in effective medicine dosage.  _Advances in Bioscience and Biotechnology_ , 1: 51–54.
  * Benet and Zia-Amirhosseini (1995) Benet, L. Z.; and Zia-Amirhosseini, P. 1995.  Basic principles of pharmacokinetics.  _Toxicologic pathology_ , 23(2): 115–23.
  * Brockman et al. (2016) Brockman, G.; Cheung, V.; Pettersson, L.; Schneider, J.; Schulman, J.; Tang, J.; and Zaremba, W. 2016.  Openai gym.  _arXiv preprint arXiv:1606.01540_.
  * Dalla et al. (2014) Dalla, M. C.; F, M.; D, L.; M, B.; B, K.; and C., C. 2014.  The UVA/PADOVA Type 1 Diabetes Simulator: New Features.  _Journal of Diabetes Science Technology_ , 8(1): 26–34.
  * Dalla, MD, and C. (2009) Dalla, M. C.; MD, B.; and C., C. 2009.  Physical activity into the meal glucose-insulin model of type 1 diabetes: in silico studies.  _Journal of Diabetes Science Technology_ , 3(1): 56–67.
  * Dasgupta and Krasowski (2020) Dasgupta, A.; and Krasowski, M. D. 2020.  Pharmacokinetics and therapeutic drug monitoring.  _Elsevier_ , 1–17.
  * FDA (2015) FDA. 2015.  FY2015 Regulatory Science Research Report:Narrow Therapeutic Index Drugs.  https://www.fda.gov/industry/generic-drug-user-fee-amendments/fy2015-regulatory-science-research-report-narrow-therapeutic-index-drugs.  Accessed: 2022-08-01.
  * Fox et al. (2020) Fox, I.; Lee, J.; Pop-Busui, R.; and Wiens, J. 2020.  Deep Reinforcement Learning for Closed-Loop Blood Glucose Control.  In _Machine Learning for Healthcare_.
  * Guengerich (2011) Guengerich, F. P. 2011.  Mechanisms of drug toxicity and relevance to pharmaceutical development.  _Drug Metabolism and Pharmacokinetics_ , 26(1): 3–14.
  * Hobbie and Roth (2007) Hobbie, R. K.; and Roth, B. J. 2007.  Exponential Growth and Decay.  In _Intermediate Physics for Medicine and Biology_ , 31–47. New York, NY: Springer.
  * Hochreiter and Schmidhuber (1997) Hochreiter, S.; and Schmidhuber, J. 1997.  Long Short-Term Memory.  _Neural Computation_ , 9(8): 1735–1780.
  * Holford (1984) Holford, N. 1984.  Drug concentration, binding, and effect in vivo.  _Pharmaceutical Research_ , 1(3): :102–5.
  * Holford (2018) Holford, N. 2018.  Pharmacodynamic principles and the time course of delayed and cumulative drug effects.  _Translational and Clinical Pharmacology_ , 26(2): 56–59.
  * Kovatchev et al. (2009) Kovatchev, B. P.; Breton, M.; Man, C. D.; and Cobelli, C. 2009.  In silico preclinical trials: A proof of concept in closed-loop control of type 1 diabetes.  _Journal of Diabetes Science and Technology_ , 3(1): 44–55.
  * Lillicrap et al. (2016) Lillicrap, T. P.; Hunt, J. J.; Pritzel, A.; Heess, N.; Erez, T.; Tassa, Y.; Silver, D.; and Wierstra, D. 2016.  Continuous Control with Deep Reinforcement Learning.  In _International Conference on Learning Representation 2016_.
  * Lin et al. (2018) Lin, R.; Stanley, M. D.; Ghassemi, M. M.; and Nemati, S. 2018.  A Deep Deterministic Policy Gradient Approach to Medication Dosing and Surveillance in the ICU.  In _40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)_ , 4927–4931.
  * Lopez-Martinez et al. (2019) Lopez-Martinez, D.; Eschenfeldt, P.; Ostvar, S.; Ingram, M.; Hur, C.; and Picard, R. 2019.  Deep Reinforcement Learning for Optimal Critical Care Pain Management with Morphine using Dueling Double-Deep Q Networks.  _Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)_ , 3960–3963.
  * Maxfield and Zineh (2021) Maxfield, K.; and Zineh, I. 2021.  Precision Dosing: A Clinical and Public Health Imperative.  _JAMA_ , 325(15): 1505–1506.
  * McCrimmon and Sherwin (2010) McCrimmon, R. J.; and Sherwin, R. S. 2010.  Hypoglycemia in type 1 diabetes.  _Diabetes_ , 59(10):2333-2339.
  * Mnih, Kavukcuoglu, and et al. (2015) Mnih, V.; Kavukcuoglu, K.; and et al., D. S. 2015.  Human-level control through deep reinforcement learning.  _Nature_ , 518(7540): 529–533.
  * Nemati, Ghassemi, and Clifford (2016) Nemati, S.; Ghassemi, M. M.; and Clifford, G. D. 2016.  Optimal medication dosing from suboptimal clinical examples: A deep reinforcement learning approach.  In _38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)_ , 2978–2981.
  * Paszke (2021) Paszke, A. 2021.  REINFORCEMENT LEARNING (DQN) TUTORIAL.  https://pytorch.org/tutorials/intermediate/reinforcement˙q˙learning.html.
  * Peschl (2020) Peschl, M. 2020.  ADRQN-PyTorch: A Torch implementation of the action-specific deep recurrent Q network.  https://mlpeschl.com/post/tiny˙adrqn/.
  * Peter and de Paula Julio (2006) Peter, A.; and de Paula Julio. 2006.  The rates of chemical reactions.  In _Atkins’ Physical chemistry_ , chapter 22, 791–823. New York: W.H. Freeman, 8 edition.
  * Ribba et al. (2020) Ribba, B.; Dudal, S.; Lavé, T.; and Peck, R. W. 2020.  Model-Informed Artificial Intelligence: Reinforcement Learning for Precision Dosing.  _Clinical pharmacology and therapeutics_ , 107(4): 853–857.
  * Sharabi et al. (2015) Sharabi, K.; Tavares, C. D. J.; Rines, A. K.; and Puigserver, P. 2015.  Molecular pathophysiology of hepatic glucose production.  _Molecular aspects of medicine_ , 46: 21–33.
  * Tejedor, Woldaregay, and Godtliebsen (2020) Tejedor, M.; Woldaregay, A. Z.; and Godtliebsen, F. 2020.  Reinforcement learning application in diabetes blood glucose control: A systematic review.  _Artificial Intelligence in Medicine_ , 104: 101836.
  * Vogenberg, Barash, and Pursel (2010) Vogenberg, F. R.; Barash, C. I.; and Pursel, M. 2010.  Personalized medicine: part 1: evolution and development into theranostics.  _Pharmacy and Therapeutics_ , 35(10): 560–76.
  * Watkins and Dayan (1989) Watkins, C.; and Dayan, P. 1989.  Q-Learning.  _Machine Learning_.
  * Weng et al. (2017) Weng, W.-H.; Gao, M.; He, Z.; Yan, S.; and Szolovits, P. 2017.  Representation and Reinforcement Learning for Personalized Glycemic Control in Septic Patients.  In _Neural Information Processing Systems_.
  * Xie (2018) Xie, J. 2018.  Simglucose v0.2.1 (2018) [Online].  https://github.com/jxx123/simglucose.
  * Zadeh, Street, and Thomas (2022) Zadeh, S. A.; Street, W. N.; and Thomas, B. W. 2022.  Optimizing Warfarin Dosing using Deep Reinforcement Learning.  _CoRR_ , abs/2202.03486.
  * Zhu, Li, and Poupart (2017) Zhu, P.; Li, X.; and Poupart, P. 2017.  On Improving Deep Reinforcement Learning for POMDPs.  _CoRR_ , abs/1704.07978.

##  Appendix A Sources of Partial Observability in the Glucose control task

Figure 7: Explaining glucose control environment. Insulin has a delayed and
prolonged effect on the blood glucose levels, which can also increase due to
Endogenous Glucose Production (EGP) under fasting. x-axis is time steps,
y-axis is respective units. In the Blood Glucose subplot, the color green
represents the target zone and orange the hyperglycemic zone. In the Blood
Glucose Gradient subplot, the color olive-green represents a reduction and
salmon-pink an increase in the blood glucose level.

As shown in figure 7, there are multiple sources of partial observability in
this environment:
1) Prolonged Effect: The effect of insulin on blood glucose is prolonged over
multiple time steps. A high dose of insulin can quickly terminate the episode
due to hypoglycemia, thereby reducing the episode length. The duration of the
effect depends on an individual’s insulin sensitivity and is proportional to
the magnitude of the dose for lower doses of insulin.
2) Constant Delay: The insulin administered at the sub-cutaneous level takes
some time to reach the bloodstream and be effective. This delay is constant
for an individual.
3) Food: An individual can consume food at any point and hence the blood
glucose can increase unexpectedly. The information of whether a meal has been
consumed is not part of the state space and is an element of surprise, thus
another source of partial observability.
4) Endogenous Glucose Production: When fasting, the human body produces
glucose by breaking down glycogen stored in the liver to meet the brain’s
energy demand (Sharabi et al. 2015).
5) Sensor Noise: The blood glucose measurements by the sensors of Continuous
Glucose Monitor (CGM) can have a little noise, making the CGM readings
slightly different from the actual blood glucose levels. However, to keep the
experiments simple, for the experiments in this paper we made noise as an
identity function, making the CGM readings same as the blood glucose.

##  Appendix B Simglucose simulator changes

We made a few changes to the simglucose simulator:

  1. 1.

Starting state: The original simulator starts all the episodes with the same
blood glucose value of 150, which is not favorable for exploration. Therefore,
we randomize the starting blood glucose in the range of 150 to 195.
Specifically, we modified the reset() method of T1DPatient() class in
simglucose/patient/t1dpatient.py.

  2. 2.

Episode termination: The original simulator terminates the episodes if the
blood glucose is below 70 or higher than 350. We redefine the termination
criteria to 70 and 200. We made this change to keep the search space small for
the experiments presented in this paper.

  3. 3.

Removing Noise: We removed noise by changing CGM = BG in the measure() method
of the CGMSensor() class found in ‘simglucose/sensor/cgm.py’ file.

  4. 4.

Reward function: We implement a reward function as described in the main
paper. The reward function is shared with the code.

  5. 5.

Action space: The original simulator allows insulin dose of upto 30 units.
Therefore, the action space is continuous between 0 and 30. We reduce the
action space to discrete values of 0 to 5 to simplify the task for the paper.

  6. 6.

Environment and Episode: we randomly select the patient adult009 as the
environment to train the RL agents presented in this paper and assume that we
do not have access to the transition dynamics. We slightly modify the
environment to redefine episodes. Our episodes begin uniformly in the blood
glucose range of (150, 195) and terminate if the blood glucose is greater than
200 (hyperglycemic death) or less than 70 (hypoglycemic death). We define
hypoglycemia as the condition when blood glucose ranges from 70 to 100, target
zone as 100 to 150 and hyperglycemia as 150 to 200. We also reduce the action
space to discrete doses of insulin of up to 5 units. The objective is to
effectively administer insulin to maximize the time spent in the target blood
glucose zone.

For more details or the exact implementation, please refer to the copy of the
updated simulator code added to the supplementary package. We will make them
public, on acceptance of the paper.

##  Appendix C Experimental Details

###  C.1 Architecture

Both move block and glucose control task use the same architecture, but only
differ in the number of neurons.

  * •

State embedding: For both the simulator (simglucose and MoveBlock), state is a
1 dimensional continuous value. We use a fully connected layer of 16 neurons
with linear activation to learn a state-embedding. We train this state-
embedding together with the rest of the Q-network. Although weights are not
shared, we use the same state-embedding architecture for all three of our
models.

  * •

Action embedding: Action is a 1 dimensional discrete value. We use a fully
connected layer with linear activation to learn an action-embedding. Effective
DQN and ADRQN use the same architecture for the action embedding. They don’t
share weights and we train them together with the rest of the Q-network and
state-embedding if applicable. MoveBlock uses 32 neurons and the glucose
control task uses 16 neurons to encode the actions.

  * •

DQN: DQN passes the state-embedding through a fully connected layer and ReLU
activation before passing it through a fully connected output layer.

  * •

Effective DQN: Effective DQN concatenates the state and action embeddings
before passing them through a fully connected layer with Relu activation
followed by a fully connected output layer.

  * •

ADRQN: The original ADRQN architecture uses convolutional layers, but since
our application has low dimensional input, we replace the convolutional layers
with linear layers similar to (Peschl 2020). It concatenates the state and
action embeddings before passing them through a LSTM unit followed by a fully
connected output layer.

MoveBlock uses 512 neurons and glucose control uses 256 neurons for the hidden
layer.

###  C.2 Hyper-parameters

We selected the hyper parameters by grid-search and present the final set of
selected hyper-parameters for the function approximation experiments in Table
2.

We use a learning rate of 0.001 for tabular q-learning and 0.005 for effective
q-learning on the tabular toy task. We choose a lambda value of 0.99 for the
effective q-learning. For both the algorithms, we use an epsilon greedy policy
that decays exponentially from 0.9 to 0.05 over a span of 10K episodes.

Table 2: Hyperparameters used to train agents on the glucose control task

Hyper-parameter | DQN | ADRQN | Eff-DQN
---|---|---|---
Glucose Control Task |  |  |
Learning Rate | 1e-6 | 0.001 | 0.001
State embedding dimension | 16 | 16 | 16
Action embedding dimension | N/A | 16 | 16
history length | N/A | N/A | 15
lambda | N/A | N/A | 0.99
hidden layer dimension | 256 | 256 | 256
seed | 1-5 | 1-5 | 1-5
replay buffer size | 100K | 100K | 100K
batch size | 512 | 512 | 512
Explore | 1000 | 1000 | 1000
Target network update rate | 1 | 1 | 1

###  C.3 Implementation details

We used the DQN implementation from the PyTorch DQN tutorial (Paszke 2021) of
PyTorch forum. We implemented effective DQN on top of this. The ADRQN
implementation published with the main paper was in C++ and Caffe. Therefore,
we have used an open-source PyTorch implementation (Peschl 2020). We self
implemented the driver and helper codes and shared them with this
supplementary package.

###  C.4 Qualitative Analysis Data generation

Table 3: Data generation for the Qualitative Analysis

seed | Starting Blood Glucose | Ending Blood Glucose
---|---|---
|  | DQN | ADRQN | Effective-DQN
100 | 191.57 | 65.60 | 65.01 | 200.26
110 | 184.02 | 63.19 | 200.66 | 200.49
120 | 162.14 | 69.08 | 208.86 | 209.82
130 | 180.80 | 62.17 | 224.56 | 224.90
140 | 160.92 | 68.61 | 214.56 | 230.39
150 | 154.27 | 66.07 | 63.10 | 243.63
160 | 155.43 | 67.66 | 66.16 | 201.68
170 | 174.21 | 60.08 | 205.03 | 200.42
180 | 193.46 | 66.84 | 68.64 | 201.52
190 | 177.90 | 61.25 | 200.37 | 202.49

Figure 8: Prolonged and Additive Action Effect in MoveBlock Environment.
_Prolonged Effect:_ Position (observed variable) does not change initially in
the absence of force (action). Unit force is exerted at t=3 followed by no
force till t=40. But since velocity (unobserved variable) keeps on changing
under the previously exerted force until friction brings it to zero, position
keeps on changing too in the absence of any further force. _Additive Effect:_
Forces exerted at t=41 and t=57 exhibit additive nature.

For the quantitative assesment, we trained 5 agents corresponding to seeds 1
to 5 for DQN, ADRQN and Effective-DQN. Out of the trained 5 agents we randomly
chose the agents with seed 3 for all the policies. For the qualitative
assessment, we fix 10 starting blood glucose levels and use them to generate
trajectories for each policy under evaluation. In table 3, we present the
seeds used to generate the episodes. For sanity check, we also present the
starting blood glucose levels corresponding to these 10 seeds. The
trajectories depend on the policy in use. We report the end of trajectory
blood glucose levels at the end of each episode.

##  Appendix D Prolonged action effect on MoveBlock toy environment

Figure 8 demonstrates the effect of prolonged actions on the states on the toy
environment. We observe that actions have a prolonged effect on the latent
variable velocity, which in turn affects the observed variable position. In
the toy environment we could control the prolonged effect by changing the co-
efficient of kinetic friction. We used a coefficient of 0.003 for the
environment used in the experiments. The reward function is same as the the
Mountain Car environment:

$r(s,a)=\begin{cases}100&\text{if $s\in S_{goal}$}\\\
-a^{2}*0.1&\text{otherwise}\\\ \end{cases}$ (5)


# Simulation of the Einstein–de Haas effect combining molecular and spin
dynamics

W. Dednam<EMAIL_ADDRESS>C. Sabater  A. E. Botha  E. B. Lombardi  J.
Fernández-Rossier  M. J. Caturla  Department of Physics, Unisa Science Campus,
University of South Africa, Johannesburg 1710, South Africa  Departamento de
Física Aplicada and Unidad asociada CSIC, Universidad de Alicante, Campus de
San Vicente del Raspeig, E-03690 Alicante, Spain.  International Iberian
Nanotechnology Laboratory, 4715-310 Braga, Portugal

###### Abstract

The spin and lattice dynamics of a ferromagnetic nanoparticle are studied via
molecular dynamics and with semi-classical spin dynamics simulations where
spin and lattice degrees of freedom are coupled via a dynamic uniaxial
anisotropy term. We show that this model conserves total angular momentum,
whereas spin and lattice angular momentum are not conserved. We carry out
simulatios of the the Einstein–de Haas effect for a Fe nanocluster with more
than 500 atoms that is free to rotate, using a modified version of the open-
source spin-lattice dynamics code (SPILADY). We show that the rate of angular
momentum transfer between spin and lattice is proportional to the strength of
the magnetic anisotropy interaction. The addition of the anisotropy allows
full spin-lattice relaxation to be achieved on previously reported timescales
of $\sim 100$ ps and for tight-binding magnetic anisotropy energies comparable
to those of small Fe nanoclusters.

###### keywords:

Spin-lattice dynamics , Molecular dynamics simulations , Magnetic anisotropy ,
Conservation of total angular momentum , Angular momentum exchange ,
Einstein–de Haas effect

††journal: Computational Materials Science

##  1 Introduction

First observed by Einstein and De Haas [1] in 1915, the so-called Einstein–de
Haas (EdH) effect involves the rotation of a magnetic object through change in
its magnetization. The converse effect, i.e. magnetization through a
rotational change, was discovered contemporaneously by Barnett [2]. Both
effects rely on the conservation of the total angular momentum and are
governed by the same gyromagnetic tensor components, satisfying the Onsager
reciprocal relations [3, 4, 5]. Thus, for any particular physical system,
experimental confirmation of either effect is thought to be sufficient to
establish conservation of total angular momentum [6]. Arguably, the EdH effect
is more widely known because of its routine use in measuring the $g$ factors
of various ferromagnetic materials: Fe, Co, Ni, etc. [7]. Recently, the EdH
effect has also been observed in a ferromagnetic thin film, deposited on a
microcantilever [8].

Despite several theoretical studies [10, 11, 12, 13], a microscopic
description of the two effects, either in terms of computer simulations or
theoretically, remains a challenging and open problem. On the theoretical
front the most recent work by Mentink et al. [13] shows that, to accurately
determine the timescale of angular momentum exchange between electronic and
lattice degrees of freedom, and to establish whether the exchange is local or
non-local, it is necessary to cast both quantum subsystems in the angular
momentum representation. In practice, this involves “dressing” every atom with
phonons that exchange angular momentum locally with the multi-electron atom’s
spins and electronic orbital angular momenta. The quantum many-body
Hamiltonian obtained in this way is fully rotationally invariant and a static
non-perturbative variational calculation shows that phonons quench the
electronic orbital angular momentum of the atoms as the coupling strength
between the lattice and electronic degrees of freedom increases. A dynamic
non-perturbative variational calculation by the same authors [13], on the
other hand, reveals the emergence of high frequency peaks in the magnetic
susceptibility as the electron-phonon coupling strength increases. These peaks
are consistent with the timescale of the ultra-fast EdH effect that was
observed recently in experiments of laser-induced demagnetization in an
isolated ferromagnetic sample [38].

Transfer of angular momentum between spin and motional degrees of freedom
ultimately occurs due to spin-orbit interactions that couple these degrees of
freedom. In the spin-lattice dynamics (SLD) model, spin-orbit interactions
lead to magnetic anisotropy terms that depend explicitly on the atomic
positions, thereby coupling spin and lattice degrees of freedom [15, 16, 17,
18, 9, 19, 20, 21, 22]. Importantly, total angular momentum must be conserved.

In this work we implement the SLD model of Ref. [17]. The model is based on a
truncated Taylor expansion of the magnetic anisotropy energy [23]. We
demonstrate that it is able to reproduce spin-lattice relaxation times $\sim
100$ ps, consistent with previously reported values [15, 17, 20]. Our
implementation of the SLD model is applicable to combined molecular and spin
dynamics of systems of arbitrary symmetry (including bulk, clusters, defect
systems, etc.) It allows the modelling of a wide range of potential new
systems at finite temperatures, for example: non-collinear magnetic domain
walls in ferromagnetic atomic-sized contacts or wires [24, 25, 26], chiral
spin textures in magnetic nanoparticles [27, 28], skyrmions [29, 30], and an
atomistic description of the Barnett [2] and Barkhausen [31] effects, amongst
others.

The remaining material is organized as follows. In Section 2, we give a brief
outline of the SLD model of ferromagnetic exchange previously implemented in
SPILADY [32] and demonstrate analytically that it conserves angular momentum.
However, we show that angular momentum exchange is not possible when only the
isotropic ferromagnetic exchange term is considered in the model. We therefore
show analytically in Section 2 that, by adding the second-order anisotropy
correction of Ref. [17] to the SLD model, not only is total angular momentum
conserved, but also exchange of angular momenta between spin and lattice is
facilitated. In Section 3, we discuss our choice of semi-classical potentials
and parameters to perform SLD simulations of ferromagnetic iron nanoclusters
and then demonstrate numerically that our implementation replicates the
Einstein–de Haas effect in a small iron nanocluster. Finally, we conclude the
paper by summarizing our findings in Section 4.

##  2 Models and analytical considerations

###  2.1 General description of the SLD model

Our implementation builds on the well-established technique known as spin-
lattice dynamics (SLD) [33, 15, 32, 20, 9]. This is a semi-classical molecular
dynamics approach which allows the simultaneous evolving of atomic spin and
position degrees of freedom. The atom-atom interactions in this model are
described by the widely-used embedded atom method (EAM) potential [34]. The
EAM potential has been implemented in SPILADY [33, 32] along with a
generalized Heisenberg model of ferromagnetic exchange:

$H=\sum_{i}^{N}\frac{\overrightarrow{p_{i}}^{2}}{2m}+V+H_{\text{s}pin}\,$ (1)

In Eq. (1), $\overrightarrow{{p}_{i}}$ is the momentum operator for the $i$th
atom and $V$ is the interatomic EAM potential energy function that models the
basic interactions between the atoms, i.e. excluding the spin interactions.
Expressed in terms of the spin angular momentum vectors of particles $i$ and
$j$, the spin Hamiltonian $H_{\rm spin}$ is a sum of the Zeeman term ($H_{\rm
Zeeman}$), the generalized Heisenberg exchange term ($H_{\rm exch}$) and the
anisotropy term ($H_{\rm anis}$):

$H_{\rm spin}=H_{\rm Zeeman}+H_{\rm exch}+H_{\rm anis}$ (2)

with

$H_{\rm
Zeeman}=-\sum_{i}\vec{M}_{i}\cdot\vec{B}_{\text{e}xt}=\sum_{i}g\mu_{B}\vec{S}_{i}\cdot\vec{B}_{\text{e}xt},$
(3)

and

$H_{\rm
exch}=-\frac{1}{2}\sum_{i}J_{ij}\left(r_{ij}\right)\left[{\vec{S}}_{i}\cdot{\vec{S}}_{j}-\left|{\vec{S}}_{i}\right|\left|{\vec{S}}_{j}\right|\right].$
(4)

The anisotropy term, $H_{\rm anis}$, is discussed in the next section. The
Zeeman term in Eq. (3) is due to an external magnetic field $\vec{B}_{\rm
ext}$, where $\vec{M}_{i}=-g\mu_{B}\vec{S}_{i}$ is the magnetic dipole moment,
$g=2.002319$ is the electron’s gyromagnetic factor and $\mu_{\rm B}$ is the
Bohr magneton. The exchange term, Eq. (4), is isotropic in form, with
$r_{ij}=|\vec{r}_{i}-\vec{r}_{j}|$ being the magnitude of the displacement
$\vec{r}_{ij}=\vec{r}_{i}-\vec{r}_{j}$ between the $i$th and $j$th atoms.
${\vec{S}}_{i}$ is the (dimensionless) atomic spin or spin angular momentum
vector at site $i$ and the quantities $J_{ij}\left(r_{ij}\right)$ are obtained
from a shifted fitting function of the form [33, 17]

$J_{ij}\left(r_{ij}\right)=J_{0}\left(1-r_{ij}/r_{c}\right)^{3}\Theta\left(r_{c}-r_{ij}\right),$
(5)

in which $J_{0}$ is the only fitting parameter,
$\Theta\left(r_{c}-r_{ij}\right)$ is the Heavyside step function, and $r_{c}$
is the cut-off radius.

Next, we discuss the conservation of total angular momentum in the Heisenberg
model of ferromagnetic exchange implemented in SPILADY.

###  2.2 Conservation of total angular momentum in SLD simulations

We briefly revisit the textbook derivation of angular momentum conservation in
a system of classical particles. We define the total angular momentum of a
system of classical particles:

$\vec{L}=\sum_{i}\vec{L}_{i}$ (6)

where $\vec{L}_{i}=m_{i}\vec{r}_{i}\times\vec{v}_{i}$ is the angular momentum
of an individual particle. The rate of change of the total angular momentum is
governed by the following equation:

$\frac{d\vec{L}}{dt}=\sum_{i}\frac{d\vec{L}_{i}}{dt}=\sum_{i}\vec{r}_{i}\times\vec{F}_{i}$
(7)

We can break down the force acting on particle $i$ into two parts, the forces
coming from external sources outside of the system and the internal forces,
$\vec{F}_{i}=\left(\vec{F}_{i,\rm{ext}}+\vec{F}_{i,\rm{int}}\right)$, where
the internal forces are given by

$\vec{F}_{i,\rm{int}}=\sum_{j\neq i}\vec{f}_{ji}$ (8)

and where $\vec{f}_{ji}$ is the force that particle $j$ exerts on $i$. By
Newton’s third law, $\vec{f}_{ji}=-\vec{f}_{ij}$. Therefore, for any pair of
particles $i$ and $j$, the contributions to the rate of change of the angular
momentum deriving from their mutual interaction only, are:

$\left.\frac{d\vec{L}_{i}}{dt}\right|_{{\rm
int},j}+\left.\frac{d\vec{L}_{j}}{dt}\right|_{{\rm
int},i}=\vec{r}_{i}\times\vec{f}_{ji}+\vec{r}_{j}\times\vec{f}_{ij}=\left(\vec{r}_{i}-\vec{r}_{j}\right)\times\vec{f}_{ji}$
(9)

Now, if the internal forces $\vec{f}_{ij}$ are central, then they do not
contribute to the rate of change of angular momentum of the system, and we are
left with only the external forces. Moreover, a force
$\vec{f}_{ij}=-\vec{\nabla}U(|\vec{r}_{i}-\vec{r}_{j}|)$ that can be obtained
as a gradient of a potential $U$, such as the EAM potential [34] used in this
work, is automatically central and, by virtue of being expressed as a gradient
of a scalar function, conservative.

This shows that, for an isolated system of classical particles, in the absence
of external forces, and with internal forces that are central, the rotational
angular momentum is conserved.

Equation (3) shows that the atomic spins in our implementation are
dimensionless. Consequently, in order to calculate the dimensionless classical
angular momentum, we define $\vec{\lambda}=\frac{\vec{L}}{\hbar}$ with
${\hbar}=6.58211899\times 10^{-16}$ eV.s in the following expressions
involving the classical angular momentum. So, substituting the Heisenberg
exchange term, Eq. (4), into Eq. (9) yields for the rate of change of the
lattice angular momentum:

$\begin{gathered}\left.\frac{d\vec{\lambda}_{i}}{dt}\right|_{{\rm
exch},j}+\left.\frac{d\vec{\lambda}_{j}}{dt}\right|_{{\rm
exch},i}=\frac{1}{\hbar}\left[\vec{r}_{i}\times\vec{f}_{ji}+\vec{r}_{j}\times\vec{f}_{ij}\right]=\frac{1}{\hbar}\left[\vec{r}_{i}\times-\vec{\nabla_{i}}H_{{\rm
exch,}j}+\vec{r}_{j}\times-\vec{\nabla_{j}}H_{{\rm exch,}i}\right]\\\
=-\vec{r}_{i}\times\frac{1}{\hbar
r_{\mathit{ij}}}\frac{dJ_{ij}(r_{\mathit{ij}})}{dr_{\mathit{ij}}}\left[{\vec{S}}_{i}\cdot{\vec{S}}_{j}-\left|{\vec{S}}_{i}\right|\left|{\vec{S}}_{j}\right|\right]\vec{r}_{ij}+\vec{r}_{j}\times\frac{1}{\hbar
r_{\mathit{ji}}}\frac{dJ_{ji}(r_{\mathit{ji}})}{dr_{\mathit{ji}}}\left[{\vec{S}}_{j}\cdot{\vec{S}}_{i}-\left|{\vec{S}}_{j}\right|\left|{\vec{S}}_{i}\right|\right]\vec{r}_{ij}\\\
=(-\vec{r}_{i}+\vec{r}_{j})\times\frac{1}{\hbar
r_{\mathit{ij}}}\frac{dJ_{ij}(r_{\mathit{ij}})}{dr_{\mathit{ij}}}\left[{\vec{S}}_{i}\cdot{\vec{S}}_{j}-\left|{\vec{S}}_{i}\right|\left|{\vec{S}}_{j}\right|\right]\vec{r}_{ij}\\\
=\frac{1}{\hbar
r_{\mathit{ij}}}\frac{dJ_{ij}(r_{\mathit{ij}})}{dr_{\mathit{ij}}}\left[{\vec{S}}_{i}\cdot{\vec{S}}_{j}-\left|{\vec{S}}_{i}\right|\left|{\vec{S}}_{j}\right|\right]\vec{r}_{ij}\times\vec{r}_{ij}=0\end{gathered}$
(10)

We therefore see that the classical or rotational angular momentum is
conserved in an isotropic model of ferromagnetic exchange.

With a view to establishing conservation of spin angular momentum in SLD
models, we need to derive an expression analogous to Eq. (9) involving the
rates of change of the angular momentum of spins $i$ and $j$. To that end, we
start with the equation of motion for the classical spin momentum of particle
$i$:

$\frac{d\vec{S}_{i}}{dt}=-\vec{S}_{i}\times\frac{1}{\hbar}\nabla_{\vec{S_{i}}}H_{\rm
spin}$ (11)

where the energy of the spin system, $H_{\rm spin}$, is in general a function
of $\vec{S}_{i}$, $\vec{S}_{j}$ and of the coordinates of the atoms and their
velocities.

We first compute the effective field associated with the first two terms of
Eq. (2), classically. For the first term, $H_{\rm Zeeman}$, we obtain:

$\nabla_{\vec{S_{i}}}H_{\rm Zeeman}=g\mu_{B}\vec{B}_{i}$ (12)

This term produces a rate of change of the spin momentum when $\vec{S}_{i}$ is
not parallel to $\vec{B}_{i}$, and, consequently, we do not consider total
angular momentum conservation in the presence of an external magnetic field
any further.

The second term, $H_{\rm exch}$, produces central forces when we compute the
associated gradients. Since individual spin magnitudes are conserved in our
model, we have for the spin torques:

$\nabla_{\vec{S_{i}}}H_{\rm exch}=-\frac{1}{2}\sum_{j}J(r_{ij})\vec{S}_{j}$
(13)

By analogy with Eq. (9), we find that for any pair of particles $i$ and $j$
the contribution to the rate of change of the spin angular momentum coming
from this term vanishes:

$\begin{gathered}\left.\frac{d\vec{S}_{i}}{dt}\right|_{{\rm
exch,}j}+\left.\frac{d\vec{S}_{j}}{dt}\right|_{{\rm exch,}i}=\\\
\frac{J(r_{ij})}{\hbar}\left(\vec{S}_{i}\times\vec{S}_{j}+\vec{S}_{j}\times\vec{S}_{i}\right)=0\end{gathered}$
(14)

In other words, isotropic interactions of the scalar form
$\vec{S}_{1}\cdot\vec{S}_{2}$ preserve the spin angular momentum pairwise.

Therefore, the Heisenberg model of ferromagnetic exchange in Eq. (4) not only
conserves total angular momentum, but does so separately within the spin and
lattice sub-reservoirs, and cannot account for angular momentum transfer
between lattice and spin degrees of freedom.

In the next section, we consider the correction $H_{\rm anis}$, which couples
the spin and lattice degrees of freedom, allowing transfer of angular momentum
between spin and lattice, while still conserving total angular momentum.

###  2.3 Addition of magnetic anisotropy

To permit exchange of angular momentum between spins and lattice, we have
implemented in the SPILADY code [32] the model of anisotropy described in Ref.
[17]. We exclude the first-order anisotropy correction of Ref. [17] since it
does not preserve time-reversal symmetry in the absence of an external field.
Therefore we only consider the second-order (uniaxial) anisotropy correction
of Ref. [17] in our implementation:

$H_{\rm
anis}=-C_{2}\sum_{i,\alpha,\beta}S_{\alpha}(i)\Lambda_{\alpha,\beta}(i)S_{\beta}(i),$
(15)

In Eq. (15),
$\Lambda_{\alpha,\beta}(i)=\frac{{\partial^{2}\rho_{i}(r_{ij})}}{{\partial
x_{i\alpha}\partial x_{i\beta}}}$, with $\alpha,\beta=1,2,3$ Cartesian
components, is the Hessian matrix element with respect to lattice site $i$, of
the semi-classical many-body, EAM electron density-like, function:

$\begin{gathered}\rho_{i}\left(r_{ij}\right)=\left\\{\begin{array}[]{ll}\sum_{j\neq
i}\left(1-\frac{r_{ij}}{r_{0}}\right)^{4}e^{\left(1-\frac{r_{ij}}{r_{0}}\right)}&\quad
r_{ij}\leq r_{0}\\\ 0&\quad r_{ij}>r_{0},\end{array}\right.\end{gathered}$
(16)

where $r_{0}$ is the cut-off distance between second and third nearest
neighbours and $C_{2}$ is a proportionality constant which determines the
relative strength of the anisotropy term [17]. 111The interested reader is
referred to Ref. [35] for full details of our implementation in SPILADY of the
algorithms described in Refs. [36, 37, 16, 17].

Thus, the contribution of the anisotropy to the effective local magnetic field
at atom i from neighbouring spins within the cut-off $r_{0}$ becomes:

$\begin{gathered}\nabla_{\vec{S_{i}}}H_{\text{anis}}=-C_{2}\sum_{j{\neq}i}\left[g\left(r_{\mathit{ij}}\right)\vec{S_{i}}+\frac{1}{r_{\mathit{ij}}}\frac{\mathit{dg}\left(r_{\mathit{ij}}\right)}{dr_{\mathit{ij}}}\left(\vec{r}_{\mathit{ij}}\cdot\vec{S_{i}}\right)\vec{r}_{\mathit{ij}}\right]\end{gathered}$
(17)

where
$g(r_{\mathit{ij}})\text{=}\frac{1}{r_{\mathit{ij}}}\frac{d\phi(r_{\mathit{ij}})}{dr_{\mathit{ij}}}$
with $r_{\mathit{ij}}\text{=}|\vec{r_{i}}-\vec{r_{j}}|$ and
$\phi(r_{\mathit{ij}})$ is the expression under the summation in Eq. (16).

Evaluating
$\frac{d\vec{S}_{i}}{dt}=-\vec{S}_{i}\times\frac{1}{\hbar}\nabla_{\vec{S_{i}}}H_{\rm
anis}$ for the spin at site $i$ yields

$\begin{gathered}\frac{d\vec{S}_{i}}{dt}=\frac{C_{2}}{\hbar}\sum_{k{\neq}i}\left[\frac{1}{r_{\mathit{ik}}}\frac{\mathit{dg}\left(r_{\mathit{ik}}\right)}{dr_{\mathit{ik}}}\left(\vec{r}_{\mathit{ik}}\cdot\vec{S_{i}}\right)\vec{S}_{i}\times\vec{r}_{\mathit{ik}}\right]\end{gathered}$
(18)

In the same way, for the spin at site $j$,

$\begin{gathered}\frac{d\vec{S}_{j}}{dt}=\frac{C_{2}}{\hbar}\sum_{k{\neq}j}\left[\frac{1}{r_{\mathit{jk}}}\frac{\mathit{dg}\left(r_{\mathit{jk}}\right)}{dr_{\mathit{jk}}}\left(\vec{r}_{\mathit{jk}}\cdot\vec{S_{j}}\right)\vec{S}_{j}\times\vec{r}_{\mathit{jk}}\right]\end{gathered}$
(19)

Finally, for the $j$th term from $\frac{d\vec{S}_{i}}{dt}$ and the $i$th term
from $\frac{d\vec{S}_{j}}{dt}$, summing gives

$\begin{gathered}\left.\frac{d\vec{S}_{i}}{dt}\right|_{j}+\left.\frac{d\vec{S}_{j}}{dt}\right|_{i}=\frac{C_{2}}{\hbar
r_{\mathit{ij}}}\frac{\mathit{dg}\left(r_{\mathit{ij}}\right)}{dr_{\mathit{ij}}}\left[\left(\vec{r}_{\mathit{ij}}\cdot\vec{S_{i}}\right)\vec{S}_{i}\times\vec{r}_{\mathit{ij}}+\left(\vec{r}_{\mathit{ij}}\cdot\vec{S_{j}}\right)\vec{S}_{j}\times\vec{r}_{\mathit{ij}}\right]\end{gathered}$
(20)

since two instances of $\vec{r}_{\mathit{ji}}$ = -$\vec{r}_{\mathit{ij}}$
occur in the second term so their minus signs cancel when the indices of $i$
and $j$ are swapped.

Next, we consider the lattice torque term, showing that it exactly compensates
the spin torque term of Eq. (20).

The expression for the “anisotropy” force on the spin at site $i$ due to all
neighbours $j$ within the cut-off of Eq. (16) is:

$\begin{gathered}\vec{f}_{{\rm anis},i}=\frac{1}{2}C_{2}\sum_{j\neq
i}\frac{1}{r_{\mathit{ij}}}\left[\frac{\mathit{dg}\left(r_{\mathit{ij}}\right)}{dr_{\mathit{ij}}}\left\\{\vec{S_{i}}{\cdot}\vec{S_{i}}+\vec{S_{j}}{\cdot}\vec{S_{j}}\right\\}\vec{r}_{\mathit{ij}}\right.\\\
-\frac{\mathit{dg}\left(r_{\mathit{ij}}\right)}{dr_{\mathit{ij}}}\left\\{\frac{\left(\vec{r}_{\mathit{ij}}{\cdot}\vec{S_{i}}\right)^{2}}{r_{\mathit{ij}}^{2}}-\frac{\left(\vec{r}_{\mathit{ij}}\cdot\vec{S_{j}}\right)^{2}}{r_{\mathit{ij}}^{2}}\right\\}\vec{r}_{\mathit{ij}}\\\
+2\frac{\mathit{dg}\left(r_{\mathit{ij}}\right)}{dr_{\mathit{ij}}}\left\\{\left(\vec{r}_{\mathit{ij}}\cdot\vec{S_{i}}\right)\vec{S_{i}}+\left(\vec{r}_{\mathit{ij}}\cdot\vec{S_{j}}\right)\vec{S_{j}}\right\\}\\\
\left.+\frac{1}{r_{\mathit{ij}}}\frac{d^{2}g\left(r_{\mathit{ij}}\right)}{d^{2}r_{\mathit{ij}}}\left\\{\left(\vec{r}_{\mathit{ij}}{\cdot}\vec{S_{i}}\right)^{2}+\left(\vec{r}_{\mathit{ij}}{\cdot}\vec{S_{j}}\right)^{2}\right\\}\vec{r}_{\mathit{ij}}\right]\end{gathered}$
(21)

Then, by substituting Eq. (21) into (9), and considering only the pair of
particles $i$ and $j$, we obtain after some algebra,

$\begin{gathered}\left.\frac{d\vec{\lambda}_{i}}{dt}\right|_{j}+\left.\frac{d\vec{\lambda}_{j}}{dt}\right|_{i}=\frac{1}{\hbar}\vec{r}_{ij}\times\vec{f}_{ji}\\\
=\frac{C_{2}}{\hbar
r_{\mathit{ij}}}\frac{\mathit{dg}\left(r_{\mathit{ij}}\right)}{dr_{\mathit{ij}}}\left[\left(\vec{r}_{\mathit{ij}}\cdot\vec{S}_{i}\right)\vec{r}_{ij}\times\vec{S}_{i}+\left(\vec{r}_{\mathit{ij}}\cdot\vec{S}_{j}\right)\vec{r}_{ij}\times\vec{S}_{j}\right],\end{gathered}$
(22)

where we have used the fact that the cross product of a vector with itself is
zero.

Clearly, Eqs. (20) and (22) compensate each other when the cross products are
commuted in Eq. (22) to appear in the same order as in Eq. (20). We have
therefore shown analytically that the model of Ref. [17] permits exchange of
angular momentum, while at the same time guaranteeing conservation of total
angular momentum.

##  3 Numerical simulation results

We validate our implementation of the anisotropy correction by studying the
phenomenon known as the Einstein–de Haas effect [1]. In its simplest version,
the EdH experiment involves using a strong magnetic field to change the
magnetization of a small freely suspended cylinder made of a ferromagnetic
material. The sudden change in magnetization in the cylinder brought on by the
strong magnetic field then manifests as a physical rotation about the axis of
the sample as it evolves towards a new equilibrium magnetization due to the
conservation of the total (spin + rotational) angular momentum. The effect is
well-known, but its microscopic origin has only been recently demonstrated in
ultrafast (femtosecond scale) laser-induced demagnetization experiments [38],
though in that particular case the electronic degrees of freedom are also
involved since the laser primarily excites the electrons [39, 40].

###  3.1 Simulation parameters for Fe

In our numerical simulations of the EdH effect, we use a small prolate
spheroid made of iron. We choose this shape because the anisotropy is uniaxial
and a positive value of $C_{2}$ will see the spins aligning preferentially
along the long axis of the spheroid. It consists of 505 Fe atoms cut out from
a bulk periodic simulation cell of body-centered cubic (BCC) Fe(001). The
particle is perfectly symmetric about its central axis along the $z$
direction, i.e., the $x$ and $y$ components of each of its constituent atoms
sum to $0$ with respect to this axis.

To describe the interatomic interactions, we use the Malerba et al. [41] EAM
potential with a per-atom cohesive energy of $-4.122$ eV/atom. This contrasts
with the $0$ eV cohesive energy of the Harmonic and Morse potentials used in
Refs. [22, 9]. The EAM potential used in this work is more realistic in that
it is able to more accurately reproduce material properties [41], including
defects and the ab-initio energy of a free Fe(001) surface.

Then, for the Heisenberg exchange term in Eq. (4), we employ the same
parameter values as in Table II of Ref. [22]: $J_{0}=0.904$ eV and
$r_{c}=3.75$ Å.

Finally, for the anisotropy term in Eq. (15), we use $C_{2}=0.025$ eV.Å2 and
$r_{0}=3.5$ Å. 222We have discovered that the units of $C_{1}$ and $C_{2}$ are
quoted incorrectly in Ref. [17], since they should be eV.Å and eV.Å2,
respectively. This value of $C_{2}$ has been estimated by determining the
difference in magnetic anisotropy energy (MAE) of the relaxed spheroid, with
spins and atoms frozen, compared to the MAE with spins perpendicular to the
first configuration, resulting in $\Delta\text{MAE}\sim 0.2$ meV/atom, in good
agreement with the value of $\sim 0.2$ meV/atom reported in Fig. 7 of Ref.
[42] for a similar-sized Fe nanocluster using a tight-binding formalism. 333We
note that our implementation does not model the weak static background
anisotropy of 3$d$ transition metals, which, for example, gives rise to their
preferred bulk magnetic easy axes. Rather, it accounts for dynamic anisotropy
due to local symmetry-breaking about a given lattice site that may dominate at
sufficiently high temperatures approaching the Debye temperature of the metal
[23]. Nevertheless, dynamic many-body effects, such as those used in this
model, play an important role [13] in fast demagnetization processes in
ferromagnetic materials that cannot be modelled by static background
anisotropy alone.

###  3.2 SLD simulation of the Einstein–de Haas effect

We simulate the EdH effect by mirroring an experiment where an initially
demagnetized sample (such as via heat shock) is allowed to relax to an
equilibrium magnetized state, while observing the rotational angular momentum
which is induced by spins flipping and precessing as the sample reaches its
equilibrium magnetization.

We therefore use an input configuration for our SLD simulation which consists
of fully randomized spins, as shown in Fig. 1 (a). The sample is then allowed
to relax for a period of $1$ ns, reaching an equilibrium magnetized state
after $\sim 100$ ps, as illustrated in Fig. 1 (b), where we see the spins have
mostly aligned parallel to the easy axis. 444See Supplemental Material at
http://link.xxx.org/supplemental/xxxx/xxxx.xxx for movies of the atoms and
spins of the spheroid in Fig. 1 for $C_{2}=0.025,0.05,0.075,0.1,0.15$ and
$0.2$ eV.Å2 showing how the nanocluster evolves to the final states depicted
in Fig. 1.

Figure 1: (a) Initial ($t=0$) spin configuration of the Fe(001) prolate
spheroid. The atoms (not shown) are placed near the sites they would occupy in
a perfectly ordered BCC lattice of Fe(001) and the spins are initially
randomly oriented. (b) After $t=100$ ps, the spins mostly point upward. The
color legend refers to the magnitude of the (dimensionless) atomic spins on
the positive z-axis, i.e., the direction of saturation magnetization of Fe,
$2.2\mu_{\rm B}$, divided by the electron’s gyromagnetic factor, $g=2.002319$,
and $1.0\mu_{\rm B}$. (c) Conservation of total average per-atom angular
momentum (red) as a sum of contributions from the atomic spins
$\left|\left<\vec{S}\right>\right|$ (green) and lattice
$-\left|\left<\vec{\lambda}\right>\right|$ (purple) for $C_{2}=0.025$ eV.Å2.
For convenience, the negative of $\left|\left<\vec{\lambda}\right>\right|$ has
been plotted so that $\left|\left<\vec{S}\right>\right|$ coincides with the
general orientation of $\left<\vec{S}\right>$ in (b). Hyperbolic tangent fits
to $\left|\left<\vec{S}\right>\right|$ (orange) and
$-\left|\left<\vec{\lambda}\right>\right|$ (yellow) are also shown. (d) The
lattice temperature (purple), initially $0$ K, rises rapidly while the spin
temperature (green), initially $\sim 12000$ K, converges equally rapidly to
the equilibrium temperature of $\sim 500$ K for $C_{2}=0.025$ eV.Å2.

In Fig. 1 (c) we plot the time evolution of the average per-atom spin and
rotational angular momentum for $C_{2}=0.025$ eV.Å2. Both the rotational
lattice and spin angular momenta are initially near $0$. As the simulation
proceeds, angular momentum is transferred from the spins as they rotate
(during the process of reaching a final equilibrium magnetized state), to
rotational lattice angular momentum, while conserving the total angular
momentum. The magnitude of the total average angular momentum, or red line, in
Fig. 1 (c) has been calculated as

$\begin{gathered}\left|\left<\vec{\lambda}+\vec{S}\right>\right|=\left({[\left<\lambda_{x}\right>+\left<S_{x}\right>]}^{2}+{[\left<\lambda_{y}\right>+\left<S_{y}\right>]}^{2}+{[\left<\lambda_{z}\right>+\left<S_{z}\right>]}^{2}\right)^{1/2}.\end{gathered}$
(23)

We note that the average per atom spin angular momentum of $\sim
0.8\left|\vec{S}\right|_{\rm initial}$ is higher than expected at the
temperature555SPILADY calculates the lattice temperature from the
equipartition of kinetic energy in three dimensions and we implemented the
generalized expression for the spin temperature developed by Nurdin et al.
[43] in SPILADY. of $\sim 500$ K shown in Fig. 1 (d), due to the model of
Refs. [33, 17, 22] which does not account for size dependence of magnetization
in nanoparticles [44, 45].

In Fig. 2 (a) we plot the difference in magnetic anisotropy energy
$\Delta\text{MAE}$ of the spin populations at $90$° to one another as a
function of $C_{2}$; it is seen that $\Delta\text{MAE}$ is linear in $C_{2}$,
verifying that the MAE is directly proportional to $C_{2}$.

Figure 2: (a) The magnetic anisotropy energy difference scales linearly with
$C_{2}$ ($R^{2}=0.993$) and the range of $\Delta\text{MAE}$ values are $\sim
1$ meV/atom [42]. (b) The dependence on $C_{2}$ of the rate of change in
angular momentum is clear, with stronger anisotropic interactions (larger
$C_{2}$) resulting in more rapid change in angular momenta. (Here,
$R^{2}=0.978$ in the case of the spins and $0.988$ in the case of the
lattice.)

Further, we verify that the strength of the anisotropy correction $C_{2}$
determines the rate of exchange of angular momentum between spins and lattice
rotational angular momentum. To this end, and for ease of comparison, we fit
the absolute values of the average per-atom lattice and spin angular momenta
to the hyperbolic tangent function $a\tanh\left(bt\right)$ using the first
$500$ ps of simulation data in each case, as illustrated in Fig. 1 (c) for the
particular case of $C_{2}=0.025$ eV.Å2. The rate of change of angular momentum
is the greatest near the origin; therefore we use the derivatives of the fits
to this function at $t=0$ as measure of the rate of exchange of angular
momentum. In Fig. 2 (b) we plot the rates of change of spin and lattice
rotational angular momenta as a function of $C_{2}$. It is seen that the rate
of exchange of angular momenta is linear in $C_{2}$, verifying that this rate
is directly proportional to the strength of the magnetic anisotropy.

Figure 3: Conservation of total angular momentum for the full simulation
duration and different values of the anisotropy coefficient $C_{2}$: (a)
$0.025$, (b) $0.075$, (c) $0.1$, (d) $0.15$ eV.Å2. Total angular momentum
(black) does not vary during the simulation, while the individual components,
the average per atom spin (green curve) and lattice rotational moment (red
curve) do vary. As the value of $C_{2}$ increases, the coordinates of the
average per-atom spin and lattice angular momenta take less time to approach
the long ($z$) axis of the spheroid via exchange of angular momentum. This is
also evident from the values of the derivatives at $t=0$ of the hyperbolic
tangent fits to the absolute values of the spin and lattice angular momenta in
Fig. 2.

In order to provide visual confirmation of the trend in Fig. 2 (b), we plot in
Fig. 3 the time evolution of the coordinates of the per-atom average angular
momenta of the spins (green), lattice (red) and their sum (black) for a range
of values of $C_{2}$. Since the anisotropy is uniaxial in this work, the green
and red traces tend to cluster near the $z$-axis, and with opposite signs. The
symmetry of the spin Hamiltonian implies that there is no preferred final
direction of the average per-atom spin and lattice angular momenta as long as
they are mutually opposite and generally point along the cluster’s easy axis
($z$ in this case). Therefore, randomly, either of the spin or the lattice
rotational angular momentum can be oriented along the positive $z$ axis (see
e.g. the traces of Fig. 3 (c) vs the other three traces); it is the relative
orientation of rotational versus spin angular momenta which is important, and
which are always opposite and equal, conserving total angular momentum.

Finally, we note that, to our knowledge, the only previous SLD simulation we
are aware of that reproduces the EdH effect (Ref. [9]) applies an external
magnetic field of $50$ Tesla throughout the entire demagnetization process of
a free face-centered cubic Co cluster. Also, only the $z$ components of the
spin and lattice angular momenta are plotted in Ref. [9], casting doubt on
whether total angular momentum is actually conserved in their EdH simulation.
(An external field can induce a non-zero overall spin torque due to individual
contributions not cancelling out in a pairwise manner.)

##  4 Conclusion

We have shown analytically that the dynamic many-body model of uniaxial
anisotropy described in Ref. [17] conserves the total angular momentum and
energy of the system. We also implemented this model in the spin-lattice
dynamics code SPILADY [32], and also make the source code available online
[46]. Our implementation reproduces the full Einstein–de Haas effect in an
isolated iron nanoparticle, consistent with the conservation of total angular
momentum by the model over a range of values of the magnetic anisotropy
coefficient $C_{2}$. The value of this coefficient not only determines the
rate at which exchange between spin and rotational angular momenta occurs, but
also how quickly the angular momenta align parallel to the long axis of the
nanoparticle for positive $C_{2}$.

The advantage of the approach followed in this work, compared to other similar
generally applicable semi-classical models, is that the realistic EAM
potential is used to describe interactions between the atoms, allowing
accurate reproduction of ab-initio free surface energies and defects. Another
advantage is that the angular momentum exchange occurs in a very stable and
consistent manner, with full spin-lattice relaxation occurring in agreement
with previously reported times $\sim 100$ ps. Moreover, neglecting shape
anisotropy, the difference in per-atom magnetic anisotropy energy between spin
populations oriented at $90$° to each other in the same frozen atomic
configuration, yields order-of-magnitude agreement, at the lower range of the
magnetic anisotropy interaction strength, $C_{2}$, with previously reported
tight-binding values for Fe nanoclusters of comparable size. The development
of dynamic anisotropy models due to local symmetry breaking at finite
temperatures is therefore crucial to modeling fast demagnetization processes.

Our implementation of anisotropy can be employed to study any single-element
ferromagnetic system in which local coupling between spins and atoms cannot be
ignored. It is also possible to extend this approach to multiple-element
systems once SPILADY, the underlying spin-lattice dynamics code used in this
implementation, makes provision for this capability. For this reason it is
likely to find applications in a significant number of situations where
uniaxial magnetization is important and coupling between lattice and spins
plays a non-negligible role.

## Acknowledgements

We acknowledge financial support from the Ministry of Science and Innovation
of Spain (grant No. PID2019-109539GB-4). This work was supported by the
Generalitat Valenciana through PROMETEO2017/139 and PROMETEO/2021/017. C. S.
gratefully acknowledges financial supports from Generalitat Valenciana with
(CDEIGENT2018/028). JFR acknowledges funding from FCT grant PTDC/FIS-
MAC/2045/2021. The SLD calculations in this paper were performed on the high-
performance computing facilities of the University of South Africa.

## References

  * [1] A. Einstein, W. Haas, Experimental proof of the existence of ampère’s molecular currents, KNAW, Proceedings, 18 I, 1915, Amsterdam 18 I (1915) 696–711.
URL http://www.dwc.knaw.nl/DL/publications/PU00012546.pdf

  * [2] S. J. Barnett, Magnetization by Rotation, Phys. Rev. 6 (1915) 239.  doi:10.1103/PhysRev.6.239.
  * [3] L. Landau, E. Lifshitz, L. Pitaevskii, Electrodynamics of Continuous Media, Course of Theoretical Physics, Butterworth-Heinemann, Oxford, 1995.
  * [4] G. E. W. Bauer, S. Bretzel, A. Brataas, Y. Tserkovnyak, Nanoscale magnetic heat pumps and engines, Phys. Rev. B 81 (2010) 024427.
  * [5] G. E. Bauer, E. Saitoh, B. J. V. Wees, Spin caloritronics, Nat. Mater. 11 (2012) 391.
  * [6] K. Fukushima, S. Pu, Z. Qiu, Eddy magnetization from the chiral Barnett effect, Phys. Rev. A 99 (2019) 032105.
  * [7] G. G. Scott, Review of gyromagnetic ratio experiments, Rev. Mod. Phys. 34 (1962) 102–109.
  * [8] T. M. Wallis, J. Moreland, P. Kabos, Einstein de Haas effect in a NiFe film deposited on a microcantilever, Appl. Phys. Lett. 89 (2006) 122502.
  * [9] M. Aßmann, U. Nowak, Spin-lattice relaxation beyond gilbert damping, J. Magn. Magn. 469 (2018) 217–223.  doi:10.1016/j.jmmm.2018.08.034.
  * [10] A. A. Kovalev, G. E. W. Bauer, A. Brataas, Current-driven ferromagnetic resonance, mechanical torques, and rotary motion in magnetic nanostructures, Phys. Rev. B 75 (2007) 014430.
  * [11] R. Jaafar, E. M. Chudnovsky, D. A. Garanin, Dynamics of the Einstein-de Haas effect: Application to a magnetic cantilever, Phys. Rev. B 79 (2009) 104410\.
  * [12] E. M. Chudnovsky, D. A. Garanin, Rotational states of a nanomagnet, Phys. Rev. B 81 (2010) 214423.
  * [13] J. H. Mentink, M. I. Katsnelson, M. Lemeshko, Quantum many-body dynamics of the Einstein-de Haas effect, Phys. Rev. B 99 (2019) 064428.
  * [14] N. Bloembergen, S. Shapiro, P. Pershan, J. Artman, Cross-relaxation in spin systems, Physical Review 114 (2) (1959) 445.
  * [15] D. Beaujouan, P. Thibaudeau, C. Barreteau, Anisotropic magnetic molecular dynamics of cobalt nanowires, Phys. Rev. B 86 (2012) 174409\.  doi:10.1103/PhysRevB.86.174409.
URL https://journals.aps.org/prb/abstract/10.1103/PhysRevB.86.174409

  * [16] D. Perera, Atomistic simulations of magnetic models with coupled translational and spin degrees of freedom, Ph.D. thesis, University of Georgia (2015).
  * [17] D. Perera, M. Eisenbach, D. M. Nicholson, G. M. Stocks, D. P. Landau, Reinventing atomistic magnetic simulations with spin-orbit coupling, Phys. Rev. B 93 (2016) 060402\.  doi:10.1103/PhysRevB.93.060402.
URL https://doi.org/10.1103/PhysRevB.93.060402

  * [18] M. Aßmann, Atomistic simulation of spin-lattice dynamics, Ph.D. thesis, University of Konstanz (2018).
  * [19] J. Hellsvik, D. Thonig, K. Modin, D. Iusan, A. Bergman, O. Eriksson, L. Bergqvist, A. Delin, General method for atomistic spin-lattice dynamics with first principles accuracy, Phys. Rev. B 99 (2018) 104302.  doi:10.1103/PhysRevB.99.104302.
  * [20] J. Tranchida, S. Plimpton, P. Thibaudeau, A. Thompson, Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics, J. Comp. Phys. 372 (2018) 406–425.  doi:10.1016/j.jcp.2018.06.042.
  * [21] P. Nieves, J. Tranchida, S. Arapan, D. Legut, Spin-lattice model for cubic crystals, Phys. Rev. B 103 (2021) 094437.  doi:10.1103/PhysRevB.103.094437.
URL https://link.aps.org/doi/10.1103/PhysRevB.103.094437

  * [22] M. Strungaru, M. O. A. Ellis, S. Ruta, O. Chubykalo-Fesenko, R. F. L. Evans, R. W. Chantrell, Spin-lattice dynamics model with angular momentum transfer for canonical and microcanonical ensembles, Phys. Rev. B 103 (2) (2021) 1–13.  doi:10.1103/physrevb.103.024429.
  * [23] M. Eisenbach, D. Perera, D. P. Landau, D. M. Nicholson, J. Yin, G. Brown, Magnetic materials at finite temperatures: Thermodynamics and combined spin and molecular dynamics derived from first principles calculations, J. Phys.: Conf. Ser. 640 (2015) 012019.
  * [24] W. Dednam, C. Sabater, O. Tal, J. J. Palacios, A. E. Botha, M. J. Caturla, Refined electron-spin transport model for single-element ferromagnetic systems: Application to nickel nanocontacts, Phys. Rev. B 102 (2020) 245415.  doi:10.1103/PhysRevB.102.245415.
URL https://link.aps.org/doi/10.1103/PhysRevB.102.245415

  * [25] J. M. Caridad, D. Mccloskey, J. F. Donegan, V. Krstić, Controllable growth of metallic nano-helices at room temperature conditions, Appl. Phys. Lett. 105 (2014) 233114.  doi:10.1063/1.4904091.
  * [26] S. Finizio, M. Foerster, M. Buzzi, B. Krüger, M. Jourdan, C. A. F. Vaz, J. Hockel, T. Miyawaki, A. Tkach, S. Valencia, F. Kronast, G. P. Carman, F. Nolting, M. Kläui, Magnetic anisotropy engineering in thin film ni nanostructures by magnetoelastic coupling, Phys. Rev. Applied 1 (2014) 021001.  doi:10.1103/PhysRevApplied.1.021001.
URL https://link.aps.org/doi/10.1103/PhysRevApplied.1.021001

  * [27] D. B. Reeves, J. B. Weaver, Approaches for modeling magnetic nanoparticle dynamics, Critical Reviews in Biomedical Engineering 42 (1) (2014) 85.
  * [28] J.-Y. Chaulaeau, et al., Electric and antiferromagnetic chiral textures at multiferroic domain walls, Nature Materials 4 (2020) 386–390.  doi:10.1038/s41563-019-0516-z.
  * [29] A. Fert, N. Reyren, V. Cros, Magnetic skyrmions: advances in physics and potential applications, Nature Reviews Materials 2 (2017) 17031.  doi:10.1038/natrevmats.2017.31.
  * [30] J. Zázvorka, F. Jakobs, D. Heinze, N. Keil, S. Kromin, S. Jaiswal, K. Litzius, G. Jakob, P. Virnau, D. Pinna, K. Everschor-Sitte, L. Rózsa, A. Donges, U. Nowak, M. Kläui, Thermal skyrmion diffusion used in a reshuffler device, Nature Nanotechnology 14 (2019) 658–661.  doi:10.1038/s41565-019-0436-8.
  * [31] B. D. Cullity, C. D. Graham, Introduction to magnetic materials, 2nd Edition, Wiley, New Jersey, 2009.
  * [32] P.-W. Ma, S. L. Dudarev, C. H. Woo, Spilady: A parallel cpu and gpu code for spin lattice magnetic molecular dynamics simulations, Comp. Phys. Comm. 207 (2016) 350 – 361.  doi:https://doi.org/10.1016/j.cpc.2016.05.017.
  * [33] P.-W. Ma, C. H. Woo, S. L. Dudarev, Large-scale simulation of the spin-lattice dynamics in ferromagnetic iron, Phys. Rev. B 78 (2) (2008) 024434.  doi:10.1103/PhysRevB.78.024434.
URL http://link.aps.org/doi/10.1103/PhysRevB.78.024434

  * [34] M. S. Daw, M. I. Baskes, Semiempirical, quantum mechanical calculation of hydrogen embrittlement in metals, Phys. Rev. Lett. 50 (17) (1983) 1285.
  * [35] W. Dednam, Atomistic simulations of competing influences on electron transport across metal nanocontacts, Ph.D. thesis, University of Alicante/University of South Africa (2019).
  * [36] M. Krech, A. Bunker, D. Landau, Fast spin dynamics algorithms for classical spin systems, Comput. Phys. Commun. 111 (1) (1998) 1–13.  doi:10.1016/S0010-4655(98)00009-5.
  * [37] D. P. Landau, S.-H. Tsai, M. Krech, A. Bunker, Improved spin dynamics simulations of magnetic excitations, Int. J. Mod. Phys. C 10 (08) (1999) 1541–1551.
  * [38] C. Dornes, Y. Acremann, M. Savoini, M. Kubli, M. Neugebauer, E. Abreu, L. Huber, G. Lantz, C. Vaz, H. Lemke, E. Bothschafter, M. Porer, V. Esposito, L. Rettig, M. Buzzi, A. Alberca, Y. Windsor, P. Beaud, U. Staub, S. Johnson, The ultrafast einstein-de haas effect, Nature 565 (01 2019).  doi:10.1038/s41586-018-0822-7.
  * [39] E. Beaurepaire, J.-C. Merle, A. Daunois, J.-Y. Bigot, Ultrafast spin dynamics in ferromagnetic nickel, Phys. Rev. Lett. 76 (1996) 4250–4253.  doi:10.1103/PhysRevLett.76.4250.
URL https://link.aps.org/doi/10.1103/PhysRevLett.76.4250

  * [40] P.-W. Ma, S. L. Dudarev, C. H. Woo, Spin-lattice-electron dynamics simulations of magnetic materials, Phys. Rev. B 85 (2012) 184301.  doi:10.1103/PhysRevB.85.184301.
URL https://link.aps.org/doi/10.1103/PhysRevB.85.184301

  * [41] L. Malerba, M. C. Marinica, N. Anento, C. Björkas, H. Nguyen, C. Domain, F. Djurabekova, P. Olsson, K. Nordlund, A. Serra, D. Terentyev, F. Willaime, C. S. Becquart, Comparison of empirical interatomic potentials for iron applied to radiation damage studies, J. Nucl. Mater. 406 (1) (2010) 19 – 38.  doi:10.1016/j.jnucmat.2010.05.017.
URL http://www.sciencedirect.com/science/article/pii/S0022311510002291

  * [42] D. Li, A. Smogunov, C. Barreteau, F. Ducastelle, D. Spanjaard, Magnetocrystalline anisotropy energy of Fe (001) and Fe (110) slabs and nanoclusters: A detailed local analysis within a tight-binding model, Physical Review B 88 (21) (2013) 1–9.  doi:10.1103/PhysRevB.88.214413.
  * [43] W. B. Nurdin, K.-D. Schotte, Dynamical temperature for spin systems, Phys. Rev. E 61 (2000) 3579–3582.  doi:10.1103/PhysRevE.61.3579.
URL https://link.aps.org/doi/10.1103/PhysRevE.61.3579

  * [44] G. Dos Santos, R. Aparicio, D. Linares, E. N. Miranda, J. Tranchida, G. M. Pastor, E. M. Bringa, Size- and temperature-dependent magnetization of iron nanoclusters, Phys. Rev. B 102 (2020) 184426.  doi:10.1103/PhysRevB.102.184426.
URL https://link.aps.org/doi/10.1103/PhysRevB.102.184426

  * [45] G. Dos Santos, R. Meyer, R. Aparicio, J. Tranchida, E. M. Bringa, H. M. Urbassek, Spin-lattice dynamics of surface vs core magnetization in Fe nanoparticles, Applied Physics Letters 119 (16) (2021) 012404.  doi:10.1063/5.0055606.
URL https://aip.scitation.org/doi/10.1063/5.0055606

  * [46] P.-W. Ma, et al., computer code SPILADY, with anisotropy corrections made by W. Dednam (2018).
URL https://github.com/wdednam/spilady/tree/anisotropy