\section{Introduction}

Over the past decades, the synthesis of superheavy nuclei (SHN) attracted much attention and has been obtained progress in experiments via massive fusion reactions. The seventh period in the periodic table was filled with the superheavy element tennessine (Ts) by using the Dubna Gas-Filled Recoil Separator (DGFRS) at the Flerov Laboratory of Nuclear Reactions in Dubna, Russia \cite{Og10}. The existence of superheavy elements (SHE) was predicted in the late 1960s by the macroscopic-microscopic theory of the atomic nucleus \cite{So66}. The synthesis of SHN is associated with testing the shell model beyond the doubly magic nucleus $^{208}$Pb, hunting for the 'island of stability', exploring the limit of the mass of atomic nucleus, providing a strong Coulomb field such as quantum electrodynamics (QED) in the super-strong electric field etc. The existence of the superheavy nucleus (SHN) ($Z\geq106$) is due to strong binding shell effect against the Coulomb repulsion. Therefore, the position of shell closure is particularly significant for the properties of SHN, i.e., half-lives of $\alpha$ decay chain and spontaneous fission, formation probability etc, which is predicted at Z=114 and N=184 by theoretical models \cite{So07,Cw05}. Attempts to synthesize elements beyond Og (Z=118) were performed via different systems, e.g., $^{64}$Ni+$^{238}$U \cite{Ho08}, $^{58}$Fe+$^{244}$Pu \cite{Og09}, $^{54}$Cr+$^{248}$Cm \cite{Ho12,Ho16} and $^{50}$Ti+$^{249}$Cf \cite{Kh13}. The mass and angle distributions of fission fragments were measured \cite{Al20}. Systematical analysis of different reactions needs to be performed for preferentially producing new SHN in experiments.

The synthesis of SHN goes back 40 years ago with the multi-nucleon transfer reactions in collisions of two actinide nuclei \cite{Hu77,Sc78}. However, the yields of the heavy fragments in strongly damped collisions were found to decrease very rapidly with increasing the atomic number and to be impossible for producing SHN because of the very low cross section. Combinations with a doubly magic nucleus or nearly magic nucleus are usually chosen owing to the larger reaction $Q$ values. Reactions with $^{208}$Pb or $^{209}$Bi based targets were firstly proposed by Oganessian et al. \cite{Og75}. The SHEs from Bh to Cn were synthesized in the cold fusion reactions at GSI (Darmstadt, Germany) with the heavy-ion accelerator UNILAC and the SHIP separator \cite{Ho00,Mu15}. Experiments on the synthesis of element Nh (Z=113) in the $^{70}$Zn+$^{209}$Bi reaction have been performed successfully at RIKEN (Tokyo, Japan) \cite{Mo04}. However, it is very difficulty to create the superheavy isotopes beyond Nh in the cold fusion reactions because of very low cross section ($\sigma<$0.1 pb). The superheavy elements from Fl (Z=114) to Og (Z=118) have been synthesized at the Flerov Laboratory of Nuclear Reactions (FLNR) in Dubna (Russia) with the double magic nuclide $^{48}$Ca bombarding actinide nuclei \cite{Og99,Og06,Og15}, in which more neutron-rich SHN were produced and identified by the subsequent $\alpha$-decay chain. The decay properties of $^{271}$Ds in the cold fusion reaction of $^{64}$Ni+$^{208}$Pb$\rightarrow ^{271}$Ds+n were identified by a gas-filled recoil separator at Institute of Modern Physics (IMP) in Lanzhou \cite{Zh12}. With constructing the new facilities in the world such as RIBF (RIKEN, Japan), SPIRAL2 (GANIL in Caen, France), FRIB (MSU, USA), HIAF (IMP, China), it would be possible to create SHNs on the 'island of stability' by using the very neutron-rich radioactive beams in the near future.

The formation dynamics of SHN in the massive fusion and multinucleon transfer reactions is complicated, which is associated with the coupling of several degrees of freedom, such the radial elongation, mass or charge asymmetry, shape configuration, relative motion energy etc. Several macroscopic models are developed for describing the fusion hindrance in massive systems, e.g., the macroscopic dynamical model \cite{Bj82}, fusion-by-diffusion (FBD) model \cite{Sw05}, dynamical models based on Langevin-type equations \cite{Za07,Li11}, dinuclear system (DNS) model \cite{Ad97,Fe06,Fe11} etc. Recently, the time-dependent Hartree-Fock (TDHF) method was also applied to investigate the quasifission and fusion-fission dynamics in the reactions of $^{48}$Ca+$^{239,244}$Pu \cite{Gu18}. Modifications of macroscopic models are still needed for self-consistently and reasonably explaining the fusion dynamics in massive systems. The production cross sections of SHEs Z=119 and 120 were estimated within the multidimensional Langevin-type equations \cite{Za15} and DNS model \cite{Na09,Fe09,Ga11,Wa12,Li18,Ad20} for different reaction systems. Systematic study on the SHN production beyond oganesson (Z=118) is needed for predicting the optimal projectile-target combinations and reaction mechanism.

In this work, the stochastic diffusion in the nucleon transfer process is taken into the DNS model via the Monte Carlo procedure. The systematic analysis on the production of new superheavy elements is performed. The paper is organized as follows. In Sec. II, a brief description of the DNS model is presented. The comparison with available data and predictions of new elements Z=119-122 are discussed in Sec. III. A summary is given in IV.

\section{Model description}

The DNS model has been applied to the the quasi-fission and fusion dynamics, multinucleon transfer reactions and deep inelastic collisions, in which the dissipation of relative motion and rotation of colliding system into the internal degrees of freedom is assumed at the touching configuration. The DNS system evolves along two main degrees of freedom to form a compound nucleus, namely, the radial motion via the decay of DNS and the nucleon transfer via the mass asymmetry $\eta=(A_{1}-A_{2})/(A_{1}+A_{2})$ \cite{Ad20b,Fe09b}. In accordance with the temporal sequence, the system undergoes the capture by overcoming the Coulomb barrier, the competition of quasi-fission and complete fusion by cascade nucleon transfer, and the formation of cold residue nuclide by evaporating $\gamma$-rays, neutrons, light charged particles and binary fission.
The production cross section of the superheavy residue is estimated by the sum of partial wave with the angular momentum $J$ at incident center of mass energy $E_{c.m.}$ as,
\begin{eqnarray}
\sigma_{ER}(E_{c.m.}) = \frac{\pi\hbar^{2}}{2\mu E_{c.m.}}\sum^{J_{max}}_{J=0}(2J+1)T(E_{c.m.},J) P_{CN}(E_{c.m.},J)W_{sur}(E_{c.m.},J)
\end{eqnarray}
Here, $T(E_{c.m.},J)$ is the penetration probability and given by a Gaussian-type barrier distribution. The fusion probability $P_{CN}$ is described by the DNS model and taking into account the competition of the quasi-fission and fission of the heavy fragment. The survival probability $W_{sur}$ is calculated with the statistical approach \cite{Ch17,Ch18}.

