Полуаналитический метод решения задачи об ансамбле ридберговских атомов Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 19, № 4, 2014

Полуаналитический метод решения задачи об ансамбле ридберговских атомов*

В. П. Жуков1'2, И. И. Бетеров3, М. П. ФЕДОРУК1'4, И. И. Рябцев3 1 Институт вычислительных технологий СО РАН, Новосибирск, Россия

2Новосибирский государственный технический университет, Россия 3Институт физики полупроводников СО РАН, Новосибирск, 4Новосибирский государственный университет, Россия e-mail: zukov@ict.nsc.ru, beterov@isp.nsc.ru, ryabtsev@isp.nsc.ru

Жуков В.П., Бетеров И.И., Федорук М.П., Рябцев И.И. Полуаналитический метод решения задачи об ансамбле ридберговских атомов // Вычислительные технологии. 2014. Т. 19, № 4. С. 32-41.

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

Ключевые слова: симплектический метод, решение обыкновенных дифференциальных уравнений, ридберговские атомы.

Zhukov V.P., Beterov I.I., Fedoruk M.P., Ryabtsev I.I. Semi-analytical method for ensemble of Rydberg atoms // Computational Technologies. 2014. Vol. 19, No. 4. P. 32-41.

The problem for determination of wave function for an ensemble of Rydberg atoms requires solution of a set of ordinaries differential equations. The number of equations increases exponentially with the number of atoms. The solution of this problem is of oscillatory type character. Calculation of the oscillation's phase and normalization of the wave function with a high degree of accuracy is very difficult for times which are much greater than the oscillations period. In the article, we suggest a very effective semi-analytic splitting method for finding of the wave function for the ensemble of Rydberg atoms. It is based on analytical (exact) solution of the problem for two states (two connected oscillators). The method gives very high accuracy for determination of the wave function (including the phases of the function) and exactly conserves particles number. It takes several hours for PC to calculate the ensemble of ten atoms with 3 levels. Earlier ensemble of only 7 atoms was calculated using Adam's type methods. The suggested semi-analytical method can be implemented to other similar problems.

Keywords: symplectic algorithm, solution of the ordinary differential equations, Ry-dberg atoms.

* Работа выполнена при финансовой поддержке Интеграционного проекта № 62 и гранта Президента РФ для государственной поддержки ведущей научной школы РФ (школа акад. Ю.И. Шокина, № НШ-5006.2014.9).

Введение

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

Многочастичная волновая функция ансамбля атомов в условиях дипольной блокады впервые была рассмотрена в работе [2] при численном моделировании вероятности одновременного возбуждения двух ридберговских атомов в ансамбле в зависимости от расстояния между ними. Моделирование проводилось для 14 двухуровневых атомов. В [3] рассматривались коллективные осцилляции населённостей в ансамбле взаимодействующих атомов, содержащем 100 двухуровневых атомов. Для решения этой задачи использовался целый ряд приближений. Так, ансамбль был разбит случайным образом на 24 ансамбля меньшего размера, а в отношении последних сделано предположение о наличии режима полной блокады, в котором взаимодействие с лазерным излучением описывается аналитически. Кроме того, было применено ограничение на число возможных возбуждений в совокупности ансамблей. Тем не менее размерность системы составляла более миллиона уравнений.

В работе [4] рассматривались адиабатические процессы при возбуждении ансамбля ультрахолодных атомов. Показано, что при возбуждении чирпованным импульсом вероятность возбуждения одного ридберговского атома не зависит от числа атомов в ансамбле. Численное моделирование с учётом резонансного диполь-дипольного взаимодействия было ограничено семью четырёхуровневыми атомами.

