\section{Introduction}

Compact stars serve as excellent astrophysical laboratories for the study of dense nuclear matter \cite{Haensel-2007,Zeldovich-71,Weinberg-72,Schutz-85,Bielich-2020}. The systematic study of pulsars and the detection of gravitational waves (GW) have already yielded significant constraints on the nuclear equation of state (EOS) \cite{Demorest-2010,Fonceca-2016,Arzoumanian-2018,Antoniadis-2013,Cromartie-2020,Linares-2018,Abbott-2017,Abbott-2018,Abbott-2019}. A question that still remains unanswered concerns the relevant degrees of freedom up to densities appearing in neutron star cores \cite{Heiselberg-2000,Weber-2005}. Compact stars could be purely hadronic, but the very dense environment indicates the possible existence of exotic forms of matter such as deconfined quarks. The latter opens up new scenarios that predict strange quark stars, composed purely of strange quark matter, or hybrid stars where a quark core is surrounded by a mantle of hadronic matter. In practice, the distinction between neutron, strange and hybrid stars is not an easy task as their radius around the observed mass region of 1.4 $M_\odot$ is rather similar. Alternative approaches that may assist identifying the phases of nuclear matter within compact stars include the study of their thermal evolution \cite{Carvalho-2016,Lyra-2023}, binary neutron star mergers \cite{Prakash-2021,Bauswein-2020,Most-2020} and phenomena related to vibration or rotation \cite{Madsen-1992,Madsen-1998,Madsen-2000,Steiner-2008,Jaikumar-2008,Alford-2010,Ranea-2018,Bratton-2022}.

The construction of hybrid EOSs often requires to describe the hadronic and quark phases separately. Depending on the dynamics of the phase transition and mainly on the speed of sound structure in quark matter a third family of compact objects may appear in the mass-radius plane. The aforementioned family of compact stars gives rise to the existence of twin stars, i.e. stars with identical mass but fairly different radius \cite{Gerlach-1968,Kampfer-1981a,Kampfer-1981b,Glendenning-2000,Schertler-2000}. Recently, the scenario of twin stars drew a lot attention, mainly because of the discovery of GW and thus the possibility of detecting them \cite{Blaschke-2013,Alvarez-Castillo,Montana-2019,Tsaloukidis-2023,Alford-2013}. Note that identifying twins would be smoking gun evidence of hadron-quark phase transtion in compact stars. In a recent study, Lyra et al. \cite{Lyra-2023} investigated the impact of compactness on the cooling of twin pairs, finding that only stars with significantly different radius exhibit considerable deviations in their thermal evolution. Furthermore, Landry and Chakravarti \cite{Landry-2022} argued on the possibility of distinguishing twins with next-generation GW detectors, through their tidal deformabilities. In the present work we study for the first time the deviation of the r-mode instability windows of twin stars and hence the differences that appear in their rotational limits.

It is well-established that relativistic stars may suffer a number of different instabilities. Among them, the r-mode instability (rotational mode) has been proposed as an explanation for the fact that neutron stars do not spin up to the theoretically allowed limit known as Kepler frequency~\cite{Andersson-1998,Friedman-1998,Friedman-1999,Andersson-2001,Andersson-2003,Kokkotas-2003,Andersson-2001-b,Lindblom-2001}. The r-modes are oscillations appearing in rotating stars, and their restoring force is the Coriolis force. In principle, the r-mode instability can only take place if the gravitational-radiation driving timescale is shorter compared to the timescales of the various dissipation mechanisms that may occur in the neutron star interior. By equating the driving and damping timescales one obtains the so-called r-mode instability window, which defines a critical frequency (maximum spin frequency for stable r-modes) as a function of temperature.

In the past decades there has been an extensive study of the r-modes (and numerous other types of oscillation) due to the possible detection of their GW \cite{Anderson-1999,Reisenegger-2003}. There are several studies predicting that accreting stars in low mass X-ray binaries (LMXBs) may be subject to long-lasting r-modes. In particular, compact stars containing exotic matter, such as deconfined quarks or hyperons, may be persistent sources of GW emission \cite{Anderson-2002,Nayyar-2006}. In addition, some authors \cite{Chugunov-2014} argue for the existence of a large unobserved population of quiescent (post-accretion) LMXBs characterised by long-lived ($\sim$10$^9$ yr) r-mode emission. Specifically, Chugunov et al. \cite{Chugunov-2014} suggested the existence of a new class of neutron stars, the so-called HOFNARS (HOt and Fast Non Accreting Rotators). Such stars retain a high temperature due to heating associated with unstable r-modes. Following the discovery of gravitational radiation from binary neutron star mergers, the search for GW signals associated with r-modes has started \cite{Caride-2019,Rajbhandari-2021,Covas-2022}. It is notable that, the absence of a detection so far has provided the opportunity to set upper limits on the GW emission and the r-mode saturation amplitude \cite{Rajbhandari-2021,Covas-2022}.

It has been shown that, the r-mode instability window of purely neutron stars is very wide to be compatible with current LMXBs data (assuming that all observed stars are stable with respect to r-modes, i.e. there are no HOFNARS). Specifically, a very strong dissipation mechanism, such as a perfectly rigid crust, is essential for the stabilization of r-modes. Numerous studies have attempted to treat this problem by considering the presence of exotic degrees of freedom in compact star cores \cite{Haskell-2012,Ofengeim-2019a,Ofengeim-2019b,Alford-2014,Pan-2006,Zheng-2006}. In particular, it has been shown that the bulk viscosity of hyperon or deconfined quark matter may be sufficient to stabilize r-modes for the frequencies and temperatures of the observed pulsars. However, it is important to comment that hyperons are expected to appear in densities of 2-3 $n_0$ (where $n_0=$ 0.16 fm$^{-3}$ is the nuclear saturation density). Thus, the fraction of the core where hyperons are present, and hence the effective damping due to their viscosity, is limited in low mass neutron stars. Subsequently, the fastest rotating pulsars can only be explained if they are massive enough \cite{Ofengeim-2019a,Ofengeim-2019b}. Similarly, it has been shown that the width of the r-mode instability window of hybrid stars is determined mainly by the amount of quark matter in the core \cite{Jaikumar-2008,Alford-2012}.

In Ref.~\cite{Papazoglou-2016}, the authors employed a set of analytical solutions of the Tolman-Oppenheimer-Volkov (TOV) equations in order to study the influence of neutron star bulk properties on the r-modes. They found that the instability window is quite sensitive to the radius of a star~\cite{Papazoglou-2016}. The latter leads to the conclusion that if twin stars do exist, their instability windows would deviate due to their radius difference. In addition, taking into account that the relevant degrees of freedom are different in the center of the two twins, the damping mechanisms that suppress the growth of the r-mode instability (bulk and shear viscosities) are going to be different as well~\cite{Madsen-1992}. This opens up a new intriguing scenario where two stars with identical mass, and similar rotational frequency and temperature profiles, may behave differently with respect to r-modes. In particular, if we assume that a star having angular velocity $\Omega_i$ and temperature $T_i$ is stable with respect to the r-modes, then any other (same mass) star with similar temperature and $\Omega\leq\Omega_i$ should be stable as well. However, this is not necessarily the case if the two stars are twins since their instability windows are expected to be different. Thus, the future detection of GW due to unstable r-modes, from multiple sources, may allow us to identify a third family of compact objects.

The motivation of the present study is twofold. Firstly, we wish to systematically study the parameters (energy density gap, crust elasticity, transition density) that affect the deviation between the r-mode instability windows of twin stars. In addition, we wish to clarify how these parameters affect the differences that appear in the spin down evolution (due to unstable r-modes) of twins. Secondly, we wish to examine if EOSs that predict a third family of compact objects are a viable solution for the explanation of current LMXBs data.

This paper is organized as it follows. Section \ref{2} is devoted to the presentation of the hadronic models employed in this work and the construction of hybrid EOSs that predict twin star configurations. In Section~\ref{3} we discuss in detail the r-mode instability formalism, while in Section~\ref{4} we present a simplified model for the spin-down of compact stars (due to unstable r-modes). In Section~\ref{5} we present our results and discuss their implications. Section~\ref{6} contains a summary of our findings.

\section{Hadron-quark phase transition} \label{2}

A hybrid EOS often results from the combination of a low density hadronic model and a high density quark EOS. The key ingredient for the construction is the matching process between the two phases. In particular, there are two widely employed methods in order to obtain hybrid EOSs: a) the Maxwell construction and b) the Gibbs construction. The main difference of the aforementioned approaches is the number of charges that are globally conserved in the system~\cite{Glendenning-1992}.~In the former case, the phase transition is abrupt (i.e.~the two phases are separate), while in the latter scenario a mixed phase is present.

In the present work we adopt the Maxwell construction for the description of the phase transition. This particular approach is the favored one, in the scenario where the surface tension $\sigma_s$ in the hadron-quark crossover, is larger than the critical value of $\mathrm{\sim 40 \;MeV\;fm^{-3}}$ and less than the highest allowed one of $\mathrm{\sim 100 \; MeV\;fm^{-3}}$, according to lattice QCD calculations~\cite{Mariani-2017}. In this case the phase transition is sharp resulting in a discontinuity in the energy density. Specifically, the energy density reads~\cite{Blaschke-2013,Alvarez-Castillo,Montana-2019,Alford-2013}
\begin{equation}
  \mathcal{E}(P) = \begin{cases}
      \mathcal{E}_{\rm HADRON}(P), & P\leq P_{\rm tr} \\
      \mathcal{E}(P_{{\rm tr}}) + \Delta \mathcal{E} + (c_s/c)^{-2}(P-P_{{\rm tr}}), & P > P_{{\rm tr}}.
   \end{cases}
   \label{1}
\end{equation}
 where $P$ stands for the pressure, $c_s$ is the speed of sound and $c$ is the speed of light. Furthermore, $P_{tr}$ and $\Delta\mathcal{E}$ denote the transition pressure and the energy density jump, respectively. It is important to comment that the first line of Eq.~(\ref{1}) refers to the hadronic phase while the second one to the quark model. We treat the quark phase using a phenomenological approach known as constant speed of sound (CSS) parametrization~\cite{Montana-2019,Alford-2013}. More precisely, the second line of Eq.~(\ref{1}) can be though as a first order Taylor expansion of the energy density around the transition pressure. Even though such a treatment lacks a rigorous theoretical basis it is widely employed as it is mimics the dynamics of the phase transition and it also allows an easy construction of EOSs predicting twin star configurations. In the present work the speed of sound is set equal to the speed of light in order to obtain EOSs consistent with the 2 $M_\odot$ constraint ~\cite{Montana-2019,Tsaloukidis-2023}.

A first order phase transition, between hadronic and quark matter, is not sufficient by itself for the appearance of a third family of compact objects. In particular, the appearance of twin stars requires the existence of an unstable region in the $M$-$R$ plane where the mass decreases with increasing central pressure. The condition that needs to be satisfied in order to obtain a third family was first studied by Seidov \cite{Seidov-1971} and it is formulated as it follows
\begin{equation}
    3 P_{\rm tr }+3\mathcal{E}_1-2\mathcal{E}_2<0,
    \label{eq2}
