Approximating Nonlinear Forces with Phase-Space Decoupling

TL;DR

This paper introduces a phase-space decoupling method using virtual coordinates to efficiently approximate nonlinear forces in beam dynamics, bridging the gap between fast linear and slow multiparticle simulations.

Contribution

The paper presents a novel phase-space decoupling technique with virtual coordinates to efficiently model nonlinear effects in beam tracking.

Findings

Significant reduction in computational cost

Effective approximation of nonlinear forces

Bridges gap between linear and multiparticle simulations

Abstract

Beam tracking software for accelerators typically falls into two categories: fast envelope simulations limited to linear beam optics, and slower multiparticle simulations that can model nonlinear effects. To find a middle ground between these approaches, we introduce virtual coordinates in position and momentum which have a cross-dependency (i.e. p=f(x) where x is an initial position and p* is a virtual projection of momentum onto the position axis). This technique approximates multiparticle simulations with a significant reduction in calculation cost.

Equations

$$\int {\rho }_x^{L}d{x}_L=\int {\rho }_x^{0}d{x}_0=\int {\rho }_p^{L}d{p}_L=\int {\rho }_p^{0}d{p}_0=N.$$

$$p_0\approx -\frac{\alpha }{\beta }{x}_0,$$

$$\int {\rho }_x^{0}d{x}_0$$

$$\frac{1}{2}\mathrm{erf}\left(\frac{\sqrt{2}{x}_0}{2{\sigma }_{x_0}}\right)$$

$$p_0=\sqrt{2}{\sigma }_{p_0}\mathrm{erfinv}\left(\mathrm{erf}\left(\frac{\sqrt{2}{x}_0}{2{\sigma }_{x_0}}\right)\right)=x_0\frac{{\sigma }_{p_0}}{{\sigma }_{x_0}}.$$

$$\int {\rho }_p^{\ast }(p_0)d{p}^{\ast }(x_0)=\int {\rho }_p^{0}d{p}_0,$$

$$N=x^{\ast }{\rho }_x^{\ast }=\int {\rho }_x^{0}d{x}_0$$

$$N=p^{\ast }{\rho }_p^{\ast }=\int {\rho }_p^{0}d{p}_0.$$

$$N^2=\mathcal{C}={\rho }_x^{\ast }{\rho }_p^{\ast }{x}^{\ast }{p}^{\ast },$$

$$\frac{d}{d{x}_0}\mathcal{C}+\frac{d}{d{p}_0}\mathcal{C}=0.$$

$$0$$

$$\displaystyle\frac{\partial{p^{*}}}{\partial{x_{0}}}=-p^{*}\frac{\rho^{0}_{x}}{\int{\rho^{0}_{x}}dx_{0}}=-\frac{p^{*}\rho^{0}_{x}}{\Upsilon_{x}},-\frac{2p^{*}\rho^{0}_{x}}{\Upsilon_{x}}\$$

$${{\rm Y}}_x\equiv \frac{1}{2}\mathrm{erf}\left(\frac{\sqrt{2}{x}_0}{2{\sigma }_x}\right).$$

$$p^{\ast }$$

$$x_L(x_0,p^{\ast })=\left\{\exp \left[-t:\mathcal{H}\left(x_0,p_0\right):\right]x_0\right\}{|}_{p_0\to p^{\ast }}.$$

$$\mathcal{H}=\frac{e}{p}\frac{k\cdot \mathfrak{Re}(x_0+i{y}_0{)}^n}{a_0^{n-1}n!}+\frac{(p_0{)}^2}{2m}.$$

$$x_1,y_1,p_{x1},p_{y1}$$

$$x_2,y_2,p_{x2},p_{y2}$$

$$p^{\ast }\approx -x_0\frac{{\sigma }_p}{{\sigma }_x}.$$

Full text

1 Introduction

The software for predicting and correcting beam dynamics in real time is largely dependent on approximations that treat the phase-space density as a single envelope [1, 2], while multiparticle tracking codes that can account for nonlinear forces are more CPU-intensive and only suitable for machine design or offline diagnostics [3, 4]. In developing the ESS Linac Simulator (ELS), our group intends to incorporate nonlinear tracking without sacrificing real-time diagnostics capability.

