Dynamics in Electron Storage Ring

The energy loss and quantum fluctuation of synchrotron radiation will dominate the dynamics in a electron storage ring.

  • Energy loss --> Damping
  • Quantum fluctuation --> Excitation

The radiation power gives:

\begin{align} P=\frac{cC_\gamma }{2\pi} \frac{E^4}{\rho^2}=\frac{e^2 c^3 C_\gamma }{2\pi} E^2B^2 \end{align}

For longitudinal maps, we care about the total radiated energy:

\begin{align} U=\oint P dt = \frac{C_\gamma E^4 }{2\pi} \oint\frac{1}{\rho^2} ds \equiv\frac{C_\gamma E^4 }{2\pi} I_2 \end{align}

As you will see below, we will express results in a series of integrals through the ring, named radiation integral. Here is the definition of the second radiation integral.

And the total photon emits per turn is:

\begin{align} N=\frac{5\pi}{\sqrt{3}}\alpha\gamma \end{align}

Damped Longitudinal Motion

Since the radiation power is energy dependent, the particle with momentum deviation $\delta$ will radiate differently than the reference particles.

In longitudinal dynamics session, we have

\begin{equation} \Delta E_{n+1}= \Delta E_{n} + eV \left(\sin \phi_n-\sin \phi_s\right) \end{equation}

Adding the term of radiation energy loss, we should have:

\begin{equation} \Delta E_{n+1}= \Delta E_{n} + eV \left(\sin \phi_n-\sin \phi_s\right) - \frac{dU}{dE}\Delta E_{n} \end{equation}

Or in small amplitude/long time scale approximation, the differential equation becomes:

\begin{equation} \frac{d\Delta E}{ds}=\frac{eV\omega_0\cos\phi_s}{2\pi}\Delta\phi-\frac{\omega_0}{2\pi}\frac{dU}{dE}\Delta E \end{equation}

Together with the phase change equation:

\begin{equation} \frac{d\Delta\phi}{dt}=h\omega_0\eta\frac{dE}{E_0} \end{equation}

The we reach a second order differential equation of the beam energy deviation with damping terms.

\begin{equation} \frac{d^2\Delta E}{dt^2}+\frac{\omega_0}{2\pi}\frac{dU}{dE}\frac{\Delta E}{dt}+\frac{eV\omega_0^2h\left(-\eta\cos\phi_s\right)}{2\pi E_0}\Delta E=0 \end{equation}

If we define the synchrotron frequency $\omega_s$ and damping coefficient of energy $\alpha_E$ as:

\begin{align} \omega_s=\sqrt{-\frac{eVh\eta\cos\phi_s}{2\pi E_0}}\omega_0 \end{align}\begin{align} \alpha_E=\frac{1}{2T_0}\frac{dU}{dE} \end{align}

The equation becomes:

\begin{equation} \frac{d^2\Delta E}{dt^2}+2\alpha_E\frac{\Delta E}{dt}+\omega_s^2\Delta E=0 \end{equation}

And the solution is approximately:

\begin{equation} \Delta E = A \exp{\left(-\alpha_E t\right)} \cos\left(\omega_s t + \phi_0\right) \end{equation}

Where $A$ and $\phi_0$ is determined by the initial condition.

Now we need to calculate the damping term carefully, to taking to account

  • Energy deviation
  • Momentum compaction
\begin{align} U&=\oint Pdt\nonumber\\ &=\frac{1}{c}\oint P\left(1+\frac{x}{\rho}\right)ds \end{align}

The derivative with respect to $E$ gives:

\begin{align} \frac{dU}{dE} &=\frac{1}{c}\oint \left(\frac{dP}{dE}\left(1+\frac{x_\beta}{\rho}\right)+ P\frac{D(s)}{\rho E_0}\right) ds \end{align}

We know that $P \sim E^2 B^2$ therefore:

\begin{align} \frac{dP}{dE} &=2P\left(\frac{1}{E_0}+\frac{D}{B_0 E_0}\frac{dB}{dx}\right) \end{align}

After ignoring the second order term:

\begin{align} \frac{dU}{dE} &=\frac{1}{c}\oint \left[\frac{2P}{E_0}+\frac{2PD}{B_0E_0}\frac{dB}{dx}+\frac{PD}{E_0\rho}\right]\Bigg|_{E=E_0} ds\\ &=\frac{2U_0}{E_0}+\frac{1}{cE_0}\oint PD\left[\frac{2}{B_0}\frac{dB}{dx}+\frac{1}{\rho}\right]\Bigg|_{E=E_0} ds \end{align}

