Wave breaking in a solar-type star

Three-dimensional simulations of internal wave breaking and the fate of planets around solar-type stars

Adrian J. Barker
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences,
Wilberforce Road, Cambridge CB3 0WA, UK
E-mail: ajb268@cam.ac.uk
Accepted 2011 February 3. Received 2011 January 26; in original form 2010 November 3

We study the fate of internal gravity waves approaching the centre of an initially non-rotating solar-type star, by performing three-dimensional numerical simulations using a Boussinesq-type model. These waves are excited at the top of the radiation zone by the tidal forcing of a short-period planet on a circular, coplanar orbit. This extends previous work done in two dimensions by Barker & Ogilvie. We first derive a linear wave solution, which is not exact in 3D; however, the reflection of ingoing waves from the centre is close to perfect for moderate amplitude waves. Waves with sufficient amplitude to cause isentropic overturning break, and deposit their angular momentum near the centre. This forms a critical layer, at which the angular velocity of the flow matches the orbital angular frequency of the planet. This efficiently absorbs ingoing waves, and spins up the star from the inside out, while the planet spirals into the star.

We also perform numerical integrations to determine the linearised adiabatic tidal response throughout the star, in a wide range of solar-type stellar models with masses in the range , throughout their main sequence lifetimes. The aim is to study the influence of the launching region for these waves at the top of the radiation zone in more detail, and to determine the accuracy of a semi-analytic approximation for the tidal torque on the star, that was derived under the assumption that all ingoing wave angular momentum is absorbed in a critical layer.

The main conclusions of this work are that this nonlinear mechanism of tidal dissipation could provide an explanation for the survival of all short-period extrasolar planets observed around FGK stars, while it predicts the destruction of more massive planets. This work provides further support for the model outlined in a previous paper by Barker & Ogilvie, and makes predictions that will be tested by ongoing observational studies, such as WASP and Kepler.

planetary systems – stars: rotation – binaries: close – hydrodynamics – waves – instabilities
pagerange: Three-dimensional simulations of internal wave breaking and the fate of planets around solar-type starsLABEL:lastpagepubyear: 2011

1 Introduction

Tidal interactions are important in determining the fate of short-period extrasolar planets and the spins of their host stars. The extent of the spin-orbit evolution that results from tides depends on the dissipative properties of the star and planet. These are usually parametrized by a dimensionless quality factor for each body, which is an inverse measure of the dissipation. This is usually defined to be proportional to the ratio of the maximum energy stored in a tidal oscillation to the energy dissipated over one cycle (e.g. Goldreich & Soter 1966).

In principal the modified quality factor111Related to the traditional by , where is the second-order potential Love number of the body. depends on tidal frequency, the internal structure of the body, and the amplitude of the tidal forcing. Unfortunately, the mechanisms of dissipation that contribute to are poorly understood. To simplify this difficult problem, nearly all studies neglect any amplitude dependence of the dissipation (except, for example Goodman & Lackner 2009). Such studies already exhibit a complicated dependence of on the tidal frequency (e.g Savonije & Papaloizou 1997; Wu 2005; Ogilvie & Lin 2004, hereafter OL04; Ogilvie & Lin 2007, hereafter OL07; Papaloizou & Ivanov 2010), and the internal structure of the body (e.g OL07; Barker & Ogilvie 2009). In this paper we extend the work of Barker & Ogilvie (2010) (hereafter BO10) to study an important nonlinear dissipation mechanism in solar-type stars.

We can decompose the tidal response of a fluid body, which in this paper we take to mean a solar-type star, into an equilibrium and a dynamical tide, defined such that the total displacement is the sum of these two displacements, i.e, that . The equilibrium tide is a quasi-hydrostatic bulge defined by


in stratified regions (Goldreich & Nicholson 1989), where is the tidal gravitational potential experienced by the body and is the gravitational acceleration. The total displacement is not well described by Eqs. 1 in convective regions (Goodman & Dickson 1998, hereafter GD98; Terquem et al. 1998, hereafter T98; OLO4). This is because in a barotropic flow (with adiabatic stratification) vorticity is conserved, so we must have , whereas , in general. The presence of a convection zone (hereafter CZ) thus implies that a dynamical tide must exist. The dynamical tide is defined as the residual response that results from the equilibrium tide not being the exact (linearised) solution to the problem, when the tidal frequency is nonzero.

The equations governing the adiabatic equilibrium and dynamical tides in linear theory are


from which it is clear that the dynamical tide is not forced directly by the tidal potential, only by the inertial terms in the equation of motion.

Dynamical tides in radiation zones of solar-type stars take the form of internal (inertia-) gravity waves (IGWs), which have frequencies below the buoyancy (or Brunt-Väisälä) frequency . These have previously been proposed to contribute to for early-type stars (e.g. Zahn 1975), where they are damped at the surface by radiative diffusion. It has also been proposed that these waves could synchronise the spin of the star with the orbit (in this case of a close-binary perturber) from the outside in (Goldreich & Nicholson, 1989). In BO10, we considered a nonlinear mechanism of tidal dissipation in solar-type stars (with radiative cores), extending an idea by GD98. The consequences of this mechanism are similar to Goldreich & Nicholson (1989), except that the star is synchronised with the orbit from the inside out.

A short-period planet excites IGWs at the top of the RZ of a solar-type star, where increases from zero with distance into the RZ from the CZ/RZ interface. There is thus a location at which , with being the planetary orbital period, at which IGWs (which have frequencies less than ) are efficiently excited. These waves propagate downwards into the radiation zone (hereafter RZ), until they reach the centre of the star, where they are geometrically focused and can become nonlinear. If their amplitudes are sufficiently large, convective overturning occurs, and the wave breaks. This has consequences for the tidal torque, and the stellar . In this paper we study this mechanism, primarily using three-dimensional numerical simulations.

