# Coarse graining a non-Markovian collisional model

Nadja K. Bernardes1111Corresponding author<EMAIL_ADDRESS>Andre R.
R. Carvalho2, C. H. Monken1, and Marcelo F. Santos3 1Departamento de Física,
Universidade Federal de Minas Gerais, Belo Horizonte, Caixa Postal 702,
30161-970, Brazil  2Centre for Quantum Computation and Communication
Technology, Department of Quantum Sciences, Research School of Physics and
Engineering, The Australian National University, Canberra, ACT 0200 Australia
3Instituto de Física, Universidade Federal do Rio de Janeiro, Caixa Postal
68528, Rio de Janeiro, 21941-972, Brazil

###### Abstract

The dynamics of systems subjected to noise is called Markovian in the absence
of memory effects, i.e. when its immediate future only depends on its present.
Time correlations in the noise source may generate non-Markovian effects that,
sometimes, can be erased by appropriately coarse graining the time evolution
of the system. In general, the coarse graining time $t_{\text{CG}}$ is taken
to be much larger than the correlation time $\tau$ but no direct relation
between them is established. Here we analytically obtain a relation between
$t_{\text{CG}}$ and $\tau$ for the dynamics of a qubit subjected to a time
correlated environment. Our results can be applied in principle to any
distribution of the environmental correlations and can be tested through a
collisional model where the qubit sequentially interacts with correlated
qutrits.

##  I Introduction

In open quantum systems two important approximations are made in order to
derive a master equation: the Born approximation, which assumes a weak system-
environment interaction, and the Markov approximation, which neglects
environmental memory effects and is justified when the time-scale
characterizing the decay of the environmental correlations is very short in
comparison to all the other significant time-scales of the problem breuer . It
is worth noticing, however, that completely uncorrelated environments are
never found in practice. A common procedure to achieve a Markovian master
equation then is to coarse grain time, i.e. the evolution of the system is
considered in time steps that are much larger than any characteristic time
scale of the environment lidar ; lidar2 , but still short enough when compared
to the time scale of changes in the state of the system. In this case, time is
discretized in intervals of finite size $t_{\text{CG}}$ and quantum master
equations in the continuous limit, such as in the Lindblad form lindblad ,
have to be considered as an approximation that is good enough if the system is
not probed in intervals smaller than $t_{\text{CG}}$. This restriction is
reasonable in many cases and such equations describe a variety of processes
like the spontaneous emission of a two-level atom or the damping of an
harmonic oscillator breuer ; wiseman .

In general, $t_{\text{CG}}$ is assumed to be much larger than the correlation
time of the environment $\tau$, but no specific information about what “much
larger” actually means is provided. From the experimental point of view, this
procedure is partially supported by the fact that any measurement apparatus is
limited by a maximum data acquisition rate. But, fundamentally, it means that
if one can probe the system fast enough, one should expect to observe non-
Markovian effects in its dynamics. This is the basis for the Zeno effect, for
example zeno . A natural question, then, regards the limits in time scales
within which the non-Markovian dynamics of a given system can be
experimentally elucidated. A more detailed relation between $t_{\text{CG}}$
and $\tau$, however, strongly depends on the specificity of the system and its
interaction with the environment and is usually difficult to achieve. It is
even more difficult to separate the role played by the correlations in the
environment and other non-Markovianity sources such as strong system-
environment interactions or the internal dynamics of the environment. In
general, all theses factors bunch together when $t_{\text{CG}}$ is established
for a given dynamics.

Recent theoretical ziman1 ; ziman2 ; giovannetti ; tomas ; ciccarello ;
vacchini ; budini ; Paternostro ; nadja  and experimental Liu ; Liu2 ; tang ;
Steve1 ; Fabio2 ; chiuri ; Xu ; Fabio1 ; nadja2 ; nadja3  works have
demonstrated how to use the so-called collisional models to simulate generic
non-unitary dynamics ziman1 ; ziman2  and, specifically, to better understand
the emergence of non-Markovian effects giovannetti ; tomas ; ciccarello ;
vacchini ; budini ; Paternostro ; nadja . Collisional models are those in
which the environment is made of particles that randomly collide with the
system, usually one at a time, generating non-unitary dynamics. In particular,
in nadja , we have designed a model where the non-Markovianity caused by
correlations in the environment can be investigated independently from the
other above mentioned factors.

In this paper, we use the results found in nadja  to show how to recover a
Markovian dynamics of a qubit by coarse graining its time evolution. In
particular, we identify in our model what “much larger” means by presenting a
clear relation between the coarse graining time and the correlation time of
the environment. We analytically derive this relation for arbitrary
correlation functions of the environment and test our results for two
particular types: a step function and an exponential decay of environment
correlations.

In Sec. II, we present our collisional model. The correlation functions are
defined in Sec. III. The dynamics is analyzed in Sec. IV and the relation
between the correlation time and the the coarse grain is established. The
limit to the continuous is established in Sec. V. We conclude in Sec. VI.

##  II The model

In classical physics, the master equation is used to describe the time-
evolution of the probability that a system undergoing stochastic evolution is
in a particular state at a given time. The quantum mechanical version of this
type of evolution for a given physical state $\rho(t)$ is given by the so-
called quantum master equation in the Lindblad form

$\displaystyle\frac{d\rho}{dt}=-i[H_{s}/\hbar,\rho]+\sum_{i}\mathcal{L}_{i}(\rho),$
(1)

where $H_{s}$ is the system Hamiltonian and
$\mathcal{L}_{i}(\rho)=-\frac{\gamma_{i}}{2}(\Gamma^{\dagger}_{i}\Gamma_{i}\rho+\rho\Gamma^{\dagger}_{i}\Gamma_{i}-2\Gamma_{i}\rho\Gamma^{\dagger}_{i})$
gives the non-unitary contribution due to the interaction with the environment
(responsible for the so-called decoherence of $\rho$). This equation can be
rewritten in time steps $dt$ as
$\rho(t+dt)=\sum_{j}W_{j}\rho(t)W^{\dagger}_{j}$ where
$W_{0}=\mathds{1}-dt(iH_{s}/\hbar+\sum_{i}\frac{\gamma_{i}}{2}\Gamma^{\dagger}_{i}\Gamma_{i})$,
$W_{i\neq 0}=\sqrt{dt\gamma_{i}}\Gamma_{i}$, $\gamma_{i}dt\ll 1$ and terms of
$O(dt^{2})$ are discarded. In order to guarantee that $\rho$ is a quantum
state at any given time, $\gamma_{i}$ has to be a positive rate for all $i$,
and $\sum_{j}W^{\dagger}_{j}W_{j}=\mathds{1}$. In this way, the trace and
positivity of $\rho$ are always preserved. Note that if we assume that $H_{s}$
is equal to zero and$\rho$ describes the state of a qubit under a random
unitary evolution with $\Gamma_{i}\equiv\sigma_{i}$ Pauli-like matrices
($\Gamma_{i}^{\dagger}=\Gamma_{i}$ and $\Gamma_{i}^{2}=1$), then Eq (1) can be
rewritten as the Completely Positive and Trace Preserving (CPTP) map

$\displaystyle\Lambda_{n+1,n}:\rho_{n}\mapsto\rho_{(n+1)}$ $\displaystyle=$
$\displaystyle\Lambda_{n+1,n}(\rho_{n})$ $\displaystyle=$
$\displaystyle(1-\sum_{i=1,2}\epsilon_{i})\rho_{n}+\sum_{i=1,2}\epsilon_{i}\sigma_{i}\rho_{n}\sigma_{i}.$

where $\epsilon_{i}\equiv\gamma_{i}dt$ gives the probability of a respective
change $\sigma_{i}$ in the system in an arbitrary time interval
$\\{ndt,(n+1)dt\\}$. We will also consider only two channels, $\sigma_{1}$ and
$\sigma_{2}$ with no loss of generality, as previously shown in nadja . This
is the scenario to be explored from now on.

This kind of time evolution, specially the discretized time version described
in Eq.(II), has been associated to collisional models Plenio3 . In this case,
each time step $n\rightarrow n+1$ corresponds to a collision with the
reservoir and the Markovian hypothesis is that the map after $m$ collisions,
that takes the system from $n$ to $n+m$, can be concatenated as the product of
$m$ intermediate CPTP maps:
$\Lambda_{n+m,n}=\Lambda_{n+m,n+m-1}\Lambda_{n+m-1,n+m-2}...\Lambda_{n+1,n}$,
for arbitrary $n$ and $m$. However, as shown in nadja , correlations in the
environmental state may lead to time dependent rates when the probabilities
$\epsilon_{i}$ depend on $n$. That happens even if the transformations
undergone by the system are equal ($\sigma_{1}=\sigma_{2}$). If, on the other
hand, $\sigma_{1}\neq\sigma_{2}$, environmental correlations also produce
“artificial” decoherence channels (extra $W$ operators) that do not derive
directly from the collisions but from their correlation, and, more important,
that may be associated to negative rates $\gamma_{i}<0$ nadja . In this case,
Eq. (II) does not describe correctly the physical process ($\rho_{n+1}$ cannot
be described as CP map applied to $\rho_{n}$) and the dynamics becomes non-
Markovian, i.e. the state of the qubit after a given collision becomes
dependent not only on its state immediately before it but also on its state at
earlier times.

In nadja  we analyzed the effects of correlations in just two consecutive
collisions because we were interested in understanding the minimum effects
that environmental correlations could produce. Physically, this situation
corresponds to a very short correlation pulse in an otherwise uncorrelated
environment. In the present work we are interested in a more realistic
scenario for the environment so that we can address the opposite question:
given that most physical environments are always somewhat correlated and that
these correlations may generate non-Markovian evolutions, we investigate the
minimum coarse graining in time that one needs to define in this model so that
any observation on the system concludes for a Markovian time evolution in the
Lindblad form. Note that to coarse grain time in our model corresponds to
allow for a certain number of collisions to occur in between any consecutive
observations of the system. In terms of the quantities already defined, it
means to establish a minimum interval of collisions $n_{\text{CG}}$ so that
$\rho_{(n+n_{\text{CG}})}=\sum_{j}W^{\prime}_{j}\rho_{n}W^{\prime\dagger}_{j}$,
for any “time” $n$, where the set $\\{W^{\prime}_{j}\\}$ shares the same
general properties of the previous $\\{W_{j}\\}$ set even though the operators
and rates will not necessarily be the same.

##  III Correlation functions

In order to better understand how environment correlation functions are
introduced in our model, it is useful to visualize its uncorrelated version in
terms of quantum trajectories. If the qubit is initially prepared in a pure
state $|\Psi_{0}\rangle$, then a quantum trajectory $Q$ that takes the system
from its initial state to its state $|\Psi_{n}\rangle$ after $n$ collisions
($T=ndt$) corresponds to a sequence
$\\{|\Psi_{0}\rangle,|\Psi_{1}\rangle,|\Psi_{2}\rangle,...,|\Psi_{n}\rangle\\}$,
where $|\Psi_{k}\rangle=U^{(k)}|\Psi_{k-1}\rangle$ and $U^{(k)}$ is drawn from
the set $\\{\sigma_{0},\sigma_{1},\sigma_{2}\\}$ ($\sigma_{0}=\mathds{1}$)
with corresponding probabilities $\\{p_{k}(0),p_{k}(1),p_{k}(2)\\}$ given by
$\\{1-(\epsilon_{1}+\epsilon_{2}),\epsilon_{1},\epsilon_{2}\\}$. A trajectory
$Q$ from $0$ to $n$ is then defined by the sequence
$Q|\Psi\rangle=\Pi_{k=1}^{n}U^{(k)}|\Psi_{0}\rangle$ which happens with
probability $P_{Q}=\Pi_{k=1}^{n}p_{k}$, where $p_{k}$ is the probability
associated to each particular realization of $U^{(k)}$. The state of the
system at $n$ is given by averaging over all possible trajectories
$\rho_{n}=\sum_{Q}Q\rho_{0}Q^{\dagger}$ and, more important, Eq. (II) holds
for any “time” $n$ when $\epsilon_{i}\ll 1$.

