Adjoint approach to beam optics sensitivity based on Hamiltonian particle dynamics

The design and optimization of complex devices and instrumentation is common to all fields of science and engineering. The design of a particular device generally depends on many parameters, and an optimum design must be found in a large, potentially infinite, dimensional parameter space. A powerful tool useful in the design process is a “sensitivity function” that quantifies how small changes in the system parameters affect a particular design metric. The purpose of the present paper is to show that, by exploiting the symplectic property¹ of Hamilton's equations, a sensitivity function can be defined that will aid in the design of charged particle beam sources and focusing systems used in accelerators and in beam-driven sources of coherent radiation such as free-electron lasers, klystrons, and traveling wave tubes. Further, we show that the sensitivity function can be computed by making a slight modification of the algorithm in the particle trajectory codes currently used to simulate these systems.

Normally, one would expect that to calculate a sensitivity function in an N-dimensional parameter space, N + 1 computations would be required to determine the gradient of the metric in that space. Since N, the number of design parameters of interest, is generally very large, the required computational burden can be so large so as to rule out a full design optimization or sensitivity analysis. However, using adjoint techniques,² this computational requirement can be reduced to as few as two computations—a “base case” and a single, especially formulated perturbed case. The latter case is referred to as the “adjoint” problem. The adjoint approach has previously been applied in circuit theory,³ electromagnetics,⁴ aerodynamics,⁵ plasma physics,⁶ and other fields.

The basic adjoint approach may be described as follows: Let the vector $X$ represent the detailed state of the system with design parameters B; the dimension of B is N. For example, X could contain trajectory information for every particle in an electron gun and B could contain information on the locations and potentials of the electrode surfaces or the position and strength of focusing magnets. We assume that small changes in X produced by small changes in B are related by a linear operator $AX0$⁠, which depends on the unperturbed system state $X0$⁠, that is, $AX0δX=δB$⁠. Now, we usually do not need all the information contained in $δX$⁠, that is, complete, detailed knowledge of the change in the system state is generally not required in the design process. Rather, what is important is the change in some particular performance metric(s) or figure(s) of merit $M$ that depends on X. The change in $M$ may be written as $δM=C†δX$⁠, where $C†≡∂M/∂X$ and the dagger symbol (⁠$†$⁠) indicates the adjoint (conjugate transpose) operator. The dimension of M is generally much smaller than that of N; in many cases of interest, the dimension of M is unity.

The key observation that provides the adjoint method its advantages is to note that the change in $M$ can be evaluated without solving for $δX$⁠. Specifically, if we solve the adjoint equation $A†X0δY=C$ for $δY$⁠, then it follows that $δM=C†⋅δX=δY†⋅δB$⁠. Thus, a single inversion of the adjoint equation to find the sensitivity function $δY†$ provides all the information needed to evaluate the change in the metric $M$ for arbitrary changes in the parameters $δB$⁠. We never need to compute $δX$⁠.

In this paper, we develop a formalism using the adjoint method to characterize the sensitivity of electron beam parameters at the exit of an electron gun to variations in the gun design parameters such as the electrode shape, position, and voltage and magnetic field distribution and alignment. A 3D rendering of such an electron gun is shown in Fig. 1. While the example considered here is that of an electron gun, the method can be extended to other beam optics structures.

Electron guns are generally designed using particle-in-cell simulation codes such as EGUN,⁷ DEMEOS,⁸ TRAK,⁹ UGUN,¹⁰ ARGUS,¹¹ AVGUN,¹² BOA,¹³ COCA,¹⁴ and MICHELLE.¹⁵ These codes solve the equations of motion for electron trajectories in self-consistent electric and magnetic fields, starting at a cathode and ending at the exit of the gun where the beam is injected into some other structure. We denote the trajectory of particle j as $(xj(t),pj(t))$⁠, where x and p are the position and momentum as functions of time since leaving the cathode; trajectory j is assumed to carry a given amount of steady current, I_j. The spatial distributions of charge and current density are

where $vjt$ is the velocity of particle j and $Tj$ is the transit time of particle j through the system, that is, the time from emission to exit. These distributions are used as sources in Maxwell's equations to update the static fields

