Использование модифицированного метода Шварцшильда для построения фазовой модели Галактики Текст научной статьи по специальности «Математика»

АСТРОНОМИЯ

УДК 524.3/4-32

ИСПОЛЬЗОВАНИЕ

МОДИФИЦИРОВАННОГО МЕТОДА ШВАРЦШИЛЬДА ДЛЯ ПОСТРОЕНИЯ ФАЗОВОЙ МОДЕЛИ ГАЛАКТИКИ*

А. А. Башаков

С.-Петербургский государственный университет, аспирант, bashakov@mail.ru

1. Введение

Основными трудностями при исследовании динамики звездных систем являются нахождение фазовых моделей и построение соответствующих дискретных систем N тел, использующихся для изучения эволюции моделей. Система N тел, моделирующая звездную систему, должна удовлетворять нескольким критериям. Один из основных критериев соответствие распределения масс в построенной модели гравитационному потенциалу исследуемой системы. В нашей работе для получения фазовых самосогласованных моделей используется модифицированный метод Шварцшильда (оригинальный метод был опубликован в 1979 году [1]) с последующим переходом к дискретной модели. С применением модифицированного метода Шварцшильда нами ранее были построены модель Галактики, основанная на двухкомпонентном потенциале Кутузова— Осипкова [2], и тестовая модель со сферическим потенциалом Пламмера [3], в которой был проверен метод перехода к дискретному представлению для исследования устойчивости методами задачи N тел. Описание модифицированного метода и полученные результаты опубликованы в статьях [4, 5].

В данной работе представлен новый алгоритм нахождения весов орбит (необходимых для построения фазовой модели), основанный на методе М^е^о-Меазиге (М2М) [6, 7], который оказался более эффективным, чем применяемые ранее. С его помощью строится фазовая, а затем дискретная модель звездной системы, основанная на потенциале Флинна и др. для Галактики [8]. Результаты моделирования показывают, что этим методом можно решать поставленные выше задачи.

* Работа выполнена при финансовой поддержке РФФИ (грант №08-02-00361а) и Совета по грантам Президента РФ для государственной поддержки молодых российских ученых и ведущих научных школ (грант №НШ-1323.2008.2).

© А.А.Башаков, 2009

2. Модификации метода Ш^варцшильда

Наша реализация метода Шварцшильда для построения фазовых моделей звездных систем описана в статьях [4, 5]. В процессе работы метод несколько раз был существенно модифицирован, в частности, были изменены алгоритм разбиения исследуемого пространства на ячейки и алгоритм вычисления весов орбит. Все это позволило существенно повысить точность метода и увеличить его быстродействие.

3. Разбиение пространства, занимаемого моделью

Разбиение исследуемого пространства модели на ячейки прямоугольной формы, использовавшееся при построении модели Галактики Кутузова—Осипкова [4], обладает рядом недостатков. В данной работе были использованы цилиндрические ячейки, при этом размеры ячеек в радиальном направлении в части экспериментов экспоненциально уменьшались к центру модели, где звездная плотность возрастает. Эффективность применения ячеек с переменной толщиной видна из следующего примера: при использовании одинаковых прямоугольных ячеек с ребром в 1 кпк для заполнения модели с радиусом диска 30 кпк и «высотой» 5 кпк необходимо 4500 ячеек. В тоже время, если использовать 4500 цилиндрических ячеек с переменной толщиной в той же модели, то можно получить разрешение порядка 0.1 кпк в наиболее важной центральной части и 0.5 кпк на периферии. При использовании цилиндрических ячеек также была предусмотрена возможность деления их по азимуту, но так как нами рассматривались только осесимметричные модели, в этом не было необходимости.

4. Определение весов орбит Ы2Ы-методом