In BO10, we derived a Boussinesq-type system of equations that is relevant for describing the dynamics of IGWs approaching the centre of a solar-type star. We then performed numerical simulations, solving these equations in two dimensional cylindrical geometry. In this paper we extend these simulations to three dimensions, in spherical geometry, and confirm that the most important results of BO10 are not affected by this extension. We first derive a linear wave solution, which represents the waves excited by planets on short-period orbits as they approach the central of the star, within which and . A weakly nonlinear analysis confirms that this solution is not exact, though the reflection of ingoing waves is close to perfect for moderate amplitude waves. We present the numerical setup and the analysis of the simulation results. We then discuss the launching region at the top of the RZ, and the possible effects of magnetic fields on the problem. Together with BO10, we provide a possible explanation for the survival of all short-period extrasolar planets observed thus far, which will be tested by ongoing observational studies, such as WASP and Kepler.

2 Tidal potential

The tidal potential experienced by a star hosting a short-period planet can be written in standard spherical polar coordinates as a sum of rigidly rotating spherical harmonics


in a non-rotating (but non-inertial) reference frame centred on the star, where is a spherical harmonic (normalised such that the integral of over solid angles is unity) and is an amplitude. Here is the frequency in that frame, related to the tidal frequency by , where is the spin angular frequency of the star. In this paper we consider the waves excited by planets on circular, coplanar orbits, relegating any studies of the waves excited by eccentric and inclined planets to future work. In this case the dominant component of the tidal potential is quadrupolar (), with , and takes the form


where is the planetary mass and is its orbital semi-major axis. Since most short-period planets orbit faster than their stars spin, as a result of stellar magnetic braking, we take the star to be (initially) non-rotating, so that , where is the orbital angular frequency of the planet.

3 Linear theory

In this section we derive a linear wave solution, starting from the adiabatic Boussinesq-type system (BO10)


where is the fluid velocity, is a buoyancy variable (proportional to the entropy perturbation) and is a modified pressure variable. This system of equation was derived in 2D but is equally valid in 3D. We note that in this model , where is a constant that measures the strength of the stable stratification at the centre. This model is valid in the innermost of the star, which contains multiple wavelengths for the gravity waves excited by short-period planets. We define , which is the pattern speed of the forcing, and equals the orbital angular frequency of the planet. We furthermore use the non-dimensionalisation that the unit of length is and the unit of time is . Linearising about hydrostatic equilibrium in 3D spherical geometry, we seek solutions steady in the frame rotating at the angular rate , proportional to , where is the azimuthal coordinate in this frame. This leads to the following equations:


To obtain the linear solution, we expand scalar quantities (i.e.  and ) in spherical harmonics , since the problem is separable in both angular coordinates. In Eqs. 1014 we have already included the -dependence of these functions (). We thus take , so that the total functions are expanded onto , and similarly for and , where from now on we drop the hats on the radial functions. The remaining velocity components are expanded onto angular functions as appropriate to satisfy Eqs. 1014. The relation


follows from incompressibility and the result


This enables us to derive the linear differential equation


whose solutions can be written in terms of Bessel functions of half-integer order (alternatively these can be written as spherical Bessel functions, or they can be reduced to elementary functions). The corresponding (total) linear solution for a standing wave in 3D can be written (where real parts are assumed to be taken)


where is an amplitude, and


Ingoing wave (hereafter IW) and outgoing wave (hereafter OW) solutions can be obtained by replacing the Bessel function of the first kind by equivalent Hankel functions of the first () and second kinds (), respectively. Note also that


Starting from Eq. 46 of BO10, we can calculate a conserved energy flux. Integrating this equation over eliminates the terms containing derivatives in and due to periodicity in , as a result of the fundamental theorem of calculus. This allows the definition of a conserved quantity proportional to the energy flux,


where the energy density per unit volume is , with the (constant) central density of the star. For linear waves, terms involving products are small, so , to this order. This leaves


Whether this is positive or negative depends on whether the wave is ingoing or outgoing.

Substituting the linear solution into Eq. 25 provides a simple expression for the flux of a single wave:


with corresponding energy flux , and angular momentum flux . We define the complex amplitudes of the IW and OW to be and , respectively. Eq. 26 follows from the orthonormality of spherical harmonics and the Wronskian of the Hankel functions


The particular wave that we will study has , and the components


where .

3.1 Criterion for overturning isentropes

A condition in 3D for isentropic overturning can be derived from considering when the radial gradient of the entropy becomes negative, i.e., when . This is equivalent to the criterion that , which can be shown by using Eq. 13. To correlate our notation with the appendix of OL07, we define the dimensionless nonlinearity parameter , such that overturning occurs if . This is defined such that the radial velocity is the real part of


in the dimensionless units that we have been using in this section (this is equivalent to Eq. 29).

Overturning is achieved at the centre when the radial velocity in the wave takes a maximum value at its innermost peak at . This value is used to compare it to the magnitude of the radial velocity achieved in numerical simulations at the onset of wave breaking. In these dimensionless units, we similarly require or . Note that there is not such a simple interpretation of the criterion on as in 2D, though the value is quantitatively very similar222In the figures below, we use a different nondimensionalisation of the velocity, by normalising it to the constant (asymptotic) radial phase velocity in the wave. this is identical to BO10, and gives values exactly 1/2 of those in the units of this section.. From the results of the 2D simulations in BO10, we expect the waves to undergo instability and break within several wave periods after these criteria begin to be satisfied.

3.2 Weakly nonlinear theory333I would like to thank Gordon Ogilvie for suggesting the calculations in this section.