Here, every term should be evaluated for particle with energy $E_0$. Then plug in the form:

$$P(E_0)=\frac{cU_0}{I_2}\frac{1}{\rho^2}$$\begin{align} \frac{dU}{dE} &=\frac{2U_0}{E_0}+\frac{U_0}{E_0 I_2}\oint \frac{D(s)}{\rho}\left[{2}K_x(s)+\frac{1}{\rho^2}\right] ds \end{align}

Therefore the damping coefficient $\alpha_E$ becomes:

\begin{align} \alpha_E&=\frac{1}{2T_0}\frac{dU}{dE} \nonumber\\ &=\frac{U_0}{2T_0E_0}\left(2+\frac{1}{ I_2}\oint \frac{D(s)}{\rho}\left[{2}K_x(s)+\frac{1}{\rho^2}\right]\right) ds\\ &\equiv\frac{U_0}{2T_0E_0}\left(2+\mathcal{D}\right) \end{align}

Here we define damping partition number $\mathcal{D}$ and the fourth radiation integral $I_4$

\begin{align} \mathcal{D}=\frac{I_4}{I_2}=\frac{\oint \frac{D(s)}{\rho}\left[{2}K_x(s)+\frac{1}{\rho^2}\right] ds}{\oint ds/\rho^2 } \end{align}

When the accelerator uses separate function magnets with identical dipoles, viz $K=0$ in dipoles and $\rho$ are piecewise constant, the damping partition simply:

\begin{align} \mathcal{D}=\frac{I_4}{I_2}=\frac{1}{\rho}\oint \frac{D(s)}{\rho}ds =\frac{\alpha_c C}{2\pi \rho} \ll 1 \end{align}
In [5]:
import numpy as np
import longitudinal as lde
import matplotlib.pyplot as plt
%matplotlib notebook

turns=3000
energy=3e9  #3GeV
mass=0.511e6   # Electron

p0=np.sqrt(energy*energy-mass*mass)

frf=500e6
rho=25
harm=1300
voltage=5e6
mcf=3.68e-4
update_eta=True
# Gamma and Betas 
gamma=energy/mass
eta=mcf-1/gamma/gamma
beta=np.sqrt(1-1/gamma/gamma)
Usr=8.85e-5*np.power(energy/1e9,4)/rho*1e9
sin_phi_s=Usr/voltage

phi_s=np.arcsin(sin_phi_s)
yf=np.sqrt(np.cos(phi_s)-(np.pi-2*phi_s)*np.sin(phi_s)/2)
height=2*np.sqrt(voltage/2/np.pi/beta/beta/energy/harm/np.abs(eta))*yf

phi_0=np.pi-phi_s
#phi_0=phi_s
delta_E=height*energy*.99
npar=2
initial_phi=np.ones(npar)*(phi_0)
initial_de=np.linspace(0,1, npar)*delta_E

result=lde.longitudinal_evolve(turns, 
                                       initial_phi,
                                       initial_de,
                                       sin_phi_s=sin_phi_s, alphac=mcf, E0_ini=energy,
                                       mass=mass, e_volt=voltage, harm=harm,
                                       update_eta=update_eta, energy_change=False, damping_turn=1000
                                      )

fig,(ax1,ax2)=plt.subplots(1,2, figsize=(8,4))
ax1.set_title(r"Longitudinal Phase Space $(\phi,\delta)$")
ax2.set_title(r"Longitudinal Phase Space $(\phi,\Delta E)$")
ax1.set_xlabel('Phase [rad]')
ax2.set_xlabel('Phase [rad]')
ax1.set_ylabel(r'$\delta=dp/p_0$')
ax2.set_ylabel(r'$\Delta E [MeV]$')
#ax1.set_xlim([-np.pi,np.pi])
ax1.set_ylim([-2*height,2*height])
ax2.set_ylim([-2*height*energy*beta*beta/1e6,2*height*beta*beta*energy/1e6])
#ax1.set_prop_cycle(plt.cycler('color', plt.cm.Oranges(np.linspace(0, 1, 8))))
ax1.plot(result[0], result[2], linestyle=None)
#ax2.set_prop_cycle(plt.cycler('color', plt.cm.Blues(np.linspace(0, 1, 8))))
ax2.plot(result[0], result[1]/1e6, linestyle=None)

