Long-term double synchronization in close-in gas giant planets

Hot Jupiters, orbiting their host stars at extremely close distances, undergo tidal evolution, with some being engulfed by their stars due to angular momentum exchanges induced by tidal forces. However, achieving double synchronization can prolong their survival. Using the MESA stellar evolution code, combined with the magnetic braking model of Matt et al. (2015), we calculate 25,000 models with different metallicity and study how to attain the conditions that trigger the long-term double synchronization. Our results indicate that massive planets orbiting stars with lower convective turnover time are easier to achieve long-term double synchronization. The rotation angular velocity at the equilibrium point ($\Omega_{\mathrm{sta}}$) is almost equal to orbital angular velocity of planet ($\mathrm{n}$) for the majority of the main sequence lifetime if a system has undergone a long-term double synchronization, regardless of their state at this moment. We further compared our results with known parameters of giant planetary systems and found that those systems with larger planetary masses and lower convective turnover time seem to be less sensitive to changes in the tidal quality factor $Q'_{_*}$. We suggest that for systems that fall on the state of $\Omega_{\mathrm{sta}} \approx n$, regardless of their current state, the synchronization will persist for a long time if orbital synchronization occurs at any stage of their evolution. Our results can be applied to estimate whether a system has experienced long-term double synchronization in the past or may experience it in the future.


INTRODUCTION
In a close-in binary system, tidal interactions between the two stars tend to align their rotational axes, circularize their orbits, and synchronize their rotations (Zahn 1977;Hut 1981;Ogilvie 2014).The importance of tidal effects on hot Jupiters was first highlighted by Rasio & Ford (1996) after the discovery of the first hot Jupiter, 51Pegasi b (Mayor & Queloz 1995).The Kepler, CoRoT, TESS and ground-based telescopes have enabled astronomers to have a deeper insight into the exoplanet systems.In the past thirty years, several hundred hot Jupiters have been discovered around main-sequence FGK stars.Generally, when the orbital period is less than 10 days, their orbital eccentricities are very low (Winn & Fabrycky 2015).According to current theories of planet formation, it is difficult to accumulate the necessary material to form close-in planets around their stars (Lin et al. 1996).Therefore, a possibility for the formation of close-in exoplanets is that in the early stages of planet formation, hot Jupiters migrate to short orbits due to the drag of the gaseous protoplanetary disk (Lin et al. 1996;Trilling et al. 1998;Ward 1997;Papaloizou 2007).
In various orbital configurations, the double synchronous orbits denote the orbital configuration where the rotational period of a star is nearly equal to the orbital period of its orbiting planet.In the case of conservation of total angular momentum, Hut (1980) proves that an equilibrium status of co-planarity, circularity and co-rotation can be attained for a binary system.However, the total angular momentum of such the system can not be conserved because the magnetic braking caused by the stellar winds takes the angular momentum out the system (Matt et al. 2015).Carone (2012) suggested that double synchronous orbits may also decay due to magnetic braking.Therefore, only F-type stars with a massive companion can maintain double synchronous states owing to the weak effect of magnetic braking.The loss of stellar angular momentum can be compensated to some extent by the tidal friction between the star and planet.In fact, the angular momentum will be transferred to the star if the orbital angular velocity of the planet is larger than the angular velocity of the star.Thus, the competition between the magnetic braking and tidal friction determines when an orbital configuration of co-rotation can occur (Damiani & Lanza 2015).
It is important to study the occurrence of the double synchronous orbit (or co-rotation) because the presence of quasi-stable states can significantly delay their tidal evolution and extend the life of the planet, which would otherwise cause the planet to fall into its host star (Damiani & Lanza 2015).There is observational evidence suggesting that close-in companion stars around F-type stars may have larger masses compared to those around hosts of lower mass.Bouchy et al. (2011) proposed an explanation for this, stating that sufficiently massive companions and fast-spinning primary stars can achieve a tidal equilibrium state with synchronized orbital motion and stellar rotation.Due to weaker magnetic braking, F-type stars can maintain fast rotation throughout their main sequence phase, allowing companion stars to survive for a longer period.Damiani & Díaz (2016) also found, by calculating the timescale for orbital decay, that F-type stars may have significantly longer durations of hosting high-mass companions compared to G-type stars, given a certain orbital period.More recently, Benni et al. (2021) discovered a short-period brown dwarf system, GPX-1, orbiting a rapidly rotating F-type young star.The authors found that the rotation period of GPX-1 is synchronized with the orbital period of the brown dwarf, and the system is in a pseudo equilibrium state and stable.In F-type stars, the loss of angular momentum due to stellar wind can be counteracted by the tidal torque induced by the close-in planet, preventing significant changes for the stellar rotation and planetary orbit.This makes it easier to maintain a double synchronization state in a longer duration.Amard & Matt (2020) found that the evolution of the rotation period of solar-like stars varies with different metallicity.Stars with low metallicity tend to rotate faster, due to their thin convective envelopes and weak magnetic braking.F-type stars, like G-type stars with low metallicity, have thin convective envelopes that result in weaker magnetic braking.This suggests that metallicity affects the balance between wind torque and tidal torque, which in turn affects orbital stability.Furthermore, the tidal interactions between the star and planet are also affected by the properties of the planet, such as planetary mass and separation.However, there is lack of sufficient investigation for the co-rotation of planets around late-type stats.It is not clear how to trigger double synchronization around late-type stars with different metallicities.Thus, these findings motivate us to investigate the phenomenon of the double synchronization further.Our aim is to explore the conditions that trigger the double synchronous.
In this paper, we focus on hot Jupiters orbiting G and F types stars.We replaced the Skumanich-type law (Skumanich 1972) with the magnetic braking scheme of Matt et al. (2015) and apply the criterion of Damiani & Lanza (2015) to determine if the star-planet system is in the state of co-rotation.We briefly described in section 2.1 the grid spacing we used in our MESA calculations.We reviewed previous studies in Section 2.2 and described the conditions for pseudo-stable equilibrium of hot Jupiters given by Damiani & Lanza (2015).In Section 2.3, we discussed the stellar spin rate Ω sta at tidal and wind torque balance based on the magnetic braking model of Matt et al. (2015).In Section 3.1, we examined the impact of stellar mass, metallicity, and initial spin rate on triggering double synchronization.In Section 3.2, we discussed the impact of planetary mass, tidal quality factor, and initial semi-major axis on triggering double synchronization.In Section 3.3, we quantitatively present the distribution of initial parameters for the occurrence of long-term double synchronization and discuss them.In Section 4, we have explored the distribution of 243 observational samples within the Darwinian diagram and the phases during which long-term double synchronization occurs.