The linear solution written down in Eqs. 1821 is not a nonlinear solution, unlike the equivalent in 2D. This can be shown by computing the nonlinear terms in the full Boussinesq-type system using the linear solution, in Mathematica, for example. We find that , in general, and similarly for the nonlinear terms in the momentum equation. This means that the reflection of the waves from the centre could be different than in 2D, since nonlinearities do not vanish for this wave. In this section we perform a weakly nonlinear analysis to determine the dominant nonlinear effects for small amplitudes. Since these nonlinearities do not vanish, this highlights the importance of numerical simulations for these waves approaching the centre. We describe the results of such simulations in § 5.

We propose a weakly nonlinear solution of the form


and similarly for the other variables, where . This form is adopted because we are interested in calculating whether the incoming wave generates harmonics through the quadratic (self-)nonlinearities. These additional waves (other than ) will escape to infinity and carry away a portion of the energy flux. Here we write from Eq. 29, for the wave above, and similarly for other variables. We substitute these expansions into the Boussinesq-type system and equate powers of . At each order we also equate coefficients of . At leading order only one mode is present, and we obtain the previously derived linear solution. After some algebra the solution at can be computed to give


which is an wave with complex amplitude , which has been computed using Mathematica and given in terms of .

For the wave described by Eq. 35, can be computed. The ratio of the energy flux in the outgoing wave to the ingoing wave can be shown to be approximately . We define a reflection coefficient


which measures the amplitude decay for a wave travelling from a radius to the centre, and back to . For perfect reflection from the centre , whereas complete absorption means that . The reflection coefficient for reflection from the centre, for a weakly nonlinear wave, can be computed from


This means that a weakly nonlinear primary wave (with ), will reflect approximately perfectly from the centre, with a reflection coefficient that is close to unity. However, a small fraction of the IW energy flux is transferred to waves with higher and -values, reinforcing the fact that Eqs. 1821 is not an exact solution, contrary to the analogous solution in 2D.

4 Numerical setup

We solve the Boussinesq-type system in three dimensions using the Cartesian pseudospectral code SNOOPY (Lesur & Longaretti, 2005), as in BO10 (see § 6 of that paper for further details). However, we modify the forcing and damping to take into account the -direction, and instead of forcing an wave in the equatorial plane, we now have


in Cartesian coordinates (with ), which is applied in the region . We study a region , where , in arbitrary units (not the same as § 3). For the solution is damped to zero by using a parabolic smoothing function, as in BO10.

We primarily use a resolution of , for which the simulations were possible to run on a single Intel Core i7 machine, utilising all 8 cores, with a typical run time of several weeks to resolve a hundred wave crossing times. We set and choose a typical IGW wavelength of , giving approximately 8 wavelengths within the box. The value of is chosen to be slightly larger than that used in most of the 2D calculations. This increases the number of grid points within a wavelenth of the primary wave, partially offsetting the reduction in resolution that results from using a smaller number of grid points per dimension than in 2D. Choosing larger wavelengths than this is found to result in unwanted effects from the proximity of the forcing region, which modifies the linear solution. The viscosity and radiative diffusion coefficients, at least one of which must be implemented in the code for numerical stability, are chosen to be and . An otherwise identical setup is used as in BO10, except using spherical geometry instead of cylindrical geometry. We have confirmed that the linear solution is well reproduced with this numerical setup.

5 Numerical results

In this section we describe the results of the numerical simulations. We analyse the results of the simulations using the IW/OW decomposition described in Appendix A, and reconstruct the solutions from the computed IW/OW amplitudes to compare with the simulation data, as described therein. From preliminary investigation, we find that is required for breaking, so a variety of simulations are performed with forcing amplitudes either side of this value. The basic results of these simulations are that the wave reflects approximately perfectly from the centre of the star if the amplitude of the wave is smaller than a certain critical value, which is found to correspond with that required for isentropic overturning. Above this value, wave breaking and critical layer formation occur. This picture is identical to that in 2D.

5.1 Low-amplitude simulations

Figure 1: In the top left panel we plot (blue) and (red) vs. , from a low-amplitude simulation in which the primary wave approximately perfectly reflects from the centre, with displayed in the top right panel. Below, are the velocity components and the buoyancy variable (blue), plotted together with the reconstructed linear solution (red), against radius. This is taken from a simulation with at , once standing waves have formed. The amplitude decay with radius can be explained as due to viscous damping.
Figure 2: 2D plots of (top) and (bottom) on the -plane for a simulation in which the wave perfectly reflects from the centre, for which This can be qualitatively compared with Fig. 3 in BO10.
Figure 3: 2D plot of on the plane for a simulation in which the wave perfectly reflects from the centre, for which . is of largest magnitude in the equatorial plane, where peaks.

For a low-amplitude simulation, with max() below the critical value for isentropic overturning (using ), we plot the variation in amplitudes of the IWs and OWs, and also the reconstructed solutions in Fig 1. We analyse the results using the method described in Appendix A, choosing a time once transients have been sufficiently damped and standing waves have formed. The decay in radius is roughly (though slightly smaller than) that which would be expected from viscous damping, by the fractional amount


where is the (constant asymptotic) radial group velocity, for a wave of the given wavelength. As in 2D, we have confirmed this explanation by running simulations without viscosity, which are found to not exhibit this decay (though these simulations eventually become numerically unstable if ). Using a smaller viscosity is also found to reduce the wavelength-scale oscillations around the mean slope. These result from the fact that the linear solution to the forced wave problem is no longer exact in the presence of viscosity. We find that increasing the number of grid points within each shell (by reducing the values of and ) slightly reduces the vertical extent of these oscillations, because this averages out the errors that result from the assumption that the inviscid linear solution is exact. However, increasing the number of grid points within each shell has negligible effect on the mean slope.

The spatial structure of the solutions in three dimensions in the -plane is very similar to that in two dimensions, as can be seen in Fig. 2 (which can be compared with Fig. 3 in BO10). In Fig. 3, we plot on the -plane. This shows that the magnitude of the azimuthal velocity peaks at , due to the latitudinal form of .

