Уравнение огибающей Капчинского—Владимирского для пучка заряженных частиц

В настоящей главе сначала рассматривается теория, а затем на Python выполняется расчёт огибающей электронного пучка в линейном ускорителе, состоящем из соленоидов и резонаторов. Выкладки следуют классическим монографиям по физике пучков заряженных частиц [32, 33]; подробности, опущенные здесь, приведены там.

Уравнения Максвелла

Объёмная плотность тока пучка \(\jmath = \rho\upsilon\), где \(\rho\) — объёмная плотность заряда, \(\upsilon\) — скорость пучка. Запишем дифференциальные уравнения Максвелла. $$ \nabla \vec{D} = 4\pi\rho, $$ $$ \nabla\times \vec{H} = \frac{4\pi\vec{\jmath}}{c}, $$ где \(\vec{D}\) — индукция электрического поля, \(\vec{H}\) — напряжённость магнитного поля, \(c\) — скорость света. Используем теорему Стокса об интегрировании дифференциальных форм, чтобы получить уравнения Максвелла, записанные в интегральной форме. $$ \oint\limits_{\partial V} \vec{D}\vec{dS} = 4\pi\int\limits_V\rho{dV}, $$ $$ \oint\limits_{\partial S} \vec{H}\vec{dl} = \frac{4\pi}{c}\int\limits_S\vec{\jmath}\vec{dS}. $$

Найдём \(D_r\) для цилиндрического пучка радиуса \(a\), заполненного зарядом с постоянной плотностью \(\rho_0\). $$ D_r = \frac{4\pi}{r}\int\limits^r_0\rho(\xi)\xi d\xi = \begin{equation*} \begin{cases} \displaystyle 2\pi\rho_0 r, r < a, \ \displaystyle \frac{2\pi\rho_0 a^2}{r}, r > a. \end{cases} \end{equation*} $$

Учитывая, что в вакууме \(D = E\), \(H = B\), \(E\) — напряжённость электрического поля, \(B\) — индукция магнитного поля, и в плоском пространстве в декартовой системе координат \(H_\alpha = \beta D_r\), где \(\displaystyle\beta = \frac{\upsilon}{c}\), радиальная компонента силы \(F_r\) из силы Лоренца \(\vec{F} = e\vec{E} + \displaystyle\frac{e}{c}\vec{\upsilon}\times\vec{B}\) записывается так. $$ F_r = eE_r - \frac{e\upsilon_z B_\alpha}{c} = eE_r(1-\displaystyle \frac{\upsilon^2}{c^2}) = \displaystyle \frac{eE_r}{\gamma^2}. $$ Полезно выразить поле через ток $$ I = \rho_0\upsilon\pi a^2, $$ тогда $$ E_r = \begin{equation*} \begin{cases} \displaystyle \frac{2 I r}{a^2\upsilon}, r < a, \ \displaystyle \frac{2 I}{r\upsilon}, r > a. \end{cases} \end{equation*} $$

Уравнения движения

Второй закон Ньютона \(\dot p_r = F_r\), используем параксиальное приближение, считая \(\gamma = const\). $$ \gamma m \ddot r = \displaystyle \frac{eE_r}{\gamma^2} = \displaystyle \frac{2 I e}{\gamma^2 a^2 \upsilon} r, $$ получено линейное уравнение, однако, поскольку движутся все частицы сразу, необходимо учесть, что \(a = a(t)\). Решение линейного уравнения, выписанного выше, можно представить как линейное преобразование фазовой плоскости, а поскольку отрезок на фазовой плоскости при невырожденном линейном преобразовании переходит в отрезок, его можно охарактеризовать одной точкой. Поэтому выберем крайнюю точку \(r = a\), которая и представляет крайнюю траекторию. $$ \gamma m \ddot r = \displaystyle \frac{2 I e}{\gamma^2 a \upsilon}. $$ Перейдём к дифференцированию по \(z\), учитывая, что \(\displaystyle dt = \frac{dz}{v}\), и тогда $$ a'' = \displaystyle \frac{e}{a}\frac{2I}{m\gamma^3\upsilon^3}. $$ Введём характерный альфвеновский ток \(I_a = \displaystyle \frac{mc^3}{e} \approx\) 17 кА, и тогда $$ a'' = \displaystyle \frac{2I}{I_a (\beta\gamma)^3} \frac{1}{a}. $$ Учтём внешнюю фокусировку, предполагая суперпозицию полей (что верно не всегда: в нелинейных средах суперпозиция не выполняется), и получим $$ a'' + k(z)a - \displaystyle \frac{2I}{I_a (\beta\gamma)^3} \frac{1}{a} = 0, $$ что напоминает уравнение огибающей $$ w'' + kw - \displaystyle \frac{1}{w^3} = 0 ,$$ где \( w =\displaystyle \sqrt \beta .\)

