Electrical control of magnetization in S/F/S junctions on a 3D topological insulator

TL;DR

This paper explores how the Josephson current in S/F/S junctions on 3D topological insulators can be used to control and stabilize complex magnetization states due to the unique spin-momentum locking of Dirac quasiparticles.

Contribution

It demonstrates the potential for magnetization control via Josephson effects in topological insulator-based junctions, revealing new ways to manipulate magnetic states.

Findings

Josephson current-phase relation strongly depends on magnetization orientation.

Voltage-driven regime can split the ferromagnet's easy-axis.

Unconventional four-fold degenerate ferromagnetic state can be stabilized.

Abstract

Strong dependence of the Josephson energy on the magnetization orientation in Josephson junctions with ferromagnetic interlayers and spin-orbit coupling opens a way to control magnetization by Josephson current or Josephson phase. Here we investigate the perspectives of magnetization control in superconductor/ferromagnet/superconductor (S/F/S) Josephson junctions on the surface of a 3D topological insulator hosting Dirac quasiparticles. Due to the spin-momentum locking of these Dirac quasiparticles a strong dependence of the Josephson current-phase relation on the magnetization orientation is realized. It is demonstrated that this can lead to splitting of the ferromagnet's easy-axis in the voltage driven regime. We show that such a splitting can lead to stabilization of an unconventional four-fold degenerate ferromagnetic state.

Equations

$$\hat{H}=\int d^2{\mathbf{r}}^{\prime }{\hat{\Psi }}^{†}({\mathbf{r}}^{\prime })\hat{H}({\mathbf{r}}^{\prime })\hat{\Psi }({\mathbf{r}}^{\prime }),$$

$$\hat{H}(\mathbf{r})=-i{v}_F(\nabla \times {\mathbf{e}}_z)\hat{\mathbf{\sigma }}+{\mathbf{h}}_{eff}(\mathbf{r})\hat{\mathbf{\sigma }}-\mu ,$$

$$\displaystyle-iv_{F}\bm{n}_{F}\hat{\nabla}\hat{g}=\Bigl{[}\varepsilon\tau_{z}-\hat{\Delta},\hat{g}\Bigr{]}_{\otimes},$$

$$-\frac{1}{2}i{v}_F{\partial }_x{f}^R+{\mathbf{h}}_{eff}{\mathbf{n}}_{\perp }{f}^R=\epsilon f^R-\Delta (x).$$

$$\displaystyle f^{R}_{\pm}=\frac{\Delta e^{\mp i\chi/2}}{\varepsilon}\exp\Bigl{[}{\frac{\mp 2i(\bm{h}_{eff}\bm{n}_{\perp}-\varepsilon)(d/2\pm x)}{v_{x}}}\Bigr{]},$$

$$\displaystyle\tilde{f}^{R}_{\pm}=-\frac{\Delta e^{\mp i\chi/2}}{\varepsilon}\exp\Bigl{[}{\frac{\mp 2i(\bm{h}_{eff}\bm{n}_{\perp}-\varepsilon)(d/2\mp x)}{v_{x}}}\Bigr{]},$$

$$j_x=-\frac{e{N}_F{v}_F}{4}\underset{-\infty }{\overset{\infty }{\int }}d\epsilon \underset{-\pi /2}{\overset{\pi /2}{\int }}\frac{d\varphi }{2\pi }\cos \varphi \times$$

$$\displaystyle\Bigl{[}(g^{R}_{+}\otimes\varphi_{+}-\varphi_{+}\otimes g^{A}_{+})-(g^{R}_{-}\otimes\varphi_{-}-\varphi_{-}\otimes g^{A}_{-})\Bigr{]},$$

$$j_s=j_c\sin (\chi -{\chi }_0),$$

$$j_c=e{v}_F{N}_{F}T\sum _{{\epsilon }_n>0}\underset{-\pi /2}{\overset{\pi /2}{\int }}d\varphi \cos \varphi \frac{{\Delta }^2}{{\epsilon }_n^2}\times$$

