[1]Haruya Kokubo

These authors contributed equally to this work.

[1]Department of Physics, Kindai University, 3-4-1 Kowakae, Higashi Osaka-
city, 577-8502, Osaka, Japan

# Impact of density inhomogeneity on the critical velocity for vortex shedding
in a harmonically trapped Bose-Einstein condensate

*<EMAIL_ADDRESS>Kenichi Kasamatsu<EMAIL_ADDRESS>*

###### Abstract

We report on a numerical study of the critical velocity for creation of
quantized vortices by a moving Gaussian obstacle in a trapped Bose-Einstein
condensate, modeled by the Gross-Pitaevskii equation. We pay attention to
impact of density inhomogeneity associated with the global inverted-parabolic
profile by a trapping potential as well as the local density suppression
around the Gaussian obstacle. When the width of the Gaussian potential is
large, the wake dynamics is significantly influenced by the nonuniformity
around the obstacle potential. The critical velocity, estimated through the
time interval between the first and second vortex emission, can be explained
by the local sound velocity by taking into account the above two
contributions. We also find that the ratio of the critical velocity to the
sound velocity at the center of the system is insensitive to the nonlinear
coefficient of the Gross-Pitaevskii equation, which supports the universal
discussion even in a inhomogeneous trapped condensate under the local density
approximation.

###### keywords:

Bose-Einstein condensate, superfluidity, critical velocity, quantized vortex

##  1 INTRODUCTION

Hydrodynamic instability is one of the key research topics in a fluid dynamics
when considering problems of a turbulent transition of a wake flow, which
occurs behind an obstacle moving in a fluid [1]. In classical hydrodynamics,
the characteristics of wake dynamics can be classified by the dimensionless
Reynolds number defined as a ratio of an inertial force and a viscous one.

