On the competing effect of wind braking and interior coupling in the rotational evolution of solar-like stars

Solar-like stars (M<1.3 Msun) lose angular momentum through their magnetized winds. The resulting evolution of the surface rotation period, which can be directly measured photometrically, has the potential to provide an accurate indicator of stellar age, and is constrained by observations of rotation periods of coeval stars, such as members of Galactic open clusters. A prominent observational feature of the mass-rotation period diagrams of open clusters is a sequence of relatively slower rotators. The formation and persistence of this slow rotators sequence across several billion years imply an approximately coherent spin-down of the stars that belong to it. In particular, the sequence is observed to evolve coherently towards longer periods in progressively older clusters. Recent observations of the 700-Myr old Praesepe and the 1-Gyr old NGC 6811 clusters, however, seem to contradict this general pattern. While the 1 Msun stars on the slow rotators sequence of the older NGC 6811 have longer periods than their counterparts in the younger Praesepe, as expected, the two sequences essentially merge at lower masses (<0.8 Msun). In other words, low-mass stars seem to have not been spinning down at all in the intervening 300 Myr. Here we show that this behavior is a manifestation of the variable rotational coupling in solar-like stars. The resurfacing of angular momentum from the interior can temporarily compensate for the loss due to wind braking at the surface. The internal redistribution of angular momentum has a steep mass dependence; as a result, the re-coupling occurs at different ages for stars of different masses. Our model explains the stalled surface spin-down of low-mass stars between Praesepe and NGC 6811, and predicts that the same behavior should be observable at other ages in other mass ranges.


Introduction
Solar-like stars (M 1.3 M ⊙ ) have outer convection zones where magnetic fields are generated through dynamo action. These stars have therefore magnetically active atmospheres, and the braking torque exerted by their magnetized winds efficiently carries away angular momentum from their surfaces (e.g., Schatzman 1962;Kraft 1967;Weber & Davis 1967;Kawaler 1988). As a result, in contrast to their more massive counterparts, solar-like stars undergo a significant rotational evolution on the main sequence. For instance, the rotation period of a star of solar mass is observed to be P rot ≈ 1 day upon reaching the zero-age main sequence (40 Myr), to be contrasted with the 26.09 days of the present Sun (4.57 Gyr). The rotation period of solar-like stars, moreover, can be derived accurately from the analysis of their photometric light curves. The availability of a directly measurable quantity undergoing such large variations on the main sequence holds great potential as an accurate age indicator ("gyrochronology": Barnes 2003; Barnes et al. 2016b).
The most stringent constraints on the rotational evolution of solar-like stars come from P rot measurements for stars in Galactic open clusters, whose age is independently known from classical methods, such as isochrone fitting (e.g., Demarque & Larson 1964). The mass-rotation period diagram of a cluster, in analogy to his classical counterpart, the color-magnitude diagram, reveals the (rotational) evolutionary state of a sample of stars of different mass at a fixed age.
While stars in early pre-main sequence clusters (age ≈ 1-10 Myr, see, e.g., Moraux et al. 2013) have a broad range of periods with an approximately uniform distribution in mass, in older clusters a clear pattern gradually emerges. A distinct sequence of relatively slower rotators is clearly visible in ≈ 70 Myr-old clusters, or older (see Figure 1 of Barnes 2003).
At its outset, the slow rotators sequence is well-defined only at the high mass end of the solar-like regime (M ≈ 1.3 M ⊙ ), while at lower masses it coexists with a less structured broad distribution of faster rotators. The fast rotators gradually disappear in clusters of increasing age, as stars of lower and lower mass converge to the slow rotators sequence. By the age of 1 Gyr, all stars down to 0.6 M ⊙ have joined the slow rotators sequence. Although observations suggest that stars of mass M < 0.6 M ⊙ eventually reach the slow rotators sequence as well, data in sufficiently old cluster and for low masses are still too scarce to draw definitive conclusions (but see Newton et al. 2018, and references therein).
The details of the formation of the slow rotators sequence are still not well understood, and only semiempirical models have been proposed so far (Barnes 2010;Brown 2014;Gondoin 2017). The observational evidences summarized above seem to imply that stars settle down on the slow rotators sequence by means of a quick, one-time transition, with a mass-dependent timescale. Since stars of the same mass do not all converge on the slow rotators sequence at once, a third parameter (besides mass and age), and/or a stochastic transition from the fast-to the slowrotator regime must be involved.
Beside the open issue of its mechanism of formation, the persistence of the slow rotators sequence (clearly observed up to at least ≈ 4 Gyr in M 67: Barnes et al. 2016a) points to an intrinsic regularity in the rotational evolution of stars once they have converged on it. The first attempts to model the evolution of the slow rotators sequence assumed its mass dependence to be factorizable from the age dependence, and that the latter follows the classical spin-down relation proposed by Skumanich (1972), P rot ∝ √ t (Barnes 2003(Barnes , 2007. Although these assumptions were later found to be inadequate (Meibom et al. 2009(Meibom et al. , 2011Barnes & Kim 2010;Lanzafame & Spada 2015, hereafter LS15), the slow rotators sequences of clusters of different ages follow a monotonic relation in the mass-rotation period diagram, with older clusters corresponding to slower periods.
The most recent observations of the Praesepe cluster (700 Myr) and NGC 6811 clusters (1 Gyr), featuring an unprecedented sampling of the low-mass regime of their slow rotators sequences, seem to contradict this (so far) wellestablished fact Douglas et al. 2019). While ≈ 1 M ⊙ stars on the slow rotators sequence of these two clusters display the expected behavior (i.e. the younger Praesepe stars rotating faster than those in the older NGC 6811), the two sequences seem to overlap below 0.8 M ⊙ (see Figure 1). In other words, low-mass stars in NGC 6811 seem to have not been spinning down significantly in the intervening 300 Myr.
In this paper, we model the evolution of the slow rotators sequence, and in particular we provide a theoretical interpretation of the puzzling new data from Praesepe and NGC 6811. We show that our model, based on our previous work (LS15), can reproduce the apparent stalled spin-down observed in NGC 6811, and that this is the result of the redistribution of angular momentum from the stellar interior to the surface, which is temporarily able to offset the effect of the magnetic wind braking. The two key ingredients of the model are the wind-braking law and the massdependence of the rotational coupling timescale originally proposed by LS15. The latter is also in remarkable agreement with the results of Somers & Pinsonneault (2016), obtained from an independent analysis of a related, but distinct problem (the lithium depletion in solar-like stars).
Our results provide valuable constraints on the mechanisms that transport angular momentum in the interior of solar-like stars, whose physical nature is still uncertain and strongly debated (e.g., Charbonneau & MacGregor 1993;Ruediger & Kitchatinov 1996;Charbonnel & Talon 2005;Spada et al. 2010;Eggenberger et al. 2019  demonstrate that the assumption of a solid-body rotation profile in stellar interiors is inadequate to obtain accurate gyrochronology relations. This paper is organized as follows: we describe the observational data used to constrain our models in Section 2; we outline the physics of our rotational evolution model in Section 3; our results are presented in Section 4; we discuss our findings in Section 5; our conclusions are summarized in Section 6.
2. Observed evolution of the slow rotators sequence from 120 Myr to 1 Gyr We constrain our models with the most up-to-date observations available for three rich, well-studied clusters, namely, the Pleiades, Praesepe, and NGC 6811 (Rebull et al. 2016;Douglas et al. 2019;Curtis et al. 2019, respectively). In the mass-rotation period diagram of each of these clusters the slow rotators sequence is clearly recognizable between M ≈ 0.3 and 1.3 M ⊙ , encompassing the entire domain of solar-like stars (see Figure 1). The Pleiades is one of the youngest clusters in which the slow rotators sequence is observed; a large fraction of stars have not yet converged onto it. After approximately half a billion year, in the Praesepe cluster, most of the stars of mass 0.6 M ⊙ are on the slow rotators sequence. Finally, in NGC 6811, at the age of 1 Gyr, all stars with measured rotation period are on the sequence.
As can be seen in Figure 1, the slow rotators sequences of Praesepe and NGC 6811 overlap below M 0.8 M ⊙ . This behavior is at odds with a simple description of the slow rotators sequence as coherently evolving in the massrotation period diagram. Indeed, stars of mass ≈ 0.7 M ⊙ seem to undergo a "stalled spin-down" epoch between 700 Myr and 1 Gyr, while their more massive counterparts spin down as expected  3. A two-zone rotational evolution model with mass-depedent coupling

Description of the model
We apply the rotational evolution model formulated by LS15 to describe the evolution of the slow rotators sequence. The main assumption of the two-zone model is that the radiative zone and the convective envelope of the star are in a state of rigid rotation, with angular velocities Ω rad and Ω env , respectively (MacGregor & Brenner 1991). These two quantities specify completely the rotational state of the star at a given time; the two zones have therefore angular momenta J rad = I rad Ω rad , and J env = I env Ω env , respectively, I rad and I env being their moments of inertia. The physics included in the model is as follows: 1. Initial conditions. We assume initial rigid rotation, The initial period is assumed to be P 0 = 8 days for all masses; this choice is in qualitative agreement with the observed rotation period distribution of very young clusters (e.g., the 1-4 Myr-old Orion Nebula Cluster, Rebull 2001; see also Moraux et al. 2013). The interaction with the circumstellar disk is taken into account according to the disk-locking hypothesis (Koenigl 1991), i.e., assuming that Ω env remains constant for the duration of the disk lifetime. We assume a circumstellar disk lifetime of 5 Myr, independent of mass (Hernández et al. 2008). 2. Wind braking. The overall rotational evolution is driven by the angular momentum loss from the surface, due to the torque imposed by the magnetized stellar wind (e.g., Schatzman 1962;Weber & Davis 1967;Kawaler 1988). We adopt the wind braking law originally formulated by LS15, which incorporates the mass dependence proposed by Barnes & Kim (2010), and follows the classical ∝ Ω 3 env rotation rate dependence of Kawaler (1988): In the equation above, I star = I rad + I env is the moment of inertia of the whole star, and τ ov is the convective overturn timescale of the convection zone. The product I star τ ov encompasses the mass dependence of our wind braking law; both quantities in the product are normalized to their values for a 1 M ⊙ model, as indicated by the subscript "⊙"; K w is an overall calibration constant, whose scaling K 0 = 1.11 × 10 47 g cm 2 s depends on the choice of the units. 3. Internal angular momentum transport. As the envelope loses angular momentum, differential rotation develops, i.e., Ω rad (t) Ω env (t); exchange of angular momentum between the two zones can also occur. During the pre-main sequence phase, the radiative interior grows at the expense of the convection zone, which comprises the whole star during the Hayashi phase (e.g., Kippenhahn et al. 2012). Once on the main sequence, angular momentum redistribution between the radiative and the convective zone can be mediated by several processes, whose relative importance is still an open issue (see Bouvier et al. 2014, for a recent review). Our model accounts for this effect phenomenologically, by introducing a constant, mass-dependent timescale τ c , over which the excess of angular momentum of the interior, ∆J ≡ I env J rad − I rad J env I rad + I env , is transferred to the envelope. The rotational coupling timescale τ c is taken to be constant along the evolution, and to scale with the stellar mass as: This scaling, derived from semi-empirical fitting of LS15, was found by these authors to be robust to the choice of the wind braking law. It was also found to be in remarkably good agreement with the independent analysis of Somers & Pinsonneault (2016).
The two-zone model equations for the the evolution of Ω rad and Ω env are: (3) where the dot denotes differentiation with respect to time; R rad and M rad are the radius and the mass of the radiative core, respectively. All the stellar structure parameters (M rad , R rad , I rad , I env , τ ov ) are derived from stellar evolution models constructed with the YREC code in its non-rotational configuration (Demarque et al. 2008).
Our model contains five parameters in total: the initial conditions are set by the initial rotation period P 0 and by the circumstellar disk lifetime τ disk ; the wind braking law contains the calibration constant K w ; the rotational coupling timescale is specified by τ c,⊙ and α.

Comparison with LS15 and re-determination of the model parameters
LS15 presented a comparison of several wind braking laws and a statistical determination of the parameters of the twozone model (P 0 , τ disk , τ c ), based on the observational constraints on the slow rotators sequence available at the time.
The best-fitting values of the model parameters were obtained by means of a Markov Chain Monte Carlo (MCMC) procedure for several different stellar masses in the range M ≈ 0.7-1.1 M ⊙ (see LS15 for details). The main conclusions of their analysis can be summarized as follows. First, the wind braking law (1) (referred to as "KB" in LS15) captures the observed shape of the slow rotators sequence sufficiently well that it can be used to describe the rotational evolution of stars of different masses with a unique value of the calibration constant K w . Second, moving from solartowards low-mass stars, it is increasingly important to take into account internal differential rotation in order to accurately reproduce the observations; the best results were obtained with a rotational coupling timescale that scales with stellar mass according to equation (2). Building on our previous work, in this paper we adopt the two-zone model (3), and we retain the functional forms of the wind braking law (1) and of the mass-dependent coupling timescale (2). However, we re-determine the parameters to take into account the new data at lower masses available for Praesepe and NGC 6811. In particular, our goal is to ascertain whether our model can reproduce the stalled spin-down of low-mass stars observed in the slow rotators sequence of NGC 6811. The parameters governing the initial conditions, P 0 and τ disk , are determined by the requirement to reproduce the slow rotators sequence of the Pleiades, and have the same values as in LS15.
The parameters K w , τ c,⊙ , and α, on the other hand, control the subsequent evolution of the slow rotators sequence from ≈ 100 Myr onward. These parameters are strongly correlated with each other (see LS15), and must therefore be re-evaluated together. The new data on Praesepe and NGC 6811, extending significantly the mass range sampled, provide the opportunity to refine the estimate of these parameters, in particular the power law exponent α.
The revised best-fitting values of the parameters K w , τ c,⊙ , and α were obtained using a Python implementation of the Nelder & Mead (1965) simplex algorithm, which is available in the function optimize.minimize, part of the SciPy 1 package. Our revised estimates are: It should be emphasized that these parameters were not obtained by LS15 as a direct result of their MCMC procedure. Rather, K LS15 w = 4.5 was obtained as an average over the entire mass range considered (≈ 0.7-1.1 M ⊙ ), while τ c,⊙ , and α were determined a posteriori from the values of τ c fitted independently at each mass. Our current estimate of K w is compatible with the LS15 value within the uncertainties, and is in excellent agreement with their estimate for a 1 M ⊙ star (3.6, cf. Table 3 of LS15). The value of τ c,⊙ is nearly unchanged (τ LS15 c,⊙ = 25 Myr). Finally, the re-calibrated α is significantly different from α LS15 = 7.3. This difference is a consequence of the extension to a wider mass range, as well as the intrinsic ill-conditioning of the best-fitting procedure, due to the strong correlations among the parameters of the model.

Modeling the evolution of the slow rotator sequence
The results of our best-fitting two-zone model are compared with the observations in Figure 2. A non-parametric fit of the slow rotators sequences, constructed consistently with the procedure employed by LS15, is also shown for each cluster (black circles with error bars). Error bars of 1.5 days for the Pleaides, and of 1 day for Praesepe and NGC 6811, represent the intrinsic scatter of the sequences, which is astrophysical in origin.
The overall fit of the slow rotators sequences of the Pleiades, Praesepe, and NGC 6811 is quite good (Figure 2). In particular, our model reproduces the apparently stalled surface spin-down of the ≈ 0.7 M ⊙ stars between 700 Myr and 1 Gyr with respect to their more massive counterparts,  Figure 2); the Sun is also plotted for reference.
which results in the merging of the slow rotators sequences of Praesepe and NGC 6811.
The key ingredient in reproducing this behavior is the scaling of the rotational coupling timescale with the mass of the star according to equation (2). This is immediately apparent when we compare the predictions of our reference model with those obtained assuming a value of τ c independent of stellar mass (all the other parameters being the same).
As Figure 3 shows, the uniform-τ c assumption produces a slow rotators sequence with a steeper slope at the age of NGC 6811 in the mass range ≈ 0.5-0.85 M ⊙ . Moreover, the difference between the slow rotators sequences at 700 Myr and at 1 Gyr (bottom panel of the Figure) is monotonically increasing towards lower masses for the uniform-τ c model. We conclude that a model with τ c independent of mass cannot achieve a satisfactory fit of the slow rotators sequences of Praesepe and NGC 6811. Such a model, in particular, fails to reproduce the stalled spin-down of the ≈ 0.7 M ⊙ stars.
On the contrary, the two-zone model implementing the scaling of τ c according to equation (2) produces a separation between the slow rotators sequences at 700 Myr and 1 Gyr that is largest at ≈ 1 M ⊙ and decreases at lower masses, thus correctly reproducing the merging of the slow rotators sequences observed in Praesepe and NGC 6811. Figure 4 shows the evolution of the surface rotation period calculated with our rotational evolution model. This alternative view further illustrates the effect of the rotational coupling. Stars of mass 0.8 M ⊙ experience a significantly reduced surface spin-down than their more massive counterparts between 700 Myr and 1 Gyr. This effect is the direct consequence of the redistribution of angular momentum from the radiative zone to the convective envelope, which temporarily offsets the angular momentum loss at the surface.

Wind braking vs. rotational coupling
The rotational evolution according to the two-zone model equations (3) is controlled by the wind braking, described by equation (1), and by the rotational coupling, characterized by equation (2). Which of these two processes dominates over the other depends on the mass of the star as well as its age. We can gain some insight into the relative importance of these two effects by comparing their timescales. The wind braking law (1) has a steep dependence on the surface rotation rate (∝ Ω 3 env ). The efficiency of the angular momentum loss therefore varies by several orders of magnitude between the zero-age main sequence, when the rotation rate is maximum for all masses, and the mature main sequence phases ( 1 Gyr). We can characterize the wind braking efficiency by the timescale: The rotational coupling timescale, on the other hand, is assumed to be constant for the entire evolution of the star. The ratio τ wb /τ c is a measure of the relative importance of the angular momentum loss and the internal angular momentum transport. This ratio is plotted in Figure 5. When τ wb /τ c < 1, the angular momentum redistribution from the interior to the surface is not efficient enough to compensate for the angular momentum lost from the surface, and significant differential rotation develops. The region of the plot in Figure 5 where this condition is realized is shaded in gray. On the contrary, when τ wb /τ c > 1, the rotational coupling is able to maintain a quasi solid-body rotation profile inside the star, while its angular momentum evolution is driven by the wind braking at the surface.
For stars of solar mass or higher, the rotational coupling catches up with the wind braking already by the age of the Pleiades, and differential rotation has essentially no chance to develop at later times. Conversely, for stars of mass ≈ 0.75 M ⊙ the angular momentum transport from the interior is still less efficient than the wind braking up to ≈ 700 Myr. For a 0.5 M ⊙ star, the rotational coupling is still relatively inefficient even after several Gyr. The evolution of the surface rotation period of stars less massive than ≈ 0.75 M ⊙ will therefore be affected by the angular momentum resurfacing from the radiative interior between the age of Praesepe and that of NGC 6811.
From Figure 5 we can also deduce that the same effect should be observed in other mass ranges at other ages. For instance, our model predicts that this angular "momentum reservoir" effect is operational for stars of mass ≈ 0.8-1.0 M ⊙ at an age of ≈ 300 Myr, and essentially up to the age of the Sun for M 0.6 M ⊙ . We look forward to testing this prediction as rotation periods for cluster stars of appropriate age and mass become available in the future.

Discussion
We have presented a simple model for the rotational evolution of solar-like stars, which contains two main physical ingredients: the magnetized wind braking at the surface, and the angular momentum transport (or "coupling") between the interior and the surface.
The essential role of the rotational de-coupling and re-coupling between the radiative interior and the convective envelope in the pre-main sequence and early main sequence evolution of solar-like stars has been recognized since the works of Endal & Sofia (1981); Stauffer et al. (1984Stauffer et al. ( , 1985; Pinsonneault et al. (1989); Soderblom et al. (1993), just to name a few. This effect has been incorporated phenomenologically in all subsequent modeling efforts, and several candidate processes have been proposed (MHD waves/instabilities, gravity waves, see, e.g., Charbonneau & MacGregor 1993;Ruediger & Kitchatinov 1996;Spruit 2002;Talon & Charbonnel 2003;Charbonnel & Talon 2005;Spada et al. 2010;Brun et al. 2011;Oglethorpe & Garaud 2013). In spite of this, the physical nature of the processes redistributing angular momentum in the interior of stars is still a major open question.
Previous works (e.g., Denissenkov et al. 2010;Gallet & Bouvier 2015) have underlined the strong mass dependence of the core-envelope rotational coupling timescale. By concentrating on the rotational evolution of the slow rotators sequence (or "I-sequence" in the terminology of Barnes 2003), LS15 quantified the mass dependence of the rotational coupling timescale. The model proposed by LS15, which incorporates such dependence, can reproduce accurately the evolution of the slow rotators sequence. The new rotational data in NGC 6811 ) provide previously unavailable constraints, by extending the range of mass with known rotation periods at 1 Gyr down to M ≈ 0.6 M ⊙ . Together with the new Praesepe data , these prompted us to re-evaluate the parameters of the LS15 model.
The re-calibrated parameters are consistent with the results of LS15 within the uncertainties of their fit to the data available at that time. Our updated model also explains quantitatively the apparent halt in the spin-down of low-mass stars (M ≈ 0.7 M ⊙ ) between 700 Myr and 1 Gyr, recently discovered by Curtis et al. (2019). Indeed, this phenomenon is a manifestation of the angular momentum redistribution from the interior of the star to its outer envelope, which can temporarily compensate for the angular momentum lost from the surface via magnetic braking.
It should be emphasized that the observed lack of spindown arises naturally in our model from simple assumptions regarding the mass-dependence of the angular momentum coupling timescale. In other words, it is not necessary to postulate a phase of weakened wind braking regime to explain the evolution of the slow rotators sequence as observed in open clusters. In addition, our model provides quantitative predictions on the duration of this suspended spin-down phase, and on the epoch at which stars of different mass are expected to experience it. These predictions will be testable as new rotational data, sampling the slow rotators sequence at low stellar masses and in clusters of different ages, become available. Barnes (2003) originally defined the I-sequence observed in open clusters as comprising stars spinning down according to the phenomenological relation P rot ∝ √ t (Skumanich 1972). Our results show that such a working definition still retains its approximate validity even in the light of the most recent observations. Indeed, the overall shape of the slow rotators sequence is captured remarkably well by the wind braking law (1), which couples the mass-dependence proposed by Barnes & Kim (2010) with a dependence on Ω 3 env , which is well-known to reproduce the Skumanich law (see, e.g., Kawaler 1988).
Departures from the Skumanich law in the slow rotators sequence have been know for a long time (e.g., Meibom et al. 2009Meibom et al. , 2011, and were already interpreted by LS15 as arising from the re-coupling of the surface with the interior. This effect manifests a strong dependence of the coupling timescale on stellar mass (cf. equation 2). Building on our previous results, we show that the halt in the spindown of low-mass stars reported by Curtis et al. (2019) is a consequence of the re-coupling occurring at the age of 1 Gyr for ≈ 0.7 M ⊙ stars.
The new observations of Curtis et al. (2019) allow us to refine the determination of the parameters in equation (2), to extend its validity to lower masses. This places valuable constraints on the unknown physical nature of the processes that transport angular momentum in the interior of solar-like main sequence stars. In particular, we confirm (within the uncertainties) the steep massdependence found by LS15, and independently recovered by Somers & Pinsonneault (2016). A successful theoretical description of the rotational coupling in solar-like stars from first principles should be able to explain the semi-empirical result of equation (2).
The deviations from a purely Skumanich spin-down that depend on the stellar mass are of practical importance for gyrochronology. Indeed, the original gyrochronological relations, based on the assumption of a purely factorable color and age dependence (where the color is used as a proxy of stellar mass), P rot = f (B − V ) g(t), were already shown to be untenable by Barnes & Kim (2010). Our model is a promising step towards the construction of more accurate gyrochronology relations, extracted from a simple physical model. This improvement comes at the price of not having a simple all-encompassing formula that can be readily inverted to derive the age from the rotation period. For convenience of use, we provide updated rotational isochrones for a wide range of ages in Appendix A.
Our model reproduces satisfactorily the slow rotators sequence as observed in the Pleiades cluster (≈ 100 Myr), and its subsequent evolution, as constrained by the recent observations of Praesepe and NGC 6811 (700 Myr and 1 Gyr, respectively), in the entire domain of solar-like stars (≈ 0.4-1.3 M ⊙ ). A two-zone model like ours, however, is obviously not applicable without modifications to fully convective stars ( 0.35 M ⊙ ). Whether such low-mass stars converge on the slow rotators sequence, and whether their evolution is similar to that of their more massive counterparts in spite of this qualitative difference in their interior structure, are still open questions.

Conclusions
We have presented a model for the rotational evolution of solar-like stars that reproduces the features of the slow rotators sequence as observed in open clusters between ≈ 100 and 1000 Myr. In particular, our model captures satisfactorily its mass-dependence in the range 0.4-1.3 M ⊙ , and the apparent stalled spin-down observed in ≈ 0.7 M ⊙ stars between 700 Myr and 1 Gyr.
The key ingredients of the model are the scalings of the wind-braking law and of the rotational coupling timescale with stellar mass. The former is well-represented by the product of the moment of inertia of the star times its convective overturn timescale; the latter follows a steep powerlaw (of exponent ≈ −5.6), with a coupling timescale ≈ 22 Myr for a 1 M ⊙ star.
Our results are a promising step towards more physically motivated gyrochronology relations, and highlight the necessity to take into account the internal transport of angular momentum in modeling the rotational evolution of solar-like stars. In addition, our model provides constraints on the currently unknown processes that transport angular momentum in the interior of solar-like stars.

Appendix A: Rotational isochrones
Selected rotational isochrones calculated with our updated two-zone model are reported in Table A.1, and plotted in Figure A.1. The Table lists the surface rotation period, in days, as a function of stellar mass and age; (B−V ) colors are also given, as calculated from a solar metallicity, 400 Myr isochrone from the YaPSI database (Spada et al. 2017).
In Figure A.2 we compare our current isochrones with those of LS15. The updated isochrones are consistent with our previous results in the range of overlap. A moderate disagreement is visible at ages 2 Gyr, reflecting the different choice of the overall calibration constant K w in equation (1)