The tidal excitation of r modes in a solar type star orbited by a giant planet companion and the effect on orbital evolution II: The effect of tides in the misaligned case

We extend the study of Papaloizou&Savonije of the tidal interactions of close orbiting giant planets with a central solar type star to the situation where the spin axis of the central star and the orbital angular momentum are misaligned. We determine the tidal response taking into account the possibility of the excitation of r modes and the effect of tidal forcing due to potential perturbations which have zero frequency in a non rotating frame. Although there is near resonance with r modes with degree l' = 1 and orders m = 1 or -1 , half widths turn out to be sufficiently narrow so that in practice dissipation rates are found to be similar to those produced by non resonant potential perturbations. We use our results to determine the evolution of the misalignment for the full range of initial inclination angles taking account of the spin down of the central star due to magnetic braking. Overall we find the rate of tidal evolution to be unimportant for a one Jupiter mass planet with orbital period of about 3.7d over a main sequence lifetime. However, it becomes significant for higher mass planets and shorter orbital periods, approximately scaling as the square of the planet mass and the inverse fourth power of the orbital period.


INTRODUCTION
We extend the study of the tidal interaction of a solar mass primary with a Jupiter mass secondary in a close circular orbit carried out by Papaloizou & Savonije (2023) (hereafter PS) to consider the case when the stellar spin and orbital angular momentum vectors are misaligned.As in PS it is assumed that turbulent viscosity (see eg. Zahn 1977;Duguid et. al. 2020) operates in the stellar convective envelope while the spin angular momentum of the planet is neglected.
A number of physical processes have been invoked in order to produce the initial distribution of alignment angles of close orbiting giant planets.These include quiescent phenomena such as disc migration that would lead to close alignment (eg Lin & Papaloizou 1986), and dynamical interactions with a range of strengths (see eg Siegel et al. 2023, and references therein) that may produce significant misalgnment.Subsequent to formation, tidal interactions can result in orbital evolution leading to synchronisation of the orbital and component spins.(see Ogilvie 2014, for a review) and so affect the observed distribution.Hence it is important to understand the extent to which this may have occurred in order to understand conditions just post formation.
In PS the spin and orbital angular momenta were assumed to be aligned.In particular the spectrum of r modes associated with spherical harmonics of degree l ′ = 3 and l ′ = 5 was studied in detail.It was found that tidal interaction is unlikely to lead to significant orbital evolution away from these r mode resonances which were found to be very narrow.However, it may become significant if resonance can be maintained during tidal evolution.This may be possible as a result of the central star being spun down through a process such as magnetic breaking.This counters the tendency of tidal interaction to spin up the star enabling the resonance to be maintained.Systems where this may have operated are Kepler 1643 and COROT 4 ( see discussion in PS and references therein).
In this paper we extend the calculations of PS to consider tidal interactions for which the spin and orbital angular momenta are misaligned.The aim is to use the results to consider the effects of tidal interactions on the spin-orbit alignment of close orbiting giant planets.This involves determining the response to tidal forcing associated with spherical harmonics of degree l ′ = 2 and azimuthal mode numbers m = ±1, as well as for m = ±2 which occurs in the aligned case.
When there is misalignment, and these forms of perturbing tidal potential are viewed in a frame that is aligned with and corotating with the central star, forcing frequencies Ω s for |m| = 1, and 2Ω s for |m| = 2 occur, with 2π/Ω s being its rotation period respectively.The associated disturbances are stationary in a non rotating frame and may be associated with a strong inertial mode response in the convective envelope (eg Papaloizou & Pringle 1981;Ogilvie & Lin 2007;Ivanov & Papaloizou 2010;Ogilvie 2014) as well as an r mode response in the radiative core (PS).
These features mean that the tidal dissipation in the misaligned case may not be simply related to that in the aligned case, an aspect we are able to investigate.
As in PS we consider a simplified model for the rotating primary star valid up to first order in Ω s .In this approximation Coriolis forces are retained but centrifugal forces neglected such that the configuration is spherically symmetric.This means that centrifugal distortion which would lead to the slow precession of the spin angular momentum vector around the total angular momentum vector is neglected.The approximation can be viewed as having introduced a background potential to cancel out the centrifugal potential which remains fixed when the model is perturbed.This situation results in the model not possessing the rigid tilt mode that exists without the approximation.
The rigid tilt mode has l ′ = 1, |m| = 1, and eigenfrequency zero in an inertial frame, and though related to spin precession is not associated with tidal dissipation.Notably both the approximate and full model are expected to have a spectrum of higher order r modes with eigenfrequencies, which when viewed in a frame rotating with the star, have magnitude → Ω s with a relative correction that → 0 as Ω s → 0 (Papaloizou & Pringle 1978;Dewberry 2023).Here we find that this spectrum does not play a dominant role on account of the narrowness of resonance widths (see Section 6).We use our results to investigate the effect of the tidal interaction in combination with magnetic breaking of the central star on the evolution of the spin-orbit alignment over the main sequence life time.
The plan of this paper is as follows.In Section 1.1 we describe the basic configuration adopted for the system consisting of a primary star and orbiting planet.The coordinate systems defining the orbit and stellar frames together with notation used is given in 1.2.The development of the perturbing tidal potential acting on the primary in terms of spherical harmonics defined in the two coordinate systems with the aid of Wigner matrices is described in Sections 2-2.3 We go on to consider the perturbation to the external gravitational potential due to the primary and use it to determine the tidal forces acting on the orbit, and equations governing its evolution, by making use of the Gauss equations in Sections 2.4-2.8.Discussion of the equations governing the evolution of the semi-major axis and the angle between the spin and orbital angular momenta, β, together with estimates of evolutionary time scales are provided in Sections 3-3.2.1.
The Numerical calculation of the tidal response of the star and the viscous dissipation in its convective envelope is outlined in Section 4, with numerical results presented in Section 5.These include calculations of the fundamental r mode resonance and the general tidal response for a range of primary rotation rates in Sections 5.1 -5.2.Numerical calculations of the evolution of the spin orbit alignment of the system incorporating the effects of magnetic braking are considered in Sections 5. 3-5.4The effects of increasing planet mass and/or decreasing the orbital period, which speed up the tidal evolution rate, are described in Section 5.5.Finally we discuss our results and conclude in Section 6.