Star-planet interaction model in MESA
In our previous work, we incorporated the process of angular momentum exchange between stars and planets into the stellar evolution code MESA, version 11554 (Guo 2023).We did not consider the differential rotation of the star, which is reasonable.Taking the Sun as an example, the error between considering and not considering the differential rotation does not exceed 5%.tectbfThe detailed calculation procedures for this are provided in Appendix A. MESA has been widely utilized in the field of stellar and planetary evolution (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019)).The values of the parameters we adopted in this study are summarized in Table 1.It is worth noting that the initial parameters for the nearly 25,000 models we calculated are all possible combinations of the six parameters listed in Table 1.

Pseudo-stable equilibrium of hot Jupiters
In this section we summarize Damiani & Lanza (2015) results.Hut (1980) found that in binary systems, when the total angular momentum of the system exceeds a critical value, i.e.  >  c 0 , the system will reach a double synchronous rotation state, which can prevent the planet from further inward migration.Subsequently, Damiani & Lanza (2015) found the conditions for maintaining a pseudostable equilibrium state when taking into account the stellar magnetic braking.They derived the conditions for tidal pseudo-equilibrium under magnetic braking and showed that the system evolves towards a circular and aligned orbit with minimum energy.The pseudo-stable condition is given by: (1) where  and Ω are the angular velocity of the planet and the star, respectively. is the semi-major axis,  * and  pl are the stellar and planetary masses.n is the mean orbital motion, and  is a function of time defined as: Due to the presence of magnetic braking, the total angular momentum of the system gradually decreases, i.e. dL dt < 0. When  = 1, the system has Ω = n, and the star and planet are in double synchronization.At this point, dΩ d ≈ 0, and the total angular momentum of the system is approximately conserved, equivalent to the case without considering magnetic braking.When the tidal torque is negligible compared to magnetic breaking,  ≈ 0. When 0 <  < 1, the system has Ω < , and dΩ d < 0 is required for the equation to hold, i.e. the tidal torque is less than the wind torque.When  > 1, the planet migrates outward under tidal action, and both tidal torque and wind torque slow down the stellar rotation.However, the equation requires dΩ dt > 0, which is clearly not true.Therefore, only when  < 1, the equilibrium state can be reached.Moreover, the equilibrium is stable only if the total angular momentum of the system  >  c and  orb >  orb,s , where (5) where ,  c ,  orb and  orb,s are the total angular momentum, critical total angular momentum, orbital angular momentum and critical orbital angular momentum respectively. * and  sp are the stellar and planetary moment of inertia.When  = 1, the equation degenerates to the case without magnetic braking.In this case, we have: In the presence of magnetic braking, it can be demonstrated that  s , defined as the maximum orbital frequency allowing for the existence of pseudo-stable states, is determined by: At  =  c 0 , it corresponds to the only average motion of rotation At  =  s the total angular momentum is When 0 <  < 1, we have 3 4 3 4  c 0 <  s <  c 0 .The corresponding conditions for the existence of a pseudo-stable state under total angular momentum can be demonstrated to satisfy  c <  s <  c 0 .
Therefore, as long as the mean motion of the orbit is  < 3 4 3 4  c 0 and  >  c 0 , the system can achieve a pseudo-stable equilibrium.The Roche limit is defined as where  p = 1.3  J and the rotational angular momentum of the planet is negligible compared to its orbital angular momentum, we can neglect  sp .

