PARALLEL ALGORITHM FOR CALCULATING THE WIGNER FUNCTION FOR A QUANTUM SYSTEM WITH A POLYNOMIAL POTENTIAL

Abstract

The article considers the construction of a parallel algorithm on the GPU computing architecture for finding the Wigner function of a quantum system with a polynomial potential. A numerical-analytical method for constructing the Wigner function, which is based on calculating the trace of the product of the density matrix and the matrix of the Weyl operator, is described. The operators are represented in the basis of a quantum harmonic oscillator, for which the Moyal equation transforms into the Liouville equation. This approach enables to visually analyze the degree of anharmonicity of the system in terms of the off-diagonal elements of the density matrix. The parallel implementation on the GPU massively parallel architecture has reduced the calculation time by two orders of magnitude compared to the single-threaded version on the x86 architecture.

Keywords

• Wigner function
• GPU computing
• quantum system with a polynomial potential

-

1 Объединенный институт ядерных исследований (ОИЯИ), Дубна, Московская область, 141980 Россия 2 Физический факультет МГУ имени М.В. Ломоносова, Москва, 119991 Россия 3 Университет «Дубна», Московская область,141980 Россия 4 Московский технический университет связи и информатики (МТУСИ), Москва, 111024 Россия

E-mail: a [e-mail] В работе рассмотрено построение параллельного алгоритма на вычислительной архитектуре GPU нахождения функции Вигнера квантовой системы с полиномиальным потенциалом. Описан численно-аналитический метод построения функции Вигнера, основанный на вычислении следа произведения матрицы плотности и матрицы оператора Вейля. Представление операторов выполнено в базисе квантового гармонического осциллятора, для которого уравнение Моэля переходит в уравнение Лиувилля. Такой подход позволяет наглядно анализировать степень ангармоничности системы по недиагональным элементам матрицы плотности. Параллельная реализация на массивно-параллельной архитектуре GPU позволила уменьшить время вычисления на два порядка в сравнении с однопоточной версией на x86 архитектуре. Ключевые слова: функция Вигнера, вычисления на GPU, квантовая система с полиномиальным потенциалом Евгений Перепёлкин, Борис Садовников, Наталья Иноземцева, Евгений Бурлаков, Римма Полякова, Павел Сысоев, Марианна Садовникова 1. Введение

Решение современных задач квантовой информатики [1-3], квантовой связи и криптографии [ 4, 5 ], задач обработки сигнала [6-8] содержит аппарат квантовой механики в фазовом пространстве. Функция Вигнера [9] играет важную роль при нахождении средних по фазовому пространству квантовых характеристик системы [ 10-20 ].

Эволюция функции Вигнера W r , p,t  описывается уравнением Моэля W t

1  1l  22l   p,r W  rU , pW    m l1 2l 1!

U r , p  2l1

W. (i.1) где U — потенциал квантовой системы. В случае гармонического осциллятора уравнение (i.1) имеет нулевую правую часть и совпадает с известным классическим уравнением Лиувилля. Указанное упрощение уравнения (i.1) позволяет эффективно получать выражения, описывающие эволюцию функции Вигнера. Для потенциалов степени выше второй решение уравнения (i.1) затруднительно. В работе [ 21 ], используя уникальность потенциала второй степени, были получены явные выражения для матричных элементов wn,k оператора Вейля в базисе гармонического осциллятора ψk  .

W  Tr ˆWˆ  , det det 1   s   W= wn,k  x, p  2  ψn  x  2  ψk  x  s  

 exp  i 2   ps   ds.  wn,k  x, p  1  k e 2x2 2p 2 Pn,k   x  i p 2  , x  i p 

 ,   где  k,n  ck cn матрица плотности;  

. Полиномы Pn,k имеют вид: m

minn,k zns zks Pn,k  z1, z2   2nk n!k !  1 2 , s0 2s s!k  s!n  s! и удовлетворяют условию ортогональности с весом  2  x  ex2 y2 :      2  x, y  Pn1,k1  x, y  Pn2 ,k2  x, y  dxdy  N n12n2 ,k1k2 n1n2 ,k1k2 ,  

