Hamiltonian preserving nonlinear optics

TL;DR

This paper introduces a method for designing nonlinear accelerator lattices that approximately preserve a Hamiltonian integral of motion, enhancing the stability and performance of particle accelerators.

Contribution

It presents a novel approach linking Hamiltonian integrators with real lens arrangements, improving nonlinear lattice design for accelerators.

Findings

Improved nonlinear lattice designs for IOTA and UMER facilities.

Application of high-order symplectic integrators to accelerator lattices.

Proposal of new lattice design using Yoshida integrator.

Abstract

In this paper we present a method of constructing a nonlinear accelerator lattice that has an approximate integral of motion that is given upfront. The integral under consideration is a Hamiltonian in normalized (canonical) coordinates that is preserved by a lattice with a given accuracy. We establish a connection between the integrator of a Hamiltonian in normalized coordinates and a real lens arrangement. We apply known algorithms of high-order symplectic integrators, to produce several nonlinear lattices and show that this approach could improve the design of the nonlinear insert considered at the IOTA and UMER facilities. We also suggest new lattice design based on the Yoshida integrator.

Tables

Table 1. Table 1: Initial conditions for the Poincare surface of section Fig. 6 for the sextupole lattice and Henon-Heiles Hamiltonian Eq.( 41 ).

set number	$q_1$	$p_1$	$q_2$	$p_2$
1	0.095	0.096	0	0.427003
2	0.15	0.0960	0	0.412958
3	-0.1	0.157	0	0.405814
4	-0.2	0.05	0	0.390086
5	0.2	0.05	0	0.403527
6	-0.12	0.005	0	0.429445

Table 2. Table 2: Initial conditions for the Poincare surface of section Fig. 9 for the octupole lattice and Henon-Heiles Hamiltonian Eq.( 45 ).

set number	$q_1$	$p_1$	$q_2$	$p_2$
1	0.02	0.06	0	0.629285
2	0.04	0.12	0	0.619676
3	0.06	0.18	0	0.603319
4	0.078	0.22	0	0.58779
5	0.1	0.3	0	0.547677
6	0.12	0.36	0	0.505862
7	0.14	0.42	0	0.451451
8	0.16	0.48	0	0.379041
9	0.18	0.54	0	0.274727
10	0.4	0.05	0	0.474025

Table 3. Table 3: Initial conditions for the Poincare surface of section Fig. 12 for the Yoshida lattice composed of three nonlinear magnets with the potential given by Eq.( 49 ).

set number	$q_1$	$p_1$	$q_2$	$p_2$
1	0	0.07967	0.1	0.54
2	0	0.41716	0.6	0.01
3	0	0.54447	0.13	0.005
4	0	0.54447	-0.13	-0.005
5	0	0.33064	0.7	0.01
6	0	0.47124	0.5	0.01
7	0	0.54597	0.04	0.04

Equations

$$\dot{X}(t)$$

$$X(0)$$

$$\dot{X}$$

$${\Psi }_{t=mh}=({\varphi }_h^{[1]}\circ {\varphi }_h^{[2]}{)}^m\text{leavevmode}m\in \mathbb{N}.$$

$$\dot{X}=-:\mathrm{H}:X.$$

$$:\mathrm{H}:g\equiv \{\mathrm{H},g\}=\sum _{i=1}^n\frac{\partial \mathrm{H}}{\partial q_i}\frac{\partial g}{\partial p_i}-\frac{\partial \mathrm{H}}{\partial p_i}\frac{\partial g}{\partial q_i}$$

$$X(t)=\exp \left(-t:\mathrm{H}:\right)X_0.$$

$$\exp \left(-t:{\mathrm{H}}_1:\right)\circ \exp \left(-t:{\mathrm{H}}_2:\right)=$$

$$\exp \left(-t:{\mathrm{H}}_1+{\mathrm{H}}_2:+\frac{t^2}{2}:\{{\mathrm{H}}_1,{\mathrm{H}}_2\}:-\mathcal{O}(t^3)\right).$$

$$\exp \left(-h:{\mathrm{H}}_1:\right)\circ \exp \left(-h:{\mathrm{H}}_2:\right)=$$

$$\exp \left(-h:{\mathrm{H}}_1+{\mathrm{H}}_2:+\mathcal{O}(h^2)\right).$$

$$\mathrm{H}=\sum _{i=1}^n\frac{q_i^2+p_i^2}{2}+V(q_1,q_2,...,q_n).$$

$$V(q_1,q_2,...,q_n)=\sum _{j=3}^l{a}_j{P}_j(q_1,q_2,...q_n),$$

$${\mathrm{H}}_1=\sum _{i=1}^n\frac{q_i^2+p_i^2}{2}$$

$${\mathrm{H}}_2=V(q_1,q_2,...,q_n)$$

$$f_1(q_1,p_1,q_2,p_2,...,q_n,p_n)=$$

$$[p_1,-q_1,p_2,-q_2,...,p_n,-q_n{]}^{\mathrm{T}}.$$

$$\left[\begin{array}{c}{\dot{q}}_i\\ {\dot{p}}_i\end{array}\right]=\left[\begin{array}{ccc}0& & 1\\ -1& & 0\end{array}\right]\left[\begin{array}{c}{q}_i\\ p_i\end{array}\right],\text{leavevmode}i\in \mathbb{Z}(1,n).$$

$$R_{\psi }^{(i)}=\left[\begin{array}{ccc}\cos (\psi )& & \sin (\psi )\\ -\sin (\psi )& & \cos (\psi )\end{array}\right].$$

$$\left[\begin{array}{c}{q}_i\\ p_i\end{array}\right]=R_{\psi }^{(i)}\left[\begin{array}{c}{q}_i^0\\ p_i^0\end{array}\right].$$

$$f_2(q_1,p_1,q_2,p_2,...,q_n,p_n)=$$

$${\left[0,-{\partial }_{q_1}V,0,-{\partial }_{q_2}V,...,0,-{\partial }_{q_n}V\right]}^{\mathrm{T}}.$$

