\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{a4margins}

\newcommand{\dd}[2]{\frac{\mathrm d #1}{\mathrm d #2}}
\newcommand{\dds}[2]{\frac{\mathrm d^2 #1}{{\mathrm d #2}^2}}
\newcommand{\pdd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pdds}[2]{\frac{\partial^2 #1}{\partial{#2}^2}}
\newcommand{\pddd}[3]{\frac{\partial^2 #1}{\partial #2\partial #3}}

\newcommand{\sq}[1]{\left(#1\right)^2}

\begin{document}

\title{Fixing the Tune of Nonscaling FFAGs in the Thin Lens Paraxial Approximation}
\author{Stephen Brooks}
\maketitle

\section{Dynamics}\label{-dynamics}

The beamline consists of drifts of length $d_n$ for $1\le n\le N$, each followed by magnetic kicks defined by $b_n(x)$.  The horizontal phase space position before the first drift is $(x_0,x'_0)$, with later positions given by
\begin{eqnarray*}
x_n &=& x_{n-1} + d_n x'_{n-1}\\
x'_n &=& x'_{n-1} + p^{-1} b_n(x_n).
\end{eqnarray*}

In physical quantities, $b = B\ell q$ or alternatively $B\ell q/p_0$ if using normalised $p$ with $p=1$ meaning the reference momentum $p_0$.  $B$ here is the from the magnetic component $B_y(x)$.

\newcommand{\dy}{\delta y}
Vertical dynamics can be analysed as small perturbations $\dy$ from the midplane, incorporating the consequence of $\nabla \cdot \mathbf B=0$:
\begin{eqnarray*}
\delta y_n &=& \dy_{n-1} + d_n \dy'_{n-1} \\
\dy'_n &=& \dy'_{n-1} - p^{-1} b'_n(x_n) \dy_n.
\end{eqnarray*}

If you consider the final positions as functions $x_N = x_N(x_0,x'_0,p)$, $x'_N=x'_N(x_0,x'_0,p)$ and similarly for $\dy_N$, the phase advances are given by

\begin{eqnarray*}
2 \cos \phi_x &=& \pdd{x_N}{x_0}+\pdd{x'_N}{x'_0} \\
2 \cos \phi_y &=& \pdd{\dy_N}{\dy_0}+\pdd{\dy'_N}{\dy'_0}.
\end{eqnarray*}

\section{Closed Orbit Evolution}\label{-closed}

The closed orbit position $(x_c,x'_c)$ at the beginning and end of the beamline (for momentum $p$) is defined by
\begin{eqnarray*}
x_N(x_c,x'_c,p) &=& x_c \\
x'_N(x_c,x'_c,p) &=& x'_c,
\end{eqnarray*}
with $\dy_c=\dy'_c=0$ by construction because only $B_y$ fields are present.  Because the $b_n$ are general nonlinear functions, these equations must be solved numerically, but once the closed orbit is known, its variation with momentum may be found by differentiation:
\begin{eqnarray*}
\pdd{x_N}{x_0}\dd{x_c}{p} + \pdd{x_N}{x'_0}\dd{x'_c}{p} + \pdd{x_N}{p} &=& \dd{x_c}{p} \\
\pdd{x'_N}{x_0}\dd{x_c}{p} + \pdd{x'_N}{x'_0}\dd{x'_c}{p} + \pdd{x'_N}{p}  &=& \dd{x'_c}{p},
\end{eqnarray*}
where all partial derivatives are evaluated at $(x_c,x'_c,p)$.  This forms a $2\times 2$ matrix system
\newcommand{\mtwo}[4]{\left[\begin{array}{cc}#1&#2\\#3&#4\end{array}\right]}
\newcommand{\vectwo}[2]{\left[\begin{array}{c}#1\\#2\end{array}\right]}
\[\mtwo{\pdd{x_N}{x_0}-1}{\pdd{x_N}{x'_0}}
{\pdd{x'_N}{x_0}}{\pdd{x'_N}{x'_0}-1}  \vectwo{\dd{x_c}{p}}{\dd{x'_c}{p}} = \vectwo{-\pdd{x_N}{p}}{-\pdd{x'_N}{p}}\]
in terms of the derivatives of the beamline functions $x_N$ and $x'_N$.  Here it is useful to note the inverse of a $2\times 2$ matrix has a relatively simple form
\[\mtwo{a}{b}{c}{d}^{-1} = \frac{1}{D}\mtwo{d}{-b}{-c}{a}\]
where $D=ad-bc$ is the determinant, here
\[ D = \left(\pdd{x_N}{x_0}-1\right)\left(\pdd{x'_N}{x'_0}-1\right) - \pdd{x_N}{x'_0}\pdd{x'_N}{x_0}. \]
This can be simplified by recalling that the overall mapping $(x_0,x'_0)\mapsto(x_N,x'_N)$ preserves phase space area, so its Jacobian determinant is unity:
\[ 1 = \pdd{x_N}{x_0}\pdd{x'_N}{x'_0} - \pdd{x_N}{x'_0}\pdd{x'_N}{x_0}. \]
This allows the product terms in $D$ to be replaced by $1$, leaving only linear terms:
\begin{eqnarray*}
D &=& 2 - \pdd{x_N}{x_0} - \pdd{x'_N}{x'_0} \\
&=& 2 - 2 \cos \phi_x,
\end{eqnarray*}
recalling the definition of the horizontal phase advance.  This means the closed orbit movement becomes singular ($D=0$) only when $\cos \phi_x =1$, that is $2\pi n$ horizontal phase advance, also called `integer tune'.

Applying the inverse matrix gives the orbit evolution explicitly as
\begin{eqnarray*}
\dd{x_c}{p} &=& \frac{1}{D}\left(\left(\pdd{x'_N}{x'_0}-1\right)\left(-\pdd{x_N}{p}\right)-\pdd{x_N}{x'_0}\left(-\pdd{x'_N}{p}\right)\right) \\
&=& \frac{1}{D}\left(\left(1-\pdd{x'_N}{x'_0}\right)\pdd{x_N}{p}+\pdd{x_N}{x'_0}\pdd{x'_N}{p}\right) \\
\dd{x'_c}{p} &=& \frac{1}{D}\left(\pdd{x'_N}{x_0}\pdd{x_N}{p}+\left(1-\pdd{x_N}{x_0}\right)\pdd{x'_N}{p}\right).
\end{eqnarray*}

The required partial derivatives can be evaluated starting with
\[\pdd{x_0}{x_0}=1 \qquad \pdd{x_0}{x'_0}=0 \qquad \pdd{x_0}{p}=0 \qquad 
\qquad \pdd{x'_0}{x_0}=0 \qquad \pdd{x'_0}{x'_0}=1 \qquad \pdd{x'_0}{p}=0,\]
and proceeding for $n\ge 1$ via the chain rule:
\begin{eqnarray*}
\pdd{x_n}{x_0} &=& \pdd{x_{n-1}}{x_0} + d_n \pdd{x'_{n-1}}{x_0} \qquad \mbox{(similarly for $x'_0$ and $p$)}\\
\pdd{x'_n}{x_0} &=& \pdd{x'_{n-1}}{x_0} + p^{-1} b'_n(x_n)\pdd{x_n}{x_0} \qquad \mbox{(similarly for $x'_0$)}\\
\pdd{x'_n}{p} &=& \pdd{x'_{n-1}}{p} - p^{-2} b_n(x_n) + p^{-1} b'_n(x_n)\pdd{x_n}{p}.
\end{eqnarray*}
For brevity it is useful to introduce the (de)focussing strengths $k_n=p^{-1} b'_n(x_n)$ and the quadratic term in inverse momentum $q_n= - p^{-2} b_n(x_n)$ to remove minus signs.

\subsection{$N=2$ case}
The simplest cell with stable focussing in both planes contains two kicks.  The expressions above can be evaluated more explicitly.
\[ \pdd{x_1}{x_0} = 1 \qquad
\pdd{x_1}{x'_0} = d_1 \qquad
\pdd{x_1}{p} = 0 \]
\[ \pdd{x'_1}{x_0} = k_1 \qquad
\pdd{x'_1}{x'_0} = 1 + k_1d_1 \qquad
\pdd{x'_1}{p} = q_1 \]
\[ \pdd{x_2}{x_0} = 1 + d_2 k_1 \qquad
\pdd{x_2}{x'_0} = d_1 + d_2(1 + k_1d_1) \qquad
\pdd{x_2}{p} = d_2 q_1 \]
\begin{eqnarray*}
\pdd{x'_2}{x_0} &=& k_1 + k_2(1 + d_2 k_1) \\
\pdd{x'_2}{x'_0} &=& 1 + k_1d_1 + k_2(d_1 + d_2(1 + k_1d_1)) \\
\pdd{x'_2}{p} &=&  q_1 + q_2 + k_2 d_2 q_1.
\end{eqnarray*}
The $x_2$ and $x'_2$ results from here can be put into the $2\times 2$ matrix system for the closed orbit movement with $p$.  The determinant of the matrix is
\begin{eqnarray*}
D &=&  2 - \pdd{x_2}{x_0} - \pdd{x'_2}{x'_0} \\
&=& 2 - (1 + d_2 k_1) - (1 + k_1d_1 + k_2(d_1 + d_2(1 + k_1d_1))) \\
&=& -d_2k_1-k_1d_1-k_2d_1-k_2d_2-k_2d_2k_1d_1 \\
&=&  - (d_1+d_2)(k_1+k_2) - d_1d_2k_1k_2.
\end{eqnarray*}

Now the expressions for the closed orbit movement itself:
\begin{eqnarray*}
\dd{x_c}{p} &=& \frac{1}{D}\left(\left(1-\pdd{x'_2}{x'_0}\right)\pdd{x_2}{p}+\pdd{x_2}{x'_0}\pdd{x'_2}{p}\right) \\
&=& D^{-1}(-(k_1d_1 + k_2(d_1 + d_2(1 + k_1d_1)))d_2 q_1 + (d_1 + d_2(1 + k_1d_1))(q_1 + q_2 + k_2 d_2 q_1)) \\ 
&=& D^{-1}(-k_1d_1d_2 q_1 + (d_1 + d_2(1 + k_1d_1))(-k_2d_2 q_1 + q_1 + q_2 + k_2 d_2 q_1)) \\
&=& D^{-1}((d_1+d_2)(q_1+q_2) + d_1d_2k_1q_2) \\
\dd{x'_c}{p} &=& \frac{1}{D}\left(\pdd{x'_2}{x_0}\pdd{x_2}{p}+\left(1-\pdd{x_2}{x_0}\right)\pdd{x'_2}{p}\right) \\
&=& D^{-1}((k_1 + k_2(1 + d_2 k_1))d_2 q_1 - d_2 k_1(q_1 + q_2 + k_2 d_2 q_1)) \\
&=& D^{-1}d_2(k_1q_1 + k_2q_1 + k_2d_2k_1q_1 - k_1q_1 - k_1q_2 - k_1k_2d_2q_1) \\
&=& D^{-1}d_2(k_2q_1 - k_1q_2).
\end{eqnarray*}

\section{Phase Advance Evolution}\label{-xphase}

To find the change of the phase advances as momentum varies, it is important to remember that the closed orbit where they are evaluated is also moving.  This gives the total derivative in $p$ the following form:
\begin{eqnarray*}
\dd{}{p}(2 \cos \phi_x) &=& \dd{}{p}\pdd{x_N}{x_0}+\dd{}{p}\pdd{x'_N}{x'_0} \\
&=& \pdds{x_N}{x_0}\dd{x_c}{p}+\pddd{x_N}{x_0}{x'_0}\dd{x'_c}{p}+\pddd{x_N}{x_0}{p}
+\pddd{x'_N}{x_0}{x'_0}\dd{x_c}{p}+\pdds{x'_N}{x'_0}\dd{x'_c}{p}+\pddd{x'_N}{x'_0}{p} \\
&=& \left(\pdds{x_N}{x_0}+\pddd{x'_N}{x_0}{x'_0}\right)\dd{x_c}{p} + 
\left(\pddd{x_N}{x_0}{x'_0}+\pdds{x'_N}{x'_0}\right)\dd{x'_c}{p}+\pddd{x_N}{x_0}{p}+\pddd{x'_N}{x'_0}{p}.
\end{eqnarray*}
Expanding the expressions for the closed orbit derivatives does not appear to yield additional simplification without substituting a specific value for $N$.
%Testing
%\[ \left(\pdds{x_N}{x_0}+\pddd{x'_N}{x_0}{x'_0}\right)\dd{x_c}{p}
%= \left(\pdds{x_N}{x_0}+\pddd{x'_N}{x_0}{x'_0}\right)\frac{1}{D}\left(\left(1-\pdd{x'_N}{x'_0}\right)\pdd{x_N}{p}+\pdd{x_N}{x'_0}\pdd{x'_N}{p}\right) \]

Now the second partial derivatives of the beamline functions are required.  For $n=0$ they are all zero and for $n\ge 1$ the chain rule gives:
\begin{eqnarray*}
\pdds{x_n}{x_0} &=& \pdds{x_{n-1}}{x_0} + d_n \pdds{x'_{n-1}}{x_0} \qquad \mbox{(similarly for all $\partial\alpha\partial\beta$)}\\
\pdds{x'_n}{x_0} &=& \pdds{x'_{n-1}}{x_0} + k'_n\sq{\pdd{x_n}{x_0}} + k_n\pdds{x_n}{x_0} \qquad \mbox{(similarly for $\alpha,\beta\in\{x_0,x'_0\}$)}\\
% \pdd{x'_n}{p} &=& \pdd{x'_{n-1}}{p} + q_n + k_n\pdd{x_n}{p}.
\pddd{x'_n}{x_0}{p} &=& \pddd{x'_{n-1}}{x_0}{p} + q'_n\pdd{x_n}{x_0} + k'_n\pdd{x_n}{x_0}\pdd{x_n}{p} + k_n\pddd{x_n}{x_0}{p} \qquad \mbox{(similarly for $\partial x_0'$)}.
\end{eqnarray*}
Here the quantities $k'_n=p^{-1} b''_n(x_n)$ and $q'_n= - p^{-2} b'_n(x_n)$ are the derivatives with respect to $x_n$ of the unprimed versions.  These second derivatives only achieve a nonzero value through coupling from the first-order derivatives.  N.B. $q'_n=-k_n/p$ so $k_iq'_j=q'_ik_j$.

\subsection{$N=2$ case}\label{-xphase-two}
\newcommand{\xe}{x^\cdot_0}
The second derivatives are evaluated below in the case of the two-element cell.  This calculation and the recurrence relation that generated it have been implemented as a Mathematica workbook to check the formulae below.
\[\pddd{x_1}{\alpha}{\beta}=0 \]
\[ \pdds{x'_1}{x_0}=k'_1 \qquad
\pddd{x'_1}{x_0}{x'_0}=k'_1d_1 \qquad
\pdds{x'_1}{x'_0}=k'_1d_1^2 \qquad
\pddd{x'_1}{x_0}{p}=q'_1 \qquad
\pddd{x'_1}{x'_0}{p}=q'_1d_1 \]
\[ \pdds{x_2}{x_0}=d_2k'_1 \qquad
\pddd{x_2}{x_0}{x'_0}=d_2k'_1d_1 \qquad
\pdds{x_2}{x'_0}=d_2k'_1d_1^2 \qquad
\pddd{x_2}{x_0}{p}=d_2q'_1 \qquad
\pddd{x_2}{x'_0}{p}=d_2q'_1d_1 \]
\begin{eqnarray*}
% \pdds{x'_{n-1}}{x_0} + k'_n\sq{\pdd{x_n}{x_0}} + k_n\pdds{x_n}{x_0}
\pdds{x'_2}{x_0} &=& k'_1 + k'_2\sq{1 + d_2 k_1} + k_2d_2k'_1 \\
\pddd{x'_2}{x_0}{x'_0} &=& k'_1d_1 + k'_2(1 + d_2 k_1)(d_1 + d_2(1 + k_1d_1)) + k_2d_2k'_1d_1 \\
\pdds{x'_2}{x'_0} &=& k'_1d_1^2 + k'_2\sq{d_1 + d_2(1 + k_1d_1)} +  k_2d_2k'_1d_1^2 \\
% \pddd{x'_n}{x_0}{p} &=& \pddd{x'_{n-1}}{x_0}{p} + q'_n\pdd{x_n}{x_0} + k'_n\pdd{x_n}{x_0}\pdd{x_n}{p} + k_n\pddd{x_n}{x_0}{p}
\pddd{x'_2}{x_0}{p} &=& q'_1 + q'_2(1 + d_2 k_1) + k'_2(1 + d_2 k_1)d_2q_1 + k_2d_2q'_1 \\
\pddd{x'_2}{x'_0}{p} &=& q'_1d_1 + q'_2(d_1 + d_2(1 + k_1d_1)) + k'_2(d_1 + d_2(1 + k_1d_1))d_2q_1 + k_2d_2q'_1d_1
\end{eqnarray*}
These may be used to evalulate the change in horizontal phase advance:
\begin{eqnarray*}
\dd{}{p}(2 \cos \phi_x) &=& \left(\pdds{x_2}{x_0}+\pddd{x'_2}{x_0}{x'_0}\right)\dd{x_c}{p} + 
\left(\pddd{x_2}{x_0}{x'_0}+\pdds{x'_2}{x'_0}\right)\dd{x'_c}{p}+\pddd{x_2}{x_0}{p}+\pddd{x'_2}{x'_0}{p} \\
&=& \left(d_2k'_1 +  k'_1d_1 + k'_2(1 + d_2 k_1)(d_1 + d_2(1 + k_1d_1)) + k_2d_2k'_1d_1\right) \dd{x_c}{p} + \\
&& \left(d_2k'_1d_1 + k'_1d_1^2 + k'_2\sq{d_1 + d_2(1 + k_1d_1)} +  k_2d_2k'_1d_1^2\right) \dd{x'_c}{p} + \\
&& d_2q'_1 + q'_1d_1 + q'_2(d_1 + d_2(1 + k_1d_1)) + k'_2(d_1 + d_2(1 + k_1d_1))d_2q_1 + k_2d_2q'_1d_1.
\end{eqnarray*}
At this point, further simplifications are possible if the $k^{(\prime)}_n$ and $q^{(\prime)}_n$ terms are replaced with their definitions in terms of derivatives  of $b_n(x_n)$ and the beam momentum $p$.  Mathematica gives:
\[
\dd{}{p}(2 \cos \phi_x) =
\frac{
\begin{array}{l}
d_1^2 \left(\left(b_2 b''_2-2 b'^2_2\right) b'^2_1+b_1 b'^2_2 b''_1\right) d_2^2 \\ - p d_1 (d_1+d_2) \left(3 b'_2 b'^2_1+\left(3 b'^2_2-(b_1+2 b_2) b''_2\right) b'_1-(2 b_1+b_2) b'_2 b''_1\right) d_2 \\ - p^2 (d_1+d_2)^2 \left(b'^2_1+2 b'_2 b'_1+b'^2_2-(b_1+b_2) \left(b''_1+b''_2\right)\right)
\end{array}
}
%/
{p^3 \left(d_1 d_2 b'_1 b'_2+p (d_1+d_2) \left(b'_1+b'_2\right)\right)},
\]
with the form of $D$ (and some powers of $p$) appearing on the denominator.  The second derivatives $b''_n$ are only linearly involved in this expression and may be separated out:
\begin{eqnarray*}
p^5D\dd{}{p}(2 \cos \phi_x) &=&
2 d_1^2 d_2^2 b_1^{\prime }{}^2 b_2^{\prime }{}^2+3 p d_1 d_2 \left(d_1+d_2\right) b_1^{\prime } \left(b_1^{\prime }+b_2^{\prime }\right) b_2^{\prime }+p^2 \left(d_1+d_2\right){}^2 \left(b_1^{\prime }+b_2^{\prime }\right){}^2\\&&-b_1^{\prime \prime } \left(p \left(d_1+d_2\right)+d_1 d_2 b_2^{\prime }\right) \left(p \left(b_1+b_2\right) \left(d_1+d_2\right)+b_1 d_1 d_2 b_2^{\prime }\right)\\&&-b_2^{\prime \prime }\left(p \left(d_1+d_2\right)+d_1 d_2 b_1^{\prime }\right) \left(p \left(b_1+b_2\right) \left(d_1+d_2\right)+b_2 d_1 d_2 b_1^{\prime }\right).
\end{eqnarray*}
In a situation where $\phi_x$ is required to be constant, the left-hand side is zero (and $D\ne0$ by construction), so this provides one linear constraint on $\{b''_1(x_1),b''_2(x_2)\}$ in terms of lower derivatives of $b_n$ at those locations and $p$.

\section{Vertical Phase Advance}\label{-yphase}
The closed orbit is $\dy_c=\dy'_c=0$ in the vertical plane by construction (there are no vertical magnetic deflections for $\dy=0$) but the phase advance given by
\[2 \cos \phi_y = \pdd{\dy_N}{\dy_0}+\pdd{\dy'_N}{\dy'_0} \]
must still be calculated in terms of the derivatives of the beamline function for the vertical plane.

The recurrence relation starts with
\[ \pdd{\dy_0}{\dy_0}=1 \qquad \pdd{\dy_0}{\dy'_0}=0 \qquad 
\qquad \pdd{\dy'_0}{\dy_0}=0 \qquad \pdd{\dy'_0}{\dy'_0}=1 \]
and for $n\ge 1$ the chain rule gives:
\begin{eqnarray*}
\pdd{\dy_n}{\dy_0} &=& \pdd{\dy_{n-1}}{\dy_0} + d_n \pdd{\dy'_{n-1}}{\dy_0} \qquad \mbox{(similarly for $\dy'_0$)}\\
\pdd{\dy'_n}{\dy_0} &=& \pdd{\dy'_{n-1}}{\dy_0} - p^{-1} b'_n(x_n)\pdd{\dy_n}{\dy_0} \qquad \mbox{(similarly for $\dy'_0$)}\\
&=& \pdd{\dy'_{n-1}}{\dy_0} - k_n\pdd{\dy_n}{\dy_0}.
\end{eqnarray*}
This is a completely linear system, which can be expressed in terms of the product of $2\times 2$ matrices containing $d_n$ and $k_n$.

The evolution of $2 \cos \phi_y$ with momentum starts simpler than the horizontal case because the position at which derivatives of $\dy^{(\prime)}_N(\dy_0,\dy'_0,p)$ are evaluated does not change with $p$:
\[ \dd{}{p}(2 \cos \phi_y) = \pddd{\dy_N}{\dy_0}{p} + \pddd{\dy'_N}{\dy'_0}{p}. \]
The values of $k_n$, however, do change with $p$ as the closed orbit shifts horizontally:
\begin{eqnarray*}
\pddd{\dy_n}{\dy_0}{p} &=& \pddd{\dy_{n-1}}{\dy_0}{p} + d_n \pddd{\dy'_{n-1}}{\dy_0}{p} \qquad \mbox{(similarly for $\dy'_0$)}\\
\pddd{\dy'_n}{\dy_0}{p} &=& \pddd{\dy'_{n-1}}{\dy_0}{p} - \dd{k_n}{p}\pdd{\dy_n}{\dy_0} - k_n\pddd{\dy_n}{\dy_0}{p} \qquad \mbox{(similarly for $\dy'_0$)}\\
\dd{k_n}{p} &=& -p^{-2}b'_n(x_n) + p^{-1}b''_n(x_n)\dd{x_n}{p} \ =\  q'_n + k'_n\dd{x_n}{p} \\
\dd{x_n}{p} &=& \pdd{x_n}{x_0}\dd{x_c}{p} + \pdd{x_n}{x'_0}\dd{x'_c}{p} + \pdd{x_n}{p}.
\end{eqnarray*}
The total derivative of $x_n$ appears for the first time because the vertical equations assume the horizontal closed orbit motion at each magnetic kick has already been determined.  Note that the closed orbit derivatives only appear multiplied by first derivatives of $\dy$, so fortunately will not become squared or cubed when the recurrence is evaluated.

\subsection{$N=2$ case}\label{-yphase-two}
The calculation has similar complexity to the one for the horizontal plane and was done in Mathematica.  The end result with second order derivatives of $b$ factored out is:
\begin{eqnarray*}
p^5D\dd{}{p}(2 \cos \phi_y) &=&
-\left(d_1 d_2 b'_1 b'_2+p (d_1+d_2) (b'_1+b'_2)\right) \left(p ((d_2+1) b'_1+b'_2)-2 d_2 b'_1 b'_2\right)\\
&&+ b''_1 \left(d_2 p+p-d_2 b'_2\right) \left(p (b_1+b_2) (d_1+d_2)+b_1 d_1 d_2 b'_2\right)\\
&&+b''_2 \left(p-d_2 b'_1\right) \left(p (b_1+b_2) (d_1+d_2)+b_2 d_1 d_2 b'_1\right).
\end{eqnarray*}

\section{Tune-Fixing Equations}\label{-fixing}
Having expressed the tune variations in the form
\[ \dd{}{p}(2 \cos \phi_x)=c_x-\sum_{n=1}^N a_{xn}b''_n(x_n) \qquad
\dd{}{p}(2 \cos \phi_y)=c_y-\sum_{n=1}^N a_{yn}b''_n(x_n), \]
where $c_{(x,y)}$ and $a_{(x,y)n}$ depend on $b_i(x_i)$ and $b'_i(x_i)$ only, setting the tune variations to zero produces a system with two linear constraints on the $b''_n(x_n)$:
\[
\left[
\begin{array}{cccc}
a_{x1}&a_{x2}&\cdots&a_{xN}\\
a_{y1}&a_{y2}&\cdots&a_{yN}
\end{array}
\right]
\left[
\begin{array}{c}
b''_1(x_1)\\
b''_2(x_2)\\
\vdots\\
b''_N(x_N)
\end{array}
\right]
=
\left[
\begin{array}{c}
c_x \\ c_y
\end{array}
\right].
\]

As the $b''_n$ are evaluated at different points, this system of coupled second-order differential equations cannot be solved as usual with $x$ as the independent variable.  Choosing to integrate with respect to $p$ allows the field derivatives to be expressed as
\[ b'_n(x_n(p)) = \dd{b_n}{p}\left(\dd{x_n}{p}\right)^{-1} = 
\dd{b_n}{p}\left(\pdd{x_n}{x_0}\dd{x_c}{p} + \pdd{x_n}{x'_0}\dd{x'_c}{p} + \pdd{x_n}{p}\right)^{-1}, \]
where the closed orbit derivatives can be found from only first order derivatives of $b$ and
\[ b''_n(x_n(p)) = \dd{b'_n}{p}\left(\dd{x_n}{p}\right)^{-1}, \]
which does not need to be expanded further provided $\dd{x_n}{p}$ is being calculated as the integration progresses.

\subsection{Initial Conditions}
The integration will start from a particular momentum $p=p_0$ and since the equation is second order, both $b_{0n}=b_n(x_n(p_0))$ and $b'_{0n}=b'_n(x_n(p_0))$ need to be defined.  Also required are closed orbit coordinates $(x_c(p_0),x'_c(p_0))$ that are consistent with the lattice at that momentum.  With the magnetic kicks fixed by assumption, for $1\le n\le N$ the dynamics satisfy
\begin{eqnarray*}
x_n &=& x_{n-1} + d_n x'_{n-1}\\
x'_n &=& x'_{n-1} + p_0^{-1} b_{0n},
\end{eqnarray*}
therefore
\begin{eqnarray*}
x'_n &=& x'_0 + p_0^{-1} \sum_{i=1}^n b_{0i} \\
x_n &=& x_0 + \sum_{i=1}^{n} d_i x'_{i-1} \\
&=& x_0 + \sum_{i=1}^{n} d_i \left(x'_0 + p_0^{-1} \sum_{j=1}^{i-1} b_{0j}\right) \\
&=& x_0 + x'_0 \sum_{i=1}^{n} d_i + p_0^{-1} \sum_{i=1}^{n-1} b_{0i} \sum_{j=i+1}^{n} d_j.
\end{eqnarray*}
If the remaining distances from kick $0\le n<N$ to the end of the beamline ($n=0$ being the start) are written as
\[ \Delta_n = \sum_{i=n+1}^{N}d_i, \]
then the conditions on the closed orbit (setting $x^{(\prime)}_0=x^{(\prime)}_N=x^{(\prime)}_{0c}$) are:
\begin{eqnarray*}
0 &=& \sum_{n=1}^N b_{0n} \\
0 &=& x'_{0c} \Delta_0 + p_0^{-1} \sum_{n=1}^{N-1} b_{0n} \Delta_n.
\end{eqnarray*}
The first of these says the magnetic kicks have to sum to zero to keep $x'$ the same as it started the beamline, while the second defines a value of $x'_{0c}$ that ensures there is no $x$ offset.  This leaves the $b_{0n}$ with $N-1$ free parameters, the $b'_{0n}$ with $N$ and $x_{0c}$ providing $1$ additional free parameter, for a total of $2N$.  The choice of $p_0$ is not counted as a free parameter because it could be used to pick a different integration starting point on exactly the same solution.

\subsection{Integration}
At any given value of $p$, the vector of $2N+2$ integration variables $(x_c,x'_c,b_n,b'_n)$ will need to have its first derivative with respect to $p$ calculated.  This can be done with the steps below.
\begin{enumerate}
\item Use the recurrence relation at the end of section \ref{-closed} to evaluate $\pdd{(x_n,x'_n)}{(x_0,x'_0,p)}$ for $n\le N$.
\item Use the other formulae in section \ref{-closed} to calculate $D$ and then $\dd{x_c}{p}$ and $\dd{x'_c}{p}$.  Because $\phi_x$ is supposed to be fixed, $D$ should be a constant to within the accuracy of numerical integration.
\item Calculate $\dd{x_n}{p}$ for all $n$ using the formula in section \ref{-yphase}.
\item The required derivatives $\dd{b_n}{p}$ are now given by $b'_n \dd{x_n}{p}$.
\item\label{i-cacoeffs} Evaluate the coefficients $c_{(x,y)}$ and $a_{(x,y)n}$ using the formulae at the end of sections \ref{-xphase-two} and \ref{-yphase-two}, or similar Mathematica output for $N>2$.
\item Using $N-2$ arbitrary functions of $p$ (if necessary), find a specific solution to the linear system for $b''_n$ given in section \ref{-fixing}.
\item The remaining derivatives $\dd{b'_n}{p}$ are given by $b''_n \dd{x_n}{p}$.
\end{enumerate}

Overall, then, the integration requires an initial momentum $p_0$ (does not affect the solution), $2N$ independent initial values and $N-2$ arbitrary functions of $p$ during integration to supply the degree of choice left in system for $b''_n$ that is underconstrained for $N\ge 3$.  Since $D$ is constant, the integer tune situation will not limit the range of integration, but two other situations may occur.  Firstly, if the $N$-vectors $[a_{xn}]$ and $[a_{yn}]$ become scalar multiples of each other, then the linear system for the $b''_n$ may not have a solution.  Secondly, if $\dd{x_n}{p}$ changes sign for any $n$, then the solution is no longer physical if different magnetic field strengths $b_n(x_n)$ are now being defined at points that $x_n(p)$ has already visited.

In step \ref{i-cacoeffs}, it may be possible to evaluate most of the recurrence for the second order derivatives with the current numerical values of $b^{(\prime)}_n$, but as the $b''_n$ are unknowns, the coefficients of those (appearing from the $k'_n$ terms) will have to be separated out.  This would allow the solution for higher numbers ($N$) of beamline elements without using Mathematica or having the resulting combinatorial explosion in the size of the equations.

\end{document}