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

Рассеяние ионов на атомах рассчитывается в предположении простого экранированного потенциала Резерфорда (потенциала Бора). Одна из целей состоит в том, чтобы проверить, что дифференциальное сечение остается конечным для рассеяния вперед, т. е. показать, что сингулярность в сечении (и фактически неинтегрируемость) при рассеянии на чисто кулоновском потенциале обусловлена его дальнодействием, т. е. отсутствием сходимости при больших прицельных параметрах.

• restart; with(plots): Digits:=11:

Вначале определим параметры: выберем единицы Бора, в которых масса электрона = 1, и рассмотрим рассеяние протона на атоме для Z2 = 10 (атом неона). Параметр потенциала Бора определялся из экспериментальных данных по рассеянию атомов неона, a = 0,52 а. е. Масса атома неона ≅ 20 масс протонов.

• M1:=1836: M2:=20*1836: mu:=M1*M2/(M1+M2);

• Z1:=1: Z2:=10: a:=0.52;

• V:=r->Z1*Z2*exp(-r/a)/r;

Зададим процедуру вычисления величины расстояния максимального сближения при заданных начальной скорости соударения и прицельном параметре. Она основана на условии перигелия, в котором потенциальная энергия максимальна вдоль траектории отталкивающего рассеяния. Для заданной скорости соударения на бесконечности (v0) и прицельного параметра b определяем полную относительную энергию рассеяния E на бесконечности (нулевой потенциал) и величину углового импульса L. Расстояние максимального сближения minR получается из решения для условия сохранения энергии (нелинейного уравнения):

• minR:=proc(v0,b) local E,L,peri; global mu,a,Z1,Z2;
• E:=mu*v0^2/2;
• L:=mu*v0*b;
• peri:=E-L^2/(2*mu*r_min^2) - V(r_min);
• fsolve(peri,r_min=0..infinity) end:

Cкорость соударения на бесконечности v в единицах Бора берем = 1, это будет промежуточным значением, соответствующим скорости протона, сравнимой с классической орбитальной скоростью электрона в 1s-состоянии водорода.

Рассматриваем быстрые и медленные столкновения (взаимодействие «меньше» и «больше»):

• v0:=0.5;

• r_clap:=minR(v0,2);

Расстояние максимального сближения получилось немного больше, чем прицельный параметр.

Нужно взять интеграл по расстоянию от величины максимального сближения до бесконечности. Похоже, что Maple может вычислить интеграл численно. Кроме того, можно уменьшить точность вычисления интеграла с помощью Digits.

• theta:=proc(v0,b) local E; global mu; E:=mu*v0^2/2;
• evalf(Pi)-2*evalf(Int(b/(r*sqrt(r^2*(1-V(r)/E)-b^2)),
  r=minR(v0,b)..infinity),Digits-1); end:
• theta(v0,0.1);

Для малого прицельного параметра и малой скорости соударения получился большой угол рассеяния. Теперь построим цикл для прицельного параметра. Надо исследовать его малые и большие значения, потому что нас интересует сравнение со случаем неэкранированного резерфордовского рассеяния.

Для малых прицельных параметров (большие углы отклонения) ожидается, что для обоих потенциалов результаты будут в согласии, так как отклонение максимально при максимальном сближении. Для больших значений b экранированный случай приводит к крошечным углам отклонения, которые становятся незначащими при b → ∞. В чисто кулоновском потенциале всегда есть отклонение, даже когда b становится бесконечным, что неправильно (граничный случай).

• db:=0.01; N:=200:

• PP:=[seq([db*i,theta(v0,db*i)],i=1..N)]:
• PP[200];

• P1:=loglogplot(PP,style=point,color=red):
• P2:=loglogplot(2*arccot(b*mu*v0^2/(Z1*Z2)),b=db..N*db,
  color=blue):
• display(P1,P2,labels=["b","theta"],axes=boxed);

Здесь видно, как численная оценка интеграла согласуется с «резерфордовским» результатом для малых прицельных параметров.

Интересно отметить, что усечение интеграла при расчете θ до конечного верхнего предела (вместо бесконечности) приводит к серьезным ошибкам для промежуточных и больших прицельных параметров. Это означает, что при численной оценке интеграла важно отображать весь диапазон интегрирования даже для короткодействующего потенциала рассеяния. Это несколько неожиданно, особенно при исследовании соотношения между θ и b путем численного решения ДУ, для чего обычно хватает конечного диапазона интегрирования.