for ax in (ax1,ax2):
    ax.set_xticks([-np.pi, -0.5*np.pi, 0., .5*np.pi, np.pi, 1.5*np.pi, 2*np.pi])
    ax.set_xticklabels([r"$-\pi$", r"$-\frac{1}{2}\pi$", "$0$", r"$\frac{1}{2}\pi$",
                     r"$\pi$", r"$\frac{3}{2}\pi$", r"$2\pi$"])
    ax.set_xlim([0,2*np.pi])
plt.tight_layout()

Damped Transverse Motion

The idea of damping in transverse motion is that the recoil of SR is always parallel and opposite to the momentum $\mathbf{P}$. The lost energy will be compensated by the RF cavity, where the momentum gain is always along the desire orbit. Therefore the transverse emittance will keep shrinking.

Vertical motion

The vertical damping is simpler than the horizontal one because the dispersion function is everywhere zero.

Since the particle does not change the net energy after one turn, we will analyze this problem in $(x,x')$ phase space. The action of the particle is given by:

\begin{align} J_y = \beta y'^2 + 2\alpha y y' + \gamma y^2 \end{align}

In the process of SR, the momentum change is along the particle 's trajectory, therefore:

\begin{align} y_\text{after SR}&=y_\text{before}\\ y'_\text{after SR}&=y'_\text{before} \end{align}

The lost energy will be replenished by the RF cavity and the momentum change $\delta p$ will along the longitudinal direction:

\begin{align} y'_n=\frac{p_y}{p+\delta p}=y'\left(1-\frac{\delta p}{p}\right) \end{align}

The change of angle gives

\begin{align} \delta y'=-y'\frac{\delta p}{p}=-y'\frac{U}{E} \end{align}

The amplitude of the betatron motion gives:

\begin{align} A_y\sim\sqrt{\beta J} \end{align}

Therefore:

\begin{align} A_y\delta A_y=2\alpha\beta yy'+\beta^2 y'\delta y' \end{align}

After taking time average over the ring, and by averaging in one betatron oscillation $\left<yy'\right>=0$:

\begin{align} \frac{\delta A_y}{A_y}=-\frac{\beta^2 y'^2}{\beta J}\frac{U}{E} \end{align}\begin{align} \frac{1}{A_y}\frac{dA_y}{dt}=-\frac{U}{2ET_0} \end{align}

Therefore the damping coefficient of vertical direction becomes:

\begin{align} \alpha_y=\frac{U}{2ET_0} \end{align}

Horizontal motion

For horizontal direction, we still have

\begin{align} x_\text{after SR}&=x_\text{before}\\ x'_\text{after SR}&=x'_\text{before} \end{align}

The orbit has extra contribution from the dispersion:

\begin{align} x=x_\beta + x_d \end{align}

As we see from the vertical plane, both the position and angle does not change therefore:

\begin{align} \delta x_\beta &=-\delta x_d= D \frac{u}{E} \\ \delta x'_\beta &=-\delta x'_d= D' \frac{u}{E} \end{align}

where $u$ is the SR energy loss in short distance $dl$

Therefore the amplitude change becomes:

\begin{align} A_x\delta A_x=\beta\gamma x_\beta \delta x_\beta+ \beta\alpha\left(x_\beta\delta x'_\beta+x'_\beta \delta x_\beta\right) + \beta^2 x'_\beta \delta x'_\beta \end{align}

Taking the average over on betatron oscillation $\left<\alpha\right>=0$ and $\left<x'_\beta\right>=0$:

\begin{align} A_x\delta A_x&=\beta\gamma x_\beta \delta x_\beta+ \beta\alpha\left(x_\beta\delta x'_\beta+x'_\beta \delta x_\beta\right) + \beta^2 x'_\beta \delta x'_\beta \\ &=x_\beta \delta x_\beta \\ &=x_\beta D \frac{u}{E} \end{align}

$u$ is calculated as:

\begin{align} u&=\frac{P\left(x_\beta\right)}{c}dl \\ &=\left(P\left(0\right)+2 x_\beta \frac{P(0)}{B}\frac{dB}{dx} \right)\left(1+\frac{x_\beta}{\rho}\right)\frac{ds}{c}\\ &=\frac{P}{c}\left(1+x_\beta \frac{2}{B}\frac{dB}{dx} +\frac{x_\beta}{\rho}\right)ds \end{align}

Then the amplitude change over the ring with averaging $\left<x_\beta\right>=0$ and $\left<x_\beta^2\right>=A_x^2/2$:

\begin{align} \frac{\delta A_x}{A_x}&= \frac{1}{2cE}\oint DP\left(\frac{2}{B}\frac{dB}{dx} +\frac{1}{\rho}\right)ds \\ &=\frac{U_0}{2E}\frac{\oint \frac{D}{\rho}\left(2K(s) +\frac{1}{\rho^2}\right)ds}{\oint ds/\rho^2} \\ &=\frac{I_4}{I_2}\frac{U_0}{2E} =\mathcal{D}\frac{U_0}{2E} \end{align}

The above result is before the considering the RF acceleration, which is same for both transverse direction. Combine these two factors (RF and dispersion):

\begin{align} \frac{\delta A_x}{A_x}=-(1-\mathcal{D})\frac{U_0}{2E} \end{align}

Or the horizontal damping coefficient reads

\begin{align} \alpha_x=(1-\mathcal{D})\frac{U_0}{2ET_0} \end{align}

Robinson Theorem

If we define:

\begin{align} \alpha_0=\frac{U_0}{2ET_0} \end{align}

Then the damping coefficients are:

\begin{align} \alpha_x&=(1-\mathcal{D})\alpha_0 \equiv \mathcal{J}_x\alpha_0\\ \alpha_y&=\alpha_0 \equiv \mathcal{J}_y\alpha_0\\ \alpha_E&=(2+\mathcal{D})\alpha_0 \equiv \mathcal{J}_E\alpha_0 \end{align}

And the Robinson theorem states:

\begin{align} \mathcal{J}_x+\mathcal{J}_y+\mathcal{J}_E=4 \end{align}

And the damping time constant:

\begin{align} \tau_{x/y/E}=\frac{2E}{\mathcal{J}_{x/y/E}U_0} T_0 \end{align}

And the damping decrements:

\begin{align} \lambda_{x/y/E}=\frac{\mathcal{J}_{x/y/E}U_0}{2E} \end{align}

The damping effect can be 'redistributed' by adjusting the damping partition $\mathcal{D}$. Besides the lattice design, the insertion device also can change the damping partition.

Quantum Excitation

The quantum effect of the SR satisfy the following properties:

  • Discrete energy
  • Instantaneous
  • Random emission, satisfy Poison distribution

Equilibrium energy spread

The longitudinal oscillation, without radiation damping, gives:

\begin{align} \Delta E = A \cos\left(\omega_s t + \phi_0\right) \end{align}

We are interested in the time average of the amplitude square of the longitudinal oscillation, $A^2$. The damping effect for $A^2$ has half of the damping time:

\begin{align} \frac{d\left<A^2\right>}{dt}=-2\frac{\left<A^2\right>}{\tau_E} \end{align}

On the other hand, the electron emits discrete photons, which is random. Assuming that, at time $t_1$, the electron radiates a photon with energy u:

\begin{align} \Delta E = A \cos\left(\omega_s t + \phi_0\right) - u \cos\left(\omega_s t -\omega_s t_1\right) \end{align}

Since $t_1$ is random, in time average, we have the amplitude change in $\Delta t$:

\begin{align} \Delta \left<A^2\right>=\left<\int n(u) u^2 du \right>_{\Delta t} \Delta t \end{align}

where $n(u)$ is the number of photon of energy $u$.

Therefore, we can combine the excitation and damping:

\begin{align} \frac{d\left<A^2\right>}{dt}=-2\frac{\left<A^2\right>}{\tau_E} +\int n(u) u^2 du \end{align}

And the equilibrium energy spread is given by $\frac{d\left<A^2\right>}{dt}=0$ and averaged over one turn

\begin{align} \sigma_E^2&=\frac{\left<A^2\right>}{2}=\frac{1}{4}\tau_E \left<\int n(u) u^2 du\right> \nonumber\\ &=\frac{1}{4}\tau_E \left<\frac{15\sqrt{3}}{8}\frac{U_{SR}}{u_c} \frac{11}{27}u_c^2\right> \nonumber\\ &=\frac{1}{4}\frac{55\sqrt{3}}{72}\frac{2E}{\mathcal{J}_{E}U_0} T_0 \left<\frac{C_\gamma E^4}{2\pi} \frac{3\hbar\gamma^3c}{2\rho^3}\right> \nonumber\\ &=\frac{1}{4}\frac{55\sqrt{3}}{72}\frac{2E}{\mathcal{J}_{E}I_2} \left< \frac{3\hbar\gamma^3c}{2\rho^3}\right>\nonumber\\ &=\frac{3}{4}\frac{55\sqrt{3}}{72}\frac{c\hbar E\gamma^3}{\mathcal{J}_{E}} \frac{\oint\rho^{-3} ds }{I_2} \end{align}

We define the third radiation integral $I_3=\oint\rho^{-3} ds$

Therefore the relative energy spread gives:

\begin{align} \frac{\sigma_E^2}{E^2} &=\frac{55\sqrt{3}}{96}\frac{\hbar}{mc}\frac{\gamma^2}{\mathcal{J}_{E}} \frac{ I_3}{I_2} \\ &\equiv C_q\frac{\gamma^2}{\mathcal{J}_{E}} \frac{ I_3}{I_2} \end{align}

For electron, $C_q$ is a constant, which is $3.83\times10^{-10}$m. Surprisingly, the energy spread does not depend on RF parameters. But the bunch length does. From the longitudinal dynamics section, we know that:

\begin{equation} \frac{\left<\delta\right>}{\left<\Delta\phi\right>}=\frac{Q_s}{h\left|\eta\right|} \end{equation}

It converts to:

\begin{equation} \sigma_t=\frac{\alpha_c}{\omega_s}\frac{\sigma_E}{E} \end{equation}

Equilibrium transverse emittance

First, we look at the horizontal direction, we know that the action $J$ of the horizontal direction is damped by SR as:

\begin{align} \frac{dJ}{dt}=-2\frac{J}{\tau_x} \end{align}

But the instantaneous emission of a energy quanta $u$ gives an prompt kick of

\begin{align} \delta x_\beta&=-D\frac{u}{E} \\ \delta x'_\beta&=-D'\frac{u}{E} \end{align}

The change of action gives:

\begin{align} \Delta J=\beta\left(x'_\beta+\delta x'_\beta\right)^2+2\alpha\left(x_\beta+\delta x_\beta\right)\left(x'_\beta+\delta x'_\beta\right)+\gamma\left(x_\beta+\delta x_\beta\right)^2 -\left(\beta {x'}_\beta^2 +2\alpha x_\beta x'_\beta+\gamma x_\beta^2 \right) \end{align}

Taking average of one turn, we have

\begin{align} \left<\Delta J\right>&=\left<\left(\beta D'^2 + 2\alpha D D' + \gamma D^2\right)\left(\frac{u}{E}\right)^2 \right> \\ &=\left<\frac{\mathcal{H}}{E^2}\left({\int n(u)u^2du}\right)^2 \right> \end{align}

Therefore the equilibrium betatron emittance gives:

\begin{align} \epsilon_x&=\frac{\left<J\right>}{2} \\ &=\frac{1}{4}\tau_x\left<\frac{\mathcal{H}}{E^2}\left({\int n(u)u^2du}\right)^2 \right> \\ &=\frac{1}{4}\frac{55\sqrt{3}}{72}\frac{2}{\mathcal{J}_{x}EI_2} \left< \mathcal{H}\frac{3\hbar\gamma^3c}{2\rho^3}\right>\\ &=\frac{55\sqrt{3}}{96}\frac{\hbar}{mc}\frac{\gamma^2}{\mathcal{J}_{x}I_2} \oint \frac{\mathcal{H}}{\rho^3}ds\\ &= C_q\frac{\gamma^2}{\mathcal{J}_{x}I_2} \left(\oint \frac{\mathcal{H}}{\rho^3}ds \right)\\ &\equiv C_q\frac{\gamma^2 I_5}{\mathcal{J}_{x}I_2} \end{align}

Here we define the fifth radiation integral $I_5=\oint\mathcal{H}/\rho^{-3} ds$

The transverse beam size, therefore, has contribution from both the equilibrium betatron motion and the energy spread/dispersion:

\begin{align} \sigma_x ^2 &= \beta_x \epsilon_x + D_x^2 \frac{\sigma_E^2}{E^2} \\ &=\frac{C_q\gamma^2}{I_2}\left(\beta_x \frac{I_5}{\mathcal{J}_x}+D_x^2\frac{I_3}{\mathcal{J}_E}\right) \end{align}

For vertical direction, follow the above treatment will give the equilibrium emittance is zero! Therefore, we have to consider higher order effect: the SR cone angle is about $\theta \sim 1/\gamma$ the random vertical betatron angle change is:

\begin{align} \delta y'\sim \frac{1}{\gamma}\frac{u}{E} \end{align}

Follow the similar treatment: we can easily find:

\begin{align} \epsilon_y = C_q\frac{ I_3}{\mathcal{J}_{y}I_2} = C_q\frac{I_3}{I_2} \end{align}

which is an extremely small number.

However, in reality, we can use betatron coupling to control the emittance of both transverse direction.

We found the relations of radiation integral and the storage ring quantities in the table:

List of Radiation Integrals Quantities in Electron Storage Ring
$I_1=\oint\frac{D}{\rho}ds$ $\alpha_c=I_1/ C$
$I_2=\oint\frac{1}{\rho^2}ds$ $U_{SR}=\frac{C_\gamma}{2\pi} E^4 I_2 $
$I_3=\oint\frac{1}{\rho^3}ds$ $\left(\frac{\sigma_E}{E}\right)=C_q\gamma^2\frac{I_3}{2I_2+I_4}$
$I_4=\oint\frac{D}{\rho}\left(\frac{1}{\rho^2}+2K\right)ds$ $\mathcal{D}=\frac{I_4}{I_2}$
$I_5=\oint\frac{\mathcal{H}}{\rho^3}ds$ $\epsilon_x=C_q\gamma^2\frac{I_5}{I_2-I_4}$

Transverse emittance optimization

A general practice is to minimized the emittance of the electron storage ring, i.e. minimize the $I_5/I_2$ since in separate magnet machine $I_4 \ll 1$. If we look at the form of $I_5$, we need to minimize $\mathcal{H}$ in dipole. Let's consider a general case here to minimize $I_5$ without any other constrains.

Consider the dipole with length $L$ and bending angle $\theta\ll 1$. The optical function on both side of the dipole can be arbitrary. Let's assume the optics at entrance of the dipole reads: $\beta=\beta_0$, $\alpha=\alpha_0$ and $\gamma_0=\left(1+\alpha_0^2\right)/\beta_0$. The dispersion function at the same location is $D_0$ and $D'_0$.

At location $s$, the twiss parameters read:

\begin{align} \beta(s)&=\frac{1}{\gamma_0}+\gamma_0\left(s-\frac{\alpha_0}{\gamma_0}\right)^2\\ \alpha(s)&=\alpha_0-\gamma_0 s \\ \gamma(s)&=\gamma_0 \end{align}

From the dipole's $3\times 3$ matrix $\mathcal{M}$

\begin{align} \mathcal{M}(s)= \left(\begin{array}{cc} 1 & s & \frac{s^2}{2\rho} \\ 0 & 1 & \frac{s}{\rho} \\ 0 & 0 & 1 \end{array}\right) \end{align}

At arbitrary location, the dispersion function gives:

\begin{align} D(s)&=\frac{s^2}{2\rho}+D_0+D'_0 s\\ D'(s)&=\frac{s}{\rho}+D'_0 \end{align}
In [2]:
import sympy
sympy.init_printing(use_unicode=True)

l,l1,t,f=sympy.symbols("L L_1 \\theta f")



beta0, tb, alpha0= sympy.symbols("\\beta_0  \\tilde\\beta \\alpha_0")
d0, td, dp0, tdp=sympy.symbols("D_0 \\tilde{d} D_p \\tilde{d_p}")

rho, s= sympy.symbols("\\rho s")
gamma0=(1+alpha0*alpha0)/beta0

beta=1/gamma0+gamma0*(s-alpha0/gamma0)*(s-alpha0/gamma0)
alpha=alpha0-gamma0*s

d=s*s/2/rho+d0+dp0*s
dp=s/rho+dp0

H=sympy.simplify(beta*dp*dp+2*alpha*d*dp+gamma0*d*d)
H0=H.subs(s,0)
Hint=sympy.integrate(H,s)
H_ave=(Hint.subs(s,l)-Hint.subs(s,0))/l
H_simp=sympy.simplify(H_ave.subs(beta0, tb*l).subs(d0,td*l*t).subs(dp0,t*tdp).subs(rho,l/t)/l/l/l*rho*rho)
#H_simplified.subs(alpha0, 10*tbeta/4)
#H_simplified.subs(alpha0, sympy.sqrt(15)).subs(tbeta,sympy.sqrt(12)/sympy.sqrt(5))

Hoptd=sympy.simplify(H_simp.subs(td, sympy.Rational(1,6)).subs(tdp, sympy.Rational(-1,2)))
sympy.simplify(Hoptd.diff(alpha0))

Hoptd.diff(tb).simplify()
#Hoptd.diff(alpha0).simplify()

Hoptd.subs(tb, 8/sympy.sqrt(15)).subs(alpha0, sympy.sqrt(15))
Out[2]:
$$\frac{\sqrt{15} \rho^{2} \theta^{2}}{180 L^{2}}$$

By defining ${\tilde{\beta}}=\beta_0/L$, ${\tilde{\gamma}}=\gamma_0 L$, ${\tilde{d}}=D_0/(L\theta)$ and ${\tilde{d}_p}=D'_0/\theta$, we found

\begin{align} \rho\frac{I_5}{I_2} &= \left<\mathcal{H}\right> \\ &=\rho\theta^3\left( \tilde{\beta}\tilde{d}_p^2+2\alpha_0\tilde{d}\tilde{d}_p+\tilde{\gamma}\tilde{d}^2 +\left(\alpha_0-\frac{\tilde{\gamma}}{3}\right)\tilde{d}\\ +\left(\tilde{\beta}-\frac{\alpha_0}{3}\right)\tilde{d}_p +\frac{\tilde{\beta}}{3}-\frac{\alpha}{4}+\frac{\tilde{\gamma}}{20} \right) \end{align}

First we can solve for the optimum dispersion function:

\begin{align} 2 \left( \begin{array}{cc} \alpha_0 & \tilde{\beta} \\ \tilde{\gamma} & \alpha_0 \\ \end{array} \right) \left( \begin{array}{cc} \tilde{d} \\ \tilde{d}_p \\ \end{array} \right) &=- \left( \begin{array}{cc} \tilde{\beta}-\frac{\alpha_0}{3} \\ \alpha_0-\frac{\tilde{\gamma}}{3} \\ \end{array} \right) \end{align}

Therefore we have $\tilde{d}=1/6$,$\tilde{d_p}=-1/2$. And it does not depend on the optics function! This actually relates to the symmetric condition. Then the ratio of the radiation integral reduces to :

\begin{align} \rho\frac{I_5}{I_2} &=\frac{\rho\theta^3}{12}\left( \tilde{\beta}-\alpha_0+\frac{4\tilde{\gamma}}{15} \right) \end{align}

After taking derivative, we have:

\begin{align} 15\tilde{\beta}&=8\alpha_0\\ 15\tilde{\beta}&=4\tilde{\gamma} \end{align}

Then we have

\begin{align} \beta&=\frac{8}{\sqrt{15}} L \\ \alpha&=\sqrt{15} \\ s^*&=\frac{\alpha_0}{\tilde\gamma}L=\frac{L}{2}\\ \beta^*&=\frac{1}{2\sqrt{15}} L \end{align}

and finally the optimized

\begin{align} \rho\frac{I_5}{I_2} = \left<H\right>_{min} &=\frac{\rho\theta^3}{12\sqrt{15}} \end{align}
In [3]:
fig, ax=plt.subplots()
ax_m=ax.twinx()
norm_s=np.linspace(0,1,1000)
d=norm_s*norm_s/2.0+1.0/6-norm_s/2.0
ds=norm_s-0.5
be_s=1/2/np.sqrt(15.0)
s_star=0.5
be=be_s+(norm_s-s_star)*(norm_s-s_star)/be_s
al=-(norm_s-s_star)/be_s
h=be*ds*ds+2*al*d*ds+d*d/be_s
print(np.average(h))
l1=ax.plot(norm_s,be, label=r'$\beta$ function')
ax.set_xlabel(r's [in $L$]')
ax.set_ylabel(r'$\beta$ [in $L$]')
ax_m.set_ylabel(r'$D$ [in $L\theta$], $\mathcal{H}$ [in $L\theta^2$]')

l2=ax_m.plot(norm_s,d, color='g',label=r'Dispersion function')
l3=ax_m.plot(norm_s,h, color='r',label=r'$\mathcal{H}$ function')
lns=l1+l2+l3
labs = [l.get_label() for l in lns]
ax.legend(lns,labs,loc='best')

0.021581253096642352
Out[3]:
<matplotlib.legend.Legend at 0x11331a630>
In [ ]: