Imports
import numpy as np
import matplotlib
from matplotlib import animation
from matplotlib import pyplot as plt
import scipy
import proplot as pplt
import psdist.visualization as psvCoupled parametric oscillators
Equations of motion
A previous post examined the analytic solutions to the equation of motion describing a parametric oscillator — an oscillator whose physical properties are time-dependent. This problem was motivated by describing the transverse oscillations of a charged particle in an accelerator. In this post, the treatment will be extended to a coupled parametric oscillator. The oscillator obeys the following equation of motion:
\[ \begin{aligned} x'' + k_{11}(s)x + k_{13}(s)y + k_{14}(s)y' &= 0, \\ y'' + k_{33}(s)y + k_{31}(s)x + k_{32}(s)x' &= 0, \end{aligned} \tag{1}\]
where the prime denotes differentiation with respect to \(s\). We will also assume that \(k_{ij}(s + L) = k_{ij}(s)\) for some \(L\).
The first possible source of coupling is the longitudinal magnetic field produced within a solenoid magnet.
If we assume that the solenoid is very long, the field within the coils points in the longitudinal direction and is approximately constant (\(\mathbf{B}_{sol} = B_0\hat{s}\)). Plugging this into the Lorentz force equation, we find:
\[ \dot{\mathbf{v}} = \frac{q}{m} \mathbf{v} \times \mathbf{B} = \frac{qB_0}{m}\left({v_y\hat{x} - v_x\hat{y}}\right). \tag{2}\]
The motion in \(x\) depends on the velocity in \(y\), and vice versa. Coupling can also be produced by transverse magnetic fields. Recall the multipole expansion of the transverse magnetic field \(\mathbf{B} = (B_x, B_y)\):
There will be nonlinear coupling terms (terms proportional to \(x^j y^k\), where \(j,k > 1\)) when \(n > 2\), but we are interested only in linear coupling. This occurs when the skew quadrupole term \(a_2\) is nonzero, which is true when a quadrupole is tilted in the transverse plane. The field couples the motion in one plane to the displacement in the other.
Solution
Let’s review the approach we took in analyzing the 1D parametric oscillator. We wrote the solution in pseudo-harmonic form with an amplitude \(\sqrt{2 J \beta(s)}\) and phase \(\mu(s)\). We then focused on the motion of a particle in \(x\)-\(x'\) phase space after each focusing period. We observed that the particle jumps around an ellipse: the area of the ellipse is constant and proportional to \(2 J\); the dimensions of the ellipse are determined by \(\beta\) and \(\alpha = -\beta' / 2\); the size of the jumps around the ellipse are determined by \(\mu\). We then wrote the symplectic \(2 \times 2\) transfer matrix \(\mathbf{M}\), which connects the initial and final phase space coordinates through one period, as
\[ \mathbf{M} = \mathbf{V}\mathbf{P}\mathbf{V}^{-1}. \tag{3}\]
\(\mathbf{V}^{-1}\), which is a function of \(\alpha\) and \(\beta\), is a symplectic transformation that deforms the ellipse into a circle while preserving its area, and \(\mathbf{P}\) is a rotation in phase space by the phase advance \(\mu\).
This is an elegant way to describe the motion with a minimal set of parameters. The question is: can we do something similar for coupled motion, in which the phase space is four-dimensional, not two-dimensional? To start, let’s track a particle in a lattice with a nonzero skew quadrupole coefficient and plot its phase space coordinates after each period. (All code is in the following collapsed cell.)
Plotting
labels = ["x", "x'", "y", "y'"]
def vector(v, ax=None, origin=(0.0, 0.0), color='black', lw=None, style='->', head_width=0.4, head_length=0.8):
props = dict()
props['arrowstyle'] = '{},head_width={},head_length={}'.format(style, head_width, head_length)
props['shrinkA'] = props['shrinkB'] = 0
props['fc'] = props['ec'] = color
props['lw'] = lw
ax.annotate('', xy=(origin[0] + v[0], origin[1] + v[1]), xytext=origin, arrowprops=props)
return ax
def animate_corner(X, vecs=None, limits=None, vec_kws=None, **kws):
if limits is None:
maxs = 1.4 * np.max(X, axis=0)
maxs[[0, 2]] = np.max(maxs[[0, 2]])
maxs[[1, 3]] = np.max(maxs[[1, 3]])
limits = [(-m, m) for m in maxs]
if vec_kws is None:
vec_kws = dict()
vec_kws.setdefault('head_width', 0.2)
vec_kws.setdefault('head_length', 0.4)
kws.setdefault('marker', '.')
kws.setdefault('mec', 'None')
kws.setdefault('lw', 0.0)
kws.setdefault('color', 'black')
kws.setdefault('ms', 5.0)
grid = psv.CornerGrid(figwidth=4.5, diag=False, limits=limits, labels=labels)
lines = [[], [], []]
old_lines = [[], [], []]
for i in range(3):
for j in range(i + 1):
ax = grid.axs[i, j]
old_line, = ax.plot([], [], alpha=0.25, **kws)
old_lines[i].append(old_line)
line, = ax.plot([], [], **kws)
lines[i].append(line)
plt.close()
def update(frame):
for ax in grid.axs:
for annotation in ax.texts:
annotation.set_visible(False)
_X = X[:frame, :]
for i in range(3):
for j in range(i + 1):
ax = grid.axs[i, j]
old_lines[i][j].set_data(X[:frame, j], X[:frame, i + 1])
lines[i][j].set_data(X[frame, j], X[frame, i + 1])
if vecs is not None:
v1 = vecs[0][frame]
v2 = vecs[1][frame]
ax = vector(v1[[j, i+1]], ax=ax, origin=(0, 0), color='blue6', **vec_kws)
ax = vector(v2[[j, i+1]], ax=ax, origin=(v1[[j, i+1]]), color='red6', **vec_kws)
grid.axs[0, 1].annotate(
'Period {}'.format(frame), xy=(0.5, 0.5),
xycoords='axes fraction', horizontalalignment='center',
)
return animation.FuncAnimation(grid.fig, update, frames=X.shape[0])Tracking
def unit_symplectic_matrix(n=2):
"""Construct 2n x 2n unit symplectic matrix.
Each 2 x 2 block is [[0, 1], [-1, 0]]. This assumes our phase space vector
is ordered as [x, x', y, y', ...].
"""
if n % 2 != 0:
raise ValueError("n must be an even integer")
U = np.zeros((n, n))
for i in range(0, n, 2):
U[i : i + 2, i : i + 2] = [[0.0, 1.0], [-1.0, 0.0]]
return U
def normalize(eigvecs):
"""Normalize transfer matrix eigenvectors."""
v1, _, v2, _ = eigvecs.T
U = unit_symplectic_matrix(4)
for i in (0, 2):
v = eigvecs[:, i]
val = np.linalg.multi_dot([np.conj(v), U, v]).imag
if val > 0:
eigvecs[:, i], eigvecs[:, i + 1] = eigvecs[:, i + 1], eigvecs[:, i]
eigvecs[:, i : i + 2] *= np.sqrt(2.0 / np.abs(val))
return eigvecs
def construct_normalization_matrix(eigvecs):
"""Construct normalization matrix from transfer matrix eigenvectors."""
v1, _, v2, _ = normalize(eigvecs).T
V = np.zeros((4, 4))
V[:, 0] = v1.real
V[:, 1] = (1.0j * v1).real
V[:, 2] = v2.real
V[:, 3] = (1.0j * v2).real
return V
def rotate_matrix(M, angle=0.0):
"""Get matrix in x-y rotated frame."""
c = np.cos(angle)
s = np.sin(angle)
R = np.array([[c, 0, s, 0], [0, c, 0, s], [-s, 0, c, 0], [0, -s, 0, c]])
return np.linalg.multi_dot([np.linalg.inv(R), M, R])
class MatrixLattice:
"""Transfer matrix representation of periodic linear system."""
def __init__(self):
self.matrices = [] # [element1, element2, ...]
self.M = None # transfer matrix
self.V = None # normalization matrix
self.eigvals = None # eigenvalues of M
self.eigvecs = None # eigenvalues of M
self.v1 = None # eigenvector 1
self.v2 = None # eigenvector 2
self.eig1 = None # eigenvalue 1
self.eig2 = None # eigenvalue 2
self.params = dict() # lattice parameters
def n_elements(self):
"""Return the number of elements in the lattice."""
return len(self.matrices)
def add(self, mat):
"""Add an element to the end of the lattice."""
self.matrices.append(mat)
self.build()
def rotate(self, phi):
"""Apply transverse rotation to all elements."""
self.M = rotate_matrix(self.M, np.radians(phi))
self.analyze()
def build(self):
"""Compute the one-turn transfer matrix."""
n_elements = self.n_elements()
if n_elements == 0:
return
elif n_elements == 1:
self.M = self.matrices[0]
else:
self.M = np.linalg.multi_dot(list(reversed(self.matrices)))
def analyze(self):
"""Compute the lattice parameters."""
self.eigvals, self.eigvecs_raw = np.linalg.eig(self.M)
self.eig1, self.eig2 = self.eigvals[[0, 2]]
self.v1_raw, self.v2_raw = self.eigvecs_raw[:, [0, 2]].T
self.V = construct_normalization_matrix(self.eigvecs_raw.copy())
self.eigvecs = normalize(self.eigvecs_raw)
self.v1, self.v2 = self.eigvecs[:, [0, 2]].T
self.Vinv = np.linalg.inv(self.V)
self.analyze_2D()
self.analyze_4D()
def analyze_2D(self):
"""Compute the 2D Twiss parameters."""
def analyze_2x2(M):
cos_phi = (M[0, 0] + M[1, 1]) / 2.0
sign = 1.0
if abs(M[0, 1]) != 0:
sign = M[0, 1] / abs(M[0, 1])
sin_phi = sign * np.sqrt(1.0 - cos_phi**2)
beta = M[0, 1] / sin_phi
alpha = (M[0, 0] - M[1, 1]) / (2.0 * sin_phi)
tune = np.arccos(cos_phi) / (2.0 * np.pi) * sign
return alpha, beta, tune
alpha_x, beta_x, tune_x = analyze_2x2(self.M[:2, :2])
alpha_y, beta_y, tune_y = analyze_2x2(self.M[2:, 2:])
self.params['alpha_x'] = alpha_x
self.params['alpha_y'] = alpha_y
self.params['beta_x'] = beta_x
self.params['beta_y'] = beta_y
self.params['tune_x'] = tune_x
self.params['tune_y'] = tune_y
def analyze_4D(self):
"""Compute the 4D Twiss parameters."""
V = self.V
beta_1x = V[0, 0] ** 2
beta_2y = V[2, 2] ** 2
alpha_1x = -np.sqrt(beta_1x) * V[1, 0]
alpha_2y = -np.sqrt(beta_2y) * V[3, 2]
u = 1.0 - V[0, 0] * V[1, 1]
nu1 = np.arctan2(-V[2, 1], V[2, 0])
nu2 = np.arctan2(-V[0, 3], V[0, 2])
beta_1y = (V[2, 0] / np.cos(nu1)) ** 2
beta_2x = (V[0, 2] / np.cos(nu2)) ** 2
alpha_1y = (u * np.sin(nu1) - V[3, 0] * np.sqrt(beta_1y)) / np.cos(nu1)
alpha_2x = (u * np.sin(nu2) - V[1, 2] * np.sqrt(beta_2x)) / np.cos(nu2)
self.params['alpha_1x'] = alpha_1x
self.params['alpha_1y'] = alpha_1y
self.params['alpha_2x'] = alpha_2x
self.params['alpha_2y'] = alpha_2y
self.params['beta_1x'] = beta_1x
self.params['beta_1y'] = beta_1y
self.params['beta_2x'] = beta_2x
self.params['beta_2y'] = beta_2y
self.params['u'] = u
self.params['nu1'] = nu1
self.params['nu2'] = nu2
self.params['tune_1'] = np.arccos(self.eig1.real) / (2.0 * np.pi)
self.params['tune_2'] = np.arccos(self.eig2.real) / (2.0 * np.pi)
def normal_form(self):
"""Return normal form of lattice transfer matrix."""
return np.linalg.multi_dot([np.linalg.inv(self.V), self.M, self.V])
def track_part(self, x, nturns=1, norm_coords=False, element_wise=False):
"""Track a single particle."""
if norm_coords:
x = np.matmul(np.linalg.inv(self.V), x)
M_oneturn = self.normal_form()
else:
M_oneturn = self.M
X = [x]
for _ in range(nturns):
if element_wise:
for M in self.matrices:
X.append(np.matmul(M, X[-1]))
else:
X.append(np.matmul(M_oneturn, X[-1]))
return np.array(X)
def fodo(k1, k2, length, fill_fac=0.5, quad_tilt=0, start='quad', n_parts=1):
"""Create a FODO lattice (focus, off, defocus, off)."""
def matrix_drift(length):
"""Drift transfer matrix."""
M = np.zeros((4, 4))
M[:2, :2] = M[2:, 2:] = [[1, length], [0, 1]]
return M
def matrix_quad(length, k, kind='qf', tilt=0):
"""Focusing quadrupole transfer matrix."""
k = np.sqrt(np.abs(k))
cos = np.cos(k * length)
sin = np.sin(k * length)
cosh = np.cosh(k * length)
sinh = np.sinh(k * length)
if kind == 'qf':
M = np.array([[cos, sin / k, 0, 0],
[-k*sin, cos, 0, 0],
[0, 0, cosh, sinh / k],
[0, 0, k * sinh, cosh]])
elif kind == 'qd':
M = np.array([[cosh, sinh / k, 0, 0],
[k * sinh, cosh, 0, 0],
[0, 0, cos, sin / k],
[0, 0, -k * sin, cos]])
if tilt:
M = rotate_matrix(M, np.radians(tilt))
return M
length_quad = 0.5 * fill_fac * length
length_quad = length_quad / n_parts
length_drift = 0.5 * (1.0 - fill_fac) * length
length_drift = length_drift / n_parts
lattice = MatrixLattice()
if start == 'quad':
for _ in range(n_parts):
lattice.add(matrix_quad(0.5 * length_quad, k1, 'qf', quad_tilt))
for _ in range(n_parts):
lattice.add(matrix_drift(length_drift))
for _ in range(n_parts):
lattice.add(matrix_quad(length_quad, k2, 'qd', -quad_tilt))
for _ in range(n_parts):
lattice.add(matrix_drift(length_drift))
for _ in range(n_parts):
lattice.add(matrix_quad(0.5 * length_quad, k1, 'qf', quad_tilt))
elif start == 'drift':
for _ in range(n_parts):
lattice.add(matrix_drift(0.5 * length_drift))
for _ in range(n_parts):
lattice.add(matrix_quad(length_quad, k1, 'qf', +quad_tilt))
for _ in range(n_parts):
lattice.add(matrix_drift(length_drift))
for _ in range(n_parts):
lattice.add(matrix_quad(length_quad, k2, 'qd', -quad_tilt))
for _ in range(n_parts):
lattice.add(matrix_drift(0.5 * length_drift))
lattice.analyze()
return lattice
# Generate a FODO lattice. The two quadrupoles are tilted in opposite directions
# by one degree to generate linear x-y coupling.
L = 5.0
k1 = 0.25
k2 = 0.25
lattice = fodo(k1, k2, L, fill_fac=0.5, quad_tilt=-1.0, start='drift')
# Set the initial particle phase space coordinates.
J1 = 0.5 * 40.0 # amplitude of mode 1
J2 = 0.5 * 10.0 # amplitude of mode 2
psi1 = np.radians(0.0) # initial phase of eigenvector 1
psi2 = np.radians(90.0) # initial phase of eigenvector 2
x1 = np.real(np.sqrt(2.0 * J1) * lattice.v1 * np.exp(1.0j * psi1)) # mode 1 contribution
x2 = np.real(np.sqrt(2.0 * J2) * lattice.v2 * np.exp(1.0j * psi2)) # mode 2 contribution
x = x1 + x2
# Track for 1000 turns.
n_turns = 1000
X1 = lattice.track_part(x1, n_turns)
X2 = lattice.track_part(x2, n_turns)
X = lattice.track_part(x, n_turns)
# Animate 45 turns in corner plot.
scale = 1.4
maxs = 1.4 * np.max(X, axis=0)
maxs[[0, 2]] = np.max(maxs[[0, 2]])
maxs[[1, 3]] = np.max(maxs[[1, 3]])
limits = [(-m, m) for m in maxs]
anim = animate_corner(X[:45, :], limits=limits)
The particle traces donut-like shapes in \(x\)-\(x'\) and \(y\)-\(y'\) instead of ellipses. Here are the shapes after 1000 periods.
Code
grid = psv.CornerGrid(diag=False, figwidth=4.5, labels=labels, limits=limits)
grid.plot_cloud(X, kind='scatter', s=0.75)
grid.axs[0, 1].annotate(
'Period 1000', xy=(0.5, 0.5), xycoords='axes fraction',
horizontalalignment='center', verticalalignment='center'
)
plt.close()
There is clearly more than one frequency present.
Code
fig, axs = pplt.subplots(ncols=2, figsize=(8.0, 1.5), sharey=False, spanx=False, space=6)
axs[0].plot(X[:150, 0], color='black', marker='.', ms=4)
axs[0].format(ylabel='x [mm]')
n = X.shape[0]
m = n // 2
f = (1.0 / n) * np.arange(m)
xf = (1.0 / m) * abs(scipy.fft.fft(X[:, 0])[:m])
axs.format(title_kw=dict(fontsize='medium'))
axs[0].format(xlabel='Period number')
axs[1].format(xlabel='Frequency')
axs[1].plot(f[1:], xf[1:], color='black')
axs[:, 1].format(ylabel='FFT amplitude')This is typical of a coupled oscillator. The motion in coupled systems can be understood as the superposition of normal modes, each of which corresponds to a single frequency. For example, consider two masses connected with a spring. There are two possible ways for the masses to oscillate at the same frequency. The first is a breathing mode in which they move in opposite directions, and the second is a sloshing mode in which they move in the same direction. The motion is the sum of these two modes. We will try to do something similar for a coupled parameteric oscillator.
Transfer matrix eigenvectors
If the phase space coordinate vector \(\mathbf{x} = (x, x', y, y')^T\) evolves according to
\[ \mathbf{x} \rightarrow \mathbf{Mx}, \tag{4}\]
where \(\rightarrow\) represents tracking through one period, it can be shown that \(\mathbf{M}\) is symplectic due to the Hamiltonian mechanics of the system. Consider the eigenvectors of \(\mathbf{M}\):
\[ \mathbf{Mv} = e^{-i\mu}\mathbf{v}. \tag{5}\]
The symplecticity condition causes the eigenvalues and eigenvectors come in two complex conjugate pairs; this gives \(\mathbf{v}_1\), \(\mathbf{v}_2\), \(\mu_1\), \(\mu_2\) and their complex conjugates. The seemingly complex motion is simplified when written in terms of the eigenvectors. We can write any coordinate vector as a linear combination of the real and imaginary components of \(\mathbf{v}_1\) and \(\mathbf{v}_2\):
\[ \mathbf{x} = \text{Re} \left\{ \sqrt{2 J_1}\mathbf{v}_1e^{-i\psi_1} + \sqrt{2 J_2}\mathbf{v}_2e^{-i\psi_2} \right\}. \tag{6}\]
\(\text{Re}\{\dots\}\) selects the non-imaginary component of \(\{\dots\}\). We have introduced two invariant amplitudes (\(J_1\) and \(J_2\)) as well as two initial phases (\(\psi_1\) and \(\psi_2\)). Applying the transfer matrix tacks on a phase to each eigenvector. Thus, what we are observing are the 2D projections of the real components of these eigenvectors as they rotate in the complex plane.
\[ \mathbf{Mx} = \text{Re} \left\{ \sqrt{2 J_1}\mathbf{v}_1e^{-i\left(\psi_1 + \mu_1\right)} + \sqrt{2 J_2}\mathbf{v}_2e^{-i(\psi_2 + \mu_2)} \right\}. \tag{7}\]
Let’s replay the animation, but this time draw a blue arrow for \(\mathbf{v}_1\) and a red arrow for \(\mathbf{v}_2\). We’ve chosen \(J_1 = 4 J_2\) and \(\psi_2 - \psi_1 = \pi/2\).
Code
anim = animate_corner(X[:45, :], vecs=[X1, X2], limits=limits)
Much simpler. Each eigenvector simply rotates at its frequency \(\mu_l\). It also explains why the amplitude in the \(x\)-\(x'\) and \(y\)-\(y'\) planes trade back and forth: it is because the projections of the eigenvectors rotate at different frequencies, sometimes aligning and sometimes anti-aligning. Because of this, the previous invariants \(J_{x,y}\) are replaced by \(J_{1,2}\). In the four-dimensional phase space, each eigenvector traces an ellipsoid, and the particle moves on a torus (represented below). The amplitudes determine the inner and outer radii, and the two phases determine the location on the surface.
Eigenvector parameterization
We are now going to introduce a set of parameters for these eigenvectors, and in turn the transfer matrix. We already have two phases, so that leaves 8 parameters. Our strategy is to observe that each eigenvector traces an ellipse in both horizontal (\(x\)-\(x'\)) and vertical (\(y\)-\(y'\)) phase space. Then, we will simply assign an \(\alpha\) function and \(\beta\) function to each of these ellipses. So, for the ellipse traced by \(\mathbf{v}_1\) in the \(x\)-\(x'\) plane, we have \(\beta_{1x}\) and \(\alpha_{1x}\), and then for the second eigenvector we have \(\beta_{2x}\) and \(\alpha_{2x}\). The same thing goes for the vertical dimension with \(x\) replaced by \(y\).
The actual eigenvectors written in terms of the parameters are
\[ \vec{v}_1 = \begin{bmatrix} \sqrt{\beta_{1x}} \\\\ -\frac{\alpha_{1x} + i(1-u)}{\sqrt{\beta_{1x}}} \\\\ \sqrt{\beta_{1y}}e^{i\nu_1} \\\\ -\frac{\alpha_{1y} + iu}{\sqrt{\beta_{1y}}} e^{i\nu_1} \end{bmatrix}, \quad \vec{v}_2 = \begin{bmatrix} \sqrt{\beta_{2x}}e^{i\nu_2} \\\\ -\frac{\alpha_{2x} + iu}{\sqrt{\beta_{2x}}}e^{i\nu_2} \\\\ \sqrt{\beta_{2y}} \\\\ -\frac{\alpha_{2y} + i(1-u)}{\sqrt{\beta_{2y}}} \end{bmatrix} \tag{8}\]
So in addition to the phases \(\mu_1\) and \(\mu_2\) we have \(\alpha_{1x}\), \(\alpha_{2x}\), \(\alpha_{1y}\), \(\alpha_{2y}\), \(\beta_{1x}\), \(\beta_{2x}\), \(\beta_{1y}\), and \(\beta_{2y}\). That’s pretty much it. There are a few other parameters we need to introduce to simplify the notation, but they are not independent. The first is \(u\), which, as noted in the figure, determines the areas of the ellipses in one plane relative to the other. The second and third are \(\nu_1\) and \(\nu_2\), which are phase differences between the \(x\) and \(y\) components of the eigenvectors (in the animation they are either \(0\) or \(\pi\)). I won’t discuss these here. The last thing to note is that the parameters reduce to their 1D definitions when there is no coupling in the lattice. So we would have \(\beta_{1x}, \beta_{2y} \rightarrow \beta_{x}, \beta_{y}\) and \(\beta_{2x}, \beta_{1y} \rightarrow 0\), and similar for \(\alpha\). The invariants and phase advances would also revert back to their original values: \(J_{1,2} \rightarrow J{x,y}\) and \(\mu_{1,2} \rightarrow \mu_{x,y}\).
Floquet transformation
These eigenvectors can also be used to construct a transformation which removes both the variance in the focusing strength and the coupling between the planes, turning the coupled parametric oscillator into an uncoupled harmonic oscillator. In other words, we seek a matrix \(\mathbf{V}\) such that
\[\mathbf{V^{-1} M V} = \mathbf{P} = \begin{bmatrix} \cos{\mu_1} & \sin{\mu_1} & 0 & 0 \\ -\sin{\mu_1} & \cos{\mu_1} & 0 & 0 \\ 0 & 0 & \cos{\mu_2} & \sin{\mu_2} \\ 0 & 0 & -\sin{\mu_2} & \cos{\mu_2} \end{bmatrix} \tag{9}\]
We can do this simply by rewriting Equation 6 as \(\mathbf{x} = \mathbf{V}\mathbf{x}_n\) with
\[ \mathbf{x}_n = \begin{bmatrix} \sqrt{2J_1}\cos{\psi_1} \\\\ -\sqrt{2J_1}\sin{\psi_1} \\\\ \sqrt{2J_2}\cos{\psi_2} \\\\ -\sqrt{2J_2}\sin{\psi_2} \end{bmatrix} \tag{10}\]
\[ \mathbf{V} = \left[ \text{Re}\{\mathbf{v}_1\}, -\text{Im}\{\mathbf{v}_1\}, \text{Re}\{\mathbf{v}_2\}, -\text{Im}\{\mathbf{v}_2\} \right] \]
Let’s observe the motion in these new coordinates \(\mathbf{x}_n\).
Code
x = x1 + x2
X1 = lattice.track_part(x1, n_turns, norm_coords=True)
X2 = lattice.track_part(x2, n_turns, norm_coords=True)
X1[:, 2:] = 0
X2[:, :2] = 0
X = lattice.track_part(x, n_turns, norm_coords=True)
anim = animate_corner(X[:45, :], vecs=[X1, X2])
The motion is uncoupled after this transformation; i.e., particles move in a circle of area \(2J_1\) in the \(x_n\)-\(x_n'\) plane at frequency \(\mu_1\), and in a circle of area \(2J_2\) in the \(y_n\)-\(y_n'\) plane at frequency \(\mu_2\).
Conclusion
The method introduced here describes a coupled parametric oscillator using the minimum number of parameters. Our physical motivation was an accelerator lattice with linear, coupled forces. There is no agreed upon method to do this among accelerator physicists, but this is the method I have used in my research. I’ve left out details which can be found in [Lebedev2010?,Willeke1989?].