Формирование медленных баров в анизотропных системах

Смирнов Антон Александрович

Формирование медленных баров в анизотропных системах

Функция распределения для обобщенно-политропных моделей использована для построения серии численных моделей анизотропных дисков. Показано, что при любой степени радиальной вытянутости орбит звёзд такие системы являются неустойчивыми относительно образования бара. Этот результат полностью отличается от результатов прошлых работ. Найдены скорости узора и амплитуды образующегося бара. Рассчитаны начальные распределения прецессий орбит. Показано, что в таких системах соблюдается все условия для возникновения неустойчивости радиальных орбит, ответственной за формирование медленного бара.

Математика

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d36645f1be77c40d58e4b

UUID: 40752d15-16e7-47db-ba97-bbeb30d3b6af

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 563

Смирнов Антон Александрович

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 8723184 bytes

Поделиться работой

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ КАФЕДРА АСТРОФИЗИКИ Дипломная работа Смирнова Антона Александровича Моделирование медленных баров в анизотропных системах Научный руководитель: доктор физ.-мат. наук Сотникова Н. Я. Рецензент: кандидат физ.-мат. наук Родионов С. А. Санкт-Петербург 2016

SAINT-PETERBURG STATE UNIVERSITY DEPARTMENT OF ASTROPHYSICS Graduation thesis Smirnov Anton Alexandrovich Modeling of slow bars in anisotropic systems Scientific supervisor: Doctor of Sciences Sotnikova N. Ya. Reviewer: Candidate of Sciences Rodionov S. A. Saint-Petesburg 2016

Содержание 1. Введение 5 2. Обобщенно-политропные модели звездных дисков 2.1. Построение модели. Нахождение потенциала обобщенно-политропной модели звездного диска . . . . . . . 2.2. Начальное распределение по координатам . . . . . . . . . . 2.3. Начальное распределение по скоростям . . . . . . . . . . . . 7 7 10 11 3. Расчет времен релаксации и оптимальных параметров сглаживания моделей 3.1. Численная релаксация в задаче N тел . . . . . . . . . . . . . 3.2. Расчет оптимальных параметров сглаживания . . . . . . . . 3.3. Оценка времен релаксации . . . . . . . . . . . . . . . . . . . 15 15 17 18 4. Численное моделирование обобщенно-политропных моделей звездных дисков 4.1. Нормировка расстояний, скоростей и времени на размер системы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2. Решение задачи N тел. Пакет N EM O . . . . . . . . . . . . 4.3. Неустойчивость обобщенно-политропных моделей звездных дисков . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4. Расчет скорости узора и амплитуды гармоники, соответствующей бару . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. Обсуждение результатов 5.1. Неустойчивость радиальных орбит. Необходимые условия ее реализации . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2. Расчет траекторий частиц в фиксированном начальном потенциале. Определение скоростей прецессии частиц, эксцентриситетов их орбит и адиабатических инвариантов . . . . . 5.2.1. Интегрирование уравнений движения . . . . . . . . . 5.2.2. Расчет скорости прецессии. Смещение апоцентра орбиты . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.3. Расчет адиабатического инварианта орбиты . . . . . 3 22 22 22 23 29 31 31 33 33 34 36

5.2.4. Расчет эксцентриситета орбиты . . . . . . . . . . . . 5.3. Эксцентричность орбит в обобщенно-политропных моделях звездных дисков . . . . . . . . . . . . . . . . . . . . . . . . . 5.4. Зависимость скорости прецесии от углового момента при фиксированном адиабатическом инварианте . . . . . . . . . . . . 5.5. Сравнение скорости узора бара и скоростей прецессии частиц 36 38 38 41 6. Выводы 42 7. Благодарности 44 4

1. Введение Из данных наблюдений известно, что примерно 50% спиральных галактик имеют характерную вытянутую крупномасштабную структуру, напоминающую “перемычку”. Эта “перемычка” называется баром. Свойства бара зависят от Хаббловского типа подстилающей галактики (Бинни и Тремейн, 2008). Как правило, у ранних типов галактик профиль поверхностной яркости бара плоский, у поздних — экспоненциальный. Обычно бары сильно вытянуты. Например, для Галактики отношение большой полуоси к малой равно примерно 3:1. Несколько связанных между собой вопросов до сих пор находятся в центре внимания астрономов. Каковы механизмы образования баров? Почему они имею такие характеристики (скорость узора, размер, профиль плотности), которые наблюдаются? По скорости узора бары разделяют на два вида: “быстрый” и “медленный”. Скорость узора “медленного” бара меньше характерных скоростей прецессии орбит звезд диска. Скорость узора “быстрого” бара может превышать максимальную скорость прецессии орбит. Размер бара определяется положением Линбладовских резонансов nr (Ω(r0 ) − Ωp ) = ±κ, (1) Ω(r0 ) = Ωp . (2) и радиусом коротации Здесь Ω(r0 ) — угловая скорость вращения звезды на радиусе r0 , κ — ее эпициклическая частота, Ωp — скорость узора, nr — порядковый номер резонанса. Линден-Беллом и Калнайсом (1972) было показано, что постоянная волна плотности индуцирует изменение углового момента в система только вблизи резонансов. Кроме того, вблизи резонансов Линдблада включается механизм затухания волн плотности, аналогичный затуханию Ландау в физике плазмы. В результате “медленный” бар ограничен внутренним Линдбладовским резонансом 2:1 (nr = 2), а “быстрый” — радиусом коротации или резонанса 4:1 (nr = 4). В центральных областях галактик орбиты звезд могут сильно отличаться от круговых. Это приводит к тому, что резонанс “расплывается”. 5

Образуется зона резонанса, в которой находятся резонансные орбиты. При фиксированном эксцентриситете орбиты удовлетворяют соотношениям (1) и (2). Изучению быстрых баров посвящено множество теоретических и численноэкспериментальных работ. Как правило, бары, наблюдаемые в галактиках, — это быстрые бары (Селвуд, 2013). Что касается медленных баров, которые могут образовываться в центральных областях звездных дисков, то работ, как теоретических, так и экспериментальных на этот счет практически нет. В качестве механизма образования медленных баров Линден-Беллом (1979) была предложена неустойчивость радиальных орбит. При определенных условиях вытянутые орбиты в центральных областях галактик можно рассматривать как медленно прецессирующие “спицы”. Под действием взаимной гравитации такие орбиты начинают “налипать” друг на друга, образуя при этом медленный бар. Эффективный способ исследования подобного рода неустойчивостей — это численное моделирование, позволяющее исследовать различные нелинейные эффекты, которые трудно учесть в теории. Однако, на этом пути тоже есть трудности. Часто оказывается, что в ходе расчетов проявляются различные численные эффекты, которые совсем не похожи на то, что наблюдается в реальности (Вайенберг & Катс, 2006, Понтзен и др., 2015). Таким образом выбор параметров модели становится одной из ключевых задач численного моделирования. С целью исследования неустойчивости радиальных орбит в настоящей работе были рассмотрены равновесные обобщенно-политропные модели звездных дисков. Они определяются профилем фазовой плотности f (E, L):  A Lβ (E − E)α , L ≥ 0 , 0 f (E, L) = 0, L < 0 . (3) ~ = |~r × V~ |), E — полная энерЗдесь L — угловой момент звезды (L = |L| гия звезды (E = Φ(r) + V 2 /2), E0 = Φ(R) — потенциальная энергия на границе системы, R — размер системы, V — скорость звезды, r = |~r| — модуль радиус-вектора звезды, α, β — безразмерные параметры, A — ко6

эффициент пропорциональности, зависящий от полной массы системы. В выражении (3) предполагается, что все звезды вращаются в одну сторону. Эти модели предоставляют широкие возможности для исследования неустойчивости радиальных орбит. Степень вытянутости орбит можно регулировать параметром β: β= Vϕ2 (r) Vr2 (r) −1 (4) где Vϕ2 (r) — среднее значение квадрата скорости звезд в азимутальном направлении на заданном радиусе r, а Vr2 (r) — в радиальном; β = 0 соответствует изотропному распределению скоростей, β = 1 — строго радиальным орбитам. Такие модели двумерного диска позволяют исследовать анизотропные системы с заданной степенью анизотропии во всем диске. Ранее подобные модели были рассмотрены в работе Поляченко и Поляченко (1994). Авторами была получена граница по неустойчивости для параметра анизотропии β ≈ −0.5. В настоящей работе полностью пересматривается данный результат. На примере моделей с высоким пространственным разрешением показано, что при любой степени анизотропии орбит β такие системы являются неустойчивыми. Найдены скорости узора и амплитуды образующегося бара. Рассчитаны начальные распределение скоростей прецессий частиц. Показано, что в таких системах соблюдаются все условия для возникновения неустойчивости радиальных орбит. 2. Обобщенно-политропные модели звездных дисков 2.1. Построение модели. Нахождение потенциала обобщенно-политропной модели звездного диска Для построения равновесной обобщенно-политропной модели звездного диска нужно найти потенциал Φ(r), отвечающий профилю фазовой плотности (3). По определению, потенциал Φ(r) можно представить в виде ин7

теграла: Z R Z 2π Φ(r) = −G 0 ρ(r0 )r0 dr0 dφ p , r2 + r02 − 2rr0 cos φ 0 (5) где ρ(r) — поверхностная плотность на расстоянии r от центра, R — размер системы. Еще одна связь потенциала Φ(r) и поверхностной плотности ρ(r) получается путем интегрирования фазовой плотности по пространству скоростей (Бисноватый-Коган, 1975): ρ(r) = ACrβ (E0 − Φ0 (r))β/2+α+1 , (6) β+1 β+1 1 , )B( , α+1) и B(x, y) — функция Эйлера. Будем где C = 2β/2 B( 2 2 2 рассматривать обобщенно-политропные модели, для которых параметры α и β связаны следующим образом: β/2 + α = 0. (7) В этом случае поверхностная плотность ρ(r) будет линейно зависеть от потенциала Φ(r): ρ(r) = ACrβ (E0 − Φ0 (r)). (8) После подстановки выражения для поверхностной плотности (8) в выражение для потенциала (5) получим интегральное уравнение, где неизвестным является функция потенциала Φ(r): R Z 2π Z Φ(r) = −ACG 0 0 r0β (E0 − Φ0 (r0 ))r0 dr0 dφ p , r2 + r02 − 2rr0 cos φ (9) Будем использовать следующую систему единиц: G = 1, A · C · G = 1, Φ(R) − Φ(0) = 1. Для упрощения записи введем замену переменных: r̃ = r/R, Ũ (r̃) = U (r) , E0 − Φ(0) 8 (10) (11)

B(r̃) = Ũ (r̃)/, (12) где R — размер системы, r — расстояние от центра, r̃ — расстояние от центра в единицах радиуса, Ũ (r̃) — нормированный потенциал, U (r) = Φ(R) − Φ(r) — относительный потенциал. Имея ввиду, что U (0) = 1, легко показать, что нормировочный множитель = 1/B(0). Таким образом, нормированный потенциал Ũ (r̃) можно выразить через B-функцию Ũ (r̃) = B(r̃) . B(0) (13) В новых переменных интегральное уравнение (9) запишется в виде: Z 1Z 2π 1 − B(r̃) = −λ 0 ξ β+1 B(ξ)dξdφ p , ξ 2 + r̃2 − 2ξr̃ cos φ 0 (14) где введено обозначение λ = Rβ+1 . Интегральное уравнение (14) решается относительно функции B(r̃) и размера системы R. Затем нормированный потенциал Ũ (r̃) находится из соотношения (13). Решение интегрального уравнения (14) представляет из себя сложную задачу. Она была аккуратна решена А. Кошкиным в дипломной работе 2006 года. В данной работе использовались найденные Кошкиным значения B(r̃) функции и размеров систем R при различных β. В работе применен алгоритм, позволяющий уточнить значения B(r̃). Имея начальное приближенное решение уравнения (14), можно интеграл в правой части взять численно. Внутренний интеграл является эллиптическим К-интегралом. Второе интегрирование выполняется гауссовыми квадратурами. После взятия обоих интегралов, получаем уточненные значения функции B(r̃). Уточнение значений B(r̃) проводилось до тех пор, пока значения справа и слева от знака равенства в уравнении (14) не совпадали с точностью до восьмого знака. Функция B(r̃) позволяет найти поверхностную плотность ρ(r̃) и кумулятивное распределение массы M (r̃): ρ(r̃) = Rβ · r̃β · 9 B(r̃) , B(0) (15)

1 M (r̃) = 2πRβ+2 λ B(0) Z r̃ ξ β+1 B(ξ)dξ, (16) 0 Величина M (r̃) равна массе, содержащейся в круге радиуса r̃. Для построения начального распределения по координатам мы будем использовать нормированное распределение массы µ(r): R r̃ β+1 ξ B(ξ)dξ M (r̃) µ(r̃) = , = R01 β+1 B(ξ)dξ Mtot ξ 0 (17) где Mtot — полная масса системы. Размеры и полные массы моделей для различных β приведены в таблице 1. β 0.0 −0.1 −0.2 −0.3 −0.4 −0.5 R 4.38 · 10−1 3.51 · 10−1 2.63 · 10−1 1.78 · 10−1 1.03 · 10−1 4.63 · 10−2 Mtot 2.40 · 10−1 1.61 · 10−1 9.89 · 10−2 5.33 · 10−2 2.36 · 10−2 7.57 · 10−3 Таблица 1: Размеры R и массы Mtot моделей при различных параметрах β. 2.2. Начальное распределение по координатам По определению µ(r̃) — монотонно возрастающая функция на отрезке [0, 1]. Поэтому существует обратная функция r̃(µ), и для генерации радиальной составляющей координаты частиц можно использовать метод обратного преобразования. Он состоит из двух шагов. 1. Задаем случайное число µ, равномерно распределенное на отрезке [0, 1]. 2. Находим расстояние r̃ по построенной зависимости r̃(µ). Угловые координаты ϕ распределяются равномерно на промежутке [0, 2π). 10

2.3. Начальное распределение по скоростям Введем переменные: x=p Vr 2(E0 − Φ(r)) ,y = Vφ . 2(E0 − Φ(r)) (18) Запишем фазовую плотность (3) в новых переменных, отбросив все, что от них не зависит: β f (x, y)dxdy = y β (1 − x2 − y 2 )− 2 dxdy. (19) C помощью плотности функции распределения (19) построим распределения частиц по скоростям (Vr , Vφ ). Плотность распределения по x получается интегрированием по всем возможным значениям y: "Z f (x) = √ # 1−x2 − β2 y β (1 − x2 − y 2 ) dy . (20) 0 Сама функция распределения F (x) будет иметь вид: F (x) = p 1 1 · (arcsin(x) + x 1 − x2 ) + . π 2 (21) Функция F (x) строго монотонна и возрастает на отрезке [−1, 1]. Это позволяет использовать метод обратного преобразования для нахождения радиальной составляющей скорости частиц. В отличие от распределения частиц по радиусу, здесь приходится численно решать уравнение вида: F (x) − p = 0, (22) где p — случайное число от нуля до единицы. Для получения распределения по y используем метод “отбора-отказа”. Для обычной модификации этого метода требуется выбрать такую константу M , что для любых x, y f (x, y) < M . Функция f (x, y) → ∞ при y → 0 при отрицательных β. Следовательно, M → ∞. Поскольку в реальной системе не может быть сингулярностей, для построения можно использовать конечное, но большое значение M . Оценим величину ошибки, связанную 11

с таким выбором. Для этого найдем, какое минимальное значение y может быть выбрано при фиксированном M . Проинтегрируем неравенство f (x, y) > M по x: f (y)dy = y β "Z √ 2 1−y − √ # β 2 (1 − x2 − y 2 ) dx dy > M · 1−y 2 Z √1−y2 − √ dxdy. (23) 1−y 2 После взятия интегралов по x приходим к выражению √ πΓ(1 − β2 ) Γ( 23 − β2 ) y β (1 − y 2 ) −β 2 > 2M. (24) Нас интересует окрестность y ≈ 0, поэтому 1 − y 2 ≈ 1. Для удобства √ введем обозначение M ? = 21 M · Γ( 23 − β2 ))( πΓ(1 − β2 ))−1 . Неравенство (24) преобразуется к виду 1 (25) y < (M ? ) β . В работе использовалось M ? = 103 и β = 0.0, −0.1, −0.2, −0.3, −0.4, −0.5. Как было отмечено выше, рассматриваемая проблема возникает только для отрицательных β. Наибольшая величина ошибки получается при β = −0.5. В этом случае y < 10−6 оказываются “отрезанными”. Эта величина оказывается порядка ошибки интегрирования уравнений движения частиц. Кроме того, y ≈ 0 соответствует угловой скорости Vϕ ≈ 0, то есть радиальным орбитам. Такое “обрезание” функции распределения может только препятствовать образованию баров из-за неустойчивости радиальных орбит. Как будет показано ниже, модели при любых β оказываются неустойчивыми, поэтому использование метода “отбора-отказа” оказывается корректным в данном случае. Опишем подробно реализацию алгоритма выбора y. 1. Выберем p случайным образом на отрезке [0, 1]. Вычислим P = p · M . 2. Найдем x из уравнения (22). √ 3. Выберем y случайным образом на отрезке [0, 1 − x2 ]. 4. Вычислим значение f (x, y) по формуле (19). 12

5. Если f (x, y) > P , то выбранная пара (x, y) подходит. Если f (x, y) < P пара (x, y) отбрасывается, и алгоритм повторяется, начиная с шага 3. От переменных (x, y) перейдем к (Vr , Vφ ): s Vr = s Vφ = 2B(r̃) · x, B(0) (26) 2B(r̃) · y. B(0) (27) Массы частиц предполагаются одинаковыми, при этом масса одной частицы Mtot , (28) m= N где Mtot — полная масса системы, N — число частиц в построенной модели. Полученные распределения по координатам для различных β приведены на рис. 1. Для проверки правильности построенных моделей были построены профили поверхностной плотности (см. рис. 2(a)). Поверхностная плотность рассчитывается двумя способами: “теоретическим” — по формуле (15) и “практическим” — суммированием масс частиц по радиальным слоям. Из рис. 2(a) видно, что полученные двумя различными способами зависимости хорошо совпадают. Для проверки построенных моделей по скоростям в соответствующих моделях был рассчитан параметр β. Из рис. 2(b) видно, что значения параметра β в слоях начинают отклоняться от выбранного значения. По-видимому, это обусловлено конечным числом частиц в слоях. 13

β = 0.0 β = −0.1 β = −0.2 β = −0.3 β = −0.4 β = −0.5 Рис. 1: Примеры получившихся моделей для различных β при N = 10000. 14

a) b) Рис. 2: Профиль поверхностной плотности системы ρ (a) и параметра анизотропии β (b) при N = 105 , β = 0.0, −0.1, −0.2, −0.3, −0.4, −0.5. a) Различными значками отложены значения плотности, рассчитанные суммированием масс частиц по радиальным слоям системы. Штрихованными кривыми показана зависимость (15) для соответствующих β. b) Точками изображены получившиеся значения параметра β в каждом слое системы, рассчитанные по формуле (4). Из-за дискретности значения β в слоях “гуляют” вокруг идеальных значений. 3. Расчет времен релаксации и оптимальных параметров сглаживания моделей 3.1. Численная релаксация в задаче N тел При исследовании систем с помощью численного моделирования возникает несколько проблем. Первая — это ограниченность вычислительных мощностей. В реальных галактиках число звезд 109 −1011 . Расчет на одном современном компьютере эволюции такого числа частиц занимал бы годы. Поэтому систему аппроксимируют меньшим числом частиц: 105 − 107 . Вторая проблема, вырастающая из первой, — это уменьшение времени релаксации системы при уменьшении количества частиц. По Чандрасекару (1939) время релаксации — это время, за которое звезда существенным образом 15

изменит свои характеристики (энергию или траекторию). За промежуток времени, равный по порядку величины времени релаксации, система “забывет” начальные условия. Такой эффект возникает из-за дискретности распределения материи в галактике. Время релаксации Trel для трехмерных систем растет с числом частиц N в системе (Бинни и Тремейн, 2008): Trel ∝ N . log(N ) (29) Галактики имеют время релаксации больше хаббловского. В численных расчетах, использующих число частиц на порядки меньше, чем в реальных галактиках, приходится делать оценку времена релаксации и учитывать ее при определении времени интегрирования модели. К сожалению, оценка (29) оказывается совсем не верна в двумерных плоских моделях. Для двумерных моделей время релаксации порядка времени пересечения системы и не зависит от числа частиц (Райбики, 1972). Это связано с тем, что в двухмерных системах роль близких сильных столкновений становится такой же, как и роль далеких. В трехмерных же системах основной вклад в отклонение траектории частицы вносят далекие слабые взаимодействия. Обычно при моделировании звездных систем модифицируют закон гравитации: 1 1 F (r) ∝ 2 7−→ F (r) ∝ 2 (30) r r + ε2 где ε — параметр сглаживания системы. Это делается по нескольким причинам. Сглаживание по закону (30) позволяет избежать расходимости вычислений при тесном сближении частиц и уменьшить флуктуации плотности, возникающие естественным образом при создании модели. Сглаживание двумерной модели позволяет рассматривать двумерный плоский диск, как трехмерный диск, имеющий конечную маленькую толщину ε. Для системы конечной, но малой толщины ε возникает зависимость времени релаксации от числа частиц (Райбики, 1972, Вайт, 1988): Trel ∝ N. 16 (31)