In the DNS model, the time evolution of the distribution probability $P(Z_{1},N_{1},E_{1}, t)$ for the fragment 1 with proton number $Z_{1}$, neutron number $N_{1}$ and excitation energy $E_{1}$ is described by the following master equations as
\begin{eqnarray}
&& \frac{d P(Z_{1},N_{1},E_{1}, t)}{dt}           \nonumber\\
  &=& \sum_{Z^{'}_{1}}W_{Z_{1},N_{1};Z^{'}_{1},N_{1}}(t)[d_{Z_{1},N_{1}}P(Z^{'}_{1},N_{1},E^{'}_{1},t)    \nonumber\\
 && - d_{Z^{'}_{1},N_{1}}P(Z_{1},N_{1},E_{1},t)] + \sum_{N^{'}_{1}}W_{Z_{1},N_{1};Z_{1},N^{'}_{1}}(t)     \nonumber\\
 && \times [d_{Z_{1},N_{1}}P(Z_{1},N^{'}_{1},E^{'}_{1},t)  - d_{Z_{1},N^{'}_{1}}P(Z_{1},N_{1},E_{1}, t)]    \nonumber\\
 && -[\Lambda^{qf}_{Z_{1},N_{1},E_{1},t}(\Theta)+\Lambda^{fis}_{Z_{1},N_{1},E_{1},t}(\Theta)]P(Z_{1},N_{1},E_{1},t).   \nonumber\\
 &&
\end{eqnarray}
Here the $W_{Z_{1},N_{1};Z^{'}_{1},N_{1}}$($W_{Z_{1},N_{1};Z_{1},N^{'}_{1}}$) is the mean transition probability from the channel($Z_{1},N_{1},E_{1}$) to ($Z^{'}_{1},N_{1},E^{'}_{1}$), [or ($Z_{1},N_{1},E_{1}$) to ($Z_{1},N^{'}_{1},E^{'}_{1}$)], and $d_{Z_{1},Z_{1}}$ denotes the microscopic dimension corresponding to the macroscopic state ($Z_{1},N_{1},E_{1}$).The sum is taken over all possible proton and neutron numbers that fragment $Z^{'}_{1}$, $N^{'}_{1}$ may take, but only one nucleon transfer is considered in the model with the relations $Z^{'}_{1}=Z_{1}\pm1$ and $N^{'}_{1}=N_{1}\pm1$. The excitation energy $E_{1}$ is determined by the dissipation energy from the relative motion and the potential energy surface of the DNS. The quasi-fission rate $\Lambda^{qf}$ and fission rate $\Lambda^{fis}$ are given by the one-dimensional Kramers formula \cite{Ad03}. The motion of nucleons in the interacting potential is governed by the single-particle Hamiltonian \cite{Fe06}, which is influenced by the local excitation energy of DNS.

The potential energy surface (PES) of the DNS is given by
\begin{eqnarray}
U(\{\alpha\}) =  B(Z_{1},N_{1})+B(Z_{2},N_{2}) - \left[B(Z,N)+V^{CN}_{rot}(J)\right]   +  V(\{\alpha\}).
\end{eqnarray}
The DNS fragments satisfy the relation of $ Z_{1}+Z_{2}=Z $ and  $ N_{1}+N_{2}=N$ with the $Z$ and $N$ being the proton and neutron numbers of composite system, respectively. The symbol ${\alpha}$ denotes the quantities of $Z_{1}$, $N_{1}$, $Z_{2}$, $N_{2}$, $J$, $R$, $\beta_{1}$, $\beta_{2}$, $\theta_{1}$, $\theta_{2}$. The $B(Z_{i},N_{i}) (i=1, 2)$ and $B(Z,N)$ are the negative binding energies of the fragment $(Z_{i},N_{i})$ and the compound nucleus $(Z,N)$, respectively. The $V^{CN}_{rot}$ is the rotation energy of the compound nucleus. The $\beta_{i}$ represent the quadrupole deformations of the two fragments and are taken the ground-state values. The $\theta_{i}$ denote the polar angles between the collision orientations and the symmetry axes of the deformed nuclei. The collision direction is sampled via the Monte Carlo method with the $\theta_{i}=\frac{\pi}{2}\xi_{i}$ with $\xi_{i}$ being the random number. It should be noticed that the angles $\theta_{i} (i=1, 2)$ of binary DNS fragments are stochastically sampled and different in the nucleon transfer process. The interaction potential between fragments $(Z_{1},N_{1})$ and $(Z_{2},N_{2})$ includes the nuclear, Coulomb, and centrifugal parts as
\begin{eqnarray}
V(\{\alpha\})=V_{N}(\{\alpha\})+\frac{J(J+1)\hbar^{2}}{2\mu R^{2}}+V_{C}(\{\alpha\})
\end{eqnarray}
$\mu$ is the reduced mass of two DNS fragments. The nuclear potential is calculated by using the double-folding method based on the Skyrme interaction force without considering the momentum and the spin dependence \cite{Ad39}. The Coulomb potential is obtained by the Wong's formula \cite{Wo40}. In the calculation, the distance $R$ between the centers of the two fragments is chosen as the minimal position of interaction pocket. The minimal path in the valley of PES is named as the driving potential and dependent on the mass asymmetry. Shown in Fig. 1 is the driving potentials for the tip-tip, waist-waist collisions and average value of random orientations in the reaction of $^{48}$Ca + $^{238}$U. The tip-tip orientation means that the interaction potential is the minimal value with varying the polar angle $\theta$ corresponding to the inner fusion barrier of $B_{fus}= 9.9$ MeV, and the waist-waist case being the maximal potential for the DNS fragments with $B_{fus}= 4.8$ MeV. In all orientations, the driving potential manifests the symmetric structure. It is obvious that the driving potential with the random collisions is basically close to the average values of the tip-tip and waist-waist collisions. The inner fusion barrier is related to the collision orientation and the waist-waist collisions undergoes a high barrier to form a compound nucleus. However, the bump structure towards the decrease of mass asymmetry hinders the quasifission process.

\begin{figure*}
\includegraphics[width=16 cm]{fig1}
\caption{The driving potentials in the reaction of $^{48}$Ca+$^{238}$U with the waist-waist, tip-tip, random collisions and fixed angles, respectively. }
\end{figure*}

In order to form a compound nucleus (CN), the DNS must have enough local excitation to overcome the internal fusion barrier. The formation probability of the compound nucleus at Coulomb barrier $B$ and angular momentum $J$ is given by the summation to the B.G. (Businaro-Gallone) point as
\begin{eqnarray}
P_{CN}(E_{c.m.},J,B) =  \frac{1}{N_{t}} \sum_{i=1}^{N_{t}} \sum^{Z_{B.G.}}_{Z_{1}=1}\sum^{N_{B.G.}}_{N_{1}=1}    P(Z_{1},N_{1},E_{1}(\theta_{1},\theta_{2}),\tau_{int}(\theta_{1},\theta_{2}))   \sin(\theta_{1}) \sin(\theta_{2})
\end{eqnarray}
Here the interaction time $\tau_{int}(E_{c.m.},J,B)$ is obtained by using the deflection function method \cite{Li41}. The excitation energy $E_{1}$ of DNS fragment ($Z_{1}, N_{1}$) is related to the collision orientation $\theta_{1}$ and $\theta_{2}$. The $N_{t}$ is the total event for Monte Carlo integration with $\theta_{i}=\frac{\pi}{2}\xi$. The fusion probability is calculated with a Gaussian-type distribution $f(B)$ as
\begin{eqnarray}
P_{CN}(E_{c.m.},J) = \int f(B) P_{CN}(E_{c.m.},J,B)dB.
\end{eqnarray}
The collision orientation influences the PES of DNS because of different interaction potential with the stochastic angle. Consequently, the formation probability of CN is related to the orientation of two DNS fragments. Shown in Fig. 2 is the dependence of the fusion probability on the excitation energy at different orientation in nucleon transfer process in the reaction of $^{48}$Ca+$^{238}$U. It is obvious that the fusion probability increases with the excitation energy of compound nucleus. The waist-waist case leads to the high fusion probability owing to the lower inner fusion barrier and a small peak towards the symmetric diffusion (quasifission path). The statistical error is included for the random collisions, in which the fusion probability lies between the fixed orientations. Usually, the tip-tip orientation is chosen in the calculation of fusion probability, which is considered to be probable nucleon transfer \cite{Fe06}.

\begin{figure*}
\includegraphics[width=14 cm]{fig2}
\caption{ A comparison of the fusion probability a function of excitation energy in the reaction of $^{48}$Ca+$^{238}$U at different collision orientations. }
\end{figure*}

Once the compound nucleus is formed by the cascade nucleon transfer, the de-excitation undergoes by emitting $\gamma$-rays, particles (n, p, d, $\alpha$ etc) and binary fission. The survival probability of heavy nuclei after evaporating particle is very important for assessment of cross section, which is usually calculated with the statistical process. The probability in the channel of evaporating the $x$-th neutron, the $y$-th proton and the $z$th-$\alpha$ is expressed as \cite{Ch17}
\begin{eqnarray}
&& W_{sur}(E^{*}_{CN},x,y,z,J)=P(E^{*}_{CN},x,y,z,J)     \nonumber  \\
&&   \times\prod^{x}_{i=1}\frac{\Gamma_{n}(E^{*}_{i},J)}{\Gamma_{tot}(E^{*}_{i},J)}\prod^{y}_{j=1}\frac{\Gamma_{p}(E^{*}_{j},J)}{\Gamma_{tot}(E^{*}_{j},J)}
\prod^{z}_{k=1}\frac{\Gamma_{\alpha}(E^{*}_{k},J)}{\Gamma_{tot}(E^{*}_{k},J)}.
\end{eqnarray}
The $E^{*}_{CN}$, $J$ and $\Gamma_{tot}$ are the excitation energy, the spin of compound nucleus and the sum of partial widths of particle evaporation, respectively. The excitation energy $E^{*}_{s}$ before evaporating the $s$-th particles is given by
\begin{eqnarray}
E^{*}_{s+1}=E^{*}_{s}-B^{n}_{i}-B^{p}_{j}-B^{\alpha}_{k}-2T_{s}
\end{eqnarray}
with the initial value $E^{*}_{1}=E^{*}_{CN}$ and $s=i+j+k$. The nuclear temperature $T_{i}$ is given by $E^{*}_{i}=aT^{2}_{i}-T_{i}$ with $a$ being the level density parameter. The widths of neutron evaporation and fission are calculated by using the Weisskopf evaporation theory. The fission barrier is evaluated from the macroscopic liquid drop model and the shell correction energy, and which is given as
\begin{equation}
B_{f}(E^{*},J)=B^{LD}_{f}+E_{shell}\exp(-E^{*}/E_{D}),
\end{equation}
where the macroscopic part $B^{LD}_{f}$ is calculated by liquid-drop model. The microcosmic shell correction energy is calculated with the Strutinsky method and obtained from Ref. \cite{Mo95}. The damping energy $E_{D}$ is associated with the level density and mass number of compound nucleus, in which the shell correction and excitation energy dependence are taken into account in the calculation \cite{Ch17}.

\section{Results and discussion}

The production rate of SHNs in massive fusion reactions is very low owing to the fusion hindrance in heavy systems, which enables the quasifission process in the binary collisions. The evaporation residue (ER) excitation functions in different channels are favorable for experimental measurements with optimal projectile-target combination and suitable beam energy. The reaction dynamics in competition of the quasifission and fusion-fission reactions, level density, separation energy of evaporated particle and fission barrier of compound nucleus influence the ER cross sections. As a test of the model, shown in Fig. 3 is the ER excitation functions in the reaction of $^{48}$Ca+$^{238}$U and compared with the data from Dubna \cite{Og43}. The ER cross sections for producing SHN strongly depend on the orientations of both DNS fragments in the nucleon transfer process. The tip-tip collisions have the lower cross sections in comparison with the waist-waist orientations and roughly lead to the two order reduction because of the higher inner fusion barrier for merging the compound nucleus. The cross sections with the stochastic selection of collision angle of two DNS fragments lie between the tip-tip and waist-waist orientation. The results with the tip-tip collisions are consistent with the available data and are chosen in the calculation. It is obvious that the maximal yield of SHN $^{283}$Cn production with 3 pb is positioned in the 3n channel via the tip-tip collisions at the excitation energy of 35 MeV. The SHN is still far away from the neutron shell closure N=184. New reaction mechanism is expected for creating the neutron-rich SHN. The pure neutron channels are dominant decay modes for surviving SHN. The cross sections with mixed channels of proton and $\alpha$ are reduced to be two order magnitude with the 3 particle channels, i.e., p2n, $\alpha$2n. The charged particle channels are of significance on the production of proton-rich actinide nuclides close to the drip line in the fusion-evaporation reactions and the multinucleon transfer dynamics \cite{Ch20}.

\begin{figure*}
\includegraphics[width=14 cm]{fig3}
\caption{ Calculated ER cross sections with different channels of (2-5)n, 1p(2-4)n and $1\alpha$(1-3)n and compared with the experimental data in the reaction of $^{48}$Ca+$^{238}$U \cite{Og43}.}
\end{figure*}

\begin{figure*}
\includegraphics[width=14 cm]{fig4}
\caption{ ER excitation functions with pure neutron channels and compared with the experimental data in the reactions of $^{48}$Ca+$^{244}$Pu, $^{248}$Cm, $^{249}$Bk, $^{249}$Cf \cite{Og43,Og44}. }
\end{figure*}

The superheavy elements from Fl (Z=114) to Og (Z=118) have been successfully synthesized with $^{48}$Ca induced reactions on actinide targets. It manifests the strong shell effect in the production and the decay chains. Shown in Fig. 4 is the ER excitation functions for producing the superheavy elements 114-118 with $^{244}$Pu, $^{248}$Cm, $^{249}$Bk and $^{249}$Cf as the targets. Different evaporation channels are distinguished by colored lines, and compared with the experimental data \cite{Og43,Og44}. It is obvious that the 3n and 4n channels are dominant decay modes for SHN production. The 2n channel is pronounced with increasing the charge number of SHN. However, the cross sections in the 4n and 5n channels decrease with becoming heavier SHN. Different with the cold fusion reactions, the maximal cross section in the $^{48}$Ca induced reactions weakly depends on the mass of ER nucleus, e.g., the value of 8 pb for $^{288}$Fl and 0.6 pb for $^{294}$Og. The difficulty is the construction of target material for experiments to synthesize new SHN. The hot fusion reactions also provide a possible way for creating the new elements at the eighth period with projectiles $^{45}$Sc, $^{50}$Ti, $^{54}$Cr, $^{58}$Fe, $^{64}$Ni etc. Besides the reaction dynamics in the formation of SHN in the hot fusion reactions, the nuclear structure effects are of importance in the evaluation of production cross section, i.e., shell effect, neutron separation energy, odd-even effect, microscopic state of compound nucleus (isomeric state) etc.

\begin{figure*}
\includegraphics[width=14 cm]{fig5}
\caption{ The evaporation residue cross sections with channels of (2-5)n of element Z=119 in collisions by different reactions.}
\end{figure*}

\begin{table}[!htb]
	\centering
	\caption{Optimal evaporation residual cross sections via different reactions leading to the formation of SHE Z=119. }
	\label{tab1}
		\begin{tabular}{ccccccc}
	Reaction systems  &$\sigma_{ER}(pb)$    &$E_{CN}^{\ast}$(MeV)       &References  \\
        \hline
            $^{249}$Bk($^{50}$Ti,3n)$^{296}119$    &0.04    &41           & \cite{Za08}   \\
            $^{249}$Bk($^{50}$Ti,4n)$^{295}119$    &0.06    &44          & \cite{Za08}  \\
            $^{254}$Es($^{48}$Ca,3n)$^{299}119$  &0.3       &35          &\cite{Za12}   \\
	    	$^{252}$Es($^{48}$Ca,4n)$^{296}$119  &0.2       &43          &\cite{Si46}   \\
			$^{254}$Es($^{48}$Ca,4n)$^{298}119$  &0.015   &41           &\cite{Si46}    \\
            $^{249}$Bk($^{50}$Ti,4n)$^{295}119$   &0.03     &36          &\cite{Si46}    \\
            $^{249}$Bk($^{50}$Ti,3n)$^{296}119$  &0.035   &27          &\cite{Wa12}  \\
            $^{249}$Bk($^{50}$Ti,4n)$^{295}119$  &0.11    &39            &\cite{Wa12}  \\
            $^{249}$Bk($^{50}$Ti,4n)$^{295}119$  &0.57   &41            &\cite{Liu11}  \\
            $^{252}$Es($^{44}$Ca,3n)$^{293}119$ &4.32  &35            &\cite{Li18}  \\
            \cline{1-4}
            $^{251}$Cf($^{45}$Sc,3n) $^{293}119$ &0.38  &37     & This work\\
            $^{249}$Cf($^{45}$Sc,3n)$^{291}119$ &0.99   &37      \\
            $^{247}$Bk($^{50}$Ti,4n)$^{293}119$ &0.024  &45    \\
            $^{249}$Bk($^{50}$Ti,4n)$^{295}119$ &0.013   &45

            \end{tabular}
\end{table}

Attempts to synthesize the superheavy elements 119 and 120 were performed at different laboratories in the world, i.e, FLNR, GSI, GANIL etc. But no decay chains were observed up to now in experiments. Theoretical predictions with various models were also made for producing the element 119 with a number of systems, e.g., $^{249}$Cf ($^{45}$Sc,xn) $^{294-x}$119, $^{251}$Cf ($^{45}$Sc,xn) $^{296-x}$119, $^{247}$Bk ($^{50}$Ti,xn) $^{297-x}$119, $^{249}$Bk ($^{50}$Ti,xn) $^{299-x}$119, $^{254}$Es($^{48}$Ca,xn)$^{302-x}119$ etc. Some results for the optimal ER cross sections in different reactions leading to the formation of Z=119 are listed in Table I. One can observe that the optimal system of $^{252}$Es ($^{44}$Ca,3n) $^{293}119$ is possible with a larger cross section of 4 pb by the DNS model \cite{Li18}. Further confirmation with different models is expected for checking the reliability of the calculation. On the other hand, the construction of target material $^{252}$Es in experiment is very difficulty. The reaction of $^{45}$Sc+$^{249,251}$Cf is also feasible for synthesizing the new element with the cross section above 0.1 pb. Difference of one order magnitude for producing the element 119 in the reaction $^{50}$Ti+$^{249}$Bk exists in the model predictions, e.g., 0.03 pb in the 4n channel at the excitation energy of 36 MeV by the FBD model \cite{Si46}, 0.57 pb at 41 MeV by the diffusion model with Langevin-type equations \cite{Liu11}, 0.11 pb in Ref. \cite{Wa12} and 0.013 pb in our calculation by the DNS model. The 3n and 4n channels are optimal way to create the new element as shown in Fig. 5. The isotopic dependence of ER cross sections is weak. More large mass asymmetry in the bombarding system is favorable to produce SHN. The synthesis of element 119 in laboratories is possible in the forthcoming experiments with the high-intensity accelerators in the world. The reliable prediction in theories is helpful for experimental managements.

\begin{figure*}
\includegraphics[width=14 cm]{fig6}
\caption{ The same as in Fig. 5, but for the production of element Z=120. }
\end{figure*}

\begin{table}[!htb]
	\centering
	\caption{The same as in table I, but for the production of SHE Z=120.}
	\label{tab2}
		\begin{tabular}{ccccccc}
			Reaction systems     &$\sigma_{ER}(pb)$     &$E_{CN}^{\ast}$(MeV)     &References \\ \hline
	       $^{249}$Cf($^{50}$Ti,4n)$^{295}$120 &0.006      &43           &\cite{Si46}   \\
		   $^{251}$Cf($^{50}$Ti,4n)$^{297}120$  &0.003      &42      \\
           $^{248}$Cm($^{54}$Cr,4n)$^{298}120$ &0.001   &35       \\
            \cline{1-4}
            $^{244}$Pu($^{58}$Fe,3n)$^{299}120$  &0.01  &36  &\cite{Si50} \\
            $^{238}$U($^{64}$Ni,3n)$^{299}120$  &0.007  &36     \\
            $^{248}$Cm($^{54}$Cr,3n)$^{299}120$ &0.076  &36   \\
            $^{249}$Cf($^{50}$Ti,3n)$^{296}120$ &0.76  &33         \\
            \cline{1-4}
            $^{249}$Cf($^{50}$Ti,3n)$^{296}120$ &0.1  &29   &\cite{Na51} \\
            $^{248}$Cm($^{54}$Cr,3n)$^{299}120$ &0.055  &30   \\
            \cline{1-4}
            $^{249}$Cf($^{50}$Ti,4n)$^{295}120$ &0.046  &43    &\cite{Za08} \\
            $^{248}$Cm($^{54}$Cr,4n)$^{298}120$ &0.028  &43   \\
            \cline{1-4}
            $^{249}$Cf($^{50}$Ti,3n)$^{296}120$ &0.06  &36   &\cite{Li11} \\
            $^{250}$Cf($^{50}$Ti,3n)$^{297}120$ &0.12  &37   \\
            $^{251}$Cf($^{50}$Ti,4n)$^{297}120$ &0.11  &38    \\
            $^{252}$Cf($^{50}$Ti,4n)$^{298}120$ &0.25  &38  \\
            \cline{1-4}
            $^{251}$Cf($^{50}$Ti,3n)$^{298}120$ &0.25  &36   &\cite{Wa12} \\
            $^{249}$Cf($^{50}$Ti,3n)$^{296}120$ &0.05  &33           \\
            $^{248}$Cm($^{54}$Cr,4n)$^{298}120$  &0.005  &42    \\
            $^{244}$Pu($^{58}$Fe,4n)$^{298}120$  &0.003  &43     \\
            $^{249}$Cf($^{50}$Ti,3n)$^{296}120$  &0.02  &31           \\
            \cline{1-4}
            $^{257}$Fm($^{40}$Ca,3n)$^{294}120$ &1.24  &48   &\cite{Li18}  \\
            $^{248}$Cf($^{46}$Ti,2n)$^{292}120$  &0.17  &34   \\
            $^{249}$Cf($^{46}$Ti,3n)$^{292}120$  &0.24  &39   \\
            $^{250}$Cf($^{46}$Ti,2n)$^{294}120$  &0.13  &36   \\
            $^{251}$Cf($^{46}$Ti,3n)$^{294}120$  &0.37  &39   \\
            \cline{1-4}
            $^{251}$Cf($^{50}$Ti,3n)$^{298}120$ &0.11  &33   &This work \\
            $^{251}$Cf($^{48}$Ti,2n)$^{297}120$ &0.25  &25         \\
            $^{244}$Pu($^{58}$Fe,3n)$^{299}120$  &0.004  &33   \\
            $^{244}$Pu($^{62}$Fe,3n)$^{303}120$  &0.0004  &31  \\
            $^{248}$Cm($^{54}$Cr,3n)$^{299}120$  &0.004  &33  \\
            $^{248}$Cm($^{52}$Cr,2n)$^{300}120$  &0.37  &25    \\
            $^{238}$U($^{64}$Ni,3n)$^{299}120$  &0.001  &31   \\
            $^{238}$U($^{62}$Ni,2n)$^{300}120$  &0.001  &27   \\
            $^{252}$Es($^{44}$Sc,3n)$^{293}120$  &3.18  &35   \\
            $^{252}$Es($^{45}$Sc,3n)$^{293}120$  &0.59  &35   \\
            $^{249}$Bk($^{49}$V,2n)$^{296}120$  &0.18  &27   \\
            $^{249}$Bk($^{51}$V,2n)$^{298}120$  &0.1  &27
 		\end{tabular}
\end{table}

The synthesis of superheavy element 120 is particularly important in understanding the shell structure in the domain of SHNs. The strongly shell effect enhances the fission barrier and $\alpha$ decay half-life, which is also favorable for producing and surviving SHN. Possible experiment is being planned at HIAF with high-intensity beams. Shown in Fig. 6 is the systematic comparison for producing SHN of Z=120 with actinide nuclides based reactions. Different panels represent different reaction systems to create the superheavy element 120, but with different neutron richness. The compound nuclei formed by different combinations are close to the neutron shell closure N=184. The channels of 2-4n are available for the SHN production with the excitation energies around 25 MeV, 30 MeV and 40 MeV, respectively. The maximal ER cross section depends on the isotopic projectiles. It is obvious that the 3n channel of the  $^{44,45}$Sc + $^{252}$Es reaction is favorable for synthesizing the new SHN owing to the large mass asymmetry in the entrance system. The smaller neutron separation energy is available for cooling the compound nucleus formed in the fusion reactions. In Table II, the production cross sections of Z=120 with the optimal channels and possible combinations are compared. Difference of model predictions exists in the calculations, e.g., the production of $^{295}$120 in collisions of $^{50}$Ti on $^{249}$Cf being 0.006 pb and 0.046 pb by the FBD model \cite{Si46} and multidimensional Langevin-type equations \cite{Za08}, respectively. Calculations support the 3n channel in collisions of $^{50}$Ti on californium isotopes is available for synthesizing the element 120 with the cross section above 0.1 pb at the excitation energy around 35 MeV. The 2n channel in the reactions of $^{49,51}$V+$^{249}$Bk is also possible for creating the new element. The position of maximal cross section is mainly determined by the odd-even effect of neutron evaporation and the energy dependence of fusion probability. It should be noticed that the proton shell closure Z=120 was predicted with the relativistic mean-field model by including the isospin dependence of the spin-orbit interaction and the effective mass \cite{Be99}. The shell correction energy calculated by the macroscopic-microscopic model is used in the calculation and the Z=114 proton shell closure is given by the approach \cite{Mo95}. The production cross section of element 120 is to be enhanced with the shell correction by the relativistic mean-field model.

\begin{figure*}
\includegraphics[width=14 cm]{fig7}
\caption{ The evaporation residue cross sections with channels of (2-5)n of element Z=121 in collisions by different combinations.}
\end{figure*}

\begin{figure*}
\includegraphics[width=14 cm]{fig8}
\caption{ The same as in Fig. 5, but for the production of element Z=122.}
\end{figure*}

As the extension of the model prediction, we analyzed the formation of superheavy elements 121 and 122 in massive fusion reactions. Shown in Fig. 7 is the ER excitation functions in the reactions of $^{54}$Cr+$^{247,249}$Bk and $^{58,64}$Fe+$^{243}$Am. In panels (a) and (b), the stable nuclide $^{54}$Cr and isotopic target nuclei are selected. It can be seen that the 2-4n channels in the $^{54}$Cr+$^{249}$Bk reaction are favorable, but with very low cross section at the level of 1 fb. The influence of stable and radioactive nuclides on the SHN production is compared in panels (c) and (d). The different structure of 2-5 channels is caused from the neutron separation energy in the cascade evaporation, i.e., the smaller separation energy for the compound nucleus $^{301}$121 resulting in the larger 2n channel probability. The production cross section of element 122 is very low in the fusion reactions as shown in Fig. 8. Both the neutron separation energy and fission barrier are sensitive to the survival of SHN. The polar angle distribution, mass, total kinetic energy spectra and excitation energy dependence of fission fragments from the SHNs are useful for extracting the fission barrier and shell evolution. More experiments are expected in the future. The systems of $^{54}$Cr+ $^{249,251}$Cf and $^{58,64}$Fe+$^{248}$Cm are chosen for synthesizing the element 122 as shown in Fig. 8. The cross section below 0.1 fb is already out the limit of experimental measurement. The new reaction mechanism is expected for creating the neutron-rich SHN and new element, i.e., the multinucleon transfer reaction, incomplete fusion with radioactive nuclide etc.

\section{Conclusions}

Within the framework of the DNS model, the SHN formation in the fusion-evaporation reactions is thoroughly investigated. The stochastic collision orientations in the nucleon transfer process are implemented into the model via the Monte Carlo approach. The calculated results are consistent with the experimental data from Dubna. The maximal cross sections of evaporation residues appear in the (2-5)$n$ evaporation channels. The yields in the 1$pxn$, 1$\alpha xn$ and 1$p$1$\alpha xn$ evaporation channels are much lower than the pure neutron evaporation. The reactions of $^{249}$Cf ($^{45}$Sc,xn)$^{294-x}$119, $^{251}$Cf($^{45}$Sc,xn)$^{296-x}$119, $^{247}$Bk($^{50}$Ti,xn) $^{297-x}$119 and $^{249}$Bk($^{50}$Ti,xn)$^{299-x}$119 are investigated for synthesizing the new element 119. It is concluded that the maximum cross sections are close to 1 pb in the 3$n$ evaporation channel for the systems and weakly depend on the isotopic target nucleus. The synthesis of the element $Z=120$ is investigated by a series of isotopic projectile nuclei bombarding on actinide targets. The optimal combination is the reaction of $^{44}$Sc+$^{252}$Es in the 3$n$ channel with the cross section of 3 pb. The production of superheavy elements 121 and 122 is obtained at the level of below 1 fb in the massive fusion reactions. New reaction mechanism is still needed to be explored for the new element production. Therefore, the synthesis of new SHN in experiments provides a good theoretical basis for the selection of collision combination.

\section{Acknowledgements}

This work was supported by the National Natural Science Foundation of China (Projects No. 12175072 and No. 11722546) and the Talent Program of South China University of Technology.


\section{\label{sec:level1}Introduction}

Remote synchrony (RS) is an interesting manifestation of the general and highly significant nonlinear phenomenon of synchronization~\cite{Kuramoto-84,*pikovsky2003synchronization,*strogatz2004sync,*Osipov-Kurths-Zhou-07}. RS implies adjusting rhythms of oscillators that do not interact directly but only through an asynchronous unit (hub).
Exploration of this effect, initially described by Bergner et al.~\cite{bergner2012RS} and further studied numerically and experimentally in Refs.~\cite{minati2015remote,*karakaya2019fading},
is crucial, e.g., for understanding functional connectivity in brain networks~\cite{vuksanovic2014functional,vlasov2017hub}.

Previous studies analyzed RS in star-like and complex networks of Stuart-Landau (SL) or phase oscillators~\cite{bergner2012RS,minati2015remote,*karakaya2019fading,gambuzza2013analysis,vlasov2017hub}. The results uncovered the role of amplitude dynamics~\cite{bergner2012RS,gambuzza2013analysis}:
RS appeared in a  network of isochronous SL units but not in its first-order phase approximation, i.e., in the Kuramoto network.
Furthermore, Vlasov and Bifone~\cite{vlasov2017hub} demonstrated that RS emerges in
networks of phase oscillators with the Kuramoto-Sakaguchi interaction~\cite{sakaguchi1986soluble}, but not in the case of zero phase shift in the sine-coupling term. Since the Kuramoto-Sakaguchi model is the first-order approximation of coupled non-isochronous SL oscillators, this result indicates the role of non-isochronicity in promoting RS.
However, the understanding of mechanisms leading to RS is yet incomplete. This paper uses a simple motif of three coupled SL oscillators to analyze the transition to RS. In contradistinction to \cite{vlasov2017hub}, we consider non-identical peripheral oscillators.
Using recent results on high-order phase reduction~\cite{gengelphasereduction}, we explain the contribution of both the non-isochronicity and amplitude dynamics and quantitatively describe the transition to RS. We demonstrate the importance of high-order phase approximation in the explanation of RS.

The paper is organized as follows. In Section~\ref{sec:model}, we introduce the
model and its second-order phase approximation. Next, we demonstrate the transition to RS in this model. In Section~\ref{sec:RS_theoretical}, we derive the condition for this transition and
in Section~\ref{sec:ref}, we present our results.
Section~\ref{sec:concl} concludes and discusses our findings.

\section{Remote synchrony in coupled Stuart-Landau oscillators}
\label{sec:model}
Consider three SL oscillators coupled in a chain as $1\xleftrightarrow{}2\xleftrightarrow{}3$. Thus, peripheral units 1 and 3 are not interacting directly but only through the central oscillator.
Let the (generally different) natural frequencies of the oscillators be $\omega_{1,2,3}$. Correspondingly, we denote the  frequencies of interacting units (observed frequencies) as $\Omega_{1,2,3}$. Following Bergner et al.~\cite{bergner2012RS}, we say that the network reaches a state of RS if, with an increase of coupling stength, $\Omega_1$ becomes equal to $\Omega_3$ while $\Omega_1\ne\Omega_2$. If all frequencies coincide, $\Omega_1=\Omega_2=\Omega_3$, then we speak about complete synchrony (CS).
We emphasize that Refs.~\cite{qin2020mediated,*qin2018stability,*nicosia2013remote} use the term RS in a different context.

In the rest of this Section, we first specify our model and present its second-order phase approximation. Next, we numerically demonstrate transitions from asynchrony to RS and CS in the full model and its phase-reduced versions.

\subsection{Model and its phase approximation}

The governing equations of the model are:
\begin{equation}\label{eqn:SLequations}
    \dot{A}_{n}=\left[1+\mathrm{i} (\omega_{n}+\alpha)\right]  A_{n}-(1+\mathrm{i}\alpha)\left|A_{n}\right|^{2} A_{n}+\varepsilon I_n \;,
\end{equation}
where $A_n \in \mathbb{C}$, $n=1,2,3$, $\omega_n$ is the natural frequency of the $n$-th oscillator,
and $\alpha$ is the non-isochronicity parameter, common for all units.
The parameter $\varepsilon$ and the terms $I_1=A_2$, $I_2=A_1+A_3$, $I_3=A_2$ describe the strength and structure of the coupling, respectively.

It is well-known that for sufficiently weak coupling, the dynamics of interacting limit-cycle oscillators reduce to that of phases.
For the coupled SL oscillators, the first-order phase approximation in $\varepsilon$ can be performed analytically because the phase of this system can be readily obtained from the state variable $A$; the reduction yields the celebrated Kuramoto-Sakaguchi phase equations \cite{sakaguchi1986soluble}.
However, phase reduction beyond the first-order approximation remains challenging and is a subject of ongoing research.  Here, we use the results of Gengel et al.~\cite{gengelphasereduction}, who
provided expressions for the second-order reduction of coupled SL oscillators~\footnote{Notice that our Eq.~(\ref{eqn:SLequations}) represents a particular case of a more general setup studied in~\cite{gengelphasereduction}}.
Let the phase of the $n$-th oscillator be $\varphi_n$. The second-order phase approximation of the system (\ref{eqn:SLequations}) reads:
\begin{widetext}
\begin{equation}
\begin{aligned}\label{eqn:second_order_phase_reduction}
\dot{\varphi}_{1}=& \omega_{1}+\varepsilon \left[\sin \left(\varphi_{2}-\varphi_{1}\right)-\alpha \cos \left(\varphi_{2}-\varphi_{1}\right)\right]\\
&+\varepsilon^{2}\left[D_{32} \cos \left(2 \varphi_{2}-\varphi_{1}-\varphi_{3}\right)+C_{32} \sin \left(2 \varphi_{2}-\varphi_{1}-\varphi_{3}\right)-D_{32} \cos \left(\varphi_{3}-\varphi_{1}\right)+C_{32} \sin \left(\varphi_{3}-\varphi_{1}\right)\right] +\mathcal{O}(\varepsilon^3)\;, \\
\dot{\varphi}_{2}=& \omega_{2}+\varepsilon \left[\sin \left(\varphi_{1}-\varphi_{2}\right)-\alpha \cos \left(\varphi_{1}-\varphi_{2}\right)+\sin \left(\varphi_{3}-\varphi_{2}\right)-\alpha \cos \left(\varphi_{3}-\varphi_{2}\right)\right] \\
&+\varepsilon^{2}\left[(D_{12}+D_{32}) \cos \left(2 \varphi_{2}-\varphi_{1}-\varphi_{3}\right)+(C_{12}+C_{32}) \sin \left(2 \varphi_{2}-\varphi_{1}-\varphi_{3}\right)\right.\\
&-\left.(D_{12}+D_{32}) \cos \left(\varphi_{1}-\varphi_{3}\right)+(C_{12}-C_{32})\sin(\varphi_1 - \varphi_3)\right] +\mathcal{O}(\varepsilon^3)\;,\\
\dot{\varphi}_{3}=& \omega_{3}+\varepsilon \left[\sin \left(\varphi_{2}-\varphi_{3}\right)-\alpha \cos \left(\varphi_{2}-\varphi_{3}\right)\right] \\
&+\varepsilon^{2}\left[D_{12} \cos \left(2 \varphi_{2}-\varphi_{3}-\varphi_{1}\right)+C_{12} \sin \left(2 \varphi_{2}-\varphi_{3}-\varphi_{1}\right) -D_{12} \cos \left(\varphi_{1}-\varphi_{3}\right)+C_{12} \sin \left(\varphi_{1}-\varphi_{3}\right)\right]+\mathcal{O}(\varepsilon^3)\;,
\end{aligned}
\end{equation}
\end{widetext}
where
\begin{equation}\label{eqn:second_order_phase_reduction_constant_C}
    C_{ij}=\frac{1+\alpha^2}{4+(\omega_i-\omega_j)^2}
\end{equation}
and
\begin{equation}\label{eqn:second_order_phase_reduction_constant_D}
    D_{ij}=\frac{1+\alpha^2}{2}\left(\frac{\omega_i-\omega_j}{4+(\omega_i-\omega_j)^2} \right)\;.
\end{equation}

Keeping in Eq.~(\ref{eqn:second_order_phase_reduction}) only the first-order terms $\sim\varepsilon$, one obtains the Kuramoto-Sakaguchi model:
\begin{equation}
\begin{aligned}\label{eqn:first_order_phase_reduction}
\dot{\varphi}_{1}=& \omega_{1}+\varepsilon \left[\sin \left(\varphi_{2}-\varphi_{1}\right)-\alpha \cos \left(\varphi_{2}-\varphi_{1}\right)\right]\;, \\
\dot{\varphi}_{2}=& \omega_{2}+\varepsilon \left[\sin \left(\varphi_{1}-\varphi_{2}\right)-\alpha \cos \left(\varphi_{1}-\varphi_{2}\right)\right.\\
&+\left.\sin \left(\varphi_{3}-\varphi_{2}\right)-\alpha \cos \left(\varphi_{3}-\varphi_{2}\right)\right] \;, \\
\dot{\varphi}_{3}=& \omega_{3}+\varepsilon \left[\sin \left(\varphi_{2}-\varphi_{3}\right)-\alpha \cos \left(\varphi_{2}-\varphi_{3}\right)\right] \;.
\end{aligned}
\end{equation}
For isochronous oscillators, $\alpha=0$, the model simplifies to the Kuramoto network.

\subsection{Remote synchrony in the full and reduced models}
\label{sec:RS_numerical}

\begin{figure}
    \includegraphics[width=\columnwidth]{Kumar_Fig1.png}
    \caption{\label{fig:compare_approximations_figures}
    (Color online) Numerically computed bifurcation diagrams illustrating the dependence of the system's observed state on the coupling strength, $\varepsilon$, and non-isochronicity parameter, $\alpha$. The oscillators' natural frequencies are (a) $\omega_1=1$, $\omega_2=\sqrt{2}$, $\omega_3=1.002$ and (b) $\omega_1=1$, $\omega_2=\sqrt{7}$, $\omega_3=1.01$.
    The white, blue (dark gray), and red (light gray) regions correspond to asynchrony, RS, and CS, respectively, upon numerical simulation of  Eq.~\eqref{eqn:SLequations}. The solid black line depicts the RS transition border as computed using the second-order phase approximation, see Eqs.~\eqref{eqn:second_order_phase_reduction}. The dashed black line shows the RS transition obtained for the first-order phase approximation, see
    Eqs.~\eqref{eqn:first_order_phase_reduction}. The diagrams demonstrate the crucial role of the non-isochronicity parameter $\alpha$. Furthermore, the diagrams clearly show the advantage of the second-order approximation.
    }
\end{figure}

This section compares and contrasts the regions of RS obtained using the SL system~\eqref{eqn:SLequations} and the phase approximations, see Eqs.~(\ref{eqn:first_order_phase_reduction},\ref{eqn:second_order_phase_reduction}). To this end, we fix the natural frequencies of all three oscillators~\footnote{In the following, we always consider frequencies of the peripheral oscillators to be close while the frequency of the hub is essentially different.}, numerically simulate the governing equations, and detect regions of asynchrony, CS and RS upon varying the coupling strength and the non-isochronicity parameter. (The description of the numerical procedures are deferred to Section \ref{sec:ref}.) This results in two-parameter bifurcation diagrams on the $\varepsilon$-$\alpha$ plane shown in Fig.~\ref{fig:compare_approximations_figures}.

Figure \ref{fig:compare_approximations_figures} provides us with two insights. Firstly, we note that the first-order approximation does not accurately reproduce the transition to RS. This approximation's failure results from not accounting for the amplitude modulation in the coupled SL oscillators. On the other hand, the second-order approximation fares well and is accurate for small and moderate coupling strengths.
Secondly, the non-isochronicity parameter essentially affects the transition to RS.  Generally, RS in the SL system \eqref{eqn:SLequations} appears for both the isochronous ($\alpha=0$) and the non-isochronous ($\alpha\ne 0$) cases. However, this feature is captured only by the second-order approximation; the first approximation does not exhibit RS for $\alpha=0$, in agreement with the results by Vlasov and Bifone \cite{vlasov2017hub}.

\section{Theoretical analysis of the phase dynamics}
\label{sec:RS_theoretical}
We use the phase equations~\eqref{eqn:second_order_phase_reduction}
to investigate the transition to RS.  It is straightforward to reduce Eqs.~\eqref{eqn:second_order_phase_reduction}
to a two-dimensional system for the phase differences:
\begin{equation}\label{eqn:define_phase_differences}
    \gamma_{13}=\varphi_1 - \varphi_3, \quad \gamma_{12}=\varphi_1 - \varphi_2 \;.
\end{equation}
The resulting equations represent the dynamics on a two-torus and can be studied using standard phase plane analysis techniques. In terms of the phase differences, the asynchronous state corresponds to an unbounded growth (or decline) of $\gamma_{13}$ and $\gamma_{12}$. Upon increasing the coupling strength, one observes RS, wherein $\gamma_{13}$ is bounded while $\gamma_{12}$ is unbounded. For transparency and brevity, we present our theory by analyzing the first-order phase equations. Then we provide the results of the same approach applied to the second-order model.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{Kumar_Fig2.png}
\caption
(Color online) Illustration of the RS transition using the phase approximation equations upon increasing the coupling strength.
(a) and (c): Phase plane representations of the phase difference dynamics at low and moderate coupling strengths, respectively. The bold blue (dark gray) and red (light gray) trajectories in panel (c) depict the stable and unstable limit cycles.
(b) and (d): Poincar\'e maps corresponding to the flows shown in panels (a) and (c), respectively, constructed using the Poincar\'e section $\gamma_{12}=2\pi$.
The birth of a stable limit cycle (fixed point) in the phase plane plot  (Poincar\'e map) corresponds to the onset of RS.
(The plots shown here are computed using the first-order approximation [Eq.~\eqref{eqn:first_order_phase_reduction}] but they are qualitatively identical for the second-order phase reduction case as well.)}
\label{fig:model}
\end{figure}

\subsection{Poincar\'e map}

The transition to RS corresponds to the appearance of a stable limit cycle (LC) on the torus. Figure 2a depicts a typical situation for the asynchronous regime at low coupling strengths. There are no attractors on the phase plane, the motion is quasiperiodic, and the phase differences $\gamma_{13}$ and $\gamma_{12}$ are unbounded. Figure 2c exemplifies the RS state once the coupling strength increases. A stable and an unstable limit cycle are born via a saddle-node bifurcation of LCs. Notice that on the LC, $\gamma_{12}$ is unbounded while $\gamma_{13}$ is bounded, which indicates the emergence of RS. Notice also that we consider $\omega_1<\omega_2$ for definiteness for the remainder of this article. Hence, the flow is from right to left. We have verified that our conclusions hold equally well for the other case.

For the following derivation, it is instructive to construct a Poincar\'e map, choosing the line $\gamma_{12}=2\pi$  as the Poincar\'e section. A trajectory that begins on this section intersects it next at $\gamma_{12}=0$, since the flow on the torus is leftwards. Thus, we have  $\gamma_{13}(0)=P(\gamma_{13}(2\pi))$, where $P(\cdot)$ denotes the Poincar\'e map. The Poincar\'e map corresponding to Figs.~2a and 2c are shown in Figs.~2b and 2d, respectively. Evidently, RS in the system equates to a stable fixed point of the Poincar\'e map. We exploit this observation to derive the condition for RS analytically.

\subsection{First-order phase dynamics}
Starting with Eqs.~\eqref{eqn:first_order_phase_reduction}, using Eq.~\eqref{eqn:define_phase_differences},
and introducing the new time $\tau=(\omega_2-\omega_1) t$,
we obtain a two-dimensional system for phase differences:
\begin{equation}
    \begin{aligned}\label{eqn:phase_differences_first_order_phase_reduction}
    \gamma'_{13}&=\nu+\widetilde{\varepsilon}\left[ -\sin\gamma_{12}-\alpha\cos\gamma_{12} - \sin(\gamma_{13}- \gamma_{12})\right. \\
    &\quad +\left.\alpha \cos (\gamma_{13} -\gamma_{12}) \right]\;,\\
    \gamma'_{12}&=-1+\widetilde{\varepsilon}\left[ -2\sin\gamma_{12} - \sin(\gamma_{12} -\gamma_{13})\right.\\
    & \quad + \left.\alpha\cos(\gamma_{12} -\gamma_{13}) \right]\;,
    \end{aligned}
\end{equation}
where
\begin{equation}
\begin{gathered}\label{eqn:my_params}
    \nu=\frac{\omega_1 - \omega_3}{\omega_2-\omega_1}\;, \quad \widetilde{\varepsilon}=\frac{\varepsilon}{\omega_2-\omega_1}\;,
\end{gathered}
\end{equation}
and $(\cdot)'$ denotes differentiation with respect to $\tau$.

To derive the Poincar\'e map $\gamma_{13}(0)=P(\gamma_{13}(2\pi))$, we divide the preceding equations to obtain:
\begin{widetext}
\begin{equation}\label{eqn:first_order_gamma13/gamma12}
\begin{aligned}
    \frac{\mathrm{d} \gamma_{13}}{\mathrm{d}\gamma_{12}} &= \frac{\nu+\widetilde{\varepsilon}\left[ -\sin\gamma_{12}-\alpha\cos\gamma_{12} - \sin(\gamma_{13}- \gamma_{12}) +\alpha \cos (\gamma_{13} -\gamma_{12}) \right]}{-1+\widetilde{\varepsilon}\left[ -2\sin\gamma_{12} - \sin(\gamma_{12} -\gamma_{13}) + \alpha\cos(\gamma_{12} -\gamma_{13}) \right]}\;.
\end{aligned}
\end{equation}
\end{widetext}
We solve Eq.~\eqref{eqn:first_order_gamma13/gamma12} with the initial condition $\gamma_{13}(2\pi)$  using a perturbation approach, for which we assume the following:
\begin{equation}\label{eqn:assumptions_sys_params}
    |\omega_1 - \omega_2| \sim \mathcal{O}(1),\quad 0<|\omega_1 - \omega_3| \ll 1, \quad \varepsilon \ll 1\;.
\end{equation}
Notice that the first pair of assumptions formally encapsulates our previous qualitative description: the peripheral oscillators are near-identical, whereas the hub oscillator is markedly different. Equivalently, in terms of the parameters present in Eq.~\eqref{eqn:phase_differences_first_order_phase_reduction}, the assumptions result in
$\widetilde{\varepsilon}\ll 1$ and $\nu \ll 1$.

The solution presented in Appendix~\ref{app:1st_order} provides the condition for the existence of the Poincar\'e map's fixed point:
\begin{equation}\label{eqn:first_order_RS_condition}
    \left| \frac{\varepsilon^2 \alpha}{(\omega_1-\omega_3)(\omega_1-\omega_2)} \right| \geq \frac{1}{2}\;.
\end{equation}
This inequality yields the necessary condition for RS in the first-order phase reduction Eqs.~\eqref{eqn:first_order_phase_reduction}. Its validity depends on the smallness of $\varepsilon$. It indicates that upon increasing the coupling strength, RS appears due to non-isochronicity. Hence, RS is impossible in a chain of three non-identical Kuramoto equations. This result agrees with the observation reported in
Ref.~\cite{bergner2012RS} and theoretical analysis in Ref.~\cite{vlasov2017hub}.

\subsection{Second-order phase dynamics}
Now, we use the same technique to construct the Poincar\'e map from the second-order phase dynamics equations. For this goal, we re-write Eqs.~\eqref{eqn:second_order_phase_reduction} in terms of phase differences and then obtain an equation for $\frac{\mathrm{d}\gamma_{13}}{\mathrm{d}\gamma_{12}}$ that is similar to Eq.~\eqref{eqn:first_order_gamma13/gamma12} but contains
additional terms proportional to $\widetilde{\varepsilon}^2$. Solving this equation by the perturbation technique (see Appendix~\ref{app:2nd_order} for details), we arrive at the following condition for RS:
\begin{equation}
\label{eqn:second_order_RS_condition}
    \left| \frac{{\varepsilon}^2 [\alpha-(\omega_1 - \omega_2)C_{12}]}{(\omega_1-\omega_3)(\omega_1-\omega_2)} \right| \geq \frac{1}{2}\;.
\end{equation}
This condition differs from the inequality \eqref{eqn:first_order_RS_condition}, derived in the first approximation, by the term $(\omega_1 - \omega_2)C_{12}$ alone.
(Notice that $C_{12}\approx C_{32}$.) This term is proportional to the amplitude of the synchronizing
terms $\sin(\varphi_3-\varphi_1)$, $\sin(\varphi_1-\varphi_3)$ in~Eqs.~\eqref{eqn:second_order_phase_reduction}.
These terms indicate the presence of an ``invisible'' coupling between oscillators 1 and 3. This coupling exists despite the absence of a physical link between the first and third units; the first-order phase reduction does not reveal it. Thus, RS is promoted by non-isochronicity and by indirect coupling through the hub.

\section{Results}
\label{sec:ref}

To validate our derivations, we compare the bifurcation diagram on the $\varepsilon$-$\alpha$ plane obtained using the various approximations against those obtained for the exact SL equations. Before discussing the plots, we briefly recall the approximations made and clarify the terminology used to distinguish between them. The results from the numerical computations using the SL system~\eqref{eqn:SLequations} will be referred to as ``exact''. If the numerical calculation used the first-order [Eq.~\eqref{eqn:first_order_phase_reduction}] or the second-order
[Eq.~\eqref{eqn:second_order_phase_reduction}] phase reduction, the corresponding result will be termed as ``NPR1'' or ``NPR2'', respectively~\footnote{We sweep the parameter space to determine the state of the system using efficient techniques. We find RS in the SL equations~\eqref{eqn:SLequations} by looking for a limit cycle solution where $\gamma_{12}$ is unbounded while $\gamma_{13}$ is bounded using a shooting method. To detect RS using the phase reduction equations exactly, we numerically construct the Poincar\'e map described in Sec. \ref{sec:RS_theoretical} by simulating Eq.~\eqref{eqn:first_order_gamma13/gamma12} (or its second-order counterpart) and check the presence of a fixed point. Regions of CS were computed using direct numerical simulations; we mark the points in the parameter space that resulted in $\Omega_1=\Omega_2=\Omega_3$ up to a tolerance of $10^{-4}$.}. Finally, the theoretical results obtained for the first-order [Eq. \eqref{eqn:first_order_RS_condition}] and second-order [Eq. \eqref{eqn:second_order_RS_condition}] phase reduction are coined as ``TPR1'' and ``TPR2'', respectively~\footnote{The borderline of the RS transition is given by the condition when inequalities (\ref{eqn:first_order_RS_condition},\ref{eqn:second_order_RS_condition}) turn to equalities.}.

As a first step, we compared the NPR1 and TPR1 borderlines of the RS transitions. We found that TPR1 very well reproduces the numerical results shown by dashed lines in Fig.~\ref{fig:compare_approximations_figures}.
This result confirms the capability of the perturbation approach to capture RS in the Kuramoto-Sakaguchi model \eqref{eqn:first_order_phase_reduction}.

\begin{figure}
    \includegraphics[width=\columnwidth]{Kumar_Fig3.png}
    \caption{\label{fig:secondorder_figures}
    (Color online) Comparison of theoretical and numerical results.
    Two-parameter bifurcation diagrams on the $\varepsilon$-$\alpha$ plane (coupling strength versus non-isochronicity) depicting the system's state. Exact domains of RS, CS, and asynchrony are shown in blue (dark gray), red (light gray), and white, respectively. The solid black line shows the RS borderline obtained numerically using the second-order phase reduction (NPR2). The dashed black line is the corresponding theoretical solution (TPR2). The oscillators' natural frequencies are (a) $\omega_1=1$, $\omega_2=\sqrt{2}$, $\omega_3=1.002$ and (b) $\omega_1=1$, $\omega_2=\sqrt{7}$, $\omega_3=1.01$.}
\end{figure}

Figure~\ref{fig:secondorder_figures} presents our main result. Here, we compare the NPR2 and TPR2 borderlines of the RS transition against the exact ones. When the frequency detuning $|\omega_1-\omega_3|$ is very small, as in Fig.~3a, all borders are practically identical for low coupling strengths. As the coupling strength $\varepsilon$ increases, the normalized coupling $\widetilde{\varepsilon}$ (see Eq.~\eqref{eqn:my_params}) is no longer small, which causes the observed deviation between the TPR2 and NPR2 borders. Note that the NPR2 border accurately reproduces the exact RS transition throughout the considered range of coupling strengths.
The bifurcation diagram for a second set of natural frequencies is presented in Fig.~3b. Again, for low coupling strengths, the agreement between the approximations and the exact solution is perfect. However, both NPR2 and TPR2 borders deviate from the exact border of the RS transition for higher values of coupling strength. This deviation occurs because $\varepsilon$ (and likewise $\widetilde{\varepsilon}$) are no longer small quantities. We mention in passing that the dynamics for higher coupling strengths is often not trivial. For instance, the transition to CS in Fig.~3b near the finger-like structure around the point $(\varepsilon=0.5,\alpha=0)$ exhibits complex, possibly chaotic, dynamics, presumably due to the effects of strong coupling. Interestingly, near this point, there exists a window of RS straddled by regions of CS on either side.

\section{Conclusions}
\label{sec:concl}
In summary, we analyzed the mechanisms of RS in a chain of three SL oscillators. We demonstrated that the RS transition is determined by the interplay of the non-isochronicity and the amplitude dynamics. The impact of the latter factor renders the standard first-order phase dynamics description of the RS phenomenon invalid. Our result emphasizes the importance of high-order phase reduction and highlights the crucial role amplitude dynamics may have in governing the behavior of networks of nonlinear oscillators.

We believe that the effect of the amplitude dynamics neglected in the first-order phase approximation and revealed by the high-order one holds for general limit-cycle oscillators.
This belief is supported by the results of numerical network reconstruction from data~\cite{Kralemann-Pikovsky-Rosenblum-11}, which demonstrated the emergence of coupling between indirectly interacting units. It will be interesting to investigate how the unit's complexity may bring about qualitatively new changes to the RS transition~\cite{lacerda2019multistable,minati2015remote,*karakaya2019fading} and if they can be explained under the present framework.

\begin{acknowledgments}
MK is grateful for the WISE scholarship by the DAAD (German Academic Exchange Service), which facilitated this work.
\end{acknowledgments}


\section{\label{intro} Introduction}
In recent studies~\cite{Kevin17,Hong18}, oscillators that sync and swarm, called {\it{swarmalators}}, have been considered, and interesting long-term states have been found. The swarmalators can move in the space under the correlated dynamics of the phase and space variables of each oscillator. The {\it{mobile}} feature of the oscillators is found to induce nonstationary states such as the active-phase wave state and the splintered phase wave state. The collective property of finite number of swarmalators is studied~\cite{Kevin18}, the implementation in robots is tried as an application~\cite{Bett20}, and various steady states by finite interaction-distance is reported~\cite{Lee21}. Mathematical property~\cite{swm2}, model extension~\cite{swm4,swm5,swm6}, and minimal modeling~\cite{swm7} are of basic interest. Application to swarming and flocking of biological colony~\cite{swm} is a primary interest. Control of engineering objects like drones or robots is also an application front~\cite{swm3}.

In the study of swarmalators, the interesting incoherent states are usually the consequence of negative-definite coupling strength in phase dynamics (the positive-definite case is also in Ref.~\cite{Hong18} where noise is added). On the other hand, the coupling characteristic among the constituents of many systems is rather complex and not so simple. For example, the interaction among the neurons in the neural network systems is given by the mixture with positive and negative ones~\cite{neuron,neuron2}, in general. As another example, Japanese frogs' calling behavior can be understood by considering the mixture of positive and negative interaction with each other~\cite{Ikkyu08}. Probably, the spin glass ~\cite{sg,spinglassXY}, to which the various interesting properties of condensed matter is attributed, can be the representative example for the non-definite sign of the couplings. Considering those features in nature, the mixed coupling with positive and negative strength deserves to be considered for understanding the collective behavior of the real systems, which motivates the present study.

In this paper, we consider a population of the oscillators that can sync and swarm, governed by the random interaction in the phase dynamics. We explore how the coupling-disorder affects the long-term states in the system. In particular, we pay attention to the possibility of the phase transition~\cite{ptStanely} in the system, and focus on whether the patterns observed in the absence of coupling disorder still appear. The effective annealing of the randomly quenched coupling strengths to their average is suggested in the mobility of swarmalators to understand the numerical results. Phase transition to the fully synchronized phase, reproducibility of the phases known in the original model of the no coupling disorder, and mixture of deformed long-term states by coupling disorders are explained from the viewpoint of the annealed couplings.

This paper is organized as follows: Section~\ref{mdl} introduces the model of coupled swarmalators with random couplings, and Sec.~\ref{pt} shows the collective behavior and the phase transition in the system. In Sec.~\ref{ls} we derive the threshold of the transition by the linear stability analysis of the fully synchronized state. Various long-term states including the nonstationary ones are shown in  Sec.~\ref{nonss}, and the system with couplings of more than one quenched randomness are understood as the combination of annealed-coupling systems in Sec.~\ref{mrc}. A brief summary is given in Sec.~\ref{sum}.

\section{\label{mdl} Model}
The generalized model of $N$-coupled oscillators that we consider here is given by
\begin{eqnarray}
\frac{d\theta_i}{dt} &=&
\frac{1}{N}\sum_{j\neq i} F(r_{ij}) K_{ij} \sin(\theta_j-\theta_i),
\label{eq:model_theta} \\
\frac{d{\bf{r}}_i}{dt} &=&
\frac{1}{N}\sum_{j\neq i}
G(r_{ij}, \theta_{ij})
\label{eq:model_r}
\end{eqnarray}
for $i=1,\cdots,N$, where $\theta_i$ and ${\bf{r}}_i$ represent the phase $(0 \leq \theta_i \leq 2\pi)$ and the position vector of the $i$th oscillator, respectively. $F(r_{ij})$ is a function for the spatial dynamics of the oscillators with $r_{ij}=|{\bf{r}}_{ij}|$ for ${\bf{r}}_{ij}\equiv{\bf{r}}_j - {\bf{r}}_i$, and $G(r_{ij},\theta_{ij})$ is a function for spatial dynamics with $\theta_{ij}\equiv\theta_j - \theta_i$. $K_{ij}$ denotes the random coupling strength between the oscillators $i$ and $j$, having the symmetric property $K_{ij}=K_{ji}$. It represents a sort of {\it{bond}} coupling between the oscillators, instead of the {\it{site}} one ($K_i$) that affects the oscillator itself. For simplicity, we randomly chose $K_{ij}$ from the two-peaks distribution
\begin{equation}
h(K_{ij})=p\delta(K_{ij}-K_p)+(1-p)\delta(K_{ij}-K_n),
\label{eq:hK}
\end{equation}
where $p$ is the probability of the positive coupling, and $K_{p} > 0$ and $K_{n} < 0$. For convenience, the ratio $Q\equiv -K_{n}/K_{p} (>0)$ has been chosen for the control parameter. With the functions $F$ and $G$ in the model, the phase variable $\theta$ and the spatial one $r$ become correlated, which means the oscillators can move around in the space. Such correlation in the dynamics of the space and the phase of the oscillators has been also studied in Ref.~\cite{Kevin17}, where the authors found five long-term states including the two nonstationary ones. Note that the phase coupling in Ref.~\cite{Kevin17} does not have any disorder. We here consider coupling disorder that is given by the random values from the distribution in Eq.~\eqref{eq:hK}.

We notice that the special case with $F(r_{ij})=1$ and $G(r_{ij},\theta_{ij})=0$ is equivalent to the mean-field $XY$ model~\cite{spinglassXY} with {\it{random}} coupling strength, governed by the Hamiltonian ${\cal{H}}=-\frac{1}{2N}\sum_{i\neq j} K_{ij}\cos(\theta_j - \theta_i)$. Equation~(\ref{eq:model_theta}) with $F=1$ and $G=0$ leads to the overdamped version of the Hamiltonian dynamics at zero temperature, and interestingly it is found that the first-order phase transition occurs at the same threshold $p_c$ as that we study here~\cite{Hong21}.

In this paper, we consider $F=1/r_{ij}$. The function $G$ consists of the attraction and repulsion forces acting on each oscillator, where the force functions are taken as the algebraic ones with a power like $r_{ij}^{-a}$ following Ref.~\cite{Kevin17}. The model is then given by
\begin{eqnarray}
{\dot{\theta}}_{i} &=& \frac{1}{N}
\sum_{j\neq i}
\frac{K_{ij}}{r_{ij}}\sin(\theta_j-\theta_i), \label{eq:theta} \\
\dot{\bf{r}}_i &=& \frac{1}{N}\sum_{j\neq i}
\left[\frac{{\bf{r}}_{ij}}{r_{ij}}(A+J\cos(\theta_j-\theta_i)) - B \frac{{\bf{r}}_{ij}}{r^2_{ij}}\right], \label{eq:x}
\end{eqnarray}
where $A$ and $B$ are respectively the parameters for the attractive force and the repulsive force. Here, we choose $A=B=1$ for convenience. $J$ is the parameter which measures how the phase similarity enhances the spatial proximity.  For example, a positive value of $J$ means {\textit{like attracts like}}, i.e., the swarmalators tend to be near the other swarmalators with similar phases. On the other hand, a negative value of $J$ means the opposite case: The swarmalators prefer to be near the others with opposite phases. To keep the attractive force always positive, the values of $J$ are constrained like $-1 \leq J \leq 1$, following the Ref.~\cite{Kevin17}. The difference from the model of Ref.~\cite{Kevin17} is the generalization of the constant $K$ to pair-dependent $K_{ij}$. For $K_{ij} > 0$, the swarmalators prefer to have the similar phase, but for $K_{ij} < 0$, the opposite tendency occurs. We are interested in the effect of random $K_{ij}$ in the formation of long-term states.

\section{\label{pt} Phase transition}
\begin{figure}
\includegraphics[width=9.0cm]{Fig1a.eps}
\includegraphics[width=9.0cm]{Fig1b.eps}
\caption{(Color Online)
(a) $R$ is shown as a function of $p$ for various values of $J$; (b) $S$ vs. $p$ for various values of $J$. Parameters: $Q=1.25$ and $N=200$.
}
\label{fig:RS}
\end{figure}
To see the collective states in the system and the effects of random coupling on the states we perform the numerical simulations on Eqs.~(\ref{eq:theta}) and (\ref{eq:x}). We first investigate the phase synchronization behavior of the oscillators, by measuring the complex order parameter defined by~\cite{Kuramoto84}
\begin{equation}
Re^{i\Theta} = \frac{1}{N} \sum_{j=1}^N e^{i\theta_j},
\label{eq:R}
\end{equation}
where $R$ measures the phase coherence, and $\Theta$ is the mean phase angle of the oscillators. For example, $R=0$ means the incoherent state where the all oscillators have the random phase $\theta_i\in[0,2\pi)$, and $R=1$ is the fully synchronized state where they all have the same one $\theta_j=\theta_s$ for all $j$'s.

Also, to see the correlation between the phase angle $\theta_j$ and the azimutal one $\phi_j=\tan^{-1}(y_j/x_j)$ for the $j$th oscillator, which is induced by the mobile feature of the oscillators in the system, we measure another order parameter defined by~\cite{Kevin17}
\begin{equation}
S_{\pm}e^{i\Psi_{\pm}} = \frac{1}{N} \sum_{j=1}^N e^{i(\phi_j \pm \theta_j)},
\label{eq:S}
\end{equation}
where $S_{\pm}$ is the magnitude, and $\Psi_{\pm}$ is its mean phase. Here, the order parameter $S$ measures the ``correlation'' between the spatial information (by the $\phi_j$) and the phase information (by the $\theta_j$). The system exhibits ``assortative and/or commutative'' correlation between the two depending on the initial conditions, where the $S_{+}$ comes from the assortative correlation, and $S_{-}$ comes from the commutative correlation. We chose $S$ by taking the maximum value of $S_{+}$ and $S_{-}$, i.e., $S={\rm{max}(S_{\pm})}$~\cite{Kevin17}.

In order to examine $R$ and $S$, we basically performed the numerical simulations using the Python or C programs. The total $5\times 10^4$ time steps with the discrete time unit $dt=0.1$ have been considered, where the first $2.5\times 10^4$ steps were discarded for the equilibrium, and then the time average of the quantities has been done for the later steps. And, the ten samples have been used for the average. Initial phases and positions were, respectively, randomly sampled from the uniform distributions over $[0,2\pi)$ and $[-1,1)\times[-1,1)$.

Figure~\ref{fig:RS} shows the behavior of $R$ and $S$ as a function of $p$ for various values of $J$. The value of $Q=1.25$ was chosen, and $K_p=1$ and $K_n=-1.25$ were chosen for convenience. We find that the system shows the phase transition behavior from the incoherent state ($R=0$) to the fully synchronized one ($R=1$) at a finite value of $p$, i.e., $p_c~\sim 0.6$. Interestingly, we find that the behavior of $R$ and its threshold $p_c$ do not depend on $J$, which can be easily understood by the absence of $J$ in the phase dynamics given by Eq.~(\ref{eq:theta}). On the other hand, the behavior of $S$ is different from that of $R$ as shown in Fig.~\ref{fig:RS}: In the regime of the fully synchronized state ($R=1$) for $p> p_c$, the order parameter $S$ is zero, which implies no correlation between the phase angle and the spatial angle. However, in the incoherent state with $R=0$ for $p<p_c$, $S$ shows finite value, which means there is some correlation between the phase angle and the spatial angle. Moreover, the value of $S$ for $p<p_c$ shows $S=1$ for a certain region of $p$, and the range of $p$ of this region depends on the value of $J$ as shown in Fig.~\ref{fig:RS}~(b): The range gets wider as $J$ increases.

To pin down the threshold $p_c$, we now estimate it numerically. We measure the value of $p$ at which the first derivative $dR/dp$ shows the maximum. Figure~\ref{fig:pc_Q} shows the behavior of $dR/dp$ as a function of $p$ for various values of $Q$, whereas the value of $p$ at which the $dR/dp$ reaches the maximum is denoted as $p^*$. The inset of Fig.~\ref{fig:pc_Q} shows $p^*$ as a function of $Q$, which shows a good consistency with the theoretical prediction $p_c=Q/(1+Q)$ (see Sec.~\ref{ls}).
\begin{figure}
\includegraphics[width=9cm]{Fig2.eps}
\caption{(Color Online)
The first derivative of $R$ over $p$, $dR/dp$, is shown as a function of $p$ for various values of $Q$, where the value of $p$ at which $dR/dp$ reaches the maximum value is denoted as $p^*$. Inset: $p^*$ vs. $Q$ is shown, where the blue solid line is the theoretical prediction given by $p_c=Q/(1+Q)$, and the value of $p^*$ is consistent with the prediction. The system size is $N=800$ and $J=0.8$, and the data have been averaged over 10 samples.
}
\label{fig:pc_Q}
\end{figure}

\section{\label{ls} Linear stability of the fully synchronized state ($R=1$)}
In this section, we investigate the linear stability of the fully synchronized state with $R=1$. Let a fully synchronized state have a common phase $\theta_i=\theta_{\rm s}$ for all $i$. In order to examine its linear stability, we consider a slightly perturbed situation such as $\theta_i =\theta_{\rm s} + \phi_i$ with small perturbation $\phi_i$ for all $i$. For this setting, Eq.~\eqref{eq:theta} is linearized to
\begin{equation}
\label{lt}
\dot \phi_i
= \frac{1}{N}\sum_{j\neq i} \frac{K_{ij}}{r_{ij}}(\phi_j-\phi_i)\,.
\end{equation}
Separating the summation according to the values of $K_{ij}$, we rewrite Eq.~\eqref{lt} as
\begin{eqnarray}
\label{lt2}
\dot \phi_i = \frac{1}{N}
&&\left( K_p \sum_{j\in C_p(i)} \frac{1}{r_{ij}}(\phi_j-\phi_i) \right.
\nonumber \\
&& \left.~~~+K_n \sum_{j\in C_n(i)} \frac{1}{r_{ij}}(\phi_j-\phi_i)\right),
\end{eqnarray}
where $C_p(i)$ is the collection of $j$ having $K_{ij}=K_p$ and $C_n(i)$ is that for $K_{ij}=K_n$.

We here note that $K_{ij}$ is not included in the spatial dynamics [see Eq.~\eqref{eq:x}]. Thus, for a fixed $i$, the fraction of such $j$'s giving $K_{ij}=K_p$ in the total $N$ would be approximately $p$, where the deviation is expected to decrease as $N$ increases. Similarly, that for $K_{ij}=K_n$ is given by $(1-p)$. This consideration motivates us to replace $\sum_{j\in C_p(i)}$ and $\sum_{j\in C_n(i)}$ with $p\sum_{j(\neq i)}$ and $(1-p)\sum_{j(\neq i)}$, respectively, in Eq.~\eqref{lt2} when $N$ is large enough. Therefore, for large $N$, Eq.~\eqref{lt2} can be cast into
\begin{equation}
\label{lt3}
\dot \phi_i =
\langle K_{ij}\rangle
\frac{1}{N}\sum_{j\neq i} \frac{1}{r_{ij}}(\phi_j-\phi_i)\,.
\end{equation}
where $\langle K_{ij}\rangle\equiv pK_p +(1-p)K_n $ is the average of $K_{ij}$ for the two-peak distribution $h(K_{ij})$.

Introducing $m_{ij} = -1/(Nr_{ij})$ for $i\neq j$ and $m_{ii}=-\sum_{j\neq i}m_{ij}$, one writes Eq.~\eqref{lt3} as
\begin{equation}
\label{lt4}
\dot \phi_i =
-\langle K_{ij}\rangle \sum_{j=1}^N m_{ij}\phi_j\,.
\end{equation}
Interestingly, the matrix $m=\{m_{ij}\}$ is positive-definite for large $N$ as follows. It holds $|m_{ij}| \sim {\cal O}(1/\sqrt{N})$ for off-diagonal terms. This is based on the numerical observation that all swarmalators reside in the two-dimensional spatial region of the ${\cal O}(1)$-area and there is no concentration of them. Thus, $r_{ij}\sim 1/\sqrt{N}$ holds, from which $|m_{ij}| \sim {\cal O}(1/\sqrt{N})$ follows. It also follows that $m_{ii}=-\sum_{j\neq i}m_{ij} \sim {\cal O}(\sqrt{N})$. Then, as $N$ increases, $m=\{m_{ij}\}$ approaches $\bar m = \{\bar m_{ij}\}$, which is given by
\begin{equation}
	\label{m-infty}
	\bar m_{ij}=\delta_{ij}\bar m_i\,,
\end{equation}
where $\delta_{ij}$ is the Kronecker delta that assigns $1$ when it assigns $i=j$ while $0$ otherwise. The matrix $\bar m$ is positive-definite because $\bar m_i\equiv\lim_{N\to\infty}m_{ii}>0$. Therefore, $m$ is also positive-definite for sufficiently large $N$.

When the positive definiteness of $m$ is considered in Eq.~\eqref{lt4}, one knows that the sign of $\langle K_{ij}\rangle$ determines whether $\phi_i$ will decay or not. This is because every eigenvalue of the positive-definite matrix is positive. Thus, $\phi_i$ decays when $\langle K_{ij}\rangle>0$, otherwise it grows when $\langle K_{ij}\rangle<0$. Since $Q\equiv -K_n/K_p$, we get $\langle K_{ij}\rangle =  K_p(p+(p-1)Q)$. Therefore, $\langle K_{ij}\rangle>0$ corresponds to $p>Q/(1+Q)$ and $\langle K_{ij}\rangle<0$ corresponds to $p<Q/(1+Q)$, which leads to
\begin{equation}
\label{pc}
p_c=Q/(1+Q)\,.
\end{equation}
For $Q=1.25$, used for the numerical data shown in Figs.~\ref{fig:RS} and \ref{fig:pc_Q}, Eq.~\eqref{pc} gives $p_c=5/9=0.55..\,$. The deviation from the numerical value $p_c\approx 0.6$ in the figures is supposed to be a finite $N$ effect. Note $p_c\approx 0.565$ of Fig.~\ref{fig:pc_Q} for $N=800$ is closer to the prediction $p_c = 5/9$ than the prediction in Fig.~\ref{fig:RS} for $N=200$. It is interesting that the threshold value shown in Eq.~(\ref{pc}) is equivalent to that for the mean-field $XY$-type oscillators with random coupling disorder~\cite{Hong21}. This means that the spatial dynamics in the system does not change the threshold value of the phase transition from the incoherent state to the fully synchronized one.

The crucial step of our argument above is the use of $\sum_{j\in C_p(i)}=p\sum_{j(\neq i)}$ and $\sum_{j\in C_n(i)}=(1-p)\sum_{j(\neq i)}$ to write Eq.~\eqref{lt3}. We here remark this is basically an annealed approximation. In fact, $K_{ij}$ is quenched by definition in the model and has never been treated in an annealed way in the numerical work. Thus the consistency between Eq.~\eqref{pc} and its numerical counterpart is interesting. Equation~\eqref{eq:x} gives a clue to understanding how the quenched $K_{ij}$ works as if annealed. The spatial dynamics governed by Eq.~\eqref{eq:x} is indifferent to $K_{ij}$. Based on this, one may expect the distribution of $K_{ij}$ for given $i$ has no position-dependence, i.e., the fractions for $K_p$ and $K_n$ are, respectively, $p$ and $(1-p)$ independent of position. This results in the annealed average of phase dynamics in Eq.~\eqref{eq:theta}. The consistency between the numerical data and the annealed average, even in the presence of quenched $K_{ij}$, means that the oscillators move to the proper places. That is, the effective annealing of $K_{ij}$ to $\langle K_{ij}\rangle$ is a characteristic property of swarmalator that can move.

We also examine that the sign of $\langle K_{ij}\rangle$ always governs the onset of the sync state in further numerical study using other $h(K_{ij})$ distributions. These observations suggest that the annealed approximation could be generally available. We remark that the effective annealing requiring proper locations of swarmalators is considered for long-term states.

\section{\label{nonss} Non-stationary states in the incoherent regime ($R=0$)}
\begin{figure}
\includegraphics[width=9cm, height=8.5cm]{Fig3.eps}
\caption{(Color Online) Behavior of $S$, $R$, $V$, and $U$ is shown as a function of $p$ for (a) $J=0.8$ and (b) $J=0.2$. Snapshots for the corresponding regime are also shown in the insets.  Parameters: $N=200$ and $Q=0.8$.}
\label{fig:RSUV}
\end{figure}
In this section, we investigate the incoherent state $(R=0)$ further. In particular, we pay attention to the possibility of the nonstationary states such as the active phase wave (APW) state and the splintered phase wave (SPW) state found in Ref.~\cite{Kevin17}. Therein, the all coupling strength among the oscillators was chosen as the negative one with no disorder like $K_{ij}=K$, which corresponds to $p=0$ in the current study. We note that, in Ref.~\cite{Kevin17}, various long-term states such as the async state, the APW state, the SPW state, and the static phase wave state were found depending on the strength of the negative coupling in the phase dynamics. Especially, the APW state and the SPW state are the nonstationary states where the oscillators are not static in both phase and space. It is thus interesting whether the dynamics of Eqs.~\eqref{eq:theta} and \eqref{eq:x} with quenched disorder $K_{ij}$ can induce such a nonstationary feature. Below, we again resort to the numerical data by the same scheme used for Fig.~\ref{fig:RS}.

To see the possibility of the nonstationary states for $p>0$, we measure the following two quantities. One is the mean velocity $V$ defined by~\cite{Kevin17,Hong18}
\begin{equation}
V = \frac{1}{N}\sum_{j=1}^N v_i,
\label{eq:V}
\end{equation}
where $v_i = \sqrt{{{\dot{x}}_i}^2+{{\dot{y}}_i}^2}$. The mean velocity measures the nonstatic property of the state. In other words, the finite value of $V$ means that the swarmalators move around in the plane in the long-term state. The other one is~\cite{Kevin17,Hong18}
\begin{equation}
U = \frac{N_{\rm{rot}}}{N},
\label{eq:U}
\end{equation}
where $N_{\rm{rot}}$ is the number of swarmalators circulating at least $2\pi$, and thus $U$ represents the fraction over the total swarmalators~\cite{Kevin17,Hong18}. Note that $U$ can play as the role of the indicator for the active phase wave state~\cite{Kevin17}. Figure~\ref{fig:RSUV}(a) and (b) show the quantities $R$, $S$, $V$, and $U$ as a function of $p$ for $J=0.8$ and $0.2$, respectively. As $p$ increases from zero, we find that the phase transition occurs from the incoherent state $(R=0)$ to the fully synchronized state ($R=1$) at a finite value of $p_c$, as shown in the behavior of $R$. Interestingly, the value of $p_c$ for two different values of $J$ are found to be same, which means the threshold does not depend on $J$. On the other hand, the behavior of $S$, $U$, and $V$ is found to differ depending on the value of $J$ (See Fig.~\ref{fig:RSUV}).

When $R=0$, note that the system shows a finite value of $S$ through the splintered and static phase wave. Also, the system shows the active phase wave state for small $p$ when $J$ is large.  This active phase wave does not occur when $J$ is small. We find that the system shows both the stationary states and the nonstationary states depending on $p$ for a given value of $Q$ and $J$.  Especially, when $J$ is large $(J=0.8)$, all five long-term states that have been reported in Ref.~\cite{Kevin17} exist depending on $p$.  It is interesting to note that the nonstationary states such as active phase wave sate and splintered phase wave also exist in the presence of the positive coupling among the swarmalators in the system ($p>0$).

\section{\label{mrc} Mixture of random couplings}
Above, we have seen that the result by random $K_{ij}$ is comparable with that by constant $K=K_{ij}$. The long-term state patterns, including the fully synchronized one, from both settings are similar for large $J$. It is remarkable that such analogy holds in the presence of finite faction for negative $K_{ij}$ that may bring about the frustration~\cite{frst,frst2} in the phase dynamics. This strongly suggests the annealed property conjectured for the onset of the sync phase at $\langle K_{ij} \rangle =0^+$ is also valid for non-zero $\langle K_{ij} \rangle$ cases. This motivates us to test the long-term states in the system with such quenched random couplings by mixture of different randomnesses. We are interested in whether the characters of long-term states attributed to each randomnesses will last or not.

Consider $N_1$ number of swarmalators in total $N$, and call them group 1 (G1). The coupling strength $K_{ij}$, therein, is assigned with $1$ for probability $p_1$ or $-1$ for $1-p_1$. Next, consider the other group, G2, of $N_2(=N-N_1)$ swarmalators, where $K_{ij}=\pm 1$ is similarly assigned using probability $p_2$. We also consider coupling between the swarmalators in both groups as assigning $1$ to $K_{ij}$ for probability $p_{\rm m}$ or $-1$ for $1-p_{\rm m}$. Below, we use $p_{\rm m}=p_2$ for simplicity. When $(N_1,p_1)=(100,1)$ and $(N_2,p_2)=(100,0)$, one may expect the sync disc in G1 and the APW state in G2 if the couplings that work between the groups do not overwhelm each character. Note the average coupling in G2, $\langle K_{ij} \rangle_2=-1$, is in the region for the APW state in the phase diagram reported in Ref.~\cite{Kevin17}. We obtained numerical data with time unit $dt=0.01$ and consider $2\times 10^5$ time steps.

Figure~\ref{fig-mix}(a) shows the numerical result.
\begin{figure}
\includegraphics[width=9.0cm]{mix-1a.eps}
\includegraphics[width=9.0cm]{mix-1b.eps}
\caption{(Color Online)
(a) Long-term state pattern for $N_1=N_2=100$, $p_1=1$, $p_2=0$, $p_{\rm m}=p_2$, and $J=1$. The inset is the trajectory of an arbitrarily chosen swarmalator in the U-shaped region, which sweeps the region. (b) Long-term state pattern for the same setting used in panel (a) except $p_2=0.4$. The black solid curve is the trajectory of a swarmalator localized there.
}
\label{fig-mix}
\end{figure}
Two patterns appear; one is a slightly deformed sync disc while the other is U-shaped one. We have observed the former vibrates a little (not demonstrated here) and the latter is an active phase state. The inset in Fig.\ref{fig-mix}(a) is the trajectory of a swarm in the U-shape. Since it sweeps the region while changing the phase, we name the state the deformed APW (dAPW). It is interesting to note that this reminds us of the ``chimera state'' observed in the oscillators with identical frequency~\cite{chimera,chimera2,chimera3,chimera4}: The partial group of the swarmalators in the system exhibits the sync cluster, but the other group of the swarmalators do not join the sync cluster, showing the U-shape APW.

We vary $N_1$ to see a possible change in pattern. One easily see the distance between the sync disc and the dAPW increases when comparing Figs.~\ref{fig-mix}(a) and \ref{fig-mix2}(a). We also see the sync disc becomes more circular and static when the distance increases.
\begin{figure}
\includegraphics[width=9cm]{mix-2.eps}
\caption{(Color Online)
Long-term state patterns composed of sync disc and dAPW for $N_1=60$, $N_2=140$, $p_1=1$, $p_2=0$ in (a) and for $N_1=140$, $N_2=60$, $p_1=0.7$, $p_2=0$ in (b).
}
\label{fig-mix2}
\end{figure}
We next test whether a sync disc by $p_1<1$ in G1 only system ($N_2=0$) will still last after G2 is involved ($N_2>0$). The interest is a validity of the annealed approximation of $K_{ij}$ in G1. Figure~\ref{fig-mix2}(b) is the result for $p_1=0.7$, which is also composed of the sync disc and the dAPW. The sync disc is almost circular and static with enough separation from the dAPW.

We finally test whether G2 can show the SPW state while G1 remains a sync disc. As the splintered state appears for the negative coupling larger than that for APW, we use $p_2=0.4$, which corresponds to $\langle K_{ij}\rangle_2 = -0.2$, the average coupling strength in G2. This value is in the region for the splintered wave in the phase diagram reported in Ref.~\cite{Kevin17}. Figure~\ref{fig-mix}(b) is the result when the other settings are same as those for Fig.~\ref{fig-mix}(a). There still appear sync disc and U-shaped region. This time, the swarmalators of different colors in the U-shaped region are not mixed, different from the dAPW case. The gradual change of the swarmalators' colors in space is a characteristic of the splintered phase wave. The black solid curve near the bottom is the trajectory of a swarm, which remains near there as time goes on. As shown, it does not sweep the U region but localized in the small area while splintering. This is always the case for all swarmalators in the U-shaped region. We thus regard the U-shaped pattern as a deformed splintered phase wave (dSPW).

In the above numerical tests, the characters by each random couplings are preserved more or less in the mixed system. There apparently appear sync disc, APW, and splintered state though deformed. We here add the result is qualitatively similar for the various values of $p_{\rm m}$ (not shown here), and this does not depend on the initial condition. These observations suggest that the long-term states, demonstrated in Secs~\ref{ls} and \ref{nonss}, are still robust more or less in the mixed system. Interestingly, the patterns survive the frustration taking place between the subsystems G1 and G2. A mixture with async state or of more than two groups was not tested this time. Further study on various interesting properties of the mixed deformed patterns will appear elsewhere.

\section{\label{sum} Summary}
We considered the population of swarmalators with random coupling strength, and explored how the coupling disorder affects the long-term states in the system. In particular, the possibility of the phase transition and the robustness of the state patterns are focused. To understand the long-term states observed in the system with quenched disorder of coupling strength, we considered the effective annealed approximation, which is mediated by the mobility of swarmalators. In the viewpoint of annealed couplings, the numerical observation in the quenched system is explained and, furthermore, such a system of the mixture of different quenched disorders is also understood.

We found that the system shows the phase transition from the incoherent state to the fully synchronized state at a certain threshold $p_c$, where the value of $p_c$ is argued in the linear stability analysis of the fully synchronized state. We also found that, in the regime of the incoherent state below the threshold, various long-term states are found to exist. Especially, the nonstationary states such as the dAPW and the dSPW besides the normal APW and SPW are discovered.

All long-term states known for the bare model in Ref.~\cite{Kevin17} are realized in the system of random coupling strength. The values of average-coupling where each state appears are similar to their counterparts in the bare model. This strongly suggests the randomly quenched coupling strengths work as if annealed. The pattern by the supposed annealed coupling is so robust that the mixture of different random couplings leads to the proper combination of each deformed patterns.

\section{Acknowledgements}
This research was supported by NRF Grant No. 2021R1A2B5B01001951 and `Research Base Construction Fund Support Program' funded by Jeonbuk National University in 2021 (H.H), and by NRF Grant No. 2018R1D1A1B07049254 (H.K.L.).