$$\exp [-\frac{2{\epsilon }_{n}d}{v_F\cos \varphi }]\cos [\frac{2h_{x}d\tan \varphi }{v_F}],$$

$${\chi }_0=2h_{y}d/v_F,$$

$$j_c=j_b\underset{-\pi /2}{\overset{\pi /2}{\int }}d\varphi \cos \varphi \times$$

$$\exp [-\frac{2\pi Td}{v_F\cos \varphi }]\cos [\frac{2h_{x}d\tan \varphi }{v_F}],$$

$$\frac{\partial \mathbf{M}}{\partial t}=-\gamma \mathbf{M}\times H_{eff}+\frac{\alpha }{M_s}\mathbf{M}\times \frac{\partial \mathbf{M}}{\partial t},$$

$$E_{tot}=E_M+E_J,$$

$$\displaystyle\frac{H_{eff,x}}{H_{F}}=\Gamma\bigg{[}\int_{-\pi/2}^{\pi/2}e^{-\tilde{d}/\cos\phi}\sin\phi\sin\bigg{(}rm_{x}\tan\phi\bigg{)}d\phi\bigg{]}\times$$

$$\displaystyle\bigg{[}1-\cos\bigg{(}\Omega_{J}t-rm_{y}\bigg{)}\bigg{]},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$$

$$\displaystyle\frac{H_{eff,y}}{H_{F}}=\Gamma\bigg{[}\int_{-\pi/2}^{\pi/2}e^{-\tilde{d}/\cos\phi}\cos\phi\cos\bigg{(}rm_{x}\tan\phi\bigg{)}d\phi\bigg{]}$$

$$\displaystyle\sin\bigg{(}\Omega_{J}t-rm_{y}\bigg{)}+m_{y},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$$

$$H_{eff,z}=0,$$

$$H_{eff,x}=A{H}_F{m}_x[1-\cos ({\Omega }_{J}t-r)],$$

$$H_{eff,y}=H_F[1+B\sin ({\Omega }_{J}t-r)],$$

$$A=r\Gamma \underset{-\pi /2}{\overset{\pi /2}{\int }}{e}^{-\tilde{d}/\cos \varphi }\sin \varphi \tan \varphi d\varphi .$$

$$\displaystyle\dot{m}_{x}=\frac{\gamma}{1+\alpha^{2}}\Bigl{[}H_{eff,y}(m_{z}-\alpha m_{x})+\alpha H_{eff,x}\Bigr{]},$$

$$\displaystyle\dot{m}_{z}=\frac{\gamma}{1+\alpha^{2}}\Bigl{[}-H_{eff,y}(m_{x}+\alpha m_{z})+H_{eff,x}\Bigr{]},$$

$$\displaystyle\dot{m}_{x}=\frac{\Omega_{F}}{1+\alpha^{2}}\Bigl{[}m_{z}-\alpha(1-A)m_{x}\Bigr{]},$$

$$\displaystyle\dot{m}_{z}=\frac{\Omega_{F}}{1+\alpha^{2}}\Bigl{[}-(1-A)m_{x}-\alpha m_{z}\Bigr{]}.$$

$$\displaystyle\dot{\Theta}=\frac{\gamma}{1+\alpha^{2}}\Bigl{[}-\alpha H_{eff,y}\sin\Theta+$$

$$\displaystyle H_{eff,x}(\sin\Phi+\alpha\cos\Theta\cos\Phi)\Bigr{]},$$

$$\displaystyle\dot{\Phi}\sin\Theta=\frac{\gamma}{1+\alpha^{2}}\Bigl{[}-H_{eff,y}\sin\Theta+$$

$$\displaystyle H_{eff,x}(-\alpha\sin\Phi+\cos\Theta\cos\Phi)\Bigr{]}.$$

Full text

I Introduction