Basic configuration
Following PS we consider a giant planet secondary of mass M p orbiting a star of mass M * taken to be a solar mass described as the primary ⋆ .Here, in an extension of the work of PS, we consider the situation where the orbital and spin angular momenta are misaligned.Though, as in PS, we assume efficient orbital circularisation potentially induced by tides assumed to act on the secondary, thus limiting consideration to circular orbits.

Coordinate system and notation
Following Ivanov & Papaloizou (2021) (herafter IP) we define three Cartesian coordinate systems, each of them with origin O at the centre of mass of the primary.The (X, ′ Y ′ , Z ′ ) system, or orbit frame, is such that the orbital angular momentum, L, defines the direction of the Z ′ axis.The X ′ and Y ′ axes lie in the orthogonal orbital plane.The (X, Y, Z) system, or stellar frame, has Z axis pointing in the direction of the stellar spin angular momentum S. We then choose the Y axis to lie in the (X ′ , Y ′ ) plane with the angle between Y and Y ′ being γ.† The angle between L and S is β.When the total angular momentum J = L + S is conserved, it is natural to consider the (X ′′ Y ′′ Z ′′ ) system, which we call the primary frame, such that Z ′′ points in the direction of J.
Naturally OZ ′ , OZ, and OZ ′′ lie in the same plane.The angle between OZ ′ and OZ ′′ is i.
The angle between the line of apsides associated with the assumed near Keplerian orbit and the X ′ axis is denoted by ϖ.Note that then the angle between the apsidal line and the Y axis, which can be chosen to be the line of nodes ‡ is given by ϖ + γ − π/2.The coordinate systems together with the angles γ, β and ϖ are illustrated in Fig. 1 of IP.
⋆ As the secondary is treated as a point mass and the tidal response is linear, results can be scaled to apply to different masses (PS) † This is defined in the sense of a positive rotation about OZ ′ from OY to OY ′ .‡ This is because OZ ′ and OZ ′′ can be generated from OZ by rotations about OY and are thus normal to it.

THE PERTURBING TIDAL POTENTIAL
The perturbing potential, U , can be expressed in spherical polar coordinates (r, θ ′ , ϕ ′ ) defined in the orbit frame (X ′ , Y ′ , Z ′ ) with origin at the centre of mass of primary(see IP and PS).As usual r is the distance to the origin, θ ′ is the angle between the radius vector and the Z ′ axis and ϕ ′ is the azimuthal angle.We have up to O(a −3 ) where a is the semi-major axis of the circular orbit that lies in the plane θ ′ = π/2, P 2 is the usual Legendre polynomial, Y 2,m (θ ′ , ϕ ′ ), denotes the usual spherical harmonic of degree l ′ = 2 and The azimuthal angle of the line joining the components is Φ.Both Φ and ϕ ′ are measured from from the X ′ axis.
For the circular orbit we set Φ = n o t + ϖ, where n o is the mean motion.The reference angle, ϖ, measured from the X ′ axis can be taken to be the longitude of the apsidal line once the orbit is slightly perturbed.Following SP we allow ϖ to vary with time on a time scale long compared to n −1 0 as the orbit, being influenced by tides is slightly non Keplerian.Thus the angle Φ increases at a mean rate given by dΦ dt where the enclosure in the angled brackets denotes a time average.Accordingly 2.1 Expressing the perturbing potential in terms of spherical harmonics defined in the stellar frame We use Wigner matrices (see e.g.Khersonskii, Moskalev & Varshalovich 1988) to express spherical harmonics defined in the orbit frame as a linear combination of spherical harmonics expressed in coordinates defined in the stellar frame.These are defined in the same way as for the orbit frame but with reference to the (X, Y, Z) system.Thus where the coefficients or Wigner matrix elements D (2) n,m (0, β, γ), being expressed using standard notation, depend on the angles β and γ defined in Section 1.2, which specify the magnitudes of the angles of rotation required to transform the (X, Y, Z) system to the (X ′ , Y ′ , Z ′ ) system.First this requires a rotation through an angle γ about the Z axis, followed by a rotation through an angle β about the original Y axis (see Section 1.2 and IP).Also we may write where d (2) n,m is an element of Wigner's (small) d-matrix and is real (see e.g.Khersonskii, Moskalev & Varshalovich 1988, and IP).

