Laser interaction with a beam of partially stripped ions¶
A. Petrenko (CERN-Novosibirsk, 2017-2019)
This notebook describes the kinematics of photon absorption and emission by a relativistic partially stripped ion beam and application of this process to ion beam cooling in a storage ring.
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
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:
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:
Hence the scattered photon has the frequency $\omega_1 = \gamma(1 +\beta\cos\theta'_1)\omega' \approx 2\gamma^2(1 +\beta\cos\theta'_1) \omega$.
The typical $\theta_1 \sim 1/\gamma$.
Small fraction of photons is emitted with $\theta_1 \sim 1$ and in this case to separate forward and backward emission it's important to know $\cos\theta_1$:
$$ \omega_1 \cos\theta_1 = \omega'\gamma(\beta + \cos\theta_1') ~\Rightarrow~ \cos\theta_1 = \frac{\beta + \cos\theta_1'}{1 +\beta\cos\theta'_1} $$Let's consider Li-like Pb example¶
Consider Li-like Pb transition $2s-2p_{1/2}$. The data from: Calculations of 2p lifetimes in the Li sequence. C. E. Theodosiou, L. J. Curtis, and M. El-Mekki. Phys. Rev. A 44, 7144 (1991):
# Excitation energy:
hw0 = 230.823 # eV (ħω0)
# lifetime:
tau_0 = 76.6e-12 # sec
g1 = 2; g2 = 2 # degeneracy factors
import numpy as np
from IPython.display import Latex
import warnings
warnings.filterwarnings('ignore')
A = 207.98 # Atomic mass
mc = A*0.931e9 # eV/c
mc2 = mc # eV
Z = 82 # Number of protons in the ion (Lead)
Ne = 3 # Number of remaining electrons (Lithium-like)
Qe = 1.60217662e-19 # elementary charge in Coulombs
gamma_p = 236/0.938 # equvalent gamma for the protons in the SPS
gamma = gamma_p*(Z-Ne)*0.938e9/mc # Lorentz factor of the ion
gamma_0 = gamma # defined for later use below
beta_0 = np.sqrt(1-1/(gamma_0*gamma_0))
# Corresponding photon energy in the lab frame will be
hc = 0.19732697e-6 # eV*m (ħc)
lambda_0 = 2*np.pi*hc/hw0 # m
hw = hw0/(2*gamma) # eV (ħω)
# max energy of scattered photons will be
hw1 = 4*gamma*gamma*hw
Latex(r"""
\begin{aligned}
\gamma &= %.2f \\
\hbar \omega' &= %.1f~\rm{eV} ~~ (\lambda' = %.1f~ \mathrm{nm}) \\
\hbar \omega~ &= %.2f~\rm{eV} ~~ (\lambda = %.1f~ \mathrm{nm}) \\
\hbar \omega_{1, \rm{max}} &= %.1f~\rm{keV} \\
\end{aligned}
""" % (gamma, hw0, 1e9*lambda_0, hw, 1e9*2*np.pi*hc/hw, hw1/1e3))
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_0 # eV/c
pz = gamma_0*mc*beta_0 # eV/c
p0 = pz # defined for later use below
kick = (px/pz)*1000 # mrad
print('(px ~ %.2f keV, pz = %.2f TeV) Typical transverse kick ~ %.0e mrad' % (px/1e3, pz/1e12, kick))
Ion beam definition¶
N_ions = 0.9e8 # Number of ions in the beam
Np = 20000 # number of ion macro-particles in this simulation
Twiss parameters in the laser-ion interaction point (see the MAD-X SPS optics output for details):
beta_x = 54.614389 # m
alpha_x = -1.535235
beta_y = 44.332517 # m
alpha_y = 1.314101
Dx = 2.444732 # m
Dpx = 0.097522
Dy = 0.0 # m
Dpy = 0.0
L = 6911.5038 # m length of the ring
Ion beam parameters:
emitt_n = 1.5e-6 # m*rad (normalized emittance)
emitt = emitt_n/gamma
sigma_z = 0.063 # m
sigma_dp = 2e-4 # relative momentum spread
Distribution in $z$ and $\delta p = \frac{\Delta p}{p}$ can be defined easily since they are not correlated:
z = np.random.normal(scale=sigma_z, size=Np)
dp = np.random.normal(scale=sigma_dp, size=Np)
Transverse coordinates are more complicated. Similar to Elegant first we generate the beam in the normalized coordinates, defined as
$$ \begin{aligned} u~ &= \frac{x}{\sqrt{\beta_x}} \\ u' &= \frac{x'}{\sqrt{\beta_x}} + \frac{x\alpha_x}{\beta_x\sqrt{\beta_x}}. \end{aligned} $$And $\sigma_u = \sqrt{\epsilon_x}$, $\sigma_{u'} = \sigma_u/\beta_x$.
$u$ and $u'$ are uncorrelated, therefore it's easy to generate their distribution:
u = np.random.normal(scale=np.sqrt(emitt), size=Np)
up = np.random.normal(scale=np.sqrt(emitt)/beta_x, size=Np)
Now we need to apply the inverse transformation from $(u,u')$ to $(x,x')$: $$ \begin{aligned} x~ &= u\sqrt{\beta_x} \\ x' &= u'\sqrt{\beta_x} - u \frac{\alpha_x}{\sqrt{\beta_x}}. \end{aligned} $$
x = u*np.sqrt(beta_x) # m
xp = up*np.sqrt(beta_x) - u*alpha_x/np.sqrt(beta_x) # rad
The same for y-plane:
u = np.random.normal(scale=np.sqrt(emitt), size=Np)
up = np.random.normal(scale=np.sqrt(emitt)/beta_y, size=Np)
y = u*np.sqrt(beta_y) # m
yp = up*np.sqrt(beta_y) - u*alpha_y/np.sqrt(beta_y) # rad
Let's plot these distributions:
import holoviews as hv
hv.extension('matplotlib')
%output backend='matplotlib' fig='png' size=120 dpi=120
%opts Points Curve [show_grid=True aspect=2]
%opts Points (alpha=0.3 s=1)
#%opts Points (alpha=0.5) [aspect=2]
dim_x = hv.Dimension('x', unit='mm', range=(-8,+8))
dim_xp = hv.Dimension('xp', unit='mrad', label="x'", range=(-0.16,+0.16))
dim_y = hv.Dimension('y', unit='mm', range=(-8,+8))
dim_yp = hv.Dimension('yp', unit='mrad', label="y'", range=(-0.16,+0.16))
dim_z = hv.Dimension('z', unit='mm', range=(-250.0,+250.0))
#dim_dp = hv.Dimension('dp', label='100%*Δp/p')
dim_dp = hv.Dimension('dp', label='100%*$\Delta p/p$', range=(-0.15,+0.15))
hv.Points((y*1e3, yp*1e3), kdims=[dim_y,dim_yp])
hv.Points((z*1e3,dp*100), kdims=[dim_z,dim_dp])
As a crosscheck we can calculate the RMS-emittance of this beam using the statistical definition of emittance $\epsilon_{RMS} = \sqrt{\left<x^2\right> \left<x'^2\right> - \left<xx'\right>^2}$.
emitt_RMS = np.sqrt(np.mean(x*x)*np.mean(xp*xp) - np.mean(x*xp)*np.mean(x*xp))
#emitt_RMS = np.sqrt(np.mean(y*y)*np.mean(yp*yp) - np.mean(y*yp)*np.mean(y*yp))
Latex("Normalized $\epsilon_{RMS}$ = %.2f mm*mrad (compared to %.2f mm*mrad set originally)."
% (emitt_RMS*1e6*gamma, emitt_n*1e6))
$\epsilon_{RMS}$ will be closer to the input value with more macro-particles per bunch.
And finally the effect of transverse dispersion function can be included as $x = x + D_x \frac{\Delta p}{p}$ and $x' = x' + D'_x\frac{\Delta p}{p}$.
x = x + Dx*dp; xp = xp + Dpx*dp
y = y + Dy*dp; yp = yp + Dpy*dp
hv.Points((x*1e3,xp*1e3), kdims=[dim_x,dim_xp]) + hv.Points((y*1e3,yp*1e3), kdims=[dim_y,dim_yp])
hv.Points((dp*100,x*1e3), kdims=[dim_dp,dim_x])
Laser beam¶
To define the laser beam we define several parameters (assuming Gaussian beam):
- laser wavelength $\lambda_l$,
- a vector to the laser focal point $\boldsymbol{a}=(a_x,a_y,a_z)$,
- a vector of laser direction of propagation $\boldsymbol{n}=(n_x,n_y,n_z)$,
- laser waist radius $w_0$,
- delay of the laser pulse with respect to the ion beam (i.e. the laser pulse reaches its focal point at the moment $t=t_l$).
We assume that the laser pulse length is much shorter than the length of the ion bunch (with this approach a long laser pulse can be represented as a sequence of several short pulses). In this case for every particle we can calculate the moment and location of its collision with the laser wavefront.
# Let the laser beam collide with the ion beam at the center:
ax = 0; ay = 0; az = 0 # m
t_l = 0.0 # sec
# i.e. some ions collide with the laser at t < 0,
# before it reaches the center of the ion beam.
# So the tracking should start at some t < 0.
We assume that the laser crosses the ion beam trajectory with an angle $\theta_l$:
#theta_l = 10.0*np.pi/180 # rad
#theta_l = 5.0*np.pi/180 # rad
theta_l = 2.6*np.pi/180 # rad
#theta_l = 1.0*np.pi/180 # rad
#theta_l = 0.1*np.pi/180 # rad
Then $n_x = 0$, $n_y = -\sin\theta_l$, $n_z = -\cos\theta_l$.
nx = 0; ny = -np.sin(theta_l); nz = -np.cos(theta_l)
And because of this angle the resonant wavelength of the laser should be shifted from $\lambda_0$ as
$$ \lambda_l = \frac{\lambda_0}{2}(1 + \beta_0 \cos\theta_l). $$lambda_l = (2*np.pi*hc/hw)*(1 + beta_0*np.cos(theta_l))/2 # m
# Shift laser wavelength for fast longitudinal cooling
lambda_l = lambda_l*(1+sigma_dp)
Latex(r"$\lambda_l$ = %.1f nm" % (lambda_l*1e9))
w_0 = 1.3e-3 # m
#w_0 = 2.0e-3 # m
#w_0 = 4.0e-3 # m
#w_0 = 0.2e-3 # m
The corresponding Rayleigh length is $$ z_R = \frac{\pi w_0^2} {\lambda_l}. $$
z_R = np.pi*w_0*w_0/lambda_l
Latex("$z_R$ = %.2f m" % z_R)
Latex(r"Divergence of the laser beam $\theta = \frac{\lambda_l}{\pi w_0} = %.2f$ mrad." %
(1e3*lambda_l/(np.pi*w_0)))
Collision of the ion beam with the laser pulse¶
The position of the laser beam center is $\boldsymbol{r}_l = \boldsymbol{a} + c(t-t_l)\boldsymbol{n}$. We can find the moment when a particle with a position $\boldsymbol{r} = \boldsymbol{r}_0 + \boldsymbol{v}t$ collides with the laser as the moment when $\boldsymbol{r}-\boldsymbol{r}_l$ is perpendicular to $\boldsymbol{n}$. Then $(\boldsymbol{r}-\boldsymbol{r}_l, \boldsymbol{n})=0$, which yields the equation
Therefore the moment of ion-photon collision is $$ t_\mathrm{col} = \frac{(\boldsymbol{r}_0 - \boldsymbol{a}, \boldsymbol{n}) + ct_l}{c - (\boldsymbol{v}, \boldsymbol{n})}. $$
c = 299792458 # m/s
vz = c*beta_0
vx = vz*xp
vy = vz*yp
t_col = ((x-ax)*nx + (y-ay)*ny + (z-az)*nz + c*t_l) / (c - (vx*nx+vy*ny+vz*nz)) # sec
For every particle we can find the location of its collision with the laser pulse as $\boldsymbol{r} = \boldsymbol{r}_0 + \boldsymbol{v}t$:
x_col = x + vx*t_col # m
y_col = y + vy*t_col # m
z_col = z + vz*t_col # m
As a test let's plot the location of all ions which have already passed the laser pulse at $t=0$:
hv.Points((z[t_col<0]*1e3, y[t_col<0]*1e3), kdims=[dim_z,dim_y])
Now we need to determine if the photon is absorbed and the ion is excited according to the laser beam intensity and angle at the location of the interaction point.
At the moment of the ion collision with the laser pulse its intensity $I$ is given by the expression
$$ I = I_0 \frac{w^2_0}{w^2} \exp \left(-\frac{2(\boldsymbol{r} - \boldsymbol{r}_l)^2}{w^2}\right),$$where $\boldsymbol{r} - \boldsymbol{r}_l$ is the distance to the laser center, $I_0$ is the maximum intensity at the focal point $\boldsymbol{a}$, and the laser beam width $w$ is evolving with the distance to the focal point as
$$ w = w_0 \sqrt{1 + \left(\frac{\boldsymbol{r}_l - \boldsymbol{a}}{z_R}\right)^2}. $$Also the radius of wavefront curvature $R$ can be written as
$$ R = |\boldsymbol{r}_l - \boldsymbol{a}| \left[ 1 + \left(\frac{z_R}{\boldsymbol{r}_l - \boldsymbol{a}}\right)^2 \right]. $$The total laser pulse energy $J_l$ can be integrated as
$$ J_l = \Delta t_l \int\limits_{S} I dS = \Delta t_l I_0 \int\limits_{0}^{\infty} \exp \left(-\frac{2r^2}{w_0^2}\right) 2 \pi rdr = \Delta t_l \frac{I_0 \pi w_0^2}{2}, $$where $\Delta t_l$ is the laser pulse duration (assuming uniform temporal profile). The same expression assuming the Gaussian temporal pulse profile can be written as
To calculate the probabilistic ion excitation with a given process cross-section we only need the density of photons (and their incoming angle which is slightly varying due to the wavefront curvature) at the location of every ion. The total number of photons per pulse
$$ N_{\hbar \omega} = \frac{J_l} {\hbar \omega}. $$While the surface density of photons (number of photons per unit surface) is proportional to the laser intensity $I$:
$$ \frac{dN_{\hbar\omega}}{dS} = \frac{I\Delta t_l}{\hbar\omega} = \frac{2N_{\hbar\omega}}{\pi w^2} \exp \left(-\frac{2(\boldsymbol{r} - \boldsymbol{r}_l)^2}{w^2}\right). $$Photon absorption cross-section¶
The cross-section of the ion excitation by a photon with a frequency $\omega'$ is given by
$$ \sigma = 2\pi^2 c r_e f_{12} g(\omega' - \omega'_0), $$where $r_e$ is classical electron radius, $f_{12}$ is the oscillator strength, $\omega'_0$ is the resonance frequency of the ion transition, $g(\omega'-\omega'_0)$ is the Lorentzian
where $\Gamma$ is the resonance width defined by the lifetime of the excited ion $\tau_0$:
$$ \Gamma = \frac{1}{\tau_0}. $$Also
$$ \Gamma = 2 r_e \omega_0'^2 f_{12}\frac{g_1}{cg_2}. $$Therefore $$ \sigma(\omega' - \omega'_0) = \frac{\sigma_0}{1 + 4\tau^2_0(\omega' - \omega'_0)^2}, $$ where $$ \sigma_0 = \frac{\lambda_0'^2 g_2}{2\pi g_1}. $$
sigma_0 = lambda_0*lambda_0*g2/(2*np.pi*g1)
Latex(r"$\sigma_0$ = %.2e $\mathrm{m}^2$ = %.1f Gbarn" % (sigma_0, sigma_0/1e-28/1e9))
Let's plot $\sigma$ as a function of relative frequency error $(\omega' - \omega'_0)/\omega'_0$
omega_0 = 2*np.pi*c/lambda_0
domega = np.linspace(-1e-7,+1e-7,201)*omega_0
sigma = sigma_0/(1 + 4*(tau_0*domega)*(tau_0*domega))
dim_sigma = hv.Dimension('sigma', label="$\sigma/\sigma_0$", range=(0,1))
dim_domega = hv.Dimension('domega', label="$(\omega-\omega_0)/\omega_0$")
hv.Curve((domega/omega_0, sigma/sigma_0), kdims=[dim_domega], vdims=[dim_sigma])
First we can estimate the energy of the laser pulse required to excite every ion in the bunch.
If the laser wavelength perfectly matches the ion level we would need one photon per every $\sigma_0$ coss-section i.e.
$$ N_{\hbar \omega_0} \sim \frac{w_0^2}{\sigma_0}. $$However since the momentum spread $\Delta p/p$ in the ion bunch is larger than the width of the resonance we need $N_{\hbar \omega_0}$ photons for all possible frequencies, i.e.
$$ N_{\hbar \omega} \sim N_{\hbar \omega_0}\frac{\Delta p}{p}\bigg/\frac{\Gamma}{\omega_0} \sim \frac{w_0^2}{\sigma_0}\frac{\Delta p}{p}\omega_0\tau_0. $$And the energy of the laser pulse should be
$$ J_l \sim \hbar\omega \frac{w_0^2}{\sigma_0}\frac{\Delta p}{p}\omega_0\tau_0. $$This expression is an order of magnitude estimate neglecting geometric factors.
N_hw = (w_0*w_0/sigma_0)*sigma_dp*omega_0*tau_0
J_l = N_hw*hw*1.602e-19 # joules
Latex(r"$N_{\hbar\omega} \sim$ %.1e, \
laser pulse energy $\sim$ %.1f mJ, \
laser average power $\sim$ %.1f W." %
(N_hw, J_l*1e3, J_l/(L/c)))
Now let's set the specific laser pulse energy:
J_l = 5e-3 # joules
N_hw = J_l/(hw*1.602e-19)
Latex(r"$N_{\hbar\omega} =$ %.1e, \
laser pulse energy = %.1f mJ, \
laser average power = %.1f W." %
(N_hw, J_l*1e3, J_l/(L/c)))
For a Gaussian (in time) beam there is a simple relation between the pulse duration $\sigma_t$ and its frequency spread $\sigma_\omega$:
$$ \sigma_t \sigma_\omega = 1. $$Suppose we would like to cool the ion beam in the momentum space. For this we would need to interact only with the ions with $\Delta p > 0$ and with a Gaussian beam one can use, for example, $\sigma_\omega/\omega = 0.5 \sigma_{\Delta p/p}$ and the laser frequency shifted by $2\sigma_\omega$ from the resonanse frequency corresponding to the equillibrium ion momentum.
sigma_w = (c*hw/hc)*sigma_dp/2
Then the pulse duration can be found as $\sigma_t = 1/\sigma_\omega$
sigma_t = 1/sigma_w
Latex(r"$\sigma_t = $ %.1f ps." % (sigma_t*1e12))
Peak power density $I_0$ in the center of this laser pulse
$$ I_0 = \sqrt{\frac{2}{\pi}} \frac{J_l}{\sigma_t \pi w_0^2}, $$I_0 = np.sqrt(2.0/np.pi)*J_l/(sigma_t*np.pi*w_0*w_0)
Latex(r"$I_0 = $ %.2f TW/$\mathrm{m}^2$ = %.2f GW/$\mathrm{cm}^2$." % ((I_0*1e-12),(I_0*1e-9*1e-4)) )
Using these estimates as a reference we can now do the detailed Monte Carlo simulation including all geometric factors.
Monte Carlo simulation of the ion bunch interaction with the laser pulse¶
Above we have already determined the moment and location of collision of every ion with the laser pulse. Using the ion excitation cross-section, laser intensity and frequency distribution we can find the probability of ion excitation after the collision with the laser pulse.
Number of excitation events for every ion can be found as
$$ N_{exc} = \frac{dN_{\hbar \omega}}{dS} \bar\sigma = \frac{2N_{\hbar\omega}}{\pi w^2} \exp \left(-\frac{2(\boldsymbol{r} - \boldsymbol{r}_l)^2}{w^2}\right) \bar\sigma, $$where $\bar\sigma$ is the excitation cross-section averaged over the laser frequency distribution. And as we've seen already $\bar\sigma \ll \sigma_0$.
In our model here we assumed that if $N_{exc} \ll 1$, then $N_{exc}$ is the probability of ion excitation. For $N_{exc} \sim 0.1$ we should take into account the saturation effects. This can be done by solving the rate equation
$$ \frac{dP_2}{dt} = m_2 P_1 - m_2 B P_2, $$where $P_2, P_1$ are the probabilities of ion to be in the excited and non-excited state ($P_1 + P_2 = 1$). The rate of excitation events $m_2 P_1$ is proportional to the population of lower level $P_1$, photon density, and the absorption cross-section, $m_2 B P_2$ is the rate of stimulated emission events ($B = g_1/g_2$). Before ion enters into the laser pulse $P_1 = 1$ and $P_2 = 0$.
To solve the rate equation we need to separate the variables:
$$ \frac{dP_2}{1 - (1 + B)P_2} = m_2 dt. $$$m_2$ depends on time during the passage of the ion through the laser pulse, but we already know the answer in the case of small $P_2 = N_{exc} \ll 1$, hence
$$ \int\limits_t m_2(t) dt = \frac{dN_{\hbar \omega}}{dS} \bar\sigma. $$(In the general case of a long laser pulse maybe we would need to keep the $P_2$ number during the slice-by-slice integration of ion motion through the laser pulse).
Therefore the result of integration
$$ N_{exc} = P_2 = \frac{1 - \exp\left[-(1+B)\int\limits_t m_2(t) dt \right]}{1+B} = \frac{1 - \exp\left[-(1+B)\frac{dN_{\hbar \omega}}{dS} \bar\sigma \right]}{1+B}. $$As we can see while $\frac{dN_{\hbar \omega}}{dS} \bar\sigma$ is small it equals $N_{exc}$, while for large $\frac{dN_{\hbar \omega}}{dS} \bar\sigma$ the result is limited by $1/(1+B)$.
More accurate model can also include the time-dependent effects of Rabi oscillations -- to be done in the future.
Excitation cross section averaged over the whole laser spectrum will depend on the energy of the ion and on the angle of collision between the ion and the laser pulse.
$$ \bar\sigma = \int\limits_{-\infty}^{+\infty} \sigma(\omega') f(\omega')d\omega' = \int\limits_{-\infty}^{+\infty} \frac{\sigma_0}{1 + 4\tau_0^2\left( \omega' - \omega_0' \right)^2} f(\omega')d\omega' $$If laser frequency distribution $f(\omega')$ is much wider than the width of the resonance line then we can simply replace $f(\omega')$ with $f(\omega_0')$ and
$$ \bar\sigma \approx \sigma_0 f(\omega_0') \int\limits_{-\infty}^{+\infty} \frac{d\omega'}{1 + 4\tau_0^2\left( \omega' - \omega_0' \right)^2} = \sigma_0 f(\omega_0') \frac{\pi}{2\tau_0}. $$The Gaussian frequency distribution in the laser pulse with the main frequency $\omega_l$ is written as
$$ f(\omega) = \frac{1}{\sigma_\omega\sqrt{2\pi}}\exp\left[-\frac{(\omega - \omega_l)^2}{2\sigma_\omega^2}\right], $$therefore
The Doppler effect $\omega' = \gamma (1 + \beta \cos\theta) \omega$ applied to $\sigma_\omega = 1/\sigma_t$ and $\omega_l$ introduces the dependency on the ion momentum
$$ \sigma_{\omega'} = \gamma \frac{1 + \beta \cos\theta} {\sigma_t}, $$$$ \omega_l' = \gamma (1 + \beta \cos\theta) \omega_l, $$where $\theta = \theta_l + \delta\theta$, and $\delta\theta$ is the additional small angle due to the anglular spread in the ion beam and due to the curvature of the laser wavefront.
For the SPS $\delta\theta \sim 0.01~$mrad and its effect is important only at the final stage of extreme ion cooling in the SPS when the energy spread is reduced from $\Delta p/p\sim10^{-4}$ to $10^{-6}$ or less. We can put $\theta = \theta_l$ for now.
In our case $\sigma_t \ll \gamma\tau_0$, which leads to $\bar\sigma \ll \sigma_0$:
Latex(r"$\sigma_t = $ %.1f ps, while $\gamma\tau_0 = $ %.1f ps." % (sigma_t*1e12, gamma_0*tau_0*1e12))
Now we can calculate the excitation probability for every ion in the bunch.
We define Doppler factor $f_D = \gamma (1 + \beta \cos\theta)$.
f_D = gamma_0*(1+dp)*(1+beta_0*np.cos(theta_l))
Then $\omega_l'$
omega_l = 2*np.pi/(lambda_l/c)*f_D
sigma_average = sigma_0*sigma_t/(tau_0*f_D) * \
np.sqrt(2*np.pi)/4 * \
np.exp(-np.power( (omega_0 - omega_l)/(sigma_w*f_D) , 2)/2)
Position of the laser center at the time of collision $\boldsymbol{r}_l = \boldsymbol{a} + c(t_{\mathrm{col}}-t_l)\boldsymbol{n}$:
x_l = ax + c*(t_col-t_l)*nx
y_l = ay + c*(t_col-t_l)*ny
z_l = az + c*(t_col-t_l)*nz
Distance to the laser center r2 = $(\boldsymbol{r} - \boldsymbol{r}_l)^2$
r2 = (x_col-x_l)*(x_col-x_l) + (y_col-y_l)*(y_col-y_l) + (z_col-z_l)*(z_col-z_l)
Now let's find the laser pulse width
$$ w = w_0 \sqrt{1 + \left(\frac{\boldsymbol{r}_l - \boldsymbol{a}}{z_R}\right)^2}. $$w = w_0*np.sqrt( 1 + ((x_l-ax)*(x_l-ax) + (y_l-ay)*(y_l-ay) + (z_l-az)*(z_l-az))/(z_R*z_R) )
And finally number of excitation events for every ion can be found as
$$ \frac{dN_{\hbar \omega}}{dS} \bar\sigma = \frac{2N_{\hbar\omega}}{\pi w^2} \exp \left(-\frac{2(\boldsymbol{r} - \boldsymbol{r}_l)^2}{w^2}\right) \bar\sigma, $$$$ N_{exc} = \frac{1 - \exp\left[-(1+B)\frac{dN_{\hbar \omega}}{dS} \bar\sigma \right]}{1+B}. $$N_exc = 2*N_hw/(np.pi*w*w)*np.exp(-2*r2/w/w)*sigma_average
N_exc = (1 - np.exp(-1*(1+g1/g2)*N_exc))/(1+g1/g2)
#%output backend='matplotlib' fig='png' size=250
#%opts Points [aspect=3]
hv.Points((z*1e3, N_exc), [dim_z,"N_exc"])
Since we treat $N_{exc}$ as the probability of ion excitation, let's randomply select excited atoms according to their $N_{exc}$ value:
For every ion select a random number from 0 to 1:
rnd = np.random.uniform(size=Np)
hv.Points((z*1e3, rnd), [dim_z,"rnd"])
And we count the ion as excited if its N_exc value is larger than the rnd value:
Excited = N_exc > rnd
Plot excited ions:
%opts Points.Excited (alpha=0.8, color='black')
%opts Points.Still (alpha=0.3)
hv.Points((z*1e3, y*1e3), [dim_z,dim_y], group='Still') * \
hv.Points((z[Excited]*1e3, y[Excited]*1e3), [dim_z,dim_y], group='Excited')
hv.Points((z*1e3, dp*100), [dim_z,dim_dp], group='Still') * \
hv.Points((z[Excited]*1e3, dp[Excited]*100), [dim_z,dim_dp], group='Excited')
print("%.1f%% of all ions are excited." % (100*len(z[Excited])/len(z)))
Photon emissions¶
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.
# random polar angle in the ion's frame:
costheta = np.random.uniform(-1,+1, size=Np)
theta1p = np.arccos(costheta)
# emitted photon energy in the lab frame:
hw1_emitted = gamma_0*(1+dp) * ( 1 + beta_0*np.cos(theta1p) )*hw0 # eV
# only excited atoms are emitting of course:
hw1_emitted = hw1_emitted*Excited
dp = dp - hw1_emitted/(gamma_0*mc)
hv.Points((z*1e3, dp*100), [dim_z,dim_dp], group='Still') * \
hv.Points((z[Excited]*1e3, dp[Excited]*100), [dim_z,dim_dp], group='Excited')
Now let's summarize the above theory into the python function to simulate the process of ion excitations and photon emissions in a turn-by-turn simulation.
X = np.matrix([
x ,
xp,
y ,
yp,
z ,
dp
])
def ExciteIons(X):
# Returns a vector of ion states after interaction with the laser: Excited / Not excited.
x = X[0].A1; xp = X[1].A1; y = X[2].A1; yp = X[3].A1; z = X[4].A1; dp = X[5].A1
vz = c*beta_0
vx = vz*xp
vy = vz*yp
t_col = ((x-ax)*nx + (y-ay)*ny + (z-az)*nz + c*t_l) / (c - (vx*nx+vy*ny+vz*nz)) # sec
x_col = x + vx*t_col # m
y_col = y + vy*t_col # m
z_col = z + vz*t_col # m
f_D = gamma_0*(1+dp)*(1+beta_0*np.cos(theta_l))
omega_l = 2*np.pi/(lambda_l/c)*f_D
sigma_average = sigma_0*sigma_t/(tau_0*f_D) * \
np.sqrt(2*np.pi)/4 * \
np.exp(-np.power( (omega_0 - omega_l)/(sigma_w*f_D) , 2)/2)
x_l = ax + c*(t_col-t_l)*nx
y_l = ay + c*(t_col-t_l)*ny
z_l = az + c*(t_col-t_l)*nz
r2 = (x_col-x_l)*(x_col-x_l) + (y_col-y_l)*(y_col-y_l) + (z_col-z_l)*(z_col-z_l)
w = w_0*np.sqrt( 1 + ((x_l-ax)*(x_l-ax) + (y_l-ay)*(y_l-ay) + (z_l-az)*(z_l-az))/(z_R*z_R) )
N_exc = 2*N_hw/(np.pi*w*w)*np.exp(-2*r2/w/w)*sigma_average
N_exc = (1 - np.exp(-1*(1+g1/g2)*N_exc))/(1+g1/g2)
rnd = np.random.uniform(size=Np)
Excited = N_exc > rnd
return Excited
Excited = ExciteIons(X)
hv.Points((z*1e3, dp*100), [dim_z,dim_dp], group='Still') * \
hv.Points((z[Excited]*1e3, dp[Excited]*100), [dim_z,dim_dp], group='Excited')
def EmitPhotons(X, Excited):
# Returns vector of emitted photons and vector X after photon emissions
x = X[0].A1; xp = X[1].A1; y = X[2].A1; yp = X[3].A1; z = X[4].A1; dp = X[5].A1
Np_Excited = len(x[Excited])
# random polar angle in the ion's frame:
cos_theta1p = np.random.uniform(-1,+1, size=Np_Excited)
theta1p = np.arccos(cos_theta1p)
# emission angle in the lab frame:
sin_theta1 = np.sin(theta1p)/(gamma_0 * (1+dp[Excited]) * ( 1 + beta_0*np.cos(theta1p) ))
theta1 = np.arcsin(sin_theta1)
# azimuthal angle of emission:
phi1 = np.random.uniform(0,2*np.pi, size=Np_Excited)
# emitted photon energy in the lab frame:
hw1_emitted = gamma_0*(1+dp[Excited]) * ( 1 + beta_0*np.cos(theta1p) )*hw0 # eV
# duration of time between the ion excitation and photon emission (in the lab frame):
dt_excited = gamma_0*np.random.exponential(scale=tau_0, size=Np_Excited) # seconds
# distance travelled by the ion between its excitation and photon emission:
ds = dt_excited*c # meters
# emitted photon location:
s_photon = z[Excited] + ds
x_photon = x[Excited] + ds*xp[Excited]
y_photon = y[Excited] + ds*yp[Excited]
p_photon = hw1_emitted # eV/c
pz_photon = p_photon*np.cos(theta1)
pt_photon = p_photon*np.sin(theta1) # photon transverse momentum
px_photon = pt_photon*np.cos(phi1)
py_photon = pt_photon*np.sin(phi1)
xp_photon = px_photon/pz_photon # rad
yp_photon = py_photon/pz_photon # rad
# The effect of photon emission on ion is noticeble only in the dp value (see earlier explanation):
dp[Excited] = dp[Excited] - hw1_emitted/(gamma_0*mc)
Photons = np.matrix([
x_photon,
xp_photon,
y_photon,
yp_photon,
s_photon,
p_photon
])
X1 = np.matrix([
x ,
xp,
y ,
yp,
z ,
dp
])
return X1, Photons
X1, Photons = EmitPhotons(X, Excited)
#%output backend='matplotlib' fig='png' size=150
#%opts Points [aspect=2]
x_photon = Photons[0].A1
xp_photon = Photons[1].A1
s_photon = Photons[4].A1
p_photon = Photons[5].A1
dim_s = hv.Dimension('s', unit='m')# , range=(-600.0,+600.0))
dim_xp_photon = hv.Dimension('xp_photon', unit='mrad', label="x'", range=(-100,+100))
hv.Points((s_photon,x_photon*1e3), [dim_s,dim_x]) + \
hv.Points((xp_photon*1e3,x_photon*1e3), [dim_xp_photon,dim_x])
dim_p_photon = hv.Dimension('p_photon', unit='keV', label="Photon energy") #, range=(0,None))
hv.Points((s_photon,p_photon/1e3), [dim_s,dim_p_photon]) + \
hv.Points((xp_photon*1e3,p_photon/1e3), [dim_xp_photon,dim_p_photon])
And now we can do the full 1-turn tracking of the ion beam.
Transverse 1-turn matrix¶
The expressions will be written for the general coupled case (in x-y) but with matrix elements simplified for the uncoupled case. 1-turn matrix in the uncoupled case can be expressed in the Twiss parametrization.
Full 1-turn matrix (excluding RF-effects) can be written as
$$ \begin{pmatrix} x \\ x' \\ y \\ y' \\ \Delta l \\ \Delta p/p \end{pmatrix}_{n+1} = \begin{pmatrix} R_{11} & R_{12} & R_{13} & R_{14} & 0 & R_{16} \\ R_{21} & R_{22} & R_{23} & R_{24} & 0 & R_{26} \\ R_{31} & R_{32} & R_{33} & R_{34} & 0 & R_{36} \\ R_{41} & R_{42} & R_{43} & R_{44} & 0 & R_{46} \\ R_{51} & R_{52} & R_{53} & R_{54} & 1 & R_{56} \\ 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ x' \\ y \\ y' \\ 0 \\ \Delta p/p \end{pmatrix}_n, $$where $n$ is the turn number.
In the notation used here $\Delta l$ is the path length difference between the trajectory of the particle with these particular coordinates and the central particle. $\Delta l$ can be used to find the longitudinal particle coordinate along the bunch $z$ (used earlier). Since particles with different momentum move with different velocity
$$ z_{n+1} = z_n - \Delta l_{n+1} + T_s \Delta \upsilon, $$where $T_s$ is the revolution period and
$$ \Delta \upsilon = \upsilon_s \frac{1}{\gamma^2}\frac{\Delta p}{p}, $$where $\upsilon_s$ is the reference ion velocity. Therefore
$$ z_{n+1} = z_n - \Delta l_{n+1} + T_s \upsilon_s \frac{1}{\gamma^2}\frac{\Delta p}{p} = z_n - \Delta l_{n+1} + \frac{L}{\gamma^2}\frac{\Delta p}{p}, $$where $T_s \upsilon_s = L$ is the ring circumference.
So we can replace $\Delta l$ with $z$ in the above matrix 1-turn expression:
$$ \begin{pmatrix} x \\ x' \\ y \\ y' \\ z \\ \Delta p/p \end{pmatrix}_{n+1} = \begin{pmatrix} R_{11} & R_{12} & R_{13} & R_{14} & 0 & R_{16} \\ R_{21} & R_{22} & R_{23} & R_{24} & 0 & R_{26} \\ R_{31} & R_{32} & R_{33} & R_{34} & 0 & R_{36} \\ R_{41} & R_{42} & R_{43} & R_{44} & 0 & R_{46} \\ -R_{51} & -R_{52} & -R_{53} & -R_{54} & 1 & -R_{56} + L\gamma^{-2} \\ 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ x' \\ y \\ y' \\ z \\ \Delta p/p \end{pmatrix}_n. $$Effects of RF-resonators should be included separately because in our case the longitudinal nonlinearity of RF is often important.
The expression looks more simple in the uncoupled case:
$$ \begin{pmatrix} x \\ x' \\ y \\ y' \\ z \\ \Delta p/p \end{pmatrix}_{n+1} = \begin{pmatrix} R_{11} & R_{12} & 0 & 0 & 0 & R_{16} \\ R_{21} & R_{22} & 0 & 0 & 0 & R_{26} \\ 0 & 0 & R_{33} & R_{34} & 0 & 0 \\ 0 & 0 & R_{43} & R_{44} & 0 & 0 \\ -R_{51} & -R_{52} & 0 & 0 & 1 & -R_{56} + L\gamma^{-2} \\ 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ x' \\ y \\ y' \\ z \\ \Delta p/p \end{pmatrix}_n. $$Since the closed orbit with some momentum offset is given by
$$ \begin{pmatrix} x \\ x' \end{pmatrix} = \begin{pmatrix} D_x \\ D'_x \end{pmatrix} \frac{\Delta p}{p}, $$we can find the $R_{16}$ and $R_{26}$ matrix elements from the vales of dispersion functions defined earlier
$$ \begin{pmatrix} D_x \\ D'_x \end{pmatrix} \frac{\Delta p}{p} = \begin{pmatrix} R_{11} & R_{12} \\ R_{21} & R_{22} \end{pmatrix} \begin{pmatrix} D_x \\ D'_x \end{pmatrix} \frac{\Delta p}{p} + \begin{pmatrix} R_{16} \\ R_{26} \\ \end{pmatrix} \frac{\Delta p}{p}. $$therefore
$$ \begin{pmatrix} R_{16} \\ R_{26} \end{pmatrix} = \begin{pmatrix} D_x \\ D'_x \end{pmatrix} - \begin{pmatrix} R_{11} & R_{12} \\ R_{21} & R_{22} \end{pmatrix} \begin{pmatrix} D_x \\ D'_x \end{pmatrix}. $$In the uncoupled case the 1-turn matrix can be written in Twiss as
$$ \begin{pmatrix} R_{11} & R_{12} \\ R_{21} & R_{22} \\ \end{pmatrix} = \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}, $$where $\mu_x = 2\pi\nu_x$ ($\nu_x$ -- betatron tune), $\alpha_x = -\beta'_x / 2$, $\gamma_x = (1 + \alpha_x^2)/\beta_x$.
And the same for the y-plane.
The matrix elements $R_{51}$ and $R_{52}$ can be found from the symplecticity of the full-turn matrix $M$ as described in Analytical Tools in Accelerator Physics by V. Litvinenko (page 43, equation L3-53):
$$ -R_{51} = R_{16}R_{21} - R_{26}R_{11}, \\ -R_{52} = R_{16}R_{22} - R_{26}R_{12}. $$# For the SPS Q26 optics the betatron tunes are:
nux = 26.12999969
nuy = 26.17999991
mux = 2*np.pi*nux
muy = 2*np.pi*nuy
R11 = np.cos(mux) + alpha_x*np.sin(mux); R12 = beta_x*np.sin(mux)
R21 = -(1+alpha_x*alpha_x)*np.sin(mux)/beta_x; R22 = np.cos(mux) - alpha_x*np.sin(mux)
R33 = np.cos(muy) + alpha_y*np.sin(muy); R34 = beta_y*np.sin(muy)
R43 = -(1+alpha_y*alpha_y)*np.sin(muy)/beta_y; R44 = np.cos(muy) - alpha_y*np.sin(muy)
R16 = Dx - (Dx*R11 + Dpx*R12)
R26 = Dpx - (Dx*R21 + Dpx*R22)
R51 = R26*R11 - R16*R21
R52 = R26*R12 - R16*R22
gamma_t = 22.77422909 # gamma transition in the ring
R56 = L/(gamma_t*gamma_t)
M = np.matrix([
[ R11, R12, 0 , 0 , 0 , R16],
[ R21, R22, 0 , 0 , 0 , R26],
[ 0 , 0 , R33, R34, 0 , 0 ],
[ 0 , 0 , R43, R44, 0 , 0 ],
[-R51, -R52, 0 , 0 , 1 , -R56+L/(gamma_0*gamma_0)],
[ 0 , 0 , 0 , 0 , 0 , 1 ]
])
for row in M:
for itm in row.A1:
print("%+.3e " % itm, end="")
print("\n")
From SPS_optics:
-4.346e-01 +3.981e+01 +0.000e+00 +0.000e+00 +0.000e+00 -3.754e-01 -4.481e-02 +1.804e+00 +0.000e+00 +0.000e+00 +0.000e+00 +3.116e-02 +0.000e+00 +0.000e+00 +1.615e+00 +4.011e+01 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 -5.566e-02 -7.633e-01 +0.000e+00 +0.000e+00 +3.036e-02 -1.918e+00 +0.000e+00 +0.000e+00 +1.000e+00 -1.318e+01 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +0.000e+00 +1.000e+00
We can quickly test that the 1-turn transformation is not changing the beam distribution:
Initial distribution:
def plot_xxp_yyp(X):
x = X[0].A1 # .A1 converts matrix into 1D-array
xp = X[1].A1
y = X[2].A1
yp = X[3].A1
return hv.Points((x*1e3,xp*1e3), kdims=[dim_x,dim_xp]) + hv.Points((y*1e3,yp*1e3), kdims=[dim_y,dim_yp])
plot_xxp_yyp(X)
# test one turn:
X1 = M*X
# test many turns
#X1 = M*M*M*M*M*M*M*M*M*M*M*M*M*M*M*M*M*M*M*X
plot_xxp_yyp(X1)
Effects of RF-resonator¶
The transverse cooling happens because RF-resonator restores only the longitudinal component $p_z$ of the ion momentum, while the photon emission reduces the full momentum $p$. So, in order to simulate the transverse radiative cooling we need to take into account the resonator effect not only on $\Delta p/p$ but also on the angles $x'$ and $y'$.
Longitudinal momentum gain of the ion after it has passed through the RF-resonator depends on the ion phase with respect to the RF:
where $E_{\rm{RF}}$ is the accelerating electric field and $\phi$ is the ion phase in the RF resonator. The resulting longitudinal momentum change:
where $V_{\rm{RF}}$ is the RF-voltage, $\upsilon_z$ is the ion velocity.
Since $x' = p_x/p_z$ and $p_x$ is not changed in the resonator the angle $x'$ after the resonator
RF-resonator frequency $f_{\rm{RF}}$ is some high harmonic $h$ of ion revolution frequency:
$$ f_{\rm{RF}} = \frac{h}{T_s}, $$Longitudinal coordinate $z$ gives the longitudinal distance from the ion to the reference particle at the moment when the reference particle arrives at the RF-phase $\phi_0$ (which is always the same). So the ion then arrives to the RF-resonator after the time
$$ \Delta T = -\frac{z}{\upsilon_z} \approx -\frac{z}{\upsilon_s}. $$Then the ion phase in the RF-resonator is
$$ \phi = \phi_0 + 2\pi f_{\rm{RF}}\Delta T = \phi_0 - 2\pi \frac{hz}{T_s\upsilon_z} \approx \phi_0 - 2\pi \frac{hz}{L}. $$# Pass beam through the RF-resonator:
def RFcavity(X, h, eVrf, phi0):
# returns vector X after RF-cavity
x = X[0].A1; xp = X[1].A1; y = X[2].A1; yp = X[3].A1; z = X[4].A1; dp = X[5].A1
p = p0 + p0*dp # eV/c
# since p*p = px*px + py*py + pz*pz, and px = pz*x', py=pz*y':
pz = p/np.sqrt(1+xp*xp+yp*yp)
px = pz*xp
py = pz*yp
v0 = c/np.sqrt(1 + mc*mc/(p0*p0))
v = c/np.sqrt(1 + mc*mc/(p*p))
vz = v*pz/p
# phase in the resonator:
phi = phi0 - 2*np.pi*h*(z/L)*v0/vz
pz = pz + eVrf*(Z-Ne)*np.cos(phi)/beta_0 # pz after RF-cavity
xp = px/pz
yp = py/pz
p = np.sqrt(px*px+py*py+pz*pz)
dp = (p-p0)/p0 # new relative momentum deviation
return np.matrix([
x ,
xp,
y ,
yp,
z ,
dp
])
Multi-turn tracking¶
First we can check if our implementation of RF-cavity is changing the transverse beam emittance (because we've included higher-order effects).
N_turns = 100000
N_plots = 11
h = 4620 # SPS 200 MHz
eVrf = 7e6 # eV
phi0 = np.pi/2
t_plots = np.arange(0,N_turns+1,int(N_turns/(N_plots-1)))
X_plots = {}
X1 = X;
for turn in range(0,N_turns+1):
if turn in t_plots:
print( "\rturn = %g (%g %%)" % (turn, (100*turn/N_turns)), end="")
X_plots[turn] = X1
X1 = RFcavity(X1, h, eVrf, phi0)
X1 = M*X1
# function to calculate beam current profile:
def get_I(z, z_bin = 5e-3, z_min=-0.4, z_max=+0.4):
# z, z_bin, z_min, z_max in meters
hist, bins = np.histogram( z, range=(z_min, z_max), bins=int((z_max-z_min)/z_bin) )
Qm = Qe*(Z-Ne)*N_ions/Np # macroparticle charge in C
I = hist*Qm/(z_bin/c) # A
z_centers = (bins[:-1] + bins[1:]) / 2
return z_centers, I
%opts Area [show_grid=True aspect=2] (linewidth=1, alpha=0.7)
z_centers, I = get_I(z)
dim_I = hv.Dimension('I', unit='A', range=(0.0,+4.0))
z_I = hv.Area((z_centers*1e3, I), kdims=[dim_z], vdims=[dim_I])
z_I
# plot z-dp + beam current profile:
def plot_z_dp(X):
z = X[4].A1
dp = X[5].A1
z_centers, I = get_I(z)
img = (hv.Points((z*1e3, dp*100), [dim_z,dim_dp]) + \
hv.Area((z_centers*1e3, I), kdims=[dim_z], vdims=[dim_I]) ).cols(1)
return img
plot_z_dp(X_plots[t_plots[-1]])
items = [(turn*L/c, (plot_z_dp(X_plots[turn]))) for turn in t_plots]
m = hv.HoloMap(items, kdims = ['t (sec)'])
m.collate()
items = [(turn*L/c, plot_xxp_yyp(X_plots[turn])) for turn in t_plots]
m = hv.HoloMap(items, kdims = ['t (sec)'])
m.collate()
And finally the full turn-by-turn Monte Carlo simulation including photon emissions¶
import time
N_turns = 1000000
#N_turns = 2000000
N_plots = 21
t_plots = np.arange(0,N_turns+1,int(N_turns/(N_plots-1)))
X_plots = {}
Excited_plots = {}
X1 = X;
T_start = time.time()
for turn in range(0,N_turns+1):
Excited = ExciteIons(X1)
if turn in t_plots:
if turn != t_plots[0]:
sec_left = ((time.time() - T_start)/turn)*(N_turns - turn)
hours_left = np.floor(sec_left/3600)
min_left = (sec_left - hours_left*3600)/60
print( "\rturn = %g (%g %%, %g h. %.0f min. left)" %
(turn, (100*turn/N_turns), hours_left, min_left), end="")
X_plots[turn] = X1
Excited_plots[turn] = Excited
X1, Photons = EmitPhotons(X1, Excited)
X1 = RFcavity(X1, h, eVrf, phi0)
X1 = M*X1
#X_plots
np.save("X_plots.npy", X_plots)
def plot_z_dp_Excited(X,Excited):
z = X[4].A1
dp = X[5].A1
txt = "%.1f%% of ions excited" % (100*len(z[Excited])/len(z))
z_centers, I = get_I(z, z_bin = 1e-3)
img = (hv.Points((z*1e3, dp*100), [dim_z,dim_dp], group='Still') * \
hv.Points((z[Excited]*1e3, dp[Excited]*100), [dim_z,dim_dp], group='Excited', label=txt) + \
hv.Area((z_centers*1e3, I), kdims=[dim_z], vdims=[dim_I]) ).cols(1)
return img
#dim_I = hv.Dimension('I', unit='A', range=(0.0,+100.0))
dim_I = hv.Dimension('I', unit='A', range=(0.0,+8.0))
items = [(turn*L/c, plot_z_dp_Excited(X_plots[turn],Excited_plots[turn])) for turn in t_plots]
m = hv.HoloMap(items, kdims = ['t (sec)'])
m.collate()
items = [(turn*L/c, plot_xxp_yyp(X_plots[turn])) for turn in t_plots]
m = hv.HoloMap(items, kdims = ['t (sec)'])
m.collate()
def plot_z_y_Excited(X,Excited):
z = X[4].A1
y = X[2].A1
txt = "%.1f%% of ions excited" % (100*len(z[Excited])/len(z))
z_centers, I = get_I(z, z_bin = 1e-3)
img = (hv.Points((z*1e3, y*1e3), [dim_z,dim_y], group='Still') * \
hv.Points((z[Excited]*1e3, y[Excited]*1e3), [dim_z,dim_y], group='Excited', label=txt) + \
hv.Area((z_centers*1e3, I), kdims=[dim_z], vdims=[dim_I]) ).cols(1)
return img
items = [(turn*L/c, plot_z_y_Excited(X_plots[turn],Excited_plots[turn])) for turn in t_plots]
m = hv.HoloMap(items, kdims = ['t (sec)'])
m.collate()
Technical info:¶
%load_ext watermark
%watermark --python --date --iversions --machine
!jupyter nbconvert --to HTML psi_beam_vs_laser.ipynb
pwd