\end{equation}
where $\mathcal{E}_1\equiv \mathcal{E}(P_{{\rm tr}})$ and $\mathcal{E}_2\equiv \mathcal{E}(P_{{\rm tr}}) + \Delta \mathcal{E}$. Thus, by reorganising Eq.~(\ref{eq2}) one obtains the minimum energy density jump for the existence of twin star configurations, which is written as
 \begin{equation}
\Delta \mathcal{E}_{\rm cr}=\frac{1}{2}\mathcal{E}_{\rm tr}+\frac{3}{2}P_{\rm tr}.
\label{eq:seidov_limit}
\end{equation}
For EOSs that predict $\Delta \mathcal{E}\geq\Delta \mathcal{E}_{\rm cr}$ two distinct stable branches appear in the $M$-$R$ plane.

The resulting hybrid EOSs ought to be consistent with neutron star observations.~For example, if one assumes that the $\sim$ 1.4 $M_\odot$ compact stars involved in GW170817~\cite{Abbott-2017} or in PSR J0030+0451~\cite{Riley-2019} are purely hadronic, then the low density sector of the EOS has to satisfy tight constraints ($\Lambda_{1.4}=190^{+390}_{-120}$ and $R_{1.4}\leq 14$ km, where $\Lambda$ denotes the dimensioless tidal deformability)~\cite{Landry-2022}. On the other hand, if these compact objects are hybrid stars, the aforementioned constraints are lifted from the hadronic part of the EOS. In the present work we adopt the GRDF-DD2 (simply DD2 from now on for practical purposes)~\cite{Typel-2018} and the NL3~\cite{Lalazissis-1997} EOSs for the description of the low density phase. It is worth commenting that both of these EOSs have been previously employed in the study of twin stars~\cite{Tsaloukidis-2023,Sen-2022}. Finally, for the description of the outer crust (in the case of the NL3 model) the well-known EOS of Baym et al.~\cite{Baym-1971} is employed.

\section{R-mode instability formalism} \label{3}
Thermodynamics and the influence of  various dissipative processes define the time evolution  of the  {\it r}-modes according to the law  $ e^{i \omega t-t/\tau}$, where $\omega$ is the real part of the frequency, given by
\begin{equation}
\omega=-\frac{(l-1)(l+2)}{l+1}\Omega. \label{omega-1}
\end{equation}
In Eq.~(\ref{omega-1}), $\Omega$ is the angular velocity of the unperturbed star~\citep{Lidblom-98}
and $l$ defines the kind of mode. In  the present study we will consider the case $l=2$.
The imaginary part $1/\tau$ is related to the effects of
gravitational radiation and the various kinds of  viscosity (shear, bulk, etc.)
\cite{Lidblom-2000a,Owen-98,Lidblom-98}. We consider the case of  small-amplitude
limit where  a mode is a driven, damped harmonic oscillator and the
exponential damping time scale  is given by
\begin{eqnarray}
\frac{1}{\tau(\Omega,T)}&=&\frac{1}{\tau_{_{GR}}(\Omega)}+
\frac{1}{\tau_{_{EL}}(\Omega,T)}+\frac{1}{\tau_{_{BV}}(\Omega,T)}\nonumber \\
&+&\frac{1}{\tau_{_{SV}}(\Omega,T)},
\label{t-1}
\end{eqnarray}
where $\tau_{_{GR}}$, $\tau_{_{EL}}$, $\tau_{_{BV}}$ and $\tau_{_{SV}}$    are the gravitational radiation time scale, the damping time scale due to viscous dissipation at the boundary layer of the rigid crust and fluid core and  the  bulk and shear viscosity dissipation times scales respectively. It is notable that there is a battle between the gravitational radiation which tends to drive the r-mode unstable and the stabilization induced by  the various  dissipation mechanisms.
The  critical angular velocity $\Omega_{\rm c}$ (or critical spin frequency $f_{\rm c}=\Omega_{\rm c}/2\pi$), corresponds to the velocity  at which the two mechanisms (amplification and damping) are balanced and it is found through the equation $1/\tau(\Omega_c)=0$.

The contribution of gravitational radiation to the imaginary part of the frequency of the mode $1/\tau_{_{GR}}$ is given by the expression
\cite{Lidblom-2000a,Lidblom-98}
\begin{eqnarray}
\frac{1}{\tau_{_{GR}}}&=&-\frac{32\pi G
\Omega^{2l+2}}{c^{2l+3}}\frac{(l-1)^{2l}}{[(2l+1)!!]^2}\left(\frac{l+2}{l+1}\right)^{2l+2}\nonumber \\
&\times&
\int_0^{R}\rho(r) r^{2l+2} dr \quad \left({\rm s}^{-1}\right), \label{tgr-1}
\end{eqnarray}
where $\rho(r)$ is the mass density profile of a star.

The bulk viscosity $\xi_{_{BV}}$ is the dominant damping mechanism at high temperatures. It originated from the variations of pressure and density due to the pulsation modes and in nucleonic matter it is given by the formula~\cite{Lidblom-98}
\begin{eqnarray}
\xi_{_{BV}}^H&=&6.0\times 10^{-59}\left(\frac{l+1}{2}\right)^2 \left(\frac{{\rm Hz}}{\Omega}  \right)^2  \nonumber\\
&\times&\left(\frac{\rho}{{\rm gr\ cm^{-3}}}\right)^{2} \left(\frac{T}{{\rm K}}\right)^{6} \quad ({\rm  gr \ cm^{-1}\  s^{-1}}).
 \label{bulk-2}
\end{eqnarray}
For quark matter, the bulk viscosity is mainly determined by the weak process $d+s\leftrightarrow{}u+s$ \cite{Madsen-1992}. Following the discussion of Refs.~\cite{Madsen-1992,Jaikumar-2008} we will use an approximate expression which is appropriate for small  oscillations of the fluid and when $2\pi T\gg\delta\mu=\mu_s-\mu_d$. Specifically,
\begin{equation}
\xi_{_{BV}}^Q = \frac{\alpha T^2}{\omega^2+\beta T^4} \quad  { (\rm g \  cm^{-1} \ s^{-1}),}
\label{bulk-quark}
\end{equation}
where
\[ \alpha T^2=6.66 \times 10^{20} \left( \frac{\mu_d}{\rm MeV}\right)^3\left( \frac{m_s}{\rm MeV}\right)^4 T_9^2 \quad  { (\rm g\ cm^{-1}\ s^{-3}),}  \]
\[ \beta T^4=3.57 \times 10^{-8} \left( \frac{\mu_d}{\rm MeV}\right)^6\left(1+ \frac{m_s^2}{4 \mu_d^2}\right)^2 T_9^4 \quad  { (\rm s^{-2}),}  \]
where $T_9=T/(10^9 {\rm K})$, $\mu_d$ is the chemical potential of the down quark and $m_s$ is the mass of the strange quark. Since our model for quark matter does not provide information about the chemical potential profiles we will rely on the approximate expression $\mu_d = 235$ MeV $(\rho/\rho_0)^{1/3}$~\cite{Madsen-1992}, which has been employed in numerous r-mode studies~\cite{Madsen-1992,Madsen-2000,Anderson-2002,Pan-2006,Zheng-2006,Lindblom-1999a}. For the strange quark mass we assume that $m_s=$ 100 MeV. Finally, the bulk viscosity timescale is given by~\cite{Lidblom-98,Vidana-2012}
\begin{eqnarray}
\frac{1}{\tau_{_{BV}}}&=&\frac{4\pi}{690}\left(\frac{\Omega}{\Omega_0} \right)^4R^{2l-2}\left(\int_0^R \rho(r) r^{2l+2} dr \right)^{-1}\nonumber\\
&\times&
\int_0^R \xi_{_{BV}} \left(\frac{r}{R}  \right)^6\left[ 1+0.86\left(\frac{r}{R}  \right)^2
\right]r^2 dr, \label{bulk-1}
\end{eqnarray}
where $\Omega_0=\sqrt{\pi G\overline{\rho}}$ and
$\overline{\rho}=3M/4\pi R^3$ is the mean density of the star.

The shear viscosity is the dominant mechanism at low temperature. This mechanism  is  due to  the momentum transport  when  particle-particle scattering processes take place. In particular, the viscosity associated with the  neutron-neutron scattering and the electron-electron scattering are given respectively by \cite{Lidblom-2000a}
\begin{equation}
\eta_{nn}=347  \left(\frac{\rho}{{\rm gr\ cm^{-3}}}\right)^{9/4} \left(\frac{T}{{\rm K}}\right)^{-2}\quad({\rm g\ cm^{-1}\ s^{-1}}),
\label{eta-nn-1}
\end{equation}
\begin{equation}
\eta_{ee}=6.0\cdot 10^6 \left(\frac{\rho}{{\rm gr\ cm^{-3}}}\right)^2 \left(\frac{T}{{\rm K}}\right)^{-2}\quad({\rm g\ cm^{-1}\ s^{-1}}). \label{eta-ee-1}
\end{equation}
For quark matter the shear viscosity is dominated by quark-quark scattering in QCD. Following Ref.~\cite{Jaikumar-2008} we have
\begin{equation}
\eta_{q}= 5\times 10^{15} \left(\frac{0.1}{\alpha_s}  \right)^{3/2}\left(\frac{\rho}{\rho_0}  \right)^{14/9}T_9^{-5/3}
\quad  { (\rm g \  cm^{-1} \ s^{-1})},
\label{shera-quark}
\end{equation}
where $a_s$ is the coupling constant for the strong interaction. In the present work we will use a typical value of $a_s=$ 0.1.
The dissipation time scale  due to the shear viscosity is given by \citep{Lidblom-98}
\begin{eqnarray}
\frac{1}{\tau_{_{SV}}}&=&(l-1)(2l+1) \left(\int_0^R \rho(r) r^{2l+2} dr \right)^{-1}\nonumber\\
&\times&
\int_0^R \eta_{_{SV}} r^{2l} dr, \quad ({\rm s^{-1}}).
\label{shear-1}
\end{eqnarray}

In the special case where the dissipation effect due to the crust has been included, the corresponding time scale is given by~\cite{Lidblom-2000a}
\begin{eqnarray}
\tau_{_{EL}}&=&\frac{1}{2\Omega}\frac{2^{l+3/2}(l+1)!}{l(2l+1)!! {\cal
C}_l}\sqrt{\frac{2\Omega R_c^2\rho_{cr}}{\eta_{cr}}}\nonumber\\
&\times&
\int_0^{R_c}\frac{\rho(r)}{\rho_{cr}}\left(\frac{r}{R_c}\right)^{2l+2}
\frac{dr}{R_c} \quad ({\rm s}).
 \label{tv-1-mew}
\end{eqnarray}
In Eq.~(\ref{tv-1-mew}), $R_c$ is the core's radius, while $\eta_{cr}$ and $\rho_{cr}$ are the viscosity and density of the fluid at the outer edge of the core, respectively. The factor ${\cal C}_l$, for $l=2$,  takes the value ${\cal C}_2=0.080411$. The expression~(\ref{tv-1-mew}) refers to  the case where the crust is rigid and consequently static in the rotating frame. However, in a more realistic case, the motion of the crust (due to the mechanical coupling with the core) induces an increase of the  timescale $\tau_{_{EL}}$  by a factor of $1/{\cal S}^2$, where ${\cal S}$ is the slippage factor defined as ${\cal S}=\Delta v/v$. In particular,  $v$ denotes the velocity of the core and  $\Delta v$ is the difference between the velocities in the inner edge of the crust and the outer edge of the core~\cite{Levin-2001}.