Уравнения огибающей для эллиптического пучка с распределением Капчинского—Владимирского с внешней фокусировкой линейными полями

Запишем распределение Капчинского—Владимирского. $$ f = A\delta(1 - \displaystyle\frac{\beta_x x'^2 + 2\alpha_x x x' + \gamma_x x^2}{\epsilon_x} - \displaystyle\frac{\beta_y y'^2 + 2\alpha_y y y' + \gamma_y y^2}{\epsilon_y} ), $$ где \(A\) — нормировочная постоянная, задаваемая полным числом частиц, а квадратичная форма \(\beta_x x'^2 + 2\alpha_x x x' + \gamma_x x^2\) и аналогичная форма по \(y\), стоящие под дельта-функцией, представляют собой инварианты Куранта—Снайдера, отнесённые к эмиттансам \(\epsilon_x\) и \(\epsilon_y\); дельта-функция помещает все частицы на поверхность постоянного инварианта. Полуоси эллипса равны $$ a = \sqrt{\epsilon_x \beta_x}, b = \sqrt{\epsilon_y \beta_y}. $$ Поле получается линейно внутри заряженного эллиптического цилиндра. $$ E_x = \displaystyle \frac{4I}{\upsilon}\frac{x}{a(a+b)}, $$ $$ E_y = \displaystyle \frac{4I}{\upsilon}\frac{y}{b(a+b)}. $$ Проверим, что \(\nabla \vec{E} = 4\pi\rho:\) $$ \displaystyle I = \rho \upsilon \pi ab, $$ $$ \nabla \vec{E} = \displaystyle \frac{4I(a+b)}{\upsilon(a+b)ab} = \displaystyle \frac{4I}{\upsilon ab} = 4\pi\rho. $$ Так как поля линейные, они складываются с полями фокусирующей линзы, рассчитанными отдельно. Подставим в уравнение огибающей \(\displaystyle a = \sqrt \epsilon_x w_x, b = \sqrt \epsilon_y w_y:\) $$ a'' + k_{xt} a - \frac{\epsilon_x^2}{a^3} = 0, $$ где \(k_{xt} = k_x - k_{xsc}\) — полная жёсткость, \(k_x\) — жёсткость линзы, а \(\displaystyle k_{xsc} = \frac{4I}{I_a (\beta\gamma)^3}\frac{1}{a(a+b)}.\) В итоге получаем систему уравнений, связанных через пространственный заряд. $$ \begin{equation*} \begin{cases} \displaystyle a'' + k_xa - \frac{4I}{I_a (\beta\gamma)^3}\frac{1}{(a+b)} - \frac{\epsilon_x^2}{a^3} = 0 , \\ \displaystyle b'' + k_yb - \frac{4I}{I_a (\beta\gamma)^3}\frac{1}{(a+b)} - \frac{\epsilon_y^2}{b^3} = 0. \end{cases} \end{equation*} $$

Количественный критерий применимости приближения ламинарности течения

Эта система учитывает два эффекта, мешающих сфокусировать пучок в точку: конечность эмиттанса и пространственный заряд. В уравнение они входят соседними членами, поэтому допускают непосредственное сравнение. При малом токе \(I\) отталкивание слабое, при большом — сильное, и тогда эмиттансом можно пренебречь, а течение считать ламинарным. Количественным критерием ламинарности является сравнение этих двух членов. $$ \displaystyle \sqrt{\frac{2I}{I_a(\beta\gamma)^3}} \gg \sqrt{\frac{\epsilon_x}{\beta_x}}. $$ Видно, что пространственный заряд сильнее сказывается там, где \(\beta_x\)-функция больше, а вблизи фокуса его влияние пренебрежимо мало.

Теорема Буша

Теорема Буша связывает угловую скорость заряженной частицы, движущейся в аксиально-симметричном магнитном поле, с магнитным потоком через окружность, проведённую через точку, где находится частица, с центром на оси.

