\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{verbatim,a4margins,amsmath,amssymb,bm,hyperref,graphicx}

\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{\rowfour}[4]{\left[\begin{array}{cccc}#1&#2&#3&#4\end{array}\right]}
\newcommand{\vectwo}[2]{\left[\begin{array}{c}#1\\#2\end{array}\right]}
\newcommand{\mxtwo}[4]{\left[\begin{array}{cc}#1&#2\\#3&#4\end{array}\right]}

\newcommand{\sq}[1]{\left(#1\right)^2}
\newcommand{\floor}[1]{\left\lfloor#1\right\rfloor}
\newcommand{\abs}[1]{\left|#1\right|}

\newcommand{\rem}[1]{}

\newcommand{\A}{\mathbf A}
\newcommand{\E}{\mathbf E}
\newcommand{\B}{\mathbf B}
\newcommand{\J}{\mathbf J}
\newcommand{\0}{\mathbf 0}

\newcommand{\f}{\mathbf f}
\newcommand{\uv}{\mathbf u}
\newcommand{\vv}{\mathbf v}
\newcommand{\w}{\mathbf w}
\newcommand{\x}{\mathbf x}
\newcommand{\y}{\mathbf y}
\newcommand{\z}{\mathbf z}

\newcommand{\s}{\mathbf s}
\newcommand{\rv}{\mathbf r}
\newcommand{\cv}{\mathbf c}

\newcommand{\e}{\mathbf e}
\newcommand{\m}{\mathbf m}
\newcommand{\M}{\mathbf M}

\newcommand{\de}{\,\mathrm d}

\newcommand{\ii}{\mathrm i}

\begin{document}

\title{Effective Strength of Sinusoidally-Varying Focussing}
\author{Stephen Brooks}
\date{April 4, 2024}
\maketitle

\section{Introduction}
Particle beam focussing and other confined physical systems often obey the equation $\ddot x=kx$, where $k<0$ implies stable motion.  A constant negative value of $k$ yields simple harmonic motion (constant focussing), while positive values give unstable (defocussing) motion.  Alternating-gradient focussing is the surprising principle that if $k(t)$ is varied periodically, the overall motion can sometimes be stable even if the average value of $k$ is zero or negative.  It is used where focussing in one direction implies defocussing in one or more other directions, such as with electromagnetic fields in free space that have to satisfy $\nabla\cdot\E=0$ and similar.

\section{Equation of Motion and Time Rescaling}
Sinusoidally-varying focussing is a natural choice in some experiments (Paul-type ion traps for example).  The equation
\[ \ddot x=k\sin(\omega t)x \]
for some constant $k$ determines the motion.  Rescaling time with $y(t)=x(\frac{t}{\omega})$ and $\ddot y(t)=\frac1{\omega^2}\ddot x(\frac{t}{\omega})$ and evaluating the original equation at time $\frac{t}{\omega}$ gives
\begin{eqnarray*}
\omega^2\ddot y(t) = \ddot x\left(\frac{t}{\omega}\right) &=& k\sin(t)x\left(\frac{t}{\omega}\right) = k\sin(t)y(t) \\
\Rightarrow \qquad \ddot y &=& \frac{k}{\omega^2}\sin(t)y.
\end{eqnarray*}
This only has one parameter, $\frac{k}{\omega^2}$.  So, without loss of generality, the remainder of this note will study behaviour of the equation
\[ \ddot x=k\sin(t)x \]
for constant $k$.

\section{Linear Dynamics}
The equation of motion is second order and linear in $x$, so values at later times satisfy
\[ \vectwo{x(t)}{\dot x(t)} = \mxtwo{a(t)}{b(t)}{c(t)}{d(t)}\vectwo{x(0)}{\dot x(0)} \]
for some functions $a,b,c,d$.  Differentiating both sides gives
\begin{eqnarray*}
\vectwo{\dot x(t)}{k\sin(t)x(t)} =
 \vectwo{\dot x(t)}{\ddot x(t)} &=& \mxtwo{\dot a(t)}{\dot b(t)}{\dot c(t)}{\dot d(t)}\vectwo{x(0)}{\dot x(0)} \\
\Rightarrow \qquad  \mxtwo{c(t)}{d(t)}{k\sin(t)a(t)}{k\sin(t)b(t)}\vectwo{x(0)}{\dot x(0)} &=& \mxtwo{\dot a(t)}{\dot b(t)}{\dot c(t)}{\dot d(t)}\vectwo{x(0)}{\dot x(0)}.
\end{eqnarray*}
This gives the time derivatives of $a,b,c,d$:
\[ \dot a=c, \quad \dot b=d, \quad \dot c=k\sin(t)a, \quad \dot d=k\sin(t) b. \]
At $t=0$ the matrix is just the identity, so the initial conditions are 
\[ a(0)=1, \quad b(0)=0, \quad c(0)=0, \quad d(0)=1. \]

\section{Effective Focussing Strength}
The equation of motion is periodic in time, so solutions can be shifted by multiples of $2\pi$.  This means once the matrix
\[ M=\mxtwo{a(2\pi)}{b(2\pi)}{c(2\pi)}{d(2\pi)} \]
is obtained, it can evolve conditions at $t=2\pi n$ to $t=2\pi(n+1)$ for any integer $n$.  Thus, the long term dynamics from $t=0$ to $t=2\pi n$ are determined by $M^n$.

For stable dynamics, the two eigenvalues of $M$ must be on the complex unit circle $e^{\pm i \phi}$ for some phase advance $\phi$.  The trace of a matrix is equal to the sum of its eigenvalues, so
\newcommand{\tr}{\mathrm{tr}}
\[ \tr(M) = a(2\pi)+d(2\pi) = 2\cos(\phi). \]
The value of this trace determines the strength of the focussing, with $\tr(M)$=$2$ ($\phi=0$) meaning no focussing and $\tr(M)$=$-2$ ($\phi=\pi$) meaning the maximum focussing possible before instability.

This note will derive an expression for $\tr(M)$ and $\phi$ as functions of the time-rescaled $k$.  This gives the long-term dynamical frequencies as a function of applied focussing strength.

It is also possible to find an `effective focussing strength', which is the constant focusing strength (simple harmonic motion) that would give the same frequency:
\[ \ddot x = -k_\mathrm{eff} x 
\quad \Rightarrow \quad x=A\sin(\sqrt{k_\mathrm{eff}}t)+B\cos(\sqrt{k_\mathrm{eff}}t)
\quad \Rightarrow \quad \sqrt{k_\mathrm{eff}}2\pi = \phi. \]

\section{Expansion in Powers of $k$}
The time derivatives of $a,b,c,d$ seem to form a coupled system but expanding in powers of $k$ makes the calculation easier.  If $a=\sum_{n=0}^\infty a_n k^n$ and similarly for $b,c,d$, then equating powers of $k$ gives
\[ \dot a_n=c_n, \quad \dot b_n=d_n, \quad \dot c_n=\sin(t)a_{n-1}, \quad \dot d_n=\sin(t)b_{n-1}, \]
with $\dot c_0=\dot d_0=0$.  Writing the integrals explicitly gives, for $n\ge1$,
\[ a_n=\int_0^t c_n\de t, \qquad b_n=\int_0^t d_n\de t, \qquad
c_n=\int_0^t\sin(t)a_{n-1}\de t, \qquad d_n=\int_0^t\sin(t)b_{n-1}\de t.\]
The initial values imply that $c_0=0$ and $d_0=1$ for all $t$ because these have derivative zero.  There is also
\[ a_0=1+\int_0^t c_0\de t=1 \qquad \mathrm{and} \qquad b_0=\int_0^t d_0\de t=t. \]
Then $c_1=1-\cos(t)$ and $d_1=\sin(t)-t\cos(t)$, and so on, alternating between $a,b$ and $c,d$.

Successive approximations to the phase advance can be calculated from
\[  2\cos(\phi) = \tr(M) = a(2\pi)+d(2\pi) = \sum_{n=0}^\infty (a_n(2\pi)+d_n(2\pi)) k^n. \]

Note that 
\[ \mxtwo{a_0}{b_0}{c_0}{d_0} = \mxtwo{1}{t}{0}{1} \]
is just the transfer matrix of a drift with no focussing force, as expected for the $k^0$ term.

\section{Start of Calculation}
We already have the $k^0$ coefficient
\[ a_0(2\pi)+d_0(2\pi) = 1+1 = 2. \]

$a_1=t-\sin(t)$

$b_1=2-t\sin(t)-2\cos(t)$

The $k^1$ coefficient is
\[ a_1(2\pi)+d_1(2\pi) = 2\pi+(-2\pi) = 0. \]

$c_2=-\frac12t+\sin(t)+\frac14\sin(2t)-t\cos(t)$

$d_2=\frac{11}8-\frac14t^2+\frac14t\sin(2t)-2\cos(t)+\frac58\cos(2t)$

$a_2=2-\frac14t^2-t\sin(t)+\frac14\sin^2(t)-2\cos(t)$

The $k^2$ coefficient is
\[ a_2(2\pi)+d_2(2\pi) = (-\pi^2)+(-\pi^2) = -2\pi^2. \]

The integrals are probably best done by computer algebra system beyond this point.  Although, approximating $\tr(M) \simeq 2-2\pi^2k^2$ already gives a reasonable approximation to the real value as shown in Figure~\ref{fig:secondorder}.

\begin{figure}[!htb]
   \centering
   \includegraphics*{order2}
   \caption{Approximation of $\tr(M)$ up to the $k^2$ term.}
   \label{fig:secondorder}
\end{figure}

Odd powers of $k$ have coefficients of zero because replacing $k\sin(t)$ by $-k\sin(t)$ in the differential equation just shifts the time axis by $\pi$, so long term behaviour should be identical on changing the sign of $k$, making the phase advance an even function of $k$.

\section{General Terms}
It appears the general terms look like $p(t)\sin(mt)+q(t)\cos(mt)$ for some polynomials $p,q$ and $m\ge0$.  The $m=0$ $\cos$ term gives the non-sinusoidal polynomial parts.  The integral of this term also has the form $P(t)\sin(mt)+Q(t)\cos(mt)$ where the new polynomials $P,Q$ can be found by a recurrence relation.  Write $p(t)=\sum_{n=0}^Np_nt^n$ and similarly for $q,P,Q$, then
\begin{eqnarray*}
& \dd{}{t}(P(t)\sin(mt)+Q(t)\cos(mt)) = \sum_{n=0}^N \dd{}{t}(P_nt^n\sin(mt)+Q_nt^n\cos(mt)) \\
& = \sum_{n=0}^N (nP_nt^{n-1}-mQ_nt^n)\sin(mt)+(nQ_nt^{n-1}+mP_nt^n)\cos(mt) \\
& = \sum_{n=0}^N ((n+1)P_{n+1}-mQ_n)t^n\sin(mt)+((n+1)Q_{n+1}+mP_n)t^n\cos(mt) \\
& = \sum_{n=0}^N p_nt^n\sin(mt)+q_nt^n\cos(mt),
\end{eqnarray*}
where by convention, $P_{N+1}=Q_{N+1}=0$.  Equating coefficients gives 
\[ p_n = (n+1)P_{n+1}-mQ_n, \qquad q_n = (n+1)Q_{n+1}+mP_n, \]
therefore
\[ P_n = \frac1m(q_n-(n+1)Q_{n+1}), \qquad Q_n = \frac1m((n+1)P_{n+1}-p_n). \]
The recurrence can be started with $P_N = \frac{q_N}m$ and $Q_N = \frac{-p_N}m$ then evaluated downwards to $n=0$.  The exception is $m=0$ when the usual polynomial integration $P_n=\frac{p_{n-1}}{n}$ and $Q_n=\frac{q_{n-1}}{n}$ should be used. 

Besides integrating, the other operation that happens in the calculation is multiplication by $\sin(t)$.  This can be dealt with by the trigonometrical product formulae
\begin{eqnarray*}
\sin(t)\sin(mt) &=& \frac12\cos((m-1)t) - \frac12\cos((m+1)t) \\
\sin(t)\cos(mt) &=& -\frac12\sin((m-1)t) + \frac12\sin((m+1)t). 
\end{eqnarray*}

Overall, a matrix entry function is represented on a computer as
\[ a_n = \sum_{m=0}^M p_{a_nm}(t)\sin(mt)+q_{a_nm}(t)\cos(mt)
 = \sum_{m=0}^M \sum_{j=0}^N p_{a_nmj}t^j\sin(mt)+q_{a_nmj}t^j\cos(mt) \]
with arrays of coefficients $(p,q)_{(a,b,c,d)_nmj}$.  For $n=0$, the upper limits are $N=1$ and $M=0$ and these can increase by one each time $n$ does.

\section{Computer Algebra Calculation}
Using the general formulae above, a computer algebra calculation gives the trace as
\begin{eqnarray*}
& \tr(M) = 2\cos(\phi) = a(2\pi)+d(2\pi) = \\
&\left(+2\right)k^{0} \\
&+\left(-2\pi^{2}\right)k^{2} \\
&+\left(-\frac{25}{8}\pi^{2}+\frac{1}{3}\pi^{4}\right)k^{4} \\
&+\left(-\frac{1169}{144}\pi^{2}+\frac{25}{24}\pi^{4}-\frac{1}{45}\pi^{6}\right)k^{6} \\
&+\left(-\frac{16824665}{663552}\pi^{2}+\frac{24329}{6912}\pi^{4}-\frac{5}{48}\pi^{6}+\frac{1}{1260}\pi^{8}\right)k^{8} \\
&+\left(-\frac{104550461873}{1194393600}\pi^{2}+\frac{25241465}{1990656}\pi^{4}-\frac{14977}{34560}\pi^{6}+\frac{5}{1008}\pi^{8}-\frac{1}{56700}\pi^{10}\right)k^{10} \\
&+\left(-\frac{1383860829361699}{4299816960000}\pi^{2}+\frac{343096621171}{7166361600}\pi^{4}-\frac{7069153}{3981312}\pi^{6}+\frac{35579}{1451520}\pi^{8}-\frac{5}{36288}\pi^{10}+\frac{1}{3742200}\pi^{12}\right)k^{12} \\
&+\left(-\frac{29289023958538918009}{23702740992000000}\pi^{2}+\frac{150905495534989}{806215680000}\pi^{4}-\frac{131109204649}{17915904000}\pi^{6}+\right.\\&\left.\frac{9427513}{83607552}\pi^{8}-\frac{10301}{13063680}\pi^{10}+\frac{1}{399168}\pi^{12}-\frac{1}{340540200}\pi^{14}\right)k^{14} + O(k^{16}).
\end{eqnarray*}
The final $k^{14}$ term is of similar magnitude to the double precision rounding error, so this should suffice for most numerical calculations.

\section{Conclusion}
If the polynomial above is donated $T(k)\simeq \tr(M) = 2\cos(\phi)$, then the phase advance per period can be calculated with $\phi \simeq \arccos(\frac12T(k))$.

The effective focussing strength satisfies $\sqrt{k_\mathrm{eff}}2\pi = \phi$ and so
\[ k_\mathrm{eff} = \left(\frac{\phi}{2\pi}\right)^2 \simeq \left(\frac{\arccos(\frac12T(k))}{2\pi}\right)^2. \]

Finally, considering the non-time-rescaled equation of motion $\ddot x=k\sin(\omega t)x$, the subtitutions $k \leftarrow \frac{k}{\omega^2}$ and $2\pi \leftarrow 2\pi/\omega$ give
\[ k_\mathrm{eff} \simeq \omega^2\left(\frac{\arccos\left(\frac12T\left(\frac{k}{\omega^2}\right)\right)}{2\pi}\right)^2. \]

\end{document}