\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}