Основная трудность при решении задачи о возбуждении ансамбля ультрахолодных атомов состоит в интегрировании системы большого числа уравнений на интервале времени, намного превышающем характерный период колебаний. При этом необходимо с высокой точностью передать фазу колебаний и обеспечить сохранение нормировки волновой функции (в других задачах — сохранение энергии). К данному классу относятся также задачи о движении планет [5-8] и движении вихревых нитей [9, 10]. Методы типа метода Рунге — Кутты, Адамса и т.п., дающие хороший результат на временах порядка нескольких периодов, для таких задач оказываются неэффективными. В этом случае применяются методы, сохраняющие фазовый объём (в том числе симплекти-ческие).

Предлагаемый в данной работе сохраняющий фазовый объём метод был использован для численного моделирования вынужденного рамановского адиабатического возбуждения в ансамбле, состоящем из 10 трёхуровневых атомов [11]. Для ансамблей из семи атомов было получено согласие с численным расчётом, проведённым с использованием системы Wolfram Mathematica 7. Кроме того, для взаимодействия с непрерывным лазерным излучением в режиме полной дипольной блокады авторы наблюдали возникновение осцилляций населённостей с частотой, пропорциональной числу взаимодействующих атомов, что согласуется с аналитическим решением [11]. При этом представленный метод оказался в десятки раз быстрее применяемых ранее (варианты методов Рунге — Кутты и Адамса). В частности, решение приведённой в настоящей работе задачи об ансамбле из 10 атомов требует нескольких часов работы на персональном компьютере.

1. Математическая постановка задачи

Ансамбль из Nat0 атомов, каждый из которых содержит Niv0 (в данной работе Niv0 = 3) энергетических уровней, описывается волновой функцией, имеющей k0 = (Niv0)Nat0 компонент Ck, представляемых уравнением [11]

llCt = (N2)V> - NÍ1)6osgn(t - t*)) Ck + ^ Vm(k,t)e-^m(k)c(nm(k)). (1)

m=1,M0(k)

Здесь член с N/(1) описывает взаимодействие с внешним постоянным полем, которое меняет знак в момент времени t*, число N^ равно числу атомов в состоянии k, находящихся на втором уровне (на первом возбуждённом энергетическом уровне), 50 —

отстройка частоты лазерного излучения от резонанса с атомными уровнями. (2)

Число N2) = 1, если в состоянии k имеются два и более атомов, находящихся на третьем (высшем) энергетическом уровне. В противном случае N^2) = 0. Этот член отвечает за коллективные взаимодействия высоковозбуждённых атомов. Представляет интерес случай, когда состояния с двумя и более атомами, находящимися на третьем уровне, могут быть исключены из рассмотрения (полная ридберговская блокада [11]). Численный код предусматривает оба варианта. В настоящей работе все примеры расчётов используют последний вариант, соответственно член с N(2 в (1) можно опустить.

Член с Vm(k,í) отвечает за переходы из состояния k в состояние nm(k). Здесь рассматриваются четыре возможных перехода системы на соседние по уровням энергии состояния. Первый — (только) один из атомов переходит с уровня 1 (основной) на уровень 2; второй — переход, обратный первому; третий — (только) один из атомов переходит с уровня 2 на уровень 3; четвёртый — переход, обратный третьему. В случае переходов 1 и 2 (рис. 1) имеем

Vm(nm(k))) = (e-(t-tci)2/á2 + e-(í-íc2)2/¿2) , (2)

в случае переходов 3, 4 —

Vm(nm(k))) = ^ (e-(í-íSl)2/¿2 + e-(t-tS2)/s) . (3)

Рис. 1. Зависимость коэффициентов Vm(nm(k))), соответствующих (2) и (3), от времени

В представляемых ниже расчётах фазы <^т(к) = 0.

В начальный момент времени все атомы находятся на основном уровне (состояние к = 1), т.е. = 0) = 1, си= = 0) = 0. При таких начальных условиях ком-

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

2. Параметризация

Поскольку число уравнений к0 = 3^4° велико, то возникает вопрос о задании коэффициентов поставленной выше задачи. Для этой цели в коде вводится массив Щу (т,к). Элемент массива равен уровню, на котором находится т-й атом в к-м состоянии. Состояние будем нумеровать следующим образом:

к =1 + (т, к) - 1)№оГ-1.

т=1

Такая нумерация позволяет определить п1у (т, к) по к. По существу слагаемое п1у (т, к) является представлением целого числа к в Л^0-м исчислении и может быть определено с помощью соответствующей процедуры.

Чтобы при реализации случая, когда состояния с двумя и более атомами, находящимися на третьем уровне, исключаются из рассмотрения, код подвергался минимальным изменениям, удобно ввести целый массив М33(к): М33(к) = 0, если состояние к содержит два и более атомов на третьем уровне, М33(к) = 1 — в противном случае.

Удобно ввести массив, описывающий переходы (¿¿г,]). Здесь ¿¿г — номер перехода (цг = 1, 2, 3,4 для переходов с первого уровня на второй, со второго на первый, со второго на третий, с третьего на второй соответственно), (цг,] = 1) — номер уровня, с которого атом уходит, (цг,] = 2) — номер уровня, на который атом приходит. Например, Ы^г (1,1) = 1, (1, 2) = 2. Поскольку число возможных переходов невелико, заполнить этот массив "вручную" несложно.

Массив ^г(Цг, Ц) позволяет легко заполнить массивы пт(к). Число пт(к) равно номеру состояния, с которым связано состояние к (т — номер ненулевого недиагонального элемента гамильтониана), М0(к) — полное число таких связей (число всех переходов из состояния к: т = 1,..., М0(к)). Пользуясь (цг,Ь), можно заполнить массив Щ(т,к), служащий для удобства программирования временных зависимостей (2), (3) при вычислении коэффициентов Ут(к) в (1): ^(т,к) = 1 или 2 в случае (2) или (3) соответственно.

Уравнения (1) устроены так, что если недиагональный коэффициент (к, к') не равен нулю, то элемент (к', к) также не равен нулю. Поэтому по данному признаку удобно разбить все компоненты по парам и ввести массивы краг (^р), траг (^р). Здесь N —