Рассмотрим заряд \(q\), движущийся в магнитном поле \(\vec B = (B_r,0,B_z)\). Приравняем \(\theta\)-составляющую силы Лоренца к производной момента импульса по времени, делённой на \(r\). $$ F_\theta = -q(\dot r B_z - \dot z B_r) = \frac{d}{rdt}(\gamma m r^2 \dot \theta). $$ Поток, пронизывающий площадь, охваченную окружностью радиуса \(r\), центр которой расположен на оси, а сама она проходит через точку, в которой расположен заряд, записывается в виде \(\psi = \int\limits^r_0 2\pi r B_z dr\). Когда частица перемещается на \(\vec{dl} = (dr,dz)\), скорость изменения потока, охваченного этой окружностью, можно найти из второго уравнения Максвелла \(\nabla \vec{B} = 0.\) Таким образом, $$ \dot\psi = 2\pi r (-B_r \dot z + B_z \dot r). $$ После интегрирования по времени отсюда получаем $$ \dot \theta = (-\frac{q}{2\pi \gamma m r^2})(\psi - \psi_0). $$

Уравнение параксиального луча

Чтобы вывести уравнение параксиального луча для системы с аксиальной симметрией при уже сделанных допущениях, приравняем силу радиального ускорения электрическим и магнитным силам внешних полей. Здесь \(\gamma = \gamma(t).\) $$ \frac{d}{dt}(\gamma m \dot r) - \gamma m r (\dot \theta)^2 = q(E_r + r\dot \theta B_z). $$ Применим теорему Буша и независимость \(B_z\) от \(r:\) $$ -\dot \theta = \frac{q}{2\gamma m}(B_z - \frac{\psi_0}{\pi r^2}). $$ Исключим \(\dot \theta\) и подставим \(\displaystyle \dot \gamma \approx \frac {\beta q E_z }{mc}\). $$ \ddot r + \frac{\beta q E_z}{\gamma m c} \dot r + \frac{q^2 B^2_z}{4\gamma ^2 m^2} r - \frac{ q^2 \psi^2_0}{4\pi^2\gamma^2 m^2} (\frac{1}{r^3}) - \frac{q E_r}{\gamma m } = 0. $$

Уравнение огибающей для круглого и эллиптического пучка

Учитывая, что $$ \dot r = \beta c r', $$ $$ \ddot r = r'' (\dot z)^2 + r'\ddot z \approx r'' \beta^2 c^2 + r' \beta' \beta c^2. $$ Если же в области пучка зарядов нет, то, разлагая в ряд Тейлора в окрестности оси и оставляя только первый член, с учётом $$\nabla \vec{E} = 0$$ получаем $$ E_r = -0.5 r E'_z \approx - 0.5 r \gamma'' mc^2/q. $$ Теперь можно окончательно записать уравнение огибающей круглого пучка радиуса \(r\), взятого с распределением Капчинского—Владимирского, при внешней фокусировке линейными полями. $$ \displaystyle r'' + \frac{1}{\beta^2\gamma} \gamma' r' + \frac{1}{2\beta^2\gamma}\gamma''r + kr - \frac{2I}{I_a (\beta\gamma)^3}\frac{1}{r} - \frac{\epsilon^2}{r^3} = 0 ; $$ и для эллиптического пучка $$ \begin{equation*} \begin{cases} \displaystyle a'' + \frac{1}{\beta^2\gamma} \gamma' a' + \frac{1}{2\beta^2\gamma}\gamma''a + k_xa - \frac{4I}{I_a (\beta\gamma)^3}\frac{1}{(a+b)} - \frac{\epsilon_x^2}{a^3} = 0 , \\ \displaystyle b'' + \frac{1}{\beta^2\gamma} \gamma' b' + \frac{1}{2\beta^2\gamma}\gamma''b + k_yb - \frac{4I}{I_a (\beta\gamma)^3}\frac{1}{(a+b)} - \frac{\epsilon_y^2}{b^3} = 0. \end{cases} \end{equation*} $$

Магнитные линзы

Фокусировка пучка осуществляется соленоидами и магнитными квадрупольными линзами.

Соленоиды

\(k_x = k_y = k_s\) — жёсткость соленоида. $$ k_s = \left ( \frac{eB_z}{2m_ec\beta\gamma} \right )^2 = \left ( \frac{e B_z}{2\beta\gamma\cdot 0.511\cdot 10^6 e \cdot \mathrm{volt}/c} \right )^2 = \left ( \frac{cB_z[\mathrm{T}]}{2\beta\gamma\cdot 0.511\cdot 10^6 \cdot \mathrm{volt}} \right )^2. $$