Для нахождения весов траекторий вместо использовавшегося в предыдущих работах симплекс-метода применялся модифицированный М2М-алгоритм [6, 7]. К числу преимуществ нового метода можно отнести возрастание скорости вычислений (что позволяет увеличить количество ячеек и орбит), постоянство исходной матрицы с временами нахождения частиц в ячейках (исключает эффекты накопления ошибок округления, возникающих в симплекс-методе, где с этой матрицей на каждом шаге производятся различные преобразования), возможность использования кинематических параметров как дополнительных условий (например, тангенциальной скорости), получение большего количества орбит в итоговом распределении (что делает его более сглаженным), поскольку метод назначает равные веса орбитам, имеющим одинаковое распределение времен/масс по ячейкам.

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

В методе используются следующие переменные: наблюдения У- — параметры исходного потенциала и параметры, полученные из наблюдений (так, в нашем случае это были массы в ячейках, найденные из уравнения Пуассона, и тангенциальные скорости в плоскости г = 0, полученные аналитически из потенциала Флинна с поправкой на дисперсию). Вычисленные значения у- —параметры модели, полученные объединением орбит с учетом весов: у- = ^ 1 т^х—, где т —искомые веса, х— —значение параметра, создаваемого *-й орбитой в ^’-й ячейке. Используются также отклонения — разности между наблюдениями и вычисленными значениями Д- = (У- — у- )/<г(У-), где ст(У-) — ошибки наблюдений.

Вычисление весов орбит производилось следующим образом: в начале все веса орбит

полагались равными 1/1, где I количество орбит. Затем они уточнялись в итерационном процессе. Каждая итерация состояла из трех шагов:

1) вычисление отклонений Д-;

2) изменение весов орбит:

3) нормирование весов: т* = т*/Ш, где Ш = ^ * т*.

Параметр е (итерационный шаг при уточнении весов) динамически изменялся в процессе работы. Общих рекомендаций по выбору его значения нет, поэтому использовался следующий подход: значение е возрастало от малых значений (10-3) в начале

(компенсируют уменьшение отклонений для поддержания скорости сходимости метода).

Так как изначально алгоритм разрабатывался для нахождения весов не орбит, а частиц, количество которых существенно превышает количество наблюдений, для устранения избыточности авторами [6] была введена энтропия $ = — ^ т* log(wj/т*), где

—начальные значения весов. В наших вариантах численного моделирования звездной системы количество орбит и количество наблюдений близки по величине, и использование энтропии излишне (коэффициент ц учета изменения энтропии брался нулевым).

По окончании работы алгоритма получены веса для каждой орбиты, причем, как показывают эксперименты, сумма абсолютных значений отклонений масс по всем ячейкам не превосходит 5-10% от общей массы построенной модели.

5. Особенности метода

После проведения первых численных экспериментов были выявлены особенности получающихся моделей: они имели большие дисперсии, порядка 80-120 км/с на расстоянии 10 кпк от центра модели, а кривая вращения на плоском участке соответствовала скорости порядка 150-180 км/с. Причина такого поведения кинематических характеристик оказалась связанной с начальной библиотекой орбит, из которой выбирался набор орбит для создания фазовой модели. Максимальное количество орбит и ячеек было ограничено сверху имеющимися вычислительными ресурсами. Нами использовались значения для них порядка 104. Таким количеством орбит необходимо было заполнить по возможности большую часть фазового пространства модели.

Для создания библиотеки орбит использовались распределения интеграла площадей J (радиус круговой орбиты Дс*г), интеграла энергии Е и энергии «вертикальных» осцилляций Ег. Особенно важным для быстро вращающихся подсистем (например, дисковая система нашей Галактики), является распределение орбит по интегралу площадей. Оно определяет положение (расстояние от оси вращения) круговых траекторий в экваториальной плоскости. Некруговые траектории с одинаковым значением интеграла площадей находятся в торообразной области, которая окружает круговую орбиту с данным значением J и размеры которой зависят от интеграла энергии Е. Чем больше энергия Е при данном интеграле площадей, тем сильнее траектория отличается от круговой. Энергия вертикальных колебаний Ег определяет, будет ли орбита расположена вблизи плоскости диска или отойдет от нее в область больших г.

