Maximizing Beam Convergence with a FODO Quadrupole Lattice
Particle accelerators rely on carefully arranged magnetic fields to keep a beam of charged particles tightly bunched as it travels down the beamline. Left alone, a beam naturally diverges due to the particles’ spread in transverse momentum. Quadrupole magnets counteract this by focusing the beam — but a magnet that focuses in one plane necessarily defocuses in the other, so accelerator physicists arrange alternating focusing and defocusing quadrupoles into a repeating structure known as a FODO cell (Focusing–Off–Defocusing–Off).
In this article we treat the strengths of the quadrupole magnets as free design parameters and solve a concrete optimization problem: find the quadrupole gradients that minimize the beam envelope (i.e., maximize beam convergence) while keeping the lattice dynamically stable. We’ll build the physics from transfer matrices, formulate the optimization, solve it in Python, and visualize the result — including a 3D rendering of the beam envelope itself.
The Physics: Transfer Matrices and Twiss Parameters
A particle’s transverse position and angle at position $s$ along the beamline are described by the vector
$$
\vec{u}(s) = \begin{pmatrix} x(s) \ x’(s) \end{pmatrix}
$$
Each beamline element (drift space, quadrupole) transforms this vector linearly:
$$
\vec{u}(s_1) = M(s_0 \to s_1), \vec{u}(s_0)
$$
For a drift of length $L$:
$$
M_{\text{drift}} = \begin{pmatrix} 1 & L \ 0 & 1 \end{pmatrix}
$$
For a quadrupole of length $l$ with focusing strength $k$ (in the plane where it focuses, $k>0$):
$$
M_{\text{focus}} = \begin{pmatrix} \cos(\sqrt{k},l) & \dfrac{1}{\sqrt{k}}\sin(\sqrt{k},l) \[6pt] -\sqrt{k}\sin(\sqrt{k},l) & \cos(\sqrt{k},l) \end{pmatrix}
$$
and in the plane where the same magnet defocuses:
$$
M_{\text{defocus}} = \begin{pmatrix} \cosh(\sqrt{k},l) & \dfrac{1}{\sqrt{k}}\sinh(\sqrt{k},l) \[6pt] \sqrt{k}\sinh(\sqrt{k},l) & \cosh(\sqrt{k},l) \end{pmatrix}
$$
Chaining these matrices through one FODO cell (QF – drift – QD – drift) gives a one-cell matrix $M$. If the lattice is periodic, the beam envelope is described by the Twiss parameters $(\beta, \alpha, \gamma)$, obtained from:
$$
\cos\mu = \frac{\mathrm{Tr}(M)}{2}, \qquad \beta = \frac{M_{12}}{\sin\mu}, \qquad \alpha = \frac{M_{11}-M_{22}}{2\sin\mu}, \qquad \gamma = \frac{1+\alpha^2}{\beta}
$$
The beam’s physical size (RMS) is then
$$
\sigma(s) = \sqrt{\beta(s),\varepsilon}
$$
where $\varepsilon$ is the beam emittance — a conserved quantity set by the source. A smaller $\beta$ means a smaller, more tightly focused beam. The lattice is stable only if $|\mathrm{Tr}(M)| < 2$; otherwise the beam amplitude grows unboundedly turn after turn.
Formulating the Optimization Problem
We treat the two quadrupole gradients $k_1$ (QF) and $k_2$ (QD) as design variables and minimize the total transverse beam envelope:
$$
\min_{k_1,,k_2} ; \sigma_x(k_1,k_2) + \sigma_y(k_1,k_2)
$$
subject to the stability constraints
$$
\left|\frac{\mathrm{Tr}(M_x)}{2}\right| < 1, \qquad \left|\frac{\mathrm{Tr}(M_y)}{2}\right| < 1
$$
This is a small but genuinely nonlinear, non-convex problem: the objective is undefined (or effectively infinite) outside the stable region, so the optimizer must navigate around a “forbidden” zone in $(k_1,k_2)$ space to find the true minimum.
Full Python Implementation
1 | import numpy as np |
Code Walkthrough
Sections 1–2 (Geometry and the quadrupole matrix) define the fixed lengths of the lattice and implement quad_matrix_thick_vec, a thick-lens transfer matrix using the exact trigonometric/hyperbolic solution of the equation of motion inside a quadrupole. Crucially, this function is written entirely with NumPy array operations (np.where, np.sqrt, np.cos/np.cosh), so it works transparently whether k is a single Python float (used during optimization) or a full 2D NumPy grid (used later for the landscape scan) — one implementation, two use cases.
Section 3 (one_turn_matrix) chains drift and quadrupole matrices in physical order — QF, drift $L_1$, QD, drift $L_2$ — by hand-multiplying the 2×2 matrix elements. Writing the multiplication out explicitly (rather than using @) is what makes the function automatically vectorize over grids later.
Section 4 (get_twiss_scalar) applies the standard Twiss formulas to a single $(k_1,k_2)$ point, returning None whenever the trace condition $|\mathrm{Tr}(M)|\geq 2$ signals an unstable lattice, or when the resulting $\beta$ is non-physical.
Section 5 (objective) is the function handed to the optimizer: it computes $\sigma_x+\sigma_y$ for a trial $(k_1,k_2)$ and returns a large penalty (1e6) whenever the configuration is unstable — this is what keeps the optimizer inside the physically allowed region.
Section 6 runs scipy.optimize.minimize with the bounded Nelder-Mead method (a derivative-free simplex algorithm, well suited here since the objective has a hard discontinuity at the stability boundary where gradients don’t exist).
Section 7 (the landscape scan) is the part of this problem that would be slow if implemented naively: evaluating the objective at every point of a 300×300 grid (90,000 configurations) with a Python-level double for loop calling get_twiss_scalar each time would carry significant per-call Python overhead. Instead, compute_beta_grid reuses the same one_turn_matrix function directly on 2D NumPy arrays, so the entire grid is evaluated with a handful of vectorized NumPy operations — no explicit loop over grid points at all. This keeps the scan fast enough to recompute interactively even at much finer resolutions.
Section 8 (propagate_twiss) steps through each lattice element in small sub-increments to trace out $\beta_x(s)$ and $\beta_y(s)$ continuously across three repeated FODO cells, using the standard beta-transport formula $\beta(s) = C^2\beta_0 - 2CS\alpha_0 + S^2\gamma_0$. This loop only touches a few hundred points, so its cost is negligible even though it isn’t vectorized.
Sections 9–10 build the two figures described below.
Running the Optimization
The script prints the optimal quadrupole strengths, the resulting beta functions, the corresponding beam sizes, and the improvement relative to the initial guess.
============================================================ Optimization result ============================================================ Optimal QF strength k1 = 2.33329 [1/m^2] Optimal QD strength k2 = 2.86340 [1/m^2] beta_x = 4.93040 m , beta_y = 0.55757 m sigma_x = 2.22045 mm sigma_y = 0.74670 mm Final objective value: 2.967155e-03 m Objective at initial guess: 7.782808e-03 Improvement over initial guess: 61.88% ============================================================
Visualizing the Results
Figure 1 shows two panels side by side. On the left, $\beta_x(s)$ and $\beta_y(s)$ are traced continuously through three repeated FODO cells for the optimized quadrupole strengths, with the QF and QD magnet locations shaded in blue and red — this is the classic “beta-function” plot used throughout accelerator physics to visualize where the beam is widest and narrowest. On the right, a contour map of $\sigma_x+\sigma_y$ across the full $(k_1,k_2)$ parameter space reveals the shape of the optimization landscape, with unstable configurations left blank and the optimum marked with a red star — this makes visible exactly how the optimizer had to thread its way around the forbidden instability region to reach the true minimum.
Figure 2 renders the beam envelope in 3D: at every position $s$ along the beamline, an ellipse of semi-axes $\sigma_x(s)$ and $\sigma_y(s)$ is drawn, and stacking these ellipses along the beam direction produces a tube-like surface. The alternating pinching and widening of this tube is a direct visualization of the strong-focusing principle — the beam is squeezed tight in one plane exactly where it’s allowed to relax in the other, and vice versa, which is precisely why alternating-gradient focusing works better than any single continuously-focusing magnet could.
Discussion
The optimized lattice sits right at the edge of the stable region but not on it — this is a hallmark of well-designed accelerator optics: pushing the focusing strength as hard as possible without crossing into instability gives the smallest possible beam envelope for a given emittance. The vectorized landscape scan makes this trade-off visually obvious: the minimum of $\sigma_x+\sigma_y$ consistently lies just inside the stability boundary, never deep in the “safe” low-$k$ region where focusing is weak, nor beyond the boundary where the beam blows up entirely.
This same transfer-matrix and Twiss-parameter framework scales directly to real accelerator lattices with dozens or hundreds of magnets; the only difference is the number of elements chained together and the number of free parameters handed to the optimizer.