A wake flow has been also studied very well in superfluids, having been
realized experimentally in cold atomic-gas Bose-Einstein condensates (BECs)
[2, 3, 4, 5, 6] or exciton-polariton condensates [7, 8, 9, 10]. The wake is
achieved in these systems by inducing a condensate flow in the presence of an
obstacle potential created by a localized external field. There have been lots
of theoretical works concerning wake superflows based on the Gross-Pitaevskii
(GP) model [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. The fundamental
interest in these works is to reveal similarity and difference from the
classical counterpart; characteristics in superfluids are mainly due to the
presence of the critical velocity of superfluidity and of vortices with
quantized circulation. These features yield a new platform to study a quantum
analogue of the von Kàrmàn vortex street [5, 16, 23] and the Reynolds number
[14].

To study the wake phenomenon in superfluids, a motion of an obstacle in a
static superfluid frame must be above a critical velocity, below which an
obstacle motion is frictionless. However, since the critical velocity is
determined by various details, e.g. the shape or the boundary condition, of
the obstacle, the systematic determination of the critical velocity is an
unresolved issue. The well-known Landau criterion of superfluidity predicts
that the critical velocity in a uniform system is given by the sound velocity
$c_{s}$ by considering the energetic instability for elementary excitations.
When there is a hard wall cylinder in incompressible flow, the local velocity
is multiplied by a factor 2 on the lateral side of the cylinder, so that the
critical velocity is expected to be $0.5c_{s}$. In the calculation of the GP
model with contributions of compressibility and quantum pressure, the critical
velocity is $0.37c_{s}$ when the size of the cylinder is sufficiently large
compared to the healing length [24, 25]. In the limit of the vanishing
obstacle size, the critical velocity approaches to $c_{s}$ [24, 25].

While the above-mentioned works have been discussed for a uniform system, the
atomic-gas BECs are essentially inhomogeneous [2, 3, 4, 5, 6]. There are two
external contributions that induce inhomogeneity, namely a trap potential and
an obstacle potential. For a BEC trapped in a harmonic potential, the global
form of the density profile is an inverted parabola. In addition, the
obstacles are usually modeled by the Gaussian form, which intrinsically
involves inhomogeneous distribution of the condensate density around the
localized obstacles. The previous experimental works [4] clarified that the
critical velocity decreases with increasing the width of the Gaussian
obstacle. One possible reason for this decrease could be the appearance of an
inhomogeneous region around the Gaussian obstacle.

In a uniform system, the GP equation is usually scaled by the bulk density and
the associated healing length, so that the results are free from the coupling
constant of the nonlinear term. For trapped systems, however, the dependence
of the coupling constant is nontrivial because of the additional length scale
due to the harmonic trap. By reducing the nonlinear coefficient, i.e., the
particle number or the (repulsive) interaction energy, the global size of the
condensate decreases while the healing length increases, the inhomogeneity
thus being reinforced. It is noticeable that the linear Schrödinger limit
gives rise to the vanishing critical velocity. Thus, it is expected that the
inhomogeneity basically plays a role of decreasing the critical velocity.

In this study, we investigate the critical velocity of the two-dimensional BEC
trapped in a harmonic potential for vortex generation by a moving Gaussian
obstacle. To focus on the inhomogeneous effects of vortex generation, we
consider the GP model with the nonlinear coefficients smaller than those
corresponding to the previous experiments [4]. Numerical simulations of the GP
equation demonstrate that the ratio of the critical velocity to the sound
velocity is smaller than previous literatures. This result suggests that the
critical velocity is influenced by the inhomogeneity of the system. We take
into account the local sound velocity due to the trap potential and the
Gaussian obstacle to determine which is responsible for the decrease in the
critical velocity.

The paper is organized as follows. In Sec.2 we introduce the GP equation with
an obstacle potential, discussing the inhomogeneous features of the system
caused by the trapping potential and the obstacle potential. In Sec.3 we
introduce the method to determine the critical velocity and analyze impact of
density inhomogeneity on the critical velocity. Finally, Sec.4 devotes to the
conclusions.

##  2 FORMULATION

Ignoring the inhomogeneity along the axial $z$-direction, we consider a two-
dimensional (2D) BEC in the $x$-$y$ space. Within a mean-field theory, a BEC
at low temperatures is described by the macroscopic wave function
$\Psi(\bm{r},t)=\sqrt{n(\bm{r},t)}e^{i\theta(\bm{r},t)}$ with the particle
density $n$ and the phase $\theta$. The dynamics of the wave function obeys
the GP equation

$i\hbar\frac{\partial}{\partial
t}\Psi=\left(-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{\rm{ext}}+g|\Psi|^{2}\right)\Psi.$
(1)

Here, $m$ is the atomic mass, $V_{\rm ext}(\bm{r})$ is the harmonic trap
potential and $g$ in the nonlinear term is the coupling constant in the 2D
system with the axial contribution being integrated out [14, 17].

###  2.1 Inhomogeneity by a trapping potential

First, we confirm the effect by the trapping potential. The stationary state
of the system can be calculated by inserting the form
$\Psi(\bm{r},t)=\Phi(\bm{r})e^{-i\mu t/\hbar}$ into Eq.(1) and by solving the
time-independent GP equation

$\left[-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{\rm
ext}-\mu+g|\Phi|^{2}\right]\Phi=0.$ (2)

A stationary density profile in a harmonic potential $V_{\rm
ext}(\bm{r})=\frac{1}{2}m\omega^{2}(x^{2}+y^{2})$ is written within the
Thomas-Fermi approximation as [26]

$n_{\rm TF}(\bm{r})=n_{0}\left(1-\frac{x^{2}}{R^{2}_{\rm
TF}}-\frac{y^{2}}{R^{2}_{\rm TF}}\right).$ (3)

Here, the Thomas-Fermi radius is $R_{\rm TF}=\sqrt{{2\mu}/{m\omega^{2}}}$ and
the density at the trap center is $n_{0}=n(0)={\mu}/{g}$. From the
normalization condition $\int d\bm{r}n_{\rm TF}(\bm{r})=N$, a relationship
between the coupling constant $g$ and the chemical potential $\mu$ is given by

$g=\frac{\pi\mu^{2}}{mN\omega^{2}}.$ (4)

Here, $N$ is a particle number in the 2D system. We regard the chemical
potential $\mu$ as a parameter to represent the coupling constant of the GP
equation instead of $g$, which is useful in the following discussion because
the relevant energy scales are often compared with $\mu$. We also confine
ourselves to the Thomas-Fermi regime $\mu\gg\hbar\omega$ or, equivalently,
$R_{\text{TF}}\gg\xi$ with the healing length $\xi=\hbar/\sqrt{2m\mu}$.

###  2.2 Inhomogeneity by an obstacle potential

We next consider a 2D trapped BEC with a moving obstacle potential at a
constant velocity. This system is described by the time-dependent GP equation

$i\hbar\frac{\partial}{\partial
t}\Psi=\left[-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{\rm{ext}}+V_{\rm{obst}}+g|\Psi|^{2}\right]\Psi.$
(5)

The obstacle potential $V_{\rm obst}(\bm{r},t)$ is given by the Gaussian form
as

$V_{\rm obst}=V_{0}{\rm exp}\\{-2[(x-x_{0}-v_{\rm
obst}t)^{2}+y^{2}]/\sigma^{2}\\},$ (6)

which moves at a constant velocity $v_{\rm obst}$ along the $x$ axis. The
parameters $V_{0}$, $x_{0}$, and $\sigma$ are the height, the initial
position, and the size of the Gaussian potential, respectively. In this work,
the length, time and energy scales are $\sqrt{\hbar/m\omega}\equiv a_{h}$,
$1/\omega\equiv\tau$ and $\hbar\omega$, respectively, by using units of a
harmonic potential. Then, the dimensionless GP equation with a moving Gaussian
obstacle is written as

$i\frac{\partial}{\partial\tilde{t}}\tilde{\Psi}=\left[-\frac{1}{2}\tilde{\nabla}^{2}+\tilde{V}_{\rm
ext}+\tilde{V}_{\rm
obst}+\pi\tilde{\mu}^{2}|\tilde{\Psi}|^{2}\right]\tilde{\Psi}.$ (7)

Here, the quantities with tildes are dimensionless and the wave function is
scaled as $\tilde{\Psi}=\Psi(\bm{r},t)(a_{h}/\sqrt{N})$. The coupling constant
$g$ is replaced by $\mu$ through the relation of Eq.(4); although Eq.(4) is
not correct in the presence of $V_{\text{obst}}$, the resultant correction is
negligibly small in our parameter setting as discussed below. When the
amplitude $V_{0}$ of the Gaussian obstacle is smaller than the chemical
potential $\mu$, the wave function can penetrate into the obstacle’s area [4].
In this work, we focus on the impenetrable regime $V_{0}>\mu$, in which the
condensate density strongly decays at the boundary of the obstacle and
vanishes inside the obstacle. The stationary solution with the Gaussian
obstacle is obtained by using the imaginary time evolution method for Eq.(7).

Figure 1: The cross sections of the stationary density $n_{0}(-0.6R_{\rm
TF},y)$ along the $y$-direction passing through the center of the Gaussian
obstacle. (a) The density profiles for $\mu/(\hbar\omega)=50$ and different
values of the Gaussian widths $\sigma/\xi=6$ (solid curve), $10$ (dashed
curve) and $14$ (dotted curve). (b) The enlarged plot of (a) around the
Gaussian obstacle. The circles, triangles, and squares correspond to the
density at $y=R_{\rm obst}$, $R_{\rm obst}+\xi$, and $R_{\rm obst}+2\xi$,
respectively, determined in Sec. 3.2. (c) A similar plot of (b), but for fixed
$\sigma/\xi=10$ and different values of $\mu/(\hbar\omega)=20$ (solid curve),
50 (dashed curve), 150 (dotted curve).

Figure 1 shows the cross section of the stationary density profile along the
$y$-direction across the central position of the Gaussian obstacle, where we
set $x_{0}=-0.6R_{\text{TF}}$ and represent the value of $\sigma$ in units of
$\xi$, following the calculation in the next section. As shown in Fig.1(a),
with increasing the Gaussian width $\sigma$ for fixed $\mu$, the size of the
density hole, created by Gaussian obstacle, is obviously increased. On the
other hands, the outer Thomas-Fermi boundary is not influenced so much by the
obstacle’s size, which supports our statement below Eq.(7). It is noteworthy
that, as clearly seen in Fig.1(b), the Gaussian tail becomes larger with
increasing $\sigma$, extending the “recovery region” of the density from zero
to the equilibrium Thomas-Fermi value. This feature indicates that the density
inhomogeneity is more apparent for a larger Gaussian obstacle and has a strong
impact to the suppression of the critical velocity as discussed below. Also,
the recovery region is extended with decreasing the chemical potential $\mu$
even for fixed $\sigma$, as seen in Fig.1(c). This is because, the smaller
$\mu$ becomes, the longer the healing length $\xi$ is.

##  3 RESULTS

We investigate the effects of inhomogeneous density to the wake dynamics and
the associated critical velocity by the Gaussian obstacle potential and the
trapping potential through the numerical calculation of Eq.(7). Since the
sound velocity depends on the local density as
$c_{s}(\bm{r})=\sqrt{gn(\bm{r})/m}$ within the local density approximation,
the sound velocity $c_{s}(\bm{r}=0)\equiv c_{s}$ at the center of the system
is used as the basis for unifying the scale of the critical velocity of the
system. The initial position of the obstacle is set as $x_{0}=-0.6R_{\rm TF}$,
and the center of the obstacle moves to the $+x$-direction with the distance
$2|x_{0}|$. The mesh sizes $\Delta_{x,y}$ are determined by
$\Delta_{x,y}=\xi/4$ and the time step is $\Delta t=5.0\times 10^{-6}\tau$. In
this study, the chemical potential $\mu$ and the width $\sigma$ of the
Gaussian potential are smaller than those of the experiment [4]. The value of
$V_{0}$ in Sec.3.2 is taken to be the same value as that of experiment for
comparison with the experiment [4].

Figure 2: The simulation result of the vortex shedding for
$\mu/(\hbar\omega)=20$, $\sigma/\xi=4.0$, $V_{0}=7\mu$ and $v_{\rm
obst}=0.4c_{s}$. The panels (a) and (b) show the time evolution of the
condensate density $n(\bm{r},t)=|\Psi(\bm{r},t)|^{2}$ and phase
$\theta(\bm{r},t)=\text{arg}[\Psi(\bm{r},t)]$, respectively. The spatial
region of the plot is $-7a_{h}\leq x,y\leq 7a_{h}$ .

Figures 2(a) and (b) show snapshots of the condensate density $n(\bm{r},t)$
and the phase $\theta(\bm{r},t)$ for $\mu/(\hbar\omega)=20$, $\sigma/\xi=4.0$,
$V_{0}/\mu=7$ and $v_{\rm obst}=0.4c_{s}$, which is expected to be above the
critical velocity. When the obstacle begins to move, the obstacle pushes
forward through the area of high density in the central region, inducing the
collective dipole motion of the entire system as well as the small-amplitude
density wave. Subsequently, vortex nucleation occurs at the lateral sides of
the moving obstacle. From the phase profile in Fig. 2(b), this is a vortex-
antivortex pair created initially in the density depleted region caused by the
obstacle potential. This is slightly different from the result reported in
Ref. [15], where two pairs of vortices appear initially inside the density
depleted region in the impenetrable regime. This difference is due to the fact
that our obstacle size is rather small $\sigma/\xi=7$ compared to
$\sigma/\xi=20$ in Ref. [15]. After the vortex pair is released from the
obstacle boundary to the finite-density region, it moves to follow the
obstacle’s motion. Again, a new vortex pair is generated inside the density
depleted region and released before the obstacle pass through the Thomas-Fermi
boundary. In our simulations, emission of vortex pairs occurs up to twice
restricted by the finite system size. It is noticeable that vortex nucleation
and emission events tend to take place in the low-density region instead of
the central high-density region.

###  3.1 Method of determining the critical velocity

Our finite size system has some problems for the determination of the critical
velocity. Here, we present our strategy to determine the critical velocity
according to the emission frequency of quantized vortices from an obstacle
[6].

When an obstacle moves through a condensate near the critical velocity $v_{\rm
obst}\gtrsim v_{c}$, vortex dipoles are periodically generated [27, 4, 28]. In
the previous studies, the phase-slip event [17] or the rapid change in a drag
force acting on an obstacle [14] have been proposed as methods to measure the
timing of vortex emission in the numerical simulations. In our case, since the
density inside the hard potential is almost vanished, the phase fluctuation
arises easily in such an extremely low-density region so that the phase slip
cannot be identified clearly. Therefore, we determine the time of vortex
emissions from the drag force in the $x$-direction caused by the obstacle:

$\langle F(t)\rangle=-\int dxdy\left(\frac{\partial}{\partial x}V_{\rm
obst}\right)|\Psi|^{2}.$ (8)

Figure 3(a) shows the time development of the expected value of the drag force
of Eq.(8). Contrary to the calculations in homogeneous systems [14], it is
difficult to obtain the exact time of vortex creation and emission, since the
clear extremes relevant to events of the vortex emission cannot be identified.
This is because the signal of the vortex emission is buried by the inevitable
phonon emission and the density compression. To emphasize the fine changes of
the time evolution of $\langle F(t)\rangle$, we take the time derivative
$d\langle F(t)\rangle/dt$. Together with the careful check of the snapshots of
the density evolution, we identify the timing of the vortex emission from the
obstacle, as shown by arrows in Fig.3(b).

A relationship between the obstacle’s velocity and the period of vortex
emission $T$ is empirically given by

$v_{\rm obst}-v_{\rm c}\propto\frac{1}{T}$ (9)

for an obstacle in a penetrable regime [6]. In our impenetrable regime, we
assume that the critical velocity can be also obtained by this relationship.
However, it is difficult to measure such a period $T$ exactly in our finite-
sized system because a number of events involving vortex emission is not
enough for measuring the period. Hence, we evaluate it approximately as
$T\simeq T_{\text{2nd}}-T_{\text{1st}}$ using the interval between the first
and second vortex emissions under the assumption that the vortices would be
generated periodically. As seen in Fig.3, the interval has a linear dependence
with respect to $v_{\text{obst}}$, which enables us to estimate the critical
velocity.

Figure 3:  (a) The time evolution of the expected value of the drag force in
the $x$-direction acting on the obstacle. It is difficult to get the definite
time of the vortex emission from this data, where the arrows indicate the
timing of the vortex emission. (b) The time derivative of the drag force,
which enhances the subtle variation of the drag force. The timing of the
vortex emission is identified through the comparison with the animation of the
dynamics. The other peaks may be caused by the density waves reflected from
the Thomas-Fermi boundary. The panels (c) and (d) show the relations between
velocity of the obstacle $v_{\rm obst}$ and the interval of the vortex
emissions for $\mu/(\hbar\omega)=20$ and $50$. The plots in (c) show the
results for $\sigma/\xi=2(\blacksquare)$, $\sigma/\xi=3(\blacktriangle)$,
$\sigma/\xi=4(\blacklozenge)$, $\sigma/\xi=5(\bullet)$ and
$\sigma/\xi=6(\blacktriangledown)$. The plots in (d) show the results for
$\sigma/\xi=4(\blacklozenge)$,$\sigma/\xi=5(\blacktriangledown)$,
$\sigma/\xi=6(\blacktriangle)$, $\sigma/\xi=8(\blacksquare)$ and
$\sigma/\xi=10(\bullet)$. The value of the critical velocity is extracted at
the crossing point of the horizontal axis at
$(T_{\text{2nd}}-T_{\text{1st}})^{-1}=0$ by interpolating the linear fit of
the data (red-dashed lines).

###  3.2 Dependence on the size of the obstacle

First, we discuss the dependence of the critical velocity on the size of the
obstacle, which was studied in the previous works [29, 4]. The experiment in
Ref. [4] observed the nontrivial decreasing behavior of the critical velocity
with respect to the size of the obstacle. To investigate whether this behavior
is relevant to the inhomogeneity of the system, we simulate the GP equation
with a height of the potential $V_{0}=7\mu$ and two different values of the
chemical potential $\mu/(\hbar\omega)=20$ and $50$. Although these values are
much smaller than the estimated value $\mu\approx 140$ corresponding to the
experiment of Ref. [4], it is hard to make simulations with such a large value
because of the numerical cost to describe accurately both the large condensate
size and the small healing length. Nonetheless, we expect that our choice of
the parameter can capture the inhomogeneous effect since the obtained critical
velocity is not so sensitive to the chemical potential, as shown below.

Figure 4: Dependence of critical velocity on the obstacle size. The height of
the Gaussian obstacle is fixed as $V_{0}/\mu=7$. Vertical axis represents a
ratio of the calculated critical velocity $v_{c}$ to the sound velocity
$c_{s}$ estimated by the chemical potential at the trap center. Filled circles
and squares are results for $\mu/(\hbar\omega)=20$ and $50$, respectively.
Solid diamonds are experimental results taken from Fig. 4 in Ref.[4]. The
thin-solid and thin-dashed curves represent $v_{c}^{\text{TF}}$ of Eq. (11)
for $\alpha=1$ and 2, respectively. Here, the choice of the parameters are (A)
$(\mu,x_{0})=(20\hbar\omega,-0.6R_{\rm TF})$, (B)
$(\mu,x_{0})=(50\hbar\omega,-0.6R_{\rm TF})$, (C)
$(\mu,x_{0})=(150\hbar\omega,-0.05R_{\rm TF})$. The Bold-solid and Bold-dotted
curves represent $v_{c}^{\text{num}}$ of Eq. (12) for $\alpha=1$ and 2,
respectively.

Plots in Fig. 4 represent the critical velocity $v_{c}/c_{s}$ of the vortex
shedding, obtained by the method in Sec. 3.1, as a function of the width
$\sigma$ of the Gaussian obstacle. Although there are variations in the data,
the critical velocity decreases monotonically from $\sim 0.4$ with increasing
$\sigma/\xi$. For $\sigma/\xi>10$ the strong perturbation caused by the
obstacle with the large sizes prevents the clear identification of the vortex
nucleation event. It is noticed that our results decrease more rapidly than
the experimental result [4]. It has been known that the critical velocity of
the vortex shedding can be determined by the sound velocity multiplied by a
certain numerical factor of $\lesssim\mathcal{O}(1)$. Since a local change of
the sound velocity occurs in an inhomogeneous system, the critical velocity is
expected to be written by the local sound velocity $c_{s}^{\rm
loc}(\bm{r})=\sqrt{gn(\bm{r})/m}$ under the local density approximation. Here,
we focus on the two factors that cause the density inhomogeneity. The first
one is a decrease in the density caused by the inverted-parabolic Thomas-Fermi
profile. This decrease in density should be considered seriously since the
position at which the vortex nucleation occurs tends to locate far from the
trap center with increasing the size of the obstacle. The second one is the
density suppression caused by the tail of the Gaussian obstacle, which becomes
more remarkable as the size of the obstacle increases as seen in Fig. 1. The
latter further brings about the decrease in the critical velocity.

It is necessary to determine the location for evaluating the local sound
velocity. Since it is difficult to determine exactly the location
corresponding to the timing of the vortex emission, we consider the local
density around the obstacle potential at the “initial position”
$x=x_{0}(=-0.6R_{\text{TF}})$. This prescription is based on the numerical
observation seen in Fig. 2 that the first vortex emission occurs near the
initial position of the obstacle. Also, we suppose that the lateral sides of
the obstacle are the points where the vortex nucleation is easy to occur. To
determine this point, we consider the effective radius $R_{\text{obst}}$ of
the obstacle, which can be estimated by solving the equation

$\displaystyle V_{0}e^{-2R_{\text{obst}}^{2}/\sigma^{2}}$
$\displaystyle=\mu(x_{0},R_{\text{obst}})$ $\displaystyle\equiv\mu-
V_{\text{ext}}(x_{0},R_{\text{obst}}),$ (10)

which corresponds to the condition determining the Thomas-Fermi boundary and
is a generalization of the discussion in Ref.[4] to the local chemical
potential. Thus, this condition provides the well-defined boundary induced by
the Gaussian obstacle. Equation (10) is solved numerically to obtain the value
of $R_{\text{obst}}$, which is slightly enhanced about 10% for
$x_{0}=-0.6R_{\text{TF}}$ from the value at the trap center. Since the vortex
nucleation is likely to occur within the scale of the healing length from the
boundary, the local sound velocities are determined by using the local density
$n(x_{0},R_{\text{obst}}+\alpha\xi)$. Here,the constant
$\alpha\sim\mathcal{O}(1)$ is an uncertain numerical factor of order unity,
including a contribution of the vortex core structure as well as a correction
of the local healing length from the bulk value. We show the results with
$\alpha=1$ and 2 in the following.

We first consider only the inhomogeneous effect caused by the Thomas-Fermi
density profile. The critical velocity, scaled by the sound velocity at the
trap center, is related to the local sound velocity $c_{\rm TF}^{\rm loc}$ as

$v_{c}^{\text{TF}}=0.5\frac{c_{\rm TF}^{\rm
loc}(\bm{r})}{c_{s}}=0.5\sqrt{\frac{n_{\text{TF}}(x_{0},R_{\text{obst}}+\alpha\xi)}{n_{\text{TF}}(0,0)}},$
(11)

where we have multiplied the factor 0.5 since the local velocity at the
lateral sides of the cylindrical object becomes twice of the background
velocity of the fluid. Several curves of the thin-solid ($\alpha=1$) and the
thin-dashed ($\alpha=2$) lines in Fig. 4 show $v_{c}^{\text{TF}}$ as a
function of $\sigma/\xi$ for $\mu/(\hbar\omega)=20$ and $50$. For
$\sigma/\xi\ll 1$ the slight decrease in the density due to the inverted
parabolic profile from the center yields $v_{c}^{\text{TF}}\simeq 0.4$ smaller
than $0.5$, which is consistent to the experimental [4] and numerical
observations. However, these results do not exhibit the rapidly decreasing
behavior obtained in the numerical results. Also, we plot the result for the
experimentally relevant parameters $\mu/(\hbar\omega)=150$ and
$x_{0}=-0.05R_{\text{TF}}$ [4]. Then, the value of $v_{c}^{\text{TF}}$ is
almost constant around 0.5, whose behavior is largely deviated from the
experimental observation.

It is clearly seen that the density around the obstacle potential is further
suppressed from the Thomas-Fermi profile, which leads to additional
suppression of the critical velocity. As a next attempt, we consider the
effect of density suppression caused by the Gaussian tail, which is more
pronounced as the size of the obstacle increases. To this end, we employ the
local density $n_{\rm num}(x_{0},R_{\text{obst}}+\alpha\xi)$ taken from the
numerical solutions of the GP equation, as indicated in Fig. 1 for example, to
evaluate the local sound velocity as

$v_{c}^{\text{num}}=0.5\frac{c_{\rm num}^{\rm
loc}(\bm{r})}{c_{s}}=0.5\sqrt{\frac{n_{\text{num}}(x_{0},R_{\text{obst}}+\alpha\xi)}{n_{\text{num}}(0,0)}}.$
(12)

The results are shown by the curves of bold-line ($\alpha=1$) and dotted-line
($\alpha=2$) in Fig.4(b) and (c). We find that the obtained curves for the
particular value of $x_{0}$ are almost independent of the values of $\mu$,
thus showing the results only for $\mu/(\hbar\omega)=50$ in the case of
$x_{0}=-0.6R_{\rm TF}$. The result indicates that including the density
suppression by the Gaussian tail can reproduce the rapid decreasing behavior
of the critical velocity with respect to $\sigma$. We see that the numerical
results are fitted better for $\alpha=2$. We also find that the decreasing
behavior obeys the power low $\sim(\sigma/\xi)^{1/3}$, which is a decay slower
than the predictions of the previous studies [15, 30]. Furthermore, we plot
the case of the experimentally relevant parameters ($\mu/(\hbar\omega)=150$
and $x_{0}=-0.05R_{\rm TF}$). Since the initial position of the Gaussian
obstacle is near the center of the condensate, the estimated critical velocity
is larger than that for $\mu/(\hbar\omega)=50$. However, the values
underestimate the experimental results.

We give two comments before enclosing this section.

  1. 1.

The value of $v_{c}^{\text{num}}$ does not depend on the chemical potential
$\mu$, while the numerical results indicate that $v_{c}$ for
$\mu/(\hbar\omega)=20$ is clearly smaller than that for
$\mu/(\hbar\omega)=50$. One possible reason of this discrepancy may be
considered as follows. The flow induced by the moving obstacle may results in
the backflow due to the reflection of the density wave from the outer Thomas-
Fermi boundary. This backflow enhances the relative velocity between the
condensate and the obstacle, resulting in the decrease in $v_{c}$. For a small
condensate with $\mu/(\hbar\omega)=20$ and a fast moving obstacle, this effect
cannot be negligible.

  2. 2.

Our result seems not to be consistent with the experimental observation in
Ref. [4]. In the experiment [4], the critical velocity was measured by seeing
a single event whether the vortices appear or not through the displacement
$\Delta x=0.1R_{\rm TF}$, which seems to be small compared with the whole
condensate size. As shown in the next section, this experimental protocol may
lead to overestimation of the critical velocity compared with that with the
method in Sec. 3.1 using the period of the vortex emission.

###  3.3 Dependence on the nonlinear coefficient

As shown in the previous section, the critical velocity, scaled by the sound
velocity at the center, is insensitive to the chemical potential i.e., the
nonlinear coefficient of the GP equation, even in the trapped system. Here, as
supplemental results in the previous section, we investigate the dependence of
the critical velocity on the nonlinear coefficient of the GP equation. We fix
$V_{0}=10\mu$ and the obstacle size $\sigma/\xi=0.6\sqrt{5}\sim 1.34$ as small
as possible. The latter choice is to reduce the effect of the Gaussian
inhomogeneity described in Sec.3.2. Figure 5(a) shows the critical velocity as
a function of $\mu$. Here, in addition to the method in Sec. 3.1, we give an
alternative estimate of the critical velocity by seeing whether the event of
the vortex emission from the obstacle occurs. In fact, the latter method
follows the protocol used in the experiment of Ref.[4]. The two methods
confirm that the critical velocity is insensitive to the value of $\mu$. These
results are evident from the density profile of the initial state in Fig.1;
the local density $n_{\text{num}}(x_{0},R_{\text{obst}}+\alpha\xi)$ decreases
as $\sigma$ increased for the fixed $\mu$, while it is almost unchanged by
varying $\mu$ for the fixed $\sigma$. The critical velocity obtained by the
latter method is always larger than the former, since this method does not
eliminate the possibility of vortex generation in the successive time
evolution. This might be a possible reason why the critical velocities in the
experiment [4] are always larger than our numerical results.

Figure 5: Dependence of critical velocity on the chemical potential $\mu$ for
$V_{0}=10\mu$ and $\sigma/\xi=0.6\sqrt{5}\sim 1.34$. In (a), the vertical axis
is a ratio of the critical velocity $v_{c}$ to the sound velocity $c_{s}$ at
the center of the system. The squares are obtained by the method in Sec.3.1,
while the circles are obtained by seeing whether the vortex emission occurs in
the simulations. (b) The time duration of the first event of the vortex
emission as a function of $\mu$.

The method in Sec. 3.1 is not applicable for small $\mu\lesssim 15$ since the
time duration for which the vortex nucleation occurs is extended, as shown in
Fig. 5(b). Also, the result in Fig. 5(b) implies that vortex nucleation is
less likely to occur as the chemical potential decreases. In a uniform system,
the dynamics is free from the chemical potential since there are no parameters
in the scaled GP equation. Thus, the tendency observed in Fig.5(b) can be
interpreted as an inhomogeneous effect, indicating that the vortex nucleation
is not easy to occur for small $\mu$ since the healing length is extended over
the system so that the local vortex nucleation may be prohibited.

##  4 Conclusion

In this study, we investigated the influence of inhomogeneity of the trapped
Bose-condensed system on the critical velocity for vortex shedding via a
Gaussian obstacle. The decrease in the critical velocity from one estimated in
a homogeneous bulk sound velocity can be explained by the local sound velocity
due to not only the inverted parabolic density profile in a trap potential but
also the density suppression by the tail of the Gaussian potential. We also
confirmed that the critical velocity is not dependent on the value of the
chemical potential, which indicate that universal discussion in a uniform
system is approximately applicable to the trapped system by properly including
the inhomogeneous feature within the local density approximation.

## References

  *   * Landau [1941] Landau, L.: Theory of the Superfluidity of Helium II. Physical Review 60(4), 356–358 (1941)
  * Raman et al. [1999] Raman, C., Köhl, M., Onofrio, R., Durfee, D.S., Kuklewicz, C.E., Hadzibabic, Z., Ketterle, W.: Evidence for a Critical Velocity in a Bose-Einstein Condensed Gas. Physical Review Letters 83(13), 2502–2505 (1999)
  * Onofrio et al. [2000] Onofrio, R., Raman, C., Vogels, J.M., Abo-Shaeer, J.R., Chikkatur, A.P., Ketterle, W.: Observation of Superfluid Flow in a Bose-Einstein Condensed Gas. Physical Review Letters 85(11), 2228–2231 (2000)
  * Kwon et al. [2015] Kwon, W.J., Moon, G., Seo, S.W., Shin, Y.: Critical velocity for vortex shedding in a Bose-Einstein condensate. Physical Review A 91(5), 053615 (2015)
  * Kwon et al. [2016] Kwon, W.J., Kim, J.H., Seo, S.W., Shin, Y.: Observation of von kármán Vortex Street in an Atomic Superfluid Gas. Physical Review Letters 117(24), 245301 (2016)
  * Lim et al. [2022] Lim, Y., Lee, Y., Goo, J., Bae, D., Shin, Y.: Vortex shedding frequency of a moving obstacle in a Bose-Einstein condensate. New Journal of Physics 24(8), 083020 (2022)
  * Amo et al. [2009] Amo, A., Lefrére, J., Pigeon, S., Adrados, C., Ciuti, C., Carusotto, I., Houdré, R., Giacobino, E., Bramati, A.: Superfluidity of polaritons in semiconductor microcavities. Nature Physics 5(11), 805–810 (2009)
  * Amo et al. [2011] Amo, A., Pigeon, S., Sanvitto, D., Sala, V.G., Hivet, R., Carusotto, I., Pisanello, F., Leménager, G., Houdré, R., Giacobino, E., Ciuti, C., Bramati, A.: Polariton Superfluids Reveal Quantum Hydrodynamic Solitons. Science 332(6034), 1167–1170 (2011)
  * Nardin et al. [2011] Nardin, G., Grosso, G., Léger, Y., Piȩtka, B., Morier-Genoud, F., Deveaud-Plédran, B.: Hydrodynamic nucleation of quantized vortex pairs in a polariton quantum fluid. Nature Physics 7(8), 635–641 (2011)
  * Lerario et al. [2017] Lerario, G., Fieramosca, A., Barachati, F., Ballarini, D., Daskalakis, K.S., Dominici, L., De Giorgi, M., Maier, S.A., Gigli, G., Kéna-Cohen, S., Sanvitto, D.: Room-temperature superfluidity in a polariton condensate. Nature Physics 13(9), 837–841 (2017)
  * Frisch et al. [1992] Frisch, T., Pomeau, Y., Rica, S.: Transition to dissipation in a model of superflow. Physical Review Letters 69(11), 1644–1647 (1992)
  * Nore et al. [2000] Nore, C., Huepe, C., Brachet, M.E.: Subcritical Dissipation in Three-Dimensional Superflows. Physical Review Letters 84(10), 2191–2194 (2000)
  * Aftalion et al. [2003] Aftalion, A., Du, Q., Pomeau, Y.: Dissipative Flow and Vortex Shedding in the Painlevé Boundary Layer of a Bose-Einstein Condensate. Physical Review Letters 91(9), 090407 (2003)
  * Reeves et al. [2015] Reeves, M.T., Billam, T.P., Anderson, B.P., Bradley, A.S.: Identifying a Superfluid Reynolds Number via Dynamical Similarity. Physical Review Letters 114(15), 155302 (2015)
  * Kwak et al. [2023] Kwak, H., Jung, J.H., Shin, Y.: Minimum critical velocity of a gaussian obstacle in a Bose-Einstein condensate. Phys. Rev. A 107(2), 023310 (2023)
  * Stagg et al. [2014] Stagg, G.W., Parker, N.G., Barenghi, C.F.: Quantum analogues of classical wakes in Bose–Einstein condensates. Journal of Physics B: Atomic, Molecular and Optical Physics 47(9), 095304 (2014)
  * Jackson et al. [1998] Jackson, B., McCann, J.F., Adams, C.S.: Vortex Formation in Dilute Inhomogeneous Bose-Einstein Condensates. Physical Review Letters 80(18), 3903–3906 (1998)
  * Winiecki et al. [2000] Winiecki, T., Jackson, B., McCann, J.F., Adams, C.S.: Vortex shedding and drag in dilute Bose-Einstein condensates. Journal of Physics B: Atomic, Molecular and Optical Physics 33(19), 4069–4078 (2000)
  * Musser et al. [2019] Musser, S., Proment, D., Onorato, M., Irvine, W.T.M.: Starting flow past an airfoil and its acquired lift in a superfluid. Phys. Rev. Lett. 123(15), 154502 (2019)
  * Kiehn et al. [2022] Kiehn, H., Singh, V.P., Mathey, L.: Superfluidity of a laser-stirred Bose-Einstein condensate. Phys. Rev. A 105(4), 043317 (2022)
  * Kunimi and Kato [2015] Kunimi, M., Kato, Y.: Metastability, excitations, fluctuations, and multiple-swallowtail structures of a superfluid in a Bose-Einstein condensate in the presence of a uniformly moving defect. Phys. Rev. A 91(5), 053608 (2015)
  * Pinsker and Berloff [2014] Pinsker, F., Berloff, N.G.: Transitions and excitations in a superfluid stream passing small impurities. Phys. Rev. A 89(5), 053605 (2014)
  * Sasaki et al. [2010] Sasaki, K., Suzuki, N., Saito, H.: Bènard–von Kàrmàn Vortex Street in a Bose-Einstein Condensate. Physical Review Letters 104(15), 150404 (2010)
  * Rica [2001] Rica, S.: A remark on the critical speed for vortex nucleation in the nonlinear Schrödinger equation. Physica D: Nonlinear Phenomena 148(3), 221–226 (2001)
  * Pham et al. [10 15, 2005] Pham, C.-T., Nore, C., Brachet, M.-É.: Boundary layers and emitted excitations in nonlinear Schrödinger superflow past a disk. Physica D: Nonlinear Phenomena 210(3), 203–226 (10 15, 2005)
  * Baym and Pethick [1996] Baym, G., Pethick, C.J.: Ground-State Properties of Magnetically Trapped Bose-Condensed Rubidium Gas. Physical Review Letters 76(1), 6–9 (1996)
  * Neely et al. [2010] Neely, T.W., Samson, E.C., Bradley, A.S., Davis, M.J., Anderson, B.P.: Observation of Vortex Dipoles in an Oblate Bose-Einstein Condensate. Physical Review Letters 104(16), 160401 (2010)
  * Kwon et al. [2015] Kwon, W.J., Seo, S.W., Shin, Y.-i.: Periodic shedding of vortex dipoles from a moving penetrable obstacle in a Bose-Einstein condensate. Physical Review A 92(3), 033613 (2015)
  * Huepe and Brachet [2000] Huepe, C., Brachet, M.-E.: Scaling laws for vortical nucleation solutions in a model of superflow. Physica D: Nonlinear Phenomena 140(1), 126–140 (2000)
  * Stießberger and Zwerger [2000] Stießberger, J.S., Zwerger, W.: Critcal velocity of superfluid flow past large obstacles in Bose-Einstein condensates. Physical Review A 62(6), 061601 (2000)


# High frequency oscillations in spin-torque nano oscillator due to bilinear
coupling

R. Arun  Department of Nonlinear Dynamics, School of Physics, Bharathidasan
University, Tiruchirapalli-620024, India  R. Gopal  Department of Physics,
Centre for Nonlinear Science and Engineering, School of Electrical and
Electronics Engineering, SASTRA Deemed University, Thanjavur 613 401, India
V. K. Chandrasekar  Department of Physics, Centre for Nonlinear Science and
Engineering, School of Electrical and Electronics Engineering, SASTRA Deemed
University, Thanjavur 613 401, India  M. Lakshmanan  Department of Nonlinear
Dynamics, School of Physics, Bharathidasan University, Tiruchirapalli-620024,
India

###### Abstract

Exchange coupling in an interfacial context is crucial for spin-torque nano
oscillator (STNO) that consists of a non-magnetic spacer which is alloyed with
a ferromagnetic material. Currently, investigations on the dynamics of the
free layer magnetization and frequency enhancement in the STNO with bilinear
coupling are still being actively pursued. In the present work, we investigate
the dynamics of the STNO in the presence of bilinear coupling but in the
absence of an external magnetic field by analyzing the associated Landau-
Lifshitz-Gilbert-Sloncewski(LLGS) equation, and consequently the impact of the
bilinear coupling on the dynamics of the magnetization of the free layer is
studied. It is observed that the frequency of the oscillations in the
magnetization component along the direction of the pinned layer polarization
can be enhanced above 300 GHz by positive bilinear coupling and up to around
30 GHz by negative bilinear coupling. We further reveal a transition from in-
plane to out-of-plane precession both for positive and negative bi-linear
couplings. We also analyze the switching of the magnetization for different
values of current and bilinear coupling. Our detailed investigations of STNO
with bilinear coupling aim at the possibilities of high-frequency devices by
considering the applied current and bilinear coupling in the absence of a
magnetic field.

##  I Introduction

A spin-polarized electrical current can impart spin angular momentum in the
ferromagnetic material, which can be used to control the magnetization state
of a magnetoresistive device called spin torque nano oscillator (STNO) slon ;
berger ; slon1 ; katine ; deng ; yama ; liu ; imai ; wu ; wolba ; eklund ;
jiang ; luba  . In particular, it is feasible to cause the oscillations or
precession of the magnetization, which is relevant for tunable microwave
devices or to reverse the magnetization that is essential for various magnetic
memory systems groll . In an STNO, two ferromagnetic layers are separated by a
thin nonmagnetic, but conductive layer called a spacer. Among the two
ferromagnetic layers, one is called the free layer, which is comparatively
thinner than the other which is the pinned layer. In the free layer the
direction of magnetization can change while it is fixed in the pinned layer.
Further, some studies also ensure that the spacer layer can promote a high
interlayer exchange coupling between its adjacent ferromagnetic layers parkin
. The bottom and top layers of the two in-plane magnetized ferromagnetic
layers are exchange-coupled via a Ruderman-Kittel-Kasuya-Yosida (RKKY)
interaction across the thin nonmagnetic spacer, whose thickness is tuned to
produce an antiferromagnetic coupling in zero applied field parkin ; ungu ;
naga ; moser ; gri .

For instance, a nonmagnetic layer typically made of Ru mck  introduces a RKKY
exchange coupling between two magnetic layers mck . The spin direction of the
ferromagnetic layers can be parallel or antiparallel to each other depending
upon the thickness of the spacer layer in magnetic multilayer systems. This
parallel or antiparallel orientation of the ferromagnetic layers can be called
collinear magnetization configuration mck ; duine . On the other hand,
obtaining a noncollinear magnetization configuration is possible due to the
competition between the interlayer coupling energy and magnetic anisotropies
of the coupled ferromagnetic layers for some structures. Recently, Nunn et al.
have reported that the influence of the exchange coupling between two
ferromagnetic layers (Fe) coupled through a nonmagnetic interlayer (Ru) is
essential in controlling the magnetic layers’ functionality nunn , and this
has now been observed in various systems. It has been explained theoretically
by several different approaches cuc ; rana ; harris ; grim ; meena ; tsu ;
arun ; kuro ; ishi .

Recent results nunn ; ishi ; kuro ; jiang ; luba ; liu  in this context show
that the presence of an exchange coupling system plays a backbone in the
emergence of many spintronic-based applications such as magnetic field
sensors, magnetic memory devices rana ; harris , magnetic resistive random
access memory (MRAM) meena  and spin-torque nano oscillators slon ; slon1 ;
berger . Based on the nanoscale size and suitability for room-temperature
operation, spin-torque oscillators (STOs) provide exciting possibilities for
these applications. However, their adjustable range and oscillation frequency
are only from 100 MHz to 10 GHz grim ; tsu . Recently, we investigated and
reported that the frequency of an STNO with bilinear and bi-quadratic
couplings can be enhanced above 300 GHz by the current arun . Also, Kurokawa
et al. kuro  have shown the oscillations of the free layer magnetization in
the components along the perpendicular directions of the pinned layer
polarization with frequencies upto 576 GHz in the presence of bilinear and
biquadratic interlayer exchange couplings in STNOs, and also with the free
layer having low transition temperature for the saturation magnetization. In
their investigation they have shown that the biquadratic coupling is essential
for the high frequency kuro .

In this connection, our present report provides a detailed study on
$\mathrm{Co\left|RuFe\right|Co}$ STNO with bilinear interlayer exchange
coupling alone between the free and pinned ferromagnetic layers and show the
existence of oscillations of the free layer magnetization in the components
along the pinned layer polarization with frequencies above 300 GHz with the
free layer having high transition temperature. This unaccompanied role of the
bilinear interlayer exchange coupling has been thoroughly researched since it
has been used in many spintronics devices ishi , and multilayer magnetic thin
films. Depending on the interfacial exchange coupling, both negative and
positive exchange couplings have been seen in ferromagnetic/ferrimagnetic
transition of metal and rare-earth alloy multilayer thin films alt ; mish  and
the role of the bilinear coupling co-efficient are experimentally studied in
Ref. nunn . However, numerical and analytical studies on the bilinear coupling
in STNO without an external magnetic field that leads to magnetization
oscillations have not been thoroughly studied in the literature gusa .

The paper is organized as follows. First, we formulate the model and the
governing LLGS equation of motion and effective magnetic field for the present
study in Sec. II. The positive and negative bilinear coupling dynamics and
expression for minimum current for oscillations are presented in Sec. III and
IV, respectively. Section V is devoted to the conclusion of the present work.

##  II Model

Figure 1: Schematic illustration of Co/ReFe/Co trilayer

The schematic picture of an STNO considered for our study, which consists of a
free layer, a spacer layer and a pinned layer, is shown in Fig.1. The
magnetization of the free layer is denoted as ${\bf M}=M_{s}{\bf m}$, where
$M_{s}$ is the saturation of the magnetization. While the magnitude of the
magnetization is fixed, its direction can change over time. The magnetization
of the pinned layer ${\bf P}=M_{s}{\bf p}$ is fixed for both magnitude and
direction. Here ${\bf m}$ and ${\bf p}$ are the unit vectors along ${\bf M}$
and ${\bf P}$, respectively. As shown in Fig.1, the positive and negative
currents correspond to the flow of electrons from the free layer to pinned
layer and vice versa, respectively. The free and pinned layers are considered
to be made up of Co. The spacer layer is a nonmagnetic conductive layer,
constituting an alloy of Ru and Fe. The magnetization dynamics described by
the LLGS equation that governs the motion of the unit vector ${\bf m}$ is
given as

$\displaystyle\frac{d{\bf m}}{dt}=$ $\displaystyle-\gamma{\bf m}\times{\bf
H}_{eff}+\alpha{\bf m}\times\frac{d{\bf m}}{dt}+\gamma H_{S}~{}{\bf
m}\times({\bf m}\times{\bf p}).$ (1)

Here, $\gamma$ and $\alpha$ are the gyromagnetic ratio and damping parameter,
respectively. The spin-torque strength is

$\displaystyle H_{S}=\frac{\hbar\eta I}{2eM_{s}V(1+\lambda{\bf m}\cdot{\bf
p}))},$ (2)