By now it is well-known that current-phase relation (CPR) in Josephson junctions with multi-layered ferromagnetic interlayers is strongly sensitive to the mutual orientation of the magnetizations in the layers Waintal2002 ; Barash2004 ; Braude2007 ; Grein2009 ; Liu2010 ; Kulagina2014 ; Moor2015_1 ; Moor2015_2 ; Mironov2015 ; Silaev2017 ; Bobkova2017 ; Rabinovich2018 . CPRs of Josephson junctions with ferromagnetic interlayers in the presence of spin-orbit coupling also depends on the magnetization orientation. This occurs primarily via the appearance of the magnetization-dependent anomalous phase shift Krive2004 ; Asano2007 ; Reynoso2008 ; Buzdin2008 ; Tanaka2009 ; Zazunov2009 ; Malshukov2010 ; Alidoust2013 ; Brunetti2013 ; Yokoyama2014 ; Bergeret2015 ; Campagnano2015 ; Konschelle2015 ; Kuzmanovski2016 . This coupling between the Josephson and magnetic subsystems leads to the supercurrent-induced magnetization dynamics Konschelle2009 ; Kulagina2014 ; Waintal2002 ; Braude2008 ; Linder2011 ; Cai2010 ; Chudnovsky2016 ; Bobkova2018 . In particular, the reversal of the magnetic moment by the supercurrent pulse Shukrinov2017 was predicted. A unique possibility of controlling the magnetization dynamics via external bias current and series of specific magnetization trajectories has been reported Shukrinov2019 . In Refs. Konschelle2009, ; Shukrinov2018, it was also reported that in the presence of spin-orbit coupling the supercurrent can cause reorientation of the magnetization easy-axis. Assuming the initial position of the easy axis along $z$-direction these works demonstrate that under the applied supercurrent stable position of the magnetization becomes between $z$− and $y$-axes depending on parameters of the system.

Here we investigate prospects of superconductor/ferromagnet/superconductor (S/F/S) Josephson junctions constructed atop a three dimensional topological insulator (3D TI) surface, which hosts Dirac quasiparticles, in the field of supercurrent-induced magnetization control. Our motivation is that these Dirac quasiparticles on the surface of the 3D TI exhibit full spin-momentum locking: an electron spin always makes a right angle with its momentum. This gives rise to a very pronounced dependence of the CPR on the magnetization direction Tanaka2009 ; Linder2010 ; Zyuzin2016 . In particular, the anomalous ground state phase shift proportional to the in-plane magnetization component perpendicular to the supercurrent direction was reported.

The second reason to study magnetization dynamics in such a system is that at present there is a great progress in experimental realization of F/TI hybrid structures. In particular, to introduce the ferromagnetic order into the TI, random doping of transition metal elements, e.g., Cr or V, has been employed Chang2013 ; Kou2013 ; Kou2013_2 ; Chang2015 . The second option, which has been successfully realized experimentally, is a coupling of the nonmagnetic TI to a high $T_{c}$ magnetic insulator to induce strong exchange interaction in the surface states via the proximity effectJiang2015 ; Wei2013 ; Swartz2012 ; Jiang2015_2 ; Jiang2016 .

Here we demonstrate that the anomalous phase shift causes the magnetization dynamics analogously to the case of a spin-orbit coupled system. However, in contrast to the spin-orbit coupled systems, where the magnetization dynamics was studied before, for the system under consideration the absolute value of the critical current also depends strongly on the magnetization orientation. It only depends on the in-plane magnetization component along the current direction. We demonstrate that such dependence, in a suitably chosen voltage-driven regime, can lead to supercurrent induced splitting of the magnetic easy axis of the ferromagnet. We show that this effect may lead to stabilization of a four-fold degenerate ferromagnetic state, which is in sharp contrast to the conventional two-fold degenerate easy-axis ferromagnetic state.

The paper is organized as follows. In Sec. II we derive a CPR for the S/F/S junction atop a topological insulator surface starting from the quasiclassical Green function formalism. This is followed by a discussion of the magnetization dynamics of such systems in Sec III. Next, in Sec IV, we discuss the stabilization of the four-fold degenerate ferromagnetic state. Finally, we conclude in Sec. V.

II Current-phase relation in a ballistic S/F/S junction on a 3D TI

The sketch of the system under consideration is presented in Fig. 1. Two conventional s-wave superconductors and a ferromagnet are deposited on top of a 3D TI insulator to form a Josephson junction.

First of all, we consider a current-phase relation of a Josephson junction. The interlayer of the junction consists of the TI conducting surface states with a ferromagnetic layer on top of it. It is assumed that the magnetization $\bm{M}(\bm{r})$ of the ferromagnet induces an effective exchange field $\bm{h}_{eff}(\bm{r})\sim\bm{M}(\bm{r})$ in the underlying conductive surface layer. The Hamiltonian that describes the TI surface states in the presence of an in-plane exchange field $\bm{h}_{eff}(\bm{r})$ reads:

[TABLE]

[TABLE]

where $\hat{\Psi}=(\Psi_{\uparrow},\Psi_{\downarrow})^{T}$, $v_{F}$ is the Fermi velocity, $\bm{e}_{z}$ is a unit vector normal to the surface of TI, $\mu$ is the chemical potential, and $\hat{\bm{\sigma}}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is a vector of Pauli matrices in the spin space. It was shown Zyuzin2016 ; Bobkova2016 that in the quasiclassical approximation $(h_{eff},\varepsilon,\Delta)\ll\mu$ the Green’s function has the following spin structure: $\check{g}(\bm{n}_{F},\bm{r},\varepsilon)=\hat{g}(\bm{n}_{F},\bm{r},\varepsilon)(1+\bm{n}_{\perp}\bm{\sigma})/2$, where $\bm{n}_{\perp}=(n_{F,y},-n_{F,x},0)$ is the unit vector perpendicular to the direction of the quasiparticle trajectory $\bm{n}_{F}=\bm{p}_{F}/p_{F}$ and $\hat{g}$ is the spinless $4\times 4$ matrix in the particle-hole and Keldysh spaces containing normal and anomalous quasiclassical Green’s functions. The spin structure above reflects the fact that the spin and momentum of a quasiparticle at the surface of the 3D TI are strictly locked and make a right angle. Following standard proceduresEilenberger1968 ; Usadel1970 it was demonstratedZyuzin2016 ; Bobkova2016 ; Hugdal2017 that the spinless retarded Green’s function $\hat{g}(\bm{n}_{F},\bm{r},\varepsilon)$ obeys the following transport equations in the ballistic limit:

[TABLE]

where $[A,B]_{\otimes}=A\otimes B-B\otimes A$ and $A\otimes B=\exp[(i/2)(\partial_{\varepsilon_{1}}\partial_{t_{2}}-\partial_{\varepsilon_{2}}\partial_{t_{1}})]A(\varepsilon_{1},t_{1})B(\varepsilon_{2},t_{2})|_{\varepsilon_{1}=\varepsilon_{2}=\varepsilon;t_{1}=t_{2}=t}$. $\tau_{x,y,z}$ are Pauli matrices in particle-hole space with $\tau_{\pm}=(\tau_{x}\pm i\tau_{y})/2$. $\hat{\Delta}=\Delta(x)\tau_{+}-\Delta^{*}(x)\tau_{-}$ is the matrix structure of the superconducting order parameter $\Delta(x)$ in the particle-hole space. We assume $\Delta(x)=\Delta e^{-i\chi/2}\Theta(-x-d/2)+\Delta e^{i\chi/2}\Theta(x-d/2)$. The spin-momentum locking allows for including $\bm{h}_{eff}$ into the gauge-covariant gradient $\hat{\nabla}\hat{A}=\bm{\nabla}\hat{A}+(i/v_{F})[(h_{x}\bm{e}_{y}-h_{y}\bm{e}_{x})\tau_{z},\hat{A}]_{\otimes}$.

Eq. (3) should be supplemented by the normalization condition $\hat{g}\otimes\hat{g}=1$ and the boundary conditions at $x=\mp d/2$. As we assume that the Josephson junction is formed at the surface of the TI, the superconducting order parameter $\Delta$ and $\bm{h}_{eff}$ are effective quantities induced in the surface states of TI by proximity to the superconductors and a ferromagnet. In this case there are no reasons to assume existence of potential barriers at the $x=\mp d/2$ interfaces and we consider these interfaces as fully transparent. In this case the boundary conditions are extremely simple and are reduced to continuity of $\hat{g}$ for a given quasiparticle trajectory at the interfaces.

