Dynamic coupling of ferromagnets via spin Hall magnetoresistance

TL;DR

This paper theoretically investigates how spin Hall magnetoresistance mediates dynamic coupling and synchronization of magnetization in ferromagnets placed on a nonmagnetic heavy metal, revealing various phase-locked states.

Contribution

It introduces a theoretical model demonstrating how spin Hall magnetoresistance causes coupled magnetization dynamics and synchronization in ferromagnetic systems.

Findings

In-phase and antiphase synchronization depending on damping and torque.

Quarter-period phase difference observed in longitudinal geometry.

Analytical relation between current, frequency, and phase difference established.

Abstract

The synchronized magnetization dynamics in ferromagnets on a nonmagnetic heavy metal caused by the spin Hall effect is investigated theoretically. The direct and inverse spin Hall effects near the ferromagnetic/nonmagnetic interface generate longitudinal and transverse electric currents. The phenomenon is known as the spin Hall magnetoresistance effect, whose magnitude depends on the magnetization direction in the ferromagnet due to the spin transfer effect. When another ferromagnet is placed onto the same nonmagnet, these currents are again converted to the spin current by the spin Hall effect and excite the spin torque to this additional ferromagnet, resulting in the excitation of the coupled motions of the magnetizations. The in-phase or antiphase synchronization of the magnetization oscillations, depending on the value of the Gilbert damping constant and the field-like torque strength, is found in the transverse geometry by solving the Landau-Lifshitz-Gilbert equation numerically. On the other hand, in addition to these synchronizations, the synchronization having a phase difference of a quarter of a period is also found in the longitudinal geometry. The analytical theory clarifying the relation among the current, frequency, and phase difference is also developed, where it is shown that the phase differences observed in the numerical simulations correspond to that giving the fixed points of the energy supplied by the coupling torque.

Equations

$$\frac{J_{\mathrm{c}x}}{J_0}=1+{\chi }^{\prime \prime }+\chi m_y^2,$$

$$\frac{J_{\mathrm{c}y}}{J_0}=-\chi m_x{m}_y-{\chi }^{\prime }{m}_z,$$

$$\begin{array}{rl}\chi =& \frac{{\vartheta }^2{\lambda }_{\mathrm{N}}}{d_{\mathrm{N}}}\left[\mathrm{Re}\frac{g^{\uparrow \downarrow }}{g_{\mathrm{N}}+g^{\uparrow \downarrow }\coth \left(d_{\mathrm{N}}/{\lambda }_{\mathrm{N}}\right)}-\frac{g^{\ast }}{g_{\mathrm{N}}}\right]{\tanh }^2\left(\frac{d_{\mathrm{N}}}{2{\lambda }_{\mathrm{N}}}\right),\end{array}$$

$$\begin{array}{r}{\chi }^{\prime }=-\frac{{\vartheta }^2{\lambda }_{\mathrm{N}}}{d_{\mathrm{N}}}\mathrm{Im}\frac{g^{\uparrow \downarrow }}{g_{\mathrm{N}}+g^{\uparrow \downarrow }\coth \left(d_{\mathrm{N}}/{\lambda }_{\mathrm{N}}\right)}{\tanh }^2\left(\frac{d_{\mathrm{N}}}{2{\lambda }_{\mathrm{N}}}\right),\end{array}$$

$$\begin{array}{rl}{\chi }^{\prime \prime }=& \frac{2{\vartheta }^2{\lambda }_{\mathrm{N}}}{d_{\mathrm{N}}}\tanh \left(\frac{d_{\mathrm{N}}}{2{\lambda }_{\mathrm{N}}}\right)\\ & -\frac{{\vartheta }^2{\lambda }_{\mathrm{N}}}{d_{\mathrm{N}}}\mathrm{Re}\frac{g^{\uparrow \downarrow }}{g_{\mathrm{N}}+g^{\uparrow \downarrow }\coth \left(d_{\mathrm{N}}/{\lambda }_{\mathrm{N}}\right)}{\tanh }^2\left(\frac{d_{\mathrm{N}}}{2{\lambda }_{\mathrm{N}}}\right),\end{array}$$

$$\frac{1}{g^{\ast }}=\frac{2}{(1-p_g^2)g}+\frac{1}{g_{\mathrm{F}}\tanh \left(d_{\mathrm{F}}/{\lambda }_{\mathrm{F}}\right)}+\frac{1}{g_{\mathrm{N}}\tanh \left(d_{\mathrm{N}}/{\lambda }_{\mathrm{N}}\right)}.$$

$${\vartheta }_{\mathrm{R}(\mathrm{I})}=\vartheta \tanh \left(\frac{d_{\mathrm{N}}}{2{\lambda }_{\mathrm{N}}}\right)\mathrm{Re}\left(\mathrm{Im}\right)\frac{g^{\uparrow \downarrow }}{g_{\mathrm{N}}+g^{\uparrow \downarrow }\coth \left(d_{\mathrm{N}}/{\lambda }_{\mathrm{N}}\right)}.$$

$$\beta =-\frac{{\vartheta }_{\mathrm{I}}}{{\vartheta }_{\mathrm{R}}},$$

$$\begin{array}{rl}{j}_{\ell x}& =\frac{\left({\chi }^{\prime \prime }+\chi m_{\ell y}^2\right)\left(1+{\chi }^{\prime \prime }+\chi m_{{\ell }^{\prime }y}^2\right)}{1-\left({\chi }^{\prime \prime }+{\chi }^{\prime }{m}_{\ell y}^2\right)\left({\chi }^{\prime \prime }+\chi m_{{\ell }^{\prime }y}^2\right)}{J}_0\\ & \simeq \left({\chi }^{\prime \prime }+\chi m_{\ell y}^2\right)J_0,\end{array}$$

$$\begin{array}{rl}{\mathbf{T}}_{\ell }^{\mathrm{L}}=& -\frac{\gamma \hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}\left(1+{\chi }^{\prime \prime }+\chi m_{{\ell }^{\prime }y}^2\right){\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right)\\ & -\frac{\gamma \hslash {\vartheta }_{\mathrm{I}}{J}_0}{2eM{d}_{\mathrm{F}}}\left(1+{\chi }^{\prime \prime }+\chi m_{{\ell }^{\prime }y}^2\right){\mathbf{e}}_y\times {\mathbf{m}}_{\ell },\end{array}$$

$$\begin{array}{rl}{\mathbf{T}}_{\ell }^{(0)}=& -\frac{\gamma \hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right)-\frac{\gamma \hslash {\vartheta }_{\mathrm{I}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{e}}_y\times {\mathbf{m}}_{\ell },\end{array}$$

$$\begin{array}{rl}& -\frac{\gamma \hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}\left({\chi }^{\prime \prime }+\chi m_{{\ell }^{\prime }y}^2\right){\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right)\\ & -\frac{\gamma \hslash {\vartheta }_{\mathrm{I}}{J}_0}{2eM{d}_{\mathrm{F}}}\left({\chi }^{\prime \prime }+\chi m_{{\ell }^{\prime }y}^2\right){\mathbf{e}}_y\times {\mathbf{m}}_{\ell },\end{array}$$

$$\begin{array}{rl}{j}_{\ell y}& =-[\frac{\left(\chi m_{\ell x}{m}_{\ell y}+{\chi }^{\prime }{m}_{\ell z}\right)}{1-\left({\chi }^{\prime \prime }+\chi m_{\ell x}^2\right)\left({\chi }^{\prime \prime }+\chi m_{{\ell }^{\prime }x}^2\right)}\\ & +\frac{\left({\chi }^{\prime \prime }+\chi m_{\ell x}^2\right)\left(\chi m_{{\ell }^{\prime }x}{m}_{{\ell }^{\prime }y}+{\chi }^{\prime }{m}_{{\ell }^{\prime }z}\right)}{1-\left({\chi }^{\prime \prime }+\chi m_{\ell x}^2\right)\left({\chi }^{\prime \prime }+\chi m_{{\ell }^{\prime }x}^2\right)}]J_0\\ & \simeq -\left(\chi m_{\ell x}{m}_{\ell y}+{\chi }^{\prime }{m}_{\ell z}\right)J_0.\end{array}$$

$$\begin{array}{rl}{\mathbf{T}}_{\ell }^{\mathrm{T}}=& -\frac{\gamma \hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right)-\frac{\gamma \hslash {\vartheta }_{\mathrm{I}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\\ & -\frac{\gamma \hslash {\vartheta }_{\mathrm{R}}\left(\chi m_{{\ell }^{\prime }x}{m}_{{\ell }^{\prime }y}+{\chi }^{\prime }{m}_{{\ell }^{\prime }z}\right)J_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_x\times {\mathbf{m}}_{\ell }\right)\\ & -\frac{\gamma \hslash {\vartheta }_{\mathrm{I}}\left(\chi m_{{\ell }^{\prime }x}{m}_{{\ell }^{\prime }y}+{\chi }^{\prime }{m}_{{\ell }^{\prime }z}\right)J_0}{2eM{d}_{\mathrm{F}}}{\mathbf{e}}_x\times {\mathbf{m}}_{\ell },\end{array}$$

$$\begin{array}{rl}\frac{d{\mathbf{m}}_{\ell }}{dt}=& -\gamma {\mathbf{m}}_{\ell }\times {\mathbf{H}}_{\ell }+\alpha {\mathbf{m}}_{\ell }\times \frac{d{\mathbf{m}}_{\ell }}{dt}+{\mathbf{T}}_{\ell }^{\mathrm{L},\mathrm{T}},\end{array}$$

$$\begin{array}{rl}& \left(1+{\alpha }^2\right)\frac{d{\mathbf{m}}_{\ell }}{dt}=-\gamma {\mathbf{m}}_{\ell }\times {\mathbf{H}}_{\ell }\\ & -\gamma \left(1+{\chi }^{\prime \prime }\right)\frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right)-\alpha \gamma {\mathbf{m}}_{\ell }\times \left({\mathbf{m}}_{\ell }\times {\mathbf{H}}_{\ell }\right)\\ & -\gamma \chi m_{{\ell }^{\prime }y}^2\frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right)\\ & -\gamma \left(1+{\chi }^{\prime \prime }\right)\left(\alpha +\beta \right)\frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times {\mathbf{e}}_y\\ & +\gamma \left(1+{\chi }^{\prime \prime }\right)\alpha \beta \frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right)\\ & -\gamma \left(\alpha +\beta \right)\chi m_{{\ell }^{\prime }y}^2\frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times {\mathbf{e}}_y\\ & +\gamma \alpha \beta \chi m_{{\ell }^{\prime }y}^2\frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right),\end{array}$$

