# -*- coding: utf-8 -*-
"""
Robinson (coupled-bunch mode 0) instabilities:
- Coupled/uncoupled Robinson dipole to octupole modes, including fast
mode-coupling, based on Bosch equations [1].
- Robinson instabilites using the Gaussian longitudinal mode-coupling
instability method [2].
References
----------
[1] Bosch, R. A., K. J. Kleman, and J. J. Bisognano. "Robinson instabilities
with a higher-harmonic cavity." Physical Review Special Topics-Accelerators
and Beams 4.7 (2001): 074401.
[2] Alves, Murilo B. "Theoretical models for longitudinal coupled-bunch instabilities
driven by harmonic cavities in electron storage rings." Physical Review Accelerators
and Beams 28.3 (2025): 034401.
"""
from __future__ import annotations
from math import factorial
from typing import TYPE_CHECKING
import numpy as np
from scipy.optimize import root, root_scalar
if TYPE_CHECKING: # pragma: no cover - hints only
from mbtrack2 import CavityResonator, Synchrotron
from pycolleff.longitudinal_equilibrium import LongitudinalEquilibrium
[docs]def omega_n(ring: Synchrotron, cavity: CavityResonator) -> float:
"""Detuned resonant angular frequency of a cavity. Near Eq. (B.3) in [1]."""
return (
2 * cavity.QL * cavity.m * ring.omega1 / (2 * cavity.QL + np.tan(-cavity.psi))
)
[docs]def phi_pm(
ring: Synchrotron, cavity: CavityResonator, Omega: float, sign: str
) -> float:
"""Upper ("+") / lower ("-") sideband phase. Near Eq. (B.4) in [1]."""
wn = omega_n(ring, cavity)
if sign == "+":
return np.arctan(2 * cavity.QL * (cavity.m * ring.omega1 + Omega - wn) / wn)
elif sign == "-":
return np.arctan(2 * cavity.QL * (cavity.m * ring.omega1 - Omega - wn) / wn)
[docs]def A_tilde(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
Omega: float,
) -> float:
"""Cavity transfer term à (dipole). Bosch Appendix B [1]."""
coef = ring.ac * ring.omega1 * I0 / (2 * ring.E0 * ring.T0)
sum_val = 0
for i, cavity in enumerate(cavity_list):
brackets = np.sin(2 * phi_pm(ring, cavity, Omega, "-")) + np.sin(
2 * phi_pm(ring, cavity, Omega, "+")
)
sum_val += cavity.m * cavity.RL * F[i] ** 2 * brackets
return coef * sum_val
[docs]def B_tilde(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
Omega: float,
bunch_length: float,
) -> float:
"""Cavity transfer term B̃ (quadrupole). Bosch Appendix B [1]."""
coef1 = ring.ac * ring.omega1 * I0 / (2 * ring.E0 * ring.T0)
coef2 = (ring.omega1 * bunch_length) ** 2
sum_val = 0
for i, cavity in enumerate(cavity_list):
brackets = np.sin(2 * phi_pm(ring, cavity, Omega, "-")) + np.sin(
2 * phi_pm(ring, cavity, Omega, "+")
)
sum_val += cavity.m**3 * cavity.RL * F[i] ** 2 * brackets
return coef1 * coef2 * sum_val
[docs]def D_tilde(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
Omega: float,
bunch_length: float,
) -> float:
"""Cavity transfer cross-term D̃. Bosch Appendix B [1]."""
coef2 = ring.ac * ring.omega1 * I0 / (2 * ring.E0 * ring.T0)
coef1 = ring.omega1 * bunch_length
sum_val = 0
for i, cavity in enumerate(cavity_list):
brackets = np.sin(2 * phi_pm(ring, cavity, Omega, "-")) - np.sin(
2 * phi_pm(ring, cavity, Omega, "+")
)
sum_val += cavity.m**2 * cavity.RL * F[i] ** 2 * brackets
return coef1 * coef2 * sum_val
[docs]def a_tilde(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
Omega: float,
) -> float:
"""Cavity transfer term ã (dipole damping). Bosch Appendix B [1]."""
coef = ring.ac * ring.omega1 * I0 / (ring.E0 * ring.T0)
sum_val = 0
for i, cavity in enumerate(cavity_list):
brackets = (
np.cos(phi_pm(ring, cavity, Omega, "-")) ** 2
- np.cos(phi_pm(ring, cavity, Omega, "+")) ** 2
)
sum_val += cavity.m * cavity.RL * F[i] ** 2 * brackets
return coef * sum_val + 2 * Omega / ring.tau[2]
[docs]def b_tilde(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
Omega: float,
bunch_length: float,
) -> float:
"""Cavity transfer term b̃ (quadrupole damping). Bosch Appendix B [1]."""
coef1 = ring.ac * ring.omega1 * I0 / (ring.E0 * ring.T0)
coef2 = (ring.omega1 * bunch_length) ** 2
sum_val = 0
for i, cavity in enumerate(cavity_list):
brackets = (
np.cos(phi_pm(ring, cavity, Omega, "-")) ** 2
- np.cos(phi_pm(ring, cavity, Omega, "+")) ** 2
)
sum_val += cavity.m**3 * cavity.RL * F[i] ** 2 * brackets
return coef1 * coef2 * sum_val + 4 * Omega / ring.tau[2]
[docs]def d_tilde(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
Omega: float,
bunch_length: float,
) -> float:
"""Cavity transfer cross-term d̃. Bosch Appendix B [1]."""
coef2 = ring.ac * ring.omega1 * I0 / (ring.E0 * ring.T0)
coef1 = ring.omega1 * bunch_length
sum_val = 0
for i, cavity in enumerate(cavity_list):
brackets = (
np.cos(phi_pm(ring, cavity, Omega, "-")) ** 2
+ np.cos(phi_pm(ring, cavity, Omega, "+")) ** 2
)
sum_val += cavity.m**2 * cavity.RL * F[i] ** 2 * brackets
return coef1 * coef2 * sum_val
[docs]def ar_B11(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
Omega: float,
bunch_length: float,
omega_r: float,
) -> float:
"""Coupled-mode damping rate. Eq. (B.11) in [1]."""
a = a_tilde(ring, cavity_list, I0, F, Omega)
b = b_tilde(ring, cavity_list, I0, F, Omega, bunch_length)
d = d_tilde(ring, cavity_list, I0, F, Omega, bunch_length)
A = A_tilde(ring, cavity_list, I0, F, Omega)
B = B_tilde(ring, cavity_list, I0, F, Omega, bunch_length)
D = D_tilde(ring, cavity_list, I0, F, Omega, bunch_length)
num = (
a * (Omega**2 - (2 * omega_r) ** 2 + B)
+ b * (Omega**2 - omega_r**2 + A)
- 2 * D * d
)
den = 2 * Omega * (2 * Omega**2 - 5 * omega_r**2 + A + B)
return num / den
[docs]def Omega_B13(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
Omega: float,
alpha_r: float,
bunch_length: float,
omega_r: float,
sign: str,
abs_val: bool = True,
return_root: bool = False,
) -> float:
"""Coupled-mode frequency (or its discriminant). Eq. (B.13) in [1]."""
a = a_tilde(ring, cavity_list, I0, F, Omega)
b = b_tilde(ring, cavity_list, I0, F, Omega, bunch_length)
d = d_tilde(ring, cavity_list, I0, F, Omega, bunch_length)
A = A_tilde(ring, cavity_list, I0, F, Omega)
B = B_tilde(ring, cavity_list, I0, F, Omega, bunch_length)
D = D_tilde(ring, cavity_list, I0, F, Omega, bunch_length)
part1 = (5 * omega_r**2 - A - B) / 2
root_val = (
(3 * omega_r**2 + A - B) ** 2 / 4
+ D**2
- d**2
+ (a - 2 * Omega * alpha_r) * (b - 2 * Omega * alpha_r)
)
if return_root:
return root_val
if abs_val:
root_val = np.abs(root_val)
if sign == "+":
return part1 + np.sqrt(root_val)
elif sign == "-":
return part1 - np.sqrt(root_val)
[docs]def robinson_damping_rate(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
Omega: float,
mu: int,
bunch_length: float,
) -> float:
"""Radiation-free Robinson damping rate of mode mu. Eq. (16) in [1]."""
coef1 = 8 * ring.ac * I0 / (ring.E0 * ring.T0)
coef2 = (
mu * (ring.omega1 * bunch_length) ** (2 * mu - 2) / (2**mu * factorial(mu - 1))
)
sum_val = 0
for i, cavity in enumerate(cavity_list):
prod = (
np.tan(-cavity.psi)
* np.cos(phi_pm(ring, cavity, Omega, "+")) ** 2
* np.cos(phi_pm(ring, cavity, Omega, "-")) ** 2
)
sum_val += cavity.m ** (2 * mu - 2) * F[i] ** 2 * cavity.RL * cavity.QL * prod
return coef1 * coef2 * sum_val
[docs]def robinson_no_coupling(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
omega_r: float,
bunch_length: float,
landau_threshold: np.ndarray,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""Robinson modes mu = 1..4 solved independently (no mode coupling).
For each mode the coherent frequency is found from the dispersion relation,
then flagged unstable if its (radiation-included) damping rate is negative and
its complex frequency shift exceeds the Landau threshold.
Parameters
----------
ring : Synchrotron
Ring parameters.
cavity_list : list of CavityResonator
Solved cavities.
I0 : float
Beam current in A.
F : numpy.ndarray
Amplitude form factor per cavity.
omega_r : float
Synchrotron angular frequency.
bunch_length : float
RMS bunch length in [s].
landau_threshold : numpy.ndarray
Per-mode thresholds from
landau_threshold.
Returns
-------
robinson : numpy.ndarray of bool, shape (4,)
Instability flag per mode.
Omega : numpy.ndarray, shape (4,)
Coherent angular frequency per mode.
converged : numpy.ndarray of bool, shape (4,)
Root-finder convergence per mode.
growth : numpy.ndarray, shape (4,)
Net growth rate (radiation included) per mode.
References
----------
Eq. (15)-(16) in [1].
"""
modes = np.array([1, 2, 3, 4])
ac_robinson = np.zeros(len(modes), dtype=bool)
Omega_modes = np.zeros(len(modes), dtype=float)
converged = np.zeros(len(modes), dtype=bool)
growth = np.zeros(len(modes), dtype=float)
for mu in modes:
Omega0 = mu * omega_r
def f(Omega):
coef1 = ring.ac * ring.omega1 * I0 / (ring.E0 * ring.T0)
coef2 = (
mu
* (ring.omega1 * bunch_length) ** (2 * mu - 2)
/ (2**mu * factorial(mu - 1))
)
sum_val = 0
for i, cavity in enumerate(cavity_list):
brackets = np.sin(2 * phi_pm(ring, cavity, Omega, "-")) + np.sin(
2 * phi_pm(ring, cavity, Omega, "+")
)
sum_val += cavity.m ** (2 * mu - 1) * cavity.RL * F[i] ** 2 * brackets
inner = (mu * omega_r) ** 2 - coef1 * coef2 * sum_val
return Omega - np.sqrt(inner)
sol = root_scalar(f, x0=Omega0, xtol=1e-6)
if sol.converged:
Omega = sol.root
Omega_modes[mu - 1] = Omega
converged[mu - 1] = True
alpha_r = robinson_damping_rate(
ring, cavity_list, I0, F, Omega, mu, bunch_length
)
alpha_r_incl_rad_damp = alpha_r + mu / ring.tau[2]
# store the net growth rate (positive => unstable), not the damping rate
growth[mu - 1] = -alpha_r_incl_rad_damp
delta_Omega = Omega - 1j * alpha_r - mu * omega_r
if alpha_r_incl_rad_damp < 0:
if np.abs(delta_Omega) > landau_threshold[mu - 1]:
ac_robinson[mu - 1] = True
else:
converged[mu - 1] = False
return ac_robinson, Omega_modes, converged, growth
[docs]def robinson_coupling(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
omega_r: float,
bunch_length: float,
landau_threshold: np.ndarray,
abs_val: bool = True,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""Coupled Robinson dipole/quadrupole modes with fast mode-coupling.
Solves the coupled dipole and quadrupole dispersion relations (Eq. B.11/B.13),
flags each as Landau-stabilised or not, and detects fast mode-coupling when the
discriminant goes negative.
Parameters
----------
ring : Synchrotron
Ring parameters.
cavity_list : list of CavityResonator
Solved cavities.
I0 : float
Beam current in A.
F : numpy.ndarray
Amplitude form factor per cavity.
omega_r : float
Synchrotron angular frequency.
bunch_length : float
RMS bunch length in [s].
landau_threshold : numpy.ndarray
Per-mode thresholds from
landau_threshold.
abs_val : bool, optional
Take the absolute value of the discriminant when solving (default True).
Returns
-------
robinson : numpy.ndarray of bool, shape (4,)
[coupled_dipole, coupled_quadrupole, fast_coupling_by_dipole,
fast_coupling_by_quadrupole].
Omega : numpy.ndarray, shape (2,)
Coherent dipole and quadrupole angular frequencies.
converged : numpy.ndarray of bool, shape (2,)
Root-finder convergence for dipole and quadrupole.
growth : numpy.ndarray, shape (4,)
Net growth rate per entry of robinson.
References
----------
Appendix B, Eq. (B.11)-(B.13) in [1].
"""
methods = ["hybr", "broyden1"]
coupled_dipole = False
coupled_quadrupole = False
fast_mode_coupling_bydip = False
fast_mode_coupling_byquad = False
Omega_dip = np.nan
Omega_quad = np.nan
growth = np.zeros(4, dtype=float)
def f_dipole(x):
y0 = x[0] - np.sqrt(
Omega_B13(
ring,
cavity_list,
I0,
F,
x[0],
x[1],
bunch_length,
omega_r,
"-",
abs_val,
)
)
y1 = x[1] - ar_B11(ring, cavity_list, I0, F, x[0], bunch_length, omega_r)
return np.array([y0, y1])
x0 = np.array([omega_r, 1 / ring.tau[2]])
for method in methods:
try:
sol = root(f_dipole, x0, tol=1e-6, method=method)
except ValueError:
sol = {"success": False}
if sol["success"]:
break
if sol["success"]:
Omega_dip = sol.x[0]
alpha_r_incl_rad_damp_dip = sol.x[1]
# store the net growth rate (positive => unstable), not the damping rate
growth[0] = -alpha_r_incl_rad_damp_dip
delta_Omega_dip = Omega_dip - 1j * alpha_r_incl_rad_damp_dip - omega_r
converged_dip = True
root_val = Omega_B13(
ring,
cavity_list,
I0,
F,
Omega_dip,
alpha_r_incl_rad_damp_dip,
bunch_length,
omega_r,
"-",
return_root=True,
)
if root_val < 0:
fast_mode_coupling_bydip = True
if alpha_r_incl_rad_damp_dip < 0:
if np.abs(delta_Omega_dip) > landau_threshold[0]:
coupled_dipole = True
else:
converged_dip = False
def f_quadrupole(x):
y0 = x[0] - np.sqrt(
Omega_B13(
ring,
cavity_list,
I0,
F,
x[0],
x[1],
bunch_length,
omega_r,
"+",
abs_val,
)
)
y1 = x[1] - ar_B11(ring, cavity_list, I0, F, x[0], bunch_length, omega_r)
return np.array([y0, y1])
x0 = np.array([2 * omega_r, 2 / ring.tau[2]])
for method in methods:
try:
sol = root(f_quadrupole, x0, tol=1e-6, method=method)
except ValueError:
sol = {"success": False}
if sol["success"]:
break
if sol["success"]:
Omega_quad = sol.x[0]
alpha_r_incl_rad_damp_quad = sol.x[1]
# store the net growth rate (positive => unstable), not the damping rate
growth[1] = -alpha_r_incl_rad_damp_quad
delta_Omega_quad = Omega_quad - 1j * alpha_r_incl_rad_damp_quad - 2 * omega_r
converged_quad = True
root_val = Omega_B13(
ring,
cavity_list,
I0,
F,
Omega_quad,
alpha_r_incl_rad_damp_quad,
bunch_length,
omega_r,
"+",
return_root=True,
)
if root_val < 0:
fast_mode_coupling_byquad = True
if alpha_r_incl_rad_damp_quad < 0:
if np.abs(delta_Omega_quad) > landau_threshold[1]:
coupled_quadrupole = True
else:
converged_quad = False
if fast_mode_coupling_bydip:
growth[2] = growth[0]
if fast_mode_coupling_byquad:
growth[3] = growth[1]
robinson = np.array(
[
coupled_dipole,
coupled_quadrupole,
fast_mode_coupling_bydip,
fast_mode_coupling_byquad,
]
)
Omega = np.array([Omega_dip, Omega_quad])
converged = np.array([converged_dip, converged_quad])
return (robinson, Omega, converged, growth)
[docs]def robinson_glmci(
ring: Synchrotron, LE: LongitudinalEquilibrium, max_azi: int = 2, max_rad: int = 0
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""
Robinson (coupled-bunch mode 0) instabilites using the Gaussian
longitudinal mode-coupling instability method [2].
This method does not distinguish between dipole Robinson, quadrupole Robisonson or
fast-mode coupling instability. A positive result is flagged as fast-mode coupling
instabilty.
Parameters
----------
ring : Synchrotron
Ring parameters.
LE : LongitudinalEquilibrium
Solved pycolleff equilibrium handle.
max_azi : int, optional
Maximum azimuthal mode number (default 2).
max_rad : int, optional
Maximum radial mode number (default 0).
Returns
-------
robinson : numpy.ndarray of bool, shape (4,)
[False, False, FMCI, FMCI] (the fast-coupling dipole/quadrupole slots).
Omega : numpy.ndarray, shape (2,)
Real parts of the first two eigenfrequencies.
converged : numpy.ndarray of bool, shape (2,)
Always [True, True] (the eigensolver does not iterate).
growth : numpy.ndarray, shape (4,)
Net growth rate placed in the fast-coupling slots.
References
----------
Alves [2].
"""
from albums.physics.ptbl import (
prepare_alves_le,
) # local: avoids a robinson<->ptbl import cycle
prepare_alves_le(ring, LE)
eigenfreq, *_ = LE.calc_mode_coupling(
w=[-300 * ring.omega1, +300 * ring.omega1],
cbmode=0,
max_azi=max_azi,
max_rad=max_rad,
use_fokker=False,
delete_m0=True,
delete_m0k0=True,
reduced=True,
)
idx = np.argmax(eigenfreq.imag)
grate = eigenfreq.imag[idx] - 1 / ring.tau[2]
FMCI = grate > 0
robinson = np.array([False, False, FMCI, FMCI])
Omega = np.array([eigenfreq.real[0], eigenfreq.real[1]])
converged = np.array([True, True])
growth = np.array([0.0, 0.0, grate, grate], dtype=float)
return robinson, Omega, converged, growth