To obtain the simplest sinusoidal form of the current-phase relation we linearize Eq. (3) with respect to the anomalous Green’s function. In this case the retarded component of the Green’s function $\hat{g}^{R}=\tau_{z}+f^{R}\tau_{+}+\tilde{f}^{R}\tau_{-}$. The anomalous Green’s function obeys the following equation:

[TABLE]

Equation for $\tilde{f}^{R}$ is obtained from Eq. (4) by $\bm{v}_{F}\to-\bm{v}_{F}$, $\Delta\to-\Delta$ and $\chi\to-\chi$.

The solution of Eq. (4) satisfying asymptotic conditions $f^{R}\to(\Delta/\varepsilon)e^{\pm i\chi/2}$ at $x\to\pm\infty$ and continuity conditions at $x=\mp d/2$ takes the form [the solution is written for $x\in(-d/2,d/2)$, the solution in the superconducting leads is also found, but it is not required for finding the Josephson current]:

[TABLE]

where the subscript $\pm$ corresponds to the trajectories ${\rm sgn}\!~{}v_{x}=\pm 1$.

The density of electric current along the $x$-axis is

[TABLE]

where $\phi$ is the angle, which the quasiparticle trajectory makes with the $x$-axis. $\varphi_{\pm}$ is the distribution function corresponding to the trajectories ${\rm sgn}\!~{}v_{x}=\pm 1$.

Here we consider the voltage-biased junction. In principle, in this case the electric current through the junction consists of two parts: the Josephson current $j_{s}$ and the normal current $j_{n}$. The Josephson current is connected to presence of the nonzero anomalous Green’s functions in the interlayer and takes place even in equilibrium. Here we assume the low applied voltage regime $eV/(k_{B}T_{c})\ll 1$. In this case the deviation of the distribution function from equilibrium is weak and can be disregarded in calculation of the Josephson current: $\varphi_{+}=\varphi_{-}=\tanh(\varepsilon/2T)$. Exploiting the normalization condition one can obtain $g^{R}_{\pm}\approx 1-f^{R}_{\pm}\tilde{f}^{R}_{\pm}/2$. Taking into account that $g^{A}_{\pm}=-g^{R*}_{\pm}$ we find the following final expression for the Josephson current:

[TABLE]

where $\varepsilon_{n}=\pi T(2n+1)$. At high temperatures $T\approx T_{c}\gg\Delta$ the main contribution to the current comes from the lowest Matsubara frequency and Eq. (8) can be simplified further

[TABLE]

where $j_{b}=ev_{F}N_{F}\Delta^{2}/(\pi^{2}T)$. Similar expression has already been obtained for Dirac materials Hugdal2017 . The normal current is due to deviation of the distribution function from the equilibrium. However, for the system under consideration, where we assume the ferromagnet to be metallic, practically all the normal current flows through the ferromagnet because in real experimental setups the TI resistance should be much larger as compared to the resistance of the ferromagnet. As for the Josephson current, it is carried by Cooper pairs and is strongly suppressed inside the ferromagnetic layer because the exchange field there is typically much larger as compared to the induced exchange field $\bm{h}_{eff}$ in the TI surface layer. Therefore, it flows through the TI surface states and we can assume that it is equal to the total electric current flowing via the TI surface states.

III Magnetization dynamics induced by a coupling to Josephson junction

The dynamics of the ferromagnet magnetization can be described in the framework of the Landau-Lifshitz-Gilbert (LLG) equation

[TABLE]

where $M_{s}$ is the saturation magnetization, $\gamma$ is the gyromagnetic ratio and $H_{eff}$ is the local effective field. The electric current flowing through the TI surface states causes spin-orbital torqueYokoyama2010 ; Yokoyama2011 ; Mahfouzi2012 ; Chen2014 due to the presence of a strong coupling between a quasiparticle spin and momentum. In principle, if the ferromagnetism and spin-orbit coupling spatially coexist, this torque is determined by the total electric current flowing through the system. However, for the case under consideration only the supercurrent flows via the TI surface states, where the spin-momentum locking takes place. Therefore, only this supercurrent generates a torque acting on the magnetization. The normal current flows through the homogeneous ferromagnet, where we assume no spin-orbit coupling. Consequently, it does not contribute to the torque.