It has been shown that the critical frequency $\Omega_{\rm c}$ is quite sensitive to the radius of a star~\cite{Papazoglou-2016}. More precisely, it has been found that for relatively low and high values of  temperature, $\Omega_{\rm c}$ scales with the radius as   $\Omega_{\rm c}\sim 1/R^{3/2}$ and $\Omega_{\rm c}\sim 1/R^{3/4}$, respectively. The latter leads to the conclusion that the r-mode instability windows of twin stars are going to be different. Furthermore, if one takes into account that the damping mechanisms in quark matter are, in principal, stronger than those in hadronic matter, then the instability window of a hybrid star is expected to be shifted to larger $\Omega_{\rm cr}$ compared to the one of its hadronic twin. Thus, there are two mechanisms which act additively and may drastically affect the instability windows of the two different branches. The above findings are essentially a strong motivation for investigating the possible identification of twin stars due to the implications of their different instability windows.

Another  crucial issue is the limitation of the instability window, at high frequencies, from the corresponding Kepler angular velocity  $\Omega_{\rm K}$ (the maximum rotation frequency of the star). To a very good approximation, the Kepler velocity is given by $\Omega_{{\rm K}}=\frac{2}{3}\Omega_0$. It is interesting that in the case of twin stars  the Kepler frequency of the hybrid branch may be even $20\%$ higher compared to the hadronic branch, due to the different radius values.  This apparent differentiation at the upper limit of the instability window can by itself be a criterion for separating the two branches. Connecting the analysis presented above with the fact that newly born compact stars are expected to rotate close to their mass shedding limit, we conclude that the spin-down evolution paths of twin stars are going to exhibit distinct deviations.

\section{Spin-down and cooling }\label{4}
We are now going to present a simplified model to describe the spin-down (due to unstable r-modes) of a hadronic or hybrid star simultaneously with its cooling. During the phase that the angular momentum is radiated away to infinity by gravitational waves, the angular velocity of a star evolves as follows~\cite{Owen-98}
\begin{equation}
\frac{d\Omega}{dt}=\frac{2\Omega}{\tau_{GR}}\frac{\alpha^2 Q}{1-\alpha^2Q},
\label{dOmegadt-1}
\end{equation}
where $\alpha$ is the dimensionless r-mode amplitude  parameter. This parameter strongly affects the r-mode evolution and usually takes values
in the large  interval  $\alpha=1-10^{-8}$. Moreover, $\alpha$ in general depends both on the viscosity (and consequently on the temperature $T$ and cooling process) and on time. However, following Ref.~\cite{Lidblom-98} we consider that $d\alpha/dt=0$.  In addition, the quantity $Q$  related with the bulk properties of a star and is defined as $Q=3\tilde{J}/2\tilde{I}$ where
\begin{equation}
\tilde{J}=\frac{1}{MR^4}\int_{0}^R \rho(r)r^6 dr, \qquad \tilde{I}=\frac{8\pi}{2MR^2}\int_{0}^R \rho(r)r^4 dr.
\label{I-J-1}
\end{equation}
Under the afforementioned assumptions, on can solve Eq.~(\ref{dOmegadt-1}) analytically and obtain~\cite{Moustakidis-2015,Routray-2021}
\begin{equation}
\Omega(t)=\left(\frac{1}{\Omega^{-6}_{in} -6{\cal C}t}\right)^{1/6},
\label{dOmegadt-2}
\end{equation}
where
\begin{equation}
{\cal C}=\frac{2\alpha^2Q}{\tilde{\tau}_{GR} (1-\alpha^2Q)}\frac{1}{\Omega_0^6}, \quad \tilde{\tau}_{GR}=\left(\frac{\Omega}{\Omega_0}\right)^{6}\tau_{GR}.
\label{Cons-1}
\end{equation}
$\Omega_{in}$ is a free parameter which corresponds to the initial angular velocity and $\tilde{\tau}_{GR}$ is the fiducial gravitational-radiation time scale.

In order to combine the concurrent processes  of the  spin-down and cooling of twin stars,  we use the standard description for the cooling of hot and young neutron stars proposed by Owen et al.~\cite{Owen-98} (where is considered that the cooling primarily due to the emission of neutrinos via the modified URCA process). In this case, the temperature drops according to the law
\begin{equation}
T(t)=\left(\frac{t}{t_c}+\left(\frac{10^9 \ {\rm K}}{T_i}  \right)^6   \right)^{-1/6} 10^9 \ {\rm K}
 \label{T-t}
\end{equation}
where  $T_i$ is the initial temperature of the star (a typical value is $T_i\simeq 10^{11} \ {\rm K}$) and $t_c$ is the cooling rate parameter ($t_c\simeq 1\ {\rm year}$~\cite{Owen-98}).

It is worth mentioning that this is a simplified model, especially concerning the cooling process. In particular, one may expect that the two twins may cool in a different way, considering the different cooling mechanisms in hadronic and quark matter. Obviously, a more elaborated study is necessary if one is interested in an accurate quantitative description of the cooling process. However, we need to comment that according to the findings of Lyra et al. \cite{Lyra-2023}, the thermal evolution of twin stars is only distinct when there is a large difference in their compactness. More precisely, in the case where there is a 10~$\%$ compactness difference (which is the case for the models employed in the present study), the thermal evolution of the two twins is nearly identical \cite{Lyra-2023}. Therefore, we conclude that the present model will provide, at least, a good qualitative picture for the evolution of twin stars.

 From the above analysis it is obvious that young, hot and rapidly rotating  twin stars follow different spin-down paths on the $f-T$ plane. This is mainly due to the following three reasons: a) the different instability  windows, b) the different Kepler frequencies (different spin frequencies at birth) and c) the different spin evolution with time.

\section{Results and discussion} \label{5}
\subsection{Mass-Radius diagrams}
 \begin{figure*}[t]
  \includegraphics[width=\textwidth,scale=1]{MR.pdf}
  \caption{Mass-Radius diagrams for the DD2 (left panel) and NL3 (right panel) EOSs. The black solid curves indicate the original EOSs. The solid (dashed) horizontal line is set to 1.4 (1.2)~$M_\odot$. The shaded areas correspond to the constraints from the analysis of the GW170817 event \cite{Abbott-2017,Abbott-2018}. Each hybrid EOS is identified from the baryon density where the phase transition occurs and by the energy density jump. The energy density gap is given in units of MeV fm$^{-3}$ which are omitted in the legend for simplicity.}
  \label{Fig1}
\end{figure*}

In order to study the differences that manifest in the r-mode instability windows and spin evolution of twin stars we constructed a set of hybrid EOSs, using the analysis presented in Section~\ref{2}. In particular, the low density phase is described by the DD2 and NL3 EOSs, where for the quark matter a phenomenological constant speed of sound model is employed. The values of the energy jump are selected in order to obtain EOSs that are consistent with the constraints from astrophysical observations. For both hadronic models, the resulting EOSs predict twin stars with mass of 1.2 or 1.4 $M_\odot$.

Fig.~\ref{Fig1} depicts the mass-radius dependence for the EOSs employed in this study.~In the left panel the hadronic phase is described using the DD2 EOS, while for the results of the right panel the NL3 model was employed. The solid black curves stand for the case where no phase transition occurs (i.e. the $M$-$R$ diagrams for the purely hadronic EOSs).~In addition, the shaded areas correspond to constraints based on the analysis of the GW170817 event~\cite{Abbott-2017,Abbott-2018}. Finally, the horizontal lines are drawn to indicate the twin configurations with 1.2 and 1.4 $M_\odot$.

As it is evident from Fig.~\ref{Fig1}, increasing the energy density jump results into a softening of the EOS. Thus, the largest values for $\Delta\mathcal{E}$ are selected so that EOSs remain consistent with the 2 $M_\odot$ constraint. Furthermore, we need to highlight that as $\Delta\mathcal{E}$ increases the radius difference of the two twins becomes larger. The latter is expected to play a critical role concerning the deviation of r-mode instability windows \cite{Papazoglou-2016}. It is important to note that our analysis does not include the limiting case where $\Delta\mathcal{E}=\Delta\mathcal{E}_{cr}$, as in such a scenario the separation of the two twins is almost negligible. In particular, if the phase transition occurs in relatively low baryon density a third family may not even appear \cite{Tsaloukidis-2023}.

\begin{figure}
  \centering
  \includegraphics[width=9.5 cm,scale=1]{eos_effect.pdf}
  \caption{Critical spin frequency $f_c$ as a function of the redshifted temperature $T^\infty$ (r-mode instability windows) for 1.4 $M_\odot$ twin stars for the DD2 (blue) and NL3 (red) EOSs. The dashed lines and solid lines correspond to the hadronic and hybrid twins, respectively. The horizontal lines denoted the Kepler frequency for its star. The value for the energy density gap is $\Delta\mathcal{E}_{cr}+$100 MeV fm$^{-3}$ for both EOSs.}
  \label{Fig2}
\end{figure}

\subsection{Qualitative analysis}
Fig.~\ref{Fig2} presents the r-mode instability windows of 1.4~$M_\odot$ twin stars for the case where $\Delta\mathcal{E}=\Delta\mathcal{E}_{cr}+$100 MeV fm$^{-3}$. The damping mechanism due to the presence of a solid crust is not included. The results for the hadronic and hybrid stars are indicated using dashed and solid curves, respectively. Additionally, the horizontal lines stand for the corresponding Kepler frequencies. Note that in the x-axis of the plot, one would not find the temperature $T$, appearing in the formalism of Section \ref{3}, but the so-called redshifted temperature which is given by $T^\infty=T\sqrt{1-2C},$
where $C=GM/Rc^2$ is the compactness of a star.

Firstly, we need to underline the sensitivity of the instability window to the employed EOS. Specifically, by comparing the instability windows of the purely hadronic configurations one finds that the predicted critical frequency is lower (in the low temperature region) when the NL3 EOS is employed. This results from the fact that the radius of a 1.4 $M_\odot$ compact star is smaller when the DD2 model is used~\cite{Papazoglou-2016}. Incidentally, the radius of the hybrid star constructed using the NL3 EOS coincides with the radius of the hadronic configuration using the DD2 EOS. The latter results into an overlap of their instability windows in the low temperature regime. However, the existence of a quark core, in the hybrid star, leads to significant differences in the critical frequencies for $T^\infty\geq$ 10$^8$ K, where the bulk viscosity plays a crucial role.

For a qualitative comparison of the r-mode instability windows of twin pairs one can divide Fig. \ref{Fig2} in three representative regions. In particular, for $T^\infty\leq$ 10$^8$ K (where the shear viscosity is the dominant dissipation mechanism), the radius difference plays a crucial role for the apparent critical frequency deviations. For 10$^8$~K~$\leq T^\infty\leq$~10$^{10}$~K, the bulk viscosity is the major damping mechanism and the trend of the $f_c(T^\infty)$ curve is altered for the hybrid twin. More precisely, the critical frequency increases and then decreases with temperature leading to a local maxima. This topological difference derives from the fact that the bulk viscosity of quark matter is not a monotononic function of temperature. Finally, for $T^\infty\geq$ 10$^{10}$ K the bulk viscosity dominates and this leads to appreciable differences in $f_c$ for the two twins. From an observational perspective, the differences that appear in the low temperature regime will lead to different limits on the spin up of accreting pulsars in LMXBs. On the other hand, the deviations in the large temperature region may affect the evolution of a rapidly rotating proto-neutron star. Apart from the differences appearing in the limiting spin frequencies of twin stars, due to unstable r-modes, we need to comment that the there is a $\sim$ 17 $\%$ difference in their Kepler frequencies. Subsequently, a young hybrid star can rotate much faster than its hadronic twin.

\subsection{Energy density jump and crust effects}

\begin{figure}
  \centering
  \includegraphics[width=8.5 cm,scale=1]{relations.pdf}
  \caption{Panel a: Radius difference between 1.4 $M_\odot$ twin stars as a function of the energy density gap, Panel b: The quark core fraction of a 1.4 $M_\odot$ hybrid star as a function of the energy density jump. }
  \label{Fig3}
\end{figure}
\begin{table}
\caption{The difference in the critical frequencies for 1.4 M$_\odot$ twin stars using the NL3 model for different values of temperature and energy density jump. The damping due to a viscous boundary layer (rigid crust) is not considered for the results presented in this table.}
\begin{center}
\setlength\tabcolsep{0pt}
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}} cccc }
    \hline
    \hline
     $\Delta\mathcal{E}$ (MeV fm$^{-3}$) & $\Delta R$ (km) & $T^\infty$ (10$^8$ K) & $\Delta f_c$ (Hz)\\
    \hline
    \multirow{3}{*}{ $\Delta\mathcal{E}_{cr}+$100} & & 1 & 81.93 \\ & 1.75 & 5 & 192.47\\ & & 10 & 82.62\\
    \hline
    \multirow{3}{*}{$\Delta\mathcal{E}_{cr}+$150} &  & 1 & 136.77 \\ & 2.40 & 5 & 306.32 \\ & & 10 & 149.006\\
    \hline
    \multirow{3}{*}{$\Delta\mathcal{E}_{cr}+$200} & & 1 & 187.07 \\ & 2.99 & 5 & 402.53\\ & & 10 & 207.14\\
    \hline
    \hline
\end{tabular*}
\end{center}
\label{t1}
\end{table}

\begin{table}
\caption{The difference in the critical frequencies for 1.4~$M_\odot$ twin stars for different values of temperature and crust elasticity. The results were obtained using the NL3 model with $\Delta\mathcal{E}=\Delta\mathcal{E}_{cr}+$200 MeV fm$^{-3}$. The results for this EOS, in the case where the crust damping mechanism is not considered, can be found in Table \ref{t1}.}
\begin{center}
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}} ccc }
    \hline
    \hline
     $S$  & $T^\infty$ (10$^8$ K) & $\Delta f_c$ (Hz)\\
    \hline
    \multirow{3}{*}{0.2}   & 1 & 102.279 \\  & 5 & 199.014\\  & 10 & 103.426\\
    \hline
    \multirow{3}{*}{1}  & 1 & 165.56 \\  & 5 & 169.475\\  & 10 & 129.472\\
    \hline
    \hline
\end{tabular*}
\end{center}
\label{t2}
\end{table}
\begin{figure*}
  \includegraphics[width=\textwidth,scale=0.8]{crust_effects.pdf}
  \caption{The effect of energy density gap in the deviation of the r-mode instability windows of 1.4~$M_\odot$ twin stars for increasing crust elasticity values. Panel a) DD2 EOS and $\Delta\mathcal{E} = \Delta\mathcal{E}_{cr}+$75 MeV fm$^{-3}$, b) DD2 EOS and $\Delta\mathcal{E} = \Delta\mathcal{E}_{cr}+$100 MeV fm$^{-3}$, c) NL3 EOS and $\Delta\mathcal{E} = \Delta\mathcal{E}_{cr}+$100 MeV fm$^{-3}$, d) NL3 EOS and $\Delta\mathcal{E} = \Delta\mathcal{E}_{cr}+$100 MeV fm$^{-3}$. Dashed (solid) lines indicate the hadronic (hybrid) twin. The dotted points correspond to observational data taken from Ref.~\cite{Gusakov-2014}.}
  \label{Fig5}
\end{figure*}
At this point we wish to systematically study the influence of certain parameters on the instability window deviations of twin stars. In particular, we are going to vary the value of the energy density gap and examine its effects. Furthermore, up to this point the only dissipative mechanisms considered in our calculations, were the bulk and shear viscocities. Now, we are also going to include the damping mechanism due to the presence of a viscous boundary layer. It is interesting that, as the aforementioned mechanism is strong and common for both twins, the critical frequency deviations due to different viscocities are expected to be less pronounced.

Firstly, we are going to investigate the importance of the energy density jump. As we mentioned, $\Delta\mathcal{E}$ is the regulator of the radius difference between twin configurations. Fig.~\ref{Fig3}a depicts the dependence of $\Delta R$ on $\Delta \mathcal{E}$ for 1.4 $M_\odot$ twin stars. Surprisingly, we find that the aforementioned quantities are connected through a linear formula. Even though the exact $\Delta R$-$\Delta \mathcal{E}$ relation is sensitive to the low density model, the slopes of the resulting fitted lines appear to be very similar. It is worth pointing out that, from the analysis presented in the previous section, an increment of $\Delta \mathcal{E}$ will result into larger deviations in the instability windows due to an increase of $\Delta R$. However, as it is evident from Fig.~\ref{Fig3}b, a larger value of $\Delta \mathcal{E}$ also results into a hybrid twin with a larger quark core fraction $x_q=R_q/R$ (where $R_q$ is the quark core radius). Hence, the damping due to quark matter's bulk viscosity is going to be even more effective. It is noteworthy that $\Delta \mathcal{E}$ and $x_q$ are also linearly dependent and that the slopes of the lines are, once again, not strongly sensitive to employed hadronic model. The relations presented in Fig.~\ref{Fig3} can be added to the other correlations found in the detailed analysis of Ref.~\cite{Sen-2022}. Finally, we need to underline that through the relations found above, the knowledge of the radius difference of twin stars, may provide important information concerning the phase transition and the interior of hybrid stars.

Fig.~\ref{Fig5} depicts the dependence of critical frequency on temperature for 1.4 $M_\odot$ twin stars constructed using different $\Delta \mathcal{E}$ values. In addition, the crust elasticity $\mathcal{S}$ is varied from 0 to 1 in order to investigate the effects of a viscous boundary layer. The circular points stand for observational data inferred from LMXBs and millisecond pulsars, while the uncertainties visualised though the error bars derive from the process of converting surface temperature measurements to estimations for the core temperature \cite{Gusakov-2014}. Tables \ref{t1}, \ref{t2} contain numerical data of $\Delta f_c$ for different energy density jump, crust elasticity and temperature values.

It is worth pointing out that, the instability window differences are more pronounced in the case where the NL3 model is employed. This results from the fact that as the NL3 model is stiffer, it allows the construction of EOSs that satisfy observational constraints even for large $\Delta \mathcal{E}$ values. Furthermore, Fig.~\ref{Fig5} illustrates the strong impact of $\Delta \mathcal{E}$ on the resulting r-mode instability window of the hybrid twin. In particular, in the case of $\Delta \mathcal{E}_{cr}$ + 200, the spin frequency difference for the two twins may reach values of $\sim$ 400 Hz (see Table \ref{t2}).

 \begin{figure*}
  \includegraphics[width=\textwidth,scale=1]{mass_effect.pdf}
  \caption{R-mode instability windows of compact stars in the mass range 1.2$-$1.9 $M_\odot$. The EOSs used are: a) DD2, $n_{tr}=$ 0.32 fm$^{-3}$, $\Delta \mathcal{E}_{cr}+$ 150, b) DD2, $n_{tr}=$ 0.35 fm$^{-3}$, $\Delta \mathcal{E}_{cr}+$ 100, c) NL3, $n_{tr}=$ 0.25 fm$^{-3}$, $\Delta \mathcal{E}_{cr}+$ 150, d) NL3, $n_{tr}=$ 0.27 fm$^{-3}$, $\Delta \mathcal{E}_{cr}+$ 150. The dotted points correspond to observational data and the are taken from Ref.~\cite{Gusakov-2014}. The (h) appearing in the legend stands for the most massive purely hadronic configuration. The value for the crust elasticity is considered to be $\mathcal{S}=$ 0.1.}
  \label{Fig6}
\end{figure*}

The most important effect when the damping due to a solid crust is included, is that the peak appearing in the instability window of the hybrid star (in a temperature region around $\sim$ 5$\times$10$^8$ K) drops down. However, depending on the selected $\Delta \mathcal{E}$ value, large $f_c$ differences for the two twins may remain. In accordance to the results present by Lyra et al. \cite{Lyra-2023} we conclude that the role of the compactness is not only critical concerning the thermal evolution of twin pairs, but it also signifficantly affects the r-mode instability window of the hybrid configuration.

\subsection{Comparison with observational data}
 \begin{figure*}
  \includegraphics[width=\textwidth,scale=1]{spin_down.pdf}
  \caption{Panel: (a) Spin frequency as a function of time for 1.4 $M_\odot$ twin stars constructed with the DD2 EOS, (b) Spin-down rate as a function of time for 1.4 $M_\odot$ twin stars constructed with the DD2 EOS, (c) Spin frequency as a function of time for 1.4 $M_\odot$ twin stars constructed with the DD2 EOS, (d) Spin-down rate as a function of time for 1.4 $M_\odot$ twin stars constructed with the NL3 EOS. In all panels two different values for $\Delta\mathcal{E}$ were used (see legends).
}
  \label{spin-down-1}
\end{figure*}

 \begin{figure*}
  \centering
  \includegraphics[width=12 cm,scale=1]{kepler_evolution.pdf}
  \caption{The spin-down evolution of 1.4 $M_\odot$ twin stars (NL3 EOS, $n_{tr}=$ 0.27 fm$^{-3}$ and $\Delta\mathcal{E}=\Delta\mathcal{E}_{cr}+$ 100 MeV fm$^{-3}$) in the frequency-temperature plane for different values of the saturation amplitude.  The initial frequencies for the twins are their corresponding Kepler frequencies. The blue (red) solid lines indicate the evolution for the hybrid (hadronic) twin. The blue and red dashed lines denote the r-mode instability window of the hybrid and hadronic star, respectively. The dotted points stand for observational data taken from Ref.~\cite{Gusakov-2014}.}
  \label{spin-down-2}
\end{figure*}

As previously mentioned, it is rather difficult to explain the observational data in the context of a purely hadronic star. In particular, the not realistic assumption of a perfectly rigid crust is essential \cite{Haskell-2012}. For that matter, several studies have investigated the r-mode instability window of compact stars containing exotic forms of matter. In a recent work, Ofengeim et al. \cite{Ofengeim-2019a,Ofengeim-2019b} examined if the existence of hyperons in the core of compact stars can lead to results compatible with current LMXBs data. What they found is that for neutron stars with $M\leq 1.9$~$M_\odot$, the bulk viscosity of hyperonic matter leads to r-mode stabilization in the $f-T^\infty$ regime where observed neutron stars appear \cite{Ofengeim-2019a,Ofengeim-2019b}.

At this point we wish to examine if the hybrid EOSs constructed in this study are in accordance to the observed spin frequencies and temperatures in LMXBs. As it is evident from Fig.~\ref{Fig5}, the r-mode instability window of the hybrid twin is always narrower. In addition, in a minimal scenario where the effects of the crust are not included the explanation of the observational data is not possible. However, depending on the energy density jump, a moderate crust elasticity value would suffice for the construction of instability windows compatible with observations. Specifically, for the NL3 model ($\Delta \mathcal{E}_{cr}+$ 200) and a relatively small crust elasticity of 0.2, most of the observed stars lay in the stable region of the $f-T^\infty$ plane. Another critical point is that, in all cases there are stars (from the dataset) that lay in the region between the $f_c(T^\infty)$ curves for the two twins. Hence, while such stars can be considered stable with respect to r-modes in the framework of the hybrid twin, they would be unstable if they were purely hadronic. The latter comment is of most importance, as the detection of GW emission, from stars laying in a $f-T^\infty$ region where r-modes are considered to be stable, would be a strong indication of hadron-quark phase transition.

In Fig.~\ref{Fig6} we present the r-mode instability windows for compact stars in the mass range 1.2$-$1.9 $M_\odot$ and a relatively low crust elasticity value $\mathcal{S}=$ 0.1. In panels a and c the twin configurations have a mass of 1.2 $M_\odot$, while in panels b and d their mass is 1.4 $M_\odot$. In the first case we find that, the bulk viscosity of quark matter is sufficient to stabilize r-modes for moderately massive compacts stars ($M \leq 1.6$~$ M_\odot$) in the whole $f-T^\infty$ range occupied by the observed stars in LMXBs. In latter case, where the phase transition occurs at higher baryon density, more massive compact star configurations (1.8 or 1.9 $M_\odot$ depending on the hadronic EOS) are essential for the explanation of current LMXBs data.

Another observation that can be made from Fig.~\ref{Fig6}, is that right after the phase transition occurs a narrowing of the instability window is evident. Then, as the mass further increases the instability window becomes wider for low temperature values ($T^\infty \leq$ 10$^8$ K). The fact that higher mass configurations have wider instability windows is a known result from previous studies \cite{Mukhopadhyay-2018}. It is interesting that while a higher mass is necessary for the stabilization of r-modes in observed stars with $T^\infty \geq$ 10$^8$ K, it fails to provide an explanation for the stars appearing in a lower temperature regime. However, the low temperature region can be cover by hybrid star configurations of lower mass. In principle, if a star slightly surpasses a critical mass, after which a phase transition occurs, then its instability window will be also slightly different from the one of the most massive purely hadronic configuration. In contrast, if the structure of the phase transition predicts the existence of a third family, then stars with mass equal or slightly larger than the aforementioned critical mass are going to exhibit considerable deviations in their r-mode instability windows. The non trivial behavior of stars having narrower instability windows compared to those of lower mass stars (for low $T^\infty$), is characteristic of an EOS predicting twin configurations.

\subsection{Spin down and thermal evolution}
In Fig.~\ref{spin-down-1} we display  the time evolution for the
frequency and the corresponding spin-down rate of 1.4~$M_\odot$ twin stars. The upper and lower panels contain results for the DD2 and NL3 EOSs, respectively. For comparison reasons, we consider the same initial frequency of 650~Hz for both twins. Furthermore, in accordance to previous studies \cite{Moustakidis-2015,Mukhopadhyay-2018}, the selected value for the r-mode saturation amplitude is $a=$ 2$\times$10$^{-7}$. The spin-down rate is slower in hybrid stars,  right after their birth. Specifically, the higher the energy density gap the lower the rate. However, after a certain amount of time the spin down rates of twin stars converge to the same value. The latter is reflected on the distinct time evolution of the frequency for the two cases. In particular, a hybrid star retains its initial rotation frequency for a longer period of time compared to its hadronic twin.

Usually, it is more convenient to study the spin down evolution of a compact star on the $f-T^\infty$ plane. The latter demands the simultaneous knowledge of the spin and thermal evolution for a star. By employing the toy model for the fall of temperature, presented in Section~\ref{4}, we intend to obtain the different evolution paths of twin stars on the $f-T^\infty$ plane. In addition, instead of considering the same initial frequency for the two twins we set as initial condition the corresponding Kepler frequencies. The results presented in Fig.~\ref{spin-down-2} were constructed using the NL3 EOS with $n_{tr}=$ 0.27 fm$^{-3}$ (hence 1.4~$M_\odot$ twins) and $\Delta\mathcal{E}=\Delta\mathcal{E}_{cr}+$ 100 MeV fm$^{-3}$. The latter EOS predicts twins with a $\sim$ 13 $\%$ difference in compactness and therefore a similar thermal evolution is not an unreasonable assumption \cite{Lyra-2023}. For the crust elasticity a low value of $\mathcal{S}=$ 0.1 was chosen. Moreover, we consider three different values of the amplitude $\alpha$ since the results are very sensitive to it. From   Fig.~\ref{spin-down-2} it is obvious that there are three main reasons which differentiate the time evolution of the two branches. The first one is the different Kepler velocities. The second one is connected to the spin down rates of the two twins, even though this effect is less pronounced. The third one is the deviation of the instability windows. In particular, the unstable region is more extended in the case of the hadronic branch. The latter is of most importance, as the r-mode instability window essentially sets the resulting frequency of a star as it comes out of the unstable region. Of course we need to stress out that, the paths presented in  Fig.~\ref{spin-down-2} can  be  improved if one considers a  more realistic cooling process for the two branches. However,  the general picture  will not change noticeably and the main conclusions of the present study are not expected to be significantly altered.

\section{Conclusion} \label{6}

The present work was dedicated to the study of twin stars and their r-mode instability windows. In particular, we have conducted a detailed investigation of the parameters that affect the deviation between the instability windows of twin stars. This is of most importance as two stars with identical mass may have different rotational limits. More precisely, two stars in the same region of the frequency-temperature diagram may behave differently with respect to r-modes. Subsequently, the future detection of (r-mode) GW emission, from stars that are considered to be stable with respect to r-modes (due to existing observations), would be a clear sign for the existence of a third family and hence of hadron-quark phase transition.

Firstly, we studied the influence of the energy density jump $\Delta\mathcal{E}$ on the deviation between the instability windows of twin stars. We found that $\Delta\mathcal{E}$ regulates the radius difference between twin configurations. In addition, hybrid twins predicted from EOSs with higher $\Delta\mathcal{E}$ exhibit larger quark core fractions. Thus, the differences in the critical spin frequencies of twins become more pronounced as the energy density jump increases. Secondly, we took into consideration the strong and common (for both twins) dissipation mechanism due to the presence of a viscous boundary layer. What we found is that, the characteristic peak appearing in the r-mode instability windows of hybrid stars (around $T^\infty$ $\sim$ 5 $\times$ 10$^8$ K) flattens as the crust elasticity increases. However, depending on the selected value of $\Delta\mathcal{E}$, considerable differences in the limiting frequency of the two twins may remain.

Furthermore, we examined if the EOSs constructed in this study (i.e. EOSs predicting a third family of compact objects) are a viable option for the explanation of current LMXBs data. We found that depending on the phase transition onset (transition density) and also the masses of stars in LMXBs, our EOSs may be compatible with the existing observational data. In particular, for EOSs that predict twin stars with 1.2 $M_\odot$, the bulk viscosity of quark matter
is adequate to stabilize r-modes for moderately massive stars ($M \leq 1.6~M_\odot$) in the whole $f-T^\infty$
region occupied by the observed stars in LMXBs. As the critical compact star mass for the phase transition to occur increases, more massive configurations are needed for the stabilization of r-modes.

Finally, we studied the differences that manifest in the spin down evolution of twin pairs. We found that the hybrid star retains its initial spin frequency for a larger period of time and this is because its spin down rate is lower compared to its hadronic twin. Furthermore, we noticed that larger $\Delta\mathcal{E}$ values result into lower spin down rates for hybrid stars. In addition, by employing a simplified thermal evolution model we evaluated the evolution paths of twins stars in the $f-T^\infty$ plane. The resulting path differences derive from: a) the fact that the Kepler frequencies (initial conditions) of twin stars are different, b) the different spin evolution, which is dependent on the bulk properties of a star, c) the different instability windows of twin stars which essentially control when and with what frequency a star is going to pass in the r-mode stable region.

 There are some other issues that a more elaborated study should take into account such as additional damping mechanisms or a more rigorous treatment of the thermal evolution. In addition, it would be interesting to explore the effects of a mixed phase (EOSs constructed with the Gibbs method). Even though such a study (already in progress) would be more complete from a quantitative point of view, we do not expect that our main conclusions will be significantly altered.
 Finally, we need to highlight that even though there are a few studies focusing on the r-mode instability and hybrid stars this is the first work dealing with the possible existence of two stars with identical mass and different r-mode instability windows. The future detection of GW  associated with unstable r-modes may finally allow us to distinguish twin stars.

\section*{Acknowledgements}
The authors would like to thank Prof. K. Kokkotas for their useful insight and comments


\section{Introduction}

Estimating how a human would move in the near future is an essential task for various applications such as surveillance \cite{surv_1, surv_2}, autonomous driving \cite{drive_1, drive_2}, and human-robot/computer-interaction \cite{hri_1}.
Many approaches have been proposed to solve this problem, often based on the motion capture datasets such as Human3.6M \cite{h36m_pami} or SMPL \cite{loper2015smpl}-based datasets such as AMASS \cite{mahmood2019amass}.
In this paper, we concern with a task whose goal is to predict a sequence of 3D pose skeletons in Human3.6M and HumanEva-I \cite{sigal2006humaneva} datasets, when a previously observed 3D pose sequence is given as an input.

Existing works on 3D skeleton motion prediction can be categorized as follows.
One line of research focuses on models for deterministic motion prediction \cite{deter_1, deter_2, s-rnn, DCT-GCN, ST-TR, 2CH-TR}.
These works aim at predicting a single motion that is most likely to be observed in the future.
Therefore, their performance is usually evaluated based on an $L2$-distance between a prediction and a ground truth.
Another line of research focuses on generative models for stochastic motion prediction \cite{hp-gan, mt-vae, Dlow, gsps}.
Their performance is evaluated based on the metrics for likelihood and diversity.
After generating a fixed number of prediction samples from a single observation,
the likelihood is measured based on the minimum distance between the prediction samples and ground truth, and the diversity is measured based on the average distance between all pairs of prediction samples.