Квадруполи

\(k_q = \displaystyle\frac{eG}{p}\) — жёсткость квадруполя, где \(G = \displaystyle\frac{\partial B_x}{\partial y} = \displaystyle\frac{\partial B_y}{\partial x}\) — градиент магнитного поля, причём \(k_x = k_q, k_y = -k_q.\) $$ k_q = \left ( \frac{eG}{m_ec\beta\gamma} \right ) = \left ( \frac{eG}{\beta\gamma\cdot 0.511\cdot 10^6 e \cdot \mathrm{volt}/c} \right ) = \left ( \frac{cG}{\beta\gamma\cdot 0.511\cdot 10^6 \cdot \mathrm{volt}} \right ). $$

Продольная динамика пучка

Продольную динамику удобно рассчитывать отдельно от поперечной, находя сначала энергию пучка как функцию \(z\), а затем подставляя её в уравнение огибающей в виде готовой функции. Скорость электрона достаточно близка к скорости света, поэтому его продольная координата \(z \approx ct\), а импульс \(p_z \approx \gamma mc\), и тогда

$$ \frac{d\gamma}{dz} \approx \frac{eE_z}{mc^2}, $$

откуда задача сводится к однократному интегрированию \(E_z(z)\) вдоль тракта. Численные значения подставляются в разделе с расчётом на Python.

Решение уравнения огибающей для эллиптического пучка с фокусирующими элементами

Уравнение огибающей эллиптического пучка с полуосями \(a, b\) и распределением Капчинского—Владимирского при внешней фокусировке линейными полями, выведенное выше, имеет следующий вид. $$ \begin{equation*} \begin{cases} \displaystyle a'' + \frac{1}{\beta^2\gamma} \gamma' a' + \frac{1}{2\beta^2\gamma}\gamma''a + k_qa - \frac{2P}{(a+b)} - \frac{\epsilon_x^2}{a^3} = 0 , \\ \displaystyle b'' + \frac{1}{\beta^2\gamma} \gamma' b' + \frac{1}{2\beta^2\gamma}\gamma''b - k_qb - \frac{2P}{(a+b)} - \frac{\epsilon_y^2}{b^3} = 0. \end{cases} \end{equation*} $$

Пусть \(\displaystyle x = \frac{da}{dz}, y = \frac{db}{dz}, \displaystyle \frac{d\gamma }{dz}\approx \frac{e E_z}{m c^2},\) тогда $$ \begin{matrix} \displaystyle \frac{dx}{dz}= - \frac{1}{\beta^2\gamma} \gamma' a' - \frac{1}{2\beta^2\gamma}\gamma''a - k_qa + \frac{2P}{(a+b)} + \frac{\epsilon_x^2}{a^3}\\ \displaystyle\frac{da}{dz} = x\\ \displaystyle \frac{dy}{dz}= - \frac{1}{\beta^2\gamma} \gamma' b' - \frac{1}{2\beta^2\gamma}\gamma''b + k_qb + \frac{2P}{(a+b)} + \frac{\epsilon_y^2}{b^3}\\ \displaystyle\frac{db}{dz} = y \end{matrix}. $$ Пусть

$$ \vec X = \begin{bmatrix} x \\ a \\ y \\ b \end{bmatrix}, $$