Balance of tidal torque and wind torque
In this section, we did not use the Skumanich-type law (Skumanich 1972), which assumes that the magnitude of the torque due to magnetic braking is  mb =  mb Ω 3 , where  mb = 1.5 × 10 −14  yr and  = 1.0 for G stars and  = 0.1 for F stars.Instead, we followed the work of Matt et al. (2015) and used their braking equation to find the equilibrium position Ω sta of the star, where the tidal and wind torques balance.This equation takes into account not only a simple constant multiplied by Ω 3 , but also the thickness of the convective envelope, which determines the magnetic braking resulting from the stellar mass and metallicity.This approach allows us to more accurately reflect the relationship between star mass and magnetic braking, rather than simply distinguishing F and G stars based on their  values.We can also discuss the equilibrium position Ω sta of stars with different metallicity.Ω sta represents the physically meaningful solution when the torque is balanced, i.e., when dΩ/d = 0. Following Goldreich (1963); Kaula (1968); Jackson et al. (2008), the time evolution of the semi-major axis can be expressed as The total angular momentum of the system: Here,  * ,  orb and  wind are the stellar angular momentum, planetary orbital angular momentum, and the angular momentum lost by the magnetic stellar wind, respectively.They are given by: Due to the relatively small variation in stellar moment of inertia during the main sequence for solar-like stars, the second term on the right side of Equation ( 17) is much smaller than the first term.For simplicity, we neglect the second term in Equation ( 17).By substituting Equation (15) into Equation ( 18) and utilizing the above equations, we can solve for Equation (21) and Equation ( 22).The negative sign on the left side of the equations represents the tidal effects, while the negative sign on the right side represents the influence of magnetic braking. (unsaturated) where  * , Ω, and Ω crit are the stellar radius, angular rotation rate, and critical angular rotation rate, respectively.The constants , , , and  are free parameters.The saturated and unsaturated states in Equations ( 19) and ( 20) represent two different states of stellar magnetic activity, both of which have a strong correlation with the Rossby number  o = 2  Ω cz .The convective turnover time  cz can be expressed as (Gossage et al. 2021): where  P () is the scale height,  c () is the convective velocity at radius , and  MLT is the convective mixing length.In Equation ( 26),  cz is the turnover timescale of the convective zone, which is defined in the location where  =  BCZ + 0.5  P (), and  BCZ is the radius of the bottom of the outer convection zone.The values of  MLT is 1.82. o ⩽  osat is the saturation region, and vice versa.The parameter  =  o / osat relies on the critical value  osat = 0.14 (Wright et al. 2018).The solar  o is around 2 (See et al. 2016).Thus we take  = 14.The constant  is set to be 0.22.In order to reproduce the current rotation period of the Sun, the values of  and  are 1.2 × 10 30  and 2.6, respectively.
Here, the equilibrium point where the tidal and wind torques are balanced corresponds to dΩ/dt = 0.By solving a high-degree onedimensional function, we can obtain Ω sta at the equilibrium point.However, unlike the cubic function in Damiani & Lanza (2015) that can be solved using Cardano's method, we can only seek numerical solutions for high-degree one-dimensional functions.We provide a method for finding the roots of high-degree one-dimensional functions in Appendix B.

The influence of stellar parameters on the triggering of double synchronization
In Figure 1, we discuss the effects of different initial stellar rotation periods, metallicity, and stellar masses of the star-planet systems on the triggering of the double synchronization mechanism.In Figure 1 (a), we show the evolution of the stellar rotation period and the planetary orbital period, the ratio of orbital angular momentum to critical orbital angular momentum for the system with  * = 1.1  ⊙ ,  pl = 8.0  J and  ini =0.05 au when the value of  ′ * is 10 6 .We also show the ratio of total angular momentum to critical total angular momentum as a function of age.For top panel: the solid colored lines represent the stellar rotation periods  rot , the dashed colored lines represent the planet orbital periods  orb , and the thin solid black line represent the Roche limit., respectively.The square markers represent the 99% position of the system's evolutionary timescale, while the pentagon markers indicate the end point of the system's evolution.The magnified inset is included to provide a clearer view for the example with [Fe/H] = -0.5 dex and  rot,ini = 3 .In this scenario, the planetary system, due to the double synchronization effect, can survive until the end of the stellar main-sequence phase.
The top panel of Figure 1 (a) shows that long-term double synchronization is possible when [Fe/H] = -0.5 dex, because the stellar magnetic braking is weak.For example, for the cases with initial rotation periods of 3.0 days and 8.0 days, long-term double synchronization can be maintained for a long time when the stellar rotation rate synchronizes with the planetary orbital rate.However, for [Fe/H] = 0 dex and +0.5 dex, for the case with initial rotation period of 8.0 days, the stellar rotation period and the planetary orbital period approach synchronization only when the planet is about to be engulfed, and for the case with initial rotation period of 3.0 days, due to the effect of magnetic braking in the early stage of evolution and the transfer of stellar rotation angular momentum to the planetary orbit, the stellar rotation period rapidly increases until the stellar rotation rate and the planetary orbital rate synchronize, but this synchronized state is immediately broken.This is because in this case, the stars with [Fe/H] = 0 dex and +0.5 dex have stronger magnetic braking than the stars with [Fe/H] = -0.5 dex, and long-term double synchronization is difficult to maintain.Moreover, long-term double synchronization requires that Ω sta ≈  = Ω, which can also be seen in Figure 1 (b), where the systems with [Fe/H] = 0 dex and +0.5 dex do not satisfy Ω sta < (red and blue dash-dotted lines), long-term double synchronization cannot be maintained, while the system with [Fe/H] = -0.5 dex satisfies Ω sta ≈ (green dash-dotted line), long-term double synchronization can be maintained.It is worth noting that from the middle and bottom panels of Figure 1 (a), all models undergo orbital decay and magnetic braking, the  orb  orb,c 0 and   c 0 ratios of the system with initial stellar rotation periods of 3.0 days and [Fe/H] = -0.5 dex always remain above 1, while the other curves will eventually lead to  orb  orb,c 0 and   c 0 ratios below 1 as they evolve.When  orb  orb,c 0 < 1, the planets fall into the star quickly and are swallowed.This is because the total angular momentum of the system is always decreasing due to the effect of stellar magnetic braking, and the higher the stellar metallicity and the thicker the convective envelope, the stronger the magnetic braking, so that the system is unstable when  orb  orb,c 0 and   c 0 are less than 1.For [Fe/H] = -0.5 dex, we find the stellar-planet system with an initial period of 3.0 days has a larger initial angular momentum, the angular momentum loss caused by magnetic braking is insufficient to decrease the total angular momentum below the critical value during the main sequence phase of the star from the bottom panel of Figure 1 (a).
Meanwhile, from the Darwin diagram in Figure 1(b), for [Fe/H] = -0.5 dex and  rot,ini = 3.0 , when  > 1, the planet migrates outward, resulting in a decrease in the orbital velocity during the evolution.Both tidal torque and wind torque contribute to the reduction of the stellar angular momentum, until the system reaches double synchronization at  = 1.On the other hand, for [Fe/H] = -0.5 dex and  rot,ini = 8.0 , when 0 <  < 1, the tidal torque causes the planet to migrate inward, leading to an increase in the orbital velocity during the evolution.Initially, the tidal torque dominates the wind torque, and the angular momentum remains approximately constant, staying nearly horizontal during the early stages of evolution.However, when the wind torque and tidal torque become comparable, the trajectory follows Ω = Ω sta .It is worth noting that the curve at Ω = 0 almost coincides with Ω =  when [Fe/H] = -0.5 dex.Therefore, as long as the wind torque remains comparable to the tidal torque, the evolution will proceed at almost constant rotation frequency.This can explain why the purple and green solid lines can evolve with al-most constant stellar rotation and orbital frequency when they reach double synchronization until stability conditions are destroyed.The purple curve gradually crosses the equilibrium point at   c 0 where tidal torque becomes greater than wind torque, becoming unstable and eventually leading to the planet being consumed.When [Fe/H] = 0 dex and [Fe/H] = +0.5 dex, the stars undergo stellar contraction during the PMS phase, resulting in a gradual decrease in radius.This leads to a reduction in the stellar moment of inertia  * .According to equations 10 and 12, it can be understood that  c 0 decreases and  c 0 increases.Therefore, the evolutionary trajectory shows a trend towards the upper left direction.For the case of [Fe/H] = 0 dex (red and light blue solid lines), as evolution proceeds, the increasing tidal torque approaches the decreasing wind torque, and evolution will follow the trajectory Ω = Ω sta for a period of time, until the planet quickly falls into the star when   c 0 > 1.For the case of [Fe/H] = +0.5 dex (blue and yellow solid lines), 99% of the planetary life is spent in the saturated phase of the star, and the equilibrium point (blue dash-dotted line) almost coincides with Ω = 0. Due to the thicker convection envelope, magnetic braking is stronger, and the wind torque dominates, causing the star to rotate slowly.As evolution proceeds until the star is in an unsaturated state, the system evolves below the blue dash-dotted line, at which point the tidal torque dominates, causing the planet to be rapidly consumed and the star to spin up.
The , the planet gradually crosses the equilibrium point, and the tidal torque is greater than the wind torque.The evolution is below the red dash-dotted line, and the planet quickly is engulfed by the star.When 0 <  < 1 (purple and light blue solid lines), the planet migrates inward, and because the tidal torque is greater than the wind torque, the evolution proceeds almost horizontally until Ω = .The evolution then follows the trajectory of Ω = Ω sta and gradually crosses the equilibrium point as   c 0 eventually being engulfed by the star.For [Fe/H] = +0.5 dex, the evolution is the same as in the case of 1.1  ⊙ during the PMS phase.The stellar contraction will cause  c 0 to decrease and  c 0 to increase, and the evolution will proceed in the upper left direction.Until the tidal torque dominates (blue dash-dotted line below), the planet will quickly fall into the star.

