Ion-photon interaction kinematics and laser cooling of ion in a storage ring¶
A. Petrenko (CERN, 2017)
This notebook describes the kinematics of photon absorption and emission by a relativistic partially stripped ion and application of this process to ion beam cooling in a storage ring. It can be executed as psi_cooling.ipynb at the MS Azure cloud service.
There is also a more detailed simulation of Li-like Pb cooling in the SPS.
In the ion's initial frame of reference:
Let's use 4-vectors $(E/c, \boldsymbol{p})$ for the Lorentz transformations. Ion's 4-momentum is $(\gamma m c, \gamma m \boldsymbol{v})$, and the photon's 4-momentum is $(\hbar \omega / c, \hbar \boldsymbol{k})$. If the $Z$-axis is aligned with the direction of ion motion, then Lorentz transformation can be written as
$$ \begin{pmatrix}E'/c \\ p_x'\\ p_y'\\ p_z' \end{pmatrix} = \begin{pmatrix} \gamma & 0 & 0 & -\beta\gamma \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ -\beta\gamma & 0 & 0 & \gamma \\ \end{pmatrix} \begin{pmatrix}E/c \\ p_x \\ p_y \\ p_z \end{pmatrix}. $$
Assuming that $k_x = -k\sin\theta$, $k_y = 0$, $k_z = -k\cos\theta$, and $k = \omega/c$ we can find the incoming photon parameters in the ion's frame of reference:
$$ \left( \begin{array}{c} 1 \\ -\sin\theta' \\ 0 \\ -\cos\theta' \end{array} \right) \frac{\omega'}{c} = \begin{pmatrix} \gamma & 0 & 0 & -\beta\gamma \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ -\beta\gamma & 0 & 0 & \gamma \\ \end{pmatrix} \left( \begin{array}{c} 1 \\ -\sin\theta \\ 0 \\ -\cos\theta \end{array} \right) \frac{\omega}{c}. $$
Since $k = \omega/c$,
$$ \omega' = ( 1 + \beta \cos\theta )\gamma \omega \approx \left( 1 + \beta - \beta\frac{\theta^2}{2} \right) \gamma \omega \approx 2 \gamma \omega.\\ $$
Incoming angular spread in the beam of $\theta \sim 1$ mrad will be translated to a frequency error of only $\sim10^{-6}$ in the ion's frame of reference. Frequency mismatch is dominated by the energy spread in the ion beam (typically $\sim 10^{-4}$). It also means that the laser beam can cross the ion trajectory with an angle of several mrad at least.
$$ \omega' \sin\theta' = \omega \sin\theta, $$
therefore $\theta'$ is very small
$$ \theta' \approx \frac{\theta}{2\gamma}. $$
Excited ion after the photon absorption:
And here is the ion after the photon emission:
Photon emission will occur in a random direction. For simplicity let's assume that the photon was emitted in the same plane ($X',Z'$) at a random angle $\theta'_1$, i.e. $k'_{1x} = k'\sin\theta'_1$, $k'_{1z} = k'\cos\theta'_1$. Then inverse Lorentz transformation gives us the emitted photon parameters in the lab frame:
$$ \left( \begin{array}{c} 1 \\ ~\sin\theta_1 \\ 0 \\ \cos\theta_1 \end{array} \right) \frac{\omega_1}{c} = \begin{pmatrix} \gamma & 0 & 0 & \beta\gamma \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ \beta\gamma & 0 & 0 & \gamma \\ \end{pmatrix} \left( \begin{array}{c} 1 \\ ~\sin\theta'_1 \\ 0 \\ \cos\theta'_1 \end{array} \right) \frac{\omega'}{c}. $$
Hence the scatterd photon has the frequency $\omega_1 = \gamma(1 +\beta\cos\theta'_1)\omega' \approx 2\gamma^2(1 +\beta\cos\theta'_1) \omega$.
$$ \omega_1 \sin\theta_1 = \omega' \sin\theta'_1 ~\Rightarrow~ \sin\theta_1 = \frac{\sin\theta'_1} {\gamma(1 +\beta\cos\theta'_1)}. $$
$$ \theta_1 \sim \frac{1}{\gamma}. $$
Let's consider some examples¶
from math import *
A = 208 # Atomic mass
mc = A*0.931e9 # eV/c
mc2 = mc # eV
Z = 82 # Number of protons in the ion (Lead)
Ne = 1 # Number of remaining electrons
gamma_p = 7000/0.938 # max. gamma for the protons in the LHC
gamma = gamma_p*(Z-Ne)*0.938e9/mc # Lorentz factor of the ion
gamma_s = gamma # defined for later use below
# Let's select the transition between the 1st and the 2nd level (Lyman-alpha line)
n1 = 1; n2 = 2
# According to the Rydberg formula (https://en.wikipedia.org/wiki/Rydberg_formula)
hw0 = (1/(n1*n1) - 1/(n2*n2))*Z*Z*13.6 # eV
# Corresponding photon energy in the lab frame will be
hw = hw0/(2*gamma)
# max energy of scattered photons will be
hw1 = 4*gamma*gamma*hw
from IPython.display import Latex
Latex(r"""
\begin{aligned}
\gamma &= %.0f \\
\hbar \omega' &= %.0f~\rm{keV} \\
\hbar \omega~ &= %.0f~\rm{eV} \\
\hbar \omega_{1, \rm{max}} &= %.0f~\rm{MeV} \\
\end{aligned}
""" % (gamma, hw0/1000, hw, hw1/1e6))
Typical transverse kick obtained by an ion due to the photon scattering is very small compared to the typical angular spread in the beam:
px = hw1/gamma # eV/c
pz = gamma*mc # eV/c
p_s = pz # defined for later use below
kick = (px/pz)*1000 # mrad
print('(px ~ %.1f MeV, pz = %.0f TeV) Typical transverse kick ~ %.0e mrad' % (px/1e6, pz/1e12, kick))
# Angular spread in the beam at the LHC interaction point for example:
emitt_n = 2.0 # mm*mrad (normalized emittance)
emitt = emitt_n/gamma
beta_x = 1.0 # m
sigma_x = sqrt(beta_x*emitt) # mm
sigma_xp = emitt/sigma_x
print('Transverse beam size = %.3f mm, angular spread = %.3f mrad' % (sigma_x, sigma_xp))
Unperturbed betatron oscillations
Ion phase-space coordinates $\boldsymbol{X} = (x, x', y, y', \Delta l, \Delta p/p)^T$, where $x' = p_x/p_z$, and $y' = p_y/p_z$. These $p_x, p_y, p_z$ are slightly different from the ones used above because now the direction of $Z$-axis is not perfectly aligned with the direction of ion motion. Now we use the local coordinate system standard for accelerator physics.
Here is the description of ion trajectory without the scattering and assuming no coupling between horizontal, vertical and longitudinal betatron oscillations:
$$ \left( \begin{array}{c} x~ \\ x' \end{array} \right)_{n+1} = \begin{pmatrix} \cos\mu_x + \alpha_x\sin\mu_x & \beta_x \sin\mu_x \\ -\gamma_x \sin\mu_x & \cos\mu_x - \alpha_x \sin\mu_x \\ \end{pmatrix} \left( \begin{array}{c} x~ \\ x' \end{array} \right)_{n}, $$
where $n$ is the turn number, $\mu_x = 2\pi\nu_x$ ($\nu_x$ -- betatron tune), $\alpha_x = -\beta'_x / 2$, $\gamma_x = (1 + \alpha_x^2)/\beta_x$.
Let's assume zero alpha-function at this location for simplicity, then
$$ \left( \begin{array}{c} x~ \\ x' \end{array} \right)_{n+1} = \begin{pmatrix} \cos\mu_x & \beta_x \sin\mu_x \\ -\beta^{-1}_x \sin\mu_x & \cos\mu_x \\ \end{pmatrix} \left( \begin{array}{c} x~ \\ x' \end{array} \right)_{n}. $$
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
nux = 0.31
mux = 2*pi*nux
N_turns = 100
x = np.zeros(N_turns + 1) # m
xp = np.zeros(N_turns + 1) # rad
x[0] = sigma_x/1e3
for n in range(0, N_turns):
x[n+1] = x[n]*cos(mux) + xp[n]*beta_x*sin(mux)
xp[n+1] = -x[n]*sin(mux)/beta_x + xp[n]*cos(mux)
plt.scatter(x*1e3, xp*1e3)
plt.xlabel("x (mm)")
plt.ylabel("x' (mrad)")
plt.grid()
plt.show()
Unperturbed synchrotron oscillations
Unperturbed longitudinal motion is described in terms of energy gain of the ion when it passes the RF-resonator as well as the ion phase with respect to the RF:
$$ E_{n+1} = E_n + e(Z-N_e)V_{\rm{RF}}\cos\phi_n, $$
where $E_n$ is the ion energy and $\phi_n$ is its phase in RF resonator at the $n$-th turn; RF system freqyency $f_{\rm{RF}}$ is some high harmonic $h$ of ion revolution frequency: $$ T_s = \frac{L}{\upsilon_s} = \frac{h}{f_{\rm{RF}}}, $$ where $L$ is the ring perimeter, and $\upsilon_s$ is the reference velocity of the ion, and $T_s$ is the revolution period. If the ion momentum is different from its reference value then the period of revolution $T$ is different from $T_s$: $$ T = \frac{L+\Delta l}{\upsilon_s + \Delta \upsilon} \approx T_s \left ( 1 + \frac{\Delta l}{L} - \frac{\Delta \upsilon}{\upsilon_s} \right ) \approx T_s \left ( 1 + \frac{\Delta l}{L} - \frac{1}{\gamma^2} \frac{\Delta p}{p} \right ). $$ The difference between the length of ion trajectory and the reference orbit length is given by the $M_{56}$ element of the 1-turn transport matrix: $$ \Delta l = M_{56} \frac{\Delta p} {p}. $$ $\Delta T = T - T_s$ in can be written as $$ \Delta T \approx T_s \left ( \frac{M_{56}}{L} - \frac{1}{\gamma^2} \right ) \frac{\Delta p}{p} = T_s \left ( \frac{1}{\gamma_t^2} - \frac{1}{\gamma^2} \right ) \frac{\Delta p}{p} = T_s\eta\frac{\Delta p}{p}. $$ Then the phase is $$ \phi_{n+1} = \phi_n + 2\pi f_{\rm{RF}}(T_s + \Delta T_{n+1}) = \phi_n + 2\pi h + 2\pi h \eta\frac{\Delta p_{n+1}}{p}, $$ $\Delta p$ and $\Delta E$ are related as $$ \frac{\Delta p}{p} = \frac{1}{\beta^2}\frac{\Delta E}{E}, $$ where $\beta = \upsilon/c$.
N_turns = 1000
dE = np.zeros(N_turns + 1) # energy deviation (eV)
phi = np.zeros(N_turns + 1) # ion phase with respect to RF
eVrf = 16e6 # eV
E_s = gamma_s*mc2 # eV
beta_s = sqrt(1-1/(gamma_s*gamma_s))
dE[0] = 3.5e-4*E_s # eV
phi[0] = pi/2
gamma_t = 55.68 # LHC
eta = 1/(gamma_t*gamma_t) - 1/(gamma_s*gamma_s)
h = 35640 # LHC
for n in range(0, N_turns):
dE[n+1] = dE[n] + eVrf*(Z-Ne)*cos(phi[n])
phi[n+1] = phi[n] + (2*pi*h*eta/(beta_s*beta_s))*dE[n+1]/E_s
#plt.scatter(range(0, N_turns+1),100*dE/E_s)
#plt.xlabel("Turn")
plt.scatter(phi/pi, 100*dE/E_s, s=2, linewidths=0)
plt.xlim(-0.5, +1.5)
plt.ylim(-0.04, +0.04)
plt.xlabel("Phase/pi")
plt.ylabel("dE/E (%)")
plt.grid()
plt.show()
The perturpation due to the photon scattering
The effect of photon absorption on the ion momentum in the lab frame can be neglected. We should take into account only the emitted photon. On every turn we will select a random direction of photon emission in the ion's frame of reference. The random polar angle $\theta_1'$ can be obtained as arccos of uniformly distributed random number from $-1$ to $+1$. Then using the above expressions we can find the frequency of the scattered photon as well as $\theta_1$ angle between the initial ion momentum and the direction of photon emission in the lab frame. Finally a random azimuthal angle (from $0$ to $2\pi$) is needed to get the direction of photon scattering in the plane which is perpendicular to the initial ion momentum.
$\theta_1'$ is named theta0 in the script below:
Longitudinal effects
The effect of photon scattering is easy to simulate in the longitudinal plane. The previous script should be modified only slightly:
N_turns = 20000
dE = np.zeros(N_turns + 1) # energy deviation (eV)
phi = np.zeros(N_turns + 1) # ion phase with respect to RF
dE[0] = 1.0e-4*E_s # eV # typical value for the LHC
phi[0] = pi/2
# number of ion-photon scatterings per turn:
N_sc = 4
# for simplicity we consider here a fixed number of photons scattered
# by the ion per one turn. In a more detailed Monte Carlo simulation
# the scattering event should have some probability proportional to
# the laser beam intensity.
eVrf = 16e6 # eV # voltage can be increased to keep the ion in the bucket.
# Other constants are already defined above.
for n in range(0, N_turns):
# photon scattering:
for sc in range(0, N_sc):
costheta = np.random.uniform(-1,+1)
# random polar angle in the ion's frame:
theta0 = acos(costheta)
gamma = (E_s + dE[n])/mc2
beta = sqrt(1-1/(gamma*gamma))
# scattered photon energy:
#\omega_1 = \gamma(1 +\beta\cos\theta'_1)\omega'
hw1 = gamma * ( 1 + beta*cos(theta0) )*hw0 # eV
dE[n] = dE[n] - hw1
# 1-turn tracking:
dE[n+1] = dE[n] + eVrf*(Z-Ne)*cos(phi[n])
phi[n+1] = phi[n] + (2*pi*h*eta/(beta_s*beta_s))*dE[n+1]/E_s
plt.scatter(phi/pi, 100*dE/E_s, s=2, linewidths=0, alpha=0.3)
plt.xlim(-0.5, +1.5)
plt.ylim(-0.04, +0.04)
plt.xlabel("Phase/pi")
plt.ylabel("dE/E (%)")
plt.grid()
plt.show()
Or in terms of longitudinal position with sespect to the synchronous particle $\Delta s = ... $
As we can see the number of photon scatterings per turn is limited by the RF voltage. With 5 scatterings per turn ion is always lost from the bucket in the case of 16 MV RF voltage (and these specific initial conditions).
The stability and efficiency can be improved if we use a laser with a higher spectral density at higher energy (i.e. a form of enhanced longitudinal cooling). For example, let's modify the script above, so that there is a small additional probaility of photon scattering for $\Delta E > 0$:
N_turns = 20000
dE = np.zeros(N_turns + 1) # energy deviation (eV)
phi = np.zeros(N_turns + 1) # ion phase with respect to RF
dE[0] = 1.0e-4*E_s # eV # typical value for the LHC
phi[0] = pi/2
eVrf = 16e6 # eV
for n in range(0, N_turns):
# photon scattering:
for sc in range(0, N_sc):
costheta = np.random.uniform(-1,+1)
# random polar angle in the ion's frame:
theta0 = acos(costheta)
gamma = (E_s + dE[n])/mc2
beta = sqrt(1-1/(gamma*gamma))
# scattered photon energy:
#\omega_1 = \gamma(1 +\beta\cos\theta'_1)\omega'
hw1 = gamma * ( 1 + beta*cos(theta0) )*hw0 # eV
dE[n] = dE[n] - hw1
if dE[n] > 0:
# +3% probability of photon scattering for dE > 0:
if np.random.uniform() < 0.03:
dE[n] = dE[n] - hw1 # let's take the previous hw1 for simplicity
# 1-turn tracking:
dE[n+1] = dE[n] + eVrf*(Z-Ne)*cos(phi[n])
phi[n+1] = phi[n] + (2*pi*h*eta/(beta_s*beta_s))*dE[n+1]/E_s
plt.scatter(phi/pi, 100*dE/E_s, s=2, linewidths=0, alpha=0.3)
plt.xlim(-0.5, +1.5)
plt.ylim(-0.04, +0.04)
plt.xlabel("Phase/pi")
plt.ylabel("dE/E (%)")
plt.grid()
plt.show()
Transverse effects
The transverse effect of laser cooling can be included after the longitudinal effects are understood. The transverse cooling is acheved because all components of ion momentum are lost due to the radiation friction force, but only the longitudinal component is restored by the RF system.
N_turns = 20000
dE = np.zeros(N_turns + 1) # energy deviation (eV)
phi = np.zeros(N_turns + 1) # ion phase with respect to RF
x = np.zeros(N_turns + 1) # m
xp = np.zeros(N_turns + 1) # rad
dE[0] = 1.0e-4*E_s # eV # typical value for the LHC
phi[0] = pi/2
x[0] = sigma_x/1e3
N_sc = 4
eVrf = 2*16e6 # eV
for n in range(0, N_turns):
# photon scattering:
for sc in range(0, N_sc):
costheta = np.random.uniform(-1,+1)
# random polar angle in the ion's frame:
theta0 = acos(costheta)
gamma = (E_s + dE[n])/mc2
beta = sqrt(1-1/(gamma*gamma))
# scattered photon energy:
#\omega_1 = \gamma(1 +\beta\cos\theta'_1)\omega'
hw1 = gamma * ( 1 + beta*cos(theta0) )*hw0 # eV
# polar angle of scattered photon in the lab frame:
sin_theta1 = sin(theta0) / ( gamma * (1 + beta*cos(theta0) ) )
cos_theta1 = ( beta + cos(theta0) ) / (1 + beta*cos(theta0) )
theta1 = atan2(sin_theta1, cos_theta1)
# take the effect of scattering into account here
# ...
# ...
# ...
dE[n] = dE[n] - hw1
# 1-turn tracking:
dE[n+1] = dE[n] + eVrf*(Z-Ne)*cos(phi[n])
phi[n+1] = phi[n] + (2*pi*h*eta/(beta_s*beta_s))*dE[n+1]/E_s
x[n+1] = x[n]*cos(mux) + xp[n]*beta_x*sin(mux)
xp[n+1] = -x[n]*sin(mux)/beta_x + xp[n]*cos(mux)
plt.scatter(phi/pi, 100*dE/E_s, s=2, linewidths=0, alpha=0.3)
plt.xlim(-0.5, +1.5)
plt.ylim(-0.04, +0.04)
plt.xlabel("Phase/pi")
plt.ylabel("dE/E (%)")
plt.grid()
plt.show()
plt.scatter(x*1e3, xp*1e3, s=2, linewidths=0, alpha=0.3)
plt.xlabel("x (mm)")
plt.ylabel("x' (mrad)")
plt.grid()
plt.show()
Some tests:¶
RF resonator changes only the $p_z$ component of ion momentum, while the photon scattering affects all components.
costheta = np.random.uniform(-1,+1)
# random polar angle in the ion's frame:
theta0 = acos(costheta)
theta0*180/pi