$$\begin{array}{rl}& \left(1+{\alpha }^2\right)\frac{d{\mathbf{m}}_{\ell }}{dt}=-\gamma {\mathbf{m}}_{\ell }\times {\mathbf{H}}_{\ell }\\ & -\gamma \frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right)-\alpha \gamma {\mathbf{m}}_{\ell }\times \left({\mathbf{m}}_{\ell }\times {\mathbf{H}}_{\ell }\right)\\ & -\gamma \left(\chi m_{{\ell }^{\prime }x}{m}_{{\ell }^{\prime }y}+{\chi }^{\prime }{m}_{{\ell }^{\prime }z}\right)\frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_x\times {\mathbf{m}}_{\ell }\right)\\ & -\gamma \left(\alpha +\beta \right)\frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times {\mathbf{e}}_y\\ & +\gamma \alpha \beta \frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right)\\ & -\gamma \left(\alpha +\beta \right)\left(\chi m_{{\ell }^{\prime }x}{m}_{{\ell }^{\prime }y}+{\chi }^{\prime }{m}_{{\ell }^{\prime }z}\right)\frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times {\mathbf{e}}_x\\ & +\gamma \alpha \beta \left(\chi m_{{\ell }^{\prime }x}{m}_{{\ell }^{\prime }y}+{\chi }^{\prime }{m}_{{\ell }^{\prime }z}\right)\frac{\hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_x\times {\mathbf{m}}_{\ell }\right).\end{array}$$

$${\mathbf{H}}_{\ell }=H_{\mathrm{K}}{m}_{\ell y}{\mathbf{e}}_y-4\pi M{m}_{\ell z}{\mathbf{e}}_z,$$

$$J_{\mathrm{c}}=\frac{2\alpha eM{d}_{\mathrm{F}}}{\hslash {\vartheta }_{\mathrm{R}}}\left(H_{\mathrm{K}}+2\pi M\right),$$

$$J^{\ast }=\frac{4\alpha eM{d}_{\mathrm{F}}}{\pi \hslash {\vartheta }_{\mathrm{R}}}\sqrt{4\pi M\left(H_{\mathrm{K}}+4\pi M\right)},$$

$$\begin{array}{rl}{J}_{\mathrm{c}}^{\mathrm{T}}& =j_{1y}+j_{3y}\\ & =-\left(\chi m_{1x}{m}_{1y}+{\chi }^{\prime }{m}_{1z}\right)J_0\\ & -\left(\chi m_{3x}{m}_{3y}+{\chi }^{\prime }{m}_{3z}\right)J_0.\end{array}$$

$$\begin{array}{rl}{J}_{\mathrm{c}}^{\mathrm{L}}& =J_0+j_{1x}+j_{2x}\\ & =J_0+\left({\chi }^{\prime \prime }+\chi m_{1y}^2\right)J_0+\left({\chi }^{\prime \prime }+\chi m_{2y}^2\right)J_0.\end{array}$$

$$\begin{array}{rl}\frac{d{\mathbf{m}}_{\ell }}{dt}\simeq & -\gamma {\mathbf{m}}_{\ell }\times {\mathbf{H}}_{\ell }-\alpha \gamma {\mathbf{m}}_{\ell }\times \left({\mathbf{m}}_{\ell }\times {\mathbf{H}}_{\ell }\right)\\ & -\frac{\gamma \hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}\left(1+{\chi }^{\prime \prime }+\chi m_{{\ell }^{\prime }y}^2\right){\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right),\end{array}$$

$$\begin{array}{rl}\frac{d{\mathbf{m}}_{\ell }}{dt}\simeq & -\gamma {\mathbf{m}}_{\ell }\times {\mathbf{H}}_{\ell }-\alpha \gamma {\mathbf{m}}_{\ell }\times \left({\mathbf{m}}_{\ell }\times {\mathbf{H}}_{\ell }\right)\\ & -\frac{\gamma \hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_y\times {\mathbf{m}}_{\ell }\right)\\ & -\frac{\gamma \hslash {\vartheta }_{\mathrm{R}}{J}_0}{2eM{d}_{\mathrm{F}}}\chi m_{{\ell }^{\prime }x}{m}_{{\ell }^{\prime }y}{\mathbf{m}}_{\ell }\times \left({\mathbf{e}}_x\times {\mathbf{m}}_{\ell }\right).\end{array}$$

$$E=-\frac{M{H}_{\mathrm{K}}}{2}{m}_{\ell y}^2+2\pi M^2{m}_{\ell z}^2.$$

$$m_x=\sqrt{1+\frac{2E}{M{H}_{\mathrm{K}}}}\mathrm{sn}\left[\frac{4\mathsf{K}(k)}{\tau (E)}t+{\phi }_0,k\right],$$

$$m_y=\sqrt{\frac{4\pi M-2E/M}{H_{\mathrm{K}}+4\pi M}}\mathrm{dn}\left[\frac{4\mathsf{K}(k)}{\tau (E)}t+{\phi }_0,k\right],$$

$$m_z=\sqrt{\frac{H_{\mathrm{K}}+2E/M}{H_{\mathrm{K}}+4\pi M}}\mathrm{cn}\left[\frac{4\mathsf{K}(k)}{\tau (E)}t+{\phi }_0,k\right],$$

$$f(E)=\frac{\gamma \sqrt{H_{\mathrm{K}}\left(4\pi M-2E/M\right)}}{4\mathsf{K}(k)},$$

$$k=\sqrt{\frac{4\pi M\left(H_{\mathrm{K}}+2E/M\right)}{H_{\mathrm{K}}\left(4\pi M-2E/M\right)}}.$$

Full text

I Introduction

Dynamic coupling of ferromagnets has been of interest in the field of magnetism. The dipole interaction has been the basic interaction to excite the coupled motion of the magnetizations, and is applied to magnetic recording Hubert and Schäfer (1998); Kudo et al. (2012). Another method to realize the coupling is to use the spin-transfer effect Slonczewski (1996); Berger (1996), where the application of an electric current to ferromagnetic/nonmagnetic multilayers results in the magnetization switching and self-oscillation Katine et al. (2000); Kiselev et al. (2003); Grollier et al. (2003); Rippard et al. (2004); Kubota et al. (2005); Tulapurkar et al. (2005); Houssameddine et al. (2007); Taniguchi and Imamura (2008); Slavin and Tiberkevich (2009); Kubota et al. (2013). The coupled dynamics through pure spin current generated in spin pumping Tserkovnyak et al. (2003); Heinrich et al. (2003) and nonlocal Kimura et al. (2006) geometries have also been observed. It should be emphasized that these couplings are strongly restricted by the characteristic length scales. For example, the dipole coupling decays according to the inverse cube detection law, whereas the spin transfer effect by a spin-polarized electric or pure spin current occurs in a system smaller than the spin diffusion length.

Recently, physical phenomena, such as the spin torque Liu et al. (2012a, b); Pai et al. (2012); Kim et al. (2012); Haney et al. (2013); Lee et al. (2013); Garello et al. (2014); Cubukcu et al. (2014) and magnetoresistance effects Huang et al. (2012); Nakayama et al. (2013); Althammer et al. (2013); Hahn et al. (2013); Y.-T.Chen et al. (2013); Liu et al. (2015); Avci et al. (2015); Cho et al. (2015); Kim et al. (2016), due to the spin Hall effect Dyakonov and Perel (1971); Hirsch (1999); Kato et al. (2004) in bilayers consisting of an insulating or metallic ferromagnet and a nonmagnetic heavy metal have attracted much attention. The latter, known as the spin Hall magnetoresistance, originates from the charge-spin conversion of an external electric current by the direct and inverse spin Hall effects, and has been observed by measuring the longitudinal and transverse electric currents, which are given by

[TABLE]

[TABLE]