where $\hbar$ is the reduced Planck’s constant ($\hbar(=h/2\pi)$), $I$ is the
current, $e$ is the electron charge, and $V$ is the volume of the free layer,
$\eta$ and $\lambda$ are the dimensionless parameters determining magnitude
and angular dependence of the spin-transfer torque.

The effective magnetic field $H_{eff}$ is given by

$\displaystyle{\bf H}_{eff}={\bf H}_{ani}+{\bf H}_{dem}+{\bf H}_{bil},$ (3)

where ${\bf H}_{ani}$ and ${\bf H}_{dem}$ is the anisotropy and the
demagnetization field, respectively. The effective field also consists of a
bilinear coupling interaction ${\bf H}_{bil}$ of interlayer exchange coupling
between the free and reference layers, the details of which are given below.
Specifically, the various interactions in (3) are given by

${\bf H}_{ani}=H_{k}m_{z}~{}{\bf e}_{z},$ (4a) ${\bf H}_{dem}=-4\pi
M_{s}m_{z}~{}{\bf e}_{z},$ (4b) ${\bf H}_{bil}=-\frac{J}{dM_{s}}~{}{\bf
e}_{x}.$ (4c)

Consequently, we have

$\displaystyle{\bf H}_{eff}=(H_{k}-4\pi M_{s})m_{z}~{}{\bf
e}_{z}-\frac{J}{dM_{s}}~{}{\bf e}_{x}.$ (5)