The influence of planetary parameter and tidal quality factor on triggering the double synchronization
In this section, we discuss the impact of planet mass, tidal quality factor, and initial orbital semimajor axis on the triggering of the double synchronization mechanism in Figure 3. Firstly, from the top panel of Figure 3(a), it can be seen that only systems with  ini = 0.06 au,  ′ * = 10 6 , and  pl = 8.0  J (orange curve) can undergo the double synchronization mechanism and maintain it.Additionally, it can be observed from the figure that systems with an initial semimajor axis of a = 0.06 au can survive the main sequence stage.The models with an initial semi-major axis of a = 0.03 au survive during the main sequence stage only under the conditions of  ′ * = 10 8 and  pl = 1.0  J .The middle panel of Figure 3 reveals that systems with a planetary mass of 1.0  J have orbital angular momentum below the critical value, while systems with 8.0  J have orbital angular momentum above the critical value.As time progresses, for systems with an initial semi-major axis of  ini = 0.03 au and 8.0  J , the orbital angular momentum will become lower than the critical value, while for systems with  ini = 0.06 au and 8.0  J , the orbital angular momentum will still be above the critical value at the end of the main sequence.The bottom panel indicates that at the beginning of the evolution, only systems with  ini = 0.03 au and 8.0  J have total angular momentum below the critical value.However, as the evolution progresses until the end of the main sequence, only systems with  ini = 0.03 au will have total angular momentum below the critical value.
In the Darwin diagram of Figure 3(b), the equilibrium point trend can generally be divided into two categories: .As  ′ * increases and  pl decreases, the equilibrium point curve corresponding to Ω = 0 is closer to Ω = 0, indicating that the wind torque plays a more dominant role compared to the tidal torque..As  ′ * decreases and  pl increases, the equilibrium point curve corresponding to Ω = 0 is closer to Ω = , indicating that the tidal torque plays a more dominant role compared to the wind torque.
For the green, light blue, and blue curves in the figure 3, the wind torque dominates over the tidal torque, and their stellar rotation speed is faster than the orbital adjustment of the planet, resulting in nearly vertical evolution.The red curve in the figure 3 quickly drops below the red dash-dotted line as it evolves, indicating that the tidal torque dominates over the wind torque, and the planet is quickly engulfed by the star, resulting in nearly horizontal evolution.The purple and yellow curves quickly evolve past  c 0 , where the stronger tidal torque causes them to evolve in an almost horizontal manner.For the orange and brown curves in the figure 3, they are located at positions where  ≥ 1.The planet migrates outward for a period of time, and then the orange curve evolves along the trajectory of Ω = Ω sta until the end of the main sequence stage when it survives.The brown curve evolves towards the equilibrium point, but due to the planet's distance from the star and weaker tidal forces, the trajectory on the Darwin diagram is short and always above the equilibrium point dash-dotted line.