From the calculation in § 3.2, we expect the effects of nonlinearity to be much weaker than the effects of viscosity for small-amplitude waves which do not cause convective overturning. Since the effects of weak nonlinearity are very small, it is difficult to quantitatively confirm the results in § 3.2, using, for example, an extension of the method described in Appendix A for multiple and values. Nevertheless, we have qualitatively confirmed the result that the reflection is coherent and nearly perfect (in that ) for amplitudes below that required for overturning the stratification. As in 2D we do not observe any instabilities that act on the waves when they have insufficient amplitude to overturn the stratification. In this case, the waves can form global modes in the RZ.

5.2 High-amplitude simulations

In high-amplitude simulations, in which the wave amplitude exceeds the overturning criterion, the wave overturns the stratification during part of its cycle and a rapid instability acts on the wave, which leads to wave breaking within wave periods. This causes the rapid (within several wave periods) deposition of primary wave angular momentum, which spins up the mean flow to (which corresponds with the orbital angular frequency of the planet) and produces a critical layer. This critical layer acts as an absorbing barrier for IWs, as is shown from Fig. 4, which plots the variation in amplitude of the IW and OW, and also the reconstructed wave solutions (which can be contrasted with Fig. 1). Once the critical layer has formed, we find . The central regions are not well described by the linear model, as we would expect. However, the region outside of is well described by the linear solution, with . In this region, , so it is reasonable to assume that the IWs are efficiently absorbed near the centre. This picture is identical to that in 2D.

The picture in 3D in the -plane is very similar to that in 2D, as can be seen in Fig. 5 (to compare with Fig. 6 from BO10). However, one noticeable difference is that the primary wave preferentially transfers its angular momentum at low latitudes, close to the equatorial plane. This can be seen in Fig. 6, where we plot the angular frequency of the flow normalised to once a critical layer has formed, in both the and planes. This is a consequence of the latitudinal form of , whose magnitude peaks at , as is illustrated in Fig. 3. The critical layer absorption is observed to continue as the wave forcing is ongoing, so this differential rotation is continually reinforced by the absorption of IWs. Since there are wave motions in the region of fluid interior of the critical layer, parts of these regions spin slightly faster than (this is not seen in Fig. 6 due to the adopted colour scale). This was also observed in the 2D simulations.

Figure 4: In the top left panel we plot and vs. in a high-amplitude simulation in which the primary wave breaks () at , with in the top right panel. Below these are the velocity components and the buoyancy variable (blue), plotted together with the reconstructed linear solution (red), against radius. The central regions are not well described by the linear model, as expected.
Figure 5: 2D plot of (top) and (bottom) on -plane for a simulation in which breaking occurs with , at .
Figure 6: 2D plot of the angular velocity of the central regions in the (left) and (right) planes, normalised to the angular pattern speed (i.e., orbital angular frequency), in a large-amplitude simulation with , at . Latitudinal differential rotation is produced by absorption of ingoing waves.

6 Wave launching region at the top of the radiation zone

The main motivation for our work is to study for solar-type stars, and in particular to connect this with the survival of close-in extrasolar planets. We have demonstrated that the fate of IGWs approaching the centre of a solar-type star, outlined in BO10, is unaffected by the extension to three dimensions. A critical layer is formed by the deposition of angular momentum through wave breaking once the waves cause isentropic overturning near the centre. For lower amplitudes, the waves reflect coherently and approximately perfectly from the centre of the star, and may form global modes. In that case, the dissipation is only efficient when the system enters a resonance with a global mode. The corresponding time-averaged dissipation rate is weak, because the system spends most of its time out of resonance (T98; GD98), resulting in negligible tidal evolution of the planetary companion. Note, however, that this neglects the possibility that passage through resonance will cause wave breaking, which should be studied in future work.

If a critical layer forms, wave absorption is efficient, and global modes (of any frequency very similar to the orbital frequency) are prevented from being set up in the RZ. The tidal torque can then be computed from assuming that the IWs are entirely absorbed. In this case a calculation along the lines of GD98 for the ingoing energy flux of the waves excited at the top of the RZ is required. This estimates the tidal torque, and thus the orbital evolution of the planetary companion. In this section, we perform numerical integrations of the linearised tidal response in an extensive set of stellar models of solar-type stars with masses in the range , throughout their main sequence lifetimes. We aim to determine the tidal torque numerically, and compare it with a simple model of the launching region at the top of the RZ, which was derived in GD98 and discussed in BO10.

6.1 Numerical computation of the linearised tidal response throughout the star

In this section we solve the linearised equations governing the adiabatic tidal response (Eqs. 2 and 3) throughout the star, computing the excitation of both the equilibrium and dynamical tides numerically. This allows us to determine the ingoing energy and angular momentum fluxes in IGWs launched at the top of the RZ, and to check the validity of approximate semi-analytic formulae for these quantities, presented in the next section. This is important because the orbital evolution of a planetary companion is determined by the ingoing angular momentum flux absorbed at the critical layer.

We solve the following coupled ODEs for the radial and horizontal displacements:


An outline of the derivation of these equations is presented in T98. Note that we are ignoring the self-gravity of the entire tidal response, which is reasonable because most of the mass of the star is concentrated near the centre. This assumption is certainly valid for the dynamical tide, and is approximately valid for the equilibrium tide. In these equations, we take the tidal potential in the frame rotating with to be equal to Eq. 5 with and , so that444Note that we are defining spherical harmonics in a standard manner, normalised so that the integral of over solid angles is unity, unlike T98.


This is the amplitude of the largest tide for a circular orbit. In this frame, the displacement field is separated into radial and horizontal (non-radial) components


The tidal response is further decomposed into an equilibrium and a dynamical tide, as defined in the introduction.

We impose a free upper boundary, i.e., take the Lagrangian pressure perturbation at , so that the Eulerian pressure perturbation . Since the relation