The torque caused by the supercurrent can be accounted for as an additional contribution to the effective field. In order to find this contribution we can consider the energy of the junction as a sum of the magnetic and the Josephson energies:

[TABLE]

where $E_{J}=E_{c}[1-\cos(\chi-\chi_{0})]$ with $E_{c}=\Phi_{0}I_{c}/2\pi$, $I_{c}=j_{c}S$ ($S$ is the junction area) and $\chi=2eVt$ in the presence of the applied voltage. is the Josephson energy. $E_{M}=-KV_{F}M_{y}^{2}/2M_{s}^{2}$ is the uniaxial anisotropy energy with the easy axis assumed to be along the $y$-axis. $V_{F}$ is the volume of the ferromagnet. The effective field $H_{eff}=-(1/V_{F})(\delta E_{tot}/\delta\bm{M})$ and takes the form:

[TABLE]

where we have introduced the unit vector $\bm{m}=\bm{M}/M_{s}$, $\tilde{d}=2\pi Td/v_{F}$ is the dimensionless junction length, $\Gamma=\Phi_{0}j_{b}Sr/2\pi KV_{F}$ is proportional to the ratio of the Josephson and magnetic energies, $r=2dh_{eff}/v_{F}$, $\Omega_{J}=2eV$ is the Josephson frequency and $H_{F}=\Omega_{F}/\gamma=K/M_{s}$.

The effective field consists of two contributions: the anisotropy field, which is directed along the easy axis, is represented by the last term in Eq. (14). The other terms are generated by the supercurrent. The same approach to study magnetization dynamics in voltage biased junctions has already been applied to systems with spin-orbit coupling in the interlayer Konschelle2009 ; Shukrinov2018 . The qualitative difference of our system based on the TI surface states from these works is that the critical current demonstrates strong dependence on the $x$-component of magnetization in our case, while it has been considered as independent on the magnetization direction earlier. This dependence leads to nonzero $H_{eff,x}\sim m_{x}$ at small $m_{x}$. This means that the easy $y$-axis can become unstable in a voltage-driven or current-driven junction, while this axis is always stable if the critical current does not depend on magnetization direction. Moreover, there is no difference for the system between $\pm m_{x}$-components of the magnetization. This leads to the remarkable fact that in a driven system the easy axis is not reoriented keeping two stable magnetization directions, as it has already been obtained before, but is split demonstrating four stable magnetization directions. In the following section we study this effect in detail.

IV Magnetization easy axis splitting

It is obvious that $m_{x}=m_{z}=0$ is an equilibrium point of Eq. (11) with $H_{eff}$ determined by Eqs. (13)-(15). Now we investigate stability of this point. In the linear order with respect to $m_{x}$ the effective field can be written as follows:

[TABLE]

where $A>0$ and $B>0$ are constants. By comparing Eqs. (16) and (13) it is seen that

[TABLE]

The LLG equation (11) in the linear order with respect to $m_{x}$ and $m_{z}$ takes the form:

[TABLE]

while $\dot{m}_{y}=0$.

One can estimate the parameter $\Omega_{F}/\Omega_{J}\sim\gamma H_{F}/(eI_{c}R_{n})$ for 3D TI Josephson junctions. Taking $I_{c}R_{n}\sim 10^{-3}V$, as has been reported for $Nb/Bi_{2}Te_{3}/Nb$ Josephson junctions Veldhorst2012 , and the easy-axis anisotropy field $H_{F}\sim 500Oe$, what was reported for Py Beach2005 ; Beach2006 , we obtain $\Omega_{F}/\Omega_{J}\sim 3\times 10^{-3}$. In the regime $\Omega_{F}/\Omega_{J}\ll 1$ the magnetization varies slowly at $t\sim\Omega_{J}^{-1}$, therefore we can average Eqs. (18) over a Josephson period thus obtaining the following system:

[TABLE]

The general solution of this system takes the form $m_{x(z)}=\sum\limits_{k=1,2}C_{k,x(z)}\exp[\lambda_{k}t]$. The equilibrium point $m_{x}=m_{z}=0$ becomes unstable under the condition ${\rm Re}\lambda_{1}>0$ or ${\rm Re}\lambda_{2}>0$. One can obtain that it is realized at $A>1$.

It is rather difficult to make accurate estimates of the numerical value of $A$ for realistic parameters. The main problem is the absence of experimental data on the value of $h_{eff}$. However, if we take $K=500J/m^{3}$ from the measurements Rusanov2004 on permalloy with very weak anisotropy, $I_{c}=10\mu A$, $v_{F}\sim 10^{5}m/c$ from Veldhorst2012 and the permalloy volume $d\times l\times w\sim 100nm\times 10nm\times 50nm$, then we can obtain $A\sim r\Gamma\sim I_{c}h_{eff}/(v_{F}eKlw)\sim 0.4-8$ for $h_{eff}\sim 10-200K$. Therefore, we can conclude that the range of $A$ values discussed in our work should be experimentally accessible.

Now we turn to study the stationary points of the magnetization and their stability. Beyond the linear approximation (with respect to $m_{x}$ and $m_{z}$) it is convenient to parametrize the magnetization as $\bm{m}=(\sin\Theta\cos\Phi,\cos\Theta,\sin\Theta\sin\Phi)$. Then from LLG equation one obtains:

[TABLE]

At $\Omega_{F}/\Omega_{J}\to 0$ effective fields $H_{eff,x,y}$ determined by Eqs. (13), (14) should be averaged over a Josephson period $H_{eff,x,y}\to\langle H_{eff,x,y}\rangle$. The stationary points are to be found as solutions of Eqs. (20) corresponding to $\dot{\Theta}=\dot{\Phi}=0$.

Fig. 2 shows vector fields in the plane $0\leq\Phi<2\pi$, $0\leq\Theta<\pi$ according to Eq. (20) at four different values of $A$. The stationary solutions are indicated by color points. The blue points correspond to unstable stationary solutions, while the red points indicate the stable magnetization directions. The Gilbert damping constant $\alpha=0.25$. We have chosen such a large unrealistic value of the Gilbert constant in order to clearly show the stability/instability of the stationary points because for $\alpha=0.01$, which is more appropriate for a realistic situation, stability/instability of a solution is not clearly seen in the figure [compare Figs. 3(a) and (b)], although in fact the topology of the vector field is not changed. Fig. 2(a) represents the regime $A<1$, when the only stable solutions $\bm{m}^{st}$ are $m_{y}^{st}=\pm 1$, what corresponds to upper and bottom horizontal lines in the figure. Panels (b) and (c) demonstrate the vector fields in the regime of not very large $A>1$. Four stable red points are clearly seen. Upon further increase of $A$ the stable points get closer to $\Theta=\pi/2$ and finally merge into two stable points at some $A_{crit}$, as it is shown in Fig. 2(d). Therefore, there exists a finite range of $1<A<A_{crit}$, where the ferromagnet has four stable magnetization directions in the voltage-biased regime considered here. From Fig.2 it is seen that all the stationary points correspond to $m_{z}=\pm 1$ or $m_{z}=0$. The stationary points $m_{z}=\pm 1$ are always unstable. Let us consider the stationary points corresponding to $m_{z}=0$, that is $\Phi=0,\pi$. It is obvious that in order to have four stable points $|m_{x}^{st}|$ and $|m_{y}^{st}|$ should be nonzero simultaneously. Substituting $m_{z}=0$ into Eq. (20) and taking into account that $\langle H_{eff,y}\rangle=H_{F}m_{y}$, we obtain that $m_{x}^{st}$ can be determined from the simple nonlinear equation $m_{x}=\langle H_{eff,x}\rangle/H_{F}$. This equation always has the solution $m_{x}=0$, but at $1<A<A_{crit}$ it also has the second nonzero solution $m_{x}^{st}$. In this situation there are four stable points $\bm{m}^{st}=(\pm|m_{x}^{st}|,\pm|m_{y}^{st}|,0)$.