номер пары (к, к'), краг — номер компоненты волновой функции (меньший из к, к1), траг — номер ненулевого элемента в строке к: пт=траг(и) (к) = к'. Параметризация остальных величин не представляет труда.

3. Метод решения

Для решения (1) перейдем к новым переменным

а- = с-е \ °

в которых гамильтониан Н (оператор в правой части (1)) не содержит диагональных элементов

г-О- = ^ Ут(к,г)а(пт(к))ехр(-р^т(к) - гДт(к)),

т=1,М°(к)

г

Дт(к) = (М^(к) - М(2))Уо1 - (Л£)(к) - ^0 У sgn(í - 1*)<И. (4)

0

Система является эрмитовой: для любого пт(к) существует т* такое, что

пт* (пт* (к)) = Ут* (пт* (к)) = Ут(к) ,

фт* (пт* (к)) = -фт(к), Дт* (пт* (к)) = -Дт(к). Соответственно сохраняется

Е |akI2 = 1. (5)

к=1,к°

Для решения (4) введём сетку по времени гп = пт, ¿п+1/2 = (п + 1/2)т. Шаг по времени т будем подбирать таким образом, чтобы г* совпадало с одним из гп. На каждом интервале (¿п+1 : гп) функции (2), (3) будем полагать постоянными и равными их значению в момент времени ¿п+1/2.

Система (4) описывает парные связи состояний к О пт(к). В случае только одной пары функций а1 и а2 на интервале (¿п+1 : гп) необходимо решать систему типа

,да1

дг ,0а2

Уе-га-гв(г-п)а2, (6)

г^ = Уега+гв(г-гп)а1 (7)

с постоянными коэффициентами У, а, в и начальными данными а1 = аП, а2 = аП при г = гп. Решение имеет вид

и е-гу+Т + 7/ ег^-т Уе-га

аг(1 = гп+1) = —--ап, + —-(е-г"+т - ег"-т) а^

и ' + V- 1 и+ + и-К ' 2'

е-гу-т — егу+т и,е-ги-т + V е™+Т

а2 (г = гп+1) = е---ап + ^-—-ап,

К и+ + Уе-га 1 v+ + 2

^ = у/У2 + в2/4 ± в/2.

Предлагаемый метод заключается в последовательном аналитическом решении уравнений (6), (7) для каждой пары без учёта действия других пар. При этом для первой пары в качестве начальных условий используются значения функций в момент времени ¿га, а для последующих — значения функций, рассчитанные на предыдущих этапах. Значения , полученные после перебора всех пар, считаются величинами (¿га+1).

Описанный метод является несимметричным относительно порядка, в котором берутся пары. Это приводит к большой ошибке. Для уменьшения ошибки будем менять порядок взятия пар на интервалах (¿га+1 : ¿га) и (¿га+1 : ¿га+2) на противоположный. Это легко реализовать с помощью массивов краг, траг.

Представленный метод сохраняет фазовый объём. Действительно, на каждом этапе (при вычислении изменений в каждой паре) имеем

д<+1 За^1 д<+1 За^1

даП даП даП даП

Очевидно, что на каждом этапе выполняется (5). При этом с точностью до ошибок округления имеем сохранение фазового объёма и выполнение (5) для полного шага по времени. Также очевидна обратимость во времени (с учётом соответствующей последовательности взятия пар).

Описываемый метод является методом расщепления. Оператор Шредингера H представляется здесь в виде суммы вообще говоря некоммутирующих операторов Шредин-гера Hj для каждой из пар и точное решение

tn+2

-i Г H (t)dt

ак (Г+2) = e ак (tn)

аппроксимируется выражением

tn + 2 tn + 2 tn+1 tn + 1

-i f Hi(tn+3/2)dt -i f Hjp(tn+3/2)dt -i f Hjp(t"+1/2)dt -i f Hi(tn+1/2)dt

e tn+1 ...e tn+1 e tn ...e tn ak (tn).

Подобная схема широко используется в работах [6, 12] и имеет второй порядок аппроксимации.

4. Результаты расчётов

Изучим влияние шага сетки т на точность предлагаемого метода на примере ансамбля

из 10 атомов при типичных для подобных задач параметрах [11], а именно: 80 = 400п,

^ = 9, 8 = //2, По = 30п, ¿С1 = ^ - 3.5, ^2 = 2^ - ^1, ^ = 40п, tsl = ^ - 5.5,

ts2 = 2^—tsl, N<^0 = 10, NI= 3. Здесь величины времени измеряются в мкс, частоты — 1

в мкс 1.

Будем рассматривать влияние шага сетки т на зависимость от времени амплитуды а и фазы ^ компоненты волновой функции, соответствующей основному состоянию я*;=1 = аег^, и действительной части компоненты волновой функции, соответствующей первому возбуждённому состоянию И,е(ае). Основное состояние, в котором все атомы находятся на нижнем (первом в нашей нумерации) уровне, с точки зрения физических приложений представляет главный интерес. Первое возбуждённое состояние, в котором один из атомов находится на втором уровне, а все остальные занимают нижний уровень, интересно в связи с тем, что соответствующая ему компонента волновой

функции колеблется с наименьшим для этой задачи периодом Tmin = 2п/80 = 0.005 мкс. При используемых начальных условиях все Nat0 компонент волновой функции ak, соответствующие первому возбуждённому состоянию, одинаковы. При численном решении между этими ak возникают различия, но они более чем на порядок меньше погрешности, связанной с конечным т.

На рис. 2 представлена зависимость амплитуды a компоненты волновой функции от времени при разных т. На графике, соответствующем всему интервалу времени, на котором ищется решение, влияние т становится трудноразличимым уже при т < 0.001 = Tmin/5. Поэтому рядом приведены более подробные графики на меньшем временном интервале. При т < 0.0002 = Tm¡n/25 различий в поведении амплитуды основного состояния нет даже в деталях. Минимальный шаг сетки т, с которым проводились расчёты, составляет 10-5.

Аналогична ситуация с фазой основного состояния ^ (рис. 3), которая неплохо передается уже при т = 0.0004. Исключение составляет диапазон времен примерно (7.5 : 10.5), в котором амплитуда основного состояния близка к 0 и соответственно фаза состояния не имеет физического смысла.

На рис. 4 показана зависимость Re(ae) от времени. Точность вычисления этой функции хотя и ниже, чем величин a и период колебаний которых намного больше, но достаточно высока уже при т < 0.0004 = Tmin/12.5. При этом даже при шаге сетки т = 0.001 мкс, соответствующем трём точкам на период, имеется количественное соответствие с точным решением. Очевидно, что несимплектические методы обеспечить такую точность при указанных шагах не в состоянии. Интересным является то обстоятельство, что в данном методе не происходит накопления фазы на больших временах. Продемонстрируем это на следующем примере. Пронумеруем максимумы функции Re(ae) в порядке увеличения t. Обозначим через tmax(Nm, т) момент времени, в который достигается Nm-й максимум сеточной функции Re^^í"^)), полученной при использовании шага т. Отличие ) от точного значения этой функции (в качестве

а б

Рис. 2. Зависимость амплитуды а компоненты волновой функции, соответствующей основному состоянию, от времени: а — на всем временном интервале, б — на малом временном интервале (сплошная линия — т < 0.0002 мкс, штрихпунктирная — т = 0.0004 мкс, штриховая — т = 0.001 мкс)

которого можно взять значение, полученное при т = 10-5) на всем интервале решения задачи характеризует накопление фазы, связанное с погрешностью метода. В качестве ^ах сеточной функции возьмём если Ке(ае(^-1,т)) < И,е(ае(^,т)) < Ке(ае(^+1,т)). Такой подход определяет положение максимума с погрешностью, примерно равной т.

На рис. 5, а изображена зависимость Д^ах(т) = ^гаах^т,т) — ^ах^т,т = 10 5)| от Nm. Видно, что за исключением двух точек Д^ах не превышает величины ~ т, причём величина ошибки не растёт. Упомянутые точки, в которых наблюдается резкий рост Д^ах, соответствуют моментам времени, когда модуль комплексной амплитуды ае (а не только её действительная часть) обращается в нуль. В эти моменты понятие фа-