respectively, where $J_{0}$ is the electric current density generated by the external electric field. The definitions of the dimensionless coefficients, $\chi$, $\chi^{\prime}$, and $\chi^{\prime\prime}$, are given below. It should be emphasized that the currents given by Eqs. (1) and (2) depend on the magnetization direction $\mathbf{m}=(m_{x},m_{y},m_{z})$. When another ferromagnet is placed onto the same nonmagnet, these currents will be converted to spin current by the spin Hall effect again, and excite spin torque on this additional ferromagnet. Then, the magnetization dynamics of two ferromagnets will be coupled through the angular dependencies of the electric current given by Eqs. (1) and (2). This coupling is unavoidable whenever several ferromagnets are placed onto the same nonmagnet, and is not restricted by the distance between the ferromagnets because it is carried by the electric current. Since the structure consisting of several ferromagnets on the nonmagnetic heavy metal will be important from the viewpoints of both fundamental physics and practical applications based on the spin Hall effect, such as magnetic random access memory, spin torque oscillators, and bio-inspired computing Locatelli et al. (2014); Grollier et al. (2016), it is of interest to clarify the role of this coupling.

In this paper, we investigate the coupled dynamics of magnetizations in ferromagnets in the presence of the spin Hall effect by solving the Landau-Lifshitz-Gilbert (LLG) equation numerically for both longitudinal and transverse geometries. In addition to the external electric current, the current contributing to the spin Hall magnetoresistance also excites the spin torque. The strength of this additional torque is estimated from the theory of the spin Hall magnetoresistance extended to the system consisting of several ferromagnets. The conventional spin torque is proportional to the spin Hall angle $\vartheta$, whereas the new torque is on the order of $\vartheta^{3}$, and therefore, its value is two orders of magnitude smaller than the conventional spin torque. Nevertheless, it is found that this additional new torque affects the phase difference of the magnetization in the self-oscillation state. The numerical simulation reveals that the in-phase or antiphase synchronization is observed in the transverse geometry. On the other hand, in addition to them, the phase difference becomes a quarter of a period in the longitudinal geometry. It is found that these phase differences depend on the values of the Gilbert damping constant and the dimensionless field-like torque strength. An analytical theory clarifying the relation among the current, frequency, and phase difference is also developed.

The paper is organized as follows. In Sec. II, we describe the system under consideration, and discuss the theoretical formula of the spin torque excited by the spin Hall effect in the presence of the several ferromagnets and the spin Hall magnetoresistance effect. In Sec. III, we study the phase differences in the synchronized state of the magnetizations for both the longitudinal and transverse geometries by solving the LLG equation numerically. In Sec. IV, the theory clarifying the relation among the current, frequency, and phase difference is developed based on the LLG equation averaged over constant energy curves. The summary of the paper is given in Sec. V.

II System description and LLG equation

In this section, we describe the system adopted in this study, and show the spin torque formulas including the coupling torques between the ferromagnets.

II.1 System description

The system we consider is schematically shown in Fig. 1(a), where three ferromagnets Fℓ ($\ell=1,2,3$) are placed onto the same nonmagnet N. We assume that the material parameters in the ferromagnets are identical, for simplicity. The external electric field, $E_{x}$, is applied to the $x$ direction, inducing the electric current density $J_{0}=\sigma_{\rm N}E_{x}$, where $\sigma_{\rm N}$ is the conductivity of the nonmagnet. The direct and inverse spin Hall effects produce electric currents in the longitudinal ($x$) and transverse ($y$) directions. These electric currents are converted to the spin current and injected into the F2 and F3 layers due to the spin Hall effect, resulting in the excitation of the spin torque. Then, the magnetization dynamics in the Fℓ ($\ell=2,3$) layer is affected by that in the F1 layer, and vice versa. We call the coupling between the F1 and F2 layers the longitudinal coupling, whereas that between the F1 and F3 layers the transverse coupling. We assume that both the ferromagnet and nonmagnet are metallic because metallic bilayers are generally used to measure the magnetization switching and oscillation by the spin Hall effect Liu et al. (2012a, b); Pai et al. (2012); Kim et al. (2012); Garello et al. (2014). Although the spin Hall magnetoresistance was originally studied for insulating ferromagnets Huang et al. (2012); Nakayama et al. (2013); Althammer et al. (2013); Hahn et al. (2013), a large spin Hall magnetoresistance in metallic system has also been reported recently Avci et al. (2015); Cho et al. (2015); Kim et al. (2016).

The dimensionless coefficients, $\chi$, $\chi^{\prime}$, and $\chi^{\prime\prime}$, in Eqs. (1) and (2) for single ferromagnets have been derived for an insulating Y.-T.Chen et al. (2013) or metallic Kim et al. (2016); Taniguchi (2016) ferromagnet, which are given by

[TABLE]

[TABLE]

[TABLE]

where $d_{\rm N}$ and $\lambda_{\rm N}$ are thickness and spin diffusion length of the nonmagnet, respectively, whereas $g_{\rm N}/S=h\sigma_{\rm N}/\left(2e^{2}\lambda_{\rm N}\right)$ with the cross section area of the ferromagnetic/nonmagnetic interface $S$. The dimensionless mixing conductance $g^{\uparrow\downarrow}=g_{\rm r}^{\uparrow\downarrow}+ig_{\rm i}^{\uparrow\downarrow}$ consists of its real and imaginary parts Brataas et al. (2000, 2001); Zwierzycki et al. (2005), and $g^{*}$ is defined as

[TABLE]

Here, $g=g^{\uparrow\uparrow}+g^{\downarrow\downarrow}$ is the sum of the conductances of the spin-up and spin-down electrons, whereas $p_{g}=(g^{\uparrow\uparrow}-g^{\downarrow\downarrow})/g$ is its spin polarization. The ferromagnetic/nonmagnetic interface resistance $r$ is related to $g$ via $g/S=(h/e^{2})/r$. We also introduce $g_{\rm F}/S=h\left(1-p_{\sigma}^{2}\right)\sigma_{\rm F}/\left(2e^{2}\lambda_{\rm F}\right)$, where $\sigma_{\rm F}$ is the conductivity of the ferromagnet and $p_{\sigma}$ is its spin polarization. The thickness and spin diffusion length of the ferromagnet are denoted as $d_{\rm F}$ and $\lambda_{\rm F}$, respectively. The term related to $g^{*}$ is neglected when the ferromagnet is an insulator Y.-T.Chen et al. (2013), i.e., $r\to\infty$. The following quantities correspond to the effective spin polarizations of the damping-like (or Slonczewski Slonczewski (1996)) torque and the field-like torque, respectively,

[TABLE]

For the later discussion, we introduce

[TABLE]

which corresponds to the ratio of the field-like torque to the damping-like torque.

The values of the parameters used in the following calculations are derived from recent experiments on the W/CoFeB heterostructure Kim et al. (2016), where $\rho_{\rm F}=1/\sigma_{\rm F}=1.6$ k$\Omega$nm, $p_{\sigma}=0.72$, $\lambda_{\rm F}=1.0$ nm, $\rho_{\rm N}=1/\sigma_{\rm N}=1.25$ k$\Omega$nm, $\lambda_{\rm N}=1.2$ nm, and $\vartheta=0.27$, whereas the thicknesses are assumed to be $d_{\rm F}=2$ nm and $d_{\rm N}=3$ nm. The interface conductances were not evaluated in Ref. Kim et al. (2016) by assuming a transparent interface. Instead, we use typical values of the interface conductances obtained from the first-principles calculations Zwierzycki et al. (2005), $r=0.25$ k$\Omega$nm2, $p_{g}=0.5$, and $g_{\rm r}^{\uparrow\downarrow}/S=25$ nm*-2*. We note that the imaginary part of the mixing conductance, $g_{\rm i}^{\uparrow\downarrow}$, is either positive or negative, depending on the material and thickness Zwierzycki et al. (2005). The sign of $g_{\rm i}^{\uparrow\downarrow}$ determines those of $\chi^{\prime}$ and $\vartheta_{\rm I}$, or equivalently, $\beta$. For example, when $g_{\rm i}^{\uparrow\downarrow}/S=1$ nm*-2*, $\chi\simeq 0.010$, $\chi^{\prime}\simeq-0.0002$, $\chi^{\prime\prime}\simeq 0.035$, $\theta_{\rm R}\simeq 0.167$, and $\theta_{\rm I}\simeq 0.002$ ($\beta\simeq-0.010$); see Appendix A. In the following calculations, we study the magnetization dynamics for several values of $\beta$.

II.2 Spin torques in longitudinal and transverse geometries

Equation (1) was derived for a system having single ferromagnet. In this case, $J_{0}$ is the external electric current density, whereas $(\chi^{\prime\prime}+\chi m_{y})J_{0}$ is the current density generated as a result of the charge-spin conversion by the direct and inverse spin Hall effects. In the longitudinal geometry in the present study, on the other hand, two ferromagnetic/nonmagnetic interfaces, i.e., F1/N and F2/N interfaces, contribute to the generation of the longitudinal current through the direct and inverse spin Hall effects, as schematically shown in Fig. 1(b). Let us denote the electric current density generated by these effects near the Fℓ/N interface as $j_{\ell x}$ ($\ell=1,2$). This current density is determined by the conservation law of the electric current, as follows. Considering in a similar manner to the case of the single ferromagnet, the electric current $J_{0}+j_{1x}$ is converted to the spin current by the spin Hall effect near the F2/N interface, and this spin current produces an additional electric current $(\chi^{\prime\prime}+\chi m_{2y}^{2})(J_{0}+j_{1x})$ by the inverse spin Hall effect. Therefore, the total longitudinal electric current density near the F2/N interface is $(1+\chi^{\prime\prime}+\chi m_{2y}^{2})(J_{0}+j_{1x})$. Similarly, the electric current density near the F1/N interface can be expressed as $(1+\chi^{\prime\prime}+\chi m_{1y}^{2})(J_{0}+j_{2x})$. These currents should be equal to the total electric current density, $J_{0}+j_{1x}+j_{2x}$, according to the conservation law of the electric current. Then, we find that $j_{\ell x}$ is given by