Further in Fig. 4 we demonstrate the full time evolution of the magnetization $\bm{m}$ obtained from the numerical solution of the LLG equation. It is seen that starting from different initial conditions it is possible to reach all four stable magnetization solutions. The results are obtained at $\Omega_{F}/\Omega_{J}=0.2$, but the averaged values of magnetization at large times are in good agreement with the results for stable points obtained in the limit $\Omega_{F}/\Omega_{J}\ll 1$, which are demonstrated in Fig. 3(a) for the same parameters $\Gamma$, $r$, $\alpha$ and $d$. Fig. 3(b) only differs from (a) by the value of $\alpha=0.25$. While the topology of the vector fields presented in panels (a) and (b) is the same, the stability/instability of all the stationary points is more clearly seen for larger values of the damping constant $\alpha$. At finite values of $\Omega_{F}/\Omega_{J}$ the magnetization oscillates around the vector trajectory presented in Fig. 3 and the amplitude of the oscillations is suppressed at $\Omega_{F}/\Omega_{J}\to 0$.

In order to show that the system under consideration can demonstrate spontaneous behavior we investigate the system evolution starting from one of unstable points $m_{z}=\pm 1$. A small noise is introduced to the system in order to allow for leaving the unstable equilibrium point. From the vector fields represented in Fig. 3(a) it is seen that at small values of $\alpha$ the system finally comes to one of the four stable states with practically equal probabilities. It is shown in Fig. 5, where different panels correspond to all the possible final states.

Fig. 6 demonstrates the behavior of the absolute values of averaged magnetization at $t\to\infty$ depending on essential parameters of the system. The dependence on $\Omega_{J}/\Omega_{F}$ is represented in Fig. 6(a). It is seen that at $\Omega_{J}/\Omega_{F}\gg 1$ $|\langle m_{i}\rangle|$ tend to constant values and, in particular, $|\langle m_{z}\rangle|\to 0$, as it follows from our analysis of stationary points of Eqs. (20).

The dependence on $\Gamma$ is plotted in Fig. 6(b). $\Gamma$ is linearly proportional to $A$. For this reason one can explicitly see in this panel the range of $A$ where four stable limiting magnetization directions exist: it corresponds to the regions, where $|\langle m_{x}\rangle|$ and $|\langle m_{y}\rangle|$ are nonzero simultaneously.

Panel (c) of Fig. 6 represents the dependence of $|\langle m_{i}\rangle|$ on the junction length. Analogously to the previous panel, the range of existence of four stable limiting magnetization directions is also clearly seen. The dependence on $r$ is qualitatively very similar to the dependence on $\Gamma$, therefore we do not represent it. Figs. 6(b)-(c) also provide the optimal range of parameters $\Gamma$ and $d$ for experimental realization of the easy axis splitting. The effect can be experimentally investigated, for example, by measuring the magnetic field pattern away from the nanomagnet.

V Conclusions

In this work we study a S/F/S Josephson junction atop a topological insulator and discuss the possibility of electrical control of magnetization in such a device. Our analysis, which combines microscopic Keldysh Green function tehnique for obtaining the Josephson current with phenomenological Landau-Lifshitz-Gilbert equations for studying magnetization dynamics, shows that the property of full spin momentum locking can lead to destabilization of a transverse easy magnetization axis $m_{y}$ in such systems in the presence of appropriate voltage or current bias. Such an instability, in its turn, results in a ferromagnet with two easy axes allowing for four stable magnetization directions instead of usual two. Switching between these states by means of voltage or current impulses is of interest from the applied point of view.

Acknowledgements

The work of I.V.B. and A.M.B. was carried out within the state task of ISSP RAS. The reported study was partially funded by the RFBR research projects 18-02-00318 and 18-52-45011-IND. Numerical calculations have been made in the framework of the RSF project 18-71-10095. K.S. thanks DST for support through INT/RUS/RFBR/P-314.