Extend single-turn tune estimation to coupled lattices

[austin-hoover]

This PR updates orbit.diagnostics.TeapotTuneAnalysisNode to handle coupled accelerator lattices.

Background

This class estimates the tunes using the Average Phase Advance (APA) method [1] over one turn. The $x-x'$ phase space coordinates are normalized such that the turn-by-turn normalized coordinates $u_x-u_x'$ jump around a circle of area $A_x = 2 \pi J_x = \pi (u_x^2 + u_x'^2)$. Here $J_x$ is the Courant-Snyder (CS) invariant, or "action". The phase $\theta_x$ is defined by

$$ \tan(\theta_x) = u_x' / u_x. $$

The tune $\nu_x$ is estimated from the average phase advance over $N$ turns:

$$ \nu_x = \frac{1}{2 \pi N} \sum_{t=1}^{N} \left( \theta_x^{(t)} - \theta_x^{(t - 1)} \right) $$

This node sets $N = 1$. The vertical pane is treated in the same way. The normalized coordinates $u_x, u_x'$ are defined as [2]:

$$ \begin{align} \tilde{x} &= x - \eta_x \delta_p, \\ \tilde{x}' &= x' - \eta_x' \delta_p, \\ \begin{bmatrix} u_x \\ u_x' \end{bmatrix} &= \sqrt{\frac{1}{\beta}} \begin{bmatrix} 1 & 0 \\\ \alpha_x & \beta_x \\ \end{bmatrix} \begin{bmatrix} \tilde{x} \\ \tilde{x}' \end{bmatrix}, \end{align} $$

where $\delta_p = (p - p_0) / p_0$ is the fractional momentum offset and $\{ \alpha_x, \beta_x, \eta_x \}$ are parameters set by the user (corresponding to Courant-Snyder parameters and dispersion).

Extension to coupled lattices

To extend the analysis to coupled lattices, we can define a $6 \times 6$ normalization matrix $\mathbf{V}^{-1}$ [3]:

$$ \mathbf{u} = \mathbf{V}^{-1} \mathbf{x}, $$

with $\mathbf{x} = [x, x', y, y', z, \delta_p]^T$ and $\mathbf{u} = [u_1, u_1', u_2, u_2', u_3, u_3']^T$. If the matrix is designed correctly, it will convert the one-turn transfer matrix $\mathbf{M}$ to the form:

$$ \mathbf{M} = \mathbf{V} \mathbf{P} \mathbf{V}^{-1}, $$

where the $\mathbf{P}$ is a block-diagonal phase advance matrix:

$$ \begin{aligned} \mathbf{P}(\nu_1, \nu_2, \nu_3) &= \begin{pmatrix} \mathbf{R}(2 \pi \nu_1) & \mathbf{0} & \mathbf{0} \\\ \mathbf{0} & \mathbf{R}(2 \pi \nu_2) & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{R}(2 \pi \nu_3) \end{pmatrix} \\ \mathbf{R}(\theta) &= \begin{pmatrix} +\cos\theta & +\sin\theta \\ -\sin\theta & +\cos\theta \end{pmatrix} \end{aligned} $$

The phase $\phi_k$ in mode $k$ is computed as in the 2D analysis ($\tan(\theta_k) = u_k' / u_k$). The tune $\nu_k$ is again estimated from the average phase advance:

$$ \nu_k = \frac{1}{2 \pi N} \sum_{t=1}^{N} \left( \theta_k^{(t)} - \theta_k^{(t - 1)} \right), $$

And the action is given by:

$$ J_{k} = \frac{u_k^2 + u_k'^2}{2}. $$

Selection of normalization matrix

In a linear lattice, we can compute the $\mathbf{V}$ matrix using the eigenvectors of the transfer matrix.

$$ \mathbf{M} \mathbf{v}_k = \lambda_k \mathbf{v}_k. $$

For an $N$-dimensional system ($2N$-dimensional phase space), the eigenvectors come in $N$ pairs: $\{ \mathbf{v}_k, \mathbf{v}_k^* \}$, where $^*$ means complex conjugate. The $\mathbf{V}$ matrix then takes the form [2]:

$$ \mathbf{V} = [Re(\mathbf{v}_1), -Im(\mathbf{v}_1), Re(\mathbf{v}_2), -Im(\mathbf{v}_2), \dots]. $$

Alternatively, if we have a collection of particles and know that the bunch's covariance matrix $\mathbf{\Sigma} = \langle \mathbf{x} \mathbf{x}^T \rangle$ is "matched" to the transfer matrix ($\mathbf{M}\mathbf{\Sigma}\mathbf{M}^T = \mathbf{\Sigma}$), we can compute the eigenvectors without the transfer matrix by solving:

$$ \mathbf{\Sigma} \mathbf{U} \mathbf{v} = i\varepsilon_k \mathbf{v}. $$

Here, $\varepsilon_k$ are the invariant emittances ("eigenemittances", "intrinsic emittances", "mode emittances"), whose product is the total emittance:

$$ \varepsilon \equiv | \mathbf{\Sigma} |^{1/2} = \prod_{k} \varepsilon_k. $$

The $\mathbf{V}$ matrix is then constructed inform the eigenvectors in the same way as above.

Implementation

• Added $6 \times 6$ normalization matrix to BunchTuneAnalysis.
• Updated the method assignTwiss to set elements of the normalization matrix.
• Added ability to set normalization matrix by hand. There are multiple ways to parameterize this matrix, so I leave it up to the user. Later, specific parameterizations could be added to mirror assignTwiss for coupled systems.
• The node stores data by adding an attribute array called ParticlePhaseAttributes to each particle. I've added methods getData to extract the data from the bunch after tracking. These attributes are also written to a file along with the particle coordinates when bunch.dumpBunch is called. (Column labels are wrong; see Issue Particle phase attribute labels printed in wrong order #78.)
• Added getTunes and getActions method to get only the tune and action data.

Tests

There are examples here. These tests calculate the tunes in a FODO lattice and compare them to those calculated from the one-turn transfer matrix. There is a test for an uncoupled FODO lattice as well as a coupled lattice. The tunes are compared to the values obtained from the transfer matrix eigenvalues $\lambda_k$:

$$ \lambda_k = \exp(-2 \pi i \nu_k). $$

References

[1] https://cds.cern.ch/record/292773/files/p147.pdf
[2] https://arxiv.org/pdf/1207.5526
[3] S. Y. Lee, Accelerator Physics

[austin-hoover]

There is also the possibility of merging this into the existing BunchTuneAnalysis class, since the two classes are using the same procedure to estimate the tunes.

[azukov]

• what should be default behavior? should the default behavior define normalization matrix from bunch or bare lattice?
• should the node automatically calculate one turn matrix? could be a separate feature.

[austin-hoover]

To clarify: currently it is up to the user to set the normalization matrix correctly for their particular case. We could add the following methods to help the user:

• setNormMatrixFromTransferMatrix(M): Set normalization matrix $V$ from transfer matrix $M$ by calculating eigenvectors. The transfer matrix is calculated externally by the user.
• setNormMatrixFromCovMatrix(S): Construct $V$ from covariance matrix $\Sigma$ by symplectic diagonalization. The covariance matrix is calculated externally by the user.

Those implement the two strategies to calculate $V$. Then the following additional methods could be defined:

• setNormMatrixFromLattice(lattice): Same as setNormMatrixFromTransferMatrix, but the transfer matix is calculated internally from the AccLattice object using PyORBIT methods.
• setNormMatrixFromBunch(bunch): Same as setNormMatrixFromCovMatrix, but the covariance matrix is calculated internally from the bunch using PyORBIT methods.
• The current method assignTwiss could be renamed to setNormMatrixFromTwiss2D.
• A flag that tells the node to recalculate $V$ on each turn, or every so often, using one of the two strategies. This would be useful if, for example, the lattice has time-dependent parameters.

[austin-hoover]

Added methods:

• setNormMatrixFromTwiss: sets normalization matrix from 2D lattice/twiss parameters.
• setNormMatrixFromTransferMatrix: sets normalization matrix from one-turn transfer matrix.
• setNormMatrixFromCovMatrix: sets normalization matrix from covariance matrix, assuming it is periodic.
• setNormMatrixFromBunch: calculations bunch covariance matrix and then calls setNormMatrixFromCovMatrix .

There are some tests of these different methods under /examples/Diagnostics/Tunes directory. They compare the average calculated tunes to the exact values from the transfer matrix.

[austin-hoover]

The user is still responsible for calculating the transfer matrix, covariance matrix, or Twiss parameters used to set the normalization matrix. The normalization matrix can be changed after each turn. I added a flag that will reset the particle phases to zero (like at the start of tracking) if the normalization matrix was updated since the last turn.

This is shown in the example test_update_norm_mat.py. The normalization matrix is changed after turns 5 and 10. Here is the output:

BunchTuneAnalysis: Adding particle phase information attribute.
BunchTuneAnalysis: Normalization matrix has been updated. Setting particle phases to zero. Tunes will be accurate after the next `analyzeBunch` call.
turn0 nux=0.50146 nuy=0.49815
turn1 nux=0.17100 nuy=0.16240
turn2 nux=0.17100 nuy=0.16240
turn3 nux=0.17100 nuy=0.16240
turn4 nux=0.17100 nuy=0.16240
turn5 nux=0.17100 nuy=0.16240
BunchTuneAnalysis: Normalization matrix has been updated. Setting particle phases to zero. Tunes will be accurate after the next `analyzeBunch` call.
turn6 nux=0.50174 nuy=0.49998
turn7 nux=0.17102 nuy=0.16242
turn8 nux=0.17098 nuy=0.16240
turn9 nux=0.17101 nuy=0.16240
turn10 nux=0.17101 nuy=0.16242
BunchTuneAnalysis: Normalization matrix has been updated. Setting particle phases to zero. Tunes will be accurate after the next `analyzeBunch` call.
turn11 nux=0.50155 nuy=0.49540
turn12 nux=0.17100 nuy=0.16240
turn13 nux=0.17100 nuy=0.16240
turn14 nux=0.17100 nuy=0.16240
turn15 nux=0.17100 nuy=0.16240
turn16 nux=0.17100 nuy=0.16240
turn17 nux=0.17100 nuy=0.16240
turn18 nux=0.17100 nuy=0.16240
turn19 nux=0.17100 nuy=0.16240

[austin-hoover]

Actually I guess we don't need to change the phases, just need the warning message.