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 Effects section, the tidal gravitational
potential (19) of an orbiting companion can be
expressed as a superposition of partial potentials
\(\PhiTlmk\). For a given &tide
namelist group appearing
in the namelist input file, gyre_tides solves for the
response of the star to each term in the superposition.
Truncating the Sums
Although the sums appearing in eqn. (19) are formally
infinite, the terms with large harmonic degree \(\ell\) and/or
orbital harmonic \(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 (if desired, minimum values can also be
set using the corresponding l_min
and k_min
parameters).
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
for a given tidal partial potential \(\tPhiTlmk\) (see eqn. 21) to be included in calculations; if \(|\yT|\) 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
associated with a specific partial tidal 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)