работы метода (когда отклонения велики) до больших значений (102) к концу работы

Использование в первоначальных вариантах одинакового числа делений по каждому из этих параметров J, Е и Ег не дало удовлетворительного результата и вело к завышенной дисперсии скоростей в диске. Количество круговых орбит получалось порядка 20 при радиусе моделируемой области ^40 кпк и круговые орбиты отстояли друг от друга на расстоянии ^2 кпк (в окрестности Солнца). Ячейки же были расположены гораздо плотнее (размер ячейки ~0.2-0.3 кпк). Чтобы заполнить ячейки, в которые не попадали круговые и близкие к круговым орбиты, приходилось привлекать орбиты с большой энергией (более «размазанные» в пространстве). Естественно, при этом дисперсия скоростей в полученной модели возрастала. Чтобы устранить описанный эффект, круговые орбиты, вокруг которых формируются области возможных движений, должны быть расположены чаще. Расстояние между ними можно оценить следующим образом: дисперсия радиальных компонент скоростей звезд в Галактике в окрестности Солнца оценивается величиной порядка 30-40 км/с, это соответствует колебаниям радиуса орбиты на величину порядка 1 кпк, следовательно, и круговые орбиты должны быть расположены не реже.

Представленные в статье результаты были получены со следующими параметрами создания библиотеки орбит: 250 значений интеграла площадей (радиусы соответствующих круговых орбит распределены равномерно на интервале от 0 до 30 кпк), 6 значений интеграла энергии для каждого значения интеграла площадей (Дс*г), 6 значений энергии вертикальных осцилляций для каждой пары интегралов площадей и энергии. При вычислении значений интеграла энергии и энергии вертикальных осцилляций использовалась геометрическая прогрессия, значение знаменателя геометрической прогрессии в разных экспериментах варьировалось от 3 до 5. Как оказалось, такой подход позволяет охватить значительную часть фазового пространства модели, и в итоговом распределении встречались как орбиты диска, так и орбиты гало, как почти круговые, так и сильно вытянутые.

Использование М2М-алгоритма для вычисления весов орбит позволяет добавить кинематические параметры в ряд характеристик, которые необходимо воспроизвести в построенной модели. В качестве таких параметров нами были использованы тангенциальные компоненты скорости в ячейках, расположенных в плоскости г = 0. Величины, которые использовались как исходные, были получены из вычисленной круговой скорости в потенциале Флинна и др. [8] с соответствующей поправкой на дисперсию. Так как движения в центральной области изучены недостаточно хорошо и кривая вращения там уже не может быть построена на основании круговой скорости, достоверность значений постепенно снижается на отрезке от 3 кпк к центру. Эффекты, вызванные этим, хорошо заметны в полученных зависимостях средних значений скоростей и их дисперсий от радиуса (левые графики на рис. 3, 4). Тангенциальная скорость в центральной области значительно снижается, в тоже время значения дисперсий резко возрастают.

6. Потенциал модели Галактики

В данной работе для построения дискретной фазовой модели использовался трехкомпонентный потенциал Флинна и др. [8]:

Ф = Фн + Фс + Фд.

Потенциал имеет следующие составляющие:

1) гало

фя = \УнНг2+г20)]

Рис. 1. Распределение масс в ячейках и круговая скорость (аналитические значения для потенциала Флинна, три линии на левом графике соответствуют трем соседним слоям по Z (толщина слоя 0.25 кпк)), значения скоростей вычислены в слое |^| < 0.25 кпк; по осям отложены расстояние от центра(кпк), масса в ячейках (в относительных единицах) и скорость (км/с).

2) центральная часть (ядро и балдж)

фс = СМо, _ СМс- •

\]Г'2+ГЪ1 \!г2+г12

3) диск

ФД = ]Г ~СМд- —

п=1 у Д2 + [ап + а/-г2 + Ъ2]2

Характеристики модели и значения параметров, взятые из работы [8], приведены в таблице 1. Потенциал имеет протяженное гало. При построении дискретной системы N тел эта часть потенциала задавалась фиксированной. Таким образом, фактически моделировались ядро, балдж и диск, то есть характеристики, представленные на рис. 1.

Таблица 1. Параметры потенциала Флинна и др. [8]

Гало V# (км/с) Го (кпк) Центр, часть мсп(ме) гСп(кпк) Ь(кпк) Диск Мпп(Ме) ап (кпк) п

220 8.5 3.0 X 10й 2.7 0.3 6.6 х 101и 5.81 1

1.6 х Ю10 0.42 -2.9 х 1010 17.43 2

3.3 х 109 34.86 3

7. Результаты

Параметры построения модели, при которых были получены представленные результаты, даны в таблице 2. Параметр сглаживания и шаг интегрирования были выбраны на основе рекомендаций [9].

Таблица 2. Параметры модели

Размеры модели К(кпк) 2:(кпк) Число ячеек по К по Z Число дроблений ПО J ПО Е ПО Ег

30 5 200 20 250 6 6

Рис. 2. Распределение масс в ячейках в моменты Т = 0 и Т ~ 6 млрд. лет, три графика соответствуют трем соседним слоям по Z (толщина слоя 0.25 кпк); по осям отложены расстояние от центра (кпк) и масса в ячейках (относительные единицы).

Для исходной библиотеки орбит было вычислено 9000 орбит. Из них для создания модели М2М-алгоритмом были отобраны 2416 орбит с весом > 10-4 и 4830 орбит с весом 2х 10-5 < < 10-4. Переход к дискретной модели был произведен по алгоритму,

разработанному и опробованному в [5], суть алгоритма в пропорциональности количества звезд на орбите и веса этой орбиты. При таком подходе на каждой траектории с Шг > 10-4 располагалось более 5 «звезд», а на траекториях с 2 х 10-5 < < 10-4 от

1 до 5 «звезд». В результате была построена дискретная модель из 50 000 тел. Распределение плотности не получилось гладким (рис. 2) — это отражение конечного числа ячеек и набора орбит.

В ходе исследования эволюции модели при помощи численного моделирования задачи N тел было обнаружено, что модель за первые 6-10 оборотов (1.5—2.5 млрд. лет) приходит в новое состояние с параметрами, близкими к начальным. В дальнейшем характеристики стабилизировались и почти не менялись. Эволюция модели показана на рис. 2—4: на левых графиках представлено начальное состояние построенной модели (сразу после перехода к дискретной модели), на правых графиках — модель примерно через 6 млрд. лет. Параметры модели отошли от начальных: так центральная область модели заметно разогревается — на расстоянии до 3 кпк дисперсии возрастают относительно значении, вычисленных при Т = 0. Также возросла толщина диска (дисперсия по Z). В числе возможных причин «диффузия» частиц в фазовом пространстве, связанная с дискретностью его заполнения орбитами; также возможна неравновесность построенной модели вследствие ряда использованных аппроксимаций.

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

8. Обсуждение и выводы

Проблему с чрезмерно большой дисперсией скоростей в первоначальных вариантах удалось решить, и в настоящий момент представленный инструментарий может использоваться для исследования различных потенциалов, применяющихся при мо-

Рис. 3. Профили скоростей в моменты Т = 0 и Т ~ 6 млрд. лет; по осям отложены расстояние от центра (кпк) и скорости (км/с).

100

80

60

40

20

100 1 1 1 .

\ е п ее

>•, \ 2

80 ' п '

Л\ . 60 ■ \\ ■

\ V \ \

\\ч

■ — - ■ 40

■ 20 ■ ^

0

10

15

20

10

15

20

Рис. 4. Профили дисперсий скоростей в моменты Т = 0 и Т ~ 6 млрд. лет; по осям отложены расстояние от центра (кпк) и дисперсии (км/с).

0

0

5

0

5

делировании Галактики и других звездных систем. Изучение эволюции построенных фазовых моделей показывает, что модифицированный метод Шварцшильда позволяет строить по заданному потенциалу звездные системы с требуемым распределением плотности. Остается не до конца решенной проблема с эволюцией центральной области и утолщением диска.

Метод имеет следующие положительные качества: сравнительно быстрая адаптация к различным потенциалам, моделирующим звездные системы; отсутствие зависимости от априорных знаний об интегралах движения (необходимо лишь заполнить фазовое пространство модели большим количеством разнообразных траекторий); в результате работы метода получаем координаты и скорости частиц равной массы, которые могут быть использованы как входные данные при исследовании задачи N тел.

Одним из минусов является необходимость избыточных вычислений: для заполнения фазового пространства модели количество орбит должно быть достаточно велико, чтобы покрыть максимально возможный объем и обеспечить необходимое разрешение внутри области модели. Однако при этом для итогового распределения отбирается лишь часть просчитанных орбит. Другим минусом является зависимость характеристик построенной модели от пространственного разрешения (размера ячеек), например, происходит увеличение дисперсии при размере ячеек более 1 кпк. Этот недостаток можно

отнести к трудностям методов «грубой силы», однако, существует некоторое значение размера ячейки, вполне достижимое, уменьшение которого уже не влияет на полученный результат. Так, в случае моделирования Галактики достаточным будет разрешение ^50 пк в центральных частях и ^200 пк в окрестности Солнца. В выполненных вычислениях с потенциалом Флинна и др. [8] для Галактики удалось достичь разрешения в центральных частях ~10 пк, в окрестности Солнца ^200 пк, на периферии модели ^500 пк.

Литература

1. Schwarzschild M. A numerical model for a triaxial stellar system in dynamical equilibrium // Astrophys. J. 1979. Vol. 232. P. 236-247.

2. Кутузов С. А., Осипков Л. П. Двухкомпонентная модель гравитационного поля Галактики // Астрон. журн. 1989. Т. 66. Вып. 5. С. 965-973.

3. Binney J., Tremaine S. Galactic Dynamics. Princeton University Press, 1987. P. 223.

4. Башаков А. А., Питьев Н. П. Построение самосогласованных моделей звездных систем методом Шварцшильда // Вестн. С.-Петерб. ун-та. 2007. Сер. 1. Вып. 3. С. 151-159.

5. Башаков А. А. Тестирование фазовых моделей звездных систем, построенных методом Шварцшильда // Вестн. С.-Петерб. ун-та. 2008. Сер. 1. Вып. 4. С. 131-143.

6. Syer D., Tremaine S. Made-to-measure N-body system // MNRAS, 1996. Vol. 282. P. 223233.

7. De Lorenzi F., Debattista V. P., Gerhard O., Sambhus N. NMAGIC: Fast parallel Implementation of a x2-Made-To-Measure Algorithm for Modeling Observational Data // MNRAS, 2007. Vol. 376. P. 71-88.

8. Flynn C., Sommer-Larsen J., Christensen P. R. Kinematics of outer stellar halo // MNRAS, 1996. Vol. 281. P. 1027-1032.

9. Родионов С. А., Сотникова Н. Я. Оптимальный выбор параметра сглаживания и шага интегрирования в N-body экспериментах // Астрономический журнал. 2005. Т. 82. №6. С. 527534.

10. Rodionov S.A., Orlov V. V. Phase models of the Milky Way stellar disc // MNRAS, 2008. Vol. 385. P. 200-214.

11. Родионов С. А., Сотникова Н. Я. Итерационыый метод построения равновесных (N-bo-dy)-моделей звездных дисков // Астрономический журнал. 2006. Т. 83. №12. С. 1095-1114.

Статья поступила в редакцию 19 марта 2009 г.