1

ч>

0.5

-0.5

10

t, мкс

15

1

0.5

-0.5

: 1;< 1 1 [Ч\ !\

Г' 11 \\ 1 ^

- и 1 '' V1 4

Vй 1 ! V \

\'у| 1 \ч- ч

и > V

7 Ц ' в 1

-1 V 1 V 1 '

- \ и 1 1 1 X"-.. ^

~ 1 - 1 |\1 1 1 1\ 1 \ * 1 ^

1 1 IV „ 1 > 1 1

I Г- \ | и ' I V1 1 1

Т * | и ^ V ч 1

"V- 1 и ^ М V ^ 1

I у | V ч 1

и 1 V 1 : \ 1 1

I ь ¡1 1 V I1 1 \\ \ 1 1

- \\ *1 V 1 1 и |' I \\ 1, \\ \ 1

~ и Ч \> | 11 V \ 1

: V ч \ 'У \ 1

- , 1 - ■ »■ ■ А 1. *.....\ ч' ........

13

13.5

1 мкс

14

Рис. 3. Зависимость фазы р компоненты волновой функции, соответствующей основному состоянию: а — на всем промежутке времени (сплошная линия — т < 0.0001, штриховая т = 0.0003), б — на малом промежутке времени (сплошная линия - т < 0.0002 мкс, штрих-пунктирная — т = 0.0004 мкс, штриховая — т = 0.001 мкс)

Ие(ае) 0.03

15 * мкс 20

0.01

Ие(ос)

0.005

-0.005

-0.01

13.99

5

13.995

_1_

t, мкс

14

Рис. 4. Зависимость Ке(ае) от времени: а — за всё время расчёта, б — на малом промежутке времени (сплошная линия т < 0.0002, квадраты — т = 0.0004, ромбы — т = 0.001 мкс)

а

а

0.02

0.01

■-j . i

0.02 г

0.01 -

1000 2000

3000

X,

0.0004

Rc(ac)

0.0002

-0.0002 г

-0.0004

2

t мкс

3.06

а

Рис. 5. Зависимость Д£тах(т) от при т = 0.001 (внизу) и т = 0.0004 (вверху) (а) и зависимость Кв(ае) от £ при т = 10-5 (сплошная линия) и т = 0.0004 (штриховая линия) в окрестности первой точки роста Д£тах(т) (б)

зы теряет физический смысл. На рис. 5, б представлена зависимость И,е(ае (¿)) при разных т в окрестности первой точки роста. Виден сбой фазы И,е(ае(£)) в непосредственной окрестности особой точки, но перед этой окрестностью и после неё наблюдается хорошее совпадение фаз.

Таким образом, можно сделать вывод о том, что предложенная схема позволяет получать достаточно точное решение (в том числе фазы колебаний) при т < Ттщ/10 на временах, на много порядков превышающих наименьший период задачи. Закон сохранения (5) выполняется в расчётах с точностью до ошибок округления (< 10-10).

Список литературы

[1] Lukin M.D., Fleischhauer M., Cote R. Dipole blockade and quantum information processing in mesoscopic atomic ensembles // Phys. Rev. Lett. 2001. Vol. 87. Paper 037901 (4 p.).

[2] Robicheaux F., Hernandez J.V. Many-body wave function in a dipole blockade configuration // Phys. Rev. A. 2005. Vol. 72. Paper 063403 (4 p.).

[3] Stanojevic J., Cote R. Many-body Rabi oscillations of Rydberg excitation in small mesoscopic samples // Ibid. 2009. Vol. 80. Paper 033418 (9 p.).

[4] Beterov I.I., Tretyakov D.B., Entin V.M. et al. Deterministic single-atom excitation via adiabatic passage and Rydberg blockade // Ibid. 2011. Vol. 84. Paper 023413 (6 p.).

[5] Чегодаева Е.А. Метод параллельного симплектического интегрирования уравнений движения малых тел Солнечной системы // Вестник ЮУрГУ. Математика, физика, химия. 2006. Вып. 7. С. 150-151.

[6] Чегодаева Е.А. Метод симплектического интегрирования уравнений движения малых тел солнечной системы // Там же. 2005. Вып. 5. С. 49-55.

[7] Emel'yanenko V. An explicit symplectic integrator for cometary orbits // Celestial Mech. and Dynam. Astronomy. 2001. Vol. 74. P. 287-295.

[8] Wisdom J., Holman M. Symplectic maps for the N-body problem // The Astronomical J. 1991. P. 1528-1538.

[9] Kuvshinov B.N., Schep T.J. Double-periodic arrays of vortices // Phys. Fluids. 2000. Vol. 12, No. 12. P. 3283-3284.

[10] Kuvshinov B.N., Schep T.J. Holtsmark distributions in point-vortex systems // Phys. Rev. Lett. 2000. Vol. 84, No. 4. P. 650-653.

[11] Beterov I.I., Saffman M., Yakshina E.A. et al. Quantum gates in mesoscopic atomic ensembles based on adiabatic passage and Rydberg blockade // Phys. Rev. A. 2013. Vol. 88(1). Paper 010303(R) (5 p.).

[12] Muslu G.M., Erbay H.A. High-order split-step Fourier schemes for generalized nonlinear Schrodinger equation // Math. and Appl. Mathematics. 1977. Vol. 57 (1). P. 1-12.

Поступила в 'редакцию 14 февраля 2014 г., с доработки — 23 июня 2014 г.