9. Symplectic Maps

The celmech.maps module implements a couple of area-preserving maps, StandardMap and EncounterMap.

9.1. Chirikov Standard Map

The Chirikov standard map is a symplectic map that depends on a single parameter, \(K\), defined by the equations

\[\begin{split}p' &= p + K \sin\theta \\ \theta' &= \theta + p'~.\end{split}\]

9.2. Comet Map

The comet map is a symplectic map that approximates the dynamics of particles on highly eccentric orbits perturbed by a planet on a circular orbit. It is described in Hadden & Tremaine (2024).

9.3. Encounter Map

The encounter map is a symplectic map that approximates the dynamics of a pair of closely-spaced planets. The map depends on three parameters , \(\epsilon,y\), and \(J\). The map is defined by the equations

\[\begin{split}x' &= x + \epsilon f(\theta;y) \\ \theta' &= \theta + 2\pi(J-x')\end{split}\]

where

\[f(\theta;y) = -\frac{4\pi}{3}\sum_{k}s_k(y)k\sin(k\theta)\]

and the \(s_k(y)\) are the resonance amplitude functions given by the sk() function.

9.4. API

class celmech.maps.CometMap(m, N, q, max_kmax=32, rtol=0.05, atol=1.49e-08, mod=True)

A class representing the comet map.

The map depends on the pericenter distance to perturber semi-major axis ratio, \(q/a_p\), the comet-perturber semi-major axis ratio, \(a/a_p\), and the perturber-star mass ratio, \(\mu\). The map is defined by the equations:

\[\begin{split}\begin{align} w' &= w + \epsilon f(\theta; q/a_p) \\ \theta' &= \theta + 2\pi\left(N + w'\right) \end{align}\end{split}\]

By default, the map is defined on the cylinder with the \(\theta\) coordinate taken modulo \(2\pi\). The parameter mod_p=True can be set to take the \(p\) coordinate modulo \(2\pi\) as well.

9. Example Usage

To initialize the map:

comet_map = CometMap(m=5e-5, q=4/3, N=20, max_kmax=32, rtol=0.05, atol=1e-8, mod=True)

param m:

Planet-star mass ratio (\(\mu\)).

type m:

float

param q:

Pericenter distance of the comet, measured in units of the perturber semi-major axis (\(q/a_p\)).

type q:

float

param N:

Center the map on an \(N:1\) mean-motion resonance (MMR).

type N:

int

param max_kmax:

Maximum order of Fourier amplitude to include before resorting to an asymptotic approximation.

type max_kmax:

int

param rtol:

Relative tolerance for Fourier amplitude calculations before using the asymptotic formula.

type rtol:

float

param atol:

Absolute tolerance for Fourier amplitude calculations before using the asymptotic formula.

type atol:

float

param mod:

If True, the \(\theta\) coordinate is taken modulo \(2\pi\). Default is True.

type mod:

bool, optional

D_QL()

Compute the quasi-linear estimate for the local diffusion coefficient given by

\[D_\mathrm{QL} = \frac{1}{2}\epsilon^2\sum_{k}k^2C_{k}(\beta)^2\]

Returns:

D_QL

Return type:

float

F(theta)

The ‘potential function’, \(F_\beta(\theta)\) from which the map’s kick function is derived.

Parameters:

theta (float) – Angle varaible at which to evaluate the potential function

Returns:

Value of the potential function

Return type:

float

F_asym(theta)

The asymptotic kick potential,

\[-\frac{1}{2} A \ln (2 (\cosh (\lambda )-\cos (\theta )))\]

Parameters:

theta (float) – Angle parameter

Returns:

value of the kick potential

Return type:

float

action(pt)

Evaluate the action zero-form,

\[\lambda(\theta, w) = 2\pi \left( \frac{w'^2}{2} - \frac{\epsilon}{2\pi} F_\beta(\theta) \right),\]

where \(w' = w - \epsilon \partial_\theta F_\beta(\theta)\). The action zero-form satisfies \(T^*(w d\theta) - w d\theta = d\lambda\), where \(T^*\) is the pullback of the map.

Parameters:

pt (array-like) – The point \((\theta, w)\) at which to evaluate the action.

Returns:

The value of the action zero-form, \(\lambda(\theta, w)\).

Return type:

float

delta_F(theta)

Difference between the kick potential and its asymptotic value.

delta_f(theta)

Difference between the kick fucntion and its asymptotic value.

f(theta)

The kick function of the map, \(-\partial_\theta F_\beta(\theta)\).

Parameters:

theta (float) – Angle varaible at which to evaluate the kick function

Returns:

Value of the kick function

Return type:

float

f_asym(theta)

The asymptotic kick function,

\[-\frac{A(\beta)}{2}\frac{\sin\theta}{\cos\theta - \cosh \lambda(\beta)}\]

Parameters:

theta (float) – Angle parameter

Returns:

value of the kick function

Return type:

float

full_map(pt)

Use version of the map, defined as:

\[\begin{split}\begin{aligned} x' &= x - 2\mu \frac{\partial F_\beta(\theta)}{\partial \theta} \\ \theta' &= \theta + \frac{2\pi}{x'^{3/2}}~. \end{aligned}\end{split}\]

where \(x\) represents the test particle’s inverse semi-major axis.

Parameters:

pt (array-like) – The point \((\theta, x)\) to map.

Returns:

pt1 – Resulting point \((\theta', x')\).

Return type:

array

get_eps_crit(tau=1, kmax=None)

Calculate the critical \(\epsilon\) parameter at which the onset of chaos is predicted based on the resonant optical depth.

Parameters:

• tau (float, optional) – Sets the resonant optical depth for the onset of chaos. The default value is 1.

• kmax (int, optional) – The maximum k value of terms to include in the optical depth calculation. By default, the critical \(\epsilon\) is calculated by including all k values, using the asympototic form of resonance widths for large k to estimate a correction.

Returns:

Critical value of epsilon.

Return type:

float

inv_partial_derivs(x0, Nmax)

Get the partial derivatives of the map up to order Nmax evaluated at point x0.

partial_derivs(x0, Nmax)

Get the partial derivatives of the map up to order Nmax evaluated at point x0.

status()

Print a summary of the current status of the map.

Return type:

None

symmetry_lines()

Return the symmetry lines of the map.

Returns:

Tuple containing three functions that parameterize the symmetry lines of the map.

Return type:

tuple

with_variational(X, dX)

Apply the map along with the tangent map to point plus variationals. In particular,

\[\begin{split}\begin{aligned} (\theta', w') &= T(\theta, w) \\ (\delta \theta',\delta w') &= DT(\theta, w) \cdot (\delta \theta,\delta w) \end{aligned}\end{split}\]

where \(T\) is the usual map and \(DT\) is the Jacobian of the map.

Parameters:

• X (array-like) – The point \(X = (\theta,w)\)

• dX (array-like) – The variational vector \((\delta\theta,\delta w)\)

Returns:

• X’ (array-like) – The new point

• dX’ (array-lke) – The new variationl vector

class celmech.maps.EncounterMap(m, J, y0, Nmax=7, mod=True)

A class representing the encounter map. The map depends on three parameters, \(\epsilon\), \(y\), and \(J\). The map is defined by the equations

\[\begin{split}\begin{align} x' &= x + \epsilon f(\theta; y) \\ \theta' &= \theta + 2\pi(J - x') \end{align}\end{split}\]

By default, the map is defined on the cylinder with the \(\theta\) coordinate taken mod \(2\pi\). The parameter mod_p=True can be set to take the \(p\) coordinate modulo \(2\pi\) as well.

See also the Chirikov standard map.

Parameters:

• m (float) – Planet-star mass ratio.

• y (float) – The eccentricity divided by the orbit-crossing eccentricity.

• J (float) – Center the map on the \(J\):J-1 MMR. For integer J, the map is centered on a first-order MMR. For rational \(J = p/q\), the map is centered on a \(q\)-th order MMR.

• mod_theta (bool, optional) – If True, the \(\theta\) coordinate is taken modulo \(2\pi\). Default is True.

• mod_p (bool, optional) – If True, the \(p\) coordinate is taken modulo \(2\pi\). Default is False.

inv_partial_derivs(x0, Nmax)

Get the partial derivatives of the map up to order Nmax evaluated at point x0.

partial_derivs(x0, Nmax)

Get the partial derivatives of the map up to order Nmax evaluated at point x0.

class celmech.maps.StandardMap(K, mod_theta=True, mod_p=False)

A class representing the Chirikov standard map. The map depends on a single parameter, \(K\) and is defind by

\[\begin{split}\begin{align} p' &=& p + K \sin\theta\\ \theta' &=& \theta + p' \end{align}\end{split}\]

By default, the map is defined on the cylinder with the \(\theta\) coordinate taken mod \(2\pi\). The parameter mod_p=True can be set to take the \(p\) coordinate modulo \(2\pi\) as well.

Parameters:

• K (float) – Map non-linearity parameter.

• mod_theta (bool, optional) – If True, the \(\theta\) coordinate is taken modulo \(2\pi\). Default is True

• mod_p (bool, optional) – If True, the \(p\) coordinate is taken modulo \(2\pi\). Default is False.

action(pt)

Evaluate The action zero-form,

\[\lambda(\theta, w) = 2\pi \left( \frac{w'^2}{2} - \frac{\epsilon}{2\pi} F_\beta(\theta) \right),\]

where \(w' = w - \epsilon \partial_\theta F_\beta(\theta)\). The action zero-form satisfies \(T^*(w d\theta) - w d\theta = d\lambda\), where \(T^*\) is the pullback of the map.

Parameters:

pt (array-like) – The point \((\theta, w)\) at which to evaluate the action.

Returns:

The value of the action zero-form, \(\lambda(\theta, w)\).

Return type:

float

inv(x)

The inverse mapping

\[\begin{split}\begin{align} \theta &=& p' - \theta' \\ p &=& p' - K \sin\theta \end{align}\end{split}\]

Parameters:

x (array-like) – The point \((\theta',p')\)

Returns:

The point \((\theta,p)\)

Return type:

array-like

inv_partial_derivs(x, Nmax)

Get the partial derivatives of the inverse map evaluated at the point x0 up to order Nmax.

Parameters:

• x (array-like) – The point at which derivatives are to be evaluated

• Nmax (int) – Maximum order of the partial derivatives

Returns:

T – The partial derivatives of the map. Writing the value of the map at a point \((x_1,x_2)\) as \(T(x_1,x_2) = (T_1(x_1,x_2),T_2(x_1,x_2))\), the entry T[i,n,m] stores

\[\frac{\partial^{(n+m)}}{\partial x_1^n \partial x_2^m} T_i\]

Note that T[:,0,0] give the value of the map.

Return type:

array, shape (2,Nmax+1,Nmax+1)

jac(x)

Evaluate the Jacobian map at \(x=(\theta,p)\), given by

\[\begin{split}DT(x) = \begin{pmatrix} 1 + K\cos\theta & 1 \\ K\cos\theta & 1 \end{pmatrix}\end{split}\]

9. Aruments

xarray

Point(s) at which to evaluate the Jacobian.

returns:

DT – Value of the Jacobain at point x.

rtype:

array

property mod_theta

Is the coordinate \(\theta\) calculated modulo \(2\pi\)?

partial_derivs(x, Nmax)

Get the partial derivatives of the map evaluated at the point x0 up to order Nmax.

Parameters:

• x (array-like) – The point at which derivatives are to be evaluated

• Nmax (int) – Maximum order of the partial derivatives to return

Returns:

T – The partial derivatives of the map. Writing the value of the map at a point \((x_1, x_2)\) as \(T(x_1, x_2) = (T_1(x_1, x_2), T_2(x_1, x_2))\), the entry T[i, n, m] stores \(\frac{\partial^{(n+m)}}{\partial x_1^n \partial x_2^m} T_i\)

Note that T[:, 0, 0] gives the value of the map.

Return type:

array, shape (2, Nmax+1, Nmax+1)

symmetry_lines()

Return the symmetry lines of the map.

Returns:

Tuple containing three functions that parameterize the symmetry lines of the map.

Return type:

tuple

with_variational(X, dX)

Apply the map along with the tangent map to a point plus variationals.

In particular,

\[\begin{split}\begin{aligned} (\theta', p') &= T(\theta, p) \\ (\delta \theta', \delta p') &= DT(\theta, p) \cdot (\delta \theta, \delta p) \end{aligned}\end{split}\]

where \(T\) is the usual map and \(DT\) is the Jacobian of the map.

Parameters:

• X (array-like) – The point \(X = (\theta, w)\).

• dX (array-like) – The variational vector \((\delta \theta, \delta w)\).

Returns:

• X’ (array-like) – The new point \((\theta', w')\).

• dX’ (array-like) – The new variational vector \((\delta \theta', \delta w')\).