Строим галактики…

Столкновение галактик EN

Столкновение галактик

Через четыре миллиарда лет Андромеда налетит на Млечный Путь. Здесь это происходит прямо сейчас: видеокарта считает притяжение каждой из сотен тысяч частиц ко всем остальным.

Что здесь видно

Тяните, чтобы облететь, колёсико — ближе-дальше. Ниже — статья с формулами.

Как это устроено

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

Масштабы и единицы

Галактики огромны и медленны. Удобно считать в единицах, где гравитационная постоянная $G = 1$, длина — килопарсек (3262 световых года), масса — $10^{10}$ масс Солнца. Тогда единица скорости — 207 км/с, а единица времени — 4,7 миллиона лет. Млечный Путь в этих единицах — диск радиусом около 15 и массой 5 внутри гало массой 100, Андромеда чуть крупнее.

Одна частица — не звезда. На сайте частица звёзд весит несколько миллионов масс Солнца, частица тёмной материи — в восемь раз больше, частица газа — в четыре раза меньше. Это «облака» из множества звёзд, которые движутся вместе. В фильме частиц в 80 раз больше и каждая во столько же раз легче.

Гравитация: каждая тянет каждую

Ускорение частицы $i$ — сумма притяжений всех остальных:

$$\mathbf a_i = G \sum_{j \ne i} m_j\, g(r_{ij})\, (\mathbf x_j - \mathbf x_i),\qquad g(r) = \frac{1}{r^3}\ \text{вдали}.$$

Вблизи сила смягчена: частица — размазанное облако, а не точка, и две частицы не могут разогнать друг друга до бесконечности. Мы используем кубический сплайн из кода GADGET-2 (Спрингел, 2005): на расстояниях больше $h = 2{,}8\,\varepsilon$ сила ньютоновская точно, а в центре потенциал как у сферы Пламмера радиуса $\varepsilon$. На сайте $\varepsilon$ — 80 пк для звёзд и 250 пк для тёмной материи; с ростом числа частиц смягчение уменьшается как $N^{-1/3}$.

Дерево Барнса–Хата

Прямая сумма по всем парам — это $N^2$: для 200 тысяч частиц сорок миллиардов слагаемых на каждый шаг. Барнс и Хат (1986) заметили, что далёкую группу частиц можно заменить одной: её полной массой в центре масс плюс поправкой на форму — квадрупольным моментом:

$$\Phi(\mathbf y) \approx -\frac{M}{|\mathbf y|} - \frac{\mathbf y^{\mathsf T} Q\, \mathbf y}{2|\mathbf y|^5},\qquad Q_{ab} = \sum_k m_k \big(3 d_a d_b - |\mathbf d|^2 \delta_{ab}\big).$$

Пространство делится восьмеричным деревом: куб — на восемь кубов, пока в ячейке не останется до 16 частиц. Узел принимается целиком, если группа частиц дальше от его центра масс, чем $l/\theta + \delta$ ($l$ — размер узла, $\delta$ — сдвиг центра масс от центра ячейки, $\theta = 1$ на сайте, $0{,}9$ в фильме), иначе раскрывается. В итоге на частицу приходится около тысячи взаимодействий вместо двухсот тысяч.

Дерево строится на видеокарте заново каждый шаг. Каждая координата переводится в 20-битное целое, биты трёх осей перемешиваются в 60-битный ключ Мортона. Частицы сортируются по ключу поразрядной сортировкой. После этого частицы любой ячейки дерева идут подряд, а её дочерние ячейки находятся двоичным поиском. Обход идёт без стека: у каждого узла есть ссылка «куда дальше», и 64 частицы группы проходят дерево вместе, принимая одинаковые решения.

Время: leapfrog

Положения и скорости продвигаются схемой «толчок–шаг–толчок»:

$$\mathbf v \mathrel{+}= \mathbf a\,\tfrac{\Delta t}{2},\qquad \mathbf x \mathrel{+}= \mathbf v\,\Delta t,\qquad \mathbf a = \mathbf a(\mathbf x),\qquad \mathbf v \mathrel{+}= \mathbf a\,\tfrac{\Delta t}{2}.$$

Схема симплектическая и обратимая во времени: ошибка энергии не накапливается, а колеблется около нуля. Поэтому орбиты не «расползаются» за тысячи шагов. Шаг на сайте — полмиллиона лет.

Как собрать галактику

Галактика должна стартовать в равновесии, иначе она начнёт пульсировать ещё до встречи. Каждая собирается из пяти частей:

  • гало тёмной материи и балдж — профиль Хернквиста $\rho = \dfrac{M a}{2\pi r (r + a)^3}$; масштаб гало получается из массы $M_{200}$ и концентрации $c = 10$ (как у гало в космологических симуляциях);
  • звёздный диск — экспоненциальный: $\Sigma(R) = \Sigma_0 e^{-R/R_d}$, по высоте — $\operatorname{sech}^2(z/z_0)$;
  • газовый диск — вдвое шире звёздного и тоньше;
  • сверхмассивная чёрная дыра в центре: 4,3 млн масс Солнца у Млечного Пути, 140 млн у Андромеды.

Скорости частиц гало и балджа берутся из функции распределения по формуле Эддингтона. Для сферической системы с изотропными скоростями она восстанавливает распределение по энергиям $\mathcal E = \Psi - v^2/2$ из плотности:

$$f(\mathcal E) = \frac{1}{\sqrt 8\,\pi^2} \frac{d}{d\mathcal E} \int_0^{\mathcal E} \frac{d\rho}{d\Psi} \frac{d\Psi}{\sqrt{\mathcal E - \Psi}}.$$

Потенциал $\Psi$ здесь общий — гало, балджа, дисков и дыры вместе. Проверка: для одиночной сферы Хернквиста численная $f(\mathcal E)$ совпадает с известной точной формулой (Хернквист, 1990) до процента.

Диск вращается со скоростью, которую дают уравнения Джинса (Хернквист, 1993). Вертикальная дисперсия — из равновесия слоя, $\sigma_z^2 = \pi G \Sigma z_0$. Радиальная задаётся параметром устойчивости Тумре:

$$Q = \frac{\sigma_R\, \kappa}{3{,}36\, G \Sigma} = 1{,}3,$$

где $\kappa$ — эпициклическая частота. Скорость вращения звёзд чуть меньше круговой $v_c$ — это асимметричный дрейф:

$$\bar v_\phi^2 = v_c^2 + \sigma_R^2 \Big(1 - \frac{\kappa^2}{4\Omega^2} - \frac{2R}{R_d}\Big).$$

Круговая скорость тонкого экспоненциального диска выражается через функции Бесселя (Фримен, 1970): $v_c^2 = 4\pi G \Sigma_0 R_d\, y^2 [I_0(y)K_0(y) - I_1(y)K_1(y)]$, $y = R/2R_d$. У нашего Млечного Пути на расстоянии Солнца выходит 230 км/с, как у настоящего. Предоставленный самому себе, такой диск за пару сотен миллионов лет обзаводится спиральными рукавами и перемычкой-баром. Настоящий Млечный Путь как раз галактика с баром.

Где Андромеда

Андромеда сейчас в 780 кпк (2,5 млн световых лет). Её скорость к нам вдоль луча зрения известна точно: −301 км/с относительно Солнца, или −109 км/с относительно центра Галактики. С поперечной скоростью сложнее: её измеряют по смещению звёзд Андромеды на небе за годы наблюдений. «Хаббл» дал 17 км/с (ван дер Марел и др., 2012), «Гайя» — 57 км/с (2019) и 80 км/с (2021) с ошибками в десятки км/с. По умолчанию стоят 17 км/с, в настройках можно поставить другую.

Диски ориентированы как настоящие. Млечный Путь вращается по часовой стрелке, если смотреть с северного полюса Галактики. Ось Андромеды получается из её наклона ($i = 77{,}5^\circ$) и позиционного угла ($37{,}7^\circ$), если учесть, что ближе к нам северо-западная сторона диска и что северо-восточная половина от нас удаляется. Выходит направление $(l, b) \approx (241^\circ, -30^\circ)$, как у ван дер Марела.

Первые три миллиарда лет, пока галактики далеко, мы считаем как задачу двух тел с протяжёнными гало и динамическим трением. Дальше, с 300 кпк, работает полная симуляция.

Столкнутся ли они на самом деле, неизвестно. Савала и соавторы (2025) учли ошибки измерений и притяжение M33 и Большого Магелланова Облака. Вероятность слияния за 10 миллиардов лет у них вышла около 50%. Здесь показан вариант, в котором слияние происходит.

Динамическое трение

Почему галактики, разлетевшись после первого пролёта, возвращаются? Массивное тело, летящее сквозь облако лёгких частиц, собирает за собой «кильватер» — сгущение, которое тянет его назад. Чандрасекар (1943) получил это торможение для однородной среды:

$$\frac{d\mathbf v}{dt} = -\frac{4\pi G^2 M \rho \ln\Lambda}{v^3}\Big[\operatorname{erf}(X) - \frac{2X}{\sqrt\pi} e^{-X^2}\Big]\mathbf v,\qquad X = \frac{v}{\sqrt2\,\sigma}.$$

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

Приливные хвосты

Звёзды почти никогда не сталкиваются: расстояния между ними в десятки миллионов раз больше их самих. Хвосты рождает прилив: ближний край диска притягивается к соседке сильнее центра, дальний — слабее. Сильнее всего вытягиваются звёзды, чьё вращение совпадает с направлением пролёта: они дольше всего «едут» рядом со встречной галактикой. Это объяснили Алар и Юри Тумре (1972), посчитав пробные частицы вокруг двух точечных масс на одном из первых компьютеров. Их модели «Антенн» и «Мышей» есть в пресетах.

Газ: гидродинамика сглаженных частиц

Газ, в отличие от звёзд, давит и сталкивается. Мы считаем его методом SPH (Люси, 1977; Гингольд и Монаган, 1977): частица газа — размазанный комок, плотность — сумма вкладов соседей с ядром Wendland C2:

$$\rho_i = \sum_j m_j W(r_{ij}, h_i),\qquad W(r,h) = \frac{21}{2\pi h^3}\Big(1 - \frac rh\Big)^4\Big(1 + \frac{4r}{h}\Big).$$

Длина сглаживания $h$ подстраивается так, чтобы соседей было около 48. Газ изотермический: температура около $10^4$ К (её поддерживает излучение звёзд), давление $P = c_s^2 \rho$, $c_s = 10$ км/с. Ускорение от давления симметрично по парам, поэтому импульс сохраняется точно:

$$\frac{d\mathbf v_i}{dt} = -\sum_j m_j \Big(\frac{P_i}{\rho_i^2} \nabla_i W(h_i) + \frac{P_j}{\rho_j^2} \nabla_i W(h_j) + \Pi_{ij}\, \nabla_i \bar W\Big).$$

$\Pi_{ij}$ — искусственная вязкость Монагана: она включается только для сближающихся пар и превращает лобовое столкновение потоков в ударную волну. Чтобы вязкость не тормозила обычное вращение диска, её гасит переключатель Балсары, который смотрит на соотношение $|\nabla\cdot\mathbf v|$ и $|\nabla\times\mathbf v|$.

Газ требует шагов мельче, чем звёзды: условие Куранта $\Delta t < 0{,}25\, h / v_{\text{sig}}$. Поэтому внутри одного гравитационного шага газ делает до восьми подшагов. Ещё две предосторожности: пол давления по Джинсу $c_{\text{eff}}^2 = \max(c_s^2,\ 3G\rho h^2)$ не даёт газу схлопнуться в комки меньше разрешения, а там, где Куранту всё равно тесно, толчок за подшаг ограничен половиной сигнальной скорости.

Звёзды из газа

Где газ плотнее примерно 0,1 атома водорода на кубический сантиметр и сжимается, рождаются звёзды. Скорость — закон Шмидта–Кенникатта: за время свободного падения $t_{ff} = \sqrt{3\pi/32G\rho}$ в звёзды превращается полтора процента газа:

$$\dot\rho_\star = \varepsilon\, \frac{\rho}{t_{ff}},\qquad p = 1 - e^{-\varepsilon \Delta t / t_{ff}}.$$

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

Как получается картинка

Свет звёзд. Частица звёзд — население одного возраста. Его цвет и светимость на единицу массы взяты по моделям Брюзуаль и Шарло (2003): молодое население голубое и в сотни раз ярче старого, жёлто-оранжевого. Между 10 млн и 10 млрд лет светимость падает примерно как $t^{-0{,}8}$. Звёзды моложе 10 млн лет ионизуют газ вокруг: светится розовая область H II (линии водорода Hα и Hβ). Её яркость спадает от центра плавно, по профилю Моффата $(1 + r^2/a^2)^{-2}$: яркий узел в несколько десятков парсек и ореол, гаснущий без края.

Пыль. Пыль идёт вместе с газом — около процента его массы. Она поглощает свет за собой, синий сильнее красного, по закону Карделли, Клейтона и Матиса (1989); нормировка $N_H/E(B{-}V) = 5{,}8 \cdot 10^{21}$ см⁻² — по Болину и др. (1978). Поэтому пылевые полосы тёмно-бурые, а свет за ними краснеет.

Сборка кадра. Каждая частица рисуется размытым диском — проекцией своего ядра сглаживания. Частицы сортируются по расстоянию от камеры и накладываются от дальних к ближним: звёзды добавляют свет, газ ослабляет всё, что за ним, отдельно в каждом цвете, $C \leftarrow E + C\,e^{-\tau}$. Часть света звёзд рисуется острой точкой: в настоящей галактике свет не гладкий, его несут яркие гиганты и скопления.

Экспозиция и цвет. Экспозиция подстраивается так, чтобы самые яркие 0,1% пикселей едва касались белого. Яркость растягивается функцией $\operatorname{asinh}$, как в цветных снимках обзора SDSS (Лаптон и др., 2004): слабое — линейно, яркое — логарифмически, без сдвига цветов. Самое яркое плавно уходит в белый, без плоской «отсечки». Свет ярче белого рассеивается в ореол, как в оптике телескопа: его яркость падает с расстоянием примерно как $1/r^2$. Поэтому яркие узлы и ядра светятся, а диск вокруг остаётся резким. Локальный контраст слегка подчёркнут, как при обработке астрофото. Цвета — «как у камеры»: свет делится на три фильтра (B, V, I), как в известных снимках «Хаббла». Глазом галактика выглядела бы бледнее и менее цветной.

Фон. Позади — далёкие галактики. Два десятка самых ярких соседей (M81 и M82, NGC 253, Центавр A, M101, Сомбреро, скопление Девы) стоят на своих настоящих местах неба и в настоящих размерах. Остальные — сочинённое по реальной статистике числа галактик глубокое поле, собранное в группы и скопления. Это декорация, а не каталог.

Что упрощено

  • Частица — не звезда, а облако в миллионы солнц на сайте и в десятки тысяч в фильме. Поэтому мелкие детали — отдельные скопления, тонкие пылевые нити — размыты по масштабу смягчения.
  • Газ изотермический. Нет охлаждения ниже $10^4$ К, нет нагрева ударными волнами до миллионов градусов, нет горячего газа в гало. При слиянии реальный газ частично нагреется и выключит звёздообразование надолго.
  • Нет обратной связи: сверхновые и ветры молодых звёзд не разгоняют газ.
  • Чёрные дыры — просто тяжёлые частицы. Они опускаются к центру и сливаются сами, но аккрецию газа и свечение квазара мы не считаем.
  • Окружения нет. Нет M33, Магеллановых Облаков и других спутников, нет космологического фона: Местная группа одна в пустоте.
  • Шаг общий для всех (кроме подшагов газа). В самых плотных ядрах орбиты считаются грубее: в настоящих кодах у каждой частицы свой шаг.
  • Цвета звёздного населения взяты для солнечной металличности. Пыль везде одного состава.

Как проверено

  • Функция распределения по Эддингтону совпадает с точной формулой Хернквиста до 1% (до 5% у самого центра), гало и балдж в своём потенциале держат вириальное равновесие $2K = |W|$ с точностью 3%.
  • Круговая скорость диска Фримена сверена с прямым суммированием по диску (0,2%) и с производной его потенциала ($10^{-5}$).
  • Положение и ориентация Андромеды сверены с данными ван дер Марела: $(l, b) = (121{,}17^\circ, -21{,}57^\circ)$, ось диска $(241^\circ, -30^\circ)$.
  • Силы дерева на видеокарте сравниваются с прямой суммой в двойной точности: медианная ошибка $2 \cdot 10^{-4}$, 99% частиц — лучше $10^{-3}$ при $\theta = 0{,}6$; при $\theta = 1$ — $8 \cdot 10^{-4}$ и $4 \cdot 10^{-3}$.
  • Плотность SPH на видеокарте совпадает с прямым подсчётом до $10^{-7}$, суммарный импульс гидродинамических сил — ноль с точностью $10^{-9}$.
  • Изолированный Млечный Путь за миллиард лет сохраняет энергию до $6 \cdot 10^{-5}$, а всё столкновение — до $2 \cdot 10^{-4}$. Первое сближение наступает через 3,9 млрд лет от сегодня — как в расчётах ван дер Марела (3,87).

Источники

Barnes J., Hut P. (1986) A hierarchical O(N log N) force-calculation algorithm. Nature 324.
Bédorf J., Gaburov E., Portegies Zwart S. (2012) A sparse octree gravitational N-body code that runs entirely on the GPU. J. Comput. Phys. 231.
Springel V. (2005) The cosmological simulation code GADGET-2. MNRAS 364.
Hernquist L. (1990) An analytical model for spherical galaxies and bulges. ApJ 356; (1993) N-body realizations of compound galaxies. ApJS 86.
Freeman K. C. (1970) On the disks of spiral and S0 galaxies. ApJ 160.
Toomre A., Toomre J. (1972) Galactic bridges and tails. ApJ 178.
Chandrasekhar S. (1943) Dynamical friction. ApJ 97.
van der Marel R. P. et al. (2012) The M31 velocity vector I–III. ApJ 753; (2019) ApJ 872.
Sawala T. et al. (2025) No certainty of a Milky Way–Andromeda collision. Nature Astronomy.
Monaghan J. J. (1992) Smoothed particle hydrodynamics. ARA&A 30; Balsara D. (1995) J. Comput. Phys. 121.
Springel V., Hernquist L. (2003) Cosmological SPH simulations: a hybrid multiphase model for star formation. MNRAS 339.
Bruzual G., Charlot S. (2003) Stellar population synthesis at the resolution of 2003. MNRAS 344.
Cardelli J., Clayton G., Mathis J. (1989) The relationship between infrared, optical, and ultraviolet extinction. ApJ 345.
Lupton R. et al. (2004) Preparing red-green-blue images from CCD data. PASP 116.