minn1,k1 minn2 ,k2  N n12n2 ,k1k2   n1 !n2 !k1 !k2 !   s0 l0 sl,n,keven/odd k  l  s 1 !!

n  l  s 1 !! s!n1  s!k1  s! l !n2  l !k2  l ! . где (i.2) (i.3) (i.4) (i.5) 1, n, k  чётные,  1, n, k  нечётные,  n,k   0, n  чётный, k  нечётный, 0, n  нечётный, k  чётный. При получении выражения (i.3) предполагалось, что волновая функция  разлагается по базису ψk  :

   x,t    cn t  ψn  x .

n0 Отметим, что в этом случае выражение (i.2) может быть записано в виде свёртки: W  CT WC,c (i.6) (i.7)

CT  c1, c2 ,..., C  c1, c2 ,...T . Используя выражения (i.2)-(i.4), в работе [22] был описан метод построения функции N Вигнера для квантовой системы с полиномиальным потенциалом UN  x   anxn степени N . n1 Пусть волновая функция  удовлетворяет уравнению Шрёдингера  2 m2 xx  UN  E   0 и допускает разложение (i.6), тогда векторы коэффициентов C  (i.7) являются собственными векторами, а спектр энергий E= собственными значениями симметричной матрицы Jn,k : JC   E=C ,  0,

Jn,k   Inl,k  l 1 Inl,k2  2 2n Inl11,k  k l1 2 In,k1

N n In0,k   al Inl,k ,

l2 , In0,k   n,k , In1,k 

(i.8) 2n In01,k  2k In0,k1 , где m 2  2a2 ;   m ;  n    n  1  − спектр энергий гармонического осциллятора.

 2  В данной работе описывается параллельный алгоритм, позволяющий эффективно производить вычисление функции Вигнера на конечно-разностной сетке по приведенным выше выражениям. 2. Параллельная реализация на графическом процессоре (GPU) Для квантовой системы, характеризуемой волновыми функциями Вигнера допускают представление: поведения волновых функций  . Как правило, с увеличением номера состояния увеличивается количество осцилляций волновой функции  , что требует дополнительных узлов разностной сетки при решении уравнения Шрёдингера и при численном интегрировании   выражения (1). Объём вычислений обусловлен не только нахождением сеточных функций Wm,n  , функции , но и необходимостью получения средних характеристик квантовой системы, требующих повторного численного интегрирования с функцией Wm,n .

Ускорение вычислений сеточных функций Wm,n и средних характеристик квантовой системы может быть достигнуто тремя способами.

Во-первых, использование явных выражений (i.3), (i.4) для матричных элементов wn,k  x, p позволяет избежать непосредственного численного интегрирования в выражении (1). Функции wn,k  x, p известны в явном виде и для фиксированной разностной сетки ( wnm,k,n ) могут быть вычислены один раз и храниться в памяти.

Во-вторых, благодаря возможности разложения волновой функции  по базису гармонического осциллятора ψ j (i.6):    cj ψ j ,

j0 ψ j  x 

1 21j j!  m 4 em2x2 H j  m x  , j  0  0, (2) (3) можно с одной стороны, вместо количества разбиений сеточной функции  контролировать количество членов NC, ряда (2), а с другой стороны, быстрый экспоненциальный спад (3) гарантирует небольшое значение NC, для приемлемой точности аппроксимации волновой функции  конечным разложением по функциям ψ j (2). Согласно (i.7), (i.8) коэффициенты разложения cj  образуют собственные вектора C  матрицы J , вид которой известен точно (i.8) и может быть вычислен один раз до начала основной процедуры вычислений. Отметим, что матрица J является симметричной, что уменьшает размер памяти для ее хранения и позволяет воспользоваться специальными численными методами при решении задачи на собственные значения.

