Tidal Forcing

This section discusses how to evaluate the stellar response (fluid displacements and perturbations) to tidal forcing, using the gyre_tides frontend. The response data can be used to calculate the secular rates-of-change of orbital elements, or to synthesize a light curve for a tidally distorted star.

Overview

As discussed in the Tidal Equations chapter, the tidal gravitational potential (14) of an orbiting companion can be expressed as a superposition of partial potentials of differing harmonic degree \(\ell\), azimuthal order \(m\) and Fourier harmonic \(k\). For each &tide namelist group appearing in its namelist input file, gyre_tides solves for the response of the star to these partial potentials with a separate calculation for every combination of \(\{\ell,m,k\}\).

Truncating the Sums

Although the sums appearing in eqn. (14) are formally infinite, the terms with large \(\ell\) and/or \(|k|\) typically produce a negligible response. gyre_tides offers a couple of approaches for truncating the sums by dropping these terms. The simplest is to set limits on the maximum values of the indices, through the l_max, k_min and k_max parameters of the &tide namelist group.

A slightly more sophisticated approach is to set these parameters to large-ish values (say, a few hundred), and then also set one or both of the y_T_thresh_abs and y_T_thresh_rel parameters. These establish a threshold on the magnitude of \(\yT\) (see eqn. 16) for a given partial potential to be included in calculations; if it does not meet this threshold, it is ignored.

Optimizing Grids

During the iterative refinement process used in setting up spatial grids, the refinement criteria are evaluated for every partial tide under consideration. If the co-rotating forcing frequency \(\sigmac\) (eqn. 17) of a specific partial potential is small compared to the dynamical frequency of the star, many levels of refinement will occur. While this is exactly what one wants in oscillation calculations (because low-frequency modes have short spatial wavelengths), it often isn’t necessary in tidal calculations because the response of a star to low-frequency forcing is the long-wavelength equilibrium tide.

One way of preventing over-refinement due to low-frequency partial potentials is to set the omega_c_thresh parameter in the &tide namelist group. This establishes a threshold on the dimensionless frequency \(\omegac \equiv \sigmac \, \sqrt{R^{3}/GM}\); partial potentials with \(|\omegac|\) below this threshold are treated as static tides (\(\omegac=0\)), and are not considered during the iterative refinement process.

An alternative approach is to avoid iterative refinement altogether, instead obtaining the spatial grid from an external file (see the 'FILE' option of the scaffold_src parameter, in the Grid Parameters section). This is the most flexible approach, but creating a grid that will adequately resolve the response to each partial potential requires some care.

Output Files

gyre_tides writes response data to summary and detail files. One detail file is created for each partial potential evaluated, and the summary file collects together global data for all partial potentials across all &tide namelist groups. The id output item can be used to determine which group a given response belongs to.

The following Python code demonstrates how the summary data might be used to evaluate the secular rates-of-change of orbital semi-major axis, eccentricity, and argument of periastron, and the stellar torque. The expression for e_dot mirrors eqn. (23) of Sun et al. (2023), and for J_dot eqn. (25) ibid.

import pygyre as pg
import numpy as np

# Read summary file from gyre_tides

s = pg.read_output('summary.h5')

# Extract the first set of responses

sg = s.group_by('id').groups[0]

Omega_orb = sg['Omega_orb']               
R_a = sg['R_a']
q = sg['q']

eps_T = R_a**3*q
    
l = sg['l']
m = sg['m']
k = sg['k']

cbar = sg['cbar']
        
Fbar = -0.5*sg['eul_Phi_ref']/(cbar*eps_T)

x = sg['x_ref']

Gbar_1 = sg['Gbar_1']
Gbar_2 = sg['Gbar_2']
Gbar_3 = sg['Gbar_3']
Gbar_4 = sg['Gbar_4']

# Evaluate secular rates-of-change

kap = np.empty(len(l))

for i in range(len(kap)):
    if k[i] == 0:
        if m[i] == 0:
            kap[i] = 0.5
        elif m[i] > 0 and m[i] <= l[i]:
            kap[i] = 1.
        else:
            kap[i] = 0.
    elif k[i] > 0:
        kap[i] = 1.
    else:
        kap[i] = 0.

# Semi-major axis (units of R per dynamical timescale)

a_dot = np.sum(4. * Omega_orb * q / R_a * 
    (R_a)**(l+3) * (x)**(l+1) * kap * Fbar.imag * Gbar_2)

# Eccentricity (units of per dynamical timescale)

e_dot = np.sum(4. * Omega_orb * q * 
    (R_a)**(l+3) * (x)**(l+1) * kap * Fbar.imag * Gbar_3)

# Argument of periastron (units of radians per dynamical timescale)

pom_dot = np.sum(4. * Omega_orb * q * 
    (R_a)**(l+3) * (x)**(l+1) * kap * Fbar.real * Gbar_1)

# Angular momentum (units of GM^2/R)

J_dot = np.sun(4. * q**2 * R_a *
    (R_a)**(l+3) * (x)**(l+1) * kap * Fbar.imag * Gbar_4)