теперь составим дифференциальное уравнение \(X' = F(X)\).

Уравнение траектории для одной частицы

Векторный потенциал в соленоиде имеет только одну компоненту и выражается через заданное на оси поле \(B_z(z)\). $$ A_\phi (z,r)=\sum_{n=0}^{\infty} \dfrac{(-1)^n}{n!(n+1)!} B_z^{(2n)}(z)\left(\dfrac{r}{2}\right)^{2n+1}. $$ Здесь \(B_z^{(2n)}\) — производная порядка \(2n\) от поля на оси, так что нулевой член даёт \(A_\phi = B_z(z),r/2\), из которого дальше и получается параксиальное приближение. Следовательно, лагранжиан в цилиндрической системе координат имеет вид $$ L = -mc^2(1-\beta^2)^{1/2} + \dfrac{e}{c}A_\phi r \dot\phi. $$ Из аксиальной симметрии следует существование интеграла движения, связанного с обобщённым импульсом. $$ P_\phi = \gamma mr^2\dot\phi + \dfrac{e}{c} r A_\phi = \gamma mr^2\dot\phi_0. $$ Продифференцировав по \(z\), получим $$ \phi' - \phi_0' = -\dfrac{e}{pc}\dfrac{A_\phi}{r}. $$ И в параксиальном приближении $$ \phi' - \phi_0' = -\dfrac{e}{2pc}B_z(z). $$ Откуда $$ \delta \phi = -\int\dfrac{e}{2pc}B_z(z)dz. $$ Тогда и уравнения траектории в приведённых координатах принимают следующий вид. $$ \begin{equation*} \begin{cases} \displaystyle \tilde x'' + \frac{1}{\beta^2\gamma} \gamma'\tilde x' + \frac{1}{2\beta^2\gamma}\gamma''\tilde x + k_s \tilde x = 0 , \\ \displaystyle \tilde y'' + \frac{1}{\beta^2\gamma} \gamma'\tilde y' + \frac{1}{2\beta^2\gamma}\gamma''\tilde y +k_s \tilde y = 0, \end{cases} \end{equation*} $$

Расчёт в Python

Подключим библиотеки и настроим оформление графиков. Построение графиков осуществляется средствами HoloViews поверх matplotlib, в векторном формате и с широкими осями, поскольку все иллюстрации главы представляют собой протяжённые графики вдоль тракта.

import numpy as np
import holoviews as hv
hv.extension('matplotlib')

%output size=200 backend='matplotlib' fig='svg'

%opts Curve Scatter Area [fontsize={'title':14, 'xlabel':14, 'ylabel':14, 'ticks':14, 'legend': 14}]

%opts Area Curve [aspect=3 show_grid=True]
%opts Area  (linewidth=1 alpha=0.25)
%opts Curve (linewidth=2 alpha=0.5)

import warnings
warnings.filterwarnings('ignore')

Зададим область счёта и шаг сетки вдоль оси ускорителя, от 0.7 до 15 метров с шагом 5 мм: на характерную ширину поля соленоида такой шаг даёт десятки точек, что позволяет безопасно интерполировать по ней.

dz = 0.005 # m
z_min = 0.7
z_max = 15
z = np.arange(z_min,z_max,dz) # m

Ускоритель представляется последовательностью элементов — соленоидов и ускоряющих зазоров — на общей оси. Каждый элемент характеризуется четырьмя параметрами: положением центра, максимумом поля, файлом с профилем поля вдоль оси и именем.

class Element:
  def __init__(self, z0, MaxField, filename, name):
    self.z0 = z0
    self.MaxField = MaxField
    self.filename = filename
    self.name = name

Первыми задаются фокусирующие соленоиды. В таблице ниже первая колонка задаёт положение центра в метрах, вторая — максимум поля на оси в теслах. Три соленоида работают на полном поле (0.03 Тл), три практически выключены (1 мТл). Профиль \(B_z(z)\) у всех соленоидов один и тот же и хранится в общем файле: соленоиды одинаковы и различаются только током.

Bz_beamline = {}
for   z0,         B0,        filename,       name in [
    # m           T                          Unique name
    [ 0.95,       0.001,    'input/fields/B_z(z).dat', 'Sol. 1' ],
    [ 2.1,        0.03,     'input/fields/B_z(z).dat', 'Sol. 2' ],
    [ 2.9077,     0.001,    'input/fields/B_z(z).dat', 'Sol. 3' ],
    [ 4.0024,     0.03,     'input/fields/B_z(z).dat', 'Sol. 4' ],
    [ 5.642,      0.03,     'input/fields/B_z(z).dat', 'Sol. 5' ],
    [ 6.760,      0.001,    'input/fields/B_z(z).dat', 'Sol. 6' ],
]:
    Bz_beamline[name] = Element(z0, B0, filename, name)

Элементы, собранные в словарь по имени, доступны напрямую: к любому из них можно обратиться по ключу.

Bz_beamline['Sol. 3'].z0

2.9077

Аналогично описываются ускоряющие резонаторы. Шесть резонаторов, каждый по −0.9 МВ/м, расположены через 1.382 метра во второй половине тракта; до седьмого метра пучок с входной энергией 2 МэВ только транспортируется, ускорение начинается дальше.

Ez_beamline = {}
for   z0,          E0,       filename,      name in [
    # m            MV/m                     Unique name
    [ 7.456,       -0.9,     'input/fields/E_z(z).dat',  'Cavity 3'],
    [ 8.838,       -0.9,     'input/fields/E_z(z).dat',  'Cavity 4'],
    [ 10.220,      -0.9,     'input/fields/E_z(z).dat',  'Cavity 5'],
    [ 11.602,      -0.9,     'input/fields/E_z(z).dat',  'Cavity 6'],
    [ 12.984,      -0.9,     'input/fields/E_z(z).dat',  'Cavity 7'],
    [ 14.366,      -0.9,     'input/fields/E_z(z).dat',  'Cavity 8'], 
]:
    Ez_beamline[name] = Element(z0, E0, filename, name)

Считывание профилей фокусирующего и ускоряющего поля

Профиль поля одного элемента хранится в текстовом файле в виде двух колонок: смещения от центра элемента и значения поля в долях максимума. Приведённая ниже функция проходит по всем элементам линейки, сдвигает профиль в положение элемента, умножает его на максимальное поле элемента и суммирует вклады в один массив на общей сетке.

FieldFiles представляет собой простейший кеш: файл читается с диска один раз, а не по разу на каждый соленоид. fill_value=(0, 0) вместе с bounds_error=False заставляет интерполятор возвращать ноль за пределами профиля вместо возбуждения ошибки; без этого было бы невозможно суммировать элементы, разнесённые по тракту.

from scipy import interpolate

FieldFiles = {} # buffer for field files
def read_field(beamline, z):
    global FieldFiles
    F = 0
    for element in beamline.values():
        if not (element.filename in FieldFiles):
            print('Reading ' + element.filename)
            FieldFiles[element.filename] = np.loadtxt(element.filename)
        M = FieldFiles[element.filename]
        z_data = M[:,0]
        F_data = M[:,1]
        f = interpolate.interp1d(
            element.z0+z_data, element.MaxField*F_data,
            fill_value=(0, 0), bounds_error=False
        )
        F = F + f(z)
    return F

Рассчитаем суммарное поле шести соленоидов.

Bz = read_field(Bz_beamline, z)

Reading input/fields/B_z(z).dat

Строка «Reading …» напечатана один раз на шесть соленоидов, следовательно, кеш работает.

Решателю дифференциальных уравнений требуется не массив, а функция: он самостоятельно выбирает точки, в которых запрашивается поле, и эти точки не обязаны совпадать с узлами сетки. Поэтому массив Bz преобразуется в функцию Bz(z).

Bz = interpolate.interp1d(z, Bz, fill_value=(0, 0), bounds_error=False)

Проверим значение Bz вблизи центра второго соленоида.

Bz(2.11111)

array(0.02984419)

Получено около 30 мТл, что близко к максимуму второго соленоида, настроенного на 0.03 Тл, как и должно быть в 11 миллиметрах от его центра.

График суммарного поля \(B_z(z)\) позволяет убедиться, что соседние соленоиды не накладываются друг на друга и что в тракте не осталось протяжённых участков без фокусировки.

dim_z  = hv.Dimension('z',  unit='m', range=(0,z_max))
dim_Bz = hv.Dimension('Bz', unit='T', label='$B_z$')

z_Bz = hv.Area((z,Bz(z)), kdims=dim_z, vdims=dim_Bz)
z_Bz

С ускоряющим полем поступаем аналогично: профиль считывается из файла, вклады шести резонаторов суммируются, результат оборачивается в интерполянт и выводится на график. Функция read_field применима как к магнитному, так и к электрическому полю.

Ez = 1*read_field(Ez_beamline, z)

Ez = interpolate.interp1d(z, Ez, fill_value=(0, 0), bounds_error=False)

dim_Ez = hv.Dimension('Ez', unit='MV/m', label=r'$E_z$')

z_Ez = hv.Area((z,Ez(z)), kdims=dim_z, vdims=dim_Ez)

z_Ez

Продольная динамика пучка: расчёт

Энергия пучка рассчитывается независимо от огибающей, а в уравнение огибающей подставляется готовая функция \(\gamma(z)\). Скорость электрона близка к скорости света, поэтому \(z \approx ct\) и \(p_z \approx \gamma mc\), откуда следует

$$ \frac{d\gamma}{dz} \approx \frac{eE_z}{mc^2}, $$

и задача сводится к однократному интегрированию \(E_z(z)\).

Энергия на входе составляет 2 МэВ; релятивистский фактор при этом равен около 4.9. mc — энергия покоя электрона в мегаэлектронвольтах, поэтому энергии и импульсы далее выражаются в тех же единицах.

mc = 0.511 # MeV, энергия покоя электрона
#p_z = 2.459 # MeV/c
E_0 = 2.000 # MeV
gamma_0 = E_0/mc + 1 # Начальный релятивистский фактор пучка
#gamma_0 = p_z/mc

Интегрирование выполняется функцией cumulative_trapezoid: она вычисляет накопленную сумму трапеций и даёт \(\gamma\) во всех точках сетки, кроме первой. Знак минус перед Ez компенсирует знак поля в резонаторах (-0.9 МВ/м), так что энергия по тракту растёт.

from scipy import integrate

gamma = gamma_0 + integrate.cumulative_trapezoid(-Ez(z)/mc, z)

dim_gamma = hv.Dimension('gamma', label=r'$\gamma$', range=(0,None))

z_gamma = hv.Curve((z,gamma), kdims=dim_z, vdims=dim_gamma)

Совместим на одном рисунке энергию и ускоряющее поле, что позволяет увидеть, где пучок набирает энергию.

(z_gamma + z_Ez).cols(1)

График \(\gamma(z)\) поднимается ступенями, по одной на резонатор, и остаётся постоянным в промежутках, где пучок движется по инерции.

Как и поля, массив gamma преобразуется в функцию. Срез z[1:] используется потому, что интегрирование по трапециям даёт на одну точку меньше, чем содержит сетка.

gamma = interpolate.interp1d(z[1:], gamma, fill_value=(gamma[0], gamma[-1]), bounds_error=False)

#hv.Curve((z,gamma(z)), kdims=dim_z, vdims=dim_gamma)

Решение уравнения огибающей для круглого пучка

Повторим уравнение огибающей для круглого пучка радиуса \(r\) с распределением Капчинского—Владимирского при внешней фокусировке линейными полями. $$ \displaystyle r'' + \frac{1}{\beta^2\gamma} \gamma' r' + \frac{1}{2\beta^2\gamma}\gamma''r + kr - \frac{2I}{I_a (\beta\gamma)^3}\frac{1}{r} - \frac{\epsilon^2}{r^3} = 0. $$

Перейдём к среднеквадратичному радиусу для распределения Капчинского—Владимирского и первеансу \(\displaystyle r_{rms} = \frac{r}{2}\) и \(\displaystyle P = \frac{2I}{I_a\beta^3\gamma^3} \). Получим $$ \displaystyle r_{rms}'' + \frac{1}{\beta^2\gamma} \gamma' r_{rms}' + \frac{1}{2\beta^2\gamma}\gamma''r_{rms} + kr_{rms} - \frac{P}{4r_{rms}} - \frac{\epsilon^2}{16{r_{rms}}^3} = 0. $$

Пусть \(\displaystyle y=\frac{dr_{rms}}{dz}, \displaystyle \frac{d\gamma }{dz}\approx \frac{e E_z}{m c^2}\), тогда $$ \begin{matrix} \displaystyle \frac{dy}{dz}= - \frac{1}{\beta^2\gamma} \gamma' r_{rms}' - \frac{1}{2\beta^2\gamma}\gamma''r_{rms} - kr_{rms} + \frac{P}{4r_{rms}} + \frac{\epsilon^2}{16r_{rms}^3}\\ \displaystyle\frac{dr_{rms}}{dz} = y \end{matrix}. $$ Пусть

$$ \vec X = \begin{bmatrix} y \ r_{rms} \ \end{bmatrix}, $$

теперь составим дифференциальное уравнение \(X' = F(X)\).

Здесь \(k\) — жёсткость соленоида.

$$ k = \left ( \frac{eB_z}{2m_ec\beta\gamma} \right )^2 = \left ( \frac{e B_z}{2\beta\gamma\cdot 0.511\cdot 10^6 e \cdot \mathrm{volt}/c} \right )^2 = \left ( \frac{cB_z[\mathrm{T}]}{2\beta\gamma\cdot 0.511\cdot 10^6 \cdot \mathrm{volt}} \right )^2. $$

Для расчёта потребуется скорость света, входящая в выражение для жёсткости соленоида.

c = 299792458 # (m/s) скорость света

Зададим параметры пучка и начальные условия на входе. Ток составляет 2 кА при альфвеновском токе 17 кА, нормализованный эмиттанс — 180 мм·мрад; на входе в тракт среднеквадратичный радиус равен 27.5 мм, а его производная — 35.4 мрад. Закомментированные значения рядом представляют собой предыдущий вариант: эти четыре числа варьируются при настройке.

I = 2000#2000.000 # A Ток в ускорителе    
Ia = 17000 # A Альфвеновский ток
emitt_n = 1.80e-4#1.80e-5 # m нормализованный эмиттанс
R_rms = 27.5e-3 #m r_rms пучка
dR_rmsdz = 35.4e-3 #50e-3 dr_rms/dz пучка ## 50/2*sqrt(2) = 35.4

Правая часть системы реализована функцией dXdz, которая получает вектор состояния X = [dr/dz, r] и координату z, а возвращает его производную. Внутри собираются все члены уравнения огибающей. Здесь P обозначает первеанс, то есть вклад пространственного заряда; emitt — геометрический эмиттанс, полученный из нормализованного делением на \(\beta\gamma\); K — жёсткость соленоида. Первые два члена расширяют пучок, третий его сжимает, а последние два описывают адиабатическое затухание при ускорении и потому требуют не только \(\gamma\), но и её производных, причём вторая вычисляется численно.

from scipy.integrate import odeint

def derivative(f, x, dx):
    return (f(x + dx) - f(x - dx)) / (2 * dx)

def dXdz(X, z):
    
    y     = X[0]
    sigma = X[1]
    
    g = gamma(z)
    
    dgdz = -Ez(z)/mc
    d2gdz2 = -derivative(lambda z:Ez(z),z, dz)/mc
    beta = np.sqrt(1-1/(g*g))
    K = ( c * Bz(z) / (2*0.511e6))**2
    emitt = emitt_n/(g*beta) # m
    P=2*I/(Ia*beta*beta*beta*g*g*g)
    dydz = P/(4*sigma) + emitt*emitt/(sigma*sigma*sigma*16) - K*sigma/(g*beta)**2 - \
           dgdz*y/(beta*beta*g) - d2gdz2*sigma/(2*beta*beta*g)
    dsigmadz = y
    
    return [dydz, dsigmadz]

Обёртка над odeint подставляет начальные условия и возвращает только среднеквадратичный радиус. Точность ослаблена до rtol = 0.001 вместо стандартных \(1.5\cdot10^{-8}\). Огибающая является гладкой функцией, лишние знаки в ней неразличимы, а расчёт заметно ускоряется; именно за счёт этого библиотеку KENV [45] удаётся вызывать тысячи раз в цикле оптимизации.

def find_sigma(z):
    X0 = [dR_rmsdz, R_rms]      
    
    sol = odeint(dXdz, X0, z, rtol = 0.001) # rtol - относительная точность (по умолчанию rtol = 1.49012e-8)
    sigma = sol[:, 1] # m

    return sigma # m

Настроим оформление кривой пучка, запустим решатель и построим огибающую.

%opts Area.Beam [aspect=3 show_grid=True] (color='red' linewidth=1 alpha=0.3)

sigma = find_sigma(z)

dim_sig = hv.Dimension('r_{rms}', label="$r_{rms}$", unit='mm', range=(0,150))

sigma_img = hv.Area((z,sigma*1e3), kdims=[dim_z], vdims=[dim_sig], group='Beam')
sigma_img

На рисунке приведён среднеквадратичный радиус пучка вдоль пятнадцатиметрового тракта: по горизонтали отложена координата в метрах, по вертикали — радиус в миллиметрах. Существенны два обстоятельства. Во-первых, положение максимума огибающей: в этом сечении пучок ближе всего к стенке вакуумной камеры, и если максимум превышает апертуру, часть пучка будет потеряна. Во-вторых, плавность кривой между соленоидами: сильные перетяжки означают, что фокусировка избыточна.

Расчёт объёмом в несколько десятков строк выполняется за доли секунды, поэтому уравнение огибающей остаётся рабочим инструментом наряду с PIC-кодами: там, где необходимо перебрать тысячи вариантов настройки соленоидов, ансамбль макрочастиц не может быть рассчитан требуемое число раз. Практическое применение рассматриваемого подхода разбирается в следующей главе.