Построение самосогласованных моделей звездных систем методом Шварцшильда Текст научной статьи по специальности «Математика»

АСТРОНОМИЯ

УДК 524.3/4-32

А. А. Башаков, Н. П. Питьев

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

1. Введение

Одной из важных задач динамики галактик является построение самосогласованных моделей звездных систем. Существует несколько основных принципов построения гравитационных моделей: возможно использование теоретических методов, численных подходов, а также их комбинаций, но универсального способа не разработано [1-6]. Один из полуэмпирических методов нахождения функции распределения фазовой звездной плотности был предложен Шварцшильдом [7]. Он основан на нахождении набора звездных орбит, распределение которых воспроизводит заданный гравитационный потенциал, в котором они вычислялись, и определения их весов. Использование этого метода не зависит от числа изолирующих интегралов движения, существующих для данной системы или заданного гравитационного потенциала, но предполагается, что система находится в стационарном или квазистационарном состоянии.

В данной работе реализован численный подход на основе метода Шварцшильда для построения самосогласованных моделей и испытан при получении кинематических характеристик звездной системы с потенциалом Кутузова—Осипкова для двухкомпонентной модели Галактики.

2. О программной реализации метода ^Шварцшильда

Для получения фазовых самосогласованных моделей был разработан на основе метода Шварцшильда программный комплекс, который для систем с конечной массой и осевой симметрией позволяет по заданному распределению плотности построить фазовую модель. Набор программ включает в себя:

• модуль, отвечающий за разбиение пространства на заданное число трехмерных ячеек и накопление в них необходимых параметров;

* Работа выполнена при финансовой поддержке РФФИ (грант №04-02-17408), ГФЕН КНР (грант №04-02-39026/19271121) и Совета по грантам Президента РФ для государственной поддержки молодых российских ученых и ведущих научных школ (грант № НШ-4929.2006.2).

© А. А. Башаков, Н. П. Питьев, 2007

• модуль генерации начальных данных для набора траекторий;

• модуль интегрирования траекторий для заполнения ячеек с заданной точностью;

• модуль, работающий с базой данных;

• модуль определения набора и весов орбит с помощью симплекс-метода;

• модуль вычисления кинематических параметров полученной модели.

Все программы были написаны на языке С/С+—+. Для получения траекторий звезд в заданном потенциале использовался метод численного интегрирования уравнений движения, предложенный Эверхартом [8]. Вычисления проводились с двойной точностью.

Кратко отметим основные пункты метода Швацшильда. Конфигурационное пространство, занимаемое звездной системой, делится на заданное количество замкнутых ячеек. Затем по определенным правилам выбираются начальные данные траекторий и строится библиотека орбит. Каждая орбита прослеживается на несколько десятков или сотен оборотов. Во время численного интегрирования уравнений движения для каждой ячейки производится накопление ряда параметров: время нахождения в ячейке, средняя скорость и т. д. Предполагается, что суммарное время, в течение которого звезда находится в ячейке, пропорционально массе, создаваемой данной орбитой в рассматриваемой ячейке,—орбита представляется как «множество звезд» равной массы, равномерно распределенных вдоль нее. При нахождении из библиотеки орбит оптимальной выборки для воспроизведения наилучшим образом плотности, соответствующей заданному потенциалу, необходимо одновременно найти «количество звезд» находящихся на каждой из выбранных траекторий или, иначе, веса орбит. Для этой цели используется некоторый вариант симплекс-метода [9, 10]. При поиске решения накладываются определенные ограничения. В частности, предполагается, что потенциал удовлетворяет условиям неотрицательной вещественной плотности и конечной полной массы. В данной работе рассмотрен вариант стационарного осесимметричного потенциала и (К, г), в котором траектории имеют два классических интеграла движения — площадей J и энергии Е.

3. Библиотека орбит