3.2. Расчет оптимальных параметров сглаживания По причинам, описанным в предыдущем разделе, двумерную модель необходимо сглаживать. Выбирать параметр сглаживания нужно так, чтобы при изменении гравитационного закона (30), система оставалась максимальна близка к своему непрерывному аналогу. Это означает, что величина (Мерит, 1996): N 1 X ~ |Fi − F~true (~ri )|2 (32) ASE = N i=1 должна быть минимальна. ASE (Average Square Error) - это величина среднего квадрата отклонения теоретической силы F~true от силы F~i , действующей на i-ую частицу в построенной модели с параметром сглаживания ε. Сила F~true вычисляется как градиент потенциала: Ftrue (ri ) = 1 dŨ (r̃i ) . R dr̃ (33) Сила F~i находится простым суммированием по всем частицам: F~i = Gm N X ~rj − ~ri , 2 + |~ 2 3/2 (ε r r i −~ j| ) j=1,j6=i (34) Величину параметра сглаживания будем измерять также в единицах радиуса ε̃ = ε/R. На рис. 3 изображена типичная зависимость ASE(ε̃) на примере моделей с β = −0.2 и β = −0.4. Оценку величины ASE можно уточнить, если усреднить ее по большому числу моделей: K 1 X ASEi , (35) M ASE = K i=1 где K - число моделей. Для нахождения величины (35) использовалось K = 100. На рис. 4 приведены рассчитанные зависимости оптимального параметра сглаживания ε̃opt от числа частиц N в диапазоне 300 ≤ N ≤ 30000 для различных β. С увеличением числа частиц оптимальный параметр сглаживания уменьшается. Зависимость ε̃ хорошо аппроксимируются степенной функцией: ε̃opt = A · N B , (36) 17

Рис. 3: Типичная зависимость ASE(ε̃). Стрелочками отмечены положение минимумов отклонений. По ним находятся оптимальные параметры сглаживания модели. где A и B зависят от параметра анизотропии β. Для моделей с большими N оптимальные значения параметра сглаживания вычисляются по найденным зависимостям (36). 3.3. Оценка времен релаксации В работах Райбики (1972) и Вайта (1988) была получена теоретическая оценка времени релаксации двумерного звездного диска с параметром сглаживания ε: σ3ε , (37) Trel ≈ 5G2 ρm где σ — дисперсия скоростей звезд, ρ — поверхностная плотность, m — масса одной частицы в системе. Поскольку m ∝ N −1 при фиксированной массе системы, то время релаксации Trel ∝ N . Оценка (37) корректна, когда поверхностная плотность не имеет резких пиков. Однако, в обобщеннополитропных моделях поверхностная плотность растет к центру и в центре имеет особенность. Поэтому корректную оценку времени релаксации приходится искать численно. Простой способ получить оценку времени релаксации — напрямую проследить, как отклоняются частицы в самосогласованном гравитационном поле (Атанасула и др., 2001). Вначале все частицы системы фиксируются. 18

Рис. 4: Зависимость оптимального параметра сглаживания ε̃opt от числа частиц N в диапазоне 300 ≤ N ≤ 30000 для различных параметров β. Штрихованными прямыми обозначены наилучшие аппроксимации зависимостей степенной функцией, ε̃opt (N ) = A · N B , для соответствующих β. Затем с какого-то характерного радиуса rh вдоль радиуса “запускаются” пробные частицы. Пробная частица будет двигаться в “замороженном” потенциале фиксированных частиц системы и отклоняться от первоначальной строго радиальной орбиты. При этом время релаксации можно оценить как hτi i Trel = , (38) hsin2 αi где τi — время движения на вычисленном отрезке траектории, α — угол между векторами начальной и конечной скорости частицы, угловые скобки обозначают усреднение по всем запускаемым частицам. С одной стороны, rh нужно выбрать так, чтобы все “запускаемые” частицы пролетали через центр, с другой — желательно, чтобы отклонения вектора скорости на каждом шаге были небольшими (в этом случае на каж19

дом шаге sin α ≈ α), иначе формула (38) не корректна (Стендишь и Акснес, 1969). В предварительных расчетах было установлено, что для этого подходит значение радиуса rh такое, что круг радиуса rh содержит половину массы системы. Начальное положение пробных частиц на окружности радиуса rh определяется случайным образом. Также было установлено, что отклонения частицы на каждом шаге зависят от ее скорости. Чтобы отклонения были маленькими, начальную скорость частиц следует задавать большой — порядка скорости убегания с радиуса rh : s vrh = 2B(r̃h ) . B(0) (39) Уравнения движения частицы: d~r = ~v , dt (40) X ~ri − ~r d~v =m , 2 + ε2 ) 32 dt ((~ r − ~ r ) i i=1 (41) N интегрируются с помощью схемы Рунге-Кутты четвертого порядка: ~k1 = dt · ~v (~rn , tn ), ~l1 = dt · ~a(~rn , tn ), ~k2 = dt · ~v (~rn + ~k1 /2, tn+1/2 ), ~l2 = dt · ~a(~rn + ~k1 /2, tn+1/2 ), ~k3 = dt · ~v (~rn + ~k2 /2, tn+1/2 ), ~l3 = dt · ~a(~rn + ~k2 /2, tn+1/2 ), ~k4 = dt · ~v (~rn + ~k3 , tn+1 ), ~l4 = dt · ~a(~rn + ~k3 , tn+1 ), 1 ~rn+1 = rn + (~k1 + 2~k2 + 2~k3 + ~k4 ), 6 1 ~vn+1 = vn + (~l1 + 2~l2 + 2~l3 + ~l4 ), 6 где n — номер шага интегрирования, dt — шаг интегрирования. Для уменьшения ошибки интегрирования в системе с оптимальным па20

раметром сглаживания εopt шаг интегрирования желательно выбрать так, чтобы частица проходила расстояние ε за несколько шагов интегрирования (Родионов, Сотникова, 2005): dt = εopt . 5vrh (42) Траектории пробных частиц прослеживаются до тех пор, пока частицы снова не достигают расстояния rh от центра. Результаты расчетов времени релаксации для различных значений числа частиц в системе N и параметров анизотропии β приведены на рис. 5. Как и ожидалось, с увеличением числа частиц время релаксации растет. Рис. 5: Зависимость времени релаксации Trel от числа частиц в системе N при различных параметрах β, рассчитанная по формуле (38). Штрихованными прямыми обозначены наилучшие аппроксимации зависимостей степенной функцией Trel (N ) = A · N B для соответствующих β. Полученные времена релаксации затем были использованы при числен21

ном расчете эволюции модели. Для оценки времени релаксации в случае большего числа частиц использовалась аппроксимация, аналогичная (36) с подходящими значениями A и B. 4. Численное моделирование обобщенно-политропных моделей звездных дисков 4.1. Нормировка расстояний, скоростей и времени на размер системы Определив время релаксации и выбрав оптимальный параметр сглаживания, можно переходить к численному моделированию системы. Для удобства вычислений перейдем к новым обозначениям: ~r → V~ ε t ~r , ε → , V~ → √ , t → 3 . R R R R2 (43) Ускорение частицы, вызываемой гравитационным притяжением других частиц в системе, в новых обозначениях записывается в виде: N X d2~ri ~rij = m . 2 + ε2 ) 32 dt2 (r ij j=1,j6=i (44) Замена переменных (43) помогает избежать ошибок в вычислениях при малых размерах системы R и упрощает отображение результатов расчетов. 4.2. Решение задачи N тел. Пакет N EM O Для моделирования использовались стандартные средства пакета N EM O (Теубен, 1995) и его расширения gyrf alcON (Дехнен, 2002). Пакет N EM O позволяет численно решать уравнения движения частиц и быстро находить силы гравитационного притяжения, действующие между частицами. “Быстро” в данном случае означает, что для расчета сил в N EM O используется алгоритм, имеющий трудоемкость меньше, чем алгоритм простого 22

суммирования. Имеется ввиду treecode-алгоритм (Барнс, Хат, 1986), имеющий сложность O(N log N ), и его комбинация с быстрым разложением по мультиполям (Грингард, Роклин, 1997), с трудоемкостью O(N ) (Дэнен, 2002). Ключевая идея этих алгоритмов — представление всей системы в виде иерархического дерева. Если частица расположена далеко от какой-то ячейки-“ветки”, то, вместо рассмотрения взаимодействия частицы с каждой частицей в ячейке, можно рассматривать взаимодействие частицы с ячейкой в целом. Параметр удаленности θ ячейки от частицы определяется как: d (45) θ= , r где d — размер ячейки, а r — расстояние между частицей и ячейкой. Фактически, θ является угловым размером ячейки, видимой из того места, где находится частица. Если угловой размер ячейки слишком большой, т. е. больше некого θmin , то рассматривается взаимодействие частицы и подъячеек данной ячейки. Экспериментально проверено, что θmin = 0.5 дает оптимальное соотношение по скорости и точности счета. Шаг интегрирования выбирался в виде ∆t = 2−n (n — натуральное число). Показатель степени n определяется из неравенства: 2−n < dt < 2−n+1 , (46) где dt — оценка шага интегрирования (42) для данной модели. Время интегрирования Tint определяется по рассчитанному для модели времени релаксации: Tint = Trel /10. Предполагается, что на временах, на порядок меньших времени релаксации, численная релаксация еще не проявляется. 4.3. Неустойчивость обобщенно-политропных моделей звездных дисков В предыдущих работах по исследованию таких моделей (Поляченко и Поляченко, 1994, Кошкин, дипломная работа, 2006) приводилась граница неустойчивости по степени анизотропии β. Было получено, что при β ≤ −0.5 модели неустойчивы, а при β > −0.5 — устойчивы. Эти результаты были получены при расчетах с неудовлетворительно малым числом 23

частиц, N ∼ 103 . Наши расчеты показали, что при любой степени вытянутости орбит β модели оказываются неустойчивыми, и образуется бар. На рис. 6 —10 представлены результаты расчетов эволюции обобщеннополитропных моделей звездных дисков для различных значений параметра анизотропии β = 0.0, −0.1, −0.2, −0.3, −0.4. Для β ≤ −0.5 неустойчивость моделей была исследована в других работах. Рис. 6: Положение частиц в квадрате −1 ≤ x ≤ 1, −1 ≤ y ≤ 1 при β = 0.0, N = 105 в различные моменты времени. 24

Рис. 7: Положение частиц в квадрате −1 ≤ x ≤ 1, −1 ≤ y ≤ 1 при β = −0.1, N = 105 в различные моменты времени. 25

Рис. 8: Положение частиц в квадрате −1 ≤ x ≤ 1, −1 ≤ y ≤ 1 при β = −0.2, N = 105 в различные моменты времени. В начальные моменты времени заметны “спицы”. 26

Рис. 9: Положение частиц в квадрате −1 ≤ x ≤ 1, −1 ≤ y ≤ 1 при β = −0.3, N = 105 в различные моменты времени. В начальные моменты времени заметны “спицы”. 27

Рис. 10: Положение частиц в квадрате −1 ≤ x ≤ 1, −1 ≤ y ≤ 1 при β = −0.4, N = 105 в различные моменты времени. 28

4.4. Расчет скорости узора и амплитуды гармоники, соответствующей бару Преобразование Фурье поверхностной плотности в координатах r и ϕ позволяет получить динамические характеристики бара: Z ∞ Z 2π ρ(r, ϕ, t) · exp {i [p ln(r) + mϕ]} rdϕdr. ρ̄(m, p, t) = 0 (47) 0 Бару соответствует гармоника с азимутальным числом m = 2 и радиальным числом p = 0. На практике плотность частиц выражается через дельта-функции. Поэтому непрерывное преобразование Фурье (47) вычисляется как сумма: N 1 X ρ̄(p, m, t) = exp {i [p ln(rj ) + mϕj ]} . N j=1 (48) На рис. 11 представлены результаты расчетов по формуле (48) при m = 2 и p ∈ [−15, 15] для моделей с β = 0.0, −0.1, −0.2, −0.3 для моментов времени, отображенных на рис. 6— 9, соответственно. В начальные моменты времени бар отсутствует, и максимум |ρ̄(p, m, t)| смещен относительно p = 0. Затем системы эволюционируют, образуется бар, и максимум смещается к p = 0. Перейдем от Фурье-образа плотности ρ̄(m, p, t) к скорости узора бара Ωp (t). Преобразование Фурье от постоянной волны плотности можно записать как (Селвуд и Атанасула, 1986): ρ̄(m, p, t) = a(p) · exp(−iωt), (49) где ω - комплексная частота собственных колебаний волны плотности. Фаза этих колебаний α(t) = −Re(ω)t определяется из соотношений:  Im(ρ̄(2, 0, t))   , tan α(t) = Re(ρ̄(2, 0, t)) Im(ρ̄(2, 0, t))   . sin α(t) = |ρ̄(2, 0, t))| 29 (50)

Рис. 11: Логарифм амплитуды Фурье-образа поверхностной плотности ln |ρ̄(p, m = 2, t)| для моделей с различными β в моменты времени, отображенные на рис. 6—9, соответственно. Жирной линией обозначаются профили Фурье-образа плотности, соответствующие конечным моментам времени. Видно, что системы эволюционируют, и максимум Фурье-образа плотности смещается к p = 0, что соответствует образованию бара. Скорость узора Ωp при произвольном m вычисляется по формуле: Ωp = 1 dα(t) . m dt (51) При численном моделировании имеется дискретный набора значений фазы 30

α. Для нахождения скорости узора бара Ωp в произвольный момент времени t фазу колебаний α(t) аппроксимируем линейной функцией на временном промежутке [t, t + ∆T ]. Величина m · Ωp будет равна угловому коэффициенту аппроксимирующей прямой. В данной работе использовалось временное окно ∆T = 35. Такая ширина временного окна приблизительно равна 1 − 2 периодам обращения бара в рассматриваемых моделях. На рис. 12(a) представлена типичная для гармоники бара зависимость α(t). На рис. 12(b) представлены рассчитанные скорости узора Ωp для различных значений параметра β. Для возможности сравнения скорости узора бара в моделях с различным β время и скорость узора преобразуются к исходным единицам, не нормированным на размер системы (обратно преобразованию (43)). Из рисунка 12(b) видно, что с уменьшением параметра β скорость узора становится меньше. 5. 5.1. Обсуждение результатов Неустойчивость радиальных орбит. Необходимые условия ее реализации Действительно ли бары в политропных моделях звездных дисков образуются благодаря неустойчивости радиальных орбит? Если это так, то радиальные, похожие на “спицы” орбиты должны найтись в системе. Кроме того, прецессирующие “спицы” должны каким-то образом изменить свою скорость прецессии Ωpr так, чтобы в итоге образовать одну большую “спицу”. Это возможно, когда орбиты прецессируют медленно: Ωpr Ω, (52) где Ω — угловая скорость вращения частицы. При этом должна существовать определенная зависимость скорости прецессии от углового момента для орбит частиц (Линден-Белл, 1979), а именно: ∂Ωpr ∂L > 0. Jf 31 (53)

a) b) Рис. 12: Скорость узора бара. a) Типичная зависимость фазы стоячей волны α(t), соответствующей бару (сплошная линия) на примере модели β = 0.0 в промежутке времени 25 ≤ T ≤ 100. Зависимость за один период хорошо аппроксимируется линейной функцией. Для примера показана аппроксимирующая прямая на момент 1 ∆α(t) времени T = 50. Скорости узора Ωp = соответствует половина ко2 ∆t эффициента наклона аппроксимирующей прямой. b) Рассчитанная скорость узора для различных β при N = 105 . Вертикальными отрезками отмечены моменты времени, соответствующие сформировавшемуся бару. Видно, что c уменьшением параметра β скорость узора бара становится меньше. Здесь Ω — угловая частота, L — угловой момент, Jf — действие, определяемое интегралом: Z 1 Jf = p~ · d~q, (54) 2π где, как обычно, p~ и ~q — обобщенные импульсы и координаты. Если орбиты прецессируют медленно, то величина Jf является адиабатическим инвариантом. Начальные значения Ωpr , Ω и зависимость Ωpr (L) можно найти, если зафиксировать распределение потенциала системы на начальный момент времени и проследить за тем, какие орбиты описывают частицы в заданном потенциале (11). При этом уравнения движения частиц приходится интегрировать численно. 32

5.2. Расчет траекторий частиц в фиксированном начальном потенциале. Определение скоростей прецессии частиц, эксцентриситетов их орбит и адиабатических инвариантов 5.2.1 Интегрирование уравнений движения Для интегрирования уравнений движения частиц использовалась схема с перешагиванием второго порядка точности: [k+1/2] ~ri [k+1] ~vi [k−1/2] = ~ri [k] [k] + δt~vi , (55) [k+1/2] (56) = ~vi + δt~ai , где δt — шаг интегрирования, k — номер шага, индекс i обозначает номер частицы, а ~ri , ~vi , ~ai соответственно радиус-вектор, скорость и ускорение частицы. Шаг интегрирования выбирался согласно (42). Ускорение частицы в фиксированном начальном потенциале [k+1/2] [k+1/2] ai dU (ri = dr ) . (57) Для сохранения второго порядка точности на первом шаге интегрирования положение частицы вычисляется по формуле: [1/2] ~ri = [0] ~ri δt [0] δt2 ~[0] + ~vi + a . 2 8 i (58) На рис. 13 и 14 представлены примеры описываемых частицами орбит различной степени вытянутости. В ходе интегрирования появляются частицы, которые удаляются от центра на расстояние большее размера системы, ri > R. Формально скорость частиц на границе системы (при ri = R) должна быть нулевой. Если бы мы рассматривали непрерывное движение, то частицы могли бы подлетать к границе сколько угодно близко. На практике, из-за дискретности шага по времени частицы могут “перелетать” границу. В наших расчетах таких частиц оказалось приблизительно 0.01% от их общего числа. Энергия частиц сохранялась с точностью до шестого знака после запятой. 33

Рис. 13: Орбиты различной степени вытянутости e в покоящейся системе отсчета (слева) и в системе отсчета, вращающейся со скоростью прецессии (справа) для модели c β = −0.2. 5.2.2 Расчет скорости прецессии. Смещение апоцентра орбиты Частицы, движущиеся в потенциале (11), будут прецессировать, их орбита будет не замкнута (cм. рис. 13 и 14, слева). Скорость прецессии и соответствующую ей эпициклическую частоту можно найти по смещению апоцентра орбиты: 2π ∆ψ Ωpr = , κ= , (59) ∆t ∆t 34

Рис. 14: Орбиты различной степени вытянутости e в покоящейся системе отсчета и в системе отсчета, вращающейся со скоростью прецессии (справа) для модели c β = −0.2. 35

где ∆ψ — угловое смещение апоцентра, ∆t — время, прошедшее между последовательным прохождением апоцентра орбиты. Положение апоцентра орбит rapo,i находится в ходе интегрирования из условий: [k−3/2] ri [k−1/2] < ri [k−1/2] , ri [k+1/2] > ri [k−1/2] , rapo,i = ri . (60) На рис. 15 показана зависимость скорость прецессии от больший полуоси орбиты. Из-за различной степени вытянутости орбит зависимость не однозначна (Страк и Линдблад, 2015). На рис. 13 и 14 справа представлены орбиты в системе отсчета, вращающейся со скоростью прецессии. Естественно, они замкнуты. 5.2.3 Расчет адиабатического инварианта орбиты Перейдем в систему отсчета, вращающуюся со скоростью прецессии орбиты. В этой системе отсчета за одну радиальную осцилляцию частица совершает половину оборота вокруг центра системы. В этом случае инвариант орбиты Jf можно представить в виде: Jf = Jr + L , 2 (61) ṙdr. (62) где Jr — радиальное действие: 1 Jr = 2π Z Интеграл (62) берется численно по двум радиальным осцилляциям частицы. 5.2.4 Расчет эксцентриситета орбиты Эксцентриситет орбиты e определяется через радиальное действие и угловой момент (Линден-Белл, 1979) из уравнения: 2 1 1 − e2 1 Jf √ = + . 2 1 − e2 2 Lz 36 (63)

Рис. 15: Распределение скорости прецессии Ωpr частиц в зависимости от большой полуоси орбиты a при различных значениях параметра анизотропии β. 37

5.3. Эксцентричность орбит в обобщенно-политропных моделях звездных дисков В настоящей работе были рассчитаны эксцентриситеты орбит частиц в обобщенно-политропных моделях звездных дисков при β = 0.0, −0.1, −0.2, −0.3 и N = 105 . На рис. 16 показаны диаграммы “количество частиц — эксцентриситет”. Из рисунка видно, что большинство орбит в обобщенно-политропных моделях звездных дисков являются сильно вытянутыми (e ≥ 0.9). Даже для изотропной модели с β = 0.0 орбит с большими эксцентриситетами, e ≥ 0.9, больше половины. Этот результат прямо противоречит тому, что утверждалось в работе Поляченко и Поляченко (1994). Таким образом, в рассмотренных моделях имеется большое количество “спиц”, которые могут быть подвержены неустойчивости радиальных орбит. На рис. 17 представлены изолинии эксцентриситета частиц в переменных скорость прецессии Ωpr — большая полуось орбиты a. Из рисунка видно, что при фиксированной большой полуоси наибольшую скорость прецессии имеют орбиты, близкие к круговым. Частицы с большими эксцентриситетами, напротив, прецессируют медленнее всего. 5.4. Зависимость скорости прецесии от углового момента при фиксированном адиабатическом инварианте Необходимым условием для неустойчивости радиальных орбит является специфическая зависимость скорости прецессии от углового момента при постоянном инварианте Jf : ∂Ωpr ∂L > 0. (64) Jf На основе значений Ωpr , Lz и Jf для каждой частицы были построены изолинии скорости прецессии в координатах действие Jf — угловой момент Lz для различных значений параметра анизотропии β при N = 105 (см. рис. 18). Изолинии строятся следующим образом: вначале все величины Ωpr , Lz и Jf сортируются по скорости прецессии и для каждой изолинии 38

Рис. 16: Диаграмма “количество частиц - эксцентриситет” при различных значениях параметра β. На графиках указано сколько частиц имеют e > 0.9 и e > 0.99. выбирается какая-то стартовая частица. Далее на график наносятся следующие за ней 103 элементов из массива полученных данных Ωpr . Около каждой изолинии показано среднее значение скорости прецессии для выбранных частиц. Пунктирная прямая Jf = Lz /2 соответствует границе круговых орбит. При постоянном значении действия Jf частицы могут двигаться только параллельно горизонтальной оси момента Lz . Из графиков на рис. 18 видно, что при таком движении выполняется условие (64). Это означает, что для всех частиц в системе выполнено условие, связывающее угловой момент и скорости прецессии, необходимое для неустойчивости 39