[TABLE]

where we neglect the higher order terms of $\vartheta$. Then, the spin torque excited on the Fℓ layer by the spin Hall effect in the longitudinal geometry is obtained by replacing the external electric current density $J_{0}=\sigma_{\rm N}E_{x}$ in the previous work Y.-T.Chen et al. (2013) with $J_{0}+j_{\ell^{\prime}x}\simeq\left(1+\chi^{\prime\prime}+\chi m_{\ell^{\prime}y}^{2}\right)J_{0}$, and is given by

[TABLE]

where $(\ell,\ell^{\prime})=(1,2)$ or $(2,1)$. The unit vector pointing in the magnetization direction of the Fℓ layer is $\mathbf{m}_{\ell}$, and $\gamma$, $M$, and $d$ are the gyromagnetic ratio, the saturation magnetization, and the thickness of the ferromagnet, respectively.

Note that the terms,

[TABLE]

in Eq. (10) are the conventional spin torques generated by the external electric current $J_{0}$ and the spin Hall effect, and are often referred to as the damping-like torque and the field-like torque, respectively. On the other hand, the terms

[TABLE]

in Eq. (10) originate from the current $j_{\ell^{\prime}x}$, and are newly introduced in this study. It should be emphasized that Eq. (12) depends on the magnetization direction of the other ferromagnet F${}_{\ell^{\prime}}$ ($\ell^{\prime}=1$ or $2$), $\mathbf{m}_{\ell^{\prime}}$, resulting in the coupling between the F1 and F2 layers.

Let us next show the spin torque formulas in the transverse geometry. We denote the electric current density flowing in the $y$ direction generated near the Fℓ/N ($\ell=1,3$) interface by the inverse spin Hall effect as $j_{\ell y}$. The conservation law of the total electric current density, $j_{1y}+j_{3y}$, gives (see also Appendix B)

[TABLE]

In addition to the conventional spin torque, Eq. (11), this transverse electric current also excites the spin torque on the other ferromagnet. In total, the spin torque acting on the Fℓ layer is given by [$(\ell,\ell^{\prime})=(1,3)$ or $(3,1)$]

[TABLE]

The last two terms represent the coupling torque between the F1 and F3 layers. Note that the direction of this coupling torque is different from that of the conventional torque because the currents $J_{0}$ and $j_{\ell y}$ flow in different directions.

In the following, the torques related to $\chi$, $\chi^{\prime}$, and $\chi^{\prime\prime}$ in Eqs. (10) and (14) are referred to as coupling torque. The ratio of these new torques to the conventional spin torque is on the order of $\chi\propto\vartheta^{2}\sim 10^{-2}$. Since the conventional spin torque due is proportional to the spin Hall angle $\vartheta_{\rm R}\propto\vartheta$, the coupling torque is proportional to $\vartheta^{3}$. Although the spin Hall angle is usually a small quantity, it will be shown that the coupling torques play a non-negligible role on the magnetization dynamics, as shown below.

II.3 LLG equation

In the following sections, we study the magnetization dynamics excited by the spin torque given by Eq. (10) or (14) both numerically and analytically. We neglect the transverse coupling when the role of the longitudinal coupling is studied, and vice versa, for simplicity, which corresponds to considering a system consisting of the F1 and F2 layers, or F1 and F3 layers. The magnetization dynamics in the Fℓ ($\ell=1,2,3$) is described by the LLG equation,

[TABLE]

where the Gilbert damping constant is denoted as $\alpha$. In the following calculations, we use the values of the material parameters, $\gamma=1.764\times 10^{7}$ rad/(Oe s) and $\alpha=0.005$, derived from the experiments Tsunegi et al. (2014). For the later discussion, it is useful to show the explicit forms of the LLG equation in the longitudinal and transverse geometries. Equation (15) for the longitudinal geometry is

[TABLE]

whereas that for the transverse geometry is

[TABLE]

III Numerical analysis of synchronization

In this section, we study the magnetization dynamics in the ferromagnets in the presence of the coupling torques by solving Eq. (15) numerically. The self-oscillation of the magnetization provides an interesting example to understand the role of the coupling torque. Note that the coupling torques are proportional to the parameters $\chi$, $\chi^{\prime}$, and $\chi^{\prime\prime}$, and their products to other parameters such as $\alpha\beta\chi$, the values of which are relatively small, as can be seen in Eqs. (16) and (17). Nevertheless, the coupling torques lead to the phase synchronization, as shown below.

The self-oscillation of the magnetization in single ferromagnets by the spin Hall effect has been observed for in-plane magnetized ferromagnets Liu et al. (2012b), and is induced by the conventional spin torque given by Eq. (11). Therefore, in the following, we assume that the magnetic field,

[TABLE]

consists of an in-plane anisotropy field $H_{\rm K}$ along the $y$ direction and a demagnetization field $4\pi M$ in the $z$ direction, which we assume as $H_{\rm K}=200$ Oe and $M=1500$ emu/c.c. Kim et al. (2016) in the following calculations. It is known for the case of the single ferromagnet Hillebrands and Thiaville (2006) that the in-plane self-oscillation can be excited when the electric current density $J_{0}$ is in the range of $J_{\rm c}<J_{0}<J^{*}$, where

[TABLE]

[TABLE]

which in this study are $J_{\rm c}\simeq 26$ and $J^{*}\simeq 33$ MA/cm2.

III.1 Transverse geometry

Let us first investigate the magnetization dynamics in the transverse geometry because this geometry provides a simple example of the coupled motion. We start with solving the LLG equation given by Eq. (17) for the F1 and F3 layers. Figure 2(a) shows an example of the trajectory of the magnetization dynamics obtained by solving Eq. (17) numerically, where $J_{0}=30$ MA/cm2. As shown, the in-plane oscillation is observed in the steady state. The initial conditions set for $\mathbf{m}_{1}$ and $\mathbf{m}_{3}$ are different as $\mathbf{m}_{1}(0)=(\cos 80^{\circ},\sin 80^{\circ},0)$ and $\mathbf{m}_{3}(0)=(\cos 95^{\circ},\sin 95^{\circ},0)$. Therefore, the dynamics of $\mathbf{m}_{1}$ and $\mathbf{m}_{3}$ near the initial time are different, as shown in Fig. 2(b), where the time evolutions of $m_{1x}(t)$ and $m_{3x}(t)$ in $0\leq t\leq 1$ ns are shown. Nevertheless, the dynamics of $\mathbf{m}_{1}$ and $\mathbf{m}_{3}$ synchronize gradually, and finally, synchronization of $m_{1x}(t)=m_{3x}(t)$ and $m_{1y}(t)=m_{3y}(t)$ is realized, as shown in Figs. 2(c) and 2(d). We emphasize here that the dynamics shown in these figures are obtained for $\beta=-0.01$.

The mutual, as well as self, synchronization of the spin torque induced magnetization oscillation by using spin waves, electric current, microwave field, or dipole coupling has been an exciting topic from the viewpoints of both nonlinear science and practical applications Kaka et al. (2005); Mancoff et al. (2005); Grollier et al. (2006); Slavin and Tiberkevich (2006); Persson et al. (2007); Ruotolo et al. (2009); Georges et al. (2008); Zhou and Akerman (2009); Urazhdin et al. (2010); Sani et al. (2013); Demidov et al. (2014); Locatelli et al. (2015); Khalsa et al. (2015); Kendziorczyk and Kuhn (2016); Tsunegi et al. (2016); Awad et al. (2017); Kudo and Morie (2017). The key quantity of the synchronization is the phase difference $\Delta\varphi$ between each magnetization to enhance the emission power from the spin torque oscillators and to investigate new practical applications such as neuromorphic architectures Locatelli et al. (2014); Grollier et al. (2016). The synchronization found here, i.e., $m_{\ell x}(t)=m_{\ell^{\prime}x}(t)$, is called the in-phase synchronization. We should emphasize here that, although the results shown in Figs. 2(b) and 2(c) are shown for one certain initial condition, the in-phase synchronizations are confirmed for the present set of the parameters even when the initial conditions are changed.

On the other hand, it was shown for the case of the current-injection locking of the spin torque oscillator that the phase difference between the magnetization oscillation and the alternative current depends on the strength of the field-like torque Zhou and Akerman (2009). The field-like torque in the present system can be either positive or negative, as mentioned above. These facts motivate us to study the magnetization dynamics for different values of $\beta$. When $\beta=+0.01$, synchronized dynamics is observed in a similar manner, but in this case, the phase difference is antiphase, i.e., $m_{1x}(t)=-m_{3x}(t)$, as shown in Figs. 2(e) and 2(f).

Figure 3 summarizes the dependences of the phase difference, $\Delta\varphi$, between the magnetizations on the field-like torque strength $\beta$ for the several values of the damping constant $\alpha$. The vertical axis in this figure represents the phase difference in the unit of the oscillation period; i.e., $\Delta\varphi=0$ corresponds to the in-phase synchronization. On the other hand, $\Delta\varphi=0.50$ means that the phase difference is half of a period, and thus, is antiphase. The algorithm evaluating $\Delta\varphi$ in the numerical simulation is summarized in Appendix A. The results indicate that the phase difference is antiphase for positive $\beta$, whereas it becomes in-phase when $\beta$ becomes smaller than a certain value, except for the narrow intermediate region near $\beta\sim-\alpha/2$, where the phase difference is in between in-phase and antiphase.

III.2 Longitudinal geometry

Next, we study the magnetization dynamics in the longitudinal geometry between F1 and F2 layers by solving Eq. (16) numerically. Figure 4(a) and 4(b) show $m_{\ell x}(t)$ and $m_{\ell y}(t)$ in a steady state, where $\beta=-0.01$. The initial conditions in these figures are $\mathbf{m}_{1}(0)=(\cos 80^{\circ},\sin 80^{\circ},0)$ and $\mathbf{m}_{2}(0)=(\cos 95^{\circ},\sin 95^{\circ},0)$. An antiphase synchronization is observed in this case. We notice, however, that the in-phase synchronization can also be realized when the initial conditions are changed. Figure 4(c) and 4(d) show such an example, where $\mathbf{m}_{1}(0)=(\cos 80^{\circ},\sin 80^{\circ},0)$ and $\mathbf{m}_{2}(0)=(\cos 85^{\circ},\sin 85^{\circ},0)$ are assumed. These numerical results indicate that both the in-phase and antiphase synchronizations are stable in this case, and whether the phase difference finally becomes in-phase or antiphase depends on the initial conditions, material parameters, and current magnitude. We also notice that the phase difference is changed when the value of $\beta$ is changed. Figures 4(e) and 4(f) show $m_{\ell x}$ and $m_{\ell y}$ for $\beta=+0.01$. In this case, the phase difference of the magnetizations is a quarter of a precession period.

The dependences of the phase difference on the field-like torque strength $\beta$ for several values of the Gilbert damping constant $\alpha$ are summarized in Fig. 5, where $\Delta\varphi=0.25$ in this figure means that the phase difference is a quarter of a period. The phase difference is found to become a quarter of the period for positive $\beta$, whereas it becomes in-phase or antiphase for negative $\beta$, depending on the initial states of the magnetizations.

III.3 Summary of numerical simulations

Let us summarize the results of the numerical simulations here. In the transverse geometry, the coupling torque induces the synchronized oscillation of the magnetizations and finally stabilizes the configuration in the in-phase or antiphase state, depending on the values of the field-like torque strength $\beta$ and the damping constant $\alpha$. The phase difference in the stable synchronized state in the longitudinal geometry also depends on the values of $\beta$ and $\alpha$, as well as the initial conditions. In this case, however, in addition to the in-phase or antiphase state, a phase difference with a quarter of a period is generated.

The phase difference can be measured from the spin Hall magnetoresistance effect. According to Eq. (13), the total electric current density in the transverse direction is

[TABLE]

Then, an oscillating current appears in the transverse direction for in-phase synchronization, whereas the transverse current becomes zero for antiphase synchronization, as shown in Fig. 6(a). The Fourier transformation of Eq. (21) for in-phase synchronization has peaks at the frequencies of $f_{n}=(2n-1)f_{0}$ ($n=1,2,3,\cdots$), where $f_{0}$ is the lowest frequency; see Appendix C. Similarly, the total electric current density in the longitudinal direction is

[TABLE]

The oscillation frequency of this current is different for the synchronizations having the phase difference of in-phase or antiphase and that of a quarter of a period, as shown in Fig. 6(b). The Fourier transformation of Eq.(22) has the peaks at $f_{n}=2nf_{0}$ for in-phase or antiphase and $f_{n}=4nf_{0}$ when the phase difference is a quarter of a period.

An interesting question regarding these numerical results is to clarify the reason why the phase difference finally becomes a certain value for a given set of the parameters, i.e., which phase difference is an attractor of the limit cycle. It is, however, difficult to answer this question directly due to the following reason. We note that the present model includes several small-valued parameters, $\alpha$, $\beta$, $\chi$, $\chi^{\prime}$, and $\chi^{\prime\prime}$, as shown in Eqs. (16) and (17), and is complicated. The torque related to the lowest order of $\beta$ in these equations is the conventional field-like torque given by $[\gamma\beta\hbar\vartheta_{\rm R}J_{0}/(2eMd_{\rm F})]\mathbf{m}_{\ell}\times\mathbf{e}_{y}$. It should be noted that the phase difference is not determined solely by this term because this lowest order term of the field-like torque does not include the coupling between the magnetizations. Similarly, the attractor is not determined solely by the lowest order term of the coupling torque, which is $\left[\gamma\hbar\vartheta_{\rm R}J_{0}/(2eMd_{\rm F})\right]\chi m_{\ell^{\prime}y}^{2}\mathbf{m}_{\ell}\times\left(\mathbf{e}_{y}\times\mathbf{m}_{\ell}\right)$ in Eq. (16) and $\left[\gamma\hbar\vartheta_{\rm R}J_{0}/(2eMd_{\rm F})\right]\chi m_{\ell^{\prime}x}m_{\ell^{\prime}y}\mathbf{m}_{\ell}\times\left(\mathbf{e}_{x}\times\mathbf{m}_{\ell}\right)$ in Eq. (17), because this torque does not include the field-like torque. Their combinations or the higher order terms including both $\beta$ and the coupling torques related to $\chi$, $\chi^{\prime}$, and $\chi^{\prime\prime}$ should be taken into account to answer the question, which is difficult due to the nonlinearity and complexity of the LLG equation.

Nevertheless, it is possible to reveal the relation between the current and frequency in the synchronized oscillation state by assuming a certain value of the phase difference between the magnetizations. The current-frequency relation has been often measured in the experiments, and therefore, it will be useful to develop a theory clarifying the role of the coupling on the current-frequency relation. In the next section, we will investigate this subject.

IV Theoretical analysis of current-frequency relation

The purpose of this section is to develop an analytical theory of the synchronization revealing the relation among the current, frequency, and the phase difference of the magnetizations in the synchronized oscillation state.

IV.1 Basis of analysis

Here, let us mention the basis of our theoretical analysis. It is difficult to solve the LLG equation exactly because of its nonlinearity. Instead, we employ the averaging technique of the LLG equation on constant energy curves Bertotti et al. (2009). This approach has been used to study the spin torque switching in thermally activated regions Apalkov and Visscher (2005); Bertotti et al. (2006); Newhall and Eijnden (2013); Taniguchi et al. (2013); Pinna et al. (2013) and spin torque oscillators Bertotti et al. (2005); Serpico et al. (2005); Ebels et al. (2008); Bazaliy and Arammash (2011); Dykman (2012); Lacoste et al. (2013); Taniguchi (2014a); Pinna et al. (2014); Taniguchi (2015), as well as the microwave assisted magnetization reversal Taniguchi (2014b); Suto et al. (2015), but has not been applied to the coupled system. This approach is valid when the magnetic energy is changed slowly compared with the oscillation period. We note that only the lowest order terms in the LLG equation is necessary to derive the current-frequency relation, as far as several parameters such as $\beta$ are small. Thus, we use the simplified LLG equation maintaining only the dominant terms. The LLG equation used in this section for the longitudinal geometry is

[TABLE]

whereas that for the transverse geometry is

[TABLE]

In the self-oscillation state, the damping torque during the precession is balanced with the spin torque, and the torque due to the magnetic field, corresponding to the first term on the right hand side of Eq. (15), becomes the dominant term determining the magnetization dynamics. The dynamic trajectory given by this field torque corresponds to constant energy curves of the energy density, $E=-M\int d\mathbf{m}_{\ell}\cdot\mathbf{H}_{\ell}$, where its explicit form is

[TABLE]

The minimum and saddle points of Eq. (25) are $\mathbf{m}_{\rm min}=\pm\mathbf{e}_{y}$ and $\mathbf{m}_{\rm saddle}=\pm\mathbf{e}_{x}$, where the corresponding energy densities are $E_{\rm min}=-MH_{\rm K}/2$ and $E_{\rm saddle}=0$. The solution of $\mathbf{m}_{\ell}$ precessing on a constant energy curve is described by the Jacobi elliptic function Byrd and Friedman (1971) as (see also Appendix C summarizing the derivations)

[TABLE]

[TABLE]

[TABLE]

where $\varphi_{0}$ is the initial phase. The period of $m_{x}$ and $m_{z}$ is $\tau(E)$, whereas that of $m_{y}$ is $\tau(E)/2$, where the precession frequency $f(E)=1/\tau(E)$ is given by

[TABLE]

where $\mathsf{K}(k)$ is the first kind of complete elliptic integral with the modulus

[TABLE]

Note that Eq. (29) reproduces the ferromagnetic resonance (FMR) frequency, $\gamma\sqrt{H_{\rm K}\left(H_{\rm K}+4\pi M\right)}/(2\pi)$, in the limit of $E\to E_{\rm min}$. Identifying $E$ and $\varphi_{0}$ corresponds to the determination of the initial condition.

The averaged technique investigates the energy change during a precession on a constant energy curve, which is obtained from the LLG equation as