Изучите взаимосвязь между прицельным параметром и углом рассеяния для различных скоростей соударения, сохраняя при этом все остальные параметры фиксированными. Влияет ли значение b на скорость соударения?

Теперь надо посмотреть на поведение дифференциального сечения при малых углах. Что происходит для экранированного потенциала с поперечным сечением рассеяния под небольшими углами?

• PP[1];

Для расчета db/dthetaсначала вычислим обратную величину dtheta/db и возьмем обратную ей, а затем вычислим dsigma/dOmega, используя формулу центральных конечных разностей на равномерной сетке значений b:

• for i from 2 to N-1 do:
  dsdO[i]:=db*i/abs((PP[i+1][2]-PP[i-1][2])/(2*db))
  /sin(PP[i][2]); od:

Чтобы правильно его построить как функцию полярного угла рассеяния (рис. 4.18.1–4.18.4), диапазон b задан как параметр:

• PPc:=[seq([PP[i][2],dsdO[i]],i=2..N-1)]:
• P3:=loglogplot(PPc,style=point,color=red):
• P4:=loglogplot(1/4*(Z1*Z2/(mu*v0^2))^2*csc(theta/2)^4,
  theta=PP[1][2]..PP[N-50][2],color=blue):
• display(P4,P3,labels=["theta","ds/dO"],axes=boxed,
  title="Bohr potential(red), Coulomb potential (blue)
  differential cross section: dsigma/dOmega");

Рис. 4.18.1

Видим, что для малых параметров численная оценка интеграла согласуется с «резерфордовским» результатом.

Рис. 4.18.2

Для диапазона показанных на рис. 4.18.2 значений θ резерфордовское сечение следует простому закону 1/θ4, так как мы находимся в режиме sin(θ) = θ. Боровское сечение продолжает расти, но не так быстро, как для малых углов. Чтобы выяснить, что на самом деле происходит, нужно рассмотреть гораздо большие значения b, т. е. гораздо меньшие полярные углы рассеяния.

Повторяем цикл для расчета с потенциалом Бора:

• db:=0.04; N:=200:

• PP:=[seq([db*i,theta(v0,db*i)],i=1..N)]:
• loglogplot(PP,style=point,color=red,labels=["b","theta"],
  axes=boxed,title="Bohr potential: theta(b)");

Рис. 4.18.3

По-видимому, есть проблема с шумом для больших b, и, вероятно, следует увеличить число значащих цифр с помощью Digits. Если это сделать, то окажется, что вычисление угла рассеяния «тормозит» из-за того, что evalf (Int) не может вернуть значение (процедура не достигает желаемой точности). Чтобы «увидеть» проблему интегрируемости, на рис. 4.18.4 показана зависимость dsigma/dtheta:

• for i from 2 to N-1 do:
  dsdO[i]:=db*i/abs((PP[i+1][2]-PP[i-1][2])/(2*db))/sin(PP[i][2]); od:
• PPc:=[seq([PP[i][2],6.283*dsdO[i]*sin(PP[i][2])],i=2..N-1)]:
• loglogplot(PPc,style=point,color=red,
  labels=["theta", "ds/dtheta"],
  axes=boxed,title="Bohr potential differential cross section");

Рис. 4.18.4

Следует учитывать, что при очень малых углах отклонения поперечное сечение нельзя рассчитать путем инвертирования производной от θ(b), потому что абсолютная ошибка в интеграле превышает его фактическое значение. Данные ниже демонстрируют провал процедуры:

• seq(dsdO[i],i=N-50..N-1);

Получились очень большие числа, порядка 1013–1015.

По-видимому, результаты для dsigma/dtheta обращаются при наименьших показанных значениях θ, что указывает на интегрируемость дифференциального сечения.

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

Проверьте режим больших углов рассеяния: насколько хорошо численный расчет для потенциала Бора согласуется с аналитическим результатом величины dsigma/dOmega для потенциала Резерфорда?

Исследуйте другой отталкивающий центральный потенциал, обладающий конечным радиусом действия. Исследуйте зависимость b-theta и дифференциальное сечение рассеяния.

Исследуйте потенциалы со степенным спадом, которые затухают быстрее, чем кулоновский потенциал, и изучите прямое дифференциальное сечение рассеяния.

Исследуйте рассеяние в конечной области для потенциала притяжения V(r).