Рис. 17: Изолинии эксцентриситета орбит e в переменных скорость прецессии Ωpr — большая полуось орбиты a. Видно, что, меньшим значениям эксцентриситета соответствует большая скорость прецессии. радиальных орбит. Более того, по крайней мере, на рис. 9 и 10 на ранних стадиях эволюции модели, еще до образования бара, отчетливо видно множество радиальных лучей, которые могут быть следствием “слипания” спицеобразных орбит. 40

Рис. 18: Изолинии скорости прецессии Ωpr для различных значений параметра анизотропии β при N = 105 . По оси ординат откладывается величина действия Jf , по оси абсцисс — угловой момент Lz . Пунктирная прямая Jf = Lz /2 соответствует круговым орбитам (Jr ≈ 0). 5.5. Сравнение скорости узора бара и скоростей прецессии частиц Сравним скорость бара, образующегося в модели, со скоростями прецессии частиц. Естественно ожидать, что если речь идет о прецессирующих “спицах”, то скорость бара должная быть близка к скоростям прецессии соответствующих частиц. На рис. 19 представлено распределение частиц по скоростям прецессии Ωpr в зависимости от большой полуось орбиты a. 41

На этом же рисунке нанесены скорости узора бара для двух моментов времени: начального, соответствующего сформировавшемуся бару, и соответствующего середине времени счета. Для возможности сравнения моделей при различных значениях параметра анизотропии β скорости прецессии и узора в разных моделях преобразуются к исходным (обратно преобразованию (43)). Из графиков на рис. 19 видно, что скорость образующегося бара оказывается больше скорости прецессии почти всех частиц. Исключение составляют частицы, движущиеся только в центральной области (большая полуось a ≤ 0.1). 6. Выводы Основной результат, полученный в данной работе — неустойчивость моделей обобщенно-политропных звездных дисков при любом параметре анизотропии системы β. Данный результат противоречит результату, полученном в работе Поляченко и Поляченко (1994). В отличие от предыдущих работ, посвященных изучению таких моделей, сделан акцент на корректной оценке параметров моделирования: оптимальном параметре сглаживания системы εopt и времени релаксации Trel при различных значениях параметра анизотропии β. В работе Поляченко и Поляченко (1994) получена граница по неустойчивости для параметра анизотропии β ≈ −0.5, однако при этом совсем не учитывалась возможная численная релаксация моделей. В работе было выбрано малое количество частиц: N ∼ 103 . При таком числе частиц численная релаксации системы наступает еще до образования бара, поэтому есть основания сомневаться в корректности результата авторов. В дипломной работе Кошкина (2006) “оптимальные” параметры сглаживания εopt были оценены из общих соображений по порядку величины. Рассчитанные в данной работе значения параметров сглаживания оказались на порядок больше. В настоящей работе получены косвенные свидетельства того, что образование бара в подобных моделях обусловлено неустойчивостью радиальных орбит. Показано, что при любых значениях параметрах анизотропии β есть большое количество спецеобразных орбит, которые из-за специфи42

Рис. 19: Скорость прецессии Ωpr частиц (точками) в сравнении со скоростями бара (прямые) в различные моменты времени. Скорость узора оказывается больше скоростей прецессии частиц, за исключением центральной области (a ≤ 0.1). ческой зависимости скорости прецессии от углового момента могут “слипаться” по механизму Линден-Белла, образуя бар. Особый интерес представляет сравнение скорости узора образующегося бара с начальным распределением скоростей прецессии частиц. По неясным причинам скорость узора бара оказывается примерно в два раза больше средней скорости прецессии частиц в начальный момент времени. Это означает, что при развивающейся неустойчивости радиальных орбит должен существовать какой-то механизм разгона бара. Этого факта не было 43

замечено авторами предыдущих работ, посвященных изучению неустойчивости радиальных орбит. В дальнейшем планируется более детальное исследование данного вопроса. 7. Благодарности Я благодарен моему научному руководителю Наталье Яковлевне Сотниковой за оказанной в работе поддержку, плодотворные обсуждения и советы. Я хотел бы выразить благодарность Петру Александровичу Тараканову за представленные технические возможности для расчетов. Список литературы 1. Атанасула и др. (E. Athanassoula at al.), A&A 376, 1135 (2001). 2. Барнес, Хат (J. Barnes, P. Hut), Nature 324, 446 (1986). 3. Бинни, Тремейн (J. Binney and S. Tremaine), Galactic Dynamics: Second Edition, Princeton University Press, 2008. 4. Бисноватый-Коган Г.С., Письма в Астрон. журн. 1, 177 (1975). 5. Вайнберг, Катс (M. Weinberg, N. Katz), MNRAS 375, 425 (2007). 6. Вайт (R. White), Astrophys. J. 330, 26 (1988). 7. Грингард, Роклин (L. Greengard, V. Rokhlin), Journal of Computational Physics 135, 280 (1997). 8. Дэнен (W. Dehnen), Journal of Computational Physics 179, 27 (2002). 9. Линден-Белл (D. Lynden-Bell), MNRAS 187, 101 (1979). 10. Линден-Белл, Калнайс (D. Lynden-Bell, A. Kalnajs), MNRAS 157, 1 (1972). 11. Мерритт (D. Merrit), Astronomical J. 111, 2462 (1996). 44

12. Поляченко В.Л., Поляченко Е.Л., Письма в Анстрон. журн. 20, 416 (1994). 13. Понтзен и др. (A. Pontzer at al.), MNRAS 451, 1366 (2015). 14. Райбики (G. Rybicki), Relaxation Times in Strictly Disk Systems (Ed. M. Lecar, IAU Colloq. 10: Gravitational N-Body Problem, 22), 1972. 15. Селвуд (J. Sellwood), Reviews of Modern Physics 86, 1 (2014). 16. Селвуд, Атанасула (J. Sellwood, E. Athanassoula), MNRAS 211, 195 (1986). 17. Родионов С.А., Сотникова Н.Я., Письма в Астрон. журн. 49, 470 (2005). 18. Стендишь, Акснес (Jr. Standish and Aksnes K.), Astrophys. J. 158, 519 (1969). 19. Страк (C. Struck), MNRAS 450, 2217 (2015). 20. Теубен (P. Teuben), The Stellar Dynamics Toolbox NEMO (Ed. R. Shaw, H. Payne, and J. Hayes, Astronomical Data Analysis Software and Systems IV, 398), 1995. 21. Чандрасекар (S. Chandrasekhar), Principles of stellar dynamics, The University of Chicago press, 1942. 45

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв