Кулоновский кристалл: заряды в ловушке

Возьмите несколько одинаковых зарядов – они отталкиваются и разлетятся в разные стороны. Но если поместить их в параболическую ловушку, которая тянет каждую частицу к центру, заряды соберутся вместе и при остывании выстроятся в аккуратные кольца-оболочки. Получается «искусственный атом»: ядро здесь заменяет сама ловушка, а электронные оболочки – кольца зарядов. Это, по сути, модель атома Томсона.

Ионы в ловушках Пауля – основа одного из типов квантовых компьютеров. Похожие структуры наблюдают в пылевой плазме: микрочастицы в газовом разряде выстраиваются в такие же кольца, и опыты с ними ставят даже на МКС. А электроны на поверхности жидкого гелия образуют вигнеровский кристалл – ту же решётку из отталкивающихся зарядов.

Как меняется число колец и то, сколько зарядов сидит на каждом кольце, когда частиц становится больше, – есть ли у этой последовательности «магические» числа. И как кристалл плавится при нагреве: среднее смещение частиц от своих узлов сначала растёт плавно, а потом скачком – это критерий Линдемана.

1. Разберитесь с физикой: одинаковые заряды отталкиваются, а параболическая ловушка тянет каждую частицу к центру – вместе они образуют «искусственный атом». Прочитайте постановку задачи и выпишите, из чего складывается энергия системы.
2. Выведите и запишите потенциал $U = \sum_i \frac{m\omega^2 r_i^2}{2} + \sum_{i \lt j} \frac{kq^2}{r_{ij}}$ и договоритесь о единицах. Кулоновское поле дальнодействующее, поэтому обрезание потенциала здесь не годится.
3. Возьмите каркас задачи MD и поменяйте в ней потенциал Леннарда-Джонса на кулоновский плюс ловушку, уберите периодические границы и обрезание потенциала. Перебор всех пар $O(N^2)$ при N до ~100 не проблема.
4. Реализуйте поиск основного состояния через вязкое трение (затухающую динамику) либо медленное охлаждение термостатом. Проверьте на двух зарядах: каждый должен стоять на расстоянии $a$ от центра, где $a^3 = \frac{kq^2}{4m\omega^2}$.
5. Проверьте на трёх зарядах: должен получиться правильный треугольник с радиусом $a^3 = \frac{kq^2}{\sqrt{3}\,m\omega^2}$. Сравните с таблицей оболочек: N = 6 – (1,5), N = 8 – (1,7), N = 9 – (2,7) (минимум).
6. Убедитесь, что программа находит именно минимум энергии: при N = 6 конфигурации (1,5) и (6) почти равны по энергии. Делайте несколько случайных стартов и отжиг – иначе легко застрять в метастабильной конфигурации.
7. Постройте таблицу конфигураций для N = 2…20 и сверьте её с рядом: N = 2…5 – одно кольцо, N = 10 – (2,8), N = 12 – (3,9), N = 15 – (5,10), N = 16 – (1,5,10) – первая третья оболочка (минимум).
8. Продвиньте расчёт оболочек до N ≈ 20–30 и займитесь плавлением: считайте среднее смещение частиц от своих узлов при нагреве. Сравнивайте энергии конфигураций, чтобы не принять метастабильное состояние за основное.
9. Расширение: сравните радиальные и угловые колебания – кольца сначала проворачиваются друг относительно друга, а потом плавятся.
10. Оформите работу: анимацию того, как заряды застывают в кольца, таблицу конфигураций для N = 2…20, график среднего смещения частиц от температуры; в выводах оговорите метастабильные состояния, в которых программа может застрять, и подготовьте доклад.