\begin{figure}
    \centering
    \includegraphics[width=0.45\textwidth]{fig/fig1.png}
    \caption{Example results when diffusion probabilistic models are used for 3D human motion prediction tasks, when observed motion is `walking'. After a single training procedure, diffusion models can be effectively used for both deterministic (Deter.) and stochastic (Sto.) motion prediction tasks.}
    \label{fig:first}
    \vspace{-3mm}
\end{figure}

However, we cannot judge which approach is always better than the other, since the efficiency would depend on the target application values.
For instance, when one needs only the most precise sample with low latency, deterministic approaches would be better.
If we say that both approaches are necessary, our next question would be whether we can propose an efficient model for both types of prediction.
To answer the question, we study the possibility of using diffusion probabilistic models \cite{ddpm, csdi} for both deterministic and stochastic 3D motion prediction tasks.

If we propose a diffusion probabilistic model \cite{ddpm, csdi} as a solution, one might ask us whether this is because we are fascinated by its performance in image generation \cite{imagen, palette}.
Frankly speaking, yes, we initiated this study out of our curiosity -- can we use diffusion probabilistic models for 3D motion prediction?
Unfortunately, our experimental results show that the diffusion model cannot perfectly replace existing state-of-the-arts for both deterministic and stochastic motion prediction tasks.
However, we found a glimpse of hope in diffusion models, due to their effectiveness in both prediction types after a single training procedure, and their ability to properly balance the trade-off between diversity and likelihood.

Figure~\ref{fig:first} shows the example results when the diffusion models are used for both deterministic and generative motion prediction tasks.
Although a diffusion model is essentially a generative model, we found that the deterministic sample with a fair performance can be obtained from the diffusion model when all randomness is excluded from its denoising process.
In addition, we found out that the diffusion models can fix the flaws of several generative methods \cite{Dlow}, which highlight the diversity of generated samples. Existing works as \cite{Dlow} claim that the likelihood of predicted samples is high when the minimum distance between samples and ground truth is low. Because of this, \cite{Dlow} can often generate the motions that are out-of-context as \cite{contexuallyplausible} pointed out. Compared to this, our diffusion models can generate prediction samples that are more likely to occur, so the generated motion does not diverge too much to be called out-of-context.

The remaining paper is constructed as follows.
After representing our literature survey in Section~\ref{sec:rel},
Section~\ref{sec:met} will explain how general diffusion models work as well as how we design ours to solve 3D motion prediction tasks.
Section~\ref{sec:exp} will show both qualitative and quantitative experiment results, and a related discussion will be also presented. Finally, this paper will end in Section~\ref{sec:con} by mentioning limitations and future works.

\section{Related Work} \label{sec:rel}

\subsection{3D Motion Prediction}
\noindent \textbf{Deterministic Models.}
The goal of deterministic 3D motion prediction is to minimize the distance between a predicted motion and ground truth.
To solve this problem, early works relying on deep neural networks \cite{s-rnn, deter_1, deter_2} often employed recurrent neural networks (RNNs) \cite{lstm, gru}, which are still well-known for their effectiveness in processing time-series data.
Among RNN-based works, a notable model is a structure RNN (S-RNN) \cite{s-rnn}, which considers the spatio-temporal information of human motion, by manually designing the high-level spatio-temporal graph to explicitly model the human body structure (i.e., spine, arm, and leg).

While S-RNN understands the human body structure based on the handcrafted network structure, there is another line of research \cite{DCT-GCN, mao2020history} that uses graph convolutional network (GCN) to overcome this manually designed spatial relationship understanding. For instance, \cite{DCT-GCN} suggested a model named DCT-GCN, where discrete cosine transform (DCT) understands the temporal information of motion, and GCN learns the spatial relationship between human body joints. DCT-GCN obtains the state-of-the-art result when evaluated on Euler-angle-based mean squared error, but its best result can be obtained when the model is separately trained for each short- or long-term prediction.

Recently, several works for deterministic motion prediction \cite{ST-TR, 2CH-TR} are based on the Transformer \cite{transformer}, which was originally suggested for language understanding problems.
Models named Spatio-Temporal Transformer (ST-TR) \cite{ST-TR} and 2-Channel Transformer (2CH-TR) \cite{2CH-TR}, understand the spatio-temporal relationship of human motion by putting the self-attention mechanism on each pose-parameter (spatial) and time (temporal) dimension.
After understanding each spatial and temporal information in parallel, outputs from both attention mechanisms are properly combined. The difference between ST-TR and 2CH-TR comes from when and how often the model combines spatial and temporal information.

\vspace{2mm}
\noindent \textbf{Generative Models.}
The goal of stochastic 3D motion prediction is to build a generative model which can sample out several future motions that are likely to happen after the observed human motion.
To solve this problem, early works \cite{hp-gan, mt-vae} employed deep generative models such as variational autoencoders (VAEs) \cite{vae} or generative adversarial networks (GANs) \cite{gan}.
For instance, \cite{mt-vae} suggested a generative model based on the conditional VAEs, and showed that VAEs can sample out several future motions that are reasonable as well as diverse.
Compared to VAE, \cite{hp-gan} showed that GANs based on the Wasserstein loss function can be effectively used in stochastic motion prediction tasks.

While these works \cite{hp-gan, mt-vae} focused on exploring the potential of using deep generative models in stochastic motion prediction tasks, another line of works \cite{Dlow, gsps} focused on sampling out as much as diverse motions that can contain the most plausible motion at the same time.
For instance, \cite{Dlow} proposed to train a post-hoc model which can be attached to the pre-trained deep generative model. This post-hoc model maps a random variable to several latent vectors of the pre-trained generative model.
Based on the \textit{diversity-promoting prior}, the post-hoc model is trained to improve the diversity between samples, which can be obtained by decoding the mapped latent vectors.

Experiments in \cite{Dlow, gsps} evaluate the likelihood of prediction samples based on the \textit{minimum} distance between the samples and the ground truth(s).
They denote the prediction samples as plausible based on the sample that is closest to the ground truth(s).
However, this can make it difficult for users to choose the most plausible motion among the prediction samples, since all samples will not be distributed near the most plausible motion.
For instance, if the observed motion is a human sitting down and drinking something, \cite{Dlow} and \cite{gsps} can produce motion samples that predict the human suddenly standing up and starting discussing something with others.
As \cite{contexuallyplausible} has pointed out, we would like to also focus on the necessity of contextually plausible and diverse motion sampling. Therefore, our paper would evaluate the likelihood of prediction also based on the mean and standard deviation of distances between the samples and ground truth.

\subsection{Diffusion Probabilistic Models}
Diffusion probabilistic models \cite{ddpm} have become a new rising star in generative models after showing excellent performance in image synthesis. Especially, its performance on text-conditioned image synthesis \cite{imagen} makes researchers as well as the public in awe.
Diffusion models consider two processes: a forward process that slowly destructs the data sample by gradually injecting the random noise, and a reverse process that learns how to reconstruct the data sample by gradually denoising the random noise.
While the advantage of diffusion models can be empirically shown based on their performances, the disadvantage is the speed of their sampling process. If the reverse process includes $1000$ times of denoising processes, it means that the data sample can be obtained after feed-forwarding the random noise to the denoising network for $1000$ times. Of course, this disadvantage can be circumvented if the application does not require the prediction samples with low latency.

Aside from image generation tasks, nowadays researchers are suggesting to use diffusion models in various generation tasks, such as text-to-speech \cite{diff_tts}, text-to-sound \cite{diff_sound}, and video \cite{ho2022video}.
Focusing on motion-related tasks like ours, several works incorporate diffusion models in text-conditioned motion generation tasks \cite{kim2022flame, zhang2022motiondiffuse}.
For the motion of intelligence agents, \cite{janner2022planning} suggests using diffusion models to sample out trajectories for properly solving a given task.
In our paper, we use diffusion models in 3D human motion prediction tasks, but to the best of our knowledge, there is no attempt yet to use diffusion models in the 3D motion prediction task.
But we believe more researchers would involve in using diffusion models to answer this question -- can diffusion models be our new savior in any kind of data generation tasks?

\section{Method} \label{sec:met}
\subsection{Preliminaries}
We will provide a short description of diffusion probabilistic models first. Note that our description relies on \cite{ddpm} and \cite{csdi}, which provide a basis for our work.
\hfill \break \newline
\textbf{Diffusion Probabilistic Model.}
Let $\mathbf{x}^0 \sim q(\mathbf{x}^0)$ denote a data point sampled from its distribution $q$. In order to learn $p_\theta(\mathbf{x}^0)$ which can model $q(\mathbf{x}^0)$, diffusion probabilistic models consider two processes.
One is a \textit{forward process} which gradually deconstructs $\mathbf{x}^0$ by injecting a subtle Gaussian noise for $K$ times, such that $\mathbf{x}^0$ can be destroyed into $\mathbf{x}^1,\ldots,\mathbf{x}^K$,  where $p(\mathbf{x}^K)=\mathcal{N}(\mathbf{0}, \mathbf{I})$.
This process can be formulated as below, which is to follow a Markov chain $q(\mathbf{x}^k|\mathbf{x}^{k-1})$ for $K$ times:
\begin{align}
    q(\mathbf{x}^{1:K}|\mathbf{x}^0)&=\prod_{k=1}^K q(\mathbf{x}^k|\mathbf{x}^{k-1}) \\ q(\mathbf{x}^k|\mathbf{x}^{k-1})&=\mathcal{N}(\sqrt{1-\beta_k}\mathbf{x}^{k-1}, \beta_k\mathbf{I}),
\end{align}
where $\beta_k$ denotes a constant for a noise level. Note that $\mathbf{x}^k$ can be sampled from $\mathbf{x}^0$ directly with a closed-form solution:
\begin{equation}
    \mathbf{x}^k = \sqrt{\alpha_k}\mathbf{x}^0 + \sqrt{1-\alpha_k}\boldsymbol{\epsilon}, \;\; \boldsymbol{\epsilon}\sim \mathcal{N}(\mathbf{0, I}),
\end{equation}
where $\hat{\alpha}_k=1-\beta_k$ and $\alpha_k=\prod_{i=1}^k \hat{\alpha}_i$.

\begin{figure*}
    \centering
    \includegraphics[width=0.975\textwidth]{fig/net.png}
    \caption{Two designs of our Transformer-based motion denoiser. Inspired by ST-TR \cite{ST-TR}, 2CH-TR \cite{2CH-TR} and CSDI \cite{csdi}, our motion denoiser processes both spatial and temporal information in series (top) or in parallel (bottom). Here, $d$ and $t$ stand for the dimension of each pose-parameter and time, and TF stands for Transformer \cite{transformer}. Note that the positional encoding also involves adding a learnable vector that represents a diffusion step $k$ as \cite{csdi} suggests.}
    \label{fig:net}
\end{figure*}

Another is a \textit{reverse process}, which goal is to obtain $\mathbf{x}^0$ starting from $\mathbf{x}^K \sim \mathcal{N}(\mathbf{0, I})$, by gradually denoising $\mathbf{x}^K$.
This process can also be formulated as following a Markov chain $p_\theta(\mathbf{x}^{k-1}|\mathbf{x}^k)$ for $K$ times:

\begin{align}
    p_\theta(\mathbf{x}^{0:K}) &= p(\mathbf{x}^K) \prod_{k=1}^K p_\theta(\mathbf{x}^{k-1}|\mathbf{x}^k), \\
    p_\theta (\mathbf{x}^{k-1}|\mathbf{x}^k) &= \mathcal{N} \big( \mathbf{x}_{t-1}; \boldsymbol{\mu}_\theta(\mathbf{x}^k, k), \sigma^2(k)\mathbf{I} \big),
\end{align}
where $p(\mathbf{x}^K)=\mathcal{N}(\mathbf{0}, \mathbf{I})$.
To obtain $\boldsymbol{\mu}_\theta$ and $\sigma$, \cite{ddpm} suggests denoising diffusion probabilistic models (DDPM), which get $\sigma^2(k) = \frac{1-\alpha_{k-1}}{1-\alpha_{k}}\beta_k$, parameterize $\boldsymbol{\mu}_\theta$ with $\theta$, and sample $\mathbf{x}^{k-1}\sim p_{\theta}(\mathbf{x}^{k-1}|\mathbf{x}^k)$ as below:
\begin{align}
    \boldsymbol{\mu}_\theta(\mathbf{x}^k, k) &= \frac{1}{\sqrt{\hat{\alpha}_k}} \bigg( \mathbf{x}^k - \frac{\beta_k}{\sqrt{1-\alpha_k}} \boldsymbol{\epsilon}_\theta(\mathbf{x}^k, k) \bigg). \\
    \mathbf{x}^{k-1} &= \boldsymbol{\mu}_\theta(\mathbf{x}^k, k) + \sigma(k)\mathbf{z}, \;\;\mathbf{z} \sim \mathcal{N}(\mathbf{0, I}).
\end{align}

In practice, $\boldsymbol{\epsilon}_\theta$ is modeled with a neural network, and it learns how much to denoise from $\mathbf{x}^k$. To train this, \cite{ddpm} suggested a simplified loss function as below:
\begin{align}
    \mathcal{L}(\theta) &= \| \boldsymbol{\epsilon} - \boldsymbol{\epsilon}_\theta(\mathbf{x}^k, k)\|^2 \nonumber \\
    &= \|\boldsymbol{\epsilon} - \boldsymbol{\epsilon}_\theta(\sqrt{\alpha_t}\mathbf{x}^0 + \sqrt{1-\alpha_t}\boldsymbol{\epsilon}, k)\|^2.
\end{align}

In a training process, $k$ is randomly sampled to obtain $\mathcal{L}(\theta)$. For more details, please refer to \cite{ddpm} and \cite{csdi}.
\vspace{2mm}

\noindent \textbf{Conditional Diffusion Model.}
A conditional score-based diffusion model for imputation (CSDI) \cite{csdi} is proposed to solve a time-series imputation problem using diffusion models.
It adds conditional information $\mathbf{x}_{co}$ to eq. (4)-(5):
\begin{align}
    p_{\theta}(\mathbf{x}^{0:K}) &= p(\mathbf{x}^K) \prod_{k=1}^K p_\theta(\mathbf{x}^{k-1}|\mathbf{x}^{k}, \mathbf{x}_{co}),\\
    p_\theta(\mathbf{x}^{k-1}|\mathbf{x}^{k}, \mathbf{x}_{co}) &= \mathcal{N}(\mathbf{x}^{k-1}; \boldsymbol{\mu}_\theta(\mathbf{x}^{k}, k | \mathbf{x}_{co}), \sigma^2(k)\mathbf{I})
\end{align}

To define $\boldsymbol{\mu}_\theta(\mathbf{x}^{k}, k | \mathbf{x}_{co})$,
eq. (6)-(7) can be rewritten by adding $\mathbf{x}_{co}$ as a condition to $\boldsymbol{\mu}_\theta$ and $\boldsymbol{\epsilon}_\theta$.
Note that $\boldsymbol{\epsilon}_\theta(\mathbf{x}^{k}, k | \mathbf{x}_{co})$ is modeled with a neural network to learn how much to denoise from $\mathbf{x}^k$ given $\mathbf{x}_{co}$.
When training the network, the same loss function as eq. (8) is used, by replacing $\boldsymbol{\epsilon}_\theta$ properly with $\mathbf{x}_{co}$ as a condition.

\subsection{Problem Formulation}
Let $\mathbf{p}_t \in \mathbb{R}^{D}$ be a 3D pose vector at time $t$, which can be denoted with various representations such as axis-angle, Euler-angle, or $xyz$-position. Here, $D=3n$ and $n$ denotes the number of joints.
A task of 3D human motion prediction can be defined as predicting future $L$ poses,
${P}_{pre}=\{\mathbf{p}_{T+1},\ldots \mathbf{p}_{T+L}\} \in \mathbb{R}^{L \times D}$,
when $T$ poses, ${P}_{obs}=\{\mathbf{p}_1,\ldots \mathbf{p}_{T}\} \in \mathbb{R}^{T \times D}$ are observed.

We utilize CSDI \cite{csdi} for obtaining $P_{pre}$ from given $P_{obs}$.
Starting from $P_{pre}^0=P_{pre}$, our forward process can obtain $P_{pre}^k$ as below:
\begin{equation}
    P_{pre}^k = \sqrt{\alpha_k} P_{pre}^0 + \sqrt{1-\alpha_k} \boldsymbol{\epsilon}, \;\; \boldsymbol\epsilon \sim \mathcal{N}(\mathbf{0, I})
\end{equation}
For a reverse process, we propose a denoiser network which models $\boldsymbol{\epsilon}_\theta(\mathbf{x}^k,k| \mathbf{x}_{co})=\boldsymbol{\epsilon}_\theta(P_{pre}^k, k | P_{obs})$.
This network is trained by minimizing
$\mathcal{L}(\theta) = \| \boldsymbol{\epsilon} - \boldsymbol{\epsilon}_\theta(P_{pre}^k, k| P_{obs}) \|^2$.

After training, we can sample $P_{pre}^0$ by repeating below reverse process for $K$ times, starting from $P_{pre}^K \sim \mathcal{N}(\mathbf{0, I})$:
\begin{equation}
    P_{pre}^{k-1} = \boldsymbol{\mu}_\theta(P_{pre}^k, k | P_{obs}) + \sigma(k)\mathbf{z}, \;\; \mathbf{z}\sim\mathcal{N}(\mathbf{0, I}),
\end{equation}%
where $\boldsymbol{\mu}_\theta(P_{pre}^k, k | P_{obs})$ is defined with $\boldsymbol{\epsilon}_\theta(P_{pre}^k, k | P_{obs})$ and properly modified version of eq. (6).
After finishing training, if our denoiser network is used for deterministic prediction in the test phase, we set $P_{pre}^K$ and $\mathbf{z}$ as zero-vectors, such that all randomness in eq. (12) can be ignored.

\subsection{Transformer-based Motion Denoiser}
Since $P_{pre}^k$ and $P_{obs}$ are time-series of human pose vectors, one can model
$\boldsymbol{\epsilon}_\theta(P_{pre}^k, k | P_{obs})$ with neural network architectures which can understand time-series data.
For example, network architectures such as RNNs \cite{lstm, gru} or Transformers \cite{transformer} can be candidates.
We empirically found out that the denoisers based on the Transformers that process both spatial and temporal information are most effective.

Figure~\ref{fig:net} shows how we design our Transformer-based denoisers in two ways.
Inspired by \cite{csdi}, the first denoiser shown on top of the figure processes both information in series. After concatenating $P_{pre}^k \in \mathbb{R}^{L \times D}$ and $P_{obs} \in \mathbb{R}^{T \times D}$ such that input can be $P_{inp}^k \in \mathbb{R}^{(T+L) \times D}$, $P_{inp}^k$ passes spatial and temporal transformer layers in series, where each layer applies self-attention to time and pose-parameter dimension.
Before passing each transformer layer, positional encoding is added to the input as \cite{transformer} suggests, with respect to pose-parameter $d\in[0, D]$ (spatial) or time $t\in[0, T]$ (temporal) dimension.
Also, the additional learnable positional encoding that projects a diffusion step $k$ into a vector space is added to the input as \cite{csdi} suggests.
Let $P_{out}^k\in \mathbb{R}^{(T+L) \times D}$ denote the output which can be obtained after $P_{inp}^k$ passing two layers.
Then, the last $L \times D$ parts from $P_{out}^k$ is obtained as $\boldsymbol\epsilon_\theta(P_{pre}^k, k|P_{obs})$, which would be used for denoising $P_{pre}^k$.

The second denoiser shown on the bottom of Figure~\ref{fig:net}, is inspired by \cite{ST-TR} and \cite{2CH-TR}, and works in parallel to understand spatio-temporal information.
After $P_{inp}^k$ passes both spatial and temporal transformer layers in parallel, two matrices with the same size as $P_{inp}^k$ are obtained, and concatenated into a 3rd-order tensor whose size is $2\times (T+L) \times D$. After this tensor passes 2-dimensional convolutional layer with $(1\times 1)$-sized kernel, the output $P_{out}^k\in \mathbb{R}^{(T+L) \times D}$ is obtained.
From $P_{out}^k$, $\boldsymbol\epsilon_\theta(P_{pre}^k, k|P_{obs})$ is obtained as same as in the first denoiser.

Note that we do not use encoder-decoder based structure, which encode a set of feature vectors from $P_{obs}$ and decode $\boldsymbol\epsilon_\theta(P_{pre}^k, k|P_{obs})$ from the encoded feature vectors and $P_{pre}^k$.
We tried various denoisers of Transformer- or RNN-based encoder-decoder, but none of them turns out to be effective.

\subsection{Implementation Details}
Our transformer-based motion denoisers have a self-attention module with 8 multi-heads and 512-dimensional query, key, and value vectors. And each temporal or spatial transformer layer shown in Figure~\ref{fig:net} consists of a single-layered transformer encoder.
To train denoisers, we set batch size as 512 and update parameters for 50,000 iterations with Adam optimizer of learning rate $0.0001$.
The diffusion step is set as $k\in[0, 20]$, with linearly scheduled noise levels $\beta_k$ that ranges between $0.001$ ($k\downarrow$) and $0.333$ ($k\uparrow$).

\section{Experiment} \label{sec:exp}
\subsection{Dataset and Metric}
\noindent \textbf{Dataset.}
We conduct our experiment for both deterministic and stochastic motion prediction tasks. For deterministic experiments, we use the Human3.6M dataset \cite{h36m_pami} and measure the Euler-angle mean square error (MSE) for evaluation as other works \cite{s-rnn, DCT-GCN, ST-TR, 2CH-TR} do. Here, with 25 fps, input observation has 50 frames, and output prediction has 25 frames.
For stochastic experiments, we preprocess Human3.6M \cite{h36m_pami} and HumanEva-I \cite{sigal2006humaneva} datasets into $xyz$-based representation as \cite{Dlow, gsps} do. Based on that, various metrics for evaluating likelihood and diversity are measured. Here, with 50 fps, an input observation has 25 frames, output prediction has 100 frames, and the number of prediction samples is 50.

\vspace{2mm}

\noindent \textbf{Metrics.}
As mentioned above, we measure the performance of our denoiser based on the Euler-angle MSE when it is used for deterministic prediction.
For stochastic prediction, we use several metrics from what \cite{Dlow} suggests to evaluate likelihood and diversity. But we propose more metrics such as aDE, sDE, aFDE, and sFDE to measure how the samples are distributed near the ground truth. Note that some of the below sentences describing metrics are borrowed from \cite{Dlow}.

(1) \textbf{Average Pairwise Distance (APD)}: average $L2$ distance between pairs from $N$ predictions $\hat{\mathbf{x}} \in \mathbb{R}^{L \times D}$, which is computed as $\frac{1}{N(N-1)}\sum_{i=1}^N \sum_{j\neq i}^N \| \hat{\mathbf{x}}_i - \hat{\mathbf{x}}_j \|_2$. This measures the diversity within $N$ predictions.
(2) \textbf{minimum Displacement Error (mDE)}: the minimum $L2$ distance between all $N$ predictions $\hat{\mathbf{x}}$ and ground truth $\mathbf{x}$, which is computed as $ \min_{\hat{\mathbf{x}}}\frac{1}{L}\|\hat{\mathbf{x}}- \mathbf{x} \|_2$.
This metric was defined as ADE in \cite{Dlow}.
(3) \textbf{average Displacement Error (aDE)}: the average $L2$ distance between all $N$ predictions $\hat{\mathbf{x}}$ and ground truth $\mathbf{x}$, which is computed as $\frac{1}{NL} \sum_{i=1}^N\|\hat{\mathbf{x}}_i- \mathbf{x} \|_2$.
(4) \textbf{standard deviation of Displacement Error (sDE)}: the standard deviation of $L2$ distances between all $N$ predictions and ground truth.
(5) \textbf{minimum Final Displacement Error (mFDE)}: the minimum $L2$ distance between final poses of $N$ predictions and ground truth, which is calculated as $\min_{\hat{\mathbf{x}}}\|\hat{\mathbf{x}}(L) - \mathbf{x}(L)\|_2$. This metric was defined as FDE in \cite{Dlow}.
(6) \textbf{average Final Displacement Error (aFDE)}: the average $L2$ distance between final poses of $N$ predictions and ground truth, which is calculated as $\frac{1}{N}\sum_{i=1}^N\|\hat{\mathbf{x}}_i(L) - \mathbf{x}(L)\|_2$.
(7) \textbf{standard deviation of Final Displacement Error (sFDE)}: the standard deviation of $L2$ distances between final poses of $N$ predictions and ground truth.

\subsection{Quantitative Results}
\noindent \textbf{Deterministic Prediction.} Table~\ref{tab:mse} compares Euler-angle MSEs when our diffusion model is used for deterministic motion prediction.
Here, bold fonts denote the best results among all approaches, and underlines denote the best results among our denoisers (series or parallel).
It is shown that the overall performance of DCT-GCN \cite{DCT-GCN} is still the best.
Among our approaches, the denoiser which understands spatial and temporal information in series is better than the parallel denoiser.
Although our models do not achieve state-of-the-art results, it is shown that our approaches are better in long-term prediction (1000ms) when compared with other transformer-based models \cite{ST-TR, 2CH-TR}.
This is a notable result, since (1) our models are originally generative ones, and (2) our models do not require additional training for deterministic prediction since ignoring all randomness in the denoising process is all they need.

\begin{table}
\centering
\caption{Average MSE errors of Deterministic Motion Prediction}
\setlength\tabcolsep{5.5pt}
\begin{tabular}{l|llllll}
\textbf{milisecond (ms)}       & \multicolumn{1}{c}{\textbf{80}} & \multicolumn{1}{c}{\textbf{160}} & \multicolumn{1}{c}{\textbf{320}} & \multicolumn{1}{c}{\textbf{400}} & \multicolumn{1}{c}{\textbf{560}} & \multicolumn{1}{c}{\textbf{1000}}  \\
\hline
S-RNN \cite{s-rnn}                          & 0.933                           & 1.166                            & 1.397                            & 1.526                            & 1.711                            & 2.139                              \\
DCT-GCN \cite{DCT-GCN}                       & {0.295}                           & \textbf{0.542}                            & \textbf{0.857}                            & \textbf{0.974}                            & \textbf{1.154}                            & \textbf{1.590}                              \\
ST-TR \cite{ST-TR}          & 0.303                           & 0.550                            & 0.901                            & 1.021                            & 1.229                            & 1.722                              \\
2CH-TR \cite{2CH-TR}                  & \textbf{0.293}                           & 0.555                            & 0.893                            & 1.016                            & 1.245                            & 1.744                              \\
\hline
\textbf{Ours (Series)} &\underline{0.325}                           & \underline{0.615}                            & \underline{0.990}                            & \underline{1.128}                            & \underline{1.309}                            & 1.721
\\
\textbf{Ours (Parallel)} & 0.350                           & 0.646                            & 1.007                            & 1.148                            & 1.317                            & \underline{1.688}
\end{tabular}
\vspace{-17mm}
\label{tab:mse}
\end{table}

\begin{table*}
\centering
\caption{Diversity and Likelihood Metrics of Stochastic Motion Prediction}
\setlength\tabcolsep{4pt}
\begin{tabular}{l|lllllll|lllllll}
                         & \multicolumn{7}{c|}{Human 3.6M \cite{h36m_pami}}                                                                                & \multicolumn{7}{c}{HumanEva-I \cite{sigal2006humaneva}}                                                                                \\
\cline{2-15}
metrics                  & \multicolumn{1}{c}{APD$\uparrow$} & \multicolumn{1}{c}{mDE$\downarrow$} & \multicolumn{1}{c}{aDE$\downarrow$} & \multicolumn{1}{c}{sDE$\downarrow$} & mFDE$\downarrow$           & aFDE$\downarrow$           & sFDE$\downarrow$           & APD$\uparrow$            & mDE$\downarrow$           & aDE$\downarrow$           & sDE$\downarrow$           & mFDE$\downarrow$           & aFDE$\downarrow$           & sFDE$\downarrow$            \\
\hline
DLow \cite{Dlow}                    & \underline{11.741}          & \underline{0.425}            & 0.968                    & 0.355                    & \underline{0.518}  & 1.387          & 0.541          & \underline{4.855}  & 0.251          & 0.585          & 0.208          & 0.268          & 0.710          & 0.255           \\
VAEs \cite{vae, Dlow}                     & 6.852                   & 0.460                    & \underline{0.720}            & \underline{0.139}            & 0.557          & \underline{1.025}  & 0.243          & 2.299          & 0.265          & 0.426          & 0.083          & 0.299          & 0.562          & 0.137           \\
GSPS \cite{gsps}                    & \textbf{14.757}         & \textbf{0.389}           & 1.206                    & 0.623                    & \textbf{0.496} & 1.554          & 0.729          & \textbf{5.825} & \textbf{0.233} & 0.655          & 0.206          & \underline{0.244}  & 0.763          & 0.268           \\
\hline
\textbf{Ours (Series)}   & 7.587                   & 0.527                    & 0.764                    & \textbf{0.132}           & 0.669          & 1.093          & \textbf{0.228} & 2.746          & 0.257          & \underline{0.383}  & \underline{0.065}  & 0.260          & \underline{0.490}  & \underline{0.130}  \\
\textbf{Ours (Parallel)} & 6.445                   & 0.477                    & \textbf{0.719}           & \underline{0.139}            & 0.584          & \textbf{1.018} & \underline{0.234}  & 1.508          & \underline{0.242}  & \textbf{0.312} & \textbf{0.037} & \textbf{0.238} & \textbf{0.385} & \textbf{0.078}
\end{tabular}
\label{tab:sto}
\vspace{-3mm}
\end{table*}

\noindent \textbf{Stochastic Prediction}
Table~\ref{tab:sto} shows the comparison of metrics for measuring the likelihood and diversity.
Here, bold fonts denote the best result and underlines denote the second best result among all approaches.
It is shown that previous works \cite{Dlow, gsps} focusing on sample diversity best perform in APD.
Also, it is shown that they are generally better in terms of mDE and mFDE.
We would like to argue here that the high diversity in prediction increases the probability of having one sample closest to the ground truth.
Then, how can we choose the most plausible result among predictions that are sampled to be diverse?

This is the same question that \cite{contexuallyplausible} also pointed out.
So in \cite{contexuallyplausible}, metrics for measuring the quality and context are proposed.
For measuring the quality, \cite{contexuallyplausible} used a pre-trained binary classifier which can discriminate the ground truths (real) from predictions (fake). If this classifier fails to discriminate the predicted motions as fake, a higher quality score is obtained.
For measuring the context, \cite{contexuallyplausible} used a pre-trained model which classifies action from motion. If it estimates that the action label of prediction is as same as the observed motion, a higher context score is obtained.

However, we were not able to use the same metric as \cite{contexuallyplausible} since its pre-trained classifiers were not openly released.
Therefore, we instead propose metrics such as aDE, sDE, aFDE, and sFDE, to measure how closely the samples are distributed near the ground truth.
Results show that our approaches generally perform better in terms of these new metrics, and the parallel denoiser performs better than the series one.
We also present the result from VAEs \cite{vae} that were implemented by \cite{Dlow}, to check how other non-diffusion generative models work.
It is shown that the overall performances of our series/parallel denoiser in diversity and likelihood are generally better than the VAEs, especially in the HumanEva-I dataset.

\begin{figure*}
    \centering
    \includegraphics[width=0.95\textwidth]{fig/qual.png}
    \caption{Deterministic (Deter.) and stochastic (Sto.) predictions from our transformer-based motion denoiser. Note that two results are given and divided based on the vertical dotted line. Predictions are obtained from observed motions labeled as `smoking' (left) and `walking' (right).}
    \label{fig:qual}
    \vspace{-3mm}
