Turbulent viscosity acting on the equilibrium tidal flow in convective stars
DDraft version July 28, 2020
Typeset using L A TEX twocolumn style in AASTeX63
Turbulent viscosity acting on the equilibrium tidal flow in convective stars
J´er´emie Vidal and Adrian J. Barker Department of Applied Mathematics, University of Leeds, Leeds, LS2 9JT, UK
Published 2020 January 16 in ApJLABSTRACTConvection is thought to act as a turbulent viscosity in damping tidal flows and in driving spinand orbital evolution in close convective binary systems. This turbulent viscosity should be reduced,compared to mixing-length predictions, when the forcing (tidal) frequency | ω t | exceeds the turnoverfrequency ω cv of the dominant convective eddies. However, two contradictory scaling laws have beenproposed and this issue remains highly disputed. To revisit this controversy, we conduct the first directnumerical simulations (DNS) of convection interacting with the equilibrium tidal flow in an idealizedglobal model of a low-mass star. We present direct computations of the turbulent effective viscosity, ν E , acting on the equilibrium tidal flow. We unexpectedly report the coexistence of the two disputedscaling laws, which reconciles previous theoretical (and numerical) findings. We recover the universalquadratic scaling ν E ∝ ( | ω t | /ω cv ) − in the high-frequency regime | ω t | /ω cv (cid:29)
1. Our results alsosupport the linear scaling ν E ∝ ( | ω t | /ω cv ) − in an intermediate regime with 1 ≤ | ω t | /ω cv (cid:46) O (10).Both regimes may be relevant to explain the observed properties of close binaries, including spinsynchronization of solar-type stars and the circularization of low-mass stars. The robustness of thesetwo regimes of tidal dissipation, and the transition between them, should be explored further in morerealistic models. A better understanding of the interaction between convection and tidal flows is indeedessential to correctly interpret observations of close binary stars and short-period planetary orbits. Keywords: binaries: close — convection — hydrodynamics — planet-star interactions — turbulence INTRODUCTIONTidal interactions determine the orbital and spin evo-lution of short-period planets and low-mass binary stars(e.g. Mazeh 2008; Zahn 2008; Ogilvie 2014). A majorweakness of tidal theory is in modeling how tidal flowsinteract with convection. Turbulent convection is be-lieved to act as an effective turbulent viscosity in damp-ing large-scale tidal flows (e.g. Zahn 1966). This mecha-nism is usually invoked to explain the circularization andsynchronization of binary systems containing low-massor solar-like main-sequence stars (e.g. Zahn 1989; Mei-bom & Mathieu 2005; Meibom et al. 2006; Van Eylenet al. 2016; Lurie et al. 2017; Triaud et al. 2017; vonBoetticher et al. 2019), and evolved stars (e.g. Verbunt &Phinney 1995; Beck et al. 2018; Price-Whelan & Good-man 2018).
Corresponding author: J´er´emie [email protected]
The turbulent viscosity, ν E , is usually estimated byneglecting the oscillatory nature of the tidal flow. Thisleads to ν E (cid:39) ν cv , where ν cv is predicted by mixing-length theory (MLT; e.g. Spiegel 1971). However, whenthe tidal frequency | ω t | is faster than the turnover fre-quency ω cv of the dominant convective eddies, ν E oughtto be reduced, as recognized initially by Zahn (1966).The magnitude of this reduction remains highly dis-puted. Two scaling laws that are based on phenomeno-logical arguments have been proposed, with either a lin-ear reduction ν E ∝ ν cv ( | ω t | /ω cv ) − (Zahn 1966, 1989),or a quadratic suppression ν E ∝ ν cv ( | ω t | /ω cv ) − (Gol-dreich & Keeley 1977; Goldreich & Nicholson 1977). Re-visiting this controversy has been attempted recently byusing direct numerical simulations (DNS). The two lawshave only been recovered in separate studies, which sup-port either the linear scaling (Penev et al. 2007, 2009b,a)or the quadratic one (Ogilvie & Lesur 2012; Braviner2016; Duguid et al. 2020). Thus, any application oftidal theory to stars (or planets) with convection zones a r X i v : . [ a s t r o - ph . S R ] J u l J. Vidal and A. J. Barker remains uncertain. Resolving this issue is essential be-fore we can apply tidal theory to interpret observationsof close binaries (e.g. Kirk et al. 2016; Lurie et al. 2017;Van Eylen et al. 2016; Triaud et al. 2017; Price-Whelan& Goodman 2018) or short-period planetary orbits (e.g.Rasio et al. 1996). For instance, circularization of sub-giant stars with orbital periods of approximately oneday could occur in either ∼ or 10 yr, depending onwhich scaling is valid (Price-Whelan & Goodman 2018).Owing to the importance of this problem to interpretobservations, we revisit this controversy using global nu-merical simulations. So far, only numerical studies us-ing local models and with simplified imposed shear flowshave been undertaken. Local DNSs may not capture thefull complexity of the tidal response existing in a globalmodel. They also could be affected by the adoptedboundary conditions. We therefore set out to gain inde-pendent physical insight from global DNSs of convectionin the presence of more realistic tidal flows. This Letteris organized as follows. We present our global modelin Section 2. Direct computations of the turbulent vis-cosity are presented in Section 3. The implications andastrophysical extrapolation of our results are presentedin Section 4, and we conclude the Letter in Section 5. DESCRIPTION OF THE TIDAL PROBLEMWe study the interplay between tides and turbulentconvection in a global model of a low-mass star (or core-less giant planet). The primary body is a full sphere ofradius R , filled with a fluid of uniform (laminar) kine-matic viscosity ν and thermal diffusivity κ . This body issubjected to tidal forcing from an orbiting companion.We model convection in the Boussinesq approximation,studying slight departures from a motionless conductionstate sustained by homogeneous internal heating. Sincemany low-mass stars are slow rotators (e.g. Nielsen et al.2013; Newton et al. 2018), and for simplicity, we neglectrotation in this study. We define the temperature per-turbation Θ and the velocity field u + U , where wesplit the flow into a background large-scale tidal flow U and a perturbation u . We use dimensionless unitsfor the simulations, adopting R as our length scale and R /ν as our timescale. Convection is then governed bytwo dimensionless numbers, the Rayleigh number Ra (which measures the strength of the convective driving)and the Prandtl number P r = ν/κ . We solve the sys-tem of equations in their weak variational form by usingthe spectral-element code Nek5000 (Fischer et al. 2007),employed previously for tidal studies (e.g. Favier et al.2014; Barker 2016; Reddy et al. 2018). Further detailsof the model are given in Appendix A. l, m − − E ( u ) − / l ≥ m ≥ Figure 1.
Simulation of uniformly heated turbulent con-vection without tides ( Ra = 10 , P r = 1). Top : three-dimensional snapshot of the velocity field component u x . Bottom : spectra of the instantaneous kinetic energy E ( u ),as a function of the spherical harmonic degree l ≥ m ≥ r ∈ [0 . , . Previous numerical studies modeled the tidal flowwith either an (ad hoc) external forcing (Penev et al.2009a) or with a background unidirectional shear flowin a shearing box (Ogilvie & Lesur 2012; Braviner 2016;Duguid et al. 2020). Here, we instead consider self-consistently the large-scale (non-wavelike) equilibrium urbulent viscosity acting on the equilibrium tidal flow in convective stars l = 2 and azimuthal order m = 2(Ogilvie 2014). In the inertial frame, the resulting (di-mensionless) flow is in the xy -plane and takes the form U = − ω t β (cid:32) sin( ω t t ) cos( ω t t )cos( ω t t ) − sin( ω t t ) (cid:33) (cid:32) xy (cid:33) , (1)where β (cid:28) ω t is the tidal (angular) frequency (twicethe orbital frequency in the absence of rotation). Theforcing amplitude β must be large enough to obtain ameasurable tidal response in the presence of convection,but large values could strongly modify the convection.The global simulations that we present here are verydemanding, because they must be run for a sufficientlylong duration to reduce turbulent noise. This inevitablyrestricts our survey of parameter space. We simu-late highly supercritical convection with Ra = 10 and P r = 1, which can be compared with the value for lin-ear onset Ra ≥ β = 0)is illustrated in Figure 1. The kinetic energy is char-acterized by a nonnegligible axisymmetric component(consistent with the flow in the top panel), and a shortinertial-like range illustrated by the Kolmogorov scaling( ∝ − /
3, bottom panel). We quantitatively estimatethe convective turnover frequency as ω cv = u rms /l E ,where u rms is the time average of the volume-averagedrms velocity and l E (cid:39) / ω cv (cid:39) . ± . TURBULENT VISCOSITYWe determine numerically the effective viscosity coef-ficient ν E , which is the leading-order component of theeffective viscosity tensor at the forcing frequency (e.g.Penev et al. 2009a). This is obtained numerically bybalancing the mean rate at which convection does workon the tidal flow with the mean rate of viscous dissipa-tion of this flow (e.g. Goodman & Oh 1997; Braviner2016). This leads to ν E = − ω t β ) ∆ T (cid:90) Tt (cid:104) u · [( u · ∇ ) U ] (cid:105) V d t, (2) | ω t | /ω cv − − − | ν E | ν E ’ l E u rms | ν E | ∝ ( | ω t | /ω cv ) − | ν E | ∝ ( | ω t | /ω cv ) − | ω t | /ω cv − − − | ν E | ν E ’ l E u rms | ν E | ∝ ( | ω t | /ω cv ) − | ν E | ∝ ( | ω t | /ω cv ) − Figure 2.
Direct measurements of the effective viscosity ν E in turbulent convection ( Ra = 10 , P r = 1), as a function of | ω t | /ω cv . Red squares: ν E >
0. Blue circles: ν E <
0. Errorbars are conservatively defined using two standard deviationsfrom the mean value.
Top : β = 10 − . Bottom : β = 5 × − .Horizontal dashed lines: expected behavior from MLT ν E (cid:39) l E u rms in the low-frequency regime ( | ω t | (cid:28) ω cv ). with (cid:104) · (cid:105) V = (1 /V ) (cid:82) V · d V the volume average and 5 ≤ ∆ T = T − t ≤
10, with t an appropriate initial timein the saturated regime. We have verified in AppendixB that the spatial average is not dominated by regionsnear the boundary, and is instead due to interactionswith turbulent flows in the bulk.Results for the effective viscosity ν E are shown in Fig-ure 2, for the tidal amplitudes β = 10 − and β = 5 × − . The former value is similar to that for a solar-massbinary in a one-day orbit. The effective viscosity de- J. Vidal and A. J. Barker creases as we increase | ω t | /ω cv . The striking feature hereis the coexistence of both heuristic scaling laws. First,we obtain an intermediate regime 1 ≤ | ω t | /ω cv (cid:46) O (10),in which ν E is consistent with the linear reduction (Zahn1966, 1989). This trend is clearer in the simulations withthe largest tidal amplitude ( β = 5 × − ), because thesignal-to-noise ratio (S/N) is lower for weaker tides (asshown by the error bars in Figure 2). Second, we clearlyobtain the quadratic law | ν E | ∝ ( | ω t | /ω cv ) − in the high-frequency regime | ω t | (cid:29) ω cv (Goldreich & Nicholson1977; Goodman & Oh 1997). The transition betweenthese two scalings is sharp, occurring when | ω t | /ω cv (cid:39) β = 5 × − (bottom panel), but appears to dependweakly on the tidal amplitude. These results demon-strate that both scaling laws are obtained in our globalmodel, which have only been found previously in sepa-rate studies in Cartesian geometry.In the low-frequency regime ( | ω t | (cid:46) ω cv ), we havebeen unable to accurately determine ν E . The ampli-tude of tidal flow (1) was too weak to give a suffi-ciently strong S/N. A crude extrapolation of our resultsis broadly consistent with MLT, which would predict ν E ∝ ν cv (cid:39) l E u rms when | ω t | → DISCUSSION4.1.
Non-Kolmogorov Turbulence?
Turbulent viscosity is often defined with a closuremodel that relates the Reynolds stress to the rate ofstrain. Ogilvie & Lesur (2012) and Duguid et al.(2020) demonstrated the viscoelastic nature of the high-frequency ( | ω t | /ω cv (cid:29)
1) tidal response, developing anasymptotic theory for the Reynolds stress. This stronglysupports the quadratic reduction for ν E . In our globalsimulations, we have confirmed the viscoelastic charac-ter of the response for high-frequency tidal forcing (seeAppendix B). However, this asymptotic theory does notstrictly apply for lower frequencies. Indeed, a linear re-duction may result from the non-Kolmogorov nature ofthe turbulence (Penev et al. 2007, 2009b,a).We illustrate in Figure 3 the frequency spectrum of theconvection. The largely non-Kolmogorov nature of theconvection is revealed by the frequency spectrum of thevolume-averaged Reynolds stress component (cid:104) u x u y (cid:105) V .The latter quantity, which is directly related to the ef-fective viscosity (see Appendix B), has a shallower decaywith frequency (in f − ) than expected from Kolmogorovtheory when | ω t | /ω cv ≤ O (10). This slope is largelyunaffected by the tidal flow, and so is a generic prop- erty of the convection in this range. The slope of thenon-Kolmorogov spectrum is similar to that reportedin Penev et al. (2007, 2009a), despite the model dif-ferences. For larger frequencies, a steeper decay is ob-served, first behaving like f − in apparent agreementwith local simulations of Rayleigh-Benard convection(Kumar & Verma 2018), and then rapidly decaying (cor-responding with a dissipation range). The transition be-tween the linear and quadratic reductions may broadlycoincide with where the shallow non-Kolmogorov scal-ing in f − ceases to be valid. Then, even though ourspectrum is still non-Kolmogorov-like, a quadratic re-duction is found for higher frequencies, in agreementwith prior asymptotic theory (Ogilvie & Lesur 2012;Duguid et al. 2020). The frequency spectrum of thethermal energy (cid:104) Θ (cid:105) V / of ν E . However, our simulationsdo not currently allow us to assess their arguments con-clusively. In summary, our new global simulations sup-port both the linear and quadratic reductions for theeddy viscosity.4.2. Astrophysical Implications
We extrapolate our findings to stellar interiors as fol-lows. MLT predicts the rms convective velocity to scaleas u rms ∝ ( Ra/P r ) / in the fully turbulent regime (e.g.Spiegel 1971), such that ν E /ν (cid:39) ν cv /ν ∝ ( Ra/P r ) / isindependent of tidal frequency when | ω t | /ω cv (cid:28)
1. Sucha frequency-independent ν E is consistent with constanttidal lag-time models (e.g. Hut 1981). Then, the effec-tive viscosity is reduced in the presence of fast tides,first with an approximately linear reduction and then aquadratic one. The transition between these two regimesoccurs when | ω t | /ω cv (cid:39) O (10). Further work is requiredto explore the robustness of the transition when Ra/P r is increased. We have also obtained statistically signifi-cant negative values of ν E for high frequencies in Figure2, which is consistent with previous local results (Ogilvie& Lesur 2012; Duguid et al. 2020). Negative valuesprobably result from (necessarily) adopting simulationparameters that are far removed from their astrophys-ical values. This phenomenon is always observed when | ν E | ≤ ν (here in the quadratic regime). As also foundin Duguid et al. (2020), the negative values occur whenwhen | ω t | lies in the dissipation range of the turbulence(see Figure 3). MLT predicts that ω cv ∝ Ra / in thefully turbulent regime, and that the inertial-like range As suggested by the referee, based on Phinney (1992). urbulent viscosity acting on the equilibrium tidal flow in convective stars f − − − | F { h u x u y i V } | ∝ f − ∝ f − f − − − − | F { h u x u y i V } | ∝ f − ∝ f − Forcing f − − − − − | F { h Θ i V } | / ∝ f − ∝ f − f − − − − − | F { h Θ i V } | / ∝ f − ∝ f − Forcing
Figure 3.
Frequency spectra of turbulent convection ( Ra = 10 , P r = 1). Spectrum of |F{(cid:104) u x u y (cid:105) V }| (top panel) and |F{(cid:104) Θ (cid:105) V / }| (bottom panel), where F is the Fourier transform and f is the (ordinary) frequency. Gray area shows theintermediate frequency range 1 ≤ | ω t | /ω cv ≤ Left : unperturbed convection.
Right : perturbed convection with β = 5 × − and | ω t | /ω cv = 14 . should extend to higher frequencies. Typical values forthe Rayleigh and Prandtl numbers in solar-like stars are Ra = 10 − and P r = 10 − − − (Hanasoge &Sreenivasan 2014), such that we expect ν cv (cid:29) ν . Thus,unrealistically large values of | ω t | may be required to getnegative values ν E ≤ τ cv = 2 π/ω cv (e.g. Terquem et al. 1998), using an estimate based onthe stellar luminosity (e.g. Price-Whelan & Goodman2018) τ cv ≈ .
37 yr ( M cv /M (cid:12) ) / ( T e / − / , (3)with M cv the mass of the convective zone, M (cid:12) the so-lar mass and T e the effective temperature (in Kelvin).Assuming a transition between the two regimes when J. Vidal and A. J. Barker
Figure 4.
Different regimes for the frequency reduction ofturbulent viscosity for solar-type binaries, with an assumedtransition at | ω t | /ω cv = 5. Red hashed zone: ν E (cid:39) ν cv . Grayzone: linear reduction. Green zone: quadratic reduction.Binaries, sorted by eccentricity e , extracted from Figure 7in Lurie et al. (2017). Top : solar-like star (1 M (cid:12) ) at 1 Gyrwith τ cv (cid:39) . M cv ≈ . M (cid:12) . Bottom : low-massstar (0 . M (cid:12) ) at 1 Gyr with τ cv (cid:39) .
25 yr and M cv ≈ . M (cid:12) . | ω t | /ω cv (cid:39) P t = 2 π/ | ω t | (cid:46) . − . τ cv ∼ . − .
25 yr, where P t = | /P s − /P orb | − / P s therotation period and P orb the orbital period (in nonsyn-chronized systems). Since low-mass stars typically havelonger timescales τ cv , the transition occurs for larger or-bital periods for these objects (see Figure 4).The range of validity of the various turbulent vis-cosity prescriptions is shown in Figure 4, for binariesgiven in Lurie et al. (2017). Both scalings are shownto be relevant for this sample. Therefore, equilibriumtide theory must be carefully applied to interpret theobservational data. Convective damping of the equi-librium tide could potentially explain the main fea-tures of this distribution. The timescale for tidal spin-synchronization, for a solar-mass binary in a circular or-bit with P orb = 10 days and P s = 15 days, is estimatedto be approximately 1 Gyr if we adopt a continuous pro-file for ν E ( | ω t | ) to fit Figure 2 (see Appendix C). Thisseems to be an efficient mechanism for P orb (cid:46)
10 d.We have also superimposed our theoretical predictionsfor spin-synchronization timescales (due to convectivedamping of the equilibrium tide), using two differentstellar models (see Appendix C) that span the majorityof the sample in Lurie et al. (2017). They suggest thatthe quadratic reduction could explain, for larger valuesof | ω t | , why some short-period binaries in Figure 4 havenot yet synchronized. Finally, the mechanism seems tooefficient to explain why some systems with P orb < P orb < P s are not synchronized. This may bedue to the young ages or high masses of these stars, orperhaps because they are affected by differential rota-tion (Lurie et al. 2017). The dynamical tide may alsobe important for some of these systems (e.g. Ogilvie &Lin 2007; Ogilvie 2014). CONCLUDING REMARKSIn this Letter, we have revisited the long-standingproblem of the interaction between tidal flows and tur-bulent convection. We have conducted the first numeri-cal simulations of turbulent convection within an ideal-ized global model of low-mass fully convective stars (orcore-less giant planets), to measure the turbulent viscos-ity acting on the large-scale equilibrium tidal flow. Wehave reconciled, for the first time and within a singleconsistent physical model, the two contradictory scalinglaws that have been proposed to describe the frequencyreduction of the effective viscosity when the tidal fre-quency exceeds the dominant convective turnover fre-quency (i.e. fast tides, Zahn 1966; Goldreich & Nichol-son 1977). Our results have confirmed the universality urbulent viscosity acting on the equilibrium tidal flow in convective stars ν E ∝ | ω t | − when | ω t | /ω cv (cid:29) ν E ∝ | ω t | − ), inan intermediate regime 1 (cid:46) | ω t | /ω cv (cid:46) O (10). Thislikely results from the non-Kolmogorov nature of theturbulence in that frequency range (e.g. Penev et al.2007, 2009b,a). This has important consequences forinterpreting astrophysical observations. Our findingsshould guide future data-driven studies to discriminatebetween these two scaling laws, for instance, when in-terpreting observations of the synchronization and cir-cularization of main-sequence binaries (Lurie et al. 2017;Triaud et al. 2017) or the circularization of evolved stars(Price-Whelan & Goodman 2018).Much further work is required before we can accu-rately model the tidal evolution of astrophysical systemsdue to this mechanism. The robustness and coexistenceof these two scaling regimes, and the transition betweenthem, should be explored further. Moreover, we haveneglected dynamical tides (e.g. Ogilvie & Lin 2007) andconsidered only circular orbits. Different tidal compo-nents could, however, be damped at different rates (e.g.Lai 2012). Simulations in spherical shells would be alsoworth exploring (e.g. Gastine et al. 2016) to model theconvective envelopes of solar-like stars. Given the im-portance of this problem, understanding the interactionbetween turbulent convection and tidal flows appearsurgent. This is necessary to correctly interpret obser-vations of close binaries (e.g. Lurie et al. 2017; Triaudet al. 2017; Price-Whelan & Goodman 2018).The validity of MLT should be also assessed usingturbulent simulations of convection. MLT predictionscould underestimate the turbulent viscosity ν cv in thelow-frequency regime (Goldman 2008). Indeed, depar-tures from MLT have been found in recent simulations ofcompressible convection (e.g. Anders et al. 2019). More- over, convection-driven turbulence is strongly affectedby rapid rotation (e.g. Gastine et al. 2016; Kaplan et al.2017), such as in giant planets or young stars. Theprescription for the turbulent viscosity from MLT (e.g.Barker et al. 2014) then ought to be modified (see Mathiset al. 2016, in the low-frequency regime). More realis-tic convection models should be considered as a long-term endeavor. Finally, by neglecting rotation, we havealso filtered out nonlinear tidal flows such as the ellip-tical (tidal) instability (e.g. Barker et al. 2016; Vidal &C´ebron 2017). They could enhance tidal dissipation forthe shortest orbital periods (Barker 2016; Vidal et al.2018, 2019), and might even modify properties of tur-bulent convection (C´ebron et al. 2010). Understandingtheir interplay with convection deserves future work.ACKNOWLEDGMENTSThis work was funded by STFC grant ST/R00059X/1.JV is grateful to S. Reddy for his support withNek5000. Numerical simulations were undertakenon ARC1/ARC2 clusters (HPC Facility, Universityof Leeds). Simulations were also performed on theDiRAC Data Intensive service at Leicester (STFCDiRAC HPC Facility). The equipment was fundedby BEIS capital funding via STFC capital grantsST/K000373/1 and ST/R002363/1 and STFC DiRACOperations grant ST/R001014/1. This work alsoused the DiRAC@Durham facility. The equipmentwas funded by BEIS capital funding via STFCcapital grants ST/P002293/1, ST/R002371/1 andST/S002502/1, Durham University and STFC opera-tions grant ST/R000832/1. DiRAC is part of the Na-tional e-Infrastructure. We acknowledge the anonymousreferee for several suggestions that have allowed us toimprove the paper. Software:
Nek5000 (https://nek5000.mcs.anl.gov/).APPENDIX A. CONVECTION MODELWe study Boussinesq thermal convection (e.g. Spiegel1971), driven by homogeneous internal heating Q T in afull sphere. We use the notation introduced in Section2. The gravitational field is g = − g r , where r is theposition vector and g is a constant. This is the leading-order component for a low-mass body that is not verycentrally condensed (e.g. Lai et al. 1993). We employ di-mensionless quantities for the simulations, adopting R as the length scale, the viscous timescale R /ν as thetimescale, and ( ν Q T R ) / (6 κ ) as the unit of temper- ature (as in Monville et al. 2019). The dimensionlessequations for u and the temperature perturbation Θ, inthe inertial frame, are ∂ u ∂t + ( u · ∇ ) u = −∇ p + ∇ u + Ra Θ r − f , (A1a) ∂ Θ ∂t + ( u · ∇ ) Θ = 1 P r (cid:2) u · r + ∇ Θ (cid:3) − Q , (A1b)where p is a dimensionless (reduced) pressure, f =( u · ∇ ) U + ( U · ∇ ) u a forcing term with U givenby (1) and Q = ( U · ∇ ) Θ. We have defined theRayleigh number Ra = α T g Q T R / (6 νκ ), where α T is J. Vidal and A. J. Barker the thermal expansion coefficient, and the Prandtl num-ber
P r = ν/κ . The nonlinear term ( U · ∇ ) U reduceshere to a pressure gradient, and thus plays no dynami-cal role within the Boussinesq approximation. We havealso neglected in Q the term ( U · ∇ ) T that shouldvanish in the limit β (cid:28) T is the background tem-perature. Equations (A1) are complemented with theincompressibility condition ∇· u = 0, and boundary con-ditions at the (dimensionless) spherical boundary r = 1.For the temperature, we employ the isothermal condi-tion Θ = 0. To avoid spurious numerical issues asso-ciated with angular momentum conservation in globalsimulations of tidal flows (as explained in Guermondet al. 2013; Favier et al. 2014), we enforce the mechanicalboundary condition u = . This is unlikely to affect the(small-scale) turbulent flows driven in the bulk withoutrotation (compared to stress-free boundary conditions).We have solved nonlinear Equations (A1) in theirweak variational form by using the spectral-element codeNek5000 (e.g. Fischer et al. 2007). The computationaldomain is decomposed into E = 3584 non-overlappinghexahedral elements. Within each element, the veloc-ity (and pressure) is represented as Lagrange polyno-mials of order N (respectively, N −
2) on the Gauss-Lobatto-Legendre (Gauss-Legendre) points. Temporaldiscretization is accomplished by a third-order method,based on an adaptive and semi-implicit scheme in whichthe nonlinear and Coriolis terms are treated explicitly,and the remaining linear terms are treated implicitly.Solutions are de-aliased following the 3 / N/ N points are used for the lin-ear terms. We have checked the numerical accuracy intargeted simulations to ensure convergence by varyingthe polynomial order from N = 7 to N = 9. We adopta time step 10 − ≤ d t ≤ × − (in dimensionlessunits, depending on the forcing frequency).For most of the simulations, we initiated the convec-tion with random noise to the temperature field and letit saturate without tides (i.e. β = 0), before switch-ing on the equilibrium tidal flow. We have checked thatinitiating the convection together with the equilibriumtidal flow does not lead to noticeably different results.We have integrated each simulation for several viscoustimescales (5 ≤ ∆ T ≤
10 in dimensionless units), corre-sponding with more than a hundred tidal periods, to ob-tain converged statistics for the effective viscosity. Thetime average in expression (2) is obtained by fitting alinear slope to the cumulative time integral (e.g. seeFigure 13 in Duguid et al. 2020), to reduce the turbu-lent noise. B. COMPLEMENTARY RESULTSThe parameters and results of the simulations behindFigure 2 are given in Table 1. We define the rms veloc-ity u rms as the time average of (2 (cid:104)E ( u ) (cid:105) V / / with E ( u ) = ( u x + u y + u z ) / l E (cid:39) /
3, which agrees with Figure 1.The latter figure indeed shows that multiple eddies spanthe radius of the body. Then, we define the turnover fre-quency as ω cv = u rms /l E . Small differences in the rmsproperties of the convection are found when the ampli-tude of the tidal flow was larger than the convectiveflow (i.e. when β | ω t | ≥ u rms ). Typically, these differ-ences are smaller than 5% for the kinetic energy and therms velocity when β ≤ × − . However, in the strongtides regime the convection can be modified more signif-icantly, which we have observed when β ≥ − , or forvery high frequencies (i.e. | ω t | /ω cv (cid:29)
100 at Ra = 10 with β = 5 × − ). Similar findings have been reportedin local simulations (e.g. Duguid et al. 2020).We illustrate in Figure 5 spatial spectra of the term u · [ u · ∇ U ] that appears in equation (2) for the effec-tive viscosity. We have computed it in the entire fluiddomain (i.e. 0 . ≤ r ≤ .
99) and omitting the bound-ary regions (i.e. 0 . ≤ r ≤ . l = 0 component in the physi-cal space) is never dominated by interactions near theboundary, but is instead due to flows in the bulk.For an alternative approach to directly using equa-tion (2), we can estimate ν E by considering the Fouriertransform of the volume-averaged component (cid:104) u x u y (cid:105) V ,as long as we account for the oscillatory nature of thetidal flow (Ogilvie & Lesur 2012). Thus, we can definethe effective viscosity (cid:98) ν E in the Fourier domain as F{(cid:104) u x u y (cid:105) V } = (cid:98) ν E ω t β F{ cos( ω t t ) } , (B2)where F denotes the Fourier transform and (cid:98) ν E is acomplex-valued quantity. The Reynolds stress and therate of strain are generally out of phase. The real part (cid:60) e ( (cid:98) ν E ) (cid:39) ν E is the turbulent viscosity (that is in phasewith the equilibrium tidal flow), whereas the imaginarypart (cid:61) m ( (cid:98) ν E ) (that is out of phase with the tidal flow) isrelated to an effective elasticity. In the regime of high-frequency tidal forcing ( | ω t | (cid:29) ω cv ), Ogilvie & Lesur(2012) and Duguid et al. (2020) demonstrated the vis-coelastic nature of the tidal response. They predictthe turbulent viscosity (cid:60) e ( (cid:98) ν E ) should scale as | ω t | − ,whereas (cid:61) m ( (cid:98) ν E ) should obey a linear reduction | ω t | − .We compute (cid:98) ν E from expression (B2) in Figure 6.The results confirm the universal nature of the vis-coelastic response, with a dominant elastic component urbulent viscosity acting on the equilibrium tidal flow in convective stars Table 1.
Table of Simulation Results for Ra = 10 and P r = 1. | ω t | (cid:104)E ( u ) (cid:105) V u rms | ω t | /ω cv ν E UpEb LwEb2 . × +2 . × +3 . × +1 . × +0 +1 . × +1 . × +1 . × +0 . × +2 . × +3 . × +1 . × +0 +1 . × +1 . × +1 . × +0 . × +2 . × +3 . × +1 . × +0 +1 . × +1 . × +0 . × +0 . × +2 . × +3 . × +1 . × +0 +1 . × +1 . × +0 . × +0 . × +2 . × +3 . × +1 . × +0 +1 . × +1 . × +0 . × +0 . × +2 . × +3 . × +1 . × +0 +7 . × +0 . × +0 . × +0 . × +2 . × +3 . × +1 . × +0 +3 . × +0 . × +0 . × +0 . × +2 . × +3 . × +1 . × +0 +2 . × +0 . × − . × − . × +2 . × +3 . × +1 . × +0 +1 . × +0 . × +0 . × − . × +2 . × +3 . × +1 . × +0 +9 . × − . × − . × − . × +3 . × +3 . × +1 . × +0 +4 . × − . × − . × − . × +3 . × +3 . × +1 . × +0 − . × − . × − . × − . × +3 . × +3 . × +1 . × +1 − . × − . × − . × − . × +3 . × +3 . × +1 . × +1 − . × − . × − . × − . × +3 . × +3 . × +1 . × +1 − . × − . × − . × − . × +4 . × +3 . × +1 . × +1 − . × − . × − . × − . × +4 . × +3 . × +1 . × +2 − . × − . × − . × − | ω t | (cid:104)E ( u ) (cid:105) V u rms | ω t | /ω cv ν E UpEb LwEb2 . × +2 . × +3 . × +1 . × +0 +1 . × +1 . × +0 . × +0 . × +2 . × +3 . × +1 . × +0 +1 . × +1 . × +0 . × +0 . × +2 . × +3 . × +1 . × +0 +1 . × +1 . × +0 . × +0 . × +2 . × +3 . × +1 . × +0 +1 . × +1 . × +0 . × +0 . × +2 . × +3 . × +1 . × +0 +8 . × +0 . × +0 . × − . × +2 . × +3 . × +1 . × +0 +7 . × +0 . × − . × +0 . × +2 . × +3 . × +1 . × +0 +6 . × +0 . × − . × +0 . × +2 . × +3 . × +1 . × +0 +4 . × +0 . × +0 . × +0 . × +2 . × +3 . × +1 . × +0 +3 . × +0 . × +0 . × +0 . × +2 . × +3 . × +1 . × +0 +3 . × +0 . × − . × − . × +3 . × +3 . × +1 . × +0 +5 . × − . × − . × − . × +3 . × +3 . × +1 . × +0 +2 . × − . × − . × − . × +3 . × +3 . × +1 . × +0 − . × − . × − . × − . × +3 . × +3 . × +1 . × +1 − . × − . × − . × − . × +3 . × +3 . × +1 . × +1 − . × − . × − . × − . × +3 . × +3 . × +1 . × +1 − . × − . × − . × − . × +4 . × +3 . × +1 . × +2 − . × − . × − . × − Note — (cid:104)E ( u ) (cid:105) V is the volume-averaged kinetic energy and u rms the rms velocity. UpEb: Upper Error bar. LwEb: LowerError bar. Top : β = 10 − . Bottom : β = 5 × − . at high frequencies. Indeed, we broadly obtain a lin-ear reduction |(cid:61) m ( (cid:98) ν E ) | ∝ | ω t | − in the high-frequencyregime. Moreover, we recover the expected scaling in | ω t | − for the turbulent viscosity |(cid:60) e ( (cid:98) ν ) E | in the high-frequency regime (in addition to a linear reduction factorin an intermediate regime), which is always smaller than |(cid:61) m ( (cid:98) ν E ) | . The effective viscosity has approximately thesame amplitude when it is calculated using (B2) or equa-tion (2), so this cross-validates our computations for the turbulent viscosity. We have checked that quantitativelysimilar results are obtained by considering the othercomponents (cid:104) u x (cid:105) V and (cid:104) u y (cid:105) V of the Reynolds stress ten-sor. This agrees with Penev et al. (2009a), who showedthat the effects of convective turbulence on a large-scaleoscillatory shear flow is well represented by an effectiveviscosity coefficient.0 J. Vidal and A. J. Barker l − − − − − − / . ≤ r ≤ . . ≤ r ≤ . . ≤ r ≤ . l − − − − − − / . ≤ r ≤ . . ≤ r ≤ . . ≤ r ≤ . Figure 5.
Power spectrum of the time-averaged and radiallyintegrated quantity − / (∆ T ( ω t β ) ) (cid:82) u · [ u · ∇ U ] d t , asa function of the spherical harmonic degree l ≥ Ra = 10 , P r =1 and β = 5 × − . Eddy viscosity (2) is given by thesquare root of the l = 0 component (in the physical space),when ∆ T is large enough to reduce the turbulent noise. Top : | ω t | /ω cv = 2 . T (cid:39) . Bottom : | ω t | /ω cv = 13 with an integration time ∆ T (cid:39) .
39. Spectrahave been computed by interpolating the data to a sphericalgrid. C.
ESTIMATES FOR TIDAL SYNCHRONIZATIONWe provide details here to compute the tidal dissi-pation timescales shown in Figure 4. For a circular andaligned orbit, the turbulent dissipation is estimated froma spherical stellar model as (e.g. see Equation (85) inRemus et al. 2012) D ν = 4 π R GM | ω t | (cid:90) α R x R ρ ∗ ν E d x R , (C3)where x R = r/R is the normalized radius, M is thestellar mass, α R is the ratio of the radius of the baseof the convective envelope to the stellar radius R , and ρ ∗ ∼ townsend/static.php?ref=ez-web) for a1 solar-mass star at 1 Gyr (assuming the metallicity Z =0 . Q (cid:48) is related to D ν by Q (cid:48) = 3 / (2 D ν ). The resulting timescale for tidalsynchronization of the stellar spin is then τ Ω = 13 πr g (cid:18) M + M M (cid:19) P orb P dyn P s D ν , (C4)where M is the mass of the companion, P orb is theorbital period, P s is the stellar rotation period, P dyn =2 π/ (cid:112) GM/R is the dynamical timescale, and r g ≈ . u cv ( x R ), the mixing length l E ( x R ) = 2 H p ( x R )with H p ( x R ) the pressure scale height, and the convec-tive frequency ω cv ( x R ) = u cv ( x R ) /l E ( x R ). Then, weassume that the frequency reduction of the turbulentviscosity ν E ( x R ) obeys the continuous profile (Figure 7) ν E = u cv (cid:96) E | ω t | /ω cv < ,ω cv / | ω t | ( | ω t | /ω cv ∈ [1 , , ω cv / | ω t | ) ( | ω t | /ω cv > , (C5)where the proportionality constant is arbitrary but ischosen here to be consistent with Figure 2. The ap-parent discontinuity, reported in Figure 2, apparentlycoincides with the rapid passage through zero of ν E inthe simulations. Since negative values of ν E may not berelevant in reality in the frequency range | ω t | /ω cv ≤ P orb = 10 days and P s ∼
15 days (for example), we would obtain τ Ω ≈ . τ Ω ≈ ν E . The urbulent viscosity acting on the equilibrium tidal flow in convective stars | ω t | /ω cv − − − | < e ( b ν E ) | ∝ ( | ω t | /ω cv ) − ∝ ( | ω t | /ω cv ) − | ω t | /ω cv − − | = m ( b ν E ) | ∝ ( | ω t | /ω cv ) − | ω t | /ω cv − − − | < e ( b ν E ) | ∝ ( | ω t | /ω cv ) − ∝ ( | ω t | /ω cv ) − | ω t | /ω cv − − | = m ( b ν E ) | ∝ ( | ω t | /ω cv ) − Figure 6.
Real and imaginary parts of the effective viscosity |(cid:60) e ( (cid:98) ν E ) | and |(cid:61) m ( (cid:98) ν E ) | , measured from the Reynolds stress tensor,as a function of | ω t | /ω cv in supercritical simulations of convection ( Ra = 10 , P r = 1). Error bars are computed by evaluatingthe noise level in the vicinity of the spike at the forcing frequency. The turbulent viscosity (relevant for tidal dissipation) ismeasured by |(cid:60) e ( (cid:98) ν E ) | , whereas |(cid:61) m ( (cid:98) ν E ) | measures the elastic component of the response. Top : β = 10 − . Bottom : β = 5 × − . resulting synchronization timescales for profile (C5) ob-tained using Equation (C4) are superimposed in Figure4. Our extrapolation indicates that convective damping of the equilibrium tide can be important in driving spinsynchronization in the sample presented in Lurie et al.(2017).REFERENCES Anders, E. H., Lecoanet, D., & Brown, B. P. 2019, ApJ,884, 65Barker, A. J. 2016, MNRAS, 459, 939Barker, A. J., Braviner, H. J., & Ogilvie, G. I. 2016,MNRAS, 459, 924Barker, A. J., Dempsey, A. M., & Lithwick, Y. 2014, ApJ,791, 13 Beck, P. G., Mathis, S., Gallet, F., et al. 2018, MNRAS,479, L123Braviner, H. J. 2016, PhD thesis, University of CambridgeC´ebron, D., Maubert, P., & Le Bars, M. 2010, GeoJI, 182,1311Duguid, C. D., Barker, A. J., & Jones, C. A. 2020,MNRAS, 491, 923 J. Vidal and A. J. Barker
Figure 7.