Here ${\bf e}_{x}$, ${\bf e}_{y}$ and ${\bf e}_{z}$ are the respective unit
vectors along the positive $x$, $y$ and $z$ directions. $H_{k}$ is the
magneto-crystalline anisotropy constant, $J$ is the coefficient of the
bilinear coupling, $M_{s}$ is the saturation magnetization and $d$ is the
thickness of the free layer. The energy density of the free layer responsible
for the effective field ${\bf H}_{eff}=-\partial E/\partial(M_{s}{\bf m})$ is
given by

$\displaystyle E=$ $\displaystyle\frac{J}{d}~{}{\bf m}.{\bf
p}-\frac{M_{s}}{2}[H_{k}-4\pi M_{s}]({\bf m}.{\bf e}_{z})^{2}.$ (6)

The pinned layer is considered to be polarized along positive $x$-direction,
i.e. ${\bf p}={\bf e}_{x}$. The material parameters are adapted as $M_{s}$ =
1210 emu/c.c., $H_{k}$ = 3471 Oe, $\eta$ = 0.54, $\lambda$ = $\eta^{2}$, $d$ =
2 nm, $A$ = $\pi\times$60$\times$60 nm2, $V$ = $Ad$, $\alpha$ = 0.005 and
$\gamma$ = 17.64 Mrad/(Oe s). Since $H_{k}<4\pi M_{s}$, the system exhibits
easy-plane anisotropy for $xy$-plane or hard axis anisotropy for $z$-axis due
to the resultant demagnetization field $-(4\pi M_{s}-H_{k})m_{z}~{}{\bf
e}_{z}$. It means that the magnetization is always pulled towards the $xy$
plane whenever it moves away from the plane with the strength directly
proportional to $m_{z}$. Therefore, before applying any current, to minimize
the energy (Eq.(6)), the magnetization of the free layer settles at (-1,0,0)
for positive bilinear coupling ($J>0$) or (1,0,0) for negative bilinear
coupling ($J<0$). This implies that the system exhibits antiferromagnetic
coupling for the positive bilinear coupling and ferromagnetic coupling for the
negative bilinear coupling between the free and pinned layers mck . It has
been shown that the magnitude and sign of the bilinear coupling coefficient
can be experimentally tuned by changing the concentration of Fe in the spacer
layer made by Ru100-xFex alloy nunn  since the oscillations are observed when
$I<0$ for the positive bilinear coupling and $I>0$ for the negative bilinear
coupling, and both the cases of the bilinear couplings are investigated
separately in the following sections.