Correlations in the reservoir are introduced by modifying the condition just
described, i.e. the joint probabilities $p_{kk^{\prime}}(i,j)$ of two (or
more) jumps of the same type become more likely to happen then for
independent, random events: $p_{kk^{\prime}}(i,i)\neq\epsilon_{i}^{2}$, where
$p_{kk^{\prime}}(i,j)$ is the joint probability that $i$-th jump occurred in
the $k$-th collision and the $j$-th jump occurred in the $k^{\prime}$-th
collision. In this work, we analyze two models of correlation functions for
the environment: a step function and an exponential time decay. In both cases
the correlation time of the reservoir is well defined as we explain later. For
the step function, we choose a certain number of collisions for which once the
first jump is of a given type, say $\sigma_{1}$, then any other jump in the
same interval must be of the same type ($p_{kk^{\prime}}(1,2)=0$ within the
interval). In this scenario, this interval is equivalent to the correlation
time of the reservoir from now on denoted by $n_{\text{cor}}$. This choice
naturally modifies the probabilities in each step and clearly introduces
memory in the dynamics. First note that since $\epsilon_{i}\ll 1$ the most
likely operation in each time step still is the no-jump $W_{0}$ associated to
applying identity to the state of the system. Second, note that the counting
of the correlated time interval always starts at the first jump. From then on
and for the duration of $n_{\text{cor}}$, for that particular trajectory any
other possible jump has to be the same.

The exponential decay model follows the same logic, the difference being that
now once a jump happens, the probability for the same type of jump to happen
in the near future decays exponentially as a function of time (number of
collisions), with the joint probability given by:
$p_{kk^{\prime}}(i,i)\propto\epsilon^{2}(1+e^{\frac{-(k^{\prime}-k)}{n_{\text{cor}}}})$,
where $i=1,2$. In this case, $1/n_{\text{cor}}$ defines the decay rate of the
reservoir correlations. In both cases, the effects of the reservoir
correlations in the dynamics of the system are analyzed through calculating
the map $\Lambda_{(n+m),n}$ that takes $\rho_{n}$ into $\rho_{(n+m)}$, for
arbitrary $n$ and number of collisions $m$, and testing the eventual
deviations from the map shown in Eq. (II).

Following the results of nadja , it suffices to investigate the configuration
in which the probabilities of the two different channels $\sigma_{1}$ and
$\sigma_{2}$ are equal ($\epsilon_{1}=\epsilon_{2}=\epsilon$), with
$\sigma_{1}=\sigma_{x}$ and $\sigma_{2}=\sigma_{z}$ (corresponding to bit flip
and dephasing). This type of decoherence is enough to encompass most of the
interesting effects the model can produce; any other configuration implies in
less non-Markovianity per time step. To illustrate the functions and confirm
that the numerical simulation is in agreement with the theoretical model, we
plotted in Fig. 1 the average correlation function $\left<\Gamma(h)\right>$
for different trajectories, where
$\Gamma(h)=\frac{\sum_{i=1}^{n-h}x_{i}x_{i+h}}{\sum_{i=1}^{n-h}|x_{i}x_{i+h}|}$
with $x_{i}=1$(-1) if the collision is $\sigma_{x}$($\sigma_{z}$).

(a)  (b)

Figure 1: Plot of the average correlation function $\left<\Gamma(h)\right>$
for $10^{7}$ trajectories and the total number of collisions of $n=100$ for
(a) the step function and (b) the exponential decay. The red line represents a
fit, where the estimated correlation length is $\bar{n}_{\text{cor}}=10.009$
in (a) and $\bar{n}_{\text{cor}}=10.28$ in (b) in agreement with the expected
value of $n_{\text{cor}}=10$.

##  IV Dynamics analysis

When $\sigma_{1}=\sigma_{x}$, $\sigma_{2}=\sigma_{z}$ and
$\epsilon_{1}=\epsilon_{2}=\epsilon$, correlations in the reservoir end up
generating non-Markovian evolutions that can be evidenced by the non-complete-
positivity divisibility of the maps nadja , i.e. CP maps for the overall
evolution over $m$ collisions that may not be decomposable into the product of
CP maps for the intermediate evolutions,
$\Lambda_{m+i,i}\neq\Lambda_{m+i,n+i}\Lambda_{n+i,i}$ ($n<m$). In this case,
finding the coarse graining time means finding the minimum $n$ beyond which
the decomposition into CP maps is always possible. In order to find
$n_{\text{CG}}$ for each model of correlation in the reservoir, we test the
maps $\Lambda_{m+i,i}$ that take $\rho_{i}$ into $\rho_{i+m}$ for different
intervals $m$. Testing the complete positivity of a map is equivalent to
testing whether the Choi matrix (also known as the dynamical matrix)
associated to this map is positive semidefinite, i.e. it has nonnegative
eigenvalues choi ; bengtsson .

(a)  (b)

Figure 2: Plot of the smallest eigenvalue $\lambda_{2n,n}$ of the dynamical
matrix of the map $\Lambda_{2n,n}$ versus the number of collisions $n$, for
$\epsilon=0.001$, $10^{8}$ trajectories and for (a) the step function and (b)
the exponential decay function. Different lines represent different
correlation lengths as it is shown in the plot. The lines are just for
illustrative reasons, since the number of collisions is a discretized
quantity.

We will analyze the case where $i=0$ and $m=2n$. In Fig. 2 the smallest
eigenvalue $\lambda_{2n,n}$ of the dynamical matrix corresponding to the map
$\Lambda_{2n,n}$ is plotted versus the number of collisions $n$ for the step
function and for the exponential decay function. If this eigenvalue is
positive, this means that $\Lambda_{2n,0}$ is divisible in
$\Lambda_{2n,0}=\Lambda_{2n,n}\Lambda_{n,0}$. We can identify here what are
the points where the curves cross the $\lambda=0$ axis, i.e. it assumes
positive values again, and consequently identify a coarse graining length
$n_{\text{CG}}$ as a function of $n_{\text{cor}}$.

By identifying a coarse grain $n_{\text{CG}}$ for the step and exponential
correlation function, we are able to show the dependence of the coarse grain
with the correlation length as presented in Fig. 3. It is possible to see that
the size of the coarse grain is increasing almost linearly with the
correlation length by a factor of approximately $2$.

This result can be explained by analyzing the maps of the resulting dynamics.
The map from one “time” $i$ to $j$ is given by
$\Lambda_{j,i}(\rho(i))=\sum_{\alpha=0,x,y,z}w_{\alpha}(j,i)\sigma_{\alpha}\rho(i)\sigma_{\alpha}$.
Note here that although in one collision, the system may undergo only a
$\sigma_{x}$ or a $\sigma_{z}$ operation, with more collisions a $\sigma_{y}$
operation may result from the combination of $\sigma_{x}$ and $\sigma_{z}$. If
the environment is uncorrelated, and to second order in $\epsilon$, this
somewhat artificial $\sigma_{y}$ channel, that is created by bunching
collisions, is proportional to the chance $\epsilon^{2}$ of having two
collisions of different types in that particular bunch (always remembering
that in each time step, the most probable collision is of the
$\sigma_{0}=\mathds{1}$ type). Truncating in $O(\epsilon^{2})$ is justifiable
for $\epsilon\ll 1$ (we will comment more on this later). In fact, for the
uncorrelated case, the probability of such $\sigma_{y}$ channel is the same as
that of an extra $\sigma_{0}$ contribution to the dynamics, since the chance
of similar or opposite collisions of type $\sigma_{x}$ and $\sigma_{z}$ is the
same. However, when correlations in the environment are introduced, these
probabilities become unbalanced and the extra $\sigma_{0}$ channel becomes
more probable then the $\sigma_{y}$ one. This unbalance is the origin of the
non-Markovian effect and the “lack” of some portion of the $\sigma_{y}$
channel reflects on a negative weight $w_{y}$ in any attempt to force the
overall map corresponding to that bunch to be a concatenation of maps, each
one representing a single collision. However, due to the finitude of
$n_{\text{cor}}$, if one takes a large enough interval
$n_{\text{CG}}>n_{\text{cor}}$, the uncorrelated events in this interval
compensate for the memory effects brought by the correlated collisions and the
dynamics coarse grained in $n_{\text{CG}}$ intervals become Markovian.
Therefore, in order to find $n_{\text{CG}}$, one needs to check the minimum
time interval beyond which all the weights, including that of $\sigma_{y}$,
become nonnegative. In this model, since the map is invertible, this is
equivalent to check if $\Lambda_{2n,n}=\Lambda_{2n,0}\Lambda_{n,0}^{-1}$ is a
CPTP map, which already imposes constraints over the weights
$w_{\alpha}(2n,n)$, i.e. $w_{\alpha}(2n,n)\geq 0$. This is exactly what was
numerically done and presented in Fig.3.

However, for small enough $\epsilon$ (in our calculations $\epsilon=0.001$),
and as long as the bunching of $M$ collisions is not too large, i.e.
$M\epsilon$ is still much smaller than one, the chance of three or more
collisions of type $\sigma_{x}$ or $\sigma_{z}$ within this bunch is very low
and one can assume at most two collisions of these types with no significant
loss of generality (in our numerics, out of $10^{8}$ generated trajectories,
less than 5% would present three or more of such collisions). This allows us
to analytically calculate a very good approximation for the weights
$w_{\alpha}$, one that becomes even better in the continuous limit of
$\epsilon\rightarrow 0$. First, let us define the weights as:

$\displaystyle w(n,0)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\binom{n}{1}p(1-p)^{n-1},$ (3) $\displaystyle
w_{y}(n,0)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\binom{n}{2}p^{2}(1-p)^{n-2}$ (4)
$\displaystyle\times\left[1-\sum_{k=1}^{n-1}f(k,n_{\text{cor}})G(k,n,p)\right],$

where $p=2\epsilon$ is the probability that in one collision the system will
undergo a $\sigma_{x}$ or $\sigma_{z}$ rotation for an uncorrelated
environment and we take $w_{x}=w_{z}=w$ since both channels are symmetric in
our model. $f(k,n_{\text{cor}})$ is the correlation function of the
environment and for the examples studied in this paper, step or exponential
decaying function, it assumes the respective forms

$\displaystyle f_{\text{step}}(k,n_{\text{cor}})$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{c}1\quad\text{if}\quad k\leq
n_{\text{cor}},\\\ 0\quad\text{otherwise},\end{array}\right.$ (7)
$\displaystyle f_{\text{exp}}(k,n_{\text{cor}})$ $\displaystyle=$
$\displaystyle\text{e}^{-\frac{k}{n_{\text{cor}}}}.$ (8)

Note, however, that at this point the model is generic enough to consider any
function $f(k,n_{\text{cor}})$.

The function $G(k,n,p)$ is related to the probability that two collisions of
the type $\sigma_{x}$ or $\sigma_{z}$ occur separated by $k$ collisions in an
interval of $n$ collisions and it is given by

$G(k,n,p)=\frac{p^{2}(n-k)(1-p)^{k-1}}{(1-p)^{n}+np-1}.$ (9)

The intermediate $\Lambda_{2n,n}$ map will be given
$\Lambda_{2n,n}(\rho(n))=\sum_{\alpha=0,x,y,z}w_{\alpha}(2n,n)\sigma_{\alpha}\rho(n)\sigma_{\alpha}$,
where the coefficients are given in terms of the coefficients from Eq. (3,4)
as

$\displaystyle w_{0}(2n,n)$ $\displaystyle=$ $\displaystyle
1-2w(2n,n)-w_{y}(2n,n),$ (10) $\displaystyle w(2n,n)$ $\displaystyle=$
$\displaystyle\frac{w(2n,0)-w(n,0)}{1-4w(n,0)},$ (11) $\displaystyle
w_{y}(2n,n)$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left[1+\frac{1-4w(2n,0)}{1-4w(n,0)}\right.$
$\displaystyle-\left.\frac{2-4w(2n,0)-4w_{y}(2n,0)}{1-2w(n,0)-2w_{y}(n,0)}\right],$

where, once again, we have considered $w_{x}(j,i)=w_{z}(j,i)=w(j,i)$ due to
the symmetry of $\sigma_{x}$ and $\sigma_{z}$ channels.