The correlation between double synchronization mechanism and initial parameters
We present in Figure 4 the initial parameter conditions for a starplanet system to sustain long-term double synchronization.It is noteworthy that the criterion for "long-term double synchronization" we define here is that the star-planet system has undergone this a significant double synchronization phase before the planets are engulfed (we assume that the star and planet maintain co-rotation for more than 10 Myr, then we consider the system to have undergone long-term double synchronization).Firstly, the density distributions of each parameter can be found from the histogram in the Figure 4. Through our grid calculations, we find that long-term double synchronization can only be sustained in systems with a star mass of 1.1 ⊙ or greater, an orbital semi-major axis of 0.02 au or greater, a planet mass of 4.0  J or greater, and a tidal quality factor  10 Q ′ * ≤ 8 (Please note that, due to the wide intervals in the parameters of planetary mass, we believe that systems with a planetary mass greater than 1.0  J are required for long-term double synchronization to occur.Additionally, it can be seen that we only find a double synchronization system if the tidal quality factor is equal to 10 8 , suggesting that long-term double synchronization is more common in systems with  ′ * ≤ 10 8 ).Furthermore, in addition to a peak occurring around an initial semimajor axis of 0.06 au, we can also see from the pairwise relationship diagrams of the initial parameters that larger stellar and planetary masses, as well as smaller values of  10 Q ′ * , [Fe/H], and the initial stellar rotation period, make it easier for double synchronous orbits to be maintained in the long term.This is because the low-metallicity stars with high mass have thinner convective envelopes that produce weaker angular momentum loss due to magnetic braking.This is reflected in the Darwin diagrams in Figures 1(b) and 2(b), where the equilibrium point corresponding to Ω = 0 is closer to the synchronous orbit Ω = n for Ω = Ω sta .On one hand, a weaker magnetic braking results in less stellar angular momentum loss, making it easier for the total angular momentum to remain above the critical value.The state of Ω ≈ n ≈ Ω sta is allowed for a long time because tidal friction between the star and planet can be ignored and the degree of orbit decay only depends on the weak magnetic braking.Moreover, larger planetary masses and smaller values of  10 Q ′ * lead to faster dissipation of tidal energy, which is similar to the previous case and is manifested by an equilibrium point at Ω = 0 that is closer to the synchronous orbit Ω =  in Figure 3(b).
When the initial stellar rotation period is less than the planetary orbital period, the planet's outward migration will inevitably evolve to the position of Ω = , implying a higher probability of long-term double synchronization.Additionally, under the same conditions, a faster initial stellar rotation period corresponds to a larger angular momentum, making it easier to achieve  >  c 0 .Regarding the initial semi-major axis, if it is too distant, the tidal torque is too weak.To achieve a balance between wind torque and tidal torque, a smaller stellar wind torque is required, resulting in fewer satisfying samples.If the semi-major axis is too close, tidal energy dissipation occurs too rapidly, triggering  >  c 0 .Moreover, a high stellar rotation rate is needed to synchronize with the orbital rate, which is evidently challenging.It is noteworthy that an intermediate orbital semi-major axis is more favorable for the long-term maintenance of tidal synchronization.As shown in Figures 1(b), 2(b), and 3(b), the system remains stable when the evolution is maintained within   c 0

DISCUSSION
It is highly significant to find the double synchronization system.Therefore, we chose 243 hot Jupiter systems and brown dwarf systems.Based on the results above, we primarily focus on close-in gas giant planetary systems with main sequence stars being F and G−type stars1 .The stellar mass ranges from 0.8  * to 1.5  * , and we consider systems with planetary masses above 0.1  J .The metallicity is within the range of +0.5 ≥ [Fe/H] ≥ −0.5, and the surface gravity is required to be logg ≥ 3.90 to ensure that the sample is in the main sequence phase.We select systems with orbital separation within 0.1 au.The parameters of these systems are provided in Table C1.For the majority of planetary systems in our sample, direct measurements of the stellar rotation period are not available.Therefore, we employ the Equation (5) from Brown (2014) to calculate the stellar rotation period, where the specific expression is given by  .blue dots represent samples that can maintain long-term double synchronization under these conditions, while blank areas represent samples that cannot maintain long-term double synchronization.The darker the color of the shaded area, the greater the number of samples at that location.The histogram in the figure represents the distribution of number densities for each initial parameter.There are a total of 425 samples in this part, which is less than 2% of our total sample size of nearly 25,000.It is important to note that the models we computed form a uniform grid as shown in Table 1.Each shaded region below is marked with blue dots, and the darkness of the shading indicates a larger number of samples in which the double synchronization state occurs.The deeper the color of the shaded region, the larger the range, representing more samples in which the double synchronization state is achieved.The addition of shaded regions is merely for a more intuitive display of the number of samples at different positions in the grid.The area between two grid points, where there is no shading, does not contain any samples.Regions with only blue dots and no shaded portions indicate very few samples at that position, possibly only one or a few.
rot = 2   *  sin  s sin  orb , the  sin  s represents the projection of stellar rotation, and  orb represents the inclination of the orbit. * is stellar radius.Assuming the system is aligned along the line of sight, the inclination of the stellar rotation axis with respect to the line of sight is denoted as  s =  orb .In order to calculate the convective turnover time( cz ), we referred to the work of See et al. (2021) and utilized the database from Amard et al. (2019) along with the maximum likelihood tool from Valle et al. (2014) to obtain  cz .When the stellar masses are higher than 1.5  * , their envelopes are radiative zones, and the influence of metallicity becomes negligible.We calculated the Rossby number  o =  rot /  , with the aim is to determine whether the star is in a saturated or unsaturated state.The distribution of these hot Jupiter systems is plotted in the Darwin diagram.
We included those samples with a stable solution regime (   c 0 and  >  c 0 ) and fast-rotating stars with Ω >  from Damiani & Lanza (2015).These samples are crucial for discussing whether a system is currently in or will experience a double synchronization state in the future.