[TABLE]

where $\mathscr{W}_{\rm s}$ is the energy change by the conventional spin torque due to the spin Hall effect, whereas $\mathscr{W}_{\alpha}$ is the dissipation due to the damping torque. The integral range is over the precession period. The explicit forms of $\mathscr{W}_{\rm s}$ and $\mathscr{W}_{\alpha}$ are given by Taniguchi et al. (2013)

[TABLE]

[TABLE]

where $\mathsf{E}(k)$ is the second kind of complete elliptic integral.

On the other hand, $\mathscr{W}_{\rm s}^{\rm L,T}$ represents the work done by the coupling torque in the longitudinal or transverse geometry, corresponding to the last term in Eq. (23) or (24). The explicit forms of $\mathscr{W}_{\rm s}^{\rm L}$ and $\mathscr{W}_{\rm s}^{\rm T}$ are shown in the following sections. For both cases, the relation between the current and frequency is clarified as follows. In the self-oscillation state, since the spin torque balances the damping torque, the following condition should be satisfied,

[TABLE]

The current density $J_{0}$ satisfying Eq. (34) is the current necessary to excite the self-oscillation on the constant energy curve of $E$, and is denoted as $J_{0}(E)$. The relation between the current and frequency in the self-oscillation state is given by this $J_{0}(E)$ and $f(E)$ given by Eq. (29). It should be emphasized that the current density $J_{0}(E)$ depends on the phase difference between the magnetizations through $\mathscr{W}_{\rm s}^{\rm L,T}$. We will therefore study the relation between the phase difference and the current-frequency relation in line with this deduction.

IV.2 Transverse geometry

In this section, we investigate the current-frequency relation in the transverse geometry. The work done by the coupling torque is defined as

[TABLE]

Before advancing the calculation, let us briefly mention the definition of the phase difference in the present approach. If the oscillation trajectory is a circle, the phase difference is easily defined, i.e., the antiphase corresponds to $\Delta\varphi=\pi$, whereas $\Delta\varphi=0$ is the in-phase. In the present case, on the other hand, the oscillation trajectory is described by the elliptic function, as shown in Eqs. (26), (27), and (28). In this case, the phase difference is defined using the elliptic integral $\mathsf{K}(k)$, where $\Delta\varphi=0$ for the in-phase synchronization, and the antiphase synchronization corresponds to $\Delta\varphi=2\mathsf{K}(k)$. It is useful to note that $\Delta\varphi=2\mathsf{K}(k)$ becomes $\pi$ in the limit of $k\to 0$, corresponding to the case that the oscillation trajectory becomes a circle. Similarly, $\Delta\varphi=\mathsf{K}(k)$ means that the phase difference is a quarter of a period.

Equation (35) for an arbitrary phase difference is evaluated by numerically calculating the integral with the solutions of $\mathbf{m}_{\ell}$ and $\mathbf{m}_{\ell^{\prime}}$ shown in Appendix D. It is, however, useful to derive the analytical solutions of Eq. (35) for specific values of the phase difference. Equation (35) for both the in-phase ($\Delta\varphi=0$) and antiphase [$\Delta\varphi=2\mathsf{K}(k)$] becomes

[TABLE]

where the double sign means the upper for the in-phase ($\Delta\varphi=0$) synchronization and the lower for the antiphase [$\Delta\varphi=2\mathsf{K}(k)$] synchronization. Equation (36) is zero at $E=E_{\rm min}$ and $E_{\rm saddle}$, and is negative (positive) for the energy density $E$ in the rage of $E_{\rm min}<E<E_{\rm saddle}$ when $\Delta\varphi=0$ [$2\mathsf{K}(k)$]. This means that the coupling torque acts as a damping (an antidamping) torque when the phase difference is in-phase (antiphase). We also note that $\mathscr{W}_{\rm s}^{\rm T}=0$ when $\Delta\varphi=\mathsf{K}(k)$; i.e., the phase difference is a quarter of a period. The calculations necessary to obtain these specific values of $\mathscr{W}_{\rm s}^{\rm T}$ are also summarized in Appendix D. We note that the sign change of $\mathscr{W}_{\rm s}^{\rm T}$ with respect to the phase difference is related to the fact that the coupling torque in the transverse geometry, Eq. (14), has the angular dependence of $m_{\ell^{\prime}x}m_{\ell^{\prime}y}\mathbf{m}_{\ell}\times\left(\mathbf{e}_{x}\times\mathbf{m}_{\ell}\right)$. Because of this angular dependence, the coupling torque acts as an anti-damping (a damping) torque when $m_{\ell x}$ and $m_{\ell^{\prime}x}$ have the opposite (same) signs, resulting in the increase (decrease) of the energy supplied to the ferromagnets by the coupling torque.

In summary, the work done by the coupling torque, $\mathscr{W}_{\rm s}^{\rm T}$, is negative and minimized at the in-phase ($\Delta\varphi=0$), zero for $\Delta\varphi=\mathsf{K}(k)$, and positive and maximized at the antiphase [$\Delta\varphi=2\mathsf{K}(k)$].

The current $J_{0}$ in the transverse geometry is defined as

[TABLE]

where $\mathcal{N}$ and $\mathcal{D}_{\rm T}$ are defined as

[TABLE]

[TABLE]

The explicit form of $\mathcal{N}$ is obtained from Eq. (33) as

[TABLE]

On the other hand, $\mathcal{D}_{\rm T}$ for the in-phase or antiphase is obtained from Eqs. (32) and (36) as

[TABLE]

where the double sign means the upper for the in-phase synchronization and the lower for the antiphase synchronization.

Figure 7(a) shows $J_{0}(E)$ as functions of the oscillation frequency $f(E)$ and the phase difference $\Delta\varphi$. The current-frequency relation in the transverse geometry can be obtained from this figure. To reveal the role of the phase difference more clearly, we show $J_{0}(E)$ for a certain value of $f(E)(=4.6{\rm GHz})$ in Fig. 7(b). Note that $J_{0}(E)$ is smaller than that in the absence of the coupling, which is shown by the dotted line, and minimized when $\Delta\varphi=2\mathsf{K}(k)$, i.e., the antiphase. This is because the work done by the coupling torque is positive and maximized at the antiphase. On the other hand, $J_{0}(E)$ is maximized at the in-phase, and is larger than that in the absence of the coupling because the work done by the coupling torque is negative and minimized at the in-phase. We notice that the phase differences observed in the numerical simulation in Sec. III, i.e., the in-phase and antiphase, correspond to $\Delta\varphi$ satisfying

[TABLE]

or equivalently,

[TABLE]

In other words, the phase differences observed in the numerical simulations correspond to those giving the extrema of $J_{0}$ ($\mathscr{W}_{\rm s}^{\rm T}$).

IV.3 Longitudinal geometry

Let us investigate the theoretical relation between the current and frequency in the longitudinal geometry. In this case, the averaged LLG equation is given by

[TABLE]

where $\mathscr{W}_{\rm s}$ and $\mathscr{W}_{\alpha}$ are given by Eqs. (32) and (33). On the other hand, $\mathscr{W}_{\rm s}^{\rm L}$ representing the energy change due to the longitudinal coupling is defined as

[TABLE]

For both the in-phase and antiphase, Eq. (45) becomes (see Appendix D)

[TABLE]

On the other hand, when the phase difference is a quarter of a period [$\Delta\varphi=\mathsf{K}(k)$], Eq. (45) becomes

[TABLE]

It should be emphasized that $\mathscr{W}_{\rm s}^{\rm L}$ is always positive for an arbitrary phase difference. This is because the coupling torque in Eq. (10) always acts as an anti-damping torque.

We can determine the current density $J_{0}(E)$ satisfying Eq. (34) in the longitudinal geometry, as in the case of the transverse geometry, by replacing $\mathcal{D}_{\rm T}$ in Eq. (37) with

[TABLE]

The explicit form of $\mathcal{D}_{\rm L}$ for both the in-phase and antiphase is obtained from Eqs. (32) and (46) as

[TABLE]

whereas that when the phase difference is a quarter of a period is obtained from Eqs. (32) and (47) as

[TABLE]

Figure 7(c) shows $J_{0}(E)$ as functions of $f(E)$ and $\Delta\varphi$. The current-frequency relation in the longitudinal geometry can be obtained from this figure. It is noted that $J_{0}(E)$ is always smaller than that in the absence of the coupling because the coupling torque in the longitudinal geometry always points to the anti-damping direction, and therefore, the work done by the coupling torque is always positive. Figure 7(b) shows $J_{0}(E)$ as a function of $\Delta\varphi$ at a certain value of $f(E)$. As shown, $J_{0}(E)$ has minima at both the in-phase ($\Delta\varphi=0$) and the antiphase [$\Delta\varphi=2\mathsf{K}(k)$], whereas it is maximized when the phase difference is a quarter of a period [$\Delta\varphi=\mathsf{K}(k)$]. We again notice that these phase differences found in the numerical simulations in Sec. III correspond to $\Delta\varphi$ satisfying $\partial J_{0}(E)/\partial\Delta\varphi=0$, or equivalently,

[TABLE]

IV.4 Phase differences in stable synchronization and fixed points of effective potential

Equation (31) describes the slow change of the magnetic energy in the oscillation state. The magnetization dynamics is regarded as a motion of a point particle in an effective potential given by its right-hand side. Equations (43) and (51) correspond to the stability conditions of the point particle in this effective potential. Therefore, the phase difference found in the numerical simulation finally converges to one of these values satisfying Eq. (43) or Eq. (51). Whether the in-phase, antiphase, or the phase difference with a quarter of a period becomes the attractor depends on the higher order terms of the small parameters, as well as the initial states of the magnetizations, as mentioned at the end of Sec. III. This discussion is beyond the scope of this paper.

IV.5 Instability threshold

At mentioned at the beginning of Sec. III, the in-plane self-oscillation for a single ferromagnet is stabilized when the current density is in the range of $J_{\rm c}<J_{0}<J^{*}$, where $J_{\rm c}$ and $J^{*}$ are given by Eqs. (19) and (20), respectively. At the end of this section, let us briefly discuss the effect of the coupling on these scaling currents.

Let us remind the readers that $J_{\rm c}$ is the current density necessary to destabilize the magnetization in equilibrium, whereas $J^{*}$ is the current necessary to overcome the energy barrier, $E_{\rm saddle}-E_{\rm min}$. These current densities are theoretically defined as Taniguchi et al. (2013)

[TABLE]

[TABLE]

It is confirmed that Eqs. (19) and (20) are reproduced by substituting Eqs. (32) and (33) in the definition of $J_{0}(E)$ in the absence of the coupling.

On the other hand, in the presence of the transverse coupling, it is confirmed from Eq. (36) that a factor $\left[1-\left(\chi/2\right)\right]$ should be multiplied to the denominator of Eq. (19) when the phase difference between the magnetizations is in-phase, whereas this factor is replaced by $\left[1+\left(\chi/2\right)\right]$ when the phase difference is antiphase. The other scaling current, $J^{*}$, is unchanged for these phase differences. In the longitudinal geometry, we see from Eqs. (46) and (47) that the factor $\left(1+\chi+\chi^{\prime\prime}\right)$ should be multiplied to the denominator of Eq. (19) when the phase difference is in-phase, antiphase, or a quarter of a period, whereas, for $J^{*}$, the factor becomes $1+\chi^{\prime\prime}+\left(\chi/2\right)\left[4\pi M/\left(H_{\rm K}+4\pi M\right)\right]$ for in-phase and antiphase, and $1+\chi^{\prime\prime}$ when the phase difference is a quarter of a period.

V Conclusion

In conclusion, the coupled magnetization dynamics in the ferromagnets through the spin Hall magnetoresistance effect was investigated. The coupling appears in both the longitudinal and transverse directions of the alignment of the ferromagnets. The in-phase or antiphase synchronization of the magnetization oscillation was found in the transverse geometry by solving the LLG equation numerically. On the other hand, in addition to them, the synchronization having the phase difference of a quarter of a period is also found in the longitudinal geometry. It was shown that these phase differences depend on the values of the damping constant and the field-like torque strength. The analytical theory revealing the relation among the current, frequency, and phase difference was also developed. It was shown that the phase differences observed in the numerical simulations correspond to that giving the fixed points of the energy supplied by the coupling torque.

Acknowledgement

The author is grateful to Takehiko Yorozu for his contributions to the analytical calculations, and to Hitoshi Kubota, Sumito Tsunegi, Yoichi Shiota, Shingo Tamaru, Tazumi Nagasawa, Kiwamu Kudo, and Yoji Kawamura for valuable discussions. The author is also thankful to Satoshi Iba, Aurelie Spiesser, Atsushi Sugihara, Takahide Kubota, Hiroki Maehara, and Ai Emura for their support and encouragement. This work was supported by JSPS KAKENHI Grant-in-Aid for Young Scientists (B) 16K17486.

Appendix A Values of parameters in numerical simulations

The exact values of the parameters used in the simulations, evaluated from the parameters found in the experiment Kim et al. (2016), are $\vartheta_{\rm R}=0.16680863$, $\beta=-0.00973617$, $\chi=0.009525272$, $\chi^{\prime}=-0.000152995$, and $\chi^{\prime\prime}=0.03516089$ for $g_{\rm i}^{\uparrow\downarrow}/S=1.0$ nm*-2*. In the main text, $\beta=-0.01$ and $\beta=+0.01$ correspond to $\beta=-0.00973617$ and $\beta=+0.00973617$, respectively. Strictly speaking, the change of the value of $g_{\rm i}^{\uparrow\downarrow}$ affects not only $\beta$ but also other quantities such as $\vartheta_{\rm R}$, $\chi$, and $\chi^{\prime}$. We, however, change the value of $\beta$ only in the numerical simulation, for simplicity, because the results do not change significantly unless $|g_{\rm i}^{\uparrow\downarrow}/g_{\rm r}^{\uparrow\downarrow}|\ll 1$. The LLG equation with these parameters is solved by using the fourth-order Runge-Kutta method from $t=0$ to $t=1$ $\mu$s with the time step of $\Delta t=20$ fs; i.e., the number of the time mesh is $N_{\rm t}=5\times 10^{7}$.

The present system has two stable states at $\mathbf{m}_{\ell}=\pm\mathbf{e}_{y}$. For convention, we assume that the magnetizations initially stay near one equilibrium, $\mathbf{m}_{\ell}=+\mathbf{e}_{y}$. For in-phase synchronizations, such as those shown in Figs. 2(c) and 2(d), the $z$ components are also synchronized with in-phase, i.e., $m_{\ell z}(t)=m_{\ell^{\prime}z}(t)$. On the other hand, for antiphase synchronization shown in, for example in Figs. 2(e) and 2(f), the $z$ components are also synchronized with antiphase, $m_{\ell z}(t)=-m_{\ell^{\prime}z}(t)$.

The algorithm evaluating the phase differences shown in Figs. 3 and 5 from the discrete numerical data is as follows. We gathered $N_{\rm i}=2^{16}=65536$ data of $\mathbf{m}_{\ell}(t)$ ($\ell=1,2,3$) from $t=\left(N_{\rm t}-N_{\rm i}+1\right)\Delta t$ to $t=N_{\rm t}\Delta t=1$ $\mu$s. Then, the averaged periods $T_{\ell}$ of the oscillation of each magnetization were evaluated from the peaks of $\mathbf{m}_{\ell}(t)$ in this time range as $T_{\ell}=\left[\left(t_{\ell,a}-t_{\ell,a-1}\right)+\cdots+\left(t_{\ell,2}-t_{\ell,1}\right)\right]/\left(N_{\ell}-1\right)=\left(t_{\ell,a}-t_{\ell,1}\right)/\left(N_{\ell}-1\right)$, where $N_{\ell}$ is the number of the peaks in $\mathbf{m}_{\ell}(t)$, whereas $t_{\ell,a}$ is the time corresponding to the $a$-th peak. Then, the phase difference is evaluated as $\Delta\varphi=\sum_{a=1}^{N^{\prime}}|t_{\ell,a}-t_{\ell^{\prime},a}|/\left(N^{\prime}\bar{T}\right)$, where $N^{\prime}={\rm min}\left[N_{\ell},N_{\ell^{\prime}}\right]$ and $\bar{T}=\left(T_{\ell}+T_{\ell^{\prime}}\right)/2$ with $(\ell,\ell^{\prime})=(1,3)$ for the transverse geometry, whereas that is $(1,2)$ for the longitudinal geometry. For the in-phase synchronization, this $\Delta\varphi$ is zero because $t_{\ell,a}=t_{\ell^{\prime},a}$. When the phase difference is antiphase, $\Delta\varphi=0.50$ because $|t_{\ell,a}-t_{\ell^{\prime},a}|=\bar{T}/2$ in this case. Similarly, $\Delta\varphi$ is $0.25$ when the phase difference is a quarter of a period.

Note that the critical current density to excite the self-oscillation, given by Eq. (19), is proportional to the damping constant $\alpha$. Therefore, the value of the current density should be increased to observe the self-oscillation when $\alpha$ is varied, as in the case of Figs. 3 and 5. In these figures, $J_{0}$ is assumed as $n\times 30$ MA/cm2 for $\alpha=n\times 0.005$ ($n=1-6$).

The numerical simulation in Fig. 2(e) indicates that the antiphase synchronization is an attractor for $\beta=+0.01$. An exception is that if the initial conditions are set to be identical, the final state becomes in-phase synchronization due to the symmetry of the LLG equation with respect to the change of $(\ell,\ell^{\prime})\to(\ell^{\prime},\ell)$. Since Eq. (43) is satisfied, the phase difference is fixed to in-phase even if it is unstable. Similar situations occur in other cases for such specific initial conditions.

Appendix B Derivation of coupling torque in transverse geometry

In a ferromagnetic/nonmagnetic bilayer, the spin current density flowing in the $i$ direction ($i=x,y,z$) with the spin polarization in the $\nu$ direction is related to the electrochemical potential $\bar{\mu}_{\rm N}$ and the spin accumulation $\delta\mu_{\rm N}$ via

[TABLE]

where $\partial_{j}\bar{\mu}_{\rm N}/e$ is the electric field in the $j$ direction, and therefore, $\sigma_{\rm N}\partial_{j}\bar{\mu}_{\rm N}/e$ is the electric current density. We assume that this equation is extended to

[TABLE]