For the reasons previously discussed regarding the origin of the $\sigma_{y}$
channel or a simple inspection of Eqs. (10-IV), $w_{y}$ is the only weight
that can assume negative values. Thus, the map $\Lambda_{2n,0}$ is divisible,
if $w_{y}(2n,n)\geq 0$. Note that this corresponds exactly to verifying when
the eigenvalues of the dynamical matrix are all positive. In fact, it is
possible to show that $w_{y}=2\lambda_{y}$ and, for small $\epsilon$, the
inequality $\lambda_{y}\geq 0$ recovers the result observed in the numerical
simulations for the coarse grain time $n_{\text{CG}}$. For example, for the
step function, if one inserts Eqs. (7,9) into Eqs. (10-IV) and then expands
the latters keeping terms up to $O(\epsilon^{2})$, the corresponding
inequality for $w_{y}(2n,n)$ simplifies to

$\displaystyle n(n-1-2n_{\text{cor}})>0,$ (13)

which has a solution for $n>2n_{\text{cor}}+1$. A similar result can be found
for the exponential function. In this way, we confirmed the previous result
encountered in the numerical simulation presented in Fig. 3.

Figure 3: Plot of the coarse grain $n_{\text{CG}}$ versus the correlation
length $n_{\text{cor}}$ for the exponential function (dots) and the step
function (squares) for $\epsilon=0.001$ and $10^{8}$ trajectories.

##  V Continuous limit

So far, we have considered a discretized model of time steps of finite size
proportional to $\epsilon$. Note, however, that in our model $\epsilon$ can be
taken as small as desired, which means that we can take the model to the
continuous limit. In order to do this, we make the following change of
variables in Eqs. (7-9): $n=T/\Delta$, $k=t/\Delta$,
$n_{\text{cor}}=\tau/\Delta$, and $\epsilon=\gamma\Delta$. The functions then
change to:

$\displaystyle f_{\text{step}}(t,\tau)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{c}1\quad\text{if}\quad t\leq\tau,\\\
0\quad\text{otherwise},\end{array}\right.$ (16) $\displaystyle
f_{\text{exp}}(t,\tau)$ $\displaystyle=$
$\displaystyle\text{e}^{-\frac{t}{\tau}},$ (17) $\displaystyle
G\left(t,T,\gamma\right)$ $\displaystyle=$
$\displaystyle\frac{4\gamma^{2}(T-t)(1-2\gamma\Delta)^{\frac{t}{\Delta}-1}\Delta}{(1-2\gamma\Delta)^{\frac{T}{\Delta}}+2\gamma
T-1}.$ (18)

In this way, with fixed $\gamma$, when $\Delta\rightarrow 0$, it follows that
$n\rightarrow\infty$ and $\epsilon\rightarrow 0$, but we guarantee that
$n\epsilon=\gamma T$. As in the Riemman-Steltjes sum, the limit of the sum
will be given by

$\lim_{\Delta\rightarrow
0}\sum_{t=\Delta}^{T-\Delta}f\left(\frac{t}{\Delta},\frac{\tau}{\Delta}\right)G\left(\frac{t}{\Delta},\frac{T}{\Delta},\gamma\Delta\right)=\int_{0}^{T}f\left(t,\tau\right)g(t,T,\gamma)dt,$
(19)

where

$g\left(t,T,\gamma\right)=\frac{4\gamma^{2}(T-t)\text{e}^{-2\gamma
t}}{\text{e}^{-2\gamma T}+2\gamma T-1}.$ (20)

Eq. (19) does not depend on the particular form of the correlation function
and, therefore, this collisional model (determined by the value of $\epsilon$)
indeed stroboscopically simulates a continuous-time quantum evolution for a
variety of different reservoir correlation functions.

##  VI Conclusions

In this work we have used a previously designed collisional model for the non-
Markovian evolution of a qubit to show that by coarse graining its dynamics in
time, i.e. defining minimal time intervals to probe the system, a Markovian
evolution can be recovered. We have also shown that our model allows one to
properly identify a relation between the correlation time of the environment
causing the non-unitary (and non-Markovian) evolution of the qubit and the
coarse graining time that restores its Markovianity. We test the model for two
types of correlation, a step function and an exponential decay and we show
that, surprisingly, depending on the correlation function, the ratio between
the coarse graining time and the correlation time, which is usually taken to
be much larger than one, can be as small as two. The model is general enough
to encompass different correlation functions, hence a broad variety of
environments, and we have derived an approximate analytical solution that fits
well our numerical results, including a limit to the continuous. Finally, the
results presented here could be also useful to identify how broad the spectral
response of a detector should be in order that memory effects caused by the
environment can be neglected.

### Acknowledgments

N.K.B, C.H.M. and M.F.S. would like to thank the support from the Brazilian
agencies CNPq and CAPES. M.F.S. would like to thank the support of FAPEMIG,
project PPM IV. This work is part of the INCT-IQ from CNPq and also of the
Australian Research Council Centre of Excellence for Quantum Computation and
Communication Technology (project number CE110001027).

## References

  * (1) see H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, New York, 2002).
  * (2) D. Bacon, D. A. Lidar, K. B. Whaley, Phys. Rev. A 60, 1944 (1999).
  * (3) D. A. Lidar, A. Shabani, R. Alicki, Chem. Phys. 268, 35 (2001).
  * (4) G. Lindblad, Commun. Math. Phys. 48, 119 (1976).
  * (5) see H. M. Wiseman and G. J. Milburn, Quantum Measurement and Control (Cambridge University Press, Cambridge, 2010).
  * (6) W. M. Itano, D. J. Heinzen, J. J. Bollinger, and D. J. Wineland, Phys. Rev. A 41, 2295 (1990).
  * (7) M. Ziman and V. Buzek, Phys. Rev. A 72, 022110 (2005).
  * (8) M. Ziman, P. Stelmachovic and V. Buzek, Open Syst. Inf. Dyn. 12, 81 (2005).
  * (9) V. Giovannetti and G. M. Palma, Phys. Rev. Lett. 108, 040401 (2012).
  * (10) T. Rybár, S. N. Filipov, M. Ziman, and V. Buzek, J. Phys. B: At. Mol. Opt. Phys. 45, 154006 (2012).
  * (11) F. Ciccarello, G. M. Palma, and V. Giovannetti, Phys. Rev. A 87, 040103(R) (2013).
  * (12) B. Vacchini, Phys. Rev. A 87, 030101(R) (2013).
  * (13) A. A. Budini, Phys. Rev. A 88, 032115 (2013).
  * (14) R. McCloskey and M. Paternostro, Phys. Rev. A 89, 052120 (2014).
  * (15) N. K. Bernardes, A. R. R. Carvalho, C. H. Monken, and M. F. Santos, Phys. Rev. A 90, 032111 (2014).
  * (16) B.-H. Liu, et al, Nat. Phys. 7, 931 (2011).
  * (17) B.-H. Liu, et al, Sci. Rep. 3, 1781 (2013).
  * (18) J.-S. T Tang, et al, Europ. Phys. Lett. 97, 10002 (2012).
  * (19) F. Fanchini, et al, Phys. Rev. Lett. 112, 210402 (2014).
  * (20) J. Jin, et al, Phys. Rev. A 91, 012122 (2015).
  * (21) A. Chiuri, C. Greganti, L. Mazzola, M. Paternostro, and P. Mataloni, Sci. Rep. 2, 968 (2012).
  * (22) J.-S. Xu, et al, Nat. Comm. 4, 2851 (2013).
  * (23) A. Orieux, et al, Sci. Rep. 5, 8575 (2015).
  * (24) N. K. Bernardes, et al, Sci. Rep. 5, 17520 (2015).
  * (25) N. K. Bernardes, et al, Sci. Rep. 6, 33945 (2016).
  * (26) R. Uzdin and R. Kosloff, New J. of Phys. 16, 095003 (2014).
  * (27) S. Lorenzo, R. McCloskey, F. Ciccarello, M. Paternostro, and G. M. Palma, Phys. Rev. Lett. 115, 120403 (2015).
  * (28) A. Rivas, S. F. Huelga, and M. B. Plenio, Rep. Prog. Phys. 77, 094001 (2014).
  * (29) M. M. Wolf and J. I. Cirac, Commun. Math. Phys. 279, 147 (2008).
  * (30) M. M. Wolf, J. Eisert, T. S. Cubitt, and J. I. Cirac, Phys. Rev. Lett. 101, 150402 (2008).
  * (31) A. Rivas, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 105, 050403 (2010).
  * (32) M.-D. Choi, Linear Algebra Appl. 10, 285 (1975).
  * (33) I. Bengtsson and K. Zyczkowski, Geometry of Quantum States: An Introduction to Quantum Entanglement (Cambridge University Press, Cambridge, UK, 2006), Sec. 10.3.


# Nonlinear features of the superconductor–ferromagnet–superconductor
$\varphi_{0}$ Josephson junction in ferromagnetic resonance region

Aliasghar Janalizadeh1 Ilhom R. Rahmonov2,3,4 Sara A. Abdelmoneim 5 Yury M.
Shukrinov2,3,4 Mohammad R. Kolahchi1 1 Department of Physics, Institute for
Advanced Studies in Basic Sciences (IASBS), P.O. Box 45137-66731, Zanjan, Iran
2 BLTP, JINR, Dubna, Moscow Region, 141980, Russia
3 Dubna State University, Dubna, 141980, Russia
4 Moscow Institute of Physics and Technology, Dolgoprudny, 141700, Moscow
Region, Russia
5 Physics department, Menofiya University, Faculty of Science, 32511, Shebin
Elkom,Egypt

(August 28, 2024)

###### Abstract

We demonstrate the manifestations of the nonlinear features in magnetic
dynamics and IV-characteristics of the $\varphi_{0}$ Josephson junction in the
ferromagnetic resonance region. We show that at small values of system
parameters, namely, damping, spin-orbit interaction, and Josephson to magnetic
energy ratio, the magnetic dynamics is reduced to the dynamics of the scalar
Duffing oscillator, driven by the Josephson oscillations. The role of
increasing superconducting current in the resonance region is clarified.
Shifting of the ferromagnetic resonant frequency and the reversal of its
damping dependence due to nonlinearity are demonstrated by the full Landau-
Lifshitz-Gilbert-Josephson system of equations, and in its different
approximations. Finally, we demonstrate the negative differential resistance
in the IV–characteristics, and its correlation with the foldover effect.

Josephson junction, Landau-Lifshitz-Gilbert equation, Duffing oscillator

##  I I. Introduction

The coupling of superconducting phase difference with magnetic moment of
ferromagnet in the $\varphi_{0}$ junction leads to a number of unique features
important for superconducting spintronics, and modern information technology
[1, 2, 3, 4, 5]. It allows to control the magnetization precession by
superconducting current and affects the current–voltage (IV) characteristics
by magnetic dynamics in the ferromagnet, in particular, to create a DC
component in the superconducting current [6, 7, 8]. A remarkable manifestation
of such coupling is the possibility to stimulate a magnetization reversal in
the ferromagnetic layer by applying current pulse through the
$\varphi_{0}$-junction [3, 9, 10, 11, 12, 13].

There are two features of our Josephson junction that come into play in our
study. One is the broken inversion symmetry in the weak link of the Josephson
junction, when the link is magnetic, which introduces an extra phase in the
current—-phase relation, preventing it from being antisymmetric. Such
Josephson junctions are named $\varphi_{0}$ junctions [1], and examples exist
such as MnSi and FeGe. Second is the nonlinear property of the system that
makes for an anomalous resonance behavior [14].

We couple such a Josephson junction to the model that describes the
magnetodynamics in thin films or heterostructure, to form the Landau-Lifshitz-
Gilbert-Josephson model (LLGJ)[15, 16, 14]. It is shown that for a particular
set of parameters, the coupled equations reduce to the dynamics of a Duffing
oscillator [14]. The cubic nonlinearity in this oscillator has applications in
describing several effects in other models too [17]. One being the resonance
effects in the antiferromagnetic bimeron in response to an alternating
current, which has applications in the detection of weak signals [15, 18, 19].

The Gilbert damping term is added phenomenologically to the Landau—-Lifshitz
model, to reproduce the damping of the precessing magnetic moment. Gilbert
damping is important in modeling other resonance features too, as its
temperature dependence affects them [20, 21], and in return in the
superconducting correlations that affect it [22]. The magnetization precession
in the ultra thin $Co_{20}Fe_{60}B_{20}$ layer stimulated by microwave voltage
under a large angle, needs modeling by Duffing oscillator too. This gets help
from the so called foldover features, again due to nonlinearity [16, 23, 24].