В-третьих, задача на собственные значения (поиск C  и E=), вычисление значений wnm,k,n и свертки (i.7) Wm,n , а также нахождение квантовых средних по функции Wm,n являются процедурами, допускающими эффективное численное распараллеливание. Подходящей вычислительной архитектурой для таких задач может быть массивно-параллельная архитектура графических процессоров GPU с использованием среды CUDA [ 23, 24 ]. Задача на собственные значения, операции по перемножению матриц имеют реализацию в бесплатных пакетах библиотек компании NVIDIA: cuBLAS, cuSPARS, cuSolver.

В качестве примера рассмотрим алгоритм параллельного вычисления матричных элементов wnm,k,n . Матричные элементы wn,k обладают свойством wn,k  wk,n , поэтому в памяти достаточно хранить только диагональ и верхнюю треугольную часть матрицы W . Размер матрицы W определяется числом NC, , которое удовлетворяет условию NC, Nx или NC,

Np . Количество разбиений N x и Np можно считать величинами одного порядка:  w1,n1,m w1,n2,m... wn,m 

1,NC,  W n,m =  0 w2n,2,m.. wn,m . (4)  3,NC,   ... 0 ...    Заметим, что при желании можно выбрать максимальное значение квантового состояния L  max (для которых будут производиться расчёты) и задать размерность матрицы (4) NC,L  NC,L .

Так как величина N N

