Intrabeam Scattering

The starting point to describe any change in emittances is

(1)
\begin{align} \Delta<g_a> = \int ds{\rm Tr}\left(G_a \frac{d\Sigma}{ds}\right) \end{align}

where $G_a$ is the invariant matrix and $\Sigma$ is the second moment matrix for the Gaussian distribution.
Thus, there are two pieces in order to understand the evolution of the emittances. One needs to know the definition of the invariants, and one needs to know the physics of the process changing the second moments- in this case IBS.

The evolution of the moment matrix by IBS is most easilly computed in the beam frame. We may also transform the invariants to the beam frame.

(2)
\begin{align} \frac{d \Sigma}{dt} = \mathcal{A}K_{ab} \end{align}

with

(3)
\begin{align} {\mathcal A} = \frac{N r_0^2 c}{8\pi^3 \beta^3 \epsilon_x \epsilon_y \epsilon_z} \end{align}

And, in the Coulomb log approximation (Bjorken-Mtingwa), we find

(4)
\begin{align} K_{ab} = Log \int d\Omega \frac{h_{ab}}{h_3} \end{align}

with $h_{ab} = {\hat r}_a {\hat r}_b - {\hat \Delta_a}{\hat \Delta_b}$, $h_3 = C_{ab}{\hat \Delta_a}{\hat \Delta_b}$, where ${\hat r}_a$ and ${\hat \Delta_a}$ are perpendicular unit vectors. For the angular integration, we have $d\Omega = \sin\theta d\theta d\phi$. We may parametrize them with

(5)
\begin{align} {\hat e}_1 = \pmatrix{\sin\theta \cos\phi\cr \sin\theta\sin\phi \cr \cos\theta},\ \ \ {\hat e}_2 =\pmatrix{\cos\theta\cos\phi\cr \cos\theta\sin\phi\cr -\sin\theta},\ \ \ {\hat e}_3 = \pmatrix{-\sin\phi \cr \cos\phi\cr 0} \end{align}

To compute $K_{ab}$ in Matlab, we need to do the 2-D numerical integral.
We use:

We connect the invariants to the beam distribution matrix via

(6)
\begin{align} \bar M_a ={1\over \epsilon_a}L^T B^T G_a B L \end{align}

We then write

(7)
\begin{align} \bar M_a = \pmatrix{A_a & B_a\cr B_a^T & C_a} \end{align}

The dispersion matrix is given by

(8)
\begin{align} \mathcal{B} = \pmatrix{1&0&0&0&0&-\eta_x\cr 0&1&0&0&0&-\eta'_x\cr 0&0&1&0&0&-\eta_y\cr 0&0&0&0&0&-\eta'_y\cr \eta'_x&-\eta_x&\eta'_y&-\eta_y&1&0\cr 0&0&0&0&0&1} \end{align}

This is for phase space ordered as (x,px,y,py,z,pz). With the ordering in AT, we need

(9)
\begin{align} \mathcal{B} = \pmatrix{1&0&0&0&-\eta_x&0\cr 0&1&0&0&-\eta'_x&0\cr 0&0&1&0&-\eta_y&0\cr 0&0&0&0&-\eta'_y&0\cr 0&0&0&0&1&0\cr -\eta'_x&\eta_x&-\eta'_y&\eta_y&0&1\cr } \end{align}

and the Lorentz transformation matrix is given by

(10)
\begin{align} L=\pmatrix{1&0&0&0&0&0\cr 0&1&0&0&0&0\cr 0&0&{1\over \gamma}&0&0&0\cr 0&0&0&{1\over P_0}&0&0\cr 0&0&0&0&{1\over P_0}&0\cr 0&0&0&0&0&{\gamma \over P_0}} \end{align}

In the uncoupled case, we have

(11)
\begin{align} C_x = \pmatrix{\beta_x &0 & -\gamma G_x \cr 0 & 0 & 0 \cr -\gamma G_x & 0 & \gamma^2 {\mathcal H}_x} \end{align}
(12)
\begin{align} C_y = \pmatrix{0&0 & 0\cr 0 & \beta_y & -\gamma G_y \cr 0&-\gamma G_y & \gamma^2{\mathcal H}_y} \end{align}
(13)
\begin{align} C_z = \pmatrix{\gamma_z\eta_x^2 & \gamma_z\eta_x\eta_y & -\alpha_z\gamma\eta_x \cr \gamma_z\eta_x \eta_y & \gamma_z\eta_y^2 & -\alpha_z\gamma\eta_y \cr -\alpha_z\gamma\eta_x & -\alpha_z \gamma\eta_y & \gamma^2 \beta_z} \end{align}

We define the IBS growth rate as ${1\over T_a} = {1\over <g_a>} \left({d<g_a>\over dt}\right)^{\rm IBS}$.

(14)
\begin{align} \frac{1}{T_a} = {\mathcal A} C^a_{ab} K_{ab} \end{align}

Once the rate is computed, the emittance may be evolved via the following equation

(15)
\begin{align} \frac{d\epsilon_a}{dt} = -\frac{2}{\tau_a}\left(\epsilon_a - \epsilon_{a0}\right) + \frac{\epsilon_a}{T_a} \end{align}

—-
The code from Demma is implemented in Mathematica via
F[NN_, Energy_, epsh0_, epsv0_, sigmas0_, sigmap0_, tauh_, tauv_,
taup_, reps_, lmax_] := Flatten[
{r0 = 2.81794092*10^(-15);
c = 2.99792458*10^8;
me = 510998.9026752615;
gamma = Energy/me;
Jeps = 2.00;
epsh = epsh0;
epsv = epsv0;
sigmas = sigmas0;
sigmap = sigmap0;
Hv = epsv0*(1 - reps)/sigmap0^2/Jeps;
eh[1] = epsh0;
ev[1] = epsv0;
sp[1] = sigmap0;
ss[1] = sigmas0;
Tpav[1] = 1;
Thav[1] = 1;
Tvav[1] = 1;
For[l = 2, l < lmax, l++,
sH = 1/Sqrt[1/sigmap^2 + h/epsh + Hv/epsv];

sigmav = Sqrt[epsv*by];

asb = Sqrt[bx*epsv/(by*epsh)];

gasb = asb^(0.021 - 0.044*Log[asb]);

Clg = Log[gamma^2*sigmav*epsh/r0/bx];
(* elegant definition of Coulomb log *)
(*sigmah=Sqrt[epsh*
bx + (dx*sigmap)^2];
ClgEle=Log[gamma^2*sigmah*epsh/r0/bx];*)
(*Clog=
NIntegrate[Clg[x],{x,0,Max[s]},MaxRecursion->20]/Max[s];*)

Clog = Apply[Plus, Clg*ds]/(Max[s] - Min[s]);

A = r0^2*c*NN*
Clog/(16*gamma^3*epsh^0.75*epsv^0.75*sigmas*sigmap^3);

Tpm1 = A*sH*gasb*(bx*by)^(-0.25);

(*Tpm1av=NIntegrate[Tpm1[x],{x,0,Max[s]},MaxRecursion->20]/Max[
s];*)
Tpm1av = Apply[Plus, Tpm1*ds]/Max[s];
Tpav[l] = Apply[Plus, Tpm1*ds]/Max[s];
hav[l] = Apply[Plus, h*ds]/Max[s];
(*HsTpav=NIntegrate[h[x]*Tpm1[x],{x,0,Max[s]},MaxRecursion->20]/
Max[s];*)
HsTpav = Apply[Plus, h*Tpm1*ds]/Max[s];
HvTpav = Apply[Plus, hy*Tpm1*ds]/Max[s];

Thm1av = sigmap^2*HsTpav/epsh;
Thav[l] = sigmap^2*HsTpav/epsh;
(*Tvm1av=sigmap^2*Hv*Tpm1av/epsv;
Tvav[l]=sigmap^2*Hv*Tpm1av/epsv;*)

Tvm1av = sigmap^2*HvTpav/epsv;
Tvav[l] = sigmap^2*HvTpav/epsv;

rtp = If[taup*Tpm1av < 1, taup*Tpm1av, 0.5];
rth = If[tauh*Thm1av < 1, tauh*Thm1av, 0.5];
rtv = If[tauv*Tvm1av < 1, tauv*Tvm1av, 0.5];
epsh1 = epsh0/(1 - rth);
epsv1 = epsv0*((1 - reps)/(1 - rtv) + reps/(1 - rth));
sigmap1 = Sqrt[sigmap0^2/(1 - rtp)];
sigmas1 = sigmas0*sigmap1/sigmap0;

eh[l] = epsh1;
ev[l] = epsv1;
sp[l] = sigmap1;
ss[l] = sigmas1;
epsh = (eh[l] + eh[l - 1])/2;
epsv = (ev[l] + ev[l - 1])/2;
sigmap = (sp[l] + sp[l - 1])/2;
sigmas = (ss[l] + ss[l - 1])/2;
rsh = epsh/epsh0;
rsv = epsv/epsv0;
rsp = sigmap/sigmap0;
];
{NN, epsh,
epsv,
sigmap*sigmas
(*sigmap,
sigmas*) (*, ClgEle,
rsh,
rsv,
rsp*)}
}]
—-
references:
Recent measurements at Cesr-TA
http://prst-ab.aps.org/abstract/PRSTAB/v16/i10/e104401
And the Bane approximation is here
http://prst-ab.aps.org/abstract/PRSTAB/v5/i8/e084403