$$\left[\begin{array}{c}{\dot{q}}_i\\ {\dot{p}}_i\end{array}\right]=\left[\begin{array}{c}0\\ -{\partial }_{q_i}V(q_1,q_2,...,q_n)\end{array}\right],\text{leavevmode}i\in \mathbb{Z}(1,n).$$

$$K_{\psi }{X}_0=[q_1^0,p_1^0-\psi {\partial }_{q_1}V,\text{leavevmode}...\text{leavevmode},q_n^0,p_n^0-\psi {\partial }_{q_n}V{]}^{\mathrm{T}},$$

$${\Psi }_h=K_h\circ R_h.$$

$$R_{h/2}\circ K_h\circ R_{h/2}=$$

$$\exp \left(-\frac{h}{2}:{\mathrm{H}}_1:\right)\circ \exp \left(-h:{\mathrm{H}}_2:\right)\circ \exp \left(-\frac{h}{2}:{\mathrm{H}}_1:\right)=$$

$$=\exp \left(-h:{\mathrm{H}}_1+{\mathrm{H}}_2:+\mathcal{O}(h^3)\right).$$

$${\Phi }_h=R_{h/2}\circ K_h\circ R_{h/2}$$

$${\varphi }_{{\gamma }_{s}h}\circ \text{leavevmode}...\text{leavevmode}\circ {\varphi }_{{\gamma }_{1}h}$$

$${\gamma }_1={\gamma }_3=\frac{1}{2-2^{1/(p+1)}},\text{leavevmode}\text{leavevmode}\text{leavevmode}{\gamma }_2=-\frac{2^{1/(p+1)}}{2-2^{1/(p+1)}}.$$

$${\Phi }_h^Y={\Phi }_{{\gamma }_{3}h}\circ {\Phi }_{{\gamma }_{2}h}\circ {\Phi }_{{\gamma }_{1}h},$$

$${\gamma }_1={\gamma }_3=\frac{1}{2-2^{1/3}},\text{leavevmode}\text{leavevmode}{\gamma }_2=-\frac{2^{1/3}}{2-2^{1/3}}.$$

$${\Phi }_h^Y=$$

Full text

I Introduction

The concept of integrable and quasi-integrable nonlinear optics has recently attracted significant attention. Initially suggested by Danilov (see Ref.Dan and references therein) and refined by Danilov and Nagaitsev in Ref.DanNag , the concept has been expanded to more realistic cases with space charge and chromaticity effects accounted for Webb ; Wal . Experimental demonstration of the integrable optics concept is currently being conducted at the IOTA facility at Fermilab IOTA ; Ant_o as well as at UMER ring at the University of Maryland Umer .

The main idea behind the integrable optics concept is a special insert of nonlinear magnets that is accommodated by a purely linear ring. The system is arranged in a way that the effective Hamiltonian for the lattice is almost time-independent and the potential produced by one nonlinear magnet warrants separation of variables, and thus a second integral of motion DanNag ; chad_r . Initial designs of the nonlinear insert considered for both experiments Ant_o ; Umer were based on an idea of approximating the smooth nonlinear potential with a certain number of nonlinear magnets (17 in the case of IOTA and 7 in the case of UMER) with their strength scaled according to a prescription derived in Ref.DanNag and placed equidistantly. Questions remain whether this number could be reduced further, and if the performance and design of the nonlinear insert could be further enhanced.

In this paper we introduce a general method of designing a nonlinear lattice based on known symplectic integration methods. After establishing a connection between the integrator of a smooth Hamiltonian in the normalized (canonical) coordinates, and a fragment of a real lattice, we demonstrate that the nonlinear lattice that preserves a given nonlinear smooth Hamiltonian could be implemented with only three nonlinear elements. As a representative example, we consider sextupole nonlinear inserts for the IOTA ring, octupole nonlinear inserts for the IOTA ring, and a toy-model FODO lattice with three nonlinear magnets introduced in Ref.DanNag .

The suggested method may be used in combination with other established tools like the normal form analysis Baz1 ; Baz2 ; TurchettiNOCE ; Dragt ; Forest1 , controls of symplectic maps HSM , and methods for the increase of integrability CaryPRAB to design lattices with high dynamic aperture.

II General approach

In this section we introduce notation, terminology and give a brief overview on symplectic integration of an autonomous Hamiltonian system. We slightly alter commonly used splitting of the flow to establish a connection between the integrator for a given Hamiltonian and a magnet arrangement of a real lattice.

For the reader’s convenience we list definitions of several terms that are used through the paper in the Appendix C.

II.1 Theoretical background

In this section we will give a brief overview of the mathematical tools and terminology that we use in the paper. For more details on the numerical integration and advanced concepts we refer the reader to the original book GInt . Basic concepts about Lie algebras and geometric methods in differential equation theory could be found in Refs. Arnold ; Arnold2 , applications of this ideas to the accelerator physics problems could be found in for example Refs. Dragt ; Forest1 .

Let us consider a nonlinear autonomous system in $\mathbb{R}^{2n}$

[TABLE]

Here, $X=[q_{1},p_{1},q_{2},p_{2},...,q_{n},p_{n}]^{\mathrm{T}}$ is a $2n$-dimensional vector of positions and momenta at a time $t$ and $X_{0}=[q_{1}^{0},p_{1}^{0},q_{2}^{0},p_{2}^{0},...,q_{n}^{0},p_{n}^{0}]^{\mathrm{T}}$ is the vector of initial conditions, $f(X)$ is a vector function called a vector field. Here and throughout the paper, the dot above the letter denotes the full time derivative. The flow, $\phi_{t}$, of the system (1) is the mapping that establishes a connection between initial condition $X_{0}$ and some point in time $\phi_{t}(X_{0})=X(t)$.

If the vector field $f(X)$ could be represented as $f(X)=f_{1}(X)+f_{2}(X)$ such that systems

[TABLE]

could be exactly integrated and the corresponding flows $\phi^{[1]}_{t}$ and $\phi^{[2]}_{t}$ could be explicitly found, then one can build an approximate flow $\Psi$ of the initial system as follows.

Let us consider a time step $h<1$. Then up to the order $\mathcal{O}(h^{2})$ the approximate flow on a time mesh with the step $h$ is

[TABLE]

Mapping $\Psi_{t=h}$ is simply an integrator of the first order by $h$ (Euler integrator) of the system (1).

The possibility of such splitting, and even its symplectic (volume preserving) property for any autonomous Hamiltonian system, could be easily seen from the following considerations.

Let $\mathrm{H}$ be a Hamiltonian and $X$ still a $2n$ vector of coordinates and momenta. Then from the Hamiltonian equations, the trajectory of the system could be found from

[TABLE]

Here ${:}\mathrm{H}{:}$ is the Lie operator with an action on $g$ defined by the Poisson bracket of $g$ with $\mathrm{H}$ as Dragt

[TABLE]

The introduced Lie operator has a simple connection to a known Liouville operator ${:}\mathrm{H}{:}=-i\widehat{\mathbf{L}}$.

The solution to the equation (5) could be represented as an exponent of the Lie operator acting on a vector of initial conditions (see for example Arnold )

[TABLE]

The operator exponent in the above expression, in general is an infinite series, however, for some special Hamiltonians this series naturally truncates at a finite number of terms leading to the exact expression for the map $\exp\left(-t{:}\mathrm{H}{:}\right)$.

Let us assume that Hamiltonian $\mathrm{H}$ could be split into $\mathrm{H}=\mathrm{H}_{1}+\mathrm{H}_{2}$ in a way that $\exp\left(-t{:}\mathrm{H}_{1}{:}\right)$ and $\exp\left(-t{:}\mathrm{H}_{2}{:}\right)$ could be evaluated explicitly. As far as the Lie operators ${:}\mathrm{H}_{1}{:}$ and ${:}\mathrm{H}_{2}{:}$ do not commute, we utilize the Baker-Campbell-Hausdorff (BCH) formula Haus to evaluate the composition

[TABLE]

If we introduce a time mesh with the step $h<1$ we will see that up to the order $\mathcal{O}(h^{2})$, the composition of the exponents is the exponent of the sum of Lie operators,

[TABLE]

This means that the composition will preserve Hamiltonian $\mathrm{H}=\mathrm{H}_{1}+\mathrm{H}_{2}$ up to the order $\mathcal{O}(h)$.

By comparing the approximate flow $\Psi_{t=h}$ Eq.(4) with Eq.(II.1) and establishing the connection $\exp\left(-h{:}\mathrm{H}_{1}{:}\right)\equiv\phi^{[1]}_{h}$, $\exp\left(-h{:}\mathrm{H}_{2}{:}\right)\equiv\phi^{[2]}_{h}$, we conclude that splitting of the flows is closely connected with the possibility of splitting the Hamiltonian. The exponent of the Lie operator is a symplectic (volume preserving map) GInt ; Arnold ; Arnold2 ; Dragt , the one-step integrator $\Psi_{t=h}$ is also symplectic as a composition of symplectic transformations is a symplectic transformation (group property).

The considerations above lead us to the following observation Tabor : a system with continuous time could be approximated with a discrete system that preserves the smooth Hamiltonian up to a given level of accuracy. It is worth mentioning that in the case of a four dimensional phase space (transverse motion in accelerator is decoupled from longitudinal) when variables in the Hamiltonian could be separated, this leads to quasi-integrability DanNag .

In the following sections, we study the case when the preserved Hamiltonian is close to the desired time independent one, and thus the dynamics are predefined by the desired Hamiltonian. We find magnet arrangements that have effective Hamiltonians of a predefined form, up to an error, that could be expressed as a power of the phase advance between nonlinear elements.

II.2 Discretization of the smooth Hamiltonian and high order integrators

The connection between a smooth Hamiltonian system and a map naturally arises from building an integrator that is essentially a discrete analog of a smooth system by definition Tabor ; GInt . The more accurate the integrator, the better it reproduces dynamics of the original system. In this section we give a brief derivation of the specific forms of the known integrators that we are going to utilize further.

We consider a smooth Hamiltonian of the form

[TABLE]

Here $V$ is the nonlinear potential of the form

[TABLE]

$P_{j}$ is a homogeneous polynomial of the degree $j$, and $a_{j}$ is a constant.

We split the Hamiltonian into a part that corresponds to linear motion

[TABLE]

and a part

[TABLE]

that combines all nonlinearities. According to Hamiltonian equations, the vector field (or simply the force) that corresponds to the Hamiltonian $\mathrm{H}_{1}$ is

[TABLE]

This leads to a system of $n$ independent pairs of equations

[TABLE]

The flow of this system is an $n$ block-diagonal matrix $R_{\psi}$ of rotations with one block given as

[TABLE]

Solutions are independent pairs of $[q_{i},p_{i}]^{\mathrm{T}}$ that are given in terms of the flow as

[TABLE]

Here $[q_{1}^{0},p_{1}^{0},q_{2}^{0},p_{2}^{0},...,q_{n}^{0},p_{n}^{0}]^{\mathrm{T}}$ is a vector of initial conditions.

The vector field that corresponds to the Hamiltonian $\mathrm{H}_{2}$ is given by

[TABLE]

Here $\partial_{q_{i}}$ stands for the partial derivative by $q_{i}$. The corresponding system of differential equations again could be written as $n$ pairs of $[q_{i},p_{i}]$, however, they are no longer completely independent.

[TABLE]

As far as the vector field keep $q_{i}$ unchanged and only modifies $p_{i}$, the flow of the system above is defined as

[TABLE]

with each partial derivative taken at the initial point $(q^{0}_{1},q^{0}_{2},...q^{0}_{n})$.

Now if we consider a time mesh, $t=mh\leavevmode\nobreak\ m\in\mathbb{N}$, with a step $h$, then the one step integrator $\Psi_{h}$ of the Hamiltonian $\mathrm{H}$ will have the from

[TABLE]

This is a well known symplectic Euler method (Fig.1 left diagram).

Next we consider a composition $R_{h/2}\circ K_{h}\circ R_{h/2}$. Keeping in mind that the exponent of the Lie operator is simply the flow: $\exp\left(-h{:}\mathrm{H}_{1}{:}\right)\equiv R_{h}$, $\exp\left(-h{:}\mathrm{H}_{2}{:}\right)\equiv K_{h}$ and using the BCH formula, we note that GInt

[TABLE]

Due to the symmetry of the composition terms of the order $\mathcal{O}(h^{2})$ cancel out $\frac{h^{2}{:}\{\mathrm{H}_{1},\mathrm{H}_{2}\}{:}}{4}+\frac{h^{2}{:}\{\mathrm{H}_{2},\mathrm{H}_{1}\}{:}}{4}=0$.

From the calculations above it is apparent that the integrator

[TABLE]

preserves the Hamiltonian $\mathrm{H}$ up to the order $\mathcal{O}(h^{2})$ and thus has higher accuracy than the simple Euler method. The integrator Eq.(23) is known GInt as Strang Strang splitting or Marchuk splitting March as well as the second order integrator introduced by Ruth Ruth1 (Fig.1 right diagram).

Having an integrator $\phi$ of given order $p$, it is often useful to build an integrator of higher order by composing $\phi$ with itself. The following theorem gives a general method of building such a composition GInt

Theorem 1

Let $\phi_{h}$ be a one step integrator of the order $p$. If $\gamma_{1}+\leavevmode\nobreak\ ...\leavevmode\nobreak\ +\gamma_{s}=1$ and $\gamma_{1}^{p+1}+\leavevmode\nobreak\ ...\leavevmode\nobreak\ +\gamma_{s}^{p+1}=0$ then the composition

[TABLE]

is an integrator of order at least $p+1$.

An important consequence of the above theorem for $s=3$, $\phi_{\gamma_{3}h}\circ\phi_{\gamma_{2}h}\circ\phi_{\gamma_{1}h}$, is the three step Yoshida integrator Forest ; Suzuki ; Yoshida with corresponding gammas given by

[TABLE]

By composing the integrator (23) ($p=2$) we arrive at the integrator of order $4$ in the form

[TABLE]

By straightforward implementation of the BCH formula one can check that the integrator above preserves the Hamiltonian $\mathrm{H}=\mathrm{H}_{1}+\mathrm{H}_{2}$ up to terms of order $\mathcal{O}(h^{4})$. To avoid negative time steps we use the identity $R_{2\pi}=\mathcal{I}$ and write the integrator (25) in the final form as

[TABLE]

We reiterate one distinctive difference between the commonly used splitting of the Hamiltonian for numerical integration and the one we used above. Commonly, the Hamiltonian is split into a part that purely depends on momentum and a part that is purely dependent on the spatial coordinates. This results in a well known integration method that is sometimes referred as the “drift-kick” method in the most simple Euler implementation. In our case, we separated the Hamiltonian into $\mathrm{H}_{1}$ - corresponding to linear motion and $\mathrm{H}_{2}$ - corresponding to a purely nonlinear “kick” (the same way as in Ref.Bella ). This splitting allows us to establish a direct connection between the integrator of the Hamiltonian in the normalized coordinates and a transformation that corresponds to a set of optical elements.

II.3 Splitting of the nonlinear potential

We consider a nonlinear potential $V$ given by Eq.(11) and recall that it is essentially a sum of several potentials $V=\sum\limits_{j=3}^{l}a_{j}P_{j}$. Each one corresponds to a a specific order of the nonlinearity. As illustrated on Fig.2 we may rewrite the Euler method (21) in the form

[TABLE]

Here $K_{h}^{(j)}$ is a flow that corresponds to the potential $a_{j}P_{i}$. Using the BCH formula one can ensure that this integrator indeed is of the order $\mathcal{O}(h)$ (preserves Hamiltonian $\mathrm{H_{1}}+\sum\limits_{j=3}^{l}a_{j}P_{j}$ up to the order $\mathcal{O}(h)$).

With a slight modification of the theorem 1 one my show that the method

[TABLE]

with $\Psi_{h}^{*}=R_{h/(l-2)}\circ K_{h}^{(3)}\circ...\circ R_{h/(l-2)}\circ K_{h}^{(l)}$ (adjoint method) is a method of the order 2 and thus preserves Hamiltonian $\mathrm{H_{1}}+\sum\limits_{j=3}^{l}a_{j}P_{j}$ up to the order $\mathcal{O}(h^{2})$ (see Ref.GInt ). Further utilization of the method (28) by substituting it into the formula (II.2) will result in increasing the order to the fourth order.

Such splitting is useful for combining different types of nonlinear lenses in the same lattice.

II.4 Connection to a real optical lattice

We assume that the longitudinal and transverse motion are decoupled. Thus we consider a four dimensional phase space. We assume as well that there is no coupling in a linear lattice (transfer matrix has a block diagonal form).

To establish a connection between integrators in normalized coordinates $\{q_{1},p_{1},q_{2},p_{2}\}$ and a real optical lattice, we recall that propagation of the particle from position $s_{0}$ to position $s_{1}$ through a linear optical channel could be described using a block diagonal transfer matrix SYL with the block of the type

[TABLE]

Here, the lower index denotes coordinate pair (either $\{x,P_{x}\}$ or $\{y,P_{y}\}$); $\psi_{x,y}=\int\limits_{s_{0}}^{s_{1}}\frac{ds}{\beta_{x,y}(s)}$ is the phase advance between position $s_{0}$ and $s_{1}$; $\mathbf{B}_{x,y}(s)$ is the corresponding block of the betatron amplitude matrix and $\mathbf{B}_{x,y}^{-1}(s)$ its inverse given by SYL

[TABLE]

Here $\beta_{x,y}(s)$, and $\alpha_{x,y}(s)=-1/2\beta_{x,y}^{\prime}(s)$ are the Twiss parameters of the linear lattice.

Now let us consider an integrator - $\Psi_{h}$, given by (23) and propagate a vector of initial conditions for one step - $h$, that corresponds to the phase advance between the points $s_{0}$ and $s_{2}$, of a linear lattice

[TABLE]

With the identity $\mathcal{I}=\mathbf{B}(s)\circ\mathbf{B}^{-1}(s)$ (here $\mathcal{I}$ is the identity matrix) equation (31) transforms as

[TABLE]

Here, $\mathrm{X}$ denotes the unnormalized state vector $\mathrm{X}_{h,0}\equiv\mathbf{B}X_{h,0}$.

We evaluate $\mathbf{B}(s_{1})\circ K_{h}\circ\mathbf{B}^{-1}(s_{1})$ further to achieve nonlinear element strength scaling with the $\beta$-function in a form

[TABLE]

with

[TABLE]

The two building blocks of the integrator are the flows $R_{h}$ and $K_{h}$, that now with the help of the Eq.(II.4), Eq.(33) and Eq.(34) could be transformed to $M(s_{2}|s_{1})$ and $N^{\beta}$ respectively. Maps $M(s_{2}|s_{1})$ and $N^{\beta}$ could be implemented with thin lenses in a real lattice. We stress that the transformation for a nonlinear kick $N^{\beta}$ is now dependent on a $\beta$ function and thus has to be properly scaled to match the linear part of the lattice. A similar result was achieved in DanNag , however, the distribution of the magnets was considered to be continuous. The results of Ref.DanNag could be considered as a limit for the present approach, in the case of infinite number of nonlinear magnets (setting the method step $h\to 0$ to zero). For the reader’s convenience the passage to the limit $h\to 0$ and recovery of the scaling derived in Ref.DanNag is outlined in the Appendix A.

II.5 Remarks

Having in hands all the necessary tools: the basic integration method that is given by Eq.(21), Theorem 1 for increasing the order, splitting strategy given by Eq.(27) and scaling of the nonlinear lens strength with respect to the linear optics from Eq.(33) and Eq.(34), one may design a nonlinear lattice that will upfront conserve any nonlinear Hamiltonian with a given accuracy. The splitting strategy in Eq.(27) allows to separate different nonlinear lenses in space, thus enabling a variety of lattices that may include sextupole, octupole and even higher order magnets at the same time if needed.

An important observation that follows from Sec.II.2 is that linear lattice configurations are not limited to the case of equal $\beta$-functions $\beta_{x}=\beta_{y}$, the $\beta$-function could differ as far as the condition $\psi_{x}=\psi_{y}$ is satisfied. This means that in order to be consistent with the integrator structure, only the phase advance in $x$ and $y$ have to be equal.

If $J=q_{1}^{2}+p_{1}^{2}+q_{2}^{2}+p_{2}^{2}$ is the action and $a_{j}$ are the amplitudes of the nonlinearities as in Eq. (11) then while action is less than unity $J<1$ and all $a_{j}<1$, higher orders in the BCH formula are suppressed in comparison to the lower orders, and the series could be thought of as a convergent series. Thus for the case $J<1$ and $a_{j}<1$, the whole scheme is stable as a first correction to the Hamiltonian that comes from the discretization error, is proportional to some positive power of action that is greater then $1$ and, on top of that, is multiplied by a small parameter $h^{p}$, where $p$ is the integrator order.

Stability over many iterations (Hamiltonian preservation property) of the integrator, or equivalently, stability of the corresponding lattice over many revolutions, is guaranteed by the following theorem Ben ; GInt

Theorem 2

Let $\mathrm{H}$ (the Hamiltonian) be an analytic function $\mathrm{H}:D\to\mathbb{R}$ (where $D\subset R^{2d}$) and $\Phi_{h}(X)$ a symplectic numerical method of the order $p$ with the step size $h$. If the numerical solution generated by the integrator stays on the compact set $K\subset D$ then there exists a $h_{0}$ such that $\mathrm{H}(X_{n})=\mathrm{H}(X_{0})+\mathcal{O}(h^{p})$ over exponentially long time intervals, $nh\leq e^{h_{0}/2h}$.

Here $d$ - is the number of degrees of freedom and $n$ - is the number of iterations.

Theorem 2 essentially states that a bounded trajectory will remain bounded for exponentially long times and addresses potential concerns regarding long term stability of the corresponding lattice. Practically this means that bounded trajectories of the discrete system will be close to the trajectories of the corresponding smooth system for any realistic time (number of revolutions).

Current consideration of an integrator based lattice has a tight connection to the problem of the resonance elimination and regularization of the particle motion, that is an established tool in creating nonlinear lattices with a large dynamic aperture (DA) CaryPRAB ; TurchettiNOCE . Indeed the route to large DA is to suppress higher order terms in the Lie exponent that could be constructed for example using the Birkhoff normal form approach TurchettiNOCE . The key difference between prior art and present studies though, is in the direction of the analysis. We suggest to perform reverse engineering of a given smooth system and to find a lattice that has the best correspondence to this smooth system in terms of the dynamics. In contrast, in the common approach, derivation and analysis of the effective Hamiltonian for a given lattice is usually the focus. Interestingly, constructing a higher order integrator and a corresponding lattice we, according to the BCH formula, eliminate low order resonances and achieve the same goal. Thus we believe that the present idea may become a useful tool in developing nonlinear lattice design strategies.

III Examples of real lens configurations

In this section we consider two examples of the nonlinear lens arrangement that correspond to the Ruth second order integrator (23) and Yoshida integrator (II.2). To provide head to head comparison in all cases of lattice tracking the integral strength of all nonlinear elements in the lattice was normalized to the unity

[TABLE]

here $f(\beta(s))$ is the amplitude scaling function that is calculated with the help of the Eq.(33) and Eq.(34), $s_{i}$ is the position of the $i$-th nonlinear element and $a_{0}$ is a dimensional multiplier. The final total strength of the nonlinear channel were chosen such that the multiplier in-font of the nonlinear potential for the effective Hamiltonian is equal to unity.

III.1 Ruth lattice for the nonlinear insert

We consider a lattice that is implemented in IOTA Ant_o and UMER Umer for the quasi-integrable and integrable optics experiments. The linear part of the lattice consists of the so called T-insert introduced in DanNag , and a drift of length $L$. The T-insert is an arrangement where part of the linear optics effectively acts as a focusing matrix in both $x$ and $y$ directions, leading to a degenerate case of equal transverse $\beta$-functions in a drift space. This configuration is shown schematically in Fig.4.

We define the total phase advance of the cell as $2\pi\nu$ where $\nu$ is the tune given by

[TABLE]

Here $k$ is the inverted focal length of the T-insert and $L$ is the length of the drift. The $\beta$-function is given by

[TABLE]

and the phase advance as a function of position inside the drift

[TABLE]

As a reference, we will use the initial design of the nonlinear insert from Ref.Ant_o where nonlinear magnets are placed with equal distance in space and nonlinear potential scaling with respect to the $\beta$-function was chosen to be

[TABLE]

Schematically this arrangement is shown in Fig.3(a).

From Eq.(II.1) and Eq.(II.2), it is apparent that in order to maintain accuracy and integrator structure, it is required that nonlinear elements are separated by the same phase advance $h$. Another important requirement is that $h<1$, as the error is proportional to the power of $h$. The linear stability criteria requires the tune, $\nu\leq 0.5$ or equivalently the phase advance of the cell in Fig.4 to be less than $\pi$. In order to have a step $h<1$ in the whole range of tunes, the best choice of the number of integrator steps per one period is $N=5$. This gives a reasonable upper bound for the integrator step $h\leq\pi/5\approx 0.62$.

As far as five phase steps are needed to fill the drift, we insert five nonlinear magnets according to the Ruth integration scheme of the second order with one step given by (23). With this lattice, the formula will read:

[TABLE]

with each nonlinear magnet $N_{j}$ strength scaled to the $\beta$-function as prescribed by Eq.(33) and Eq.(34). Schematically this lattice is shown in Fig.3(b). The length of each drift $O_{j}$ is calculated according to Eq.(III.1). It is worth mentioning that in order to maintain the second order of the integrator, the first and the last drifts, $O_{1}$, in the lattice should correspond to half of the phase advance, $h/2$, between the nonlinear elements.

III.1.1 Sextupole channel

For the first illustration of the Ruth lattice Fig.3(b) based on the second order integrator, we consider the sextupole magnet as a nonlinear element. This type of nonlinear insert was fist proposed in Ref.AntHH where authors suggested experiment for real world testing of regular and chaotic motion at IOTA ring based on Henon-Heiles system.

The transverse part of the nonlinear potential of a thin sextupole is given by

[TABLE]

here $a^{(3)}$ is the strength of the sextupole. The corresponding smooth Hamiltonian for this system is a well studied Henon-Heiles Hamiltonian HenHal ; Tabor

[TABLE]

A system with this Hamiltonian is not integrable and exhibits chaotic motion. Nevertheless, trajectories that correspond to the level set $\mathrm{H}>1/12$ may be chaotic but they are still bounded up to the level set $\mathrm{H}=1/6$. This is a consequence of the fact that the equipotential lines are closed up to the level set $\mathrm{H}=1/6$ Tabor .

For comparison we build a lattice where sextupoles are placed with equal distance in accordance with the previous method of discretization from Ref.Ant_o (Fig.3(a)), as well as scale the sextupole strength with the $\beta$-function according to Eq.(39) following the Ref.DanNag

[TABLE]

Here $s_{i}$ is the physical position of the $i$-th thin sextupole in the lattice.

In case of the Ruth lattice, nonlinear magnets are placed with equal phase advance between the magnets, with the first and last step being half of this phase advance. Positions in real space are calculated with Eq.(III.1). The sextupole strength scaling with the $\beta$-function, with Eq.(33), Eq.(34) reads

[TABLE]

Here, as before, $s_{i}$ - is the physical position of the $i$-th thin sextupole in the lattice. We notice that the scaling law for the Ruth lattice is different from the original scaling (42) suggested in Ref.DanNag .

To produce the comparison we chose the following parameters: Hamiltonian Eq.(41) level set $\mathrm{H}=0.1$ and the tune of the linear lattice $\nu=0.3344$ (slightly above the characteristic resonance $1/3$ of the sextupole).

First we build a Poincare surface of section $q_{2}=0$ for the equidistant lattice with the sextupole strength scaling given by (42) and compare it to the Poincare surface of a section of a smooth Hamiltonian in Fig.5. We observe that tracking of the equidistant lattice failed as most of the initial conditions rapidly escaped to infinity. Value of the Hamiltonian Eq.(41) calculated from the tracking deviates significantly from the initial value of 0.1 and rapidly grows as it could be seen from the bottom panel in Fig.5.

Next we demonstrate that the Ruth lattice indeed preserves the Hamiltonian. We build a Poincare surface of section $q_{2}=0$ for the Ruth lattice and compare it to the Poincare surface of a section of a smooth Hamiltonian in Fig.6. From the upper plots in Fig.6 it is apparent that the topology of both surfaces of the section coincide quite well and Hamiltonian (Fig.6 bottom panel) is conserved with high accuracy.

To produce Poincare surface of section in both cases (Ruth and equidistant lattice) several sets of initial conditions were tracked for $10^{6}$ iterations. Every point that was closer than $10^{-3}$ to the surface of section was projected on this surface. To produce Poincare surface of section of a smooth Hamiltonian 6th order symplectic integrator was used. Once the phase trajectory passed though the surface, a step back was made, and integration step was decreased. The procedure was repeated until the distance from the closet point to the surface of section was below $10^{-12}$. This point was considered to be located on the surface of section.

It is worth to mention that Poincare surfaces of section for a smooth system and for the Ruth lattice are not identical. Slight difference that is observed in the upper right and upper left plots Fig.6 are the result of a high value of the parameter $h\approx 0.42$ - equivalent integrator step in phase for the Ruth lattice. From Fig.6 it is clearly seen that mean value of the effective Hamiltonian in case of a lattice is lower and this results in a different mean phase space trajectory that would closer correspond to the phase trajectory of a smooth system at a lower level set.

In order to produce another head to head comparison of the Ruth lattice and the equidistant lattice, we pick a much lower level set of the Hamiltonian Eq.(41) $\mathrm{H}=0.00332$ and track both lattices for $1.5\times 10^{4}$ iterations. We again picked the tune of the linear lattice $\nu=0.3344$ slightly above the characteristic resonance $1/3$ of the sextupole for both Ruth lattice and equidistant lattice. Projection of the trajectory on the $(q_{1},q_{2})$ plane is presented in Fig.7 for both cases. We observe that in case of the equidistant lattice the projection looks irregular and spread in comparison to the Ruth lattice. We also observe that while the Ruth lattice preserves energy as expected (Fig.7 bottom left panel), for the equidistant lattice fluctuations in Hamiltonian are extremely high (Fig.7 bottom right panel).

III.1.2 Octupole channel

Next, we move to a more practical application and consider an octupole as the nonlinear element. It is stated in Ref.Ant_o that implementation of an octupole channel with one invariant of motion could provide large betatron frequency spread, and thus paves the way to overcome fast coherent beam instabilities. As a suitable invariant, per the initial suggestion in Ref.DanNag , for both experiments described in Ref.Ant_o and in Ref.Umer , the Hamiltonian in normalized coordinates was considered.

The transverse part of the potential for a thin octupole is given by

[TABLE]

The corresponding smooth system is again a Henon-Heiles type system with the Hamiltonian

[TABLE]

Using the same reasoning as for the original Henon-Heiles system, one may show that due to the fact that all equipotential lines are closed up to the level set $\mathrm{H}=1/4$ motion remain bounded, however as before it may be also chaotic.

According to the prescription of Ref.DanNag ; Ant_o , the octupole strength now should be scaled with the $\beta$-function according to Eq.(39) as

[TABLE]

and again placed with equal distance between the magnets. Here $s_{i}$ - is the physical position of the $i$-th thin octupole in the lattice. In the case of the Ruth lattice, octupole strength scaled with the $\beta$-function according to Eq.(33) and Eq.(34) as

[TABLE]

We notice that the scaling law for the Ruth lattice is different and has the same trend as in Eq.(43): the power of $\beta(s)$ is less by $1$ in comparison to the Ref.DanNag .

As before, we build a Poincare surface of section $q_{2}=0$ for the equidistant lattice and Ruth lattice by tracking several sets of initial conditions for $10^{6}$ iterations and compare it to the Poincare surface of section of a smooth Hamiltonian Eq.(45) in Fig.8 and Fig.9 correspondently. We utilized the same technique as in Sec.III.1.1 to produce these surface of sections. For the comparison we considered level set $\mathrm{H}=0.2$ with $\mathrm{H}$ defined by Eq.(45). We picked the tune of the linear lattice $\nu=0.2234$ slightly below the characteristic resonance $1/4$ of the octupole. By comparing upper left and upper right plots in Fig.9 we again observe that the topology of both surfaces of section coincide quite well and the Hamiltonian Eq.(45) (Fig.9 bottom panel) is conserved with high accuracy, as expected. Tracking of the equidistant lattice with the potential scaling given by Eq.(46) again failed as part of the considered initial conditions rapidly escaped to infinity as illustrated in Fig.8.

On can notice that as in the case of a Ruth lattice with the sextupoles Poincare surfaces of section for a smooth system and for the Ruth lattice with octupoles do not coincide completely. This again could be explained by a high value of the parameter $h\approx 0.28$ - equivalent integrator step in phase.

In order to highlight the significant difference between the equidistant lattice and Ruth lattice we track two sets of initial conditions near the separatrix $\mathrm{H}=1/4$. In Fig.10, purple color corresponds to the level set $\mathrm{H}=0.2232$ slightly below the separatrix and the black color to the level set $\mathrm{H}=0.2575$ - slightly above the separatrix. From Fig.10 we conclude, that the motion is quite different and the Hamiltonian in case of the equidistant lattice is conserved poorly. We also note that for the equidistant lattice, projection of the trajectories splits into regions and form an island-like structure. It is worth mentioning that recently such structure of the projection was observed experimentally at IOTA NK with the help of the SyncLight system NKs , though the number of nonlinear magnets used in the experiment was 17 and not 5 as in the present numerical example.

III.2 Yoshida lattice

A linear lattice that mimics the Yoshida integration scheme should correspond to the modified integrator structure given by Eq.(II.2). It is schematically depicted in Fig.3(c). As seen from Eq.(II.2), only three nonlinear elements are needed, and the physical distances between them could be quite large. Locations of the nonlinear elements are defined by the integrator coefficients, $\gamma_{1}$ and $\kappa_{1}$, and relative strength of the magnetic field by $\gamma_{1}$ and $\kappa_{2}$. Positions of the magnets also depend on the choice of $h$ - the integrator step in phase that could be an arbitrary number less then unity. The total phase advance of the linear lattice should be $4\pi+h$ as follows again from Eq.(II.2).

To illustrate the Yoshida lattice we picked the simplest example where linear optics is implemented with FODO cells without bends. For this example we treat quadrupoles as infinitely thin lenses. The phase advance of one cell is chosen as prescribed by the first rotation in Eq.(II.2) - $h\gamma_{1}/2$

For the specific example, linear optics is implemented with 30 FODO cells (see Fig.11) with the phase advance of one cell equal to $0.14\pi$. The total phase advance of the channel is $4.208\pi$ and corresponding step of the integrator (II.2) is $h=0.208\pi$.

It was suggested in Ref.DanNag to design a machine with the following effective Hamiltonian

[TABLE]

here $V_{\mathrm{Dr}}(q_{1},q_{2})$ is the Darboux potential (a solution to the Bertrand-Darboux partial differential equation) given by

[TABLE]

More details on this potential and corresponding nonlinear magnet could be found in the recent publication by Mitchell in Ref.chad_r .

Hamiltonian Eq.(48) is completely integrable and thus has infinite region of stable and regular particle motion (infinite dynamic aperture). The first integral of motion is the Hamiltonian itself and second integral of motion is given by DanNag ; Bella

[TABLE]

With $\xi(q_{1},q_{2})$ and $\eta(q_{1},q_{2})$ given as in Eq.(49) and

[TABLE]

We have considered thin nonlinear lens with the potential given by Eq.(49) as a nonlinear element for the Yoshida lattice (red blocks in Fig.11) and calculated corresponding transformations of these lenses according to the prescription of Eq.(33) and Eq.(34). Exact values of the $\beta$-function as well as the $\beta$-plot in Fig.11 were calculated using “OptiMX” software optim .

We build a Poincare surface of section $q_{1}=0$ for the Yoshida lattice by tracking several sets of initial conditions for $5\times 10^{5}$ iterations and compare it to the Poincare surface of section of a smooth Hamiltonian given by Eq.(48). In Fig.12 we show results for both the tracking of Yoshida lattice and a smooth system. For the comparison we took the level set $\mathrm{H}=0.15$. The same technique as described in Sec.III.1.1 was utilized to produce these plots.

As expected we see that the topology of the Poincare surface of section is the same, and moreover it looks like they are almost identical. We also observe that the Hamiltonian Eq.(48) as well as the second integral of motion Eq.(50) are conserved within $\approx 3\%$ accuracy.

Less amount of points on the right top plot Fig.12 in comparison to the left top plot Fig.12 is explained by a difference in methods when constructing the Poincare surface of section as well as relatively low number of iterations for the Yoshida lattice tracking.

This illustration explicitly demonstrates that the nonlinear integrable channel introduced in Ref.DanNag could be implemented with as low as just three nonlinear elements with the field that corresponds to the Darboux potential given by Eq.(49). It is worth mentioning that nonlinear magnets in the Yoshida lattice could be any nonlinear lens, be it sextupole, octupole or any other potential that could be implemented using magnetic coils.

IV Conclusion

We have presented a new approach to the nonlinear optics design strategy based on symplectic integration schemes. We demonstrated how one may utilize known symplectic integrators to produce optics configurations that preserve given Hamiltonians with any desired degree of accuracy. As relevant examples, we considered the Ruth integrator of the second order and the Yoshida integrator of the forth order to produce optics lattices.

In the presented examples we have demonstrated that the current design of the octupole and nonlinear channels under active study at the IOTA facility could be improved with just a few modifications. We showed as well that with just three nonlinear magnets one can produce a lattice that will still conserve the given Hamiltonian and thus have a large region of stable and regular particle motion.

We would like to highlight that the suggested method of lattice design is independent of the choice of nonlinear potential and thus could be applied for any nonlinear lens, be it conventional sextupole and octupole lens, or nonlinear lens with the Darboux potential DanNag . It may even be expanded to the case of electron lenses el1 ; el2 with some modifications.

The approach that we introduced may also serve as an initial seed for the conventional lattice design and could be incorporated into the lattice optimization workflow.

Appendix A Recovery of the Hamiltonian from the Ref.DanNag

In this section we recover original scaling of the nonlinear potential with the $\beta$-function from the Ref.DanNag . First we observe that if $h\to 0$ and $\beta_{x}=\beta_{y}=\beta$ then substitution

[TABLE]

is exact and the flow that is given by Eq.(33) could be rewritten as

[TABLE]

with the additional substitution for the $\partial_{x,y}U$ defined in Eq.(34).

It is apparent that Hamiltonian that corresponds to this flow (53) is exactly

[TABLE]

As far as the flow that is given by Eq.(II.4) is generated by a Hamiltonian

[TABLE]

the final Hamiltonian that generates the flow

[TABLE]

in the case of $h\to 0$ is exactly the one from the Ref.DanNag and is given by

[TABLE]

It is worth to mention that differences in scaling coefficient that were mentioned in Sec.III.1.1 and Sec.III.1.2 vanish in the limit of $h\to 0$. Indeed, as mentioned above, substitution $h=ds/\beta$ is exact in this case and thus there is no difference in using $h$ or $ds$ in therms of an integrator. However, switching to $ds$ introduces additional factor $1/\beta$ as one may notice from Eq.(53).

Appendix B Cross-check with the MADX

To produce such a comparison we consider Yoshida lattices that is described in Sec.III and schematically presented in Fig.11 - 30 FODO cells with the phase advance of one cell equal to $0.14\pi$. The total phase advance of the channel is $4.208\pi$ and corresponding step of the integrator (II.2) is $h=0.208\pi$. A sextupole was picked as a nonlinear element instead of a magnet with the Darboux potential. Consequently, the invariant of motion that this lattice is aimed to conserve is a Henon-Heiles Hamiltonian in normalized coordinates that is given by Eq.(41).

In Fig.13 we present Poincare surface of sections ($q_{2}=0$) for the smooth Hamiltonian Eq.(41) and Yoshida lattice with sextupoles. We observe, that as before lattice indeed conserves given Hamiltonian as prescribed. For the MADX input initial conditions that are listed in Table 1 were transformed to the real coordinate space with the help of the betatron amplitude matrix given by Eq.(II.4) and linear lattice functions that were calculated initially using OptiMX optim and then confirmed by MADX.

MADX input file as well as Wolfram Mathematica post-processing script could be found in supplemental materials.

Appendix C Definition of some terms used in the paper

Bounded Motion - the motion of is called bounded if for the initial conditions $X_{0}$ and a state vector $X$ there exist some constant $C$ such that $|X|\leq C$.

Dynamic Aperture - maximum phase space volume that contains bounded trajectories of the system.

Chaotic motion - a random motion of a deterministic system.

Integrable system - a Hamiltonian system with $2n$ dimensional phase space that has the maximal number of independent Poisson commuting invariants (including the Hamiltonian itself) equals to $n$.

Quasi-integrable system - a Hamiltonian system that has approximate Poisson commuting invariants that are close to the invariants of some integrable system.

Smooth Hamiltonian - a Hamiltonian of a system with continuous time.

Effective Hamiltonian - a smooth Hamiltonian that captures the dynamics of a discreet system.

Appendix D Parameters for the tracking

In this Appendix in order for the reader to be able to reproduce tracking results we present sets of initial conditions that were usec to produce Fig.5,Fig.6, Fig.8, Fig.9 and Fig.12. Note that coordinates are given in normalized (canonical) coordinates that are connected to the real coordinates through the betatron amplitude matrix SYL .

Acknowledgements.

This work was supported by the U.S. National Science Foundation under Award No. PHY-1549132, the Center for Bright Beams and under Award No. PHY-1535639; The author is grateful to Sergei Nagaitsev for the introduction to the problem and several fruitful discussions, to Alexander Valishev for couple of useful remarks, to Ivan Bazarov for sharing Ref.GInt , to Gerard Andonian for reading the manuscript.