where $ϕx$ and $Ax$ are the electric and magnetic potentials, respectively. The process of launching particles, integrating trajectories, accumulating sources, and updating fields is repeated until a converged result is obtained. An example of a diode with planar symmetry is shown in Figs. 1 and 2. In this case, the beam (shown in green in Fig. 2) is emitted from a cathode on the left, electrostatically focused by a pair of electrodes and extracted on the right through a drift space. The 2D calculation involved following the trajectories of 600 simulation particles through a domain consisting of a grid of 120 000 triangles. Typically, the process of computing trajectories, accumulating charge densities, and calculating self-consistent fields was iterated 200 times to reach a highly accurate converged solution.

Let us suppose that we have a reference design in which particle trajectories have been traced in combined electric and magnetic fields, and the fields have been calculated self-consistently with the particle trajectories, as in Fig. 2. This represents our base system state X₀. We consider a perturbation of this state such that trajectories $δxj,δpj$⁠, charge and current densities $(δρ,δj)$⁠, and potentials $(δϕ,δA)$ all change. The equations of motion can be written in terms of a perturbed $q(δΦ−vj⋅δA)$ and an unperturbed Hamiltonian, H(p, x). The equations for the perturbed trajectories are

We now consider two distinct perturbations: one (⁠$δX$⁠) caused by changing the parameters of the system (⁠$δB$⁠) and the other that will be related to the adjoint solution (⁠$δY$⁠). We refer to these as the “true” and “adjoint” perturbations and label them with superscripts X and Y, respectively.

We dot Eqs. (3a) and (3b) for each of the perturbation types (X or Y) with the change in the position (3a) and momentum (3b) for the other perturbation type, subtract, multiply by the beamlet current for each trajectory, sum over trajectories, and integrate over time. The result is

The left hand side of (4) is the current-weighted sum of the change, from the cathode to a point near the exit, of the symplectic area defined by the pairs of trajectories. We will assume that the true trajectories (X) are unperturbed at the cathode, and consequently, the left-hand side of (4) only depends on perturbed trajectories near the exit. We will choose the perturbations of the adjoint trajectories (Y) near the exit such that the left side becomes the change in the metric $δM$ whose sensitivity function we seek.

Before doing this, we note that the evaluation point is the time of flight of the unperturbed orbit, T_j. We want to replace this with an evaluation at, not just near, the exit of the device. The perturbed orbits will have different times of flight, given by $Tj+δTj(X,Y)$⁠, to the exit of the domain. We can express perturbed quantities at the unperturbed time of flight in terms of their values at the exit of the domain and the rate of change of the corresponding unperturbed coordinate

Using Hamilton's equations for the unperturbed trajectory, the left-hand side of (4) becomes

where the z components of the perturbed displacements now vanish. We note that the same result can be obtained by initially making a canonical transformation in which z replaces t as the independent variable. Here and below, we will consider only perturbations that leave the Hamiltonians unchanged, i.e., $δHj(X,Y)=0$⁠.

The right side of (4) can be expressed in terms of the perturbed charge and current density. For example

A similar relation applies to the vector potential terms. The result is

The right-hand side of (5) will become the sensitivity function. Its expression depends on the boundary conditions that apply to the solutions of the field Eq. (2).

We assume that electrostatic conditions are applied at the simulation boundary where either the potential or its normal derivative is specified for the true solution and the perturbed potential $δΦ(Y)$ vanishes for the adjoint solution. For the magnetostatic fields, we assume that the true system is perturbed by a change in current density in coils or in the magnetization of permanent magnets, $δjm(x)$⁠; we also assume that all magnetic fields approach zero at infinity. We then apply Green's theorem to the integral on the right side of (5) to obtain

Equation (6) is our central mathematical result. It is a form of Green's theorem extended to calculations using a mixture of fields and discrete particles. It states that the weighted sum of symplectic areas at the exit of the gun is given in terms of the perturbed electrostatic potential on the boundary and the perturbed current density creating the applied magnetic field; note that a displacement $Δx$ of the boundary is equivalent to a perturbed potential $δΦX|B=−Δx·∇Φx$⁠. These quantities are weighted by the sensitivity functions $n⋅∇δΦ(Y)$ and $δA(Y)$⁠; once these are known, Eq. (6) may be used to compute the symplectic areas for any perturbation $δϕX, δjmX$⁠. We show below how the symplectic area may be related to a figure of merit of interest by the suitable choice of the adjoint particle conditions $δxjY, δpjY$ at the exit plane. We pick the final coordinates and integrate the trajectories backward in time through the gun, determining the self-consistent fields iteratively, exactly as in the forward integration for the unperturbed case.

We consider three examples of figures of merit for the gun of Fig. 1: (a) the average vertical displacement of the beam, (b) the RMS width of the beam, and (c) the emittance of the beam. In the first case, we pick as the final perturbed coordinates $(δxj(Y)=0, δp⊥j(Y)=λe⊥)L$⁠. That is, the perturbed momentum perturbation at the exit is a constant in the direction of the unit vector, $ê⊥$⁠. In this case, the left-hand side of (6) becomes $−λ∑jIjê⊥⋅δxj(X)|L/I=−λê⊥⋅⟨δxj(X)⟩$⁠, where the angular brackets indicate a current weighted average. Thus, the left side of (6) is the average displacement in the direction of the unit vector. Both sides of Eq. (6) are proportional, in principle, to the constant $λ$⁠, and it can then be divided out or set to unity to obtain the desired sensitivity. This would be the case if in practice we solved linearized Eq. (3). However, this would require extensive programming to implement in the existing codes. Instead, we solve the nonlinear equations and select $λ$ to be sufficiently small, essentially taking a numerical derivative. This will be made clearer in the second example.

To determine the sensitivity of the RMS beam width, we pick as final perturbed coordinates $(δxj(Y)=0, δp⊥j(Y)=λx⊥j)L$⁠, where $λ$ is again a constant, to be chosen so that the adjoint perturbations are sufficiently small. Then, the left side of (6) becomes essentially the change in the RMS radius of the beam for the true solution

We thus can use the self-consistent adjoint fields as the sensitivity function for the RMS beam radius.

To test this approach, we will compute the change in the RMS beam thickness in two ways, for the device in Fig. 1, due to a change in potential on the anode. The first way is to directly simulate the device with a range in applied voltages. This will allow us to determine the range in voltage in which the changes in the RMS thickness are linear in the change in applied voltage. Shown in Fig. 3 with round blue symbols is the change in the mean squared beam width divided by the change in voltage as a function of the relative change in voltage. The line appears to be straight as would be expected for a numerically obtained derivative and shows an intercept of about 5.96 × 10⁻⁸ m²/V as the change in voltage decreases to zero. This is to be compared with the value obtained from the adjoint formulation shown by the curve with the red squares in Fig. 3. Here, we have computed the sensitivity function in the form of $−n⋅∇ϕ(Y)/λ$ obtained from the backward integration of the particles with final conditions specified according to (7). The sensitivity function is defined on the boundaries of the simulation domain and is shown as arrows in Fig. 2. This emphasizes that the sensitivity function can be used to determine the change in the mean squared beam thickness for any change in the anode shape or potential, not just the uniform change in potential considered here. We compute the expected change in the mean square beam width for a one volt change in anode potential by doing the surface integral in (6) numerically. The results are then plotted in Fig. 3 as a function of the constant $λ$⁠. For $λ>0.001$⁠, the curve shows expected linear behavior with an intercept that is within 1% of the value obtained from the true calculations. For extremely small values of $λ<0.001$⁠, the accuracy of the numerical differentiation is lost due to a number of limitations associated with numerical differentiation.

As a third example, we may determine the sensitivity of the beam emittance to changes in electrode placement. To do so, we consider changes in the quantity

where we have assumed that the average of the transverse positions and momenta vanish. The change of the emittance is then given by

Thus, if we choose

and

then the left-hand side of (6) becomes $λδε2/2$⁠, and so, the change in beam emittance due to changes in electrode potential or position may be computed using Eq. (6).

Finally, we generalize to figures of merit of the form

that is, figures of merit based on a current weighted sum of a function of phase space coordinates at the exit of the simulation domain. The deviation in this figure of merit that results from changes in electrode shapes, positions, and potentials or magnet strengths and locations may be written as

Thus, if we select

and

then (6) gives the change in this figure of merit due to changes in the focusing geometry

A sensitivity function for design of electron beam optics has been introduced. The sensitivity function, which gives the change in some metric for arbitrary changes in electrode potentials and shapes and arbitrary changes in coil or magnet locations, can be calculated with a few runs of the codes currently used to design electron guns. While the present situation is static, the approach can be generalized to time varying situations as well.

This work was supported by DARPA Contract No. HR0011-16-C-0080 with Leidos, Inc. The views, opinions, and/or findings expressed are those of the authors and should not be interpreted as representing the official views or policies of the Department of Defense or the U.S. Government. Approved for Public Release, Distribution Unlimited.