Computer studies of a dependence of equilibrium state structure on a number of particles for a two-dimensional system of charged particles confined in a disk potential Текст научной статьи по специальности «Физика»

Discrete & Continuous Models

& Applied Computational Science__2024, 32 (1) 99-105

ISSN 2658-7149 (Online), 2658-4670 (Print) http://journals.rudn.ru/miph

Research article

UDC 538.911, 538.915, 519.6

PACS 07.05.Tp, 02.70.Ns, 71.90.+q, 21.60.Jz, 24.60.Lz, 27.80.+W

DOI: 10.22363/2658-4670-2024-32-1-99-105 EDN: BXKOPH

Computer studies of a dependence of equilibrium state structure on a number of particles for a two-dimensional system of charged particles confined in a disk potential

Eduard G. Nikonov1,2'3, Rashid G. Nazmitdinov1,2, Pavel I. Glukhovtsev2

1 Joint Institute for Nuclear Research, 6 Joliot-Curie St, 141980, Dubna, Russian Federation

2 Dubna State University, 19 Universitetskaya St, Dubna, 141980, Russian Federation

3 HSE University, 34 Tallinskaya St, Moscow, 123458, Russian Federation

(received: February 12, 2024;revised: March 10, 2024;accepted: March 24, 2024)

Abstract. The problem of finding equilibrium configurations of one-component charged particles, induced by external electrostatic fields in planar systems, is a subject of active studies in fundamental as well in experimental investigations. In this paper the results of numerical analysis of the equilibrium configurations of charged particles (electrons), confined in a circular region by an infinite external potential at its boundary are presented. Equilibrium configurations with minimal energy are searched by means of special calculation scheme. This computational scheme consists of the following steps. First, the configuration of the system with the energy as close as possible to the expected energy value in the ground equilibrium state is found using a model of stable configurations. Next, classical Newtonian molecular dynamics is used using viscous friction to bring the system into equilibrium with a minimum energy. With a sufficient number of runs, we obtain a stable configuration with an energy value as close as possible to the global minimum energy value for the ground stable state for a given number of particles. Our results demonstrate a significant efficiency of using the method of classical molecular dynamics (MD) when using the interpolation formulas in comparison with algorithms based on Monte Carlo methods and global optimization. This approach makes it possible to significantly increase the speed at which an equilibrium configuration is reached for an arbitrarily chosen number of particles compared to the Metropolis annealing simulation algorithm and other algorithms based on global optimization methods.

Key words and phrases: Thomson atom, Wigner crystal, molecular dynamics

1. Introduction

The question of how charged particles arrange themselves in a restricted planar geometry attracted continuous attention for many decades (for a review see [1]). Modern technology allows us to study such phenomena on the same scale, from Bose condensates with some thousand atoms to quantum dots with a few electrons, providing rich information about specific features of correlation effects in mesoscopic systems (see, for example, [2, 3]). However, finding the exact analytical equilibrium charge distribution (the one that makes the body an equipotential) is not a simple problem. The existence of the symmetry for considered system may simplify the task. Thomson was the first to suggest an instructive solution for interacting electrons, reducing the 3D harmonic oscillator confinement to a circular (2D) harmonic oscillator [4]. He developed an analytical approach, which enables us to trace a self-organization for a small number of electrons (n < 10) in a family of rings (shells) with a certain number of electrons in each shell.

Nowadays, many ideas and concepts introduced in condensed matter physics can be realized and analyzed with high accuracy as a function of particle number and boundary properties. In fact, nanotechnology gave rise to emergence of lateral quantum dots creation which properties where not so obvious in ninetieth. Assuming a circular symmetry of such quantum dots, the first attempt to understand the distribution of electrons in such confined systems were based on the Monte Carlo

© Nikonov E. G., Nazmitdinov R. G., Glukhovtsev P. I., 2024

This work is licensed under a Creative Commons Attribution 4.0 International License https://creativecommons.Org/licenses/by-nc/4.0/legalcode