follows from the non-radial equation of motion, this relates our variables at the surface of the star. We take an IW BC at the inner boundary at , where we match and onto the analytic solution for an IW derived in §3,


where . We take to be a free parameter, so that these equations relate the ratio of and . This solution is quite accurate when since is negligible in this region. This BC is meant to represent the IW absorption at a critical layer.

We solve Eqs. 4041 using data interpolated from a stellar model at points required by a 4th/5th order adaptive step Runge-Kutta integrator, using a cubic spline interpolation. In particular, the coefficients of , , in Eqs. 4041 are singular at the origin, so we first multiply these quantities by before interpolating their values to the locations required by the ODE integrator, using the stellar model parameters. We then divide by , after the interpolation. This is done to get the correct behaviour for small .

Our method of solution is a shooting method to an intermediate fitting point (Press et al. 1992), which we take to be the CZ/RZ interface, where we enforce continuity of the solution. The freely specifiable initial conditions for each ODE integration are chosen to be at the surface, and at the inner boundary. We use as our starting “freely specifiable” estimate at the surface, which is an accurate approximation since at .

The Eulerian pressure perturbation for the dynamical tide is


from the horizontal component of Eq. 3. The radial energy flux at each radius is


which follows from manipulating Eqs. 2 and 3 to derive an energy equation.

From Eq. 3, we derive the relation


which can be averaged over solid angle to give


This is telling us that the dynamical tide is forced by the equilibrium tide. As part of the validation of our numerical code, we have confirmed that this is accurately satisfied from the numerical solutions computed with an ingoing BC. This should adequately convince ourselves that the code is able to accurately compute the energy flux in the dynamical tide.

6.2 Results

Figure 7: Typical values of the real parts of the radial (, in blue solid lines) and horizontal (, in red dashed lines) velocity components, in . We use Model S of the current Sun, and consider the tidal perturber to be a d, , planet orbiting the current Sun.

As an illustration, we present the results of our integrations for a fiducial case with a d, planet orbiting the current Sun (for which we use Model S, described in Christensen-Dalsgaard et al. 1996) in Fig. 7. We plot the real parts of (solid blue) and (dashed red) throughout the star, which represent typical values of the radial and horizontal velocity components. The radial wavelength of the waves decreases as the waves propagate deeper into the RZ, where increases. The large increase in the velocity amplitude near the centre is evident, as predicted from the linear solutions in § 3. This can be compared with a similar calculation in T98, displayed in their Fig. 1, which our code correctly reproduces for the given planetary orbital period when a regularity condition in applied at the centre. The only difference between our calculations and theirs is that we use an IW BC, whereas they allow the waves to perfectly reflect from the centre.

In Fig. 8 we plot the ingoing energy flux, normalised to both the semi-analytic prediction of GD98 (red dashed lines), and a revised expression derived in Appendix B (blue solid lines), which will be discussed in the next section. This illustrates that the ingoing energy flux oscillates about its final asymptotic value, which it eventually approaches deep in the RZ.

Figure 8: Numerically computed , normalised to the prediction from the analytic theory in Eq. 51 (solid blue line) and GD98 Eq. 13 (red dashed line). Eq. 51 overestimates by less than , and is thus a good approximation of the ingoing energy flux.

6.3 Semi-analytical calculation of the ingoing energy flux

In this section we compare the numerically computed with the semi-analytic estimate of GD98, and present a slight refinement which is found to be appropriate for solar-type stars.

Figure 9: normalised to (solid blue line), together with our local approximation (red dashed line). This approximation is reasonable over a region of size . We also plot an arbitrarily scaled profile of , to compare its radial profile with the profile of , for d.

In the launching region at the top of the RZ, there is a location at which at , near to which it follows from the dispersion relation (Eq. 4 in BO10) that the radial wavenumber of gravity waves vanishes. Near this turning point we can approximate the solution in this region by assuming a functional form for . GD98 take , in which case the problem in the launching region reduces to the solution of Airy’s differential equation for . In this case GD98 derive their Eq. 13 for the resulting ingoing energy flux (Eq. 92 in BO10, which we denote herein as ). This should match the asymptotic numerical computation of . We find that a slightly better approximation is to take over the launching region555However, no simple physical arguments for such a profile have been found, which is perhaps to be expected since stellar models contain complicated combinations of physics., since this is valid over from the interface, as we illustrate in Fig. 9. The slope of the curve can be obtained through fitting to the profile of in the stellar model. This allows a slightly more accurate calculation of based on the stellar model than is obtained through direct application of (where the gradient is not uniquely defined). Nevertheless, the differences are only a factor of two at most. In Appendix B we calculate the solution in the launching region and the resulting , for any power law profile of with a positive exponent, using the framework of OL04. Using our approximation, the radial extent of the launching region .

For the case , we obtain, for a non-rotating background,


We express , and compute by solving GD98 Eq. 3 in the CZ, using a linear shooting method, with the boundary conditions that . This matches the numerically computed asymptotic value of for the current Sun quite well, to within (see Fig. 8). The remaining discrepancy is due to the slight variation in background density, and the magnitude of the equilibrium tide, over the launching region. We note that the main change in the energy flux from modifying the profile of in the transition region is to change its frequency dependence.

As can be seen from Fig. 8, the numerically computed asymptotic value of differs from by a factor , even in the case of short-period forcing. We would expect to be correct to within a factor even if the slope varies by an order of magnitude, since it is only raised to the power . However, our slight refinement improves this estimate by for d.

We have performed numerical integrations for several different forcing frequencies (planetary orbital periods), for which we find for fixed stellar and planetary properties (other than the orbital period). If we take into account the fact that , and then consider a fixed tidal potential, we find that , which is slightly different from the power law dependence in Eq. 51. The discrepancy most likely results from the variation in the background density, and the magnitude of the equilibrium tide (which forces the dynamical tide), within this region.