The inverse transformation
The inverse transformation to (4) expresses spherical harmonics in the stellar frame in terms of those in the orbit frame.The (X ′ , Y ′ , Z ′ ) system converts to to the (X, Y, Z) system through a rotation −β about the Y ′ axis followed by a rotation through an angle −γ about the Z ′ axis.Thus the transformation corresponding to the inverse of ( 4) is given by Corresponding to (5) we have The inverse transformation (6) also follows from the big D and small d matrices being unitary.
2.3 The perturbing tidal potential in the stellar frame Applying the Wigner transformation (4) with (5) the perturbing tidal potential (3) is then defined with coordinates defined in the stellar frame.One obtains a sum of Fourier modes with time dependence through a factor, exp(−imn o t), in the form We also note the following useful relations.From the definition of spherical harmonics we have ) and the Wigner matrices satisfy D (2) n,m ) * .These relations are easily shown to imply that the sum in equation ( 9) is real.Accordingly it also follows from equation ( 9) the primary´s tidal response can be assembled by summing individual responses to the real parts of harmonically varying tidal potentials of the form with δ i,j denoting the Kronnecker delta.As the β dependence is only through the factor d2 n,m we find it convenient to remove it from the response calculation by considering the tidal response to potentials of the form and then reintroduce it through applying a response amplitude factor when required.Note that, where the integral is taken over the volume of the primary.
For |R| >> |r| we perform an expansion in inverse powers of |R| and spherical harmonics.
The dominant term then takes the form of a quadrupole as (see PS) where R = |R|.We may now assemble the response to the tidal potential given by equation ( 9) using linear superposition with the result The component of the of the specific torque in the Z direction is then 2.5 The potential perturbation in the orbit frame and its time average at the location of We make use of ( 6) -( 8) to transform the spherical harmonics expressed in terms of coordinates associated with the stellar frame that appear in equation ( 15) to a linear combination of spherical harmonics expressed in terms of coordinates defined in the orbit frame.Equation (15 ) specifying where we denote Setting θ ′ = π/2, ϕ ′ = Φ and R = a in equation( 17) we obtain the potential perturbation at the location of the companion as Neglecting any possible slow variation of ϖ, the time average of, ψ ′ , evaluated at M p found by integrating around the orbit is Noting the connection between T z and ψ ′ indicated by equation( 16), the time average of T z is given by

The force per unit mass acting on M p
We write the components of the force F acting on M p resulting from ψ ′ in the orbit frame as These act in the radial direction, the azimuthal direction and the direction perpendicular to the orbital plane or Z ′ direction respectively.We have These are readily obtained from equation( 17) from which we obtain where

Gauss equation for the evolution of the semi-major axis
The Gauss equations express the rate of change of orbital elements in terms of the components of the force per unit mass F, assumed to act on M p .Because our coordinate system is a non inertial accelerating reference frame with origin at the centre of mass of the primary, and the force acting on M p produces an equal and opposite reaction acting on the primary, the force per unit mass, F, determining the evolution of the orbit is obtained from F by dividing by the reduced mass, For a solar mass star with a planetary mass companion this correction, though implemented, is negligible.
The Gauss equation governing the rate of change of the semi-major axis is where, e, is the orbital eccentricity which we set to zero for circular orbits.Setting e = 0 in (26), noting the above remarks about determining F, making use of equation ( 24), and then taking the time average, we obtain Using the form of B n,m given by ( 15) together with the relations given in Section 2.3 and noting that the overlap integral is such that, , it follows that the expression ( 27) is real and can be written in the form We remark that because the Wigner matrix with real elements, are equal to unity for any possible n.This means that for a spherical non rotating star for which Q j,2 does not depend on j, ⟨da/dt⟩ does not depend on β as expected.
The upper panel shows the resonance curve (with resonance frequency ω 0 and full resonance half width, D) associated with the kinetic energy response (erg) for the nr = 0 r-mode with n = −1 and m = −2.The lower panel shows the corresponding resonance curve obtained for the viscous dissipation rate (erg/s).Note that the subscripts n, and m have been removed from the forcing frequency ω f .

Gauss equation for the evolution of the orbital inclination
The Gauss equation for the rate of change of the inclination is We recall that the orbital inclination i is the angle between the orbital angular momentum L and the total angular momentum J of the system.Substituting F ⊥ , which is given by equation ( 25) into equation ( 29), taking into account the comments at the beginning of Section 2.7, and setting Noting that Φ = n o t + ϖ and taking the time average we obtain Using the expressions given with equation ( 15) to express the result in terms of the overlap integrals and performing the sum over m, we find In this case we note that because the Wigner matrix with real elements, d 2 n,m , is unitary the sums (2) j,m are zero for n ̸ = m.This means that for a spherical non rotating star for which Q j,n does not depend on j, for any possible n, we have ⟨di/dt⟩ = 0 as expected.

DISCUSSION OF THE EVOLUTION EQUATIONS FOR SEMI-MAJOR AXIS AND INCLINATION
Following the discussion in Section 5.4 of PS, we note that the overlap integrals appearing in equations ( 28) and ( 32) can be related to the mean rate of change of the kinetic energy associated with forcing due to the real part of the corresponding potential given by ( 12).This relation is given by a minor adaption of equation ( 30) of PS2 in the form This quantity is expected to be negative definite.
Substituting into equation ( 28) we obtain Here we recall that ω f,j,m = −mn o + jΩ s .Similarly from equation ( 32) governing the mean rate of change of inclination we obtain We note that the β dependence in equations ( 34) and ( 35) is contained entirely within the Wigner small d matrix elements.In the case of, a, the rate of evolution is determined by the five rate of kinetic energy changes (dE kin.,j,2 /dt) for j = −2, −1, 0, 1, 2. In the case of the inclination, the evolution also depends on these together with (dE kin,j,0 /dt) for j = 1 and j = 2, Contributions from negative values of j are related to these and in this case there is no contribution from j = 0.