The distribution of planetary systems in the Darwin diagram
In Figure 5, the locations of all the systems in the Darwin diagram, revealing two notable features.This result is consistent with the theoretical model predictions in Figure 4.In addition, the systems with large convective turnover time are generally closer to the position of Ω = 0 in the diagram.This can be explained by the fact that Gtype stars have stronger magnetic braking and slower rotation rates than F-type stars, resulting in lower total angular momentums.Under the same conditions, a larger planetary mass implies a larger orbital angular momentum.The figure also shows that all the planetary mass, stellar mass, and total angular momentum decrease sharply when the orbital velocity is high.This can be attributed to two aspects: on the one hand, when the planet is closer to the star, its orbital angular momentum is smaller than that of a more distant planet under the same conditions.On the other hand, the tidal forces in systems composed of large-mass planets and stars are stronger, leading to faster orbital decay and rapid engulfment of the planet.This makes it difficult to find such systems in close-in orbits.
Finally, for systems with / c 0 > 2, the planetary mass tends to be smaller.This can be attributed to two reasons: Firstly, Damiani & Lanza (2015) pointed out that for systems with / c 0 > 1.5, planets with a mass of about 3.0  J cannot maintain a stationary state, regardless of the magnetic braking efficiency of their host stars.This suggests that massive planets in close proximity find it difficult to remain stable, as tidal forces can lead to their rapid engulfment.Secondly, a larger planetary mass is usually associated with a larger critical frequency ( c 0 ), which to some extent favors smaller values of / c 0 for systems with larger planetary masses.

The Phases of long-term double synchronization Experienced by Planetary Systems
Through the analysis of the equilibrium states of planetary systems, we can infer whether they can maintain long-term double synchronization.As these systems evolve, the planets in a state of long-term double synchronization may eventually stay at this state, while systems currently non-synchronized could potentially enter into such a state with the evolution.
We show the relationship between Ω/, planetary mass and effective temperature in Figure 6(a) for our selected observational samples.It can be seen from the figure that Ω/ and  pl have a positive correlation, and the smaller the planetary mass and Ω/, the lower the effective temperatures of the stars.It can also be found that in the region where Ω/ ≥ 0.3, the effective temperatures are above 6000 K.This is because the stars in the low-temperature regions have stronger magnetic braking because of the thicker convective envelopes, and their rotation rates are often slower than those in the high-temperature regions.This leads to two consequences.First, there is a smaller Ω/ for stars in the low-temperature.Second, the planets with large mass can be quickly engulfed by their host star with lower effective temperatures due to a rapid transfer of the planetary angular momentum to star.This also makes it very difficult for us to find outward-migrating (Ω/ ≥ 1) samples in the low-temperature region.
We also show the relationship between stellar mass, metallicity and convective turnover time in Figure 6(b).We can clearly see that the convective turnover time decreases with the increase of stellar mass, and there is weak correlation between metallicity and convective turnover time.At a given stellar mass, the convective turnover times increase with increasing metallicity.
In Figure 7, we present the relationship between Ω/ and Ω sta /Ω for the 425 theoretical models that have experienced long-term double synchronization and 243 observational data with tidal quality factors  ′ * of 10 5 , 10 7 , and 10 9 .We find that on a logarithmic scale, those systems that have experienced a long-term double synchronization state are situated near the curve of Ω sta = .This means that for systems that have undergone long-term double synchronization states, the equilibrium state's Ω sta is very close to the orbital rate n for most of the time in life.The figure also shows that some discrete points with smaller  pl / cz values are below the curve of Ω sta = .Considering the distribution range of Ω sta /Ω, we assume points with (Ω sta − )/Ω ≥ 0.01 as discrete points, which is less than 0.5% of the total sample.Moreover, the discrete points with Ω/ < 1 occur when the planet is almost to be engulfed (in the last 1% of the entire evolutionary stage), where tidal torques dominate.The discrete points with Ω/ > 1 are at the early phase of the evolution when the planet starts to migrate outward (in the first 1% of the entire evolutionary stage).For points with (Ω sta − )/Ω ≤ 0.01, we consider as the situation of Ω sta ≈ , and this sample can be understood as falling on the curve Ω sta = .The figure illustrates that long-term double synchronization states require larger  pl / cz , which corresponds to shorter convective turnover times and larger planetary masses, consistent with the findings presented in Section 3. From the figure, we can also learn that even though these systems that have undergone long-term double synchronization states cannot guarantee to maintain this state throughout the entire main sequence of the star, most of them can still maintain Ω sta ≈  in the evolution process.This provides a new approach to discover long-term double synchronization systems.Using the condition of Ω sta = , we cannot only find systems that are currently in this state, but also the systems that have experienced this state in the past (Ω/ < 1), and systems that may experience this state in the future (Ω/ > 1).For example, it is evident that there are three systems with larger  pl / cz values that fall on the curve Ω sta =  regardless of the  ′ * value, though the planetary system's positions are closer to the curve of Ω sta =  if  ′ * is smaller.These three systems are CoRoT-3, HAT-P-57, and KELT-21.Among them, CoRoT-3 is near Ω/=1, and this system may currently be in a long-term double synchronization state.HAT-P-57 and KELT-21 are located in the area where Ω/ > 1, and as they evolve, they will maintain a long-term double synchronization state if these two systems can achieve orbital synchronization during the main sequence.
In the regions of Ω/ < 1, many points with  ′ * = 10 5 and some with  ′ * = 10 7 also fall on the curve.This indicates that if the initial state allows the system to experience orbital synchronization, then the synchronization state that occurred in the past can also be maintained for a long time.If we can accurately measure these  ′ * values in the future, for planets that fall on the curve Ω sta = , we can more accurately determine the stages of long-term double synchronization states experienced by the system.The vast majority of systems that fall below this curve can be excluded, even though many of them are in stable states and near synchronous orbits on the Darwin diagram.However, this synchronous orbit will be immediately broken as they evolve.
Orbital evolution can lead to changes in Ω sta , which raises the question of whether systems currently at Ω sta ≈  will evolve to Ω sta < , and if systems at Ω sta <  will evolve to reach Ω sta ≈ .To illustrate this issue, we searched for systems within an orbital period range of 0.2 to 10.0 days in which some maintain Ω sta <  while the rest are Ω sta = .We have identified five systems exhibiting this particular characteristic, which only manifests when the tidal quality factor  ′ * = 10 5 .Figure 8 shows the relationship between  orb and  sta , where HAT-P-56 and EPIC 246851721 are currently in the state of Ω/ > 1, indicating outward orbital migration.Interestingly, EPIC 246851721 is in a saturated state ( o < 0.14), and as  orb increases,  sta deviates further from Ω sta ≈ .However, as the stellar rotation period increases through evolution, the system will transit to an unsaturated state before orbital synchronization, where systems are at Ω sta < , this prevents long-term double synchronization.For HAT-P-56, in cases where orbital period  orb is greater than 8.0 days,  sta gradually moves away from Ω sta = .If synchronization occurs within  orb < 8 days, this long-term synchronized state will also be experienced.TOI-628, HATS-39, and HATS-41 are currently in a state where Ω/ < 1, and HATS-41's orbit is almost synchronized.However, HATS-41 is found to be at Ω sta < , indicating that this synchronization will not be maintained in the long term.The situations for TOI-628 and HATS-39 are more uncertain.If orbital synchronization occurred during a past evolutionary stage at Ω sta ≈ , the long-term double synchronization would be experienced.Otherwise, even if orbital synchronization has occurred, it will not be sustained in a long time.

CONLUSION
We used MESA for modeling and considered Matt et al. (2015) magnetic braking model based on Damiani & Lanza (2015).In this case, we were able to discuss in detail the relationship between magnetic braking and stellar mass and metallicity, and explore the different effects of stellar mass and metallicity on system stability.Due to the complexity of the model, we solved for the stellar rotation rate Ω sta at the balance between the wind torque and tidal torque for various parameter values.We extensively discussed the relationship between the stellar rotation period, planetary orbital period,  orb  orb,c 0 and   c 0 as a function of time for various initial parameters.We found that long-term double synchronization is more easily achieved for higher-mass stars, lower metallicity, smaller  ′ * , larger planetary mass, moderate orbital distance, and faster initial stellar rotation rate.In addition, we discussed the evolution of Ω sta as a function of   c 0 for various parameter values.We found that for lower metallicity, higher stellar mass, higher planetary mass, and smaller tidal quality factor  ′ * , the equilibrium curve of Ω sta is closer to Ω = n, and vice versa.When the equilibrium curve of Ω sta coincides with Ω = n, the system will maintain long-term double synchronization under the stable solutions proposed by Hut (1981)  We calculated uniform grids for nearly 25,000 models and found that 425 systems can undergo and maintain a long-term double synchronous phase.It is believed that in certain conditions, a star-planet system can maintain a double synchronous state for a long time, even if the planet is a poor metal star.Generally, when the planetary mass is larger and the tidal quality factor is smaller, tidal interactions tend to be stronger, and planet engulfment occurs faster.However, when the planetary orbit distance is moderate and the stellar rotation period is equal to the planetary orbit period, weaker magnetic braking can maintain double synchronization in stars with thinner convective envelopes.This mechanism may make these planets easier to survive.We conducted a quantitative analysis of double synchronization and found that: (i) For stars, systems that maintain double synchronization for a long time tend to occur when the star mass is 1.1  ⊙ or higher.And as their mass increases, the upper limit of metallicity is also higher.As for the initial rotation period of the star, if it is less than the orbital period of the planet, the exchange of angular momentum between them will be faster, and therefore a faster stellar rotation period also tends to maintain double synchronization more easily.
(ii) For planets, in most cases, when the initial semi-major axis  is around 0.06 au, the planet mass  pl ⩾ 4.0  J , and  10 Q ′ * ⩽ 7, double synchronization will be more easier to occur.The recent studies indicate that the logarithm of the tidal quality factor  10 Q ′ * for solar-like stars with orbital periods less than 15.0 days falls within the range of 5.7 to 7.5.(Patel et al. 2023) This means that when the convective layer of the host star is thin and the stellar rotation period is close to the planetary orbit period, the system is easier to be in a long-term double synchronization, making it easier for close-in and massive planets to survive during the main sequence of the star.
We surveyed a sample of 243 hot Jupiter and brown dwarf systems, with the primary stars being F-type and G-type stars that have relatively complete observational parameters.We calculated the convective turnover timescale,  cz , for each star and, in conjunction with theoretical models, found that for the majority of the main sequence stage, the equilibrium position Ω sta of long-term double synchronization systems is close to the current orbital rate n, regardless of whether the system is in a co-rotating state.We also discovered that in the regions near Ω/ ≥ 1, the systems CoRoT-3, HAT-P-57, and KELT-21 are positioned at Ω sta ≈  regardless of the value of  ′ * .Our work provides theoretical constraints for identifying observational samples in a double synchronization state, and we can find not only systems that may currently be in a long-term double synchronization state but also predict systems that have experienced or may experience this phenomenon in the past and future.As the tidal quality factor  ′ * is measured more precisely in the future, more and more systems that have gone through this phase will be discovered.
software:MESA R11554 (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019) )    and C as constants.The current accepted values for these constants are: A 14.713 degrees/day, B = -2.396degrees/day, C = -1.787degrees/day Finally, we utilized the aforementioned equation to calculate the normalized moments of inertia under two scenarios (sun rotation with differential rotation, using the equatorial rotation rate in the solid body model).The results are presented in the table A1.We find that the impact of differential rotation results in a difference of normalized moments of inertia of just under 5%, a negligible effect.gular momentum loss rate by an order of magnitude compared to G stars.In contrast, our model employs the Matt et al. (2015) magnetic braking model, which utilizes the convective turnover timescale to describe variations in stellar mass.Consequently, in the DL15 model, the angular momentum loss due to magnetic braking varies discontinuously with stellar mass, potentially resulting in a reversed trend in the calculated  sta .For systems with masses of 1.3  ⊙ and 10.0  J , both models yield  sta values equal to  orb when  orb is less than 5.0 days.However, when  orb exceeds 5 days, the  sta values obtained by the DL15 model gradually surpass  orb .Finally, the overall trend of our model is consistent with that of the DL15 model, namely, longterm double synchronization tends to occur in systems with larger planetary masses around F stars.However, for real observed systems, the magnetic braking scheme of the DL15 model cannot accurately provide the angular momentum loss rate for each star.Therefore, in the search for long-term double synchronous systems, employing a more sophisticated magnetic braking model such as Matt et al. (2015) can provide better constraints on stellar parameters.
Middle panel: the solid colored lines represents  orb  orb,c 0 , and the thin dashed black line represents the  orb =  orb,c 0 .In the bottom panel, the solid colored lines represent   c 0 , and the thin dashed black line represents  =  c 0 .In Figure 1 (b), we present the Darwin diagram.The solid colored lines represent the locus of the star-planet system for different metallicities and initial stellar rotation periods.The dashed and solid black lines represent the locus for Ω =  and Ω = 0, respectively.The colored dashdotted lines represent the locus when the rate of magnetic braking angular momentum loss is equal to the rate of planetary orbital transfer angular momentum to the star ( Ω = 0).The vertical solid black line represents the position of the Roche limit.Thin horizontal dashed black line and thin vertical dashed black line indicate the values  =  c 0 and   c Figure 2 is similar to Figure 1, and the only difference is that the stellar mass is 1.3  ⊙ .For a 1.3  ⊙ star, both planetary systems with [Fe/H] = -0.5 dex and [Fe/H] = 0 dex can trigger the mechanism of long-term double synchronization.Similar to the case of 1.1  ⊙ , in the system with an initial rotation period of 3.0 days and [Fe/H] = -0.5 dex, the planets can still survive until the end of the stellar main sequence.In the Darwin diagram of Figure 2(b), we can see that the Ω = 0 equilibrium points for [Fe/H] = -0.5 dex and [Fe/H] = 0 dex almost coincide with Ω = .When  > 1 (green and red solid lines), the planet migrates outward, and / c 0 decreases.When Ω = , green and red solid lines is located exactly at Ω = Ω sta in the Figure 2(b), and the torque from the wind is equal to the tidal torque.The evolution will follow the trajectory of Ω = Ω sta .For the case of [Fe/H] = 0 dex, when   c