For , which is the relevant regime for which this process is potentially important for the survival of close-in planets (BO10), there is a few per-cent variation in background parameters over a lengthscale . This means that our calculation of differs from the numerically computed value of , by an amount that increases as is made smaller to a maximum of when d. Nevertheless, this discrepancy is small, and our semi-analytic estimate is a good approximation to for all cases that we have modelled. This is used to provide an estimate of , and thus the tidal torque, in § 7.

6.4 Variation between different stellar models

Eq. 51, matches the numerically computed value to within a few per-cent for a variety of solar-type stars. This is because it generally arises that , with , when the launching region is a few percent of the stellar radius. This occurs for the waves excited by planets in short-period orbits. We have confirmed this by computing in a number of stellar models with masses in the range , and ages that represent the range of main-sequence ages expected for these stars. These were computed using ASTEC (Christensen-Dalsgaard 2008).

The main uncertainty in these models is the profile of within the launching region, especially since the slope within of the interface is not well constrained by theory or observations. Helioseismic observations are not yet able to constrain the stratification within this region, owing to the lack of observed g-modes (Ellis 1984; Appourchaux et al. 2010). In addition, the relevant physics included in the stellar models is also uncertain, particularly as a result of changes to the compositional gradient from convective overshoot. The inclusion of helium settling tends to make the interface profile sharper, and the inclusion of turbulent diffusion, that results from convective overshoot, and often parameterised using a simplified 1D model of this process, tends to smooth out the profile near the interface (Jørgen Christensen-Dalsgaard, private communication). However, these changes are small and occur only within a region smaller than the launching region for d. This means that uncertainties in the observations, and the physics, at the interface are unlikely to significantly change for d, which are the planets whose survival could be threatened by our mechanism (as we discuss below). In addition, is only raised to the -2/5 power in our model.

Figure 10: In the top panel we plot in a star at Gyr, with our fitted solutions in the bottom panel. This model has a “bump” in the profile, but this occurs over a region smaller than , and so does not reduce the accuracy of our analytic model (from Eq. 71) in the launching region (or the corresponding energy flux), shown in the bottom panel.

For Model S, the compositional gradient () is unimportant in the launching region, and primarily results from temperature gradients. As the star evolves, the settling of elements heavier than , produces a compositional gradient at the top of the RZ, which can produce “bumps” in the profile of . In Fig. 10 we plot an example of the profile of for a model (with ) at Gyr, together with the best fit solution in the launching region. For d, is generally much larger than these “bumps”, so the wave launching process does not notice such departures from a smooth stratification profile, and Eq. 51 remains a good approximation for the energy flux. However, if the frequency is sufficiently low ( d), the radial extent of the launching region can become comparable with the size of these “bumps”, and the numerical solution can depart appreciably from our analytical model. Planets in such orbits are very unlikely to be affected by stellar tides at such orbital distances, so we do not consider such effects worthy of further consideration. To summarise, we have confirmed that Eq. 51 is a good approximation for the energy flux for planets on orbits of a few days throughout the range of solar-type stars in our study.

6.5 Are 1D linear hydrodynamic calculations reasonable?

The calculations done in § 6 are performed under the assumption of linearity. We have demonstrated in § 5 and in BO10, that this is not a valid assumption near the centre of a star. A simple estimate of the nonlinearity in IGWs in the launching region, for want of a better measure, compares the radial displacement to the radial wavelength, for which , for a hot Jupiter orbiting the Sun on a one-day orbit (this estimate can be obtained from Fig. 7). This value increases linearly with the mass of the planet, but we are still in the linear regime even if we have a close-binary perturber, which indicates that linearity is likely to be a good approximation in the launching region (and throughout the RZ, except for the central regions).

However, it remains to be seen whether for the 1D calculations is the same as that in 2D or 3D numerical simulations of realistic tidal forcing in a model including both a CZ and a RZ. Such simulations as Rogers et al. (2006) could be performed of the whole star subject to tidal forcing. The turbulent convection in these simulations will produce a spectrum of waves at the top of the CZ, in addition to those excited by tidal forcing. It may be that the interaction of these waves reduces the ingoing energy flux. Alternatively, the profile of in the transition region could be modified by realistic modelling of convective overshoot, altering the strength of the stratification within a few percent of a pressure scale height from the interface. This would affect if the overshoot region is comparable with the size of the launching region, though probably not significantly, as we discussed in the previous section.

A toroidal magnetic field in the launching region could also affect the amplitudes of these waves, and the value of . Rogers & MacGregor (2010) find that when the IGW frequency is approximately equal to the Alfvén frequency (), strong wave reflection occurs. This behaviour follows from the dispersion relation for IGWs in the presence of a magnetic field (e.g. Kumar et al. 1999), and could have important consequences for , if the magnetic field is sufficiently strong. However, for in the launching region, we require , when d, which is close to being ruled out on empirical grounds, from measurements of the solar oblateness (Friedland & Gruzinov 2004). If , little attenuation of wave energy over the non-magnetic case is found by Rogers & MacGregor (2010), so this seems unlikely to affect in our case.

In the innermost wavelength, a strong magnetic field would be able to reflect IGWs before they reach the centre, if in this region. This requires a toroidal field of strength , or a poloidal (radial) field of strength . Since it is unlikely that such fields could exist in the RZ (Friedland & Gruzinov, 2004), a magnetic field will probably not affect the reflection of IGWs (excited by short-period planets) in the RZ. It is therefore appropriate to ask whether the waves will break on reaching the centre.

7 Tidal and the breaking criterion

The modified tidal quality factor of the star that results from critical layer absorption can be obtained from Eq. 51 (a slight modification from BO10, resulting from the adopted profile of in the launching region), evaluating to