in the transverse geometry, where $J_{{\rm s}z\nu,{\rm N}}^{(\ell)}$ is the spin current density flowing near the Fℓ/N interface in the $z$ direction with the spin polarization in the $\nu$ direction. The spin accumulation obeys the diffusion equation, and the boundary conditions of the diffusion equation are given by the spin current density at the boundaries. Using Eq. (55), the solution of the spin accumulation is given by Taniguchi (2016)

[TABLE]

where we assume that the nonmagnet is in the region of $-d_{\rm N}\leq z\leq 0$. The spin current density, $J_{{\rm s}z\nu}^{{\rm F}_{\ell}/{\rm N}}$, at the Fℓ/N interface is given by Y.-T.Chen et al. (2013); Taniguchi (2016)

[TABLE]

where the vector notation in boldface represents the direction of the spin polarization, whereas the spatial direction of Eq. (57) is defined as the positive direction, i.e., from the nonmagnet to the ferromagnet.

On the other hand, the electric current density in the nonmagnet flowing in the $i$ direction is given by

[TABLE]

In the present case, the electric current density near the Fℓ/N interface flowing in the $y$ direction becomes

[TABLE]

Substituting Eqs. (56) and (57) into Eq. (59), and averaging along the $z$ direction as $\overline{J_{{\rm c}y,{\rm N}}^{(\ell)}}=(1/d_{\rm N})\int_{-d_{\rm N}}^{0}J_{{\rm c}y,{\rm N}}^{(\ell)}dz$, we find that

[TABLE]

The conservation law of the electric current along the $y$ direction requires that $\overline{J_{{\rm c}y,{\rm N}}^{(\ell)}}=j_{\ell y}+j_{\ell^{\prime}y}$. Solving this equation for $\ell=1$ and $3$, we obtain Eq. (13). The spin torque is defined from Eq. (57) as

[TABLE]

Appendix C Analytical solution of magnetization on a constant energy curve

In this appendix, we shows the derivation of the analytical solution of the magnetization on a constant energy curve. For simplicity, we remove the subscript $\ell(=1,2,3)$ distinguishing the ferromagnets. The magnetization dynamics on a constant energy curve is described by the Landau-Lifshitz (LL) equation

[TABLE]

The magnetic field, Eq. (18), is related to the magnetic energy density $E$ via $E=-M\int d\mathbf{m}\cdot\mathbf{H}$, as mentioned in the main text. Using the relation $m_{x}^{2}+m_{y}^{2}+m_{z}^{2}=1$, Eq. (25) is rewritten as

[TABLE]

This equation indicates that $m_{z}$ and $m_{x}$ can be expressed as $m_{z}=v^{\prime}\cos u$ and $m_{x}=\left(v^{\prime}/v\right)\sin u$, respectively, where $v$ and $v^{\prime}$ are defined as $v^{2}=H_{\rm K}/\left(H_{\rm K}+4\pi M\right)$ and $v^{\prime}=\left(2E/M+H_{\rm K}\right)/\left(H_{\rm K}+4\pi M\right)$. Then, $du/dt=\left(du/d\sin u\right)\left(d\sin u/dt\right)=\left(1/\cos u\right)\left[d\left(v/v^{\prime}\right)m_{x}/dt\right]=\left(v/m_{z}\right)\left(dm_{x}/dt\right)$, which becomes, from Eq. (62),

[TABLE]

Introducing new variable $w=\sin u$, this equation gives

[TABLE]

where the modulus $k$ is given by Eq. (30). The modulus monotonically varies from [math] to $1$ by changing the energy density $E$ from its minimum $E_{\rm min}$ to saddle $E_{\rm saddle}$. We also notice that $\left(H_{\rm K}+4\pi M\right)v\sqrt{1-v^{\prime 2}}=\sqrt{H_{\rm K}\left(4\pi M-2E/M\right)}$. Equation (65) indicates that $w$ is given by

[TABLE]

where ${\rm sn}(u,k)$ is the Jacobi elliptic function, and $\varphi_{0}$ is the initial phase determined by the initial condition. Using the relations ${\rm sn}^{2}(u,k)+{\rm cn}^{2}(u,k)=1$ and ${\rm dn}^{2}(u,k)=\sqrt{1-k^{2}{\rm sn}^{2}(u,k)}$, we find that the solution of $\mathbf{m}$ on the constant energy curve is given by Eqs. (26), (27), and (28).

The peak frequencies of the Fourier transformation of Eq. (21) in the transverse geometry are discussed as follows. We note that Eq. (21) for in-phase synchronization is proportional to $m_{\ell x}(t)m_{\ell y}(t)$ and $m_{\ell z}(t)$. Substituting the following formulas Byrd and Friedman (1971),

[TABLE]

[TABLE]

[TABLE]

to Eqs. (26), (27), (28), where $q=\exp\left[-\pi\mathsf{K}\left(\sqrt{1-k^{2}}\right)/\mathsf{K}(k)\right]$, it is found that the peak frequencies of Eq. (21) appear at $f_{n}=(2n-1)f_{0}$ ($n=1,2,3,\cdots$), where the lowest frequency $f_{0}$ is given by Eq. (29).

In the longitudinal geometry, Eq. (22) is proportional to $m_{1y}^{2}+m_{2y}^{2}$. When the phase difference of the magnetizations is in-phase or antiphase, it becomes $2m_{1y}^{2}$. In this case, using the formula Kiper (1984)

[TABLE]

it is found that the Fourier transformation of Eq. (22) has the peaks at $f_{n}=2nf_{0}$. On the other hand, when the phase difference is a quarter of a period, Eq. (22) is proportional to $g(u)\equiv{\rm dn}^{2}(u,k)+{\rm dn}\left[u+\mathsf{K}(k),k\right]={\rm dn}^{2}(u,k)+\left[\left(1-k^{2}\right)/{\rm dn}^{2}(u,k)\right]$. We notice that $g\left[u+\mathsf{K}(k)\right]=g(u)$, indicating that the Fourier transformation of Eq. (22) in this case has the peaks at $f_{n}=4nf_{0}$.

Appendix D Details of calculations of Eqs. (35) and (45)

Equations (35) and (45) can be calculated by substituting the solution of $\mathbf{m}_{\ell}$ on a constant energy curve to the integrals. As emphasized in the main text, the phase difference $\Delta\varphi$ between the magnetizations is an important quantity. According to Eqs. (26), (27), and (28), we set $\mathbf{m}_{\ell}$ and $\mathbf{m}_{\ell^{\prime}}$ as

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

The value of $\Delta\varphi$ varies in the rage of $0\leq\Delta\varphi<4\mathsf{K}(k)$. $\Delta\varphi=0$ corresponds to the in-phase synchronization, whereas $\Delta\varphi$ is $2\mathsf{K}(k)$ for the antiphase synchronization.

The analytical formulas of Eqs. (35) and (45) for the in-phase and antiphase synchronizations can be obtained as follows. First, since the elliptic functions satisfy ${\rm sn}\left[u+2\mathsf{K}(k),k\right]=-{\rm sn}(u,k)$, ${\rm cn}\left[u+2\mathsf{K}(k),k\right]=-{\rm cn}(u,k)$, and ${\rm dn}\left[u+2\mathsf{K}(k),k\right]={\rm dn}(u,k)$, $\mathscr{W}_{\rm s}^{\rm T}$ has the same magnitude but different sign for $\Delta\varphi=0$ and $\Delta\varphi=2\mathsf{K}(k)$, whereas $\mathscr{W}_{\rm s}^{\rm L}$ is the same for the in-phase and antiphase. Therefore, it is sufficient to calculate $\mathscr{W}_{\rm s}^{\rm T}$ and $\mathscr{W}_{\rm s}^{\rm L}$ for the in-phase case. In this case, it is unnecessary to distinguish $\mathbf{m}_{\ell}$ and $\mathbf{m}_{\ell^{\prime}}$. Next, it should be noted that Eq. (35) includes the following two integrals,

[TABLE]

[TABLE]

By replacing the integral variable from $u$ with $x={\rm sn}(u,k)$, and noting that $du=dx/\sqrt{\left(1-x^{2}\right)\left(1-k^{2}x^{2}\right)}$, these integrals are calculated as

[TABLE]

[TABLE]

Using these integrals, Eq. (36) is obtained. On the other hand, Eq. (45) includes the following three integrals,

[TABLE]

[TABLE]

[TABLE]

Using these integrals, Eq. (46) is obtained.

When the phase difference is a quarter of a period ($\Delta\varphi=\mathsf{K}(k)$), the relations, ${\rm sn}\left[u+\mathsf{K}(k),k\right]={\rm cn}(u,k)/{\rm dn}(u,k)$, ${\rm cn}\left[u+\mathsf{K}(k),k\right]=-\sqrt{1-k^{2}}{\rm cn}(u,k)/{\rm dn}(u,k)$, ${\rm dn}\left[u+\mathsf{K}(k),k\right]=\sqrt{1-k^{2}}/{\rm dn}(u,k)$, are used to evaluate Eqs. (35) and (45). In this case, we notice that all the integrands in Eq. (35) become odd functions of $t$, and the integrals over $0\leq t\leq\tau$ are zero. Therefore, $\mathscr{W}_{\rm s}^{\rm T}=0$ for $\Delta\varphi=\mathsf{K}(k)$. On the other hand, the following integrals are necessary to obtain Eq. (47),

[TABLE]

[TABLE]

[TABLE]