Библиотека орбит должна содержать достаточно широкий спектр траекторий, расположенных во всем предполагаемом объеме моделируемой звездной системы, и включать орбиты как близкие к круговым, так и имеющие достаточно большие амплитуды движения по К и г (соответственно, расстояния от оси и от плоскости симметрии системы). С учетом этих требований, выбор начальных данных траекторий при создании библиотеки орбит производился следующим образом. Исходя из размеров модели, определялся интервал значений интеграла площадей J (от = 0 до

увеличивающимся шагом от Ітіп к 1"тах) выбиралось заданное количество NJ значений интеграла площадей (її, І2,..., JNJ), а по ним находился набор соответствующих значений радиусов круговых орбит К і, 1 < і < (рис. 1). Для каждого значения іі определялся интервал возможных значений интеграла энергии Е, где Етіп соответствовало круговому движению при данном значении интеграла площадей и Етах — энергии, при которой контур области возможных движений достигает границ модели. Из найденного интервала (Етіп, Етах) выбирался дискретный набор NE значений интеграла энергий (Еї, Е2,..., Е^в) (рис. 1). Для каждого значения интеграла энергии Е^ вычислялся промежуток [Ктіп, Ктах], соответствующий пересечению левого и пра-

тах

этом интервале по определенному правилу (с постепенно

вого края контура нулевой скорости (Ej = .]‘2/2К'2 — и (К, г)) с плоскостью г = 0, и на найденном интервале К определялся заданный набор Жд значений Кк начальных («стартовых») точек траекторий для библиотеки орбит, соответствующих выбранной паре значений 3^ и Ej. Координаты начальных точек траекторий лежали в плоско-

Такой вариант выбора начальных данных не охватывает всех траекторий, возможных при данной паре значений 3^ и Ej, в частности, не включает траектории, которые всегда пересекают плоскость симметрии системы наклонно. Рассмотрение различных вариантов генерации начальных данных и влияния способа разбиения по значениям 3^ и Ej предполагается в дальнейшем.

В итоге получается набор из К орбит с заданными наборами 3, E, К. Время счета каждой траектории определялось по стабилизации значения доли времени пребывания траектории в каждой из ячеек, пересекаемой данной траекторией, ко всему времени. Если изменение относительного времени во всех ячейках, через которые проходила траектория, становилось ниже заданного уровня (этот параметр задавался), то предполагалась достаточная стабилизация среднего времени пребывания, и происходил переход к вычислению следующей орбиты. Для каждой траектории в процессе интегрирования уравнений движения накапливались необходимые характеристики орбиты, и при завершении интегрирования траектории накопленные параметры размещались в базе данных.

4. Аппроксимация заданной плотности симплекс-методом

Для нахождения выборки и весов орбит, с помощью которых плотность, соответствующая заданному потенциалу, моделируется наилучшим образом, использовался симплекс-метод, в котором решение системы из N линейных уравнений с К неизвестными ищется с учетом ряда ограничивющих условий [9, 10]. В нашем случае N равнялось числу пространственных ячеек, К — числу траекторий в библиотеке орбит и налагались следующие условия на веса т траекторий:

сти симметрии 2 = 0, вектор скорости перпендикулярен к плоскости (К = 0, г > 0),

начальные прямоугольные координаты и скорости имели вид

(1)

Е

Диапазон

изменения

1

^"Круговая орбита

Рис. 1. Контуры областей возможного движения и распределение значений Jг,Ej.

К

т —> тах, ЛШ < Р, т > 0, г Е {1,К},

(2)

/

где Л =

ац . .. а1к

массив, содержащий суммарное время aji пребывания

\ ад 1 ... адк

I т1

— столбец с

\ тк

траектории г в соответствующей пространственной ячейке ], Ш

(

искомыми весами орбит, Р =

Р1

величины, пропорциональные пространствен-

\ РД

ной плотности из уравнения Пуассона Ди = —4пОр. Точность аппроксимации распределения плотности зависела от задаваемого параметра, контролирующего в процессе счета каждой орбиты стабилизацию относительного времени пребывания траектории во всех из пересекаемых ею ячейках к полному времени счета.

5. Определение кинематических характеристик модели

В процессе интегрирования каждой траектории при пересечении пространственных ячеек накапливались следующие данные: время нахождения орбиты в ячейке, скорость и квадрат скорости пересечения ячейки с учетом времени пребывания внутри нее в прямоугольных (Х,у, г) или цилиндрических (К, в, г) координатах.

По окончании интегрирования г-й орбиты для каждой 3-й ячейки, которую она пересекла, вычислялись суммарное время Т3 нахождения орбиты г в ячейке ], средняя скорость Уц пересечения ячейки ПО компонентам, дисперсия скорости ИУц = (Ту. . в ячейке по компонентам. Эта информация, структурно связанная с г-й орбитой и 3-й ячейкой модели, заносилась в базу данных.

После завершения интегрирования траекторий и создания библиотеки орбит определялась выборка, моделирующая задаваемое потенциалом распределение плотности (п. 4). Затем производилось сравнение найденного распределения с заданным, и если общая масса модели отличалась от исходной в рассматриваемом объеме меньше 5-10%, то проводился анализ кинематических характеристик полученной модели. Для каждой 3-й ячейки определялся поднабор орбит 1(о), пересекавших рассматриваемую ячейку.

По этому поднабору 1(3) с учетом весов траекторий в данной 3-й ячейке вычислялись

по каждой компоненте х, у, г или К, в, г дисперсии оу6 и средние скорости Vj:

У?' = ^3^3’ (3)

3 *е*(3)

Ъ = ^7 Л аЬ^^з> (4)

3 *ё1(3)

а

^ ^ , • Тз , (6)

*ё£(3)

где т —вес орбиты, Тз —суммарное время нахождения г-й орбиты в 3-й ячейке.

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

6. Потенциал Кутузова—Осипкова

В качестве пробного потенциала была взята двухкомпонентная модель Галактики Кутузова—Осипкова [11, 12]. В ней потенциал имеет вид

и (Я, г) = Цо[а^а(£) + Ъ^ь(т)],

а + Ъ = 1,

где ио = и(0,0) —потенциал в центре, а и Ъ — веса компонентов, и <^>ь их безразмер-

ные потенциалы. Компонента а, сплюснутая и более протяженная, моделирует распределение в диске и гало:

-1 -1

<Ра(0 = а а - 1 + Vх! +

компонента Ъ с быстрым падением плотности представляет ядро и балдж:

^ь(т) = (1 + Ат)

-1

Здесь

+ 2(1 -е) |^£2 + С2 - < П = Я/то,

Vх??2 + С2

С = г/то,

Я и г —цилиндрические координаты, то — единица длины.

Модель содержит два масштабных (то, ио) и пять структурных (а, к, е, А, а) параметров, определяющих ход круговой скорости, сплюснутость и веса компонент. Значения принятых параметров, приведенных в таблице 1, оценка которых была выполнена по наблюдательным данным Галактики, взяты из работы [12].

Таблица 1. Параметры потенциала

2

2

С

т

т

С/о, 10ь км2/с2 3,81

Го, КПК 2,54

а 0,37

а 4,17

лг 0,69

£ 0,11

А 6,60

Г рафик круговой скорости в плоскости симметрии модели при данных параметрах потенциала Кутузова—Осипкова представлен на рис. 2.

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

Был получен ряд моделей, отличающихся размерами по Я и г, наборами траекторий, величиной и количеством ячеек, точностью аппроксимации исходного распределения плотности. На рис. 3-6 в качестве примера приводится один из найденных вариантов.

крс

Рис. 2. Кривая круговой скорости в двухкомпонентной модели Осипкова—Кутузова.

Пространство было разбито на ячейки со стороной 1.4 кпк, протяженность модели по Д равно 19.6 кпк, по г — ±18.2 кпк. Была учтена симметрия модели относительно оси г и плоскости г = 0. В библиотеку орбит было включено 1728 траекторий со спектром значений Ji(NJ = 12), Е^ (Же = 12) и Д^(Жд = 12), описанным в п. 3. Каждая орбита в зависимости от достижения стабилизации относительного времени пребывания во всех ячейках вдоль траектории считалась на интервале от нескольких десятков оборотов (орбиты с малой амплитудой по Д или г) до 100-150 оборотов (траектории с большой амплитудой по Д или г). Точность интегрирования контролировалась по сохранению интегралов площадей и энергии, относительное изменение J и Е было меньше 10-8—10-9. С условием аппроксимации распределения плотности не хуже 5% симплекс-методом было найдено решение, включающее 264 орбит (15%) с ненулевыми весами. Этот набор с учетом полученных весов орбит достаточно хорошо приближает распределение массы в диске, балдже и гало. На рис. 3-4 в логарифмической шкале показаны профили плотности по Д и г, соответствующие исходному потенциалу и смоделированные по найденной выборке орбит. Значения плотности нормированы к максимальному значению в районе центра модели. В большей части рассматриваемых интервалов по Д и г кривые практически совпадают, различие начинается на периферии модели, где плотность существенно уменьшается и отличается от значений в центре примерно на два порядка.

Кривая вращения построенной модели (рис. 5) близка к плоской на уровне 180— 190 км/с, примерно от радиуса 4 кпк до периферии и проходит существенно ниже кривой круговой скорости (рис. 2). Это связано с большими дисперсиями скоростей по Д, в и г (рис. 6) полученной модели — модель заметно «горячее», чем показывают, например, данные наблюдений в нашей Галактике. Максимальные значения дисперсии по

density density

R axis (kpc)

Рис. 3. Сравнение профилей заданной и смоделированной плотностей в плоскости = 0.

Z axis (kpc)

Рис. 4. Сравнение профилей заданной и смоделированной плотностей вдоль оси z.

кт/в кт/в

крс

Рис. 5. Кривая вращения.

крс

Рис. 6. Ход дисперсий скоростей по Я, в и г.

всем компонентам оказались в районе центральной области модели, с удалением от центра значения дисперсий уменьшаются и сближаются. В данной модели до расстояния примерно Д=12 кпк наибольшее значение дисперсии скоростей получилось по компоненте z (рис. 6). Вероятно, это связано с тем, что в данном варианте для библиотеки орбит были рассмотрены только траектории, стартующие перпендикулярно плоскости симметрии модели, и отсутствуют траектории с наклонными начальными скоростями. Характер кинематических свойств построенной модели существенно зависит от использованного набора траекторий.

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

Summary

A. A. Bashakov, N. P. Pityev. Constructing self-consistent models of stellar systems by Schwarzschild’s method.

The set of programs for constructing self-gravitating models of stellar systems in a given distribution function of density using Schwarzschild’s method is developed. As a trial gravitational potential the two-component model by Kutuzov and Ossipkov has been considered for which the set of phase models and kinematic parameters have been found.

Литература

1. Binney J. J., Tremaine S. Galactic dynamics. Princeton, 1987. 733 p.

2. Hernquist L. N-body realization of compound galaxies // Astron. and Astrophys. Suppl. Ser. 1993. Vol. 86. N 2. P. 389-400.

3. Kuijken K., Dubinski J. Nearly self-consistent disc-buldge-halo models for galaxies // Mon. Not. Roy. Astron. Soc. 1995. Vol. 277. P. 1341-1353.

4. Binney J. J., Ossipkov L.P. // Astrophysics at Frontier of Centuries / Edited by N. S. Kardashev, R. D. Dagkesamansky, Yu. A. Kovalev (Moscow: Janus-K, 2001). P. 181-183.

5. Widrow L. M., Dubinski J. Equilibrium disc-buldge-halo models for the Milky Way and Andromeda galaxies // Astrophys. J. 2005. Vol. 631. P. 838.

6. Родионов С. А., Сотникова Н. Я. Итерационный метод построения (N-body)-моделей звездных дисков // Астрон. журн., 2006. Т. 83. Вып. 12. С. 1-19.

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

8. Everhart E. Implicit single-sequence methods for integrating orbits // Celest. Mech., 1974. Vol. 10. P. 35-55.

9. Ашманов С. А. Линейное программирование. М., 1981.

10. Васильев Ф.П., Иваницкий А. Ю. Линейное программирование. М: Факториал Пресс, 2003.

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

12. Кутузов С. А., Осипков Л. П. Оценка параметров двухкомпонентной модели Галактики интервальным методом // Вопросы небесной механики и звездной динамики. Алма-Ата, 1990. С. 110-116.

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