with the exact value depending on the stellar model adopted. However, this value is found to vary by only a factor of 5 throughout the range of main-sequence stars in our study, for a given orbit. As the star evolves on the main-sequence, the position of the interface moves inwards towards higher density material, which slightly increases . This means that tends to decrease with main-sequence age, though only by up to a factor of 5. The dissipation that results from wave absorption at a critical layer thus becomes more effective as the star evolves.

The waves break if


where is a coefficient666Which was written as in BO10 Eq. 103; see that paper for further details. which is unity for the current Sun but varies with the stellar model, and . This expression is valid if we are sufficiently far from resonance with a global mode. Note, however, that if we are close to a resonance, the central amplitude may be much larger than these estimates would predict, which would make breaking more likely. Note also that this criterion becomes much easier to satisfy as the star evolves, since increases with main-sequence age (see Fig. 1 in BO10). If a planet exceeds this criterion, then its orbital migration is rapid (; see Eq. 105 in BO10) once critical layer absorption occurs, and would therefore threaten the survival of any short-period planet that causes wave breaking with d. Very massive perturbers, with larger orbital moments of inertia than the spin moment of inertia of the RZ of the star, may be able to synchonise the spin of the RZ with the orbit, and prevent further inward migration. This may be an explanation for the synchronisation of close-binary stars (e.g. Mazeh 2008).

7.1 The survival of short-period planets

Our most important result is that this mechanism can potentially explain the survival of all short-period extrasolar planets around solar-type stars777Contained in the catalogue at http://www.exoplanet.hanno-rein.de/ or http://exoplanet.eu/catalog.php with masses in the range . From using the closest fit stellar model to each of these stars, we find that no planet clearly satisfies Eq. 53 with a period d (planets much further out may satisfy the criterion, though this tidal effect is then unimportant). All of these planets are insufficently massive (or orbit sufficiently young stars, with low values of ) that they are unlikely to cause wave breaking at the centre of their hosts. The dominant mechanism of tidal dissipation in these stars is therefore likely to be damping of the equilibrium tide by turbulent convection. This is likely to be relatively inefficient for these periods, since the turbulent viscosity must be reduced when the orbital period is shorter than the convective timescale (see e.g. Zahn 1966; Goldreich & Nicholson 1977; Goodman & Oh 1997; Penev et al. 2007). This could be an important explanation for the survival of these planets.

This mechanism can also explain why the most massive short-period planets (which still have smaller or comparable moments of inertia in the orbit as the RZ of the star), such as WASP-18 b (Hellier et al., 2009), WASP-14 b (Joshi et al., 2009) or CoRoT-14 b, are exclusively found around F-stars, which have convective cores, and in which critical layer formation induced by wave breaking at the centre is unable to operate. In addition, these stars have been found to have weaker tidal dissipation in their outer CZs than G-stars (Barker & Ogilvie, 2009). Note, however, that planets with much larger masses may have sufficient orbital angular momentum to be able to synchronise their stars and reach an approximate tidal equilibrium state, neglecting stellar magnetic braking, which would prevent orbital decay (e.g. Levrard et al. 2009), even for planets in orbit around a solar-type star. Consequently, we make the prediction that fewer massive planets in the range , in orbits with d, that satisfy Eq. 53, will be found around solar-type stars.

8 Conclusions

In this paper, we have presented the results of a set of 3D simulations of IGWs approaching the centre of a solar-type star. The aim of this work is to determine the importance of a nonlinear mechanism of tidal dissipation in solar-type stars, first proposed by GD98, continuing the investigation of BO10. We have confirmed that the main results of BO10 are unaffected by the extension to 3D by performing numerical hydrodynamical simulations, using a Boussinesq-type model solved with a pseudospectral method.

We first derived a linear wave solution in 3D, and found that nonlinearities do not vanish for this wave, unlike the 2D solution, which is exact. Nevertheless, these waves are found to reflect approximately perfectly for moderate-amplitudes, a result which we have qualitatively confirmed in numerical simulations of moderate-amplitude tidal forcing.

The general picture for high-amplitude forcing is that IGWs break within the innermost wavelengths of a star, if they reach the centre with sufficient amplitude to overturn the stratification. If this occurs, they form a critical layer, which we have confirmed from the simulations efficiently absorbs ingoing wave angular momentum. This results in the star being spun up to the orbital angular frequency of the planet, from the inside out. This could be very important to the survival of massive planets in short-period orbits around solar-type stars.

One noticeable difference in 3D is that the absorption of ingoing waves results in the formation of latitudinal differential rotation. This is perpetually reinforced by critical layer absorption. Instabilities may act on this rotation profile, which could homogenise the horizontal angular momentum distribution. These include shear instabilities, which can be linear (Watson, 1981) or nonlinear instabilities, that set in at a critical Reynolds number (Richard & Zahn, 1999). These have growth times comparable to the tidal period, and could transfer angular momentum latitudinally. There are also doubly diffusive instabilities (Goldreich & Schubert 1967; Knobloch & Spruit 1982), or magnetic instabilities, such as the magnetorotational instability (Balbus & Hawley, 1994). However, these mechanisms are unlikely to be able to prevent the critical layer absorption, and thus prevent the tidal engulfment of a short-period planet.

In these simulations we artificially forced waves at the top of the computational domain. In reality, the waves are launched at the top of the RZ. In this paper we also solved the linearised equations for the adiabatic tidal response, in order to study the influence of the launching region in more detail. We provide an expression for the energy flux Eq. 51, which has been checked to agree very well with the numerically computed energy flux for all orbits for which this effect is potentially important, and in a wide range of solar-type stellar models. We then presented the breaking criterion that must be satisfied for the waves to break, slightly refining BO10.

Finally, in § 7 we presented the tidal that results from this process, together with a breaking criterion that determines when these waves break. This allowed us to outline a possible explanation for the survival of all short-period extrasolar planets, that is not in conflict with current observations, and makes predictions which will be tested by ongoing observational studies such as WASP and Kepler.