The consequences of the nonlinear nature of the coupled set of LLGJ system of
equations in the weak coupling regime was demonstrated recently in Ref. [14].
We showed in this regime, where the Josephson energy is small compared to the
magnetic energy, the $\varphi_{0}$ Josephson junction is equivalently
described by a scalar nonlinear Duffing equation. An anomalous dependence of
the ferromagnetic resonant frequency (FMR) with the increase of the Gilbert
damping was found. We showed that the damped precession of the magnetic moment
is dynamically driven by the Josephson supercurrent, and the resonance
behavior is given by the Duffing spring. The obtained results were based on
the numerical simulations. The role of dc superconducting current, and the
state with negative differential resistance (NDR) in IV-characteristic were
not clarified. Also, the effects of the Josephson to magnetic energy ratio and
the spin-orbit coupling (SOC) were not investigated at that time.

In the present paper, we study the nonlinear aspects of the magnetic dynamics
and IV-characteristics of the $\varphi_{0}$ Josephson junction in the
ferromagnetic resonance region. We compare description of the anomalous
damping dependence (ADD) exhibited by full LLGJ system of equations with
approximated equations and demonstrate the Duffing oscillator features in the
small parameter regime. Effects of the Josephson to magnetic energy ratio, and
the spin-orbit coupling on the ADD, referred to earlier as the $\alpha$-effect
[14] are demonstrated. By deriving the formula which couples the dc
superconducting current and maximal amplitude of magnetization we discuss the
correlation of superconducting current and the negative differential
resistance in the resonance region. Finally, we discuss the experimentally
important features by emphasizing the details of the magnetization dynamics
and the IV-characteristics of the $\varphi_{0}$ junction.

We have shown that in the limit of small system parameters; that is, the
Josephson to magnetic energy ratio $G$, the damping $\alpha$, and the spin-
orbit coupling $r$, the dynamics is given by the Duffing spring [14]. We focus
on the shift in resonance and the effects of nonlinear interactions. We give
semi-analytic models to explain our results in various limits.

The paper is organized as follows. In Section II we outline the theoretical
model and discuss the methods of calculations. The ferromagnetic resonance and
effects of system parameters on the anomalous damping dependence are
considered in Subsection A of Section III. In Subsection B we present
analytical description of the dynamics and IV-characteristics of the
$\varphi_{0}$ junction at small system parameters. Manifestation of the
negative differential resistance in IV-characteristics through the foldover
effect is discussed. We compare the description of the anomalous damping
dependence by full LLGJ system of equation with approximated equation, and
show how the Duffing oscillator captures the nonlinearities in the small
parameter regime in Subsection C. We present results on the critical damping
and derive the formula which couples the dc superconducting current and
maximal amplitude of magnetization in the ferromagnetic layer. Finally, in
Section IV we concludes the paper.

##  II II. Models and Method

The following section is closely related to our work in [13]. The
$\varphi_{0}$ junction [25, 12, 6] that we study is shown in Fig.1. The
current-phase relation in $varphi_{0}$ junction has the form
$I_{s}=I_{c}\sin{(\varphi-\varphi_{0})}$, where $\varphi_{0}=rM_{y}/M_{0}$,
$M_{y}$ denotes the component of magnetic moment in $\hat{y}$ direction,
$M_{0}$ is the modulus of the magnetization. The physics of $\varphi_{0}$
Josephson juncton is determined by system of equations which consists of
Landau-Lifshits-Gilbert (LLG), resistively capacitively shunted junction
(RCSJ) model expression with current-phase relation ($I_{s}$) described above,
and Josephson relation between phase difference and voltage.

The dynamics of the magnetic moment $\mathbf{M}$ is described by the LLG
equation [26]

Figure 1: Schematic view of SFS $\varphi_{0}$ Josephson junction. The external
current applied along x direction, ferromagnetic easy axis is along z
direction.

$\frac{d\mathbf{M}}{dt}=-\gamma\mathbf{M}\times\mathbf{H}_{eff}+\frac{\alpha}{M_{0}}\left(\mathbf{M}\times\frac{d\mathbf{M}}{dt}\right),$
(1)

where $\mathbf{M}$ is the magnetization vector, $\gamma$ is the gyromagnetic
relation, $\mathbf{H}_{eff}$ is the effective magnetic field, $\alpha$ is
Gilbert damping parameter, $M_{0}=|\mathbf{M}|$.

In order to find the expression for the effective magnetic field we have used
the model developed in Ref.[6], where it is assumed that the gradient of the
spin-orbit potential is along the easy axis of magnetization taken to be along
${\hat{z}}$. In this case the total energy of the system can be written as

$E_{\text{tot}}=-\frac{\Phi_{0}}{2\pi}\varphi
I+E_{s}\left(\varphi,\varphi_{0}\right)+E_{M}\left(\varphi_{0}\right),$ (2)

where $\varphi$ is the phase difference between the superconductors across the
junction, $I$ is the external current,
$E_{s}\left(\varphi,\varphi_{0}\right)=E_{J}\left[1-\cos\left(\varphi-\varphi_{0}\right)\right]$,
and $\displaystyle E_{J}=\Phi_{0}I_{c}/2\pi$ is the Josephson energy. Here
$\Phi_{0}$ is the flux quantum, $I_{c}$ is the critical current,
$r=l\upsilon_{so}/\upsilon_{F}$ $l=4hL/\hbar\upsilon_{F}$, $L$ is the length
of $F$ layer, $h$ is the exchange field of the $F$ layer,
$E_{M}=-K\mathcal{V}M^{2}_{z}/(2M_{0}^{2})$, the parameter
$\upsilon_{so}/\upsilon_{F}$ characterizes a relative strength of spin-orbit
interaction, $K$ is the anisotropic constant, and $\mathcal{V}$ is the volume
of the ferromagnetic ($F$) layer.

The effective field for LLG equation is determined by

$\displaystyle{\bf H_{eff}}$ $\displaystyle=$
$\displaystyle-\frac{1}{\mathcal{V}}\frac{\partial E_{tot}}{\partial{\bf M}}$
(3) $\displaystyle=$
$\displaystyle\frac{\Omega_{F}}{\gamma}\bigg{[}Gr\sin\bigg{(}\varphi-r\frac{M_{y}}{M_{0}}\bigg{)}{\bf\widehat{y}}+\frac{M_{z}}{M_{0}}{\bf\widehat{z}}\bigg{]}$

where $\Omega_{F}=\gamma K/M_{0}$ is frequency of ferromagnetic resonance and
$G=E_{J}/(K\mathcal{V})$ determines the ratio of Josephson energy to magnetic
one.

In order to describe the full dynamics $\varphi_{0}$ junction the LLG
equations should be supplemented by the equation for phase difference
$\varphi$, i.e. equation of RCSJ model for bias current and Josephson relation
for voltage. According to the extended RCSJ model, which takes into account
derivative of $\varphi_{0}$ phase shift, the current flowing through the
system in underdamped case is determined by

$\displaystyle I=\frac{\hbar
C}{2e}\frac{d^{2}\varphi}{dt^{2}}+\frac{\hbar}{2eR}\left[\frac{d\varphi}{dt}-\frac{r}{M_{0}}\frac{dM_{y}}{dt}\right]$
(4) $\displaystyle+I_{c}\sin{\left(\varphi-\frac{r}{M_{0}}M_{y}\right)}.$

where $I$ is the bias current, $C$ and $R$ are the capacitance and resistance
of Josephson junction respectively. The Josephson relation for voltage is
given by :

$\frac{\hbar}{2e}\frac{d\varphi}{dt}=V.$ (5)

We note that in the framework of RCSJ–model the displacement current is
proportional to the first derivative of voltage (or second derivative of phase
difference). From the other hand, the magnetization dynamics plays role of the
external force and first order derivative of $\varphi_{0}$ is a source of
external current for JJ. This was demonstrated in Ref.[27, 25] where the
authors included the first derivative of $\varphi_{0}$ as the source of the
electromotive force. Voltage is determined by the phase difference, and does
not depend on $\varphi_{0}$. From this point of view, in the framework of RCSJ
model the external current source cannot modify the expression for
displacement current. That’s why we do not include the second derivative of
$varphi_{0}$ in our model.

Using (1), (3), (4) and (5) we can write the system of equations, in
normalised variables, which describes the dynamics of $\varphi_{0}$ junction

$\begin{array}[]{llll}\displaystyle\dot{m}_{x}=\frac{\omega_{F}}{1+\alpha^{2}}\\{-m_{y}m_{z}+Grm_{z}\sin(\varphi-
rm_{y})\\\ -\alpha[m_{x}m_{z}^{2}+Grm_{x}m_{y}\sin(\varphi-
rm_{y})]\\},\vspace{0.2 cm}\\\
\displaystyle\dot{m}_{y}=\frac{\omega_{F}}{1+\alpha^{2}}\\{m_{x}m_{z}\\\
-\alpha[m_{y}m_{z}^{2}-Gr(m_{z}^{2}+m_{x}^{2})\sin(\varphi-
rm_{y})]\\},\vspace{0.2 cm}\\\
\displaystyle\dot{m}_{z}=\frac{\omega_{F}}{1+\alpha^{2}}\\{-Grm_{x}\sin(\varphi-
rm_{y})\\\ -\alpha[Grm_{y}m_{z}\sin(\varphi-
rm_{y})-m_{z}(m_{x}^{2}+m_{y}^{2})]\\},\vspace{0.2 cm}\\\
\displaystyle\dot{V}=\frac{1}{\beta_{c}}[I-V+r\dot{m}_{y}-\sin(\varphi-
rm_{y})],\\\ \vspace{0.05 cm}\\\ \displaystyle\dot{\varphi}=V\end{array}$ (6)

where $m_{x,y,z}=M_{x,y,z}/M_{0}$ and satisfy the constraint
$\sum_{i=x,y,z}m_{i}^{2}(t)=1$, $\beta_{c}=2eI_{c}CR^{2}/\hbar$ is McCumber
parameter. In order to use the same time scale in the LLG and RCSJ equations
in this system of equations we have normalized time to the $\omega^{-1}_{c}$,
where $\omega_{c}=\frac{2eI_{c}R}{\hbar}$, and
$\omega_{F}=\Omega_{F}/\omega_{c}$ is the normalized frequency of
ferromagnetic resonance $\Omega_{F}=\gamma K/M_{0}$. Bias current is
normalized to the critical current $I_{c}$ and voltage $V$ – to the
$V_{c}=I_{c}R$. The system of equations (6), is solved numerically using the
fourth-order Runge-Kutta method(see Ref.[14]).

##  III III. Results and Discussion

###  III.1 A. Effect of system parameters on the anomalous damping dependence

ADD of the FMR frequency with increasing $\alpha$ was discussed in Ref. [14].
It was found that the resonance curves demonstrate features of Duffing
oscillator, reflecting the nonlinear nature of Landau-Lifshitz-Gilbert-
Josephson (LLGJ) system of equations. There is a critical damping value at
which anomalous dependence comes into play. This critical value depends on the
system parameters. Here we present the details of such transformation from
usual to anomalous dependence with variation in spin-orbit coupling and ratio
of Josephson to magnetic energy.

To investigate the effect of damping, we calculate the maximal amplitude of
magnetization component $m_{y}$ taken at each value of the bias current based
on the LLGJ system of equations (6). In Fig.2 we show the voltage dependence
of maximal amplitude $m_{y}^{max}$ in the ferromagnetic resonance region at
different damping parameter and small values of Josephson to magnetic energy
ratio G=0.05 and spin-orbit coupling $r=0.05$. We found that the ferromagnetic
resonance curves demonstrate the different forms. An increase in damping shows
a nonuniform change in the resonant frequency: it is approaching the
$\omega_{F}$ instead of moving away with increase in $\alpha$. We stress that
this happens at small $G$ and $r$. We consider that such behavior can be
explained by the nonlinear nature of the LLGJ system of equations. There is a
manifestation of subharmonics of the FMR in Fig.2 at $\omega=1/2,1/3,1/4$.

Figure 2: Maximal amplitude of magnetization $m_{y}-$component at each values
of bias current and voltage along IV-characteristics of the $\varphi_{0}$
junction in the ferromagnetic resonance region for various $\alpha$. Inset
enlarges the main maximum. Parameters:
$\beta_{c}=25$,G=0.05,r=0.05,$\omega_{F}=0.5$.

We usually expect the resonance peak to move away from resonance as the
$\alpha$ increases. Figure 2 shows that this normal effect is accompanied with
an anomalous behaviour as can be seen in the inset to this figure, where the
resonance peak approaches $\omega_{F}$ as $\alpha$ increases [14].

The manifestation of FMR in IV-characteristics of the $\varphi_{0}$ junction
at three values of damping parameter is demonstrated in Fig. 3. The strong
deviation of the IV-curve is observing at $\alpha=0.01$, which is
characteristic value for many magnetic materials. This fact indicates that ADD
can be observed experimentally by measuring IV-characteristics in wide
interval of the damping parameter.

Figure 3: Part of the IV characteristic of the $\varphi_{0}$ junction at
$G=0.05,r=0.05$ and different values of Gilbert damping. The numbers show
$\alpha$ value. Inset shows the total IV-characteristic and arrow indicates
the resonance region

Interesting features of ADD appear by a variation of spin-orbit coupling. As
it was demonstrated in Ref.[28], an increase in SOC leads to the essential
change in IV-characteristics and magnetization precession in the ferromagnetic
resonance region. The nonlinearity is going stronger and the state with
negative differential resistance appears at large SOC.

Figure 4(a) demonstrates results of numerical simulations of $m_{y}^{max}$
dependence on $\alpha$ at different values of SOC parameter $r$.

Figure 4: (a) Demonstration of ADD at different values of SOC parameter $r$ at
$G=0.05$. Numbers indicate: 1 - $r=0.05$; 2 - $r=0.1$; 3 - $r=0.5$; Arrows
show critical $\alpha$ value, corresponded to the reversal in the $\alpha$
dependence (b) Demonstration of ADD at different values of the Josephson to
magnetic energy ratio $G$ at $r=0.05$. Numbers indicate: 1 - $G=0.01$; 2 -
$G=0.1$; 3 - $G=1$.

It shows two specific features of ADD. First, with an increase in $r$, the
critical value of $V_{peak}$ is decreasing (the curve moves away from
$\omega_{F}$). The second important feature is an increasing of
$\alpha_{crit}$ which is demonstrated in this figure by arrows.

Another model parameter which affects the phenomenon discussed in the present
paper is the ratio $G$ of Josephson to magnetic energies. Figure 4(b)
demonstrates the results of numerical simulations of $m_{y}^{max}$ dependence
on $\alpha$ at different values of G.

Similar to the effect of $r$, increasing $G$ also causes the value of
$\alpha_{crit}$ to increase. By changing the volume of the ferromagnetic
layer, the ferromagnetic energy and consequently the value of G can be changed
[6]. For small G, i.e. a situation where the magnetic energy is much larger
than the Josephson energy, the magnetic layer receives less energy, and its
amplitude decreases in the y direction, and also the maximum value of the
oscillation frequency is closer to the magnetic frequency, $\omega_{F}$.

###  III.2 B. Dynamics and IV-characteristics of the $\varphi_{0}$ junction at
small system parameters

As it was discussed in Refs.[6, 29, 30], in case of $G,r,\alpha<<1$ and
$m_{z}\approx 1$, first three equations of the system (6) can be simplified.
Taking into account $\varphi=\omega_{J}t$ and neglecting quadratic terms of
$m_{x}$ and $m_{y}$, we get

$\left\\{\begin{array}[]{ll}\displaystyle\dot{m}_{x}=\omega_{F}[-m_{y}+Gr\sin(\omega_{J}t)-\alpha
m_{x}]\vspace{0.2 cm}\\\ \displaystyle\dot{m}_{y}=\omega_{F}[m_{x}-\alpha
m_{y}],\end{array}\right.$ (7)

This system of equations can be written as the second order differential
equation with respect to $m_{y}$

$\displaystyle\ddot{m}_{y}+2\alpha\omega_{F}\dot{m}_{y}+\omega_{F}^{2}m_{y}=\omega_{F}^{2}Gr\sin\omega_{J}t.$
(8)

Corresponding solution for $m_{y}$ has the form

$m_{y}(t)=\frac{\omega_{+}-\omega_{-}}{r}\sin(\omega_{J}t)-\frac{\gamma_{+}+\gamma_{-}}{r}\cos(\omega_{J}t),$
(9)

where

$\omega_{\pm}=\frac{Gr^{2}\omega_{F}}{2}\frac{\omega_{J}\pm\omega_{F}}{\Omega_{\pm}},$
(10)

and

$\gamma_{\pm}=\frac{Gr^{2}\omega_{F}}{2}\frac{\alpha\omega_{J}}{\Omega_{\pm}}.$
(11)

with $\Omega_{\pm}=(\omega_{J}\pm\omega_{F})^{2}+(\alpha\omega_{J})^{2}$ (see
Ref.[6] and corresponded Erratum[31]).

When the Josephson frequency $\omega_{J}$ is approaching the ferromagnetic one
$\omega_{F}$, $m_{y}$ demonstrates the damped ferromagnetic resonance.
Differential resistance in the resonance region is decreasing and it is
manifested in the IV–characteristic as a resonance branch [7].

Taking into account $rm_{y}<<1$, we rewrite expression for superconducting
current as

$\displaystyle\displaystyle I_{s}(t)$ $\displaystyle=$
$\displaystyle\sin(\omega_{J}t-rm_{y}(t))$ (12) $\displaystyle=$
$\displaystyle\sin(\omega_{J}t)-rm_{y}\cos(\omega_{J}t)$

Using solution (9) we can obtain

$\displaystyle\displaystyle I_{s}(t)$ $\displaystyle=$
$\displaystyle\sin\omega_{J}t-\frac{\omega_{+}-\omega_{-}}{2}\sin
2\omega_{J}t$ (13) $\displaystyle+$
$\displaystyle\frac{\gamma_{+}+\gamma_{-}}{2}\cos 2\omega_{J}t+I_{0}(\alpha)$

where

$\displaystyle I_{0}=\frac{\gamma_{+}+\gamma_{-}}{2}.$ (14)

This superconducting current explains the appearance of the resonance branch
in the IV–characteristic. The generated current $I_{0}$ can be expressed
through the amplitude of $m_{y}$ and SOI parameter $r$

$\displaystyle I_{0}=\frac{r}{2}m_{y}^{max}(\omega_{J}),$ (15)

with $m_{y}^{max}(\omega_{J})$ being the frequency response of $m_{y}$.

At small model parameters $\alpha<<Gr<<1$ of SFS $\varphi_{0}$ Josephson
junction, the states with a negative differential resistance appear in the IV-
characteristics in the FMR region. Due to the nonlinearity, the resonance peak
is asymmetric. Increasing of nonlinearity leads to the bistability (foldover
effect). A natural questions appears if the states with a negative
differential resistance are the origin of the foldover and ADD. In order to
clarify this question, we show in Fig.5 a part of IV–characteristics of
$\varphi_{0}$ junction together with IV-characteristic of SIS junction in the
ferromagnetic resonance region and numerically calculated superconducting
current through this junction. The total IV–characteristics are demonstrated
in the inset to this figure.

Figure 5: IV-characteristic of $\varphi_{0}$ and SIS junctions and calculated
average superconducting current through the $\varphi_{0}$ junction.

We see the correlation of the foldover effect in superconducting current
(blue) with the NDR part of IV-curve. The peak in the superconducting current
and minimum of IV-curve are at the same voltage value. So, both effects
reflect the nonlinear features of the ferromagnetic resonance in $\varphi_{0}$
junction. But in contrast to the foldover and ADD effects, which start
manifesting themselves at relatively small deviations from linear case, the
nonlinearity in case of the NDR plays a more essential role.

We note that in the resonance region for considered limit of model parameters
the $m_{y}$ amplitude are coupled to the value of superconducting current (see
Eq.(15)). We stress the importance of the performed analysis demonstrating the
analytical coupling of time independent superconducting current and
magnetization, reflecting the Duffing oscillator features of the $\varphi_{0}$
junction.

As it is well known, the states with negative differential resistance appear
in IV-characteristics of Josephson structures in different physical
situations. In particular, the nonlinear superconducting structures being
driven far from equilibrium exhibit NDR states [32]. The NDR states plays an
essential role in the applications related to the THz radiation emission [33].
A detailed explanation for different types of negative differential resistance
(NDR) in Josephson junction (i.e., N-shaped and S-shaped) is introduced
in[34]. The authors emphasize that the nonlinear behavior of the Josephson
junction plays a key role in the NDR feature. In our case, the NDR states
appear as a result of system’s nonlinearity at small values of $\varphi_{0}$
junction parameters, such as SOC, ratio of Josephson to magnetic energy and
Gilbert damping. We demonstrate these effects here by by presented results of
detailed investigation of the NDR state at different system parameters and
discuss possibility of their control near the ferromagnetic resonance.

Figure 6 shows the effect of the spin-orbit coupling on the IV-characteristic
at $G=0.05$ and $\alpha=0.01$. We see the NDR feature which is getting more
pronounced with increase in $r$. A further increase in $r$ leads to the jump
down in voltage and then practically linear growth of IV-characteristic.

Figure 6: Enlarged parts of the IV-curves in the resonance region at different
values of SOC parameter $r$ at $\alpha=0.01$ and $G=0.05$. The numbers
indicate the increasing of $r$ from 0.09 to 0.15 with an increment 0.01. Inset
demonstrates total IV-characteristic ar $r=0.15$

An interesting question is concerning the effect of Gilbert damping. Results
of IV-characteristics simulations in the resonance region at a certain range
of the damping parameter $\alpha$ at $G=0.05$ and $r=0.13$ are shown in
Fig.7(a). In this case, the most pronounced characteristic appears at
$\alpha=0.01$. At chosen set of parameters $G=0.05$ and $r=0.13$, the range of
$\alpha$ with pronounced NDR features is $0.01\leqslant\alpha<0.014$.

The maximal amplitude $m^{max}_{y}$ as a function of voltage is shown in the
Fig.7(b). Based on the results, presented in Fig.7(a) and Fig.7(b), we came to
the important conclusion that the foldover effect (bistability) and NDR state
have strong correlations and have the same origin related to the nonlinearity
at small system parameters.

Figure 7: (a) Enlarged parts of the IV-curves at different values of $\alpha$;
(b) Voltage dependence of $m^{max}_{y}$ at different $\alpha$; (c)
$\alpha$-dependence of the resonance curve maximum in current ($I_{peak}$) and
voltage ($V_{peak}$). The numbers indicate the value of $\alpha$ from 0.01 to
0.018 in (a) and from 0.01 to 0.02 in (b) by an increment 0.001. Results are
obtained at $r=0.13$, and $G=0.05$.

But anomalous damping dependence does not show a one-to-one correlation with
either negative differential resistance or foldover effect. The resonance peak
positions of $m_{y}^{max}$ in bias current $I_{peak}$, and in voltage
$V_{peak}$ as the functions of $\alpha$ are demonstrated in Fig.7(c).
According to our results, we can divide $\alpha$ interval into two regions
(see Fig.7(c)). Region I introduces the values of $\alpha$ where the NDR
feature is present, while in the Region II it disappears. In the Region II the
foldover effect (bistability) disappears as well, but ADD is realized.

###  III.3 C. Duffing oscillator features of $\varphi_{0}$ junction and
critical damping

The system of equations (6) is nonlinear and very complex, so in order to
provide analytical study of dynamics of the $\varphi_{0}$ junction we need to
derive an approximated equation for some limited values of model parameters.
In Ref. [14] was shown that the resonance curves demonstrate features of
Duffing oscillator, reflecting the nonlinear nature of Landau-Lifshitz-
Gilbert-Josephson (LLGJ) system of equations. In this section, we present
analytical approach to describe the nonlinear dynamics of
$\varphi_{0}$-junction and compare analytical results followed from
approximated Duffing equation with numerical simulations of total system of
equations (6). We show that in the limit of $\alpha<<G,r<<1$, we arrive to the
Duffing oscillator. We start by first three equations of the (6) for
magnetization components

$\displaystyle\frac{\dot{m}_{x}}{\omega_{F}}$ $\displaystyle=$ $\displaystyle-
m_{y}m_{z}+Grm_{z}\sin(\varphi-rm_{y})-\alpha m_{x}m_{z}^{2}$
$\displaystyle\frac{\dot{m}_{y}}{\omega_{F}}$ $\displaystyle=$ $\displaystyle
m_{x}m_{z}-\alpha m_{y}m_{z}^{2}$ (16)
$\displaystyle\frac{\dot{m}_{z}}{\omega_{F}}$ $\displaystyle=$ $\displaystyle-
Grm_{x}\sin(\varphi-rm_{y})+\alpha m_{z}(m_{x}^{2}+m_{y}^{2})$

Simplifying this system of equations by the same procedure as it was done in
Ref.[14] we can write equation for $m_{y}$ as

$\displaystyle{\ddot{m}_{y}}+2\xi\alpha\dot{m}_{y}+\xi^{2}(1+\alpha^{2})m_{y}-\xi^{2}(1+\alpha^{2}-\alpha^{4})m_{y}^{3}=\xi^{2}\gamma\sin\varphi.$
(17)

Finally, by neglecting the $\alpha^{2}$ and $\alpha^{4}$ terms which are much
smaller than 1, we come to well known Duffing equation,

$\displaystyle{\ddot{m}_{y}}+2\omega_{F}\alpha\dot{m}_{y}+\omega_{F}^{2}m_{y}-\omega_{F}^{2}m_{y}^{3}=\omega_{F}^{2}Gr\sin\varphi.$
(18)

In the small parameter range, this Duffing equation can describe the dynamics
of $m_{y}$. We will have the full dynamics, once we consider the coupling with
the Josephson equation,

$\displaystyle{\ddot{\varphi}}+{1\over{\beta_{c}}}[{\dot{\varphi}}-r{\dot{m}}_{y}+\sin(\varphi-
rm_{y})]={1\over{\beta_{c}}}I.$ (19)

The system of equations, (18) and (19) can replace the LLGJ equations in the
limit of $G,r<<1$ and $G,r>>\alpha$.

Taking into account $\varphi=\omega_{J}t$ we can right analytically obtained
frequency response for equation (18)

$\displaystyle(m_{y}^{max})^{2}=\frac{\big{(}Gr\big{)}^{2}}{\big{[}\omega^{2}-1+\frac{3}{4}(m_{y}^{max})^{2}\big{]}^{2}+\big{(}2\alpha\omega\big{)}^{2}}$

where $\omega=\omega_{J}/\omega_{F}$. From Eq. (III.3) we get

$\displaystyle(m_{y}^{max})^{6}$ $\displaystyle+$
$\displaystyle\frac{8}{3}(\omega^{2}-1)(m_{y}^{max})^{4}$ (21)
$\displaystyle+$
$\displaystyle\bigg{(}\frac{4}{3}\bigg{)}^{2}\bigg{[}(\omega^{2}-1)^{2}+\big{(}2\alpha\omega\big{)}^{2}\bigg{]}(m_{y}^{max})^{2}$
$\displaystyle-$ $\displaystyle\bigg{(}\frac{4}{3}Gr\bigg{)}^{2}=0.$

This equation allows to determine analytically frequency dependence of the
$m^{max}_{y}$ amplitude. To find it we solve the equation (21) by the Newton
method. Results of analytical calculations (blue dots) corresponded to (21)
and numerical one (red doted line) corresponded to the full system of equation
(6) are demonstrated in Fig.8.

Figure 8: Numerically (curve 1) and analytically (curve 2) calculated
amplitude dependence of $m_{y}$.

We can see that they are close to each other which proves the correctness of
the chosen approximation. Both curves demonstrate an asymmetric resonance
peak, which is common for Duffing oscillator. When a role of the cubic term is
getting larger, we observe a bistability of the resonance curve, which is
usually called a foldover effect. Note that the foldover effect can be also
achieved by the damping decreasing; i.e., by the decreasing of dissipative
term in (18), we can increase the influence of the cubic term in this
equation.

The comparison of analytically and numerically calculated superconducting
current as a function of the Josephson frequency is demonstrated in Fig. 9. We
note that in our normalization $V=\omega_{J}$.

Figure 9: Numerically calculated superconducting current for SFS junction
(plot 1) and analytical $I_{0}$ (plot 2) and superconducting current for SIS
junction (plot 3).

We can see the manifestation of the asymmetric resonance peak in the frequency
dependence of superconducting current. So, the approximated system of
equations 7 reflects one of the main feature of Duffing oscillator.

Figure (10) compares anomalous damping dependence of the resonance peak of
$m_{y}^{max}(V)$ calculated numerically according to the full LLGJ system of
equations (6) with calculated numerically according to the generalized Duffing
model (equations (17, 19)). We see that in the damping parameter interval
[0.001 – 0.2] the coincidence of the dependences is enough good.

Figure 10: The $\alpha$-dependence of the resonance maximum of
$m_{y}^{max}(V)$ in the damping parameter interval [0.001 – 0.12]. Green
squares show results calculated numerically according to the full system of
equations (6), blue circles show results calculated numerically according to
the generalized Duffing and Josephson equations (17,19). The dashed line
connects the symbols to guide eyes. Solid line show analytical
$\alpha$-dependence calculated according to the Eq. (22). All calculation have
been done at $\beta_{c}=25$, G=0.05, r=0.05, $\omega_{F}=0.5$.

Using equation (18) with $\varphi=\omega_{J}t$, we can find (see Supplementary
materials LABEL:si:1) a relation between position of the resonance peak in
$m_{y}^{max}(V)$ dependence and damping

$\displaystyle\omega_{peak}=\sqrt{\frac{1-3\alpha^{2}}{2}+\frac{1}{2}\sqrt{(1-\alpha^{2})^{2}-12(\frac{Gr}{4\alpha})^{2}}}$
(22)

where $\omega_{peak}=\frac{\omega_{J,peak}}{\omega_{F}}$ determines the
position of the resonance peak.

Equation (22) allows to find the formula for critical damping $\alpha_{crit}$
which is an important parameter determining the reversal point in damping
dependence of the resonance peak of $m_{y}^{max}(V)$ .

Taking into account equation (22) we can write equation with respect of
$Gr/(4\alpha)$ (See supplementary materials LABEL:si:1).

$\displaystyle 9\bigg{(}\frac{Gr}{4\alpha_{crit}}\bigg{)}^{4}$
$\displaystyle+$ $\displaystyle
3\alpha_{crit}^{2}(10\alpha_{crit}^{2}-1)\bigg{(}\frac{Gr}{4\alpha_{crit}}\bigg{)}^{2}$
$\displaystyle-$ $\displaystyle 2\alpha_{crit}^{4}(\alpha_{crit}^{2}-1)^{2}=0$

Using approximation $10\alpha_{crit}^{2}<<1$ and $\alpha_{crit}^{2}<<1$ it
gives (see Supplementary Materials)

$\displaystyle\alpha_{crit}\approx\frac{1}{2}\sqrt{\sqrt{\frac{3}{2}}Gr}$ (24)

Table 1: A comparison between the numerical and analytical values of
$\alpha_{crit.}$ at different values of $G$ and $r$.

G | r | Gr | $\alpha_{crit.},numerics$ | $\alpha_{crit.},analytics$
---|---|---|---|---
0.01 | 0.05 | 0.0005 | 0.0100 | 0.0123
0.05 | 0.05 | 0.0025 | 0.0300 | 0.0276
0.05 | 0.10 | 0.0050 | 0.0400 | 0.0391
0.05 | 0.30 | 0.0150 | 0.0700 | 0.0677
0.05 | 0.50 | 0.0250 | 0.0900 | 0.0874
0.10 | 0.05 | 0.0050 | 0.0391 | 0.0391
0.60 | 0.05 | 0.0300 | 0.0950 | 0.0958
0.70 | 0.05 | 0.0350 | 0.1000 | 0.1035
1.00 | 0.05 | 0.0500 | 0.1200 | 0.1237

Figure 11 presents comparison of numerical and analytical results
$\alpha_{crit}$ versus $Gr$.

Figure 11: Numerical calculations according to Eq. (6) (squares), analytical
according to Eq. (III.3)(solid line) and approximated analytical according to
Eq. (24) (dashed line).

As we see, it shows a good agreement of numerical and analytical results of
calculations at small product of Josephson to magnetic energy ratio and spin-
orbit interaction.

##  IV IV. CONCLUSIONS

The understanding of the nonlinear features of magnetization dynamics in
superconductor-ferromagnet-superconductor Josephson junction and their
manifestation in the IV-characteristics has implications for superconductor
spintronics, and modern information technology. In $\varphi_{0}$ junctions the
nonlinear features can affect the control of magnetization precession by
superconducting current and external electromagnetic radiation [28].

Here, using numerical and analytic approaches, we have demonstrated that at
small system parameters, namely, the damping, spin-orbit interaction and
Josephson to magnetic energy ratio in $\varphi_{0}$ junction, magnetic
dynamics is reduced to the dynamics of the scalar Duffing oscillator, driven
by the Josephson oscillations. We have clarified the role of increasing
superconducting current in the resonance region leading to the foldover effect
in the ferromagnet magnetization. We have demonstrated the parameter
dependence of the anomalous ferromagnetic resonant shifting with anomalous
damping dependence due to nonlinearity of the full LLGJ equation and in its
different approximations. We have derived the analytical expression for
critical damping value. Also, we demonstrated appearance of negative
differential resistance in the IV-characteristics and the correlation with
occurrence of the foldover effect in the magnetization of ferromagnet.

We have stressed that the manifestation of negative differential resistance is
related to the nonlinear features of the system[35, 34]. It was demonstrated
that in the small model parameters case the equation for magnetic subsystem
takes form of Duffing equation where nonlinearity manifest itself as the cubic
term. We have shown that the appearance of negative differential resistance in
the I-V curve is related to the appearance of foldover in the
$m_{y}^{max}$-$V$ curve.

We believe that the experimentally measured IV-characteristics of
$\varphi_{0}$ junction with manifestations discussed in detail in the present
paper, would allow close investigations of its nonlinear features important
for superconductor electronics and spintronics.

##  V Supplementary

In supplementary material are presented the details of calculations for Eq.22
and Eq.24.

##  VI Funding

Numerical simulations were funded by Project No. 18-71-10095 of the Russian
Science Foundation. The presented results concerning the calculations of DC
superconducting current in the section V are supported by the Russian Science
Foundation in the framework of project 22-42-04408. A.J. and M.R.K. are
grateful to IASBS for financial support.

## References

  * Buzdin [2008] Buzdin, A. _Physical Review Letters_ 2008, _101_ (10), 107005\.
  * Linder and Robinson [2015] Linder, J., Robinson, J. W. _Nature Physics_ 2015, _11_ (4), 307-315.
  * Bobkova et al. [2018] Bobkova, I., Bobkov, A., Silaev, M. _Physical Review B_ 2018, _98_ (1), 014521.
  * Bobkova et al. [2020] Bobkova, I., Bobkov, A., Rahmonov, I., Mazanik, A.; Sengupta, K., Shukrinov, Y. M. _Physical Review B_ 2020, _102_ (13), 134505\.
  * Szombati et al. [2016] Szombati, D., Nadj-Perge, S., Car, D., Plissard, S., Bakkers, E., Kouwenhoven, L. _Nature Physics_ 2016, _12_ (6), 568–572.
  * Konschelle and Buzdin [2009] Konschelle, F., Buzdin, A. _Physical Review Letters_ 2009, _102_ (1), 017001.
  * Shukrinov et al. [2019] Shukrinov, Y. M., Rahmonov, I., Sengupta, K. _Physical Review B_ 2019, _99_ (22), 224513.
  * Shukrinov [2022] Shukrinov, Y. M. _Physics-Uspekhi_ 2022, _65_ (4), 317.
  * Linder and Yokoyama [2011] Linder, J., Yokoyama, T. _Physical Review B_ 2011, _83_ (1), 012501\.
  * Hoffman et al. [2012] Hoffman, S., Blanter, Y. M., Tserkovnyak, Y. _Physical Review B_ 2012, _86_ (5), 054427.
  * Eschrig [2015] Eschrig, M. _Reports on Progress in Physics_ 2015, _78_ (10), 104501.
  * Shukrinov et al. [2017] Shukrinov, Y. M., Rahmonov, I., Sengupta, K., Buzdin, A. _Applied Physics Letters_ 2017, _110_ (18), 182407.
  * Mazanik et al. [2020] Mazanik, A., Rahmonov, I., Botha, A., Shukrinov, Y. M. _Physical Review Applied_ 2020, _14_ (1), 014003.
  * Shukrinov et al. [2021] Shukrinov, Y. M., Rahmonov, I., Janalizadeh, A., Kolahchi, M. _Physical Review B_ 2021, _104_ (22), 224511.
  * Shen et al. [2020] Shen, L., Xia, J., Zhang, X., Ezawa, M., Tretiakov, O. A., Liu, X., Zhao, G., Zhou, Y. _Physical review letters_ 2020, _124_ (3), 037202\.
  * Azovtsev et al. [2021] Azovtsev, A. V., Nikitchenko, A. I., Pertsev, N. A. _Physical Review Materials_ 2021, _5_ (5), 054601.
  * Zhu et al. [2020] Zhu, J., Zhang, T., Yang, Y., Huang, R. _Applied Physics Reviews_ 2020, _7_ (1), 011312.
  * Wang et al. [1999] Wang, G., Chen, D., Lin, J., Chen, X. _IEEE Transactions on industrial electronics_ 1999, _46_ (2), 440–444.
  * Almog et al. [2007] Almog, R., Zaitsev, S., Shtempluck, O., Buks, E. _Physical Review Letters_ 2007, _98_ (7), 078103.
  * Zhao et al. [2016] Zhao, Y., Song, Q., Yang, S.-H., Su, T., Yuan, W., Parkin, S. S., Shi, J., Han, W. _Scientific reports_ 2016, _6_ (1), 1–8.
  * Yao et al. [2018] Yao, Y., Song, Q., Takamura, Y., Cascales, J. P., Yuan, W., Ma, Y., Yun, Y., Xie, X., Moodera, J. S., Han, W. _Physical Review B_ 2018, _97_ (22), 224414.
  * Silaev [2020] Silaev, M. _Physical Review B_ 2020, _102_ (14), 144521.
  * Nayfeh and Mook [1979] Nayfeh, A., Mook, D. _Nonlinear Oscillations. John Willey & Sons_,  1979\.
  * Chen et al. [2009] Chen, W., De Loubens, G., Beaujour, J.-M., Sun, J., Kent, A. _Applied Physics Letters_ 2009, _95_ (17), 172513.
  * Guarcello and Bergeret [2020] Guarcello, C., Bergeret, F. _Physical Review Applied_ 2020, _13_ (3), 034012.
  * Landau and Lifshitz [1955] Landau, L., Lifshitz, E. _Phys. Rev_ 1955, _100_ , 1243.
  * Rabinovich et al. [2019] Rabinovich, D. S., Bobkova, I. V., Bobkov, A. M., Silaev, M. A. _Phys. Rev. Lett._ 2019, _123_ , 207001. doi:10.1103/PhysRevLett.123.207001.
  * Abdelmoneim et al. [2022] Abdelmoneim, S., Shukrinov, Y. M., Kulikov, K., ElSamman, H., Nashaat, M. _Phys. Rev. B_ 2022, _106_ , 014505.
  * Shukrinov and Rahmonov [2020] Shukrinov, Y. M., Rahmonov, I. _Physics of Particles and Nuclei_ 2020, _51_ (4), 816–822.
  * Shukrinov et al. [2020] Shukrinov, Y. M., Rahmonov, I., Botha, A. _Low Temperature Physics_ 2020, _46_ (9), 932–938.
  * Konschelle and Buzdin [2019] Konschelle, F., Buzdin, A. _Physical Review Letters_ 2019, _123_ (16), 169901.
  * Pedersen et al. [2014] Pedersen, N. F., Filatrella, G., Pierro, V., Sørensen, M. P. _Physica C: Superconductivity and its Applications_ 2014, _503_ , 178–182.
  * Kadowaki et al. [2008] Kadowaki, K., Yamaguchi, H., Kawamata, K., Yamamoto, T., Minami, H., Kakeya, I., Welp, U., Ozyuzer, L., Koshelev, A., Kurter, C., Gray, K., Kwok, W.-K. _Physica C: Superconductivity and its applications_ 2008, _468_ (7-10), 634–639.
  * Filatrella et al. [2014] Filatrella, G., Pierro, V., Pedersen, N. F., Sorensen, M. P. _IEEE Transactions on Applied Superconductivity_ 2014, _24_ (6), 1–7.
  * Nagel et al. [2008] Nagel, J., Speer, D., Gaber, T., Sterck, A., Eichhorn, R., Reimann, P., Ilin, K., Siegel, M., Koelle, D., Kleiner, R. _Physical Review Letters_ 2008, _100_ , 217001.


11institutetext:  T-5, Center for Nonlinear Studies
Los Alamos National Laboratory
Los Alamos, New Mexico 87545
<EMAIL_ADDRESS>

# A Reflection on the Structure and Process of the Web of Data

Marko A. Rodriguez

11footnotetext: Rodriguez, M.A., “A Reflection on the Structure and Process of
the Web of Data,” Bulletin of the American Society for Information Science and
Technology, volume 35, number 6, pages 38–43, ISSN:1550-8366, LA-UR-09-03724,
September 2009.

The Web community has introduced a set of standards and technologies for
representing, querying, and manipulating a globally distributed data structure
known as the Web of Data. The proponents of the Web of Data envision much of
the world’s data being interrelated and openly accessible to the general
public. This vision is analogous in many ways to the Web of Documents of
common knowledge, but instead of making documents and media openly accessible,
the focus is on making data openly accessible. In providing data for public
use, there has been a stimulated interest in a movement dubbed Open Data
opendata:miller2008 . Open Data is analogous in many ways to the Open Source
movement. However, instead of focusing on software, Open Data is focused on
the legal and licensing issues around publicly exposed data. Together, various
technological and legal tools are laying the groundwork for the future of
global-scale data management on the Web. As of today, in its early form, the
Web of Data hosts a variety of data sets that include encyclopedic facts, drug
and protein data, metadata on music, books and scholarly articles, social
network representations, geospatial information, and many other types of
information. The size and diversity of the Web of Data is a demonstration of
the flexibility of the underlying standards and the overall feasibility of the
project as a whole. The purpose of this article is to provide a review of the
technological underpinnings of the Web of Data as well as some of the hurdles
that need to be overcome if the Web of Data is to emerge as the defacto medium
for data representation, distribution, and ultimately, processing.

Technically, on the Web of Data, Uniform Resource Identifiers (URI) are used
to identify resources uri:berners2005 . For example, depending on what is
being modeled, a URI can denote a city, a protein, a music album, a scholarly
article, a person, etc. In general, any thing can be assigned a URI. An
example URI is http://www.lanl.gov#marko. This URI denotes the author of this
article, Marko. The URI has information pertaining to the what (marko), where
(www.lanl.gov), and how (http) of a resource. The URI is more general than the
URL of common knowledge as URIs are not required to resolve to retrievable
digital objects such as documents and media. Instead, URIs can denote abstract
concepts such as the person Marko, the class of dogs, or the notion of
friendship. Finally, the space of all URIs is an inherently distributed and
theoretically infinite space. This makes the URI space fit to represent
massive amounts of data distributed world wide. A convenient consequence of
this space is that the Web of Data can emerge atop it. However, while URIs can
denote things, they can not denote how things relates to each other. Relating
URIs is necessary in order to give greater meaning and context to datum.
Moreover, relating URIs is necessary to create the “web” aspect of the Web of
Data.

The Resource Description Framework (RDF) is a standardized data model for
linking URIs in order to create a network/graph of space of all URIs
rdfintro:miller1998 . RDF also supports the linking of URIs to primitive
literals such as strings, integers, floating point values, etc. An example RDF
statement to denote the fact that “Marko knows Fluffy” is
$\langle$http://www.lanl.gov#marko, http://xmlns.com/foaf/0.1/knows,
http://www.lanl.gov#fluffy$\rangle$. In order to make long URIs more readable,
namespace prefixes are generally used. With namespace prefixes, the previous
statement can be represented as
$\langle\texttt{lanl:marko},\texttt{foaf:knows},\texttt{lanl:fluffy}\rangle$.
All RDF statements have this three component form, where there exists a
subject (lanl:marko), a predicate (foaf:knows), and an object (lanl:fluffy).
As such, an RDF statement is also known as a triple. A URI can be involved in
multiple statements. For example, it is possible to state that, while being
known by Marko, Fluffy is also a dog, 5 years old, and lives in Santa Fe, New
Mexico. Data on Santa Fe and data on Fluffy become merged when statements
involving their two URIs are joined (directly or indirectly through multiple
links). The Web of Data becomes powerful when seemingly different data sets
are interlinked. The fact that Fluffy lives in Santa Fe automatically connects
data about Fluffy to geographic and encyclopedic data about Santa Fe, New
Mexico—its geospatial coordinates, nearby cities, culture, population, etc. As
more and more statements are added to the Web of Data, the Web of Data serves,
in a sense, as a global database of interlinked heterogeneous data. The
combination of both the URI and RDF has moved the World Wide Web beyond a Web
of Documents to that of a Web of Data, where every minutia of information can
be represented and interlinked for consumption by both man and machine.

RDF’s original use case has evolved beyond that of a logic-language for
knowledge representation and reasoning on the Semantic Web
webinterpret:rodriguez2009 . As the foundational technology of the Web of
Data, RDF can be seen as a general-purpose data model. It can be used to model
formal knowledge (the Semantic Web), graph/networks (the Giant Global Graph),
and software and abstract virtual machines (the Web of Programs) to name a
few. In many ways, URIs and RDF afford a memory structure analogous to the
local memory of a physical machine except that this memory structure is
distributed over physical machines world wide. Each physical machine stores
and manages a subset representation of the full Web of Data. RDF can be stored
on a physical machine in many ways. A simple, straightforward way is to
represent RDF statements in a file—an RDF document. A common misconception is
that RDF and RDF/XML are one in the same. RDF is a data model that has various
serialized representations with RDF/XML being one such serialization. Other
popular serializations include N3 and N-Triple. Thus, there are many types of
RDF documents. For the small-scale exposure of RDF data, an RDF document
suffices. For the large-scale exposure and processing of RDF data, an RDF
repository known as a triple store or graph database is usually the chosen
solution. The expanded use of RDF has been greatly facilitated by the
continued increase in the capacity and speed of RDF triple stores. Modern
high-end RDF triple stores can hold and process on the order of 10 billion
triples. Example high-end triple stores include Neo4j111Neo4j is available at
http://neo4j.org/. and AllegroGraph222AllegroGraph is available at
http://www.franz.com/agraph/allegrograph/.. What has been the sole territory
of relational database technologies may soon be displaced by the use of RDF
and the triple store. Moreover, because RDF is the common data model utilized
by triple stores, it is possible to integrate data sets across different
triple stores—across different RDF data providers. This integration is
conveniently afforded by the URI and RDF as Web standards and is a function
foreign to the relational database domain. With the Web of Data, no longer is
information isolated in individual inaccessible data silos, but is instead
exposed in an open and interconnected environment—the Web environment. The
means to integrate RDF data across different RDF data sets is explained next.

##  1 Linked Data and a Distributed Data Structure

In an effort to provide a seamless integration between the data provided by
different RDF data providers, the Linked Data community is focused on
developing the specifications and tools for linking RDF data sets into a
single, global “web of data” linkeddata:bizer2008 . Two RDF data sets link
together when one data set uses a URI maintained by another. For example,
suppose the URI lanl:fluffy minted and maintained by the Los Alamos National
Laboratory (LANL). As previously explained, this URI is denoting something in
the world—namely Fluffy. However, it is possible for someone other than LANL
to express statements about Fluffy. Assume that the Rensselaer Polytechnic
Institute (RPI) mints their own URI rpi:fluffy to denote Fluffy, where rpi is
the namespace prefix that resolves to http://www.rpi.edu#. At this point, both
lanl:fluffy and rpi:fluffy denote the same thing—they both denote the same
real-world object known as Fluffy. This idea is diagrammed in Figure 1a, where
the dashed lines identify which worldly things the URIs stand in reference to.
In order to link the LANL data set with the RPI data set, LANL can add the RDF
statement
$\langle\texttt{lanl:fluffy},\texttt{owl:sameAs},\texttt{rpi:fluffy}\rangle$
to its data set. This statement states that both lanl:fluffy and rpi:fluffy
denote the same real-world thing, Fluffy. This idea is diagrammed in Figure
1b. Given this statement, it is possible to traverse from lanl:fluffy (LANL)
to rpi:fluffy (RPI) and thus, migrate from the LANL data set to the RPI data
set. When two data sets denote the same thing, they can be linked.

Figure 1: Two RDF repositories utilize different URIs to denote the same
thing. By making it explicit that two URIs denote the same thing (i.e.
$\langle\texttt{lanl:fluffy},\texttt{owl:sameAs},\texttt{rpi:fluffy}\rangle$),
it is possible to merge data sets together. This merging of data sets is what
creates the Web of Data.

The Linked Data community is interested in both unifying RDF data sets as well
as specifying the behaviors associated with URI resolution. A Linked Data-
compliant data provider should return data when a URI is dereferenced—when a
representation of the resource being identified by the URI is requested. More
specifically, when a URI is dereferenced, a collection of statements
associated with that URI should be returned in some RDF serialization such as
RDF/XML. Given the example above, if a machine dereferences lanl:fluffy, it
will get the statement $\langle$lanl:fluffy, owl:sameAs, rpi:fluffy$\rangle$
returned to it. In other words, LANL returns all the RDF statements for which
lanl:fluffy is the subject of the triple (i.e. the outgoing edges from
lanl:fluffy). Now, the machine knows that lanl:fluffy and rpi:fluffy denote
the same thing. Thus, if it wants to know what RPI has stated about Fluffy, it
will dereference rpi:fluffy. Upon doing so, it should get the statement
$\langle$rpi:fluffy, rdf:type, rpi:Dog$\rangle$ returned to it. To the
machine, the Web of Data is one expansive interlinked web of URIs. To the
underlying servers, the Web of Data is broken up into multiple RDF subgraphs
(multiple data providers) and linked together when one data provider
references a URI minted and maintained by another data provider. It is noted
that resolving a Linked Data-compliant store’s URI is one way of getting data
from the Web of Data. For more complicated data gathering situations, many RDF
data providers expose SPARQL end-points to their triple stores. SPARQL is a
query language similar to SQL, but focused on graph queries as opposed to
table queries sparql:prud2004 .

An interesting consequence of the Web of Data is that it can greatly shift the
role of application and data providers. Currently, web applications are
required to maintain their own data source. For example, Amazon.com maintains
its database of books, Springer its database of journal articles, and iTunes
its database of music metadata. In order for users to utilize this data in
interesting ways, these same data providers must provide a front-end
application to interact with the data. In this way, data providers and
application developers are one in the same entity. This idea is diagrammed in
Figure 2a, where each application utilizes its own back-end database to
provide its front-end application with data. With the Web of Data, this model
is significantly altered. On the Web of Data, application providers and data
providers are cleanly separated. Data providers can provide and interlink
book, article, and music data on the Web of Data and application providers can
develop software to utilize this data for different end-user services—book
recommendations, citation analysis, and music metadata population. Moreover,
this same data can be utilized by multiple different application developers
and thus, this can yield many ways for the end-user to interact with the Web
of Data. In other words, Amazon.com’s data may be more efficiently presented
and processed if it was open for any developer to create a front-end
application for it. This idea is diagrammed in Figure 2b. The clean separation
between data and application providers is already taking place as plenty of
interlinked heterogeneous data currently exists on the Web of Data. A few
examples are provided here. Book data can be found at Amazon.com’s RDF
BookMashup and the RDF representation of Project Guttenburg. Scholarly data is
provided by the Digital Bibliography and Library Project (DBLP), ACM, IEEE,
amongst many others. Finally, various music data sets exist such as
MusicBrainz and AudioScrobbler. This data is leveraged, as mentioned
previously, by resolving URIs. For example, if one were to dereference this
URI in a standard web browser
http://rdf.freebase.com/ns/guid.9202a8c04000641f800000000001a49d, what is
returned is a set of RDF statements (as an RDF/XML document) linking this URI
to other URIs and literals. Accessible, interlinked, structured data is the
point of the Web of Data. An ecology of applications leveraging this data may
greatly advance applications and algorithms for processing data as no longer
are application developers burdened by the “cold start” problem of requiring
large-amounts of data to initiate a successful service faith2:rodriguez2009 .
No longer will consumers be confined to use certain web applications as no
longer are applications and data so tightly coupled.

Figure 2: a.) The typical web application requires its own data source on
which to provide its service. b.) On the Web of Data, application providers
are cleanly separated from the data providers.

##  2 Linked Process and a Distributed Process Infrastructure

While the Web of Data and the efforts expended by the Linked Data community
have provided a path towards global-scale data management, this model is
lacking one important component: an infrastructure for data processing. A
significant hurdle to overcome for this community is that of distributed
processing on this distributed data structure. Traversing the Web of Data is
not quite the same as traversing the Web of Documents. For the human, it is
reasonable to traverse from URI to URI exploring the Web of Data in a manner
similar to how the Web of Documents is traversed. That is, a human, using
their Web browser, can resolve URIs and view the RDF data returned. Moreover,
various human-friendly RDF browsers exist (usually in the form of a Web
browser plugin) to make it easy for humans to view and traverse the data on
the Web of Data. However, for a machine (i.e. an application, an algorithm),
the Web of Data can be traversed much faster than what a human can do by
manually clicking from URI to URI. Moreover, there will be orders of magnitude
more resources and links on the Web of Data than what is found on the Web of
Documents. While a machine can crawl and pull the data to its local
environment for processing, this becomes inefficient when the data
requirements span large parts of the Web of Data. Again, note that every time
a URI is dereferenced, the resolving server prepares an RDF subgraph and
returns it (over the wire) to the requesting machine. Thus, “traversing” the
Web of Data requires data to be migrated to the traversing machine and
processed remotely from the data source. This architecture is analogous to the
current Web of Documents whereby “traversing” the Web of Documents pulls HTML
documents and media to the requesting machine. For human consumption, this is
necessary as data/documents must be rendered where the human is physically
located—remote from the data source. For a machine (a virtual machine) its
physical location need not be a factor in how data is consumed and processed.
Thus, for processing large parts of a distributed data structure, a more
efficient mechanism would be to migrate the process between data providers so
that information is not pulled over the wire, but instead, processed where the
data is maintained. In other words, an efficient mechanism for processing the
Web of Data would be to move the process to the data, not the data to the
process rodriguez:gpsemnet2009 .

