# 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)
```