do directly. However, we may simplify matters by noting that

$\displaystyle \frac{(p_y - a_y)^2}{2(1+p_z)} =  \exp(:-\displaystyle\int_{0}^{y} a_y dy':)\frac{p_x^2}{2(1+p_z)}$

which may be proven by simply carrying out the definition of the Lie exponential. Next, we apply the well-known Lie algebraic identity

$\exp(:A:)\exp(:B:)\exp(:-A:) = \exp(:\exp(:A:)B:)$

and find that

$$\begin{gathered}
I_{y/2} =  \exp(:-\frac{\Delta s}{2} H_y:) \\
=  \exp(:-\displaystyle\int_{0}^{y} a_y dy':)\exp(:-\frac{\Delta s}{2}\frac{p_x^2}{2(1+p_z)}:)\exp(:\displaystyle\int_{0}^{y} a_y dy':).\end{gathered}$$

The argument of the middle exponential is a simple Hamiltonian, and thus may be evaluated by analytically solving Hamilton’s equations. The outer exponentials may be evaluated by applying the definition of the Lie exponential operator. Doing so gives

$\exp(:\pm\displaystyle\int_{0}^{y} a_y dy':)  :
\begin{cases}
p_y \rightarrow p_y \pm a_y       \\
p_x \rightarrow p_x \pm \frac{\partial}{\partial x}\displaystyle\int_{0}^{y} a_y dy
\end{cases}$

with all other variables being unchanged. Furthermore, evaluating $I_{z}$ by applying the definition of the Lie exponential gives

$\exp(:\Delta s a_z:)  :
\begin{cases}
p_y \rightarrow p_y  + \Delta s \frac{\partial a_z}{\partial y}       \\
p_x \rightarrow p_x + \Delta s \frac{\partial a_z}{\partial x}
\end{cases}$

with all other variables being unchanged.

For our model, all possible terms are already accounted for in the Bmad symplectic integrator except for the constant term. Thus, in order to make the appropriate changes to the program, we calculated 