For the Web of Documents, the search engine philosophy of “download and index”
has made it possible for end users to find information in a more efficient
manner than by simply surfing and bookmarking. With modern commercial triple
stores scaling to the order of 10 billion triples, centralized indexing
repositories will have to contend with not only the volume of data, but also
the computational complexities of analyzing such data in sophisticated ways.
The Web of Data provides a much richer machine processable data structure than
what is provided by HTML and thus, antiquated keyword search simply does not
take significant advantage of what the Web of Data is providing. The future of
the Web of Data will be rife with algorithms from many schools of
thought—formal logic, graph analysis, object-oriented programming, etc.
webinterpret:rodriguez2009 . Many of these algorithms will compute across
various underlying stores of the Web of Data and will require a distributed
Turing complete infrastructure to do so. For any algorithm of sufficient
complexity, there is simply too much data to pull over the wire and thus, the
Web of Data in its current form greatly reduces what is possible. This is an
unfortunate state of affairs. Given the potential role of the Web of Data as
the defacto medium for interconnecting data, a distributed computing
environment is necessary. The Linked Data community needs a parallel Linked
Process effort. In a sense, data providers already expose their processors for
public use by way of their SPARQL-endpoints. SPARQL serves as an on-site data
processing language. However, this language, being a query language, is not
sufficient for representing complex algorithms. What is needed is a framework
that is more general-purpose and respective of the three following basic
requirements:

  1. 1.

safe: applications must not be able to destroy the integrity of the open
processor or its data set when using this infrastructure.

  2. 2.

efficient: applications must run faster in this infrastructure than what is
possible when pulling the required data over the wire.

  3. 3.

easy to use: application developers must be be able to utilize common
programming languages and packages and be relatively blind to the underlying
infrastructure.

Developing a distributed process infrastructure that accounts for these three
factors will ultimately drive its adoption. With the widespread adoption of
such a processing infrastructure by RDF data providers, the Web of Data will
reach a new level of functionality. No longer will the Web of Data be only a
database serving data over the wire to third-party applications, but instead,
a distributed computing environment supporting complex algorithms that can
leverage rich data in ways not previously possible in the history of
computing. The unification of Linked Data and Linked Process in many ways is
similar to “cloud computing.” However, with the integration of data sets and
hardware processors world wide, this “cloud” will be much richer and more
decentralized than what exists with other cloud providers. In this form, the
Web of Data will afford the world a democratization of both data and process
and may perhaps enjoy a frenzied adoption similar to what has occurred with
its predecessor, the Web of Documents.

##  3 Conclusion

The Web of Data provides an infrastructure that supports an instantiation of a
distributed graph of Web resources. This distributed graph is created by many
data providers representing and interrelating their data. What emerges from
this collective effort is a publicly accessible global database that can be
leverage by both man and machine to any end they deem appropriate. However,
the current instantiation of the Web of Data lacks one crucial component: a
distributed processing infrastructure. For the Web of Documents of common
knowledge, the solution to the issue of processing the vast amount of
information has been to literally download the entire Web and index and
process it in a single environment. While the content on the Web of Documents
is distributed, the means by which the information on the Web of Documents is
analyzed is not. The Web of Data need not fall into this same model. With the
nearly limitless ways in which RDF data can be processed, it would be a
disappointment if the data on the Web of Data was left solely to centralized
repositories to store, index, and provide query functionality. Beyond
disappointment, it would reduce the potential utility the Web of Data would
have given a distributed process infrastructure. By extending the work of the
Linked Data community with Linked Process, the Web of Data may one day rise to
become the defacto medium for representing and processing data much like the
Web of Documents is the defacto medium for storing and sharing documents.

## References

  * (1) Tim Berners-Lee, Roy T. Fielding, Day Software, Larry Masinter, and Adobe Systems.  Uniform Resource Identifier (URI): Generic Syntax, January 2005.
  * (2) Christian Bizer, Tom Heath, Kingsley Idehen, and Tim Berners-Lee.  Linked data on the web.  In Proceedings of the International World Wide Web Conference, Linked Data Workshop, Beijing, China, April 2008.
  * (3) Eric Miller.  An introduction to the Resource Description Framework.  Bulletin of the American Society for Information Science and Technology, 25(1):15–19, November 1998.
  * (4) Paul Miller, Rob Styles, and Tom Heath.  Open Data commons: A license for Open Data.  In Proceedings of the Workshop on Linked Data on the Web, New York, NY, April 2008. ACM Press.
  * (5) Eric Prud’hommeaux and Andy Seaborne.  SPARQL query language for RDF.  Technical report, World Wide Web Consortium, October 2004.
  * (6) Marko A. Rodriguez.  Data Management in the Semantic Web, chapter Interpretations of the Web of Data.  Nova Publishing, 2009.
  * (7) Marko A. Rodriguez.  Emergent Web Intelligence, chapter General-Purpose Computing on a Semantic Network Substrate.  Springer-Verlag, Berlin, DE, 2009.
  * (8) Marko A. Rodriguez and Jennifer H. Watkins.  Faith in the algorithm, part 2: Computational eudaemonics.  In Juan D. Velásquez, Robert J. Howlett, and Lakhmi C. Jain, editors, Proceedings of the International Conference on Knowledge-Based and Intelligent Information & Engineering Systems, Lecture Notes in Artificial Intelligence. Springer-Verlag, 2009.