(MC) calculations for charged particles (ions and electrons) confined by 2D parabolic and hard-wall potentials (see, e.g., [5-7]). The results of these calculations confirmed the predictions based on the Thomson model for N < 52. The next step in the attempt to find the analytical description of the distribution of charged particles in the disc geometry have been done in [8]. In this paper the basic principles of self-organization of one-component charged particles, confined in disk have been proposed. A system of equations was derived, which allows to determine equilibrium configurations for an arbitrary, but finite, number of charged particles that are distributed over several rings. This approach reduces significantly the computational effort in minimizing the energy of equilibrium configurations and demonstrates a remarkable agreement with the values provided by molecular dynamics calculations. This paper gave a new impetus to activity in finding the bridge between the distribution of finite number of confined charged particles to their continuous limit (e.g., [9-15]).

From the analysis, based on MC and MD calculations for a relatively small number of charged particles, it follows that the number of stable configurations grows very rapidly with the number of particles. There are many local minima that have energies very close to the global minimum. These metastable states with lower (or higher) symmetry are found with much higher probabilities than the true ground state [9,10]. This picture is akin to a liquid-solid transition, when a rapid cooling gives rise to a glass-like disordered solid rather than a crystal with lower energy. Therefore, the description of this phenomenon requires the development of various not only analytical approaches but as well the effective numerical recipes with growing number of confined particles.

In this paper, we present a new approach for numerical-analytical analysis of the equilibrium configurations of charged particles (say, electrons) confined in a disk geometry. Using a model of stable configurations, which takes into account the interaction between shells of charged particles [8, 11], we obtain functional dependencies of the total number of particles of the system on the number of rings and those of the energy of the equilibrium configuration on the total number of particles. These dependencies make it possible to significantly simplify the search for the absolute minimum of the system for a given total number of one component charged particles.

2. System description

The physical formulation of the problem can be as follows. A system of similarly charged particles is given; they are located in a region with a cylindrical confining potential at the boundary. The configurations of the particles are determined by the Hamiltonian, in which the potential energy of the interparticle interaction dominates over the kinetic energy. It is necessary to find a stable configuration of N particles with the minimal energy inside the given region.

Thus, we consider two-dimensional system on a plane consisting of identically charged particles with mutual Coulomb interaction in the confining disk potential with the radius R. The Hamiltonian of such a system can be written as follows:

N N 1 N

H = E VextOi) + « 2 + E Tf (1)

1=1

¡<j

where rt = |rj| is the distance to the center of the region bounded by the potential V"ext(r), a = e2/4 nz0zr is a quantity characterizing the interaction of charges in the medium, and Tt is the particle kinetic energy. The confining potential Vext(r) is defined as follows:

VextM = { 0 r<R; (2)

oo, r >■ R.

To avoid a large number of metastable states (local minima), the system is considered at temperatures that are close to zero; at these temperatures, the potential energy dominates over the kinetic one. As a result, it is possible to rewrite the function of the total energy of the system (1) as follows:

N N 1

H = EVext(rd + a E r^TT. (3)

r* t*. _t*.

i=i ',}=i |£i N!

i<i

The problem is to find a stable configuration of the particle inside the given region with the minimally possible energy.

Using his model, Thomsom obtained equation (4) for analytical calculation of the coordinates of particles in the equilibrium state for one ring

= S

¡=1 J=i+1 sin - (li-jl)

anSn 4r

(4)

Sn =S

1

1 sin - k

n

Here, En(r) is the Coulomb energy of n particles with the charge e uniformly distributed along a circle with the radius r, a = e2l4ne0Er is the quantity characterizing the interaction of charged particles in the medium. Without loss of generality, electrons with the charge e are considered as charged particles below. An original approach to the calculation of equilibrium configurations and the corresponding energy was formulated in [8,11]; in this approach, the interaction between shells consisting of charged particles is also taken into account in addition to the energy of one ring (one shell). As a result of solving the (5)

% = 0, i = 2,...,p.

* = r? S

j=i+i

njE\

-T8sn. + n s

'rjE

[M)2] *[№?]

}=i

l)

(5)

Here, K = X-1(E = X{) are the complete elliptic integrals of the first (second) kind: Xp(x) = f^i2 dt(1 — x sin2 t)p/2; rt is the value of the ith optimal radius for the given stable configuration of charged particles; and nt is the number of particles in the ith shell.

As preliminary analysis showed, the solution of these equations makes it possible to significantly reduce the amount of computational work by finding the state that is closest to the equilibrium configuration. It is worth noting that this approach allows us to almost exactly determine both the equilibrium configuration and the total energy of the equilibrium state for N < 52.

In this paper, we develop a new approach and methods for calculating the coordinates of particles and the energy of the equilibrium configuration for an arbitrary finite number of particles using the obtained analytic dependencies of the distribution of particles and the energy of the equilibrium configuration on the total number of particles in the system for N < 1000.

3. Computational scheme

To further increase the efficiency and reduce the calculation time, the following modification of the traditional approach based on the molecular dynamics method [16] is proposed in this paper. The computational scheme consists of the following steps.

1. We use interpolation formulas (6) for calculation initial particle distribution

Nt (N) = atN 2 —bt, (6)

i 1 2 3 4 5 6

ai 2.7948 1.3439 1.1323 1.0127 0.9482 0.8517

bt 3.9444 7.2999 10.845 14.850 19.128 21.732

Formulas (6) were obtained by interpolating solutions [17,18] to equations (3) of the equilibrium configuration model [8,11].

2. Using the initial distribution of particles inside the circular region obtained in the previous step, we run ab initio calculations using Newtonian molecular dynamics. The dissipation of energy for cooling the system upon reaching the ground state with a minimum of potential and zero kinetic energy is modeled by adding to the equations of motion a viscous friction force

proportional to the speed (7), which realizes the corresponding outflow of energy from the system.

mri'' = -ViV(ri)-bfrl', Vie{1,2,3,...,N}. (7)

Here v(ij) = Vext(ri) + a ^^=1 \ . includes the external confining potential plus the Coulomb

terms, and the friction coefficient bf is the parameter controlling the quenching of velocities.

For calculations we used quenched velocity Verlet algorithm. Let at some moment of time t the coordinates rt(t) and velocities vt(t) of the particles be given, then at the moment of time t + At, where At is some fixed time step, we can obtain the coordinates rt(t +At) and velocities \i(t +At) using the following formulas (8), (9), (10).

fi(t) = -W,V(r,(t)). (8)

rit+At) = r () + V(t)At + ft) - bfVi(t))At2. (9)

(2 - bfAt) vi(t) + (f-(t) + f-(t + At)]At

v(t + At) =-2+bf^t-—. (10)

This is the so-called quenched molecular dynamics method. 3. The final stage of calculations includes a fairly large sequence of runs of molecular dynamics calculations described at the previous step, which is due to the existence of a number of metastable states that exponentially increases with the number of particles near the point of global minimum of energy.

In order to estimate the minimum number of runs Nruns to obtain the ground state energy value Emin as close as possible to the global minimum, the following procedure is used. First, a certain number of runs are performed depending on the number of particles in the system. Then the probability of finding the

minimum P(Emin < Eva\) is calculated. Then the estimate for the minimum number of runs has the following form Nruns = 1/P(Emin < Eval).

For a sufficiently large number of runs Nruns, in accordance with the Central Limit Theorem, the energy values of metastable states calculated in each run are distributed in accordance with the Gaussian distribution. Typically, the value Eavr -3a is used as Eval. In this case, for a system consisting, for example, of 1000 particles, Eavr -3a = 736980.1734 in reduced units (Eavr is an average value for the general sample for all runs. And a is the corresponding standard deviation.). At a given value of Eval, the minimum number of runs Nruns to obtain the global minimum energy must be greater than 741.

However, we use a different version of probability estimation and, accordingly, a minimum number of runs to practically guarantee obtaining the minimum energy value. This value was obtained using the asymptotic formula (11) for the minimum energy of the ground state, obtained by the authors [17, 18] based on an analysis of solutions to the equations of the equilibrium configuration model [8,11].

Emin(N)=^4N2 -2N3 + (11)

In this case, the estimate for Eval for 1000 particles will be equal to 736978.54, which is approximately equal to Eavr - 4a. The corresponding minimum number of runs Nruns to obtain the global minimum energy must be greater than 25237. This value for the number of runs Nruns looks more realistic for the case of the initial configuration of particles distributed randomly within the disk.