The relation between ⟨di/dt⟩ and ⟨dβ/dt⟩
The angles i and β are related by the conservation of the total angular momentum, J = L + S.
Setting J = |J|, L = |L|, S = |S|, and considering the components of J, perpendicular and parallel to L, we have (see IP) Differentiating the above expressions with respect to time while using the fact that J is conserved, making use of ( 34 3.2 Estimating ⟨dβ/dt⟩ when n o /Ω s is large We are primarily interested in the situation when the orbital period is significantly shorter than the rotation period such as for a hot Jupiter with an orbital period of a few days and a solar type star with rotation period an order of magnitude larger. In this case we suppose that the dominant contribution to the righthand side of equation ( 38) comes from the forcing terms with m = 0 and hence zero forcing frequency, j = ±1, and j = ±2.
We remark that zero frequency forcing with j = 0 corresponds to the production of a static tide which is not expected to result in orbital evolution.
Forcing with m = ±2 occurs with an associated forcing frequency a factor ∼ n o /Ω s larger in magnitude than is the case for m = 0. Thus contributions from such terms are expected to be smaller than those arising from m = 0 by at least this factor given comparable energy dissipation rates.In addition the terms on the right hand side of (38) ∝ d (2) j,1 d (2) j,2 cancel when summed over j in the limit Ω s → 0. Thus these terms, derived from m = ±2 forcing, are expected to be smaller in magnitude than terms derived from zero frequency forcing by a factor ∼ (n o /Ω s ) 2 for the comparable energy dissipation rates we find.
Thus we neglect contributions arising from terms with m = ±2 and write In this case by making use of ( 33) and noting that Q −j,0 = Q * j,0 , and dE kin.−n,−m /dt = dE kin.n,m /dt, we can write (39) as Making use of the small d Wigner matrix elements given in appendix A, we can write this as This compares with a −1 ⟨da/dt⟩ evaluated in the limit β = 0 by making use of equation ( 34).We find which also applies in the limit Ω s → 0. Equations ( 41) and (42 give two equations for β and a, or equivalently L. These must be complemented by making use of the expression for the square of the magnitude of the total angular momentum which has been assumed to be conserved in the form J 2 = L 2 + S 2 + 2LS cos β = (S + L cos β) 2 + L 2 sin 2 β.This can be used to specify ,S, and so complete the system.
We may write where the moment of inertia I = kM * R 2 * and we have assumed M p /M * ≪ 1. Making use of ( 43) equation ( 41) can be written Specialising to the case when the orbital angular momentum significantly exceeds the spin angular momentum which is for the most part appropriate here, only the second term in the first set of brackets in (44) need be retained.Thus we have.
We can compare the time scales for the evolution of the semi-major axis and the inclination angle by considering 3.2.1 Dependence on β and an estimate of tidal evolution time scales.
From equation (45) β has stationary points at β = 0, π/2, and π.The first and the last of these are stable while the one corresponding to a polar orbit is unstable.Even so evolution in its neighbourhood will be slowed down suggesting a possible accumulation of systems with polar orbits if partially effective tides operate on a uniform initial distribution of misalignments (see Section 6 below).Note that if we take into account the magnitude of the spin angular momentum in comparison to the orbital angular momentum and use equation ( 44) the unstable point is shifted to be slightly in excess of π/2 when this is small.
If the energy dissipation rates are comparable as found here3 and n o ≫ Ω s then, in the limit of small β, we have the simple result For Ω s /Ω c = 0.004, a/R s = 10, k = 0.08, and M p /M * = 0.001, This yields a ratio of 59.Thus we expect the inclination angle to evolve about 60 time faster than the semi-major axis for small β.On the other hand Equation (42) gives |a(⟨da/dt⟩) −1 | ≈ 2.5 × 10 12 (|dE kin,2,2 /dt|/6 × 10 23 cgs) −1 y.
From the results presented below ( see also PS) a characteristic value of |dE kin,2,2 /dt)| ∼ 6 × 10 23 cgs.This indicates that |β(⟨dβ/dt⟩) −1 | ∼ 4 × 10 10 y suggesting that for the parameters adopted the evolution of β will be be modest over the lifetime of a solar type star (but see Section 5.5 below for the scaling with planet mass and orbital period ) .From the above considerations we  evolved to the stage when X c = 0.40, when its radius had become R * = 1.0051R ⊙ and its luminosity L * = 1.01L ⊙ , these parameters being quite similar to those of the current Sun.The calculations were performed with version r22.11.1 of the MESA stellar evolution code (Paxton et. al. 2015).
We use this stellar model to calculate the tidal evolution of a solar type star and a hot Jupiter in close orbit.The convective envelope extends from r/R * = 0.727 to r/R * = 0.999.The moment of inertia of the star is I = 7.0413 × 10 53 g cm 2 and the critical stellar spin rate Ω c = 6.2176 × 10 −4 s −1 .When comparing the tidal response of this model to a companion in a circular orbit with a = 10R s to that obtained for a more evolved star (with X c = 0.20), we found very similar  viscous dissipation rates while the stellar moment of inertia differed by only a few percent.Hence we adopted the response for the stellar model with X c = 0.4 to apply for further main sequence stellar evolution.
As in PS we consider a Jupiter mass companion in a circular orbit.The numerical procedures followed to obtain the tidal response in the star (generated by the orbiting planet) are described in detail in PS and references therein.In solar type stars the tidal torque on the star is directly related to the viscous dissipation of the kinetic energy of the tidal oscillations generated in the convective envelope.Radiative dissipation in the radiative regions of the star is expected to be much smaller and can be neglected (PS).The artificial viscous damping introduced in the radiative core to deal  with very short wavelength gravity waves, that would otherwise be unresolved (see below), is not included in the computation of the the stated viscous dissipation rates on account of its resolution dependence.
Following PS a turbulent viscosity was assumed.The kinematic viscosity ν(r) was taken from Duguid et. al. (2020) as We introduced a thin layer with artificial viscosity in the transition layer between the radiative core and the inner boundary of the convective envelope at r/R * = r cb /R * = 0.727 by extrapolating the kinematic viscosity ν cb at this boundary inwards according to ν(r) = ν cb exp{− ((r cb − r)/(c * H P )) 2 } up to the adopted minimum artificial viscosity of ν min = 10 9 cgs in the radiative core, where H P is here taken to be the pressure scale height at the convective boundary and the scaling factor is c * = 0.05.In all cases, except for the forcing with n = −2, the adopted minimum core viscosity of ν min = 10 9 cgs gives sufficient damping to enable the determination and calculation of the l ′ = 1 resonance for forcing with However, forcing with (n = −2, m = 0) requires ν min = 10 10 cgs in the radiative core to damp grid oscillations of the displacement vector in the core.We comment that in the latter case ν min /(R 2 * Ω c ) ∼ 3.25 × 10 −9 .Thus for the smallest angular velocity considered, namely, Ω s = 0.002Ω c , the Ekman number ν min /(R 2 * Ω s ) ∼ 1.6 × 10 −6 .As indicated above the effect of this viscosity is not included in the presented total rate of dissipation which is dominated by effects in the convection zone (see also PS).
As in PS the viscous force is derived from the viscous stress tensor for compressible flow Σ i,j expressed in spherical coordinates, whereby ρ −1 ∇ • Σ i,j gives the viscous force per unit mass The viscous dissipation rate in the convective envelope follows by calculating4 where the volume integral is taken over the convective envelope (up to r = 0.995R * thus omitting the low mass superadiabatic surface layer).
When applying the frozen convection approximation, as was implemented in PS, the tidal dissipation rate can become unrealistically large in the outer convective envelope where non adiabatic effects play a significant role in the evolution of the entropy variations induced by the tidal perturbation (see Bunting et. al. 2019).This effect was much less significant for the calculations of PS which were mainly focused on r mode resonances and their neighbourhoods, with associated perturbations located mainly in the radiative core.
In the outer convective envelope the time scale associated with convective energy becomes short enough to result in the frozen convection approximation becoming invalid.Unlike the situation when frozen, the convection is able to adjust and smooth out the rapid entropy variations that would be produced by rapid variations of the divergence of the perturbed radiative flux and so reduce the amplitude of the response (see Bunting et. al. 2019).Taking this effect properly into account requires a rigorous treatment of convection that is currently unavailable.
To proceed we adopt a heuristic procedure based on physical arguments.We. assume that the convective flux perturbation , F ′ c , obeys a simple relaxation equation where F ′ rad and F ′ eq − F ′ rad , are respectively the radiative flux and the equilibrium convective flux perturbations the relaxation process leads to, with τ c being the convective time scale.The quantity eq is also the total flux the system would relax to in the limit of vanishing convective time scale.For a forcing frequency, ω f,n,m , we thus have We assume that rapid variations in F ′ rad produced by density and temperature variations occur while the system relaxes towards a slowly varying total flux, F ′ eq , such as would be expected for a near isentropic convection zone.Accordingly we assume that its divergence may be neglected.
In addition we also neglect variation of τ c .Accordingly (51) for the divergence of the total energy flux perturbation becomes This procedure leads to the application of a complex reduction factor ϵ c in equation.( 17) of ). Notably this is unity in the limit of long convective time scale corresponding to frozen convection while producing, as expected from the above arguments, a reduction factor for fluctuations in the radiative flux when the convective time scale is short.