Figure 1 .Figure 2 .
Figure 1.The panel (a) shows the relationship between age and various parameters, including stellar rotation period, orbital period, ratio of orbital angular momentum to critical orbital angular momentum, and ratio of total angular momentum to critical total angular momentum.The panel (b) shows the Darwin diagram of the evolution of planetary system angular momentum under the effects of tidal dissipation and magnetic braking.(The specific meaning of the curves in the figure can be found in the first paragraph of Section 3.1.)

Figure 3 .
Figure 3. Designations are the same as those in Figure 1.In panel (a), The star has a mass of 1.1  ⊙ and the metallicitiy [Fe/H] = -0.5 dex and the initial stellar rotation period is 3.0 days.In panel (b), the vertical solid black line and vertical dotted black line represent the position of the Roche limit of 1.0  J and 8.0  J , respectively.

Figure 4
Figure 4. blue dots represent samples that can maintain long-term double synchronization under these conditions, while blank areas represent samples that cannot maintain long-term double synchronization.The darker the color of the shaded area, the greater the number of samples at that location.The histogram in the figure represents the distribution of number densities for each initial parameter.There are a total of 425 samples in this part, which is less than 2% of our total sample size of nearly 25,000.It is important to note that the models we computed form a uniform grid as shown in Table1.Each shaded region below is marked with blue dots, and the darkness of the shading indicates a larger number of samples in which the double synchronization state occurs.The deeper the color of the shaded region, the larger the range, representing more samples in which the double synchronization state is achieved.The addition of shaded regions is merely for a more intuitive display of the number of samples at different positions in the grid.The area between two grid points, where there is no shading, does not contain any samples.Regions with only blue dots and no shaded portions indicate very few samples at that position, possibly only one or a few.

Figure 5 .Figure 6 .Figure 7 .
Figure 5.In the figure, the circle markers represent observed systems, where the color of the circle markers represents the convective turnover time of the stars, and the size of the circle markers represents the planetary mass.The black dashed line and solid line represent the trajectories for Ω =  and Ω = 0, respectively.

Figure 8 .
Figure 8.The figure shows the relationship between orbital period  orb and stellar rotation period corresponding to torque balance  sta for five observational systems: HAT-P-56, TOI-628, EPIC 246851721, HATS-39, and HATS-41.The black bold line indicates the position where  orb = sta , and all the colored solid lines represent the relationship between  sta and  orb for each system.The cyan dashed line is our hypothetical relationship for EPIC 246851721 when it is in an unsaturated state.The circle and pentagram symbols denote the current positions of the systems.

Figure A1 .
Figure A1.The figure depicts the relationship between the position from the center to the interior of the Sun and its density.The blue dots represent asteroseismic data provided byDziembowski et al. (1994), while the green solid line represents the polynomial fitting results we obtained.

Figure C1 .
Figure C1.Similar to Figure 8, we use dash-dotted lines to represent the results of our model and dashed lines to represent the results of the DL15 model.

Table 1 .
All possible initial parameter combination.

Table A1 .
The normalized moments of inertia

Table C1 :
This work considers the stellar and planetary parameters in the considered systems.The full table can be found online in a machine readable format.