As an initial approximation, the particle distribution obtained by solving equations of the model of equilibrium configurations (5) for a certain number of rings is taken. After that, calculations are started using the quenched molecular dynamics method under the condition of a gradual decrease in the system temperature. When the zero temperature is reached, the calculations of the time evolution of the system are considered to be completed, after which the energy of the resulting equilibrium configuration of particles in the system is calculated.

For example, for the above-mentioned system of 1000 particles, we fixed the distribution of particles on the outer ring in the initial configuration N1 = 276 in accordance with the formula (6). As a result, we needed only 500 runs to obtain the minimum ground state energy equal to 736979.7283. The best published value known to us is 736977.7079. This value can be obtained, according to our estimates, for a number of runs no less than 25237 [14].

4. Conclusions

The effectiveness of the approach proposed in this article for calculating the global minimum energy for the ground state of the above-mentioned systems of charged particles is based on the following modifications of the traditional molecular dynamics method.

First, the initial configuration for molecular dynamics calculations was calculated using analytic solutions of the equations of the equilibrium configuration model [8,11]. Secondly, to cool the system, energy dissipation was carried out due to viscous friction forces. For this purpose, the Verlet velocity scheme with quenching was used. Finally, to estimate the number of molecular dynamics runs with different initial conditions, we used an asymptotic formula [17,18] for the ground state energy, obtained from an analysis of solutions to the equations of the equilibrium configuration model [8, 11]. As a result, the computational efficiency in terms of computation time has increased by more than two orders of magnitude.

The algorithms and programs developed by us can be used to numerically study the stability of systems of charged particles in various fields of physics, chemistry, molecular biology, and nanotechnology, including, for example, the study of nano-objects, such as quantum dots.

Funding: The work was supported within the framework of the Joint Institute for Nuclear Research project no. 05-6-11182014/2023 with the use of resources of the HybriLIT heterogeneous cluster of the Laboratory of Information Technologies, Joint Institute for Nuclear Research. Work by E.G.N. partially supported by the Basic Research Program of the National Research University "Higher School of Economics".

References

1. Bowick, M. J. & Giomi, L. Two-dimensional matter: order, curvature and defects. Advances in Physics 58, 449-563. doi:10.1080/ 000 187 3 0 9 0 30431 66 (2009).

2. Saarikoski, H., Reimann, S. M., Harju, A. & Manninen, M. Vortices in quantum droplets: Analogies between boson and fermion systems. Rev. Mod. Phys. 82, 2785-2834. doi:10. 1103/ RevModPhys.82.2785 (3 Sept. 2010).

3. Birman, J., Nazmitdinov, R. & Yukalov, V. Effects of symmetry breaking in finite quantum systems. Physics Reports 526,1-91. doi:10.1016/j.physrep.2012.11.005 (2013).

4. Thomson, J. XXIV. On the structure of the atom: an investigation of the stability and periods of oscillation of a number of corpuscles arranged at equal intervals around the circumference of a circle; with application of the results to the theory of atomic structure. Phil. Mag. 7, 237-265. doi:10.1080/14786440409463107 (39 1904).

5. Lozovik, Y. E. & Mandelshtam, V. A. Classical and quantum melting of a Coulomb cluster in a trap. Physics Letters A 165, 469-472. doi:10.1016/0375-9601(92)90349-Q (1992).

6. Bolton, F. & Rossler, U. Classical model of a Wigner crystal in a quantum dot. Superlattices and Microstructures 13,139-146. doi:10.1006/spmi.1993.1026 (1993).

7. Bedanov, V. M. & Peeters, F. M. Ordering and phase transitions of charged particles in a classical finite two-dimensional system. Phys. Rev. B 49, 2667-2676. doi:10.1103/PhysRevB.49.2667 (4 Jan. 1994).

8. Cerkaski, M., Nazmitdinov, R. G. & Puente, A. Thomson rings in a disk. Phys. Rev. E 91, 032312. doi:10.1103/PhysRevE.91.032312 (2015).

9. Amore, P. Comment on "Thomson rings in a disk". Phys. Rev. E 95, 026601. doi:10 . 1103/ PhysRevE.95.026601 (2 Feb. 2017).

10. Puente, A., Nazmitdinov, R. G., Cerkaski, M. & Pichugin, K. N. Reply to "Comment on'Thomson rings in a disk'". Phys. Rev.E 95, 026602. doi:10.1103/PhysRevE.95.026602 (2 Feb. 2017).