\end{figure*}

\subsection{Qualitative Results}
Figure~\ref{fig:qual} shows two example results from our transformer-based motion denoiser. Predictions on the left of the dotted line are obtained from the motion observation labeled as `smoking'. It is shown that the deterministic prediction is similar to the ground truth, while the stochastic predictions show the diversity between samples. But note that still the context of `smoking' looks remained in all samples. This phenomenon is also observed from the predictions on the right, which are obtained from the motion observation of `walking'. While its deterministic prediction resembles the ground truth, the stochastic predictions are diverse and contain the context of `walking'.
For better visualization, please refer to our supplementary video.

\section{Conclusion} \label{sec:con}
In this work, we study the potential of diffusion probabilistic models for 3D human motion prediction tasks.
We propose two types of diffusion models based on the transformers, which understand the motion's spatial and temporal information in series or parallel.
Since the diffusion model is originally a generative model, its main usage would be for the stochastic motion prediction task.
But once it is trained, we show that it can also be used in deterministic prediction if all randomness in its denoising process is ignored.

To show the effectiveness of diffusion models in both deterministic and stochastic motion prediction tasks, we conduct experiments based on various metrics.
Results from deterministic prediction show that the diffusion model is not superior to the state-of-the-art. But it is shown that our long-term (1000ms) prediction performance is better than other transformer-based approaches.
When it comes to evaluating stochastic predictions, it is conventional to suggest metrics measuring both likelihood and diversity.
However, we claim that the conventional metrics for measuring the likelihood do not represent how much the samples are distributed near the plausible motion, since they measure the \textit{minimum} distance between samples and ground truth.
Therefore, we suggest additional metrics to measure the mean and standard deviation of that distances, and the results show that our diffusion models can properly balance the trade-off between diversity and likelihood.

Although our results would provide nice answers to our first question -- can we use diffusion probabilistic models for 3D motion prediction? -- the most concerning disadvantage of a diffusion model is its sampling frequency.
Since our diffusion model requires a $K=20$ number of denoising processes to obtain prediction samples, this might occur a bit high latency. To overcome this issue, one might consider recent works for efficient sampling \cite{salimans2021progressive}, which would be our future work, such that efficient 3D human motion prediction can be made for various real-time applications.

\section*{ACKNOWLEDGMENT}
This work was supported by the Institute of Information \& communications Technology Planning \& Evaluation (IITP) grant funded by the Korea government (MSIT) (No.2020-0-01336, Artificial Intelligence Graduate School Program (UNIST)),
and funded by Marie Sklodowska-Curie Action Horizon 2020 (Grant agreement No. 955778) for project `Personalized Robotics as Service Oriented Applications' (PERSEO).

\addtolength{\textheight}{-4cm}