As in BO10, our simulations do not show any instabilities that act on the waves when they have insufficient amplitude to overturn the stratification. This result suggests that to obtain any efficient tidal dissipation, the waves must have sufficient amplitude to satisfy Eq. 53, so that they overturn the stratification. However, it may be that weaker parametric instabilities operate for waves with lower amplitudes (suggested by GD98). A detailed stability analysis of the 2D exact wave solution written down in BO10 is ongoing, and should shed light on this matter.

In our explanation in the previous section, we ignored the potential effects of the evolution of either the stellar eigenfrequencies or the tidal frequency, that could result in a passage through a resonance with a global mode of oscillation. If this occurs, and simple estimates suggest that this is very likely at least once in the lifetime of a system, then the waves may break, even if the amplitude of the tide is small before the passage through resonance. Future work is required to study this problem in detail.

The simulations performed in this paper used a Boussinesq-type model, which is only valid in the central by radius of a solar-type star. For future work, it is possible to perform simulations using an anelastic code, which could take into account the radial variation in stellar properties throughout the whole RZ (e.g. Rogers & Glatzmaier 2005). We have also neglected the rotation (or possible differential rotation) of the RZ, which is reasonable since we are considering the waves excited by short-period planets in slowly rotating stars. However, the inclusion of rotation (and differential rotation) would certainly be interesting to study for more rapidly rotating stars. One consequence could be that the latitudinal differential rotation observed in the simulations in 3D may be enhanced by rotation, because inertia-gravity waves can be confined in the equatorial regions.

In addition, we have so far considered a planet on a circular, coplanar orbit. Since many planets have been observed with eccentric or inclined orbits, this makes a study of the circularisation and potential spin-orbit alignment through wave breaking and the resulting critical layer absorption, seem an attractive avenue of future research.


I would like to thank STFC for a research studentship, Gordon Ogilvie for his insight and for many useful discussions, Jørgen Christensen-Dalsgaard for providing a set of stellar models, and finally the referee for several helpful suggestions which have improved the readability of the manuscript.

Appendix A IW/OW decomposition

In § 3 we derived an analytic solution for the waves near the centre. We can use this to deconstruct the numerical solution into a single IW and OW, as in BO10. Doing this enables us to quantify the amount of angular momentum absorbed as these waves approach and/or reflect from the centre. However, the reduction in resolution per dimension required for the transition to three dimensions means that an alternative method is preferable, which includes the information at several points, to reduce numerical error. The approach we use is now described.

At every point with spatial coordinates , we have four pieces of information, namely, , , and . Thus, it is possible to compute the (complex) IW and OW amplitudes and for the wave at a single point. This is, however, computationally expensive, since the routines for computing Bessel functions are relatively slow (typically several hundred times slower than a square root). In addition, this method would be subject to potentially significant round-off errors.

We fit the simulation output to a linear model, corresponding to a an IW and an OW. This problem is an overdetermined system of linear equations in terms of the unknown wave amplitudes, which we write as


where is a vector of data variables at each grid point, of size , where is the total number of grid points. is the number of spherical harmonics for which we compute the wave amplitudes, which is usually taken to be one. is a vector of size , whose components are the (complex) IW/OW amplitudes (for each and value), whose values are to be determined. The matrix contains the IW and OW radial functions and spherical harmonic functions, evaluated at the selected grid points using the linear solution in § 3, and has size .

We use a method of least squares to fit our model to the data (Press et al., 1992). This finds the best fit between the linear model data and simulation data, and computes the solution for which the sum of the squared residuals is its least value. This problem has a unique solution if the columns of the matrix are linearly independent, which is true in our case, since , . The solution is found by solving the normal equations


We take into account radial variations in the amplitudes by splitting up the region inside the forcing region into a set of concentric spherical shells of thickness , after removing an inner region of a few grid cells. This approach assumes the solution is locally independent of , hence we can ignore radial derivatives of the amplitudes within each shell. This is not valid in regions where the solution varies rapidly. In addition, we speed up computation by stepping over the grid points in each Cartesian direction, by factors , chosen to take values between 1 and 10. We always ensure that sufficient grid points are available in each shell to accurately compute the amplitudes.

A Matlab routine which reads in SNOOPY data and solves the linear least squares problem has been written, and tested with various combinations of analytic IW/OW solutions. We compute the reflection coefficient . For perfect standing waves, , and . If the IW is entirely absorbed at the centre, then .

The main disadvantage of our approach is that components also contribute to the amplitudes. This problem can be ignored if we trust the computed values of only where the solution is well described by the linear solution, i.e., far from the wave breaking and forcing regions. We therefore reconstruct , and plot it along the -axis, along the line , along the line , and along the line . All components of the solution are not zero for all radii along these lines when (this requires shifting the phase of the linear model solution when the wave phase is nonzero in the simulation data). The simulation variables and IW/OW solutions are plotted using this method in Figs. 1 and 4. Comparing the reconstructed solution with the simulation data allows us to check whether the deconstruction has worked.

Appendix B Analytic calculation of in the launching region

In this section, we adopt the notation of OL04, to avoid reproducing the results in that paper, since we are simply extending the results of their §4.4. In the launching region, we want to match the solutions of the RZ and CZ. Near the boundary , we assume , where , for the moment left unspecified. The characteristic radial extent of the transition region is of order , where , and is the ratio of the spin frequency to the dynamical frequency, which is much smaller than unity for slowly rotating bodies. We write


where is an inner variable, which is of order unity in the launching region. We can write


where . If we pose the perturbation expansions


then we obtain the linearised system of equations 111-116 from OL04, except that we replace by in the buoyancy equation (Eq. 115). A linearised equation for the modified pressure () in the launching region can now be derived, for which


where all coefficients are evaluated at . This equation can be solved using the separation of variables


The operator contains only angular derivatives, and , where