NUMERICAL RESULTS
All numerical calculations of the stellar tidal response are calculated for a fixed semi-major axis a = 10 R * , corresponding to an orbital period P orb = 3.697d.The forcing frequency for (n, m) forcing in the rotating stellar frame is ω f,n,m = −mn o + nΩ s .For m = 0 forcing the tidal response depends on Ω s only.In this case the response for other values of a can be obtained by applying the scaling factor (a/R * ) 3 for the tidal perturbations and the factor (a/R * ) 6 for the viscous dissipation rate.When m ̸ = 0 and n o does not correspond to the above orbital period the same scaling factors may be applied to obtain results for the value of a corresponding to the specified n o .

r mode resonance
Toroidal or r mode resonances are expected to occur when the forcing frequency ω f,n,m = −mn o + nΩ s is close to 2nΩ s /l ′ (l ′ + 1), where l ′ is an integer (see PS) with the relative deviation → 0 as Ω s → 0. This means that the precise location of a resonance can be found for a small value of n o .The case of interest here has l ′ = 1, with n = ±1 together with m = ±2.Without loss of generality we focus on the case n = −1, and m = −2.The resonant frequencies are then very close to −Ω s which would be expected for a rigid tilt mode for which to within a constant of The dominant l ′ = 1 r-mode resonance was calculated for a spin rate Ω s /Ω c = 4 × 10 −3 by forcing with n = −1 and m = −2 and zooming in on the resonance by searching the forcing frequency for which the kinetic energy in the star becomes maximal (the procedure followed by PS).The system's semi-major axis was kept fixed at a = 10R s noting the possibility of later scaling.
The properties of the resonance are illustrated in Fig. 1 which shows the resonance curves for the kinetic energy and the viscous dissipation.As in PS these curves are fitted by the functional form I 0 /(1 + ((ω f − ω 0 )/D) 2 ), Where the fitted parameters I 0 , ω 0 , and D are indicated in Fig. 1.
Contour plots for the tidal displacement components at the resonance frequency are shown in Fig.
2. It will be seen that the displacement is mainly toroidal with l ′ = 1 for which ξ θ is independent of θ and ξ ϕ ∝ cos θ.Thus for the horizontal components of the displacement there are no nodes in θ.This corresponds to the dominant or fundamental mode with n r = 0 in the notation of PS.It is the mode with eigenfrequency in the rotating frame closest to Ω s in magnitude.However, it is not close to a rigid tilt mode, which does not exist for our model, and for which ξ θ and ξ ϕ would be ∝ r and thus have no nodes.Fig. 2. indicates at least one node in both Re(ξ θ ) and Im(ξ ϕ ).
Note that the resonance width is much smaller than the magnitude of the deviation of |ω 0 | from Ω s .This has the consequence that the calculations below, are typically significantly into resonance wings.Here we remark that, although not included in the viscous dissipation calculation, the calculated resonance width is significantly increased by the presence of artificial viscosity in the radiative core.As the mode is predominantly located in the radiative core, this causes the profiles illustrated in Fig. 1 to be significantly broadened.Although the structure of the centre of the resonance, which in any case is likely to be affected by nonlinear effects (see PS) is modified, the wings and consequently the discussion below are to a very good approximation unaffected.6 5.2 General tidal response for rotation rates such that 0.003Ω c < Ω s < 0.010Ω c As required for an application of equations ( 42) and ( 44) to determine the evolution of a and β, calculations of the tidal response and consequent viscous dissipation rates are presented with forcing frequencies with (n, m) = (−1, 0), (−2, 0) and (−2, −2).Discrete values of Ω s /Ω c between 0.003 and 0.010 corresponding to rotation periods between 39.0d and 11.7d were adopted, see tables 1-4.
For other values of Ω s , and the corresponding ω f,n,m , the viscous dissipation rates were obtained by cubic spline interpolation/extrapolation, They are plotted in Figs. 3 and 4.These results enable the Runge-Kutta time integration of the evolution equations of the system ( Press et al. 1996) (see the next section).
The responses to forcing with (n, m) = (−1, 0) for Ω s = 3.0×10 −3 Ω c .and Ω s = 5.0 × 10 −3 Ω c are illustrated by the contour plots in Figs. 5 and 6 respectively.An r-mode response in the radiative core may potentially be excited.Though, as explained above the calculations presented here are significantly into the resonance wing.In spite of this Figs.Table 3. Results for forcing of the star with (n, m) = (−2, −2) are tabulated: the annotation is the same as for table 1 tions still have a strong toroidal component.The viscous dissipation rate as a function of Ω s /Ω c for (n, m) = (−1, 0) plotted in the lower panel of Fig. 3 indicates that resonance is approached as Ω s decreases.
Similar rates of dissipation are obtained for forcing with the component of the tidal perturbation with (n, m) = (−2, 0).Contour plots illustrating the response for Ω s = 4.0 × 10 −3 Ω c .are provided in Fig. 7.In this case there is a strong toroidal component with l ′ = 3 as expected.
The forcing frequencies ω f,−2,−2 for (n, m) = (−2, −2) lie outside the inertial range and result in the excitation of g-modes of high radial order in the radiative core as indicated in the contour plots in Fig. 8.In this case the perturbation is mainly spheroidal with l ′ = 2 for which ξ ϕ ∝ sin θ.
For r/R s < 0.3 the response could not be resolved and it was accordingly artificially damped.)( Ėkin,2,2 ) (assuming β = 0) is found to be significantly smaller than the adopted (Skumanich 1972) magnetic spin down rate ṠMB = −6.606× 10 47 Ω 3 s for all spin rates considered.

Magnetic braking and tidal evolution of the system
The viscous dissipation rates associated with the obliquity tides with (n, m) = (−1, 0) and (−2, 0) tend to drive the evolution towards alignment of L and S. We remark that the energy dissipation rates provided in tables 1 -3 and Figs. 3 and 4 may be inserted into equation ( 45) in order to obtain the evolution rate of β.
Since P orb is taken to be much smaller than P spin the generated tide with (n, m) = (−2, −2) causes the stellar spin rate to increase.However, standard Skumanich magnetic braking (Skumanich 1972) of the stellar spin is expected to dominate that process (see table 4).Thus this spin-down can in principle drive the system closer to an r-mode resonance with l ′ = 1.The dominance of magnetic braking has the consequence that the tidal evolution is almost entirely driven by the responses with (n, m) = (−1, 0) and (−2, 0).These are associated with a forcing frequencies −Ω s and −2Ω s respectively, and thus the relationship to the r mode resonance discussed above, namely being located in its wings, does not depend on the orbital period.In addition the affect on the evolution of changing the latter, to a very good approximation, can be taken into account by scaling the amplitude of the tidal response.This is in the same manner as changing the planet mass.In this way our calculations below undertaken for an initial orbital period of 3.697d can be extended to apply to other values.A consequence is that the rate of evolution changes.This aspect is discussed further in Section 5.5.
Unfortunately, the physics of magnetic braking is complicated and not well understood.Magnetic braking in evolved solar type stars generally follows the observationally derived Skumanich expression (Skumanich 1972) whereby the stellar spin angular momentum S decreases as ∝ t −1/2 due to magnetic braking with Ṡ = ṠMB = f M B Ω 3 s ≡ −κS, defining the quantity κ.The con- stant f M B = −6.60 × 10 47 (cgs) is adapted to obtain the current solar spin period P spin = 27 d.
Kepler observations of open clusters of known ages show solar type stars that follow Skumanich spin-down but there are also fast rotating stars that are not consistent with the standard magnetic braking expression (see eg. Gossage et al. 2023).