##  III Dynamics for the positive bilinear coupling

In the absence of current the equilibrium state of the unit magnetization
vector ${\bf m}$ for the positive bilinear coupling is ${\bf S}_{1}$ =
(-1,0,0) since the field due to the interaction ${\bf H}_{bil}$ acts along the
negative $x$-direction. This is confirmed in Figs.2(a) and 2(b), where the
time evolution of $m_{x}$ and $m_{y}$ are plotted for $J$ = 0.756 mJ/m2 and
0.352 mJ/m2, respectively, for different initial conditions. In both these
figures 2(a) and 2(b), we can observe that the magnetization finally reaches
the state ${\bf S}_{1}$. These numerical results coincide well with the
experimental results obtained by Nunn $et~{}al$ nunn , where the same system
exhibits antiparallel configuration between the magnetizations of the free and
pinned layers for $J$ = 0.756 mJ/m2 and 0.352 mJ/m2 corresponding to Ru32Fe68.

Figure 2: (a) Time evolution of $m_{x}$ and $m_{y}$ for (a) $J$ = 0.756 mJ/m2
and (b) $J$ = 0.352 mJ/m2.

When the current is applied, depending upon the magnitude of the current, the
system exhibits three different dynamics for ${\bf m}$. (i) When
$|I|<|I_{min}|$, the unit magnetization vector ${\bf m}$ stays in the state
${\bf S}_{1}$ where it was existing already. (ii) When
$|I_{min}|<|I|<|I_{max}|$, the vector ${\bf m}$ exhibits continuous
precession. (iii) When $|I|>|I_{max}|$ the vector ${\bf m}$ moves away from
(-1,0,0) and settles into the state S2 (near (0,0,$\pm$1)) for small $J$
($<$2.8 mJ/m2) or settles into the state ${\bf S_{3}}$=(1,0,0) for large $J$
($>$2.8 mJ/m2). Hence the states ${\bf S}_{1}$, ${\bf S}_{2}$ and ${\bf
S}_{3}$ are associated with the currents when $|I|<|I_{min}|$, $|I|>|I_{max}|$
for $J$ ($>$2.8 mJ/m2) and $|I|>|I_{max}|$ for $J$ ($<$2.8 mJ/m2),
respectively. The critical value of the positive bilinear coupling strength
$J_{c}$ = 2.8 mJ/m2 is derived in Eq.(10). Here, $I_{min}$ and $I_{max}$ are
the minimum and maximum currents, respectively, between which oscillations can
be exhibited.

Figure 3: (a) Trajectory of ${\bf m}$ during $t$ = 299-300 ns when $I$ = -0.5
mA (red), -1 mA (blue), -1.5 mA (brown), -2.5 mA (̱black), -3.25 mA (magenta)
and -4 mA (black point) for $J$ = 0.4 mJ/m2. (b) Trajectory of ${\bf m}$
during $t$ = 299-300 ns when $I$ = -2 mA (red), -2.1 mA (blue), -2.2 mA
(black), -2.3 mA (magenta), -2.35 (orange) and -3 mA (black point) for $J$ = 7
mJ/m2. (c) Time evolution of $m_{x}$ when $J$ = 0.4 mJ/m2 and $I$ = -1.5 mA.
(d) Time evolution of $m_{x}$ when $J$ = 7 mJ/m2 and $I$ = -2.3 mA.

To confirm the precession of ${\bf m}$, oscillations of $m_{x}$ and tunability
of the frequency by current, Eq.(1) is numerically solved by adaptive step
size Runge-Kutta-4 method. The initial condition of ${\bf m}$, for the
numerical simulation, is randomly chosen near the state ${\bf S}_{1}$. When a
negative current is applied with the magnitude $|I_{min}|<|I|<|I_{max}|$, the
magnetization which was in the ${\bf S}_{1}$ state moves away from it due to
the spin-transfer torque. This is due to the fact that the incoming electrons
in the free layer, which were spin polarized along the positive $x$-direction,
always move the magnetization to align with the positive $x$-direction. Once
the magnetization moves away from the state ${\bf S}_{1}$ by STT, continuous
precession is achieved due to the balance between the damping (due to the
effective field) and the STT. The trajectories of ${\bf m}$ (after transition
and between $t$ = 299 ns and $t$ = 300 ns) in continuous precession at
different currents for a low value of $J$(= 0.4 mJ/m2) and the time evolution
of $m_{x}$ corresponding to $J$ = 0.4 mJ/m2 and $I$ = -1.5 mA are plotted in
Figs.3(a) and (c), respectively. Similarly, the trajectories of ${\bf m}$ in
the same duration for a high value of $J$(= 7.0 mJ/m2) and the time evolution
of $m_{x}$ corresponding to $J$ = 7 mJ/m2 and $I$ = -2.3 mA are plotted in
Figs.3(b) and (d), respectively. We can observe from Fig.3(a) that the
trajectory corresponding to the current $I$ = -0.5 mA (red) exhibits in-plane
precession around the $x$-axis due to the field from positive bilinear
coupling. The direction of the precession is clockwise as seen from the
positive $x$-axis. When the strength of the current is increased further to
$I$ = -1 mA (blue), the trajectory of the magnetization slightly transforms as
shown in Fig.3(a). It seems that the trajectory has been folded along the
negative $x$-axis. The magnetization gets close to the positive $x$-axis when
it reaches the $xy$-plane. This is due to the fact that the resultant
demagnetization field becomes weaker when the magnetization gets closer to the
$xy$-plane. Therefore the STT, which always moves the ${\bf m}$ towards the
positive $x$-axis, becomes stronger and moves the magnetization towards the
positive $x$-axis as much as possible. Once the magnetization crosses the
$xy$-plane, the magnetization moves away from the positive $x$-axis. This is
due to the fact that the resultant demagnetization field rotates the
magnetization from negative to positive $y$-axis in the northern hemisphere
and from positive to negative $y$-axis in the southern hemisphere. When the
current is further increased to -1.5 mA (brown), the magnetization shows a
transition from the in-plane precession to out-of-plane precession around the
$z$-axis as shown in the Fig.3(a). This is because an increase of curent
increases the magnitude of the STT and consequently the projection of ${\bf
m}$ in the $xy$-plane crosses the positive $x$-axis before the ${\bf m}$
reaches the $xy$-plane. Therefore the bilinear exchange coupling field and the
resultant demagnetization field along with the STT precess the magnetization
within the northern hemisphere continuously. The out-of-plane precessions may
symmetrically take place in the southern or northern hemisphere. Further
increment in the current to -2.5 mA (black) and -3.25 mA (magenta) makes the
concentric trajectories of ${\bf m}$ around the equilibrium magnetization
state where the ${\bf m}$ settles when $|I|>|I_{max}|$, with $I_{max}$ = - 3.4
mA for $J$ = 0.4 mJ/m2. The black point in Fig.3(a) corresponds to the
equilibrium state at which the unit vector ${\bf m}$ settles for $I$ = -4 mA
when $J$ = 0.4 mJ/m2. This equilibrium state can be identified as follows: The
LLGS equation given by Eq.(1) is transformed into spherical polar coordinates
using the transformation equations
$m_{x}=\sin\theta\cos\phi,~{}m_{y}=\sin\theta\sin\phi,~{}m_{z}=\cos\theta$ as

$\displaystyle\frac{d\theta}{dt}~{}=~{}$
$\displaystyle\frac{\gamma}{1+\alpha^{2}}\Bigg{\\{}-\frac{J}{dM_{s}}(\alpha\cos\theta\cos\phi-\sin\phi)$
$\displaystyle-\alpha(H_{k}-4\pi M_{s})\sin\theta\cos\theta$ $\displaystyle-
H_{S0}\frac{(\alpha\sin\phi+\cos\theta\cos\phi)}{(1+\lambda\sin\theta\cos\phi)}\Bigg{\\}}=P(\theta,\phi),$
(7) $\displaystyle~{}\frac{d\phi}{dt}~{}=~{}$
$\displaystyle\frac{\gamma\csc\theta}{1+\alpha^{2}}\Bigg{\\{}\frac{J}{dM_{s}}(\cos\theta\cos\phi+\alpha\sin\phi)$
$\displaystyle+(H_{k}-4\pi M_{s})\sin\theta\cos\theta$
$\displaystyle+H_{S0}\frac{(\sin\phi-\alpha\cos\theta\cos\phi)}{(1+\lambda\sin\theta\cos\phi)}\Bigg{\\}}=Q(\theta,\phi).$
(8)

Here, $\theta$ and $\phi$ are the polar and azimuthal angles, respectively,
$H_{S0}=\hbar\eta I/2eM_{s}V$. The equilibrium state is obtained from the
equations $P(\theta^{*},\phi^{*})=0$ and $Q(\theta^{*},\phi^{*})=0$, where
$\phi^{*}$ is numerically observed as $\phi^{*}\approx 0$. This leads us to
derive the relation

$\displaystyle\sin\theta^{*}=J/(dM_{s}(4\pi M_{s}-H_{k})).$ (9)

Therefore, the equilibrium state ${\bf S}_{2}$ for ${\bf m}$ when
$|I|>|I_{max}|$ is given by ${\bf
S}_{2}\approx(\sin\theta^{*},0,\pm\cos\theta^{*})$, where $\sin\theta^{*}$ is
as given above.

However, when the magnitude of the current is increased much further than
$|I_{max}|$, the equilibrium state will slightly move away from the state
${\bf S}_{2}$ and if the magnitude of the current is extremely large
($|I|>>|I_{max}|$), i.e above $\sim$100 mA, then the magnetization will settle
in the state ${\bf S}_{3}$ = (1,0,0).

Figure 4: (a) Frequency tunability by current for different values of bilinear
coupling $J$ (given in mJ/m2) from 0.1 mJ/m2 to 12 mJ/m2. (b) Frequency
against bilinear coupling for different values of current $I$.

From Eq.(9), we can understand that the value of $\theta^{*}$ becomes $\pi/2$
when $J=dM_{s}(4\pi M_{s}-H_{k})$. It means that the equilibrium state ${\bf
S}_{2}$ of the magnetization moves towards the state ${\bf S}_{3}$ = (1,0,0)
as the strength of the positive bilinear coupling $J$ increases and reaches
(1,0,0) when $J\rightarrow J_{c}$, where

$\displaystyle J_{c}=dM_{s}(4\pi M_{s}-H_{k})=2.8~{}{\rm mJ/m^{2}}.$ (10)

Similarly, the magnetization precession for the high strength of bilinear
coupling ($J$ = 7.0 mJ/m2) is also investigated by plotting the trajectories
for the currents $I$ = -2 mA (red), -2.1 mA (blue), -2.2 mA (black), -2.3 mA
(magenta), -2.35 mA (orange) and -3 mA (black point) in Fig.3(b). Unlike the
case of low bilinear coupling as shown in Fig.3(a), there is no transition
from in-plane to out-of-plane precession while increasing the magnitude of the
current and the magnetization exhibits in-plane precession only around the
$x$-axis. This can be reasoned as follows: When the strength of the bilinear
coupling field is strong due to large $J(>0)$, the STT and the resultant
demagnetization field are dominated by this bilinear coupling field.
Therefore, the rotations due to the resultant demagnetization field and the
approach of the magnetization towards the positive $x$-axis due to the STT are
not exhibited. When the current is increased further, the trajectory moves
from the negative to positive $x$-axis and settles into the equilibrium state
${\bf S}_{3}$ when $|I|>|I_{max}|$, where $I_{max}$ = -2.35 mA for $J$ = 7.0
mJ/m2. The equilibrium state for the current -3 mA is shown by the black point
in the Fig.3(b).

To confirm the oscillations the time evolutions of the component $m_{x}$ are
plotted in Fig.3(c) for $J$ = 0.4 mJ/m2, $I$ = -1.5 mA and in Fig.3(d) for $J$
= 7.0 mJ/m2, $I$ = -2.3 mA. The frequencies of the oscillations are 16 GHz and
163 GHz, respectively.

The frequencies of the oscillations, of $m_{x}$ are plotted against the
current for different values of bilinear coupling strengths (given in mJ/m2)
from 0.1 mJ/m2 to 12 mJ/m2 in Fig.4(a) and against bilinear coupling for
different values of current in Fig.4(b). From Fig.4(a), we can understand that
when the bilinear coupling coefficient is low, the frequency decreases up to
some critical current $I_{c}$ and then increases. This change in the frequency
from decrement to increment is attributed to the transition of magnetization
precession from the in-plane to out-of-plane as discussed earlier with
reference to Fig.3(a). In Fig.4(a), the existence of $I_{min}$ and $I_{max}$
is evident, and the range of current for the oscillations
($|I_{max}|-|I_{min}|$) confirms the wide frequency tunability by the current.
The magnitude of $I_{c}$ slightly decreases with the increase of $J$. Also, we
can observe that when $J$ is large ($\geq$2.9 mJ/m2) the frequency decreases
with the increase in the magnitude of the current up to $I_{max}$ and the
$I_{c}$ does not exist. This is due to the nonexistence of out-of-plane
precession, as shown in Fig.3(b). From Fig.4(a) it is observed that the
tunability range ($|I_{max}|-|I_{min}|$) decreases and increases with $J$ when
the strength of $J$ is small and large, respectively. At a given current, the
frequency increases with the magnitude of bilinear coupling. Also, it is
confirmed that the frequency can be enhanced up to 300 GHz for $J$ = 12.0
mJ/m2 and even above when $J$ is increased further.

Similarly, the frequency is plotted against $J$ for different values of the
current in Fig.4(b). Due to the nonexistence of out-of-plane precession at
large strengths of $J$, the discontinuity appears in the frequency while
increasing the value of $J$ as shown in Fig.4(b). From Fig.4(b) we can observe
that the frequency almost linearly enhances with $J$. The frequency range is
around 30 GHz and 300 GHz when the values of $J$ are small and large,
respectively. The enlargement of frequency and switching time can be
essentially attributed to the large value of the bilinear coupling strength
$J$, which causes the system to behave more like a layered antiferromagnet vol
; vigo ; roy ; fira ; khym . The large value of $J$ in our system is possibly
due to Nunn et al.’s recently proposed RuFe spacer layer nunn . The current
density corresponding to the frequency 299.6 GHz when $I$ = -3.35 mA can be
obtained as 2.96$\times 10^{7}$ A/cm2 for the cross sectional area
$A=\pi\times 60\times 60$ nm2. Also, it is visible that the magnitude of the
current can increase the range of $J$ for which the oscillations are possible.

Figure 5: Frequency dependence on current and different ranges of bilinear
coupling coefficient. The open circles are the minimum critical current
$I_{max}$, for the onset of the oscillations, obtained from Eq.(14). S1 =
(-1,0,0), S2 = $(\sin\theta*,0,\pm\cos\theta*)$ and S3 = (1,0,0) are the
equilibrium states.

Figs.5(a) and (b) summarize the dependence of the frequency on current and $J$
while $J$ is below and above 2.3 mJ/m2, respectively. The white color region
is nonoscillatory region. From Figs.5(a) & (b), we can see that the magnitude
of the current above which the oscillations occur ($|I_{min}|$) linearly
increases with $J$. The value $I_{min}$ for $J>$0 can be derived as follows:

The nature of the stability of an equilibrium state which is represented by
polar coordinates can be identified from the following Jacobian matrix by
using Eqs.(7) and (8)

$\displaystyle\mathcal{J}=\begin{pmatrix}{\left.\frac{dP}{d\theta}\right|}_{(\theta^{*},\phi^{*})}&{\left.\frac{dP}{d\phi}\right|}_{(\theta^{*},\phi^{*})}\\\
{\left.\frac{dQ}{d\theta}\right|}_{(\theta^{*},\phi^{*})}&{\left.\frac{dQ}{d\phi}\right|}_{(\theta^{*},\phi^{*})}\end{pmatrix}.$
(11)

The equilibrium state $(\theta^{*},\phi^{*})$ will be stable only when the
system is dissipative about it. It will be dissipative if and only if the
trace of the matrix $\mathcal{J}$ becomes negative,

$\displaystyle{\rm Tr}(\mathcal{J})<0.$ (12)

We knew that when $|I|<|I_{c}^{min}|$ and $J>0$ the magnetization settles at
${\bf S}_{1}$, i.e, ($\pi/2,\pi$) in polar coordinates. Therefore specific set
of values $(\theta^{*},\phi^{*})=(\pi/2,\pi)$ satisfies Eq.(12). The trace of
the matrix corresponding to ($\pi/2,\pi$) is given by

$\displaystyle{{\rm
Tr}(\mathcal{J})|}_{(\theta^{*},\phi^{*})}=\frac{\gamma}{1+\alpha^{2}}\left[-\frac{2J\alpha}{dM_{s}}+(H_{k}-4\pi
M_{s})\alpha-\frac{2H_{S0}}{1+\lambda}\right].$ (13)

The minimum critical current $I_{min}$ (for $J>0$), below which the ${\bf
S}_{1}$ is stable can be derived from Eqs.(12) and (13) as

$\displaystyle I_{min}=\frac{eA\alpha(\lambda-1)}{d\hbar\eta}\left[2J+(4\pi
M_{s}-H_{k})dM_{s}\right]$ (14)

and it has been plotted as open circles in Figs.5(a) and (b), which matches
well with the numerical results and confirms the validity of the numerical
results. From Fig.5(a) and (b) we can observe that value of $I_{max}$
decreases with $J$ at lower strengths of $J$ and increases (almost linearly)
with $J$ at higher strengths of it. Fig.5(b) evidences that the range of
current which exhibits oscillations increases with $J$ while $J$ is large. In
the case of positive current, the STT always moves the magnetization to be
aligned with the negative x-direction. Therefore the positive current does not
move the magnetization from the state (-1,0,0), where it existed already
before the application of the current, and therefore no precession is
exhibited.

Figure 6:  Magnetization switching for different values of (a) $J$ when $I$ =
-2.5 mA and (b) $I$ when $J$ = 4.5 mJ/m2. Time to switch from (-1,0,0) to
(1,0,0) with respect to (c) current and (d) bilinear coupling.

We can observe in Figs.5(a) and 5(b) that the magnetization settles into the
equilibrium states S2 and S3, respectively when $I>I_{max}$. It indicates a
transition from S2 to S3 while increasing the strength of the positive
bilinear coupling. As discussed in Eq.(10), the transition occurs at $J$ = 2.8
${\rm mJ/m^{2}}$.

From Fig.5(b), we can observe that when the magnitude of the current is above
the magnitude of $I_{max}$, the magnetization will settle into the state ${\bf
S}_{3}$ from ${\bf S}_{1}$ for the positive bilinear coupling. This indicates
the existence of current-induced magnetization switching from the negative to
positive $x$-direction. The corresponding switchings of $m_{x}$ from -1 to +1
for different values of bilinear coupling when $I$ = -2.5 mA and current when
$J$ = 4.5 mJ/m2 are plotted in Figs.6(a) and (b), respectively. From Fig.6(a)
we can observe that the switching times for $J$ = 3.0, 4.5 and 6.0 mJ/m2 are
4.42, 6.01 and 9.42 ns, respectively. Hence, the switching time increases with
the magnitude of the positive bilinear coupling. On the other hand, from
Fig.6(b) we can understand that the switching times for the currents $I$ =
-2.0, -2.5 and -3.0 mA are 9.88, 6.01 and 3.892 ns, respectively. This implies
that the switching times reduce with the increase of the magnitude of the
current. The variation of the switching time against current and the strength
of the bilinear coupling for different values of $J$ and $I$ are plotted in
Figs.6(c) and (d), respectively. Figs.6(c) and (d) confirm the decrement and
increment of the switching time with the increase in the magnitude of current
and positive bilinear coupling, respectively. Since the field due to the
positive bilinear coupling acts along the negative $x$-direction, the
enhancement in the magnitude of the negative current can quickly reverse the
magnetization from negative to positive $x$-direction as shown in Fig.6(c).
Similarly, when the strength of the positive bilinear coupling increases, its
corresponding field along the negative $x$-direction increases, and
consequently the magnetization takes much time to reverse from the negative to
positive $x$-direction by the application of negative current as confirmed in
Fig.6(d). The above current-induced magnetization switching has spin torque
magnetic random access memory applications and is much more efficient than the
field-induced switching. The field-free switching may help produce magnetic
memory devices with low power consumption and greater device density Loca ;
Ral .

As observed from Figs.5, when the current $I$ is kept constant and the
strength of the positive bilinear coupling $J$ is increased, the magnetization
reaches the equilibrium state ${\bf S}_{2}$ via out-of-plane precession (see
Fig.3(a)). When $J$ is increased further, the equilibrium state of the
magnetization ${\bf S}_{2}$ becomes (1,0,0) as $J\rightarrow J_{c}$ (see Eq.10
in the revised manuscript). After the magnetization reaches the state ${\bf
S}_{3}$ it continues to settle there without showing any oscillations until
the further increase in $J$ is strong enough to move away the magnetization
from the state ${\bf S}_{3}$ against the STT due to the incoming spin
polarized electrons. As observed in Fig.4(b) and Figs.5, the gap between the
offset of oscillations of ${\bf m}$ when reaching ${\bf S}_{2}$ and the onset
of oscillations when emanating from ${\bf S}_{3}$ increases with the magnitude
of the current. This is due to the fact that the strength of the STT which
tends to keep the magnetization along the positive $x$-direction increases
with the magnitude of current and consequently the strength of the bilinear
coupling is required to be high enough to regain the oscillations from the
equilibrium state ${\bf S}_{3}$.

##  IV Dynamics for the negative bilinear coupling

Figure 7: (a-b) Frequency tunability by current for different values of
negative bilinear coupling $J$. (c) The magnetization trajectory when $I$ = 1
mA (red), 2 mA (blue), 10 mA (brown), 20 mA (black), 36 mA (magenta) and 37 mA
(black point) for $J$ = -0.1 mJ/m2. (d) The frequency variation against
negative bilinear coupling for different values of current.

In the presence of negative bilinear coupling the magnetization will initially
be oriented at ${\bf S}_{3}$ since the field due to the negative bilinear
coupling ${\bf H}_{bil}$ acts along the positive $x$-direction. The
magnetization continues to be settled at ${\bf S}_{3}$ until the current $I$
is increased to $I_{min}$. The STT, due to the positive current, will always
move the magnetization to be aligned with the negative $x$-direction. When
$I>I_{min}$, the magnetization is moved away from ${\bf S}_{3}$, and the
system shows continuous precession for the vector ${\bf m}$. The frequency of
the oscillations of $m_{x}$ is plotted against low values of current in
Fig.7(a) and high values of current in Fig.7(b) for different values of the
negative bilinear coupling (given in mJ/m2). From Fig.7(a), we can understand
that similar to the case of the positive bilinear coupling, the frequency
decreases with current up to a critical value $I_{c}$ and then increases with
current. Similar to the previous case, this increment in frequency after
decrement is attributed to the transition from in-plane to out-of-plane
precession. This is verified by plotting the trajectories of the vector ${\bf
m}$ corresponding to $I$ = 1 mA (red) and 2 mA (blue) for $J$ = -0.1 mJ/m2 in
Fig.7(c). Since the field, due to negative bilinear coupling, acts along the
positive $x$-direction, the magnetisation trajectory corresponding to $I$ = 1
mA (red) has been folded along the positive $x$-axis and exhibits in-plane
precession.

When the current increases to 2 mA (blue), the magnetization transforms from
in-plane precession to out-of-plane precession in the northern hemisphere.
However, the out-of-plane precession may also be symmetrically placed in the
southern hemisphere. The explanation behind this transition is similar to
those discussed in the case of positive bilinear coupling. The out-of-plane
precessions corresponding to the currents $I$ = 10 mA (brown), 20 mA (black)
and 36 mA (magenta) for $J$ = -0.1 mJ/m2 also are plotted in Fig.7(c). From
Fig.7(a), we can understand that when the strength of the negative bilinear
coupling is relatively high, the frequency shows only an increment with the
current. This is because at higher values of negative bilinear coupling, the
unit magnetization vector ${\bf m}$ exhibits out-of-plane precession instead
of exhibiting any transition from in-plane to out-of-plane precession. In
Fig.7(b), the frequency is plotted up to large values of current for different
values of $J$. The frequency increases with current and reaches its maximum.
For small values of $J$, the frequency increases to its maximum and then
decreases. Fig.7(b) shows that there is a maximum current $I_{max}$ above
which oscillations are not possible. For the currents above $I_{max}$, the
magnetization settles into ${\bf S}_{1}$ without showing any precession.

In Fig.7(b) we can observe the discontinuities for frequencies near $I_{max}$
upto $J\approx$ -0.4 mJ/m2, where the system exhibits multistability i.e the
magnetization may precess continuously or settle at ${\bf S}_{1}$. It is
confirmed in Fig.7(c) by precession for $I$ = 36 mA (magenta) and equilibrium
state ${\bf S}_{1}$ for $I$ = 37 mA (black point). In Fig.7(b) it is observed
that the discontinuities in the frequencies have disappeared above $J$ = -0.4
mJ/m2. This is because the magnetization does not settle at ${\bf S}_{1}$
below $I_{max}$. The magnetization exhibits three different nature of
equilibrium states for $|J|>\sim 0.4$ and $I>I_{max}$. When the current is
increased near above $I_{max}$, the magnetization settles near poles at ${\bf
S}_{2}$. When $I$ is increased further the unit vector ${\bf m}$ settles into
${\bf S}_{2}$ or ${\bf S}_{1}$. If the current is increased further to
extremely large values, the magnetization settles into ${\bf S}_{1}$. The
range of the current in which the oscillations are possible
($I_{max}-I_{min}$) also increases (decreases) with $|J|$ when $|J|$ is small
(large).

From Figs.7(a) and (b), it is observed that the frequency can be reached
around 30 GHz by increasing the current and the magnitude of the negative
bilinear coupling. In Fig.7(d), the frequency is plotted against the negative
bilinear coupling for different values of the currents. It seems that the
frequency increases almost linearly with the increase in the magnitude of
negative bilinear coupling coefficient. Also, at a given $J$, the frequency
increases with the magnitude of the current.

Figure 8: Dependence of the frequency on $J$ and $I$. (a) $I<$ 90 mA and (b)
$I<$3 mA.

The dependence of the frequency on the negative bilinear coupling and current
is plotted for the large values of current in Fig.8(a) and small values of
current in Fig.8(b). The white background corresponds to the non-oscillatory
region. From Fig.8(a) we can observe that the value of $I_{max}$ increases up
to -0.33 mJ/m2 and then decreases abruptly. From the bright green and red
regions in Fig.8(a) we can understand that the frequency can be maintained
constant while increasing the current at fixed $J$. Also, it is clearly
visible that the tunability range of the frequency by current drastically
reduces after $\sim$-0.3 mJ/m2. This is different from the case of positive
bilinear coupling where the ocillatory region ($|I_{max}|-|I_{min}|$) can be
expanded with the increase of $J$. For currents above $I_{max}$, three
different regions are identified for ${\bf m}$ as shown in Fig.8(a). The three
different regions for equilibrium states ${\bf S}_{1}$, ${\bf S}_{2}$ and
${\bf S}_{1}$/${\bf S}_{2}$ for the current above $I_{max}$ are indicated in
Fig.8(a). To see the minute variation of frequency in the low current region,
Fig.8(b) is plotted for currents upto 3 mA. Fig.8(b) confirms the decrement
and increment in frequency with current when $|J|<1$ mJ/m2. Also, the
frequency at a given current increases with the strength of the negative
bilinear coupling.

The minimum current $I_{min}$ for $J<$0 is similarly derived as in the
previous case for positive bilinear coupling. When $I<I^{min}$ and $J<0$, the
state ${\bf S}_{3}$ becomes stable and the magnetization settles into ${\bf
S}_{3}$, corresponding to $(\pi/2,0)$ in polar coordinates. The trace of the
matrix $J$ corresponding to the state $(\pi/2,0)$ is derived as

$\displaystyle{{\rm
Tr}(\mathcal{J})|}_{(\pi/2,0)}=\frac{\gamma}{1+\alpha^{2}}\left[\frac{2J\alpha}{dM_{s}}+(H_{k}-4\pi
M_{s})\alpha+\frac{2H_{S0}}{1-\lambda}\right].$ (15)

From the condition (12) and Eq.(15), we can derive the minimum current (for
$J<0$) below which the equilibrium state ${\bf S}_{3}$ is stable as

$\displaystyle
I_{min}=-\frac{eA\alpha(1+\lambda)}{d\hbar\eta}\left[2J+(H_{k}-4\pi
M_{s})dM_{s}\right].$ (16)

Eq.(16) is plotted in Fig.8(b) as open circles and matches well with the
numerical results. This confirms the validity of the numerical results.

If the current is negative, the STT always moves the magnetization towards the
positive x-direction. Therefore the magnetization does not move from the state
${\bf S}_{3}$, where it was already existing before applying the current, by
the negative current, and no precession is exhibited.

Figure 9:  Magnetization switching (negative bilinear coupling) for different
values of (a) $J$ when $I$ = 33.5 mA and (b) $I$ when $J$ = -0.05 mJ/m2. Time
to switch from ${\bf S}_{1}$ to ${\bf S}_{3}$ with respect to (c) current and
(d) bilinear coupling.

Similar to the case of positive bilinear coupling, magnetization switching can
also be identified for negative bilinear coupling. As discussed in Fig.8(a)
when a current corresponding to the region of equilibrim state ${\bf S}_{1}$
is applied the magnetization will switch from ${\bf S}_{3}$ to ${\bf S}_{1}$.
In Figs.9(a) and (b) the component $m_{x}$ is plotted to confirm the switching
from positive to negative $x$-direction for different values of $J$ when $I$ =
33.5 mA and for different values of $I$ when $J$ = -0.05 mJ/m2, respectively.
The variation of the switching time against current and the coupling is
plotted in Figs.9(c) and (d), respectively. From Figs.9(a) and (c), we can
understand that similar to the positive bilinear coupling, the switching time
decreases with the increase in the magnitude of the current. Fig.9(d) confirms
that there is no definite relationship between the switching time and the
negative bilinear coupling. The switching time variation against the magnitude
of the coupling is not smooth like in the case of positive bilinear coupling.

##  V Conclusion

In conclusion, we have investigated the dynamics of
$\mathrm{Co\left|RuFe\right|Co}$ STNO using the LLGS equation and identified
high-frequency oscillations in the magnetization of the free layer due to the
presence of bilinear coupling. The obtained orientations of the magnetization
of the free layer with that of the pinned layer in the absence of current
match well with the experimental results. A transition in the precession of
the magnetization from in-plane precession to out-of-plane precession while
increasing the current is observed for both positive and negative bilinear
coupling cases. However, the transition does not occur at higher strengths of
the bilinear coupling. Only an in-plane precession for the positive bilinear
coupling and an out-of-plane precession for the negative bilinear coupling are
exhibited. A wide range of frequency tunability by the current is observed for
both cases of bilinear coupling. While the frequency is enhanced upto 30 GHz
by the negative bilinear coupling, the positive bilinear coupling enhances the
frequency upto and above 300 GHz. This high frequency has been shown for the
oscillations of the magnetization vector (free layer) along the pinned layer
polarization and with the free layer having high transition temperature for
the saturation magnetization. The range of the current in which the frequency
can be tuned increases with the strength of the positive bilinear coupling
corresponding to the in-plane precession.

Oscillations are exhibited for the positive (negative) bilinear coupling when
the current is applied in the negative (positive) direction. Also,
oscillations are possible only when the current is between $I_{min}$ and
$I_{max}$. When $|I|<|I_{max}|$, the magnetization settles into (-1,0,0) for
$J>0$ and (1,0,0) for $J<0$. If the strength of the positive bilinear coupling
is large, then the magnetization settles into (1,0,0) for all the magnitudes
of the current above $|I_{max}|$. On the other hand, if the strength is small,
it settles near poles (${\bf S}_{2}$) when $|I|>|I_{max}|$ or into (1,0,0)
when $|I|>>|I_{max}|$. If the bilinear coupling is negative, there are three
regions corresponding to the equilibrium states ${\bf S}_{2}$, ${\bf S}_{1}$
(or) ${\bf S}_{2}$ and ${\bf S}_{1}$ above $I_{max}$ depending upon the values
of $I$ and $J$. The magnetization switching induced by the current alone is
identified for both of the bilinear couplings. It is observed that the
switching time reduces with the increase in the magnitude of the current for
both cases of the bilinear coupling.

We have also analyzed the expressions for the minimum currents to achieve the
oscillations for both the positive and negative bilinear couplings. We have
shown that they match well with the numerically obtained results. We have also
proved that the bilinear coupling is sufficient for the high-frequency
oscillations among two interlayer exchange couplings, namely bilinear and
biquadratic couplings. We wish to point out that this study has been carried
out for the temperature $T$ = 0 K. However, the free layer we have considered
is perpendicular magnetic anisotropic one and this is normally robust against
thermal noise Tudu . We believe that our detailed study on bilinear coupling
can be helpful in applications related to microwave generation with high-
frequency enhancement and magnetic memory devices.

## Acknowledgement

The works of V.K.C. and R. G are supported by the DST-SERB-CRG Grant No.
CRG/2020/004353 and they wish to thank DST, New Delhi for computational
facilities under the DST-FIST programme (SR/FST/PS-1/2020/135) to the
Department of Physics. M.L. wishes to thank the Department of Science and
Technology for the award of a DST-SERB National Science Chair under Grant No.
NSC/2020/00029 in which R. Arun is supported by a Research Associateship.

## References

  * (1) J.C. Slonczewski, Phys. Rev. B 39, 6995 (1989); J. Magn. Magn.Mater. 159, L1 (1996).
  * (2) L. Berger, Phys. Rev. B 54, 9353 (1996).
  * (3) J. C. Slonczewski, Phys. Rev. B 71, 024411 (2005).
  * (4) J. A. Katine, F. J. Albert, R. A. Buhrman, E. B. Myers, and D. C. Ralph, Phys. Rev. Lett. 84, 3149 (2000).
  * (5) K. Deng, X. Li, and B. Flebus, Phys. Rev. B 107, L100402 (2023).
  * (6) T. Yamaguchi, S. Tsunegi, K. Nakajima, and T. Taniguchi, Phys. Rev. B 107, 054406 (2023).
  * (7) C. Liu, Y. Kurokawa, N. Hashimoto, T. Tanaka, and H. Yuasa, Scientific Reports 13, 3631 (2023).
  * (8) Y. Imai, S. Tsunegi, K. Nakajima, and T. Taniguchi, Phys. Rev. B 105, 224407 (2022).
  * (9) H. T. Wu, Lei Wang, Tai Min, and X. R. Wang, Phys. Rev. B 104, 014305 (2021).
  * (10) B. Wolba, O. Gomonay, and V. P. Kravchuk, Phys. Rev. B 104, 024407 (2021)
  * (11) A. J. Eklund, M. Dvornik, F. Qejvanaj, S. Jiang, S. Chung, J. Akerman, and B. G. Malm, Phys. Rev. B 103, 214433 (2021).
  * (12) B. Jiang aand W. Zhang, Journal of Magnetism and Magnetic Materials 490, 165470 (2019).
  * (13) M. V. Lubarda, M. Kuteifan, C.-H. Lambert, V. Lomakin, E. E. Fullerton, and S. Mangin, Phys. Rev. B 102, 014405 (2020).
  * (14) J. Grollier, V. Cros, A. Hamzic, J. M. George, H. Jaffres, A. Fert, G. Faini, J. B. Youssef, and H. Legall, Appl. Phys. Lett. 78, 3663 (2001).
  * (15) S. S. P. Parkin, N. More, and K. P. Roche, Phys. Rev. Lett. 64, 2304 (1990).
  * (16) J. Unguris, R. J. Celotta, and D. T. Pierce, J. Appl. Phys. 75, 6437 (1994).
  * (17) T. Nagasima, A. Oguri, and H. Ishii, Journal of Magnetism and Magnetic Materials 177-181, 1008 (1998).
  * (18) A. Moser, A. Berger, D. T. Margulies, and Eric E. Fullerton, Phys. Rev. Lett. 91, 097203 (2003).
  * (19) S. V. Grigoriev, D. Lott, Yu. O. Chetverikov,A. T. D. Grunwald, R. C. C. Ward, and A. Schreyer Phys. Rev. B 82, 195432 (2010).
  * (20) T. McKinnon and E. Girt, App. Phys. Lett 113, 192407 (2018).
  * (21) R. A. Duine, K. J. Lee, S. S. P. Parkin, and M. D. Stiles, Nat. Phys. 14, 217 (2018).
  * (22) Z. R. Nunn, C. Abert, D. Suess, E. Girt, Sci. Adv. 6, eabd8861 (2020).
  * (23) L. Cuchet, B. Rodmacq, S. Auffret, R. C. Sousa, I. L. Prejbeanu, and B.Dieny, Sci. Rep. 6, 21246 (2016).
  * (24) B. Rana and Y. Otani, Commun. Phys. 2, 90 (2019).
  * (25) V. G. Harris, IEEE Trans. Magn. 48, 1075 (2012).
  * (26) J. S. Meena, S. M. Sze, U. Chand, and T.-Y. Tseng, Nanoscale Res. Lett. 9, 526 (2014).
  * (27) E. Grimaldi, A. Dussaux, P. Bortolotti, J. Grollier, G. Pillet, A. Fukushima, H. Kubota, K. Yakushiji, S. Yuasa, and V. Cros, Phys. Rev. B 89, 104404 (2014).
  * (28) S. Tsunegi, H. Kubota, K. Yakushiji, M. Konoto, S. Tamaru, A. Fukushima, H. Arai, H. Imamura, E. Grimaldi, R. Lebru, Appl. Phys. Express 7, 063009 (2014).
  * (29) R. Arun, R. Gopal, V. K. Chandrasekar, and M. Lakshmanan, J. Appl. Phys. 132, 094301 (2022).
  * (30) Y. Kurokawa, K. Yamada, T. Taniguchi, S. Horiike, T. Tanaka, and H. Yuasa, Sci. Rep. 12, 10849 (2022).
  * (31) Y. Ishikuro, M. Kawaguchi, T. Taniguchi and M. Hayashi, Phys. Rev. B 101, 014404 (2020).
  * (32) B. Altuncevahir and A.R. Koymen J. Magn. Magn. Mater. 261, 424 (2003).
  * (33) S. K. Mishra, F. Radu, H. A. Dürr, and W. Eberhardt, Phys. Rev. Lett. 102, 177208 (2009).
  * (34) D. Gusakova, D. Houssameddine, U. Ebels, B. Dieny, L. Buda-Prejbeanu, M. C. Cyrille, and B. Delaet, Phys. Rev. B 79, 104406 (2020).
  * (35) I. Volvach, A.D. Kent, E.E. Fullerton, and V. Lomakin, Phys. Rev. Applied 18, 024071 (2022).
  * (36) H. Vigo-Cotrina and A.P. Guimaraes, J. Magn. Magn. Mater. 497, 166009 (2020).
  * (37) P. E. Roy, R. M. Otxoa, and J. Wunderlich Phys. Rev. B 94, 014439 (2016).
  * (38) I. Firastrau, L. D. Buda-Prejbeanu, B. Dieny, and U. Ebels, J. Appl. Phys. 113, 113908 (2013).
  * (39) R. Khymyn, I. Lisenkov, V. Tiberkevich, B. A. Ivanov, and A. Slavin, Sci. Rep. 7, 43705 (2017).
  * (40) N. Locatelli, V. Cros and J. Grollier, Nat. Mat. 13, 11 (2014).
  * (41) D. C. Ralph and M. D. Stiles, J. Magn. Magn. Mater. 320, 1190 (2008).
  * (42) B. Tudu and A. Tiwari, Vaccum 146, 329 (2017).