Expanding on prior works [5, 6], we introduce a multiparticle optimization method by using “virtual” phase-space coordinates that allow for the independent calculation of position and momentum densities. It is important to classify these coordinates as non-physical: they are derived assuming a cross-dependency exists (e.g. $p^{*}=f(x_{0})$, where $p^{*}$ and $x_{0}$ are virtual momentum and real position, respectively).

When used in conjunction with standard techniques for multiparticle tracking, these virtual coordinates allow for the use of monovariate polynomials, which present a significant reduction in the number of required calculations when compared with the bivariate polynomials normally needed.

Although the results to follow consider only multiparticle tests, these techniques may be integrated into beam-envelope simulations by building contour maps from low-particle-count samples.

2 Theory

For simplicity, we will begin by considering only transverse motion along a single axis with an initial Gaussian distribution (though we will test the resulting approximation in 2D on various initial distributions).

The starting point for this method, outlined in Ref. [6], involves taking particle count $N$ as an invariant as the position and momentum envelopes $\rho_{x}$ and $\rho_{p}$ evolve:

[TABLE]

where $\rho_{x}^{L}$ and $\rho_{x}^{0}$ are the respective final and initial position densities, and likewise for the momentum densities.

To exploit this identity, $\rho_{x}^{L}$ must be independent of $p_{0}$. This can be accomplished using the approximation [7, 8]

[TABLE]

where $\alpha$ and $\beta$ are the well-known Twiss parameters. Unfortunately, this approximation is only valid if the previous history of the beam is linear (thus maintaining elliptical phase space densities). Since we want an algorithm that remains accurate for iterated nonlinear kicks (which develop irregular density profiles), a new approximation is needed.

Proceeding under the constraint that the initial distribution is Gaussian in both $x$ and $p$, we can solve Eq. 1 for $x_{0}$ and $p_{0}$.

[TABLE]

which yields

[TABLE]

This solution does produce a bigaussian phase-space ellipse, but is unsuitable for approximating $p_{0}$ with distributions of an irregular shape.

We continue by guessing that a solution exists for $p_{0}~{}=~{}f(x_{0})$ for irregularly shaped distributions. Denoting these solutions as $p^{*}$ and $\rho_{p}^{*}$ for momentum and momentum density, respectively, we have

[TABLE]

where it is critical to note that $p^{*}$ is solely dependent on $x_{0}$, while $\rho_{p}^{*}$ remains a function of $p_{0}$ (and likewise for $x^{*}$ and $\rho^{*}_{x}$). With the left-hand-side integrand and integration variables decoupled, it follows (for both $x^{*}$ and $p^{*}$):

[TABLE]

Exploiting particle-count invariance again, and squaring $N$, we can assert that

[TABLE]

and thus

[TABLE]

We can then simplify, treating all $p_{0}$-dependent terms as $f(p_{0})$:

[TABLE]

Then, dividing by $\int{\rho^{0}_{p}}dp_{0}$, the expanded $f(p_{0})$ terms become zero and we have:

[TABLE]

By reusing Eqn. 6, all $p_{0}$ dependence can be eliminated, leaving

[TABLE]

where the second solution can be obtained integrating by parts, and, in the case of a Gaussian initial distribution, the placeholder in the denominator is defined as

[TABLE]

Thus, in contrast with Eq. 4, we have an expression where $\frac{\partial{p}}{\partial{x}}$ is no longer constant.

We now check the following approximation:

[TABLE]

Where we normalize $\mathcal{D}$ using Eq. 4; setting to $p^{*}\approx-p_{0}$ near $|x_{0}|=0$ , leaving

[TABLE]

which can be shown numerically to agree with Eqn. 11 for $|x_{0}|\lesssim 6\ \sigma_{x}$.

At this point, the updated particle postition $x_{L}$ can be calculated using an exponential Lie-operator method [9]:

[TABLE]

where $t$ is elapsed time in the lab frame and the Hamiltonian for a normal multipole magnet in the transverse plane is

[TABLE]

Here, $e$, $p$, $m$, and $a_{0}$ are the fundamental charge, reference longitudinal momentum, particle mass, and magnet-pole raidus, respectively; $n=3,4,5...$ for sextupoles, octupoles, decapoles, etc; and $k$ has units of $[\rm{T}\cdot\rm{m}^{-1}]$. In the following sections, longitudinal momentum is normalized to 1 GeV/c and $a_{0}$ is set to 20 mm unless otherwise noted.

In implementing Eq. 15, $\mathcal{H}$ must be calculated symbolically first for each element. Then, $p^{*}(x_{0})$ and $x_{0}$ are substituted in at each step, reducing the bivariate $x_{L}(x_{0},p_{0})$ to a monovariate $x_{L}(x_{0},\sigma_{x},\sigma_{p})$, where $\sigma_{x}$ and $\sigma_{p}$ remain constant for a given timestep.

Although an analogous $x^{*}(p_{0})$ can be derived, it is not useful in practice. Specifically, in calculating Eq. 16 in 2D for position and momentum – $x_{L}(x_{0},y_{0},p^{*}_{x},p^{*}_{y})$ and $p_{xL}(x^{*},y^{*},p_{x0},p_{y0})$ – the resulting $x_{L}$ expression is dependent on $\sigma_{x}$ and $\sigma_{px}$, while $p_{L}$ is dependent on $\sigma_{x}$, $\sigma_{y}$, $\sigma_{px}$, and $\sigma_{py}$, rendering it computationally inefficient. Other schema involving alternate forms such as $p_{L}(x_{0},p^{*})$ have been checked, but the following is found to be most stable, with notable performance gains:

[TABLE]

Where the $D$ subscript denotes a drift space of at least five times the kick length. This effectively limits the technique to a thin-lens approximation. Such drift spaces can be reserved for incorporating space-charge effects, leading to a comparable number of calculation steps using Eq. 2 versus a standard nonlinear beam-physics code.

3 Multiparticle Simulation

Figure 1 compares the accuracy of multiparticle transformations following Eq. 2 with and without using $p^{*}$. Also shown are tests for $p^{*}=0$ and a “naive” approximation, where $p^{*}\approx p_{0}$ from Eq. 4 is used in the low $|x_{0}|$ limit of Eq. 14:

[TABLE]

To emphasize visible discrepancies, the results shown have their $\sigma$ values updated after each timestep by taking a new standard deviation. However, if mean absolute deviations are taken instead, an improved matching with the baseline can be observed.

Both the naive and null-momentum approximations fail at $\sigma_{p}\gtrsim\sigma_{x}$ (i.e. at energies exceeding 1 GeV). Figure 2 illustrates such a case for 2D Gaussian proton distributions with a kinetic energy of 8 GeV passing through an octupole magnet. In the 2D case, beam parameters were derived relativistically from Twiss parameters and $B_{0}$ by normalizing the kinetic term in Eq. 16 to the beam’s average kinetic energy then verified against Tracewin [10].

4 Conclusion

For the non-null $p^{*}$ approximations, performance improves with increasing particle count, with increasing magnetic pole count, and particularly with increased order of Lie-transform series truncation (Fig. 3). Since trajectory variations are negligible beyond a 6th-order truncation in most cases, the average reduction in CPU overhead using Eq. 14 is roughly 15%.

Similar results were obtained for sextupoles, decapoles, and high-order magnets, as well as with waterbag distributions, despite the assumption of a Gaussian shape in deriving Eq. 14. At low energies (or specifically, any low $\frac{\sigma_{p}}{\sigma_{x}}$ ratio), all three approximations tested have essentially identical results, with escalating performance in the following order: $p^{*}~{}=~{}(f[\operatorname{erf}])$;  $p^{*}~{}=-x_{0}\cdot\sigma_{p}/\sigma_{x}$;  $p^{*}~{}=~{}0$  .

For all the approximations tested, trajectories only became unstable in cases where the momentum of the baseline exceeded ${\sim}100\sigma_{p}$. Thus, the major limitation to this technique is its large drift–kick ratio requirement.