x p определяет максимальный объём вычислений, то логично произвести распараллеливание именно по индексам n, m , задав двумерную сетку нитей и блоков nx, np . Заголовок функции-ядра (kernel-function) имеет вид: __global__ void W_Matrix (int n, int k, double *WM...) {int nx = threadIdx.x + blockIdx.x * blockDim.x; int np = threadIdx.y + blockIdx.y * blockDim.y; // (x, p) parallelism ... }

Каждая вычислительная нить (thread) находит свой матричный элемент в соответствующей ее номеру точке фазового пространства (см. рис. 1). Функция-ядро «W_Matrix» запускается порядка 12 NC2,L раз c «host», чтобы произвести вычисления всех элементов матрицы (4). Результаты вычислений можно хранить в глобальной памяти GPU (global memory) с целью их переиспользования другими функциями-ядрами.

Рис. 1. Распределение вычислительной нагрузки по потокам 3. Заключение

Время работы параллельной реализации алгоритма (GPU GTX1050) вычисления функции Вигнера на два порядка меньше времени выполнения однопоточной версии (CPU класса Core i7). 4. Благодарность

Работа выполнена при поддержке гранта РФФИ No. 18-29-10014 и Междисциплинарной научно-образовательной школы МГУ им. М.В. Ломоносова «Фотонные и квантовые технологии. Цифровая медицина» Литература [1] R. P. Rundle, Todd Tilma, J. H. Samson, V. M. Dwyer, R. F. Bishop, and M. J. Everitt, General approach to quantum mechanics as a statistical theory, Phys. Rev. A 99, 012115 – Published 16 January 2019 [2] Ievgen I.Arkhipov, Artur Barasiński, Jiří Svozilík, Negativity volume of the generalized Wigner function as an entanglement witness for hybrid bipartite states, Sci Rep 8, 16955 (2018) [4] Alberto Casado, Santiago Guerra, José Plácido, From Stochastic Optics to the Wigner Formalism: The Role of the Vacuum Field in Optical Quantum Communication Experiments, Atoms 2019, 7, 76

References

Е.Е. Перепёлкин1 , 2 ,4, Б.И. Садовников2, Н . Г. Иноземцева3 , 4, Е . В. Бурлаков2 ,4, Р.В. Полякова1, a, П.Н. Сысоев2, М .Б. Садовникова2 [3] Andersen , U. , Neergaard-Nielsen , J., van Loock, P. et al. Hybrid discrete- and continuousvariable quantum information . Nature Phys 11 , 713 - 719 ( 2015 ).

[5] Casado , A. ; Guerra , S. ; Plácido , J. Wigner representation for experiments on quantum cryptography using two-photon polarization entanglement produced in parametric down-conversion . J.

Phys. B At . Mol. Opt. Phys . 2008 , 41 , 045501 .

Cohen L. Time-Frequency Analysis . Prentice Hall, Englewood Cliffs, 1995 [7] Zayed , A. A New Perspective on the Two-Dimensional Fractional Fourier Transform and Its Relationship with the Wigner Distribution . J Fourier Anal Appl 25 , 460 - 487 ( 2019 ) [8] Claasen , T.A.C.M. , Mecklenbräuker , W.F.G. : TheWigner distribution-a tool for timefrequency signal analysis. II: discrete-time signals, part 2 . Philips J. Res . 35 , 276 - 300 ( 1980 ) [9] E.P. Wigner , On the quantum correction for thermodynamic equilibrium , Phys. Rev . 40 ( June 1932 ) 749 - 759 [10] Meystre P. , Walls D.F. Nonclassical Effects in Quantum Optics . American Institute of Physics, New York, 1991 [11] Kimble H.J. , Walls D .F. // J. Opt. Soc. Am. B . 1987 . V. 4A0 ) [12] Knight P. , Loudon R . // J. Mod. Opt. 1984 . V. 34F ) [13] Giacobino E. , Fabre С // Appl. Phys. B. 1992 . V. 55C ) [14] Wu L.A. , Xiao М. , Kimble H.J. Squeezed states of light from an optical parametric oscillator // J. Opt. Soc. Am. B . 1987 . V. 4. P. 1465 - 1475 .

[15] Schleich W.P. , Wheeler J.A. Oscillations in the Photon Distribution of Squeezed States and Interference in Phase Space // Nature (London) 326 , 574 - 577 A987a) [16] Smithey D.T. , Beck M. , Raymer M.G. , Faridani A. Measurement of the Wigner distribution and the density matrix of a light mode using optical homodyne tomography: application to squeezed states and the vacuum // Phys . Rev. Lett . 1993 . V. 70. P. 1244 - 1247 .

[17] Breitenbach G. , Schiller S. , Mlynek J . Measurement of the quantum states of squeezed light / / Nature. 1997 . V. 387. P. 471 - 475 .

[18] Royer A . Measurement of the Wigner function // Phys. Rev. Lett . 1985 . V. 55. P. 2745 - 2748 .

[19] Vogel K. , Risken H . Determination of quasiprobability distributions in terms of probability distributions for the rotated quadrature phase // Phys . Rev. A. 1989 . V. 40. P. 2847 - 2849 .

[20] Radon J. fiber die Bestimmung von Funktionen durch ihre Integralwerte langs gewisser Mannigfaltigkeiten // Ber. Verh. Sachs. Akad. Wiss. Leipzig, Math.-Nat. Kl. 1917 . V.69. P. 262 - 277 .

[21] Perepelkin E.E. , Sadovnikov B.I. , Inozemtseva N.G. , Burlakov E.V. , Explicit form for the kernel operator matrix elements in eigenfunction basis of harmonic oscillator , Journal of Statistical Mechanics: Theory and Experiment , ( 2020 ) № 023109 [22] Perepelkin E.E. , Sadovnikov B.I. , Inozemtseva N.G. , Burlakov E.V. , Wigner function of a quantum system with polynomial potential , Journal of Statistical Mechanics: Theory and Experiment , ( 2020 ) № 053105 [23] Н. Г . Иноземцева, Е.Е. Перепёлкин, Б.И. Садовников, Оптимизация алгоритмов задач математической физики для графических процессоров, Физический факультет МГУ имени М . В. Ломоносова Москва , 2012 , ISBN 978-5-8279-0107-5 , 256 с.

[24] Е. Е . Перепёлкин, Б.И. Садовников, Н.Г. Иноземцева, Вычисления на графических процессорах (GPU) в задачах математической и теоретической физики, Серия «Классический учебник МГУ» , URSS Москва , 2019 , ISBN 978-5-9710-6490-9 , 240 с.