Numerical calculation of the evolution of the system
We follow the tidal evolution of the system by calculating the rate of change of a given by equation ( 42) and the rate of change of β given by equation ( 45).We remark that including magnetic braking is problematic as the angular momentum is no longer conserved as has been assumed above.Here we deal with this issue in a numerical treatment of the evolution by adopting an approach based on operator splitting (Strang 1968).Up to now, L, and S obey equations of the form, dL/dt = T tidal , dS/dt = −T tidal , and dJ/dt = d(L + S)/dt = 0, where the last of these gives the conservation of angular momentum and T tidal is a tidal torque.In order to incorporate spin down we modify the system to read The first system leads to equations ( 42) and ( 45) with total angular momentum conservation as given above.The second system is such that S preserves its direction while decreasing in magnitude.Thus β does not change.Also L and hence L do not change Accordingly the right hand sides of ( 42) and ( 45) are identically zero and we now have We remark that although it's magnitude changes J remains coplanar with L and S throughout . 7 Now let the operator O 1 (∆t) advance the first system (53) through a time step, ∆t, and the operator O 2 (∆t) advance the second system (55) through the time step ∆t.The splitting procedure advances a step, 2∆t, correct to second order, by applying the sequence of operators, (Strang 1968).This procedure can be seen to be equivalent to solving the system of three equations (42)), (45), and (56) directly. 8We implement it using an adaptive step size controlled fifth order Runge-Kutta ( Press et al. 1996) subroutine RKQC to perform the first step and subroutine RK4 to perform the second step and advance J using (56).
In Fig. 9 we show the time evolution of the angle β for a solar type star and a planet with mass M P = M J applying the Skumanich magnetic braking with f M B = −6.60 × 10 47 cgs starting from an inclination β = 40 • .It is found that with this induced stellar spin-down rate only a small reduction of ∼ 2 • in β occurs.This rate of evolution is consistent with the discussion in Section 3.2.1.As expected from the discussion in that Section the orbital period is found to decrease very slightly to 3.683 d.The system has not shifted sufficiently close to resonance with the l ′ = 1 r-mode spectrum for larger changes to β to occur during the Main Sequence phase.
7 Using this result we could start by taking, J 0 , the initial value of J to define the Z ′′ axis of our coordinate system rather than J itself..Note that L, S.J 0 and J are all coplanar.The inclination i will then be respect to J 0 and this can be shifted to be respect to J by taking into account the angle between J and J 0 .This approach yields identical conclusions as that of the splitting approach.
8 .The same primary model is assumed throughout.Thus changes as a result of stellar evolution that were found to be small are neglected 5.5 Effect of increasing M p or decreasing P orb .
In the context of these calculations we note that from Fig. 9 that for a rotation period, P rot = 33.6d, the angle β decreases by ∼ 1 • in 10 9 y.This is again consistent with the estimate in Section 3.2.1.
However, it is important to stress that these results apply to a planet with mass, M p , equal to one Jupiter mass with an orbital period of 3.697d.From equation ( 45) it follows that for fixed M * , the tidal evolution rate is ∝ M 2 p /P 4 orb .Thus it will be increased by ∼ one order of magnitude if M p is increased by a factor of ∼ 3. Evolution of β would then become significant over the main sequence life time were the rotation period to be maintained at ∼ 34d.Given that negligible change to the orbital period is still expected during the tidal evolution, similar changes in β would be expected were the initial orbital period reduced by a factor ∼ 1.75.
In order to illustrate the above discussion we calculated the the evolution of β for a system for which the planet mass M P = 3M J The results are illustrated in Fig. 10.In this case β decreased from 40 • to 24.5 • during the Main Sequence stage.The orbital period decreased to 3.651 d.The corresponding evolution of the spin period (right panel in Fig. 10) shows that the stellar rotation period attained the value P spin ≃ 32 d.
In this context we remark that Attia et al. (2023) note that the statistics of hot Jupiter missalingments indicate increased significance of tidal effects for higher masses and shorter orbital periods.
While stressing the uncertainties resulting from the use of a simplified stellar model as well as the crude treatment of convection, our results indicate that while it might have significant effects in some circumstances, the obliquity tide is unlikely to produce strongly aligned hot Jupiter systems overall.

DISCUSSION
In this paper we investigated the tidal interaction between a giant planet on a circular orbit around a solar type primary star.Extending the work of PS, we considered the situation when the orbital and spin angular momenta were misaligned.We obtained equations governing the inclination angle between the spin and orbital angular momenta β and the semi-major axis which depended on the energy dissipation rates due to tidal perturbations associated with forcing frequencies, ω f,n,m = nΩ s − mn o , with (n, m) = (−1, 0), (−2, 0), and (−2, −2).We focused initially on the case where M p was one Jupiter mass and the orbital period was 3.697d.
For the first of these for which (n, m) = (−1, 0), the perturbation is stationary in the non rotating frame and the spectrum of r modes with l ′ = 1 which has nearby eigenfrequencies may potentially affect the response.These modes have eigenfrequencies close to the frequency of a putative rigid tilt mode.However, we recall that such a mode does not exist for the model adopted apart from in the limit Ω s = 0.The properties of the fundamental r mode ( with eigenfrequency closest to that of the putative tilt mode) were investigated in Section 5.1.We found that the resonance width was extremely small such that even the small frequency mismatch associated with the forcing resulted in a response far into the wings.This was found to be the case for stellar rotation rates varying between 11.8d and 39.0d (see table 1 in conjunction with Fig. 1).This can be understood in the following simple manner.Consider the case illustrated in Fig. 1.The dimensionless resonant width is determined by the dissipation rate and, even though it increased by artificial viscosity, is ∼ 7.5 × 10 −7 ω 0 ( see Fig. 1 and PS).Whereas the relative frequency separation from the resonant frequency −Ω s is ∼ (Ω s /Ω c ) 2 ∼ 1.6 × 10 −5 This has the consequence that the full tidal response was essentially non resonant with the energy dissipation rates associated with all relevant values of ω f,n,m considered being comparable.
Given this we estimated tidal evolution time scales from the governing equations in Section 3.2.1 obtaining |a(⟨da/dt⟩) −1 | ≈ 2.5 × 10 12 y and |β(⟨dβ/dt⟩) −1 | ∼ 4 × 10 10 y taking M p to be one Jupiter mass.The time scale for changing the orbital semi-major axis was thus estimated as a factor ∼ 60 longer than that for changing β.These estimates were later confirmed by numerical.calculations of the orbital evolution that in addition took into account magnetic braking in Section 5.4.
From equation (45) the evolution of β is towards 0 for β < π/2 and towards π for β > π/2 with an unstable stationary point at β = π/2.Slow evolution in the neighbourhood of that point could result in a relative accumulation of systems in near polar orbits.This could occur if tides are effective without an initial preference arising through the formation process.Albrecht et al. (2021) and Attia et al. (2023) have found statistical evidence for such an accumulation.However, Siegel et al. (2023) using a different approach do not find strong support for this at present.This should be resolved by future work.
The spin up rate of the central star induced by tides was found to be very much less in magnitude than the estimated spin down rate arising from magnetic braking (see table 4).Thus it is a reasonable approximation to consider the evolution of β and Ω s .to occur while the orbit remains fixed.
It can be seen from integrating equation ( 45) that this has the consequence that if β 0 < π/2 and β f are the initial and final values of β for the tidal evolution during the main sequence lifetime, for a given primary a function G(β 0 , β f ) which approaches the form (β f − β 0 ))/β 0 as β f → β 0 scales as M 2 p /P 4 orb .Notably in this approximation the scaling with P orb enables scaling to different orbital periods without the need for further tidal response calculations 9 .
Thus although we estimated from our calculations that for a one Jupiter mass planet and a primary rotation period ∼ 33.6d, β would change by about 1 • in 10 9 y, this would increase by about an order of magnitude for M p ∼ 3 Jupiter masses.Thus the decrease of β of from 40 • by about 2 • we found for a one Jupiter mass planet over a main sequence lifetime as illustrated in Fig. 9 increases to ∼ 15 • for 3 Jupiter masses.This is consistent with the finding of evidence supporting the increased efficacy of tides for larger planetary masses and shorter orbital periods by Attia et al. (2023).
The alignment distribution of close in giant planets is potentially determined by a multitude of processes affecting individual objects in different ways (see eg. Siegel et al. 2023;Wright et al. 2023;Wu et al. 2023).Some planets may undergo quiescent disc migration with only modest dynamical interactions leading to modest spin-orbit misalignment.Others may undergo more violent interactions leading to larger misalignments.It has been suggested that tidal interactions may be responsible for greater alignment of planets around cool stars but not hotter stars beyond the Kraft break on account of the lack of an envelope convection zone (eg.Albrecht et al. 2012).
However, strong dynamical interactions may be preferred for stars beyond the Kraft break (Wright et al. 2023;Wu et al. 2023).In addition the distribution of warm Jupiter misalignments for which tidal effects are expected to be ineffective indicates a quiescent formation process can occur.
Although not obviously required ab initio, the discussion of Attia et al. (2023) indicates some influence of tides, though this seems to be modest.In support of that view the results presented here tend to indicate a potentially significant influence of tides, but only for giant planets with relatively large masses and short orbital periods.
However, the limitations and uncertainties associated with our results need to be emphasised.These relate to the use of a simplified stellar model that neglected centrifugal distortion, being equivalent to one immersed in a fixed background potential designed to cancel out the centrifugal potential.The tidal forcing with (n, m) = (−1, 0) may be more strongly affected by resonance with r modes leading to faster tidal evolution in a more realistic model.In addition there are significant uncertainties associated with the effective viscosity arising from convection.These are issues to be addressed in future work. 9The case β 0 > π/2 can be considered through the mapping β → π − β.

2. 4
The perturbation to the external gravitational potential due to the primary After separating out the exponential factor exp(−im (n o t + ϖ + γ)), following PS the perturbation to the external gravitational potential produced by the forcing potential d (2) n,m (β)U n,m 1 at position vector R in the stellar frame is

Figure 5 .Figure 6 .
Figure 5. Contour plots in the primary's meridional plane ϕ = 0 showing the response to forcing with, U −1,0 , (n, m) = (−1, 0)), for Ωs = 3.0 × 10 −3 Ωc.For these values of n and m the forcing frequency is −Ωs corresponding to a time independent perturbation in the non rotating frame.This calculation was done for a/R * = 10.The response for other values can be obtained by applying the scaling factor (a/R * ) 3 .The Cartesian coordinates along the two axes indicate the relative radius r/R * .The vertical colour bars on the right indicate the local value of sign(|ξx| 1 4 , ξx), where ξx is the component of the displacement vector illustrated.
) whereby the convective mixing length L mx = α|H P | is here scaled by the parameter α = 2.The local pressure scale height |H P (r)| and the local convective velocity v c (r) are taken from the MESA input stellar model.Any mismatch of the timescale of the forced oscillations (P osc = 2π/ω f,n,m ) and that of the turbulent convection (τ c = 1/ |N 2 |), where N is the Brunt-Väisälä frequency, is taken into account by the term in the denominator raised to the power s = 2.

Figure 9 .Figure 10 .
Figure 9.The left panel shows the time evolution of the angle β (in degrees) between S and L driven by Skumanich magnetic braking with f M B = −6.60 × 10 47 (cgs).The right panel shows the corresponding evolution of the stellar spin period in days.
dL/dt = T tidal , dS/dt = −T tidal − κS, and dJ/dt = d(L + S)/dt = −κS.(54) Here we have introduced the spin down torque −κS.To proceed by operator splitting we split (54) into two systems.The first system is given by (53) and the second by dL/dt = 0, dS/dt = −κS, and dJ/dt = d(L + S)/dt = −κS, mn o , is the forcing frequancy in the non rotating frame which corresponds to the the forcing frequency, ω f,n,m = nΩ s − mn o , in the frame corotating with the star.

Table 1 .
Results for forcing of the star with (n, m) = (−1, 0): the columns moving from left to right contain the dimensionless spin rate Ωs/Ωc, the kinetic energy, the Dissipation rate − Ėkin,n,m ≡ − Ėkin,−n−,m , the dimensionless forcing frequency ω f,n,m /Ωc and spin period in days.

Table 2 .
5 and 6 indicate that the perturba-Results for forcing of the star with (n, m) = (−2, 0) are tabulated : the annotation is the same as for table1