11. Nazmitdinov, R. G., Puente, A., Cerkaski, M. & Pons, M. Self-organization of charged particles in circular geometry. Phys. Rev.E 95, 042603. doi:10.1103/PhysRevE.95.042603 (2017).

12. Ash, B., Chakrabarti, J. & Ghosal, A. Static and dynamic properties of two-dimensional Coulomb clusters. Phys. Rev.E 96, 042105. doi:10.1103/PhysRevE.96.042105 (4 Oct. 2017).

13. Jami, P., Ash, B. & Ghosal, A. Significance of the nature of disorder on the universal features of the spatio-temporal correlations of two-dimensional Coulomb-clusters. Eur. Phys.J. B 96,45. doi:10.1140/epjb/s10051-023-00510-5 (4 2023).

14. Amore, P. & Zarate, U. Thomson problem in the disk. Phys. Rev. E108, 055302. doi:10. 1103/ PhysRevE.108.055302 (5 Nov. 2023).

15. Batle, J., Vlasiuk, O. & Ciftja, O. Correspondence between Electrostatics and Contact Mechanics with Further Results in Equilibrium Charge Distributions. Annalen derPhysik, 2300269. doi:10. 1002/andp.202300269 (Nov. 2023).

16. Frenkel, D. & Smit, B. Understanding Molecular Simulation: From Algorithms to Applications 664 pp. (Academic, New York, 2001).

17. Nikonov, E. G., Nazmitdinov, R. G. & Glukhovtsev, P. I. Molecular dynamics studies of equilibrium configurations of equally charged particles in planar systems with circular symmetry. Computer research and modeling 14, 609-618. doi:l0.20537/2076-7633-2022-l4-3-609-618 (2022).

18. Nikonov, E. G., Nazmitdinov, R. G. & Glukhovtsev, P. I. On the Equilibrium Configurations of Charged Ions in Planar Systems with Circular Symmetry. Journal of Surface Investigation: X-ray, Synchrotron and Neutron Techniques 17, 235-239. doi:l0.ll34/Sl02745l0230l0354 (2023).

To cite: Nikonov E.G., Nazmitdinov R. G., Glukhovtsev P. I., Computer studies of a dependence of equilibrium state structure on a number of particles for a two-dimensional system of charged particles confined in a disk potential, Discrete and Continuous Models and Applied Computational Science 32 (1)(2024)99-105.D01:10.22363/2658-4670-2024-32-1-99-105.

Information about the authors

Nikonov, Eduard G.—Doctor of Physical and Mathematical Sciences, Head of Sector MLIT JINR (e-mail: e.nikonov@jinr.ru, phone: +7(496)2164722, ORCID: https://orcid.org/0000-0001-7162-0344, ResearcherlD: C-4841-2016, Scopus Author ID: 6603099928)

Nazmitdinov, Rashid G.—Doctor of Physical and Mathematical Sciences, Leading Researcher BLTP JINR (e-mail: rashid@theor.jinr.ru, ORCID: https://orcid.org/000-0003-0489-7858, ResearcherlD: G-4860-2016,

Scopus Author ID: 7004352313)

Glukhovtsev Pavel I.—Master's degree student of Department of distributed information computing systems of Dubna State University (e-mail: pavelgl2018@gmail.com, ORCID: https://orcid.org/0009-0005-6424-4455)

УДК 538.911, 538.915, 519.6

PACS 07.05.Tp, 02.70.Ns, 71.90.+q, 21.60.Jz, 24.60.Lz, 27.80.+W DOI: 10.22363/2658-4670-2024-32-1-99-105

EDN: BXKOPH

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

Э. Г. Никонов1,23, Р. Г. Назмитдинов1,2, П. И. Глуховцев2

1 Объединенный иституут ядерных исследований, ул. Жолио-Кюри, д. 6, Дубна, 141980, Российская Федерация

2 Государственный университет «Дубна», ул. Университетская, д. 19, Дубна, 141980, Российская Федерация

3 Национальный исследовательский университет «Высшая школа экономики», ул. Таллинская, д. 34, Москва, 123458, Российская Федерация

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

Ключевые слова: атом Томсона, кристалл Вигнера, молекулярная динамика