CC BY
35
УДК 524.3/.4 Вестник СГ16ГУ. Сер. 10, 2008, вып. 1
С. А. Кутузов, Н. В. Распопова
РЕЛЬЕФ ПОЛЯ СИЛ И ОРБИТЫ В МОДЕЛИ ГАЛАКТИКИ *>
Введение. Гравитационное поле в моделях галактик обычно описывается с помощью потенциала как функции координат (а также времени в нестационарных моделях). Выражения для этих функций весьма разнообразны. Однако их геометрические образы (графики зависимости от координат, эквипотенциали и т. п.) на вид могут отличаться незначительно, а то и быть неразличимыми. Между тем соответствующие поля сил различаются гораздо заметнее, а распределения масс даже разительно. Это связано с применением операций дифференцирования к потенциалу: напряженность поля есть градиент потенциала, а плотность масс пропорциональна дивергенции напряженности.
С нолем сил непосредственно связаны свойства галактических орбит, так как напряженность определяет правую часть уравнений движения. В общем случае зависимость силы от координат описывается гиперповерхностью в четырехмерном пространстве. В случае систем с ротационной или сферической симметрией размерность пространства уменьшается на единицу или две. В поверхности силы можно выделять такие элементы рельефа как гребни и впадины. Сила, препятствующая потоку звезд от центра наружу, на гребнях максимальна. Контуры изодинам гребней определяют динамическое ядро системы. В ядре в общем случае сила почти линейна по координатам, что обусловливает независимые колебания звезд, т. е. ящикообразные орбиты, или фигуры Лиссажу [1]. Динамическую толщину сплюснутых систем определяют гребни компонента силы, перпендикулярного плоскости системы.
В данной работе исследуется связь рельефа поля сил со свойствами орбит в семействе моделей гравитационного поля галактик с ротационной симметрией.
1. Функции описания и единицы. Все фигурирующие в настоящей работе переменные являются безразмерными. Переход к размерным величинам выполняется умножением на соответствующие размерные единицы. Основными служат единицы длины, потенциала и массы:
г, Ф, М = г <1/(3, (1)
где (7 - гравитационная постоянная. Кроме того, для наглядности и удобства сравнения конкретных моделей будем использовать еще и безразмерные единицы, характерные для поля тяготения и распределения масс. Для систем с ротационной и зеркальной симметрией применим цилиндрическую систему координат Д, 9, г.
Одну из возможностей определения безразмерной единицы длины дает круговая скорость У(В), получаемая из потенциала г):
У2{К) = -Л
д-ф(Я, г)
Ж
(2)
г=0
В рассматриваемых моделях это унимодальная функция. Положение ее максимума гс можно принять в качестве «динамической» единицы длины.
"' Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант №04-02-17447), а также программы поддержки ведущих научных школ России (грант X» НШ 4929.2006.2).
© С. А. Кутузов, II. В. Распопова, 2008
Значение максимальной круговой скорости может служить «динамической» единицей скорости
VC = V(rc) . (3)
Еще одной безразмерной «динамической» единицей скорости может быть скорость освобождения в определенной точке системы ve, например в месте максимума круговой скорости или в центре системы:
vesc(R, z) = у/2z), ve = vesc(0, 0). (4)
В дальнейшем будем (как это было предложено в [2]) задавать потенциал аналитическим выражением со свободными параметрами и с одним аргументом Зависимость последнего от координат определяют эквипотенциали:
Ч> = <Ж), = /(Д, г). (5)
Случай произвольного потенциала z) можно свести к предыдущему, положив
í2 = ш z), <р = е-
2. Рельеф гравитационного поля. В цилиндрической системе отсчета сила F имеет два отличных от нуля компонента Fr и F,, которые являются компонентами градиента потенциала (5):
г dip ,df dip ,д f
Fr=M = ipM> = (6)
где
dip
<S(0 = •
Уровенные поверхности модуля силы
г) = , (7)
являющегося функцией йиг, будем называть изодинамами. В трехмерном евклидовом пространстве изодинамы в анализируемых моделях суть поверхности вращения. Поэтому ограничимся рассмотрением изолиний на меридиональной плоскости, получаемых как сечения соответствующих поверхностей, и будем называть их также изодинамами. Они аналогичны изодинамам напряженности магнитного поля Земли, а также изогипсам, изображаемым на географических картах. Нас интересуют экстремумы модуля силы Р (7). Максимумы дают гребни, а минимумы - впадины рельефа. Заметим, что значение силы вдоль гребня может изменяться, т. е. гребень в общем случае не является изодинамой. Тогда на гребне можно выделять верхние и нижние точки. Наглядным примером служит однородный шар радиуса Ь с потенциалом
(1-г2/(ЗЬ2), г О; <р(г) = ^ г2 = Д2 + г2. (8)
(2Ь/(Зг), Ь <С г ^ оо .
Ему соответствует известное выражение силы
2r/(362), OíC r^b;
F(r) =
2Ь/(Зг2), b^r^oo.
Сила имеет максимальное значение на границе, т. е. положение гребня совпадает с граничной сферой.
Для вычисления экстремумов F можно применять различные способы. Возможны следующие алгоритмы:
A. Задаем набор значений R, и для каждого значения находим такое z, при котором максимально
B. Задаем набор значений z, и для каждого значения находим такое R, при котором максимально
C. Задаем набор значений системоцентрических широт /?, подставляем в (7) выражения R = г cos (j и z = г sin /3, а затем для каждого значения находим такое г, при котором F2 максимально на заданном луче. Контур экстремума получается в полярных координатах.
Здесь сила предполагается дифференцируемой функцией. Простейшим численным методом нахождения экстремума, не требующим вычисления производных, является деление отрезка пополам [3]. Считаем, что гребни, пики и впадины включают центр системы. Если же рельеф поля сил сложный и внутри квадранта меридиональной плоскости имеются изолированные гребни, пики и впадины, то описанные способы требуют доработки.
3. Семейство моделей. Метод раздельного задания эквипотенциалей и закона потенциала был предложен в [2]. В работе [4] перечислены десять моделей, которые представляют собой частные случаи рассмотренного в [2] семейства моделей. Это модели П. П. Паренаго, Г. М. Идлиса, Г. Г. Кузмина (две), А. Тоомре, X. Пламмера, М. Миямото и Р. Нагаи, Д. Линден-Белла, М. Энона, Г. Г. Кузмина и Ю.-И. К. Велт-манна. Построенные там эквипотенциали были обобщены позднее [5], но впоследствии почти не использовались. Согласно (5), эквипотенциали определяются заданием функции /(Л, z). Были рассмотрены два трехпараметрических семейства эквипотенциалей с параметрами 7, е, fi:
1(7, е, /х) : f(R, z) = r2 + 2fi(q - ер + e2);
[11(7, e, fx) : f(R, z) = 7'2 + 2/z(l - e)(q - es),
где __
p = ^/s2 + 2q + z2 +e2 , q = \Je2s2 + г2 , s = \/l + 7R2 .
Здесь мы ограничимся частным случаем 7 = 0, когда оба семейства совпадают. Эквипотенциали задаются функцией
?=г2 + 2ц(1 - е)(Я - е), ц = \/е2 + г2 . (9)
Заметим, что в модели Сато [6] размерный потенциал можно привести к следующему виду:
Ф(0= СЯ= Г» +2(1 -£)(„-£).
Г у/ 1 - +
Очевидно, что эти эквипотенциали - частный случай анализируемых нами моделей при ц = 1. В сферической системе с е = 1 потенциал имеет разрыв в центре, где £ = 0. Система вырождается в материальную точку с кеплеровым потенциалом. Примем иной закон потенциала:
<Ж) = --, Р = а - 1, (Ю)
р + у/1 + <2
который в плоскости г = 0 совпадает с потенциалом Кузмина-Маласидзе [7]. Таким образом, принятое семейство моделей имеет четыре безразмерных параметра е, /х, а, /с, которые являются структурными. Оно дает широкие возможности моделирования различных галактик и звездных скоплений. Потенциал нормирован так, что в центре Р(0) = 1.
Асимптотика потенциала с учетом принятых единиц (1) приводит к выражению для безразмерной массы модели
АЛ а у/К
Распределение массы в модели описывает пространственная плотность. При заданном потенциале она легко получается из уравнения Пуассона, которое приобретает вид [2,
5]
¿V а?
Дифференцирование потенциала дает функции <£>'(£)> <¿>"(0 0) ():
4тге(Д, г) = -Р{Я, 2)^(0 - *)¥>"(0 , ¥>"(0 =
V'«) =
/ м У3
\2а/ а — Р<р ' . За - 2/3<р
п,( к ^ 5 2Р</>
Они зависят лишь от закона потенциала. Функции же Р(Д, г), (?(/?, г) определяются заданием эквиденсит:
Р(Д, г) - У2у1/(Д, г), <?(Д, г) = + '
Оператор набла берется в цилиндрической системе координат. Беря частные производные от выражения (9), получаем [5]
Р{Я, г) + , <2(11, г) = 4Я2 + 4г2 и;2 , ю = 1 + и}—^-.
<Г 7
Полная сила Р вычисляется с помощью компонентов (6):
Р\ = Р2=4^'2(0г2и,2. (11)
Отсюда видно, что в случае е = 1 поле сил сферически симметрично, так как, по (9), £ = г. В случае е = 0 сила испытывает скачок в плоскости г = 0, поскольку q = \г\. В плоскости находится бесконечно тонкий диск, погруженный в непрерывное гало; безразмерная масса диска равна /л [5].
4. Рельеф поля сил в конкретных моделях. Были вычислены и построены графики рельефа силы (11) на меридиональной плоскости для нескольких наборов значений параметров семейства моделей (таблица). В таблице приведены единицы длины гс и скорости ус, определенные через круговую скорость, как это предлагалось в (2), (3). Параметр к означает отношение массы ядра, заключенной внутри соответствующего контура гребня сил, к полной массе модели. Для рассматриваемого набора параметров оно составляет всего 9%.
Некоторые наборы значений параметров модели
№ Рисунок £ а К Гс Vc к
1 1, а; 2, а 1 0.5 2 1 2.20 0.59 0.089
2 1, б; 2, б 0.1 0 2 1 2.20 0.59 -
На рис. 1 передний край поверхности изображает зависимость силы от Лв плоскости г = 0. Центральная часть соответствует впадине с Р = 0. Она окружена замкнутым гребнем переменной высоты (нижайшие точки находятся в плоскости 2 = 0). Получается как бы кратер.
При е = 0 (рис. 1, б\ 2, б) кратер исчезает, и на его месте появляется глобальный пик. Рост ц повышает края кратера и делает его внешний склон более крутым. Рост а раздвигает границы кратера, а внешний склон делает более пологим.
Рис. 1. Рельеф сил в меридиональной плоскости. Объяснение в тексте.
Рельеф описывает карта изодинам (рис. 2). Для значений силы, близких к величине на вершине гребня FmatX, изодинамы не разрываются. Далее, с уменьшением силы каждая изодинама разрывается на две дуги, начиная со значения 0.73Fmax^ Внешняя и внутренняя дуги соответствуют 0.43Fmax.
5. О траекториях и орбитах. Уравнения движения. Целесообразно различать понятия «траектория» и «орбита». Траектория - это непрерывная линия, которую фактически описывает движущееся тело (звезда, материальная точка) по отношению к заданной системе отсчета. За конечный интервал времени траектория приобретает конечную длину. Орбита - это потенциально возможная траектория на (полу)бесконечном интервале времени, соответствующая данному полю сил и начальным (краевым) условиям. Траекторию можно наблюдать, а орбиту рассчитывать a priori. Например, для кеплерова движения барицентрические орбиты являются
О 0.5
Рис. 2. Изодинамы. Объяснение в тексте.
коническими сечениями, а траектории - их дугами. Возможно многократное покрытие орбиты траекторией.
Для построения траектории материальной точки (изображающей звезду) будем решать задачу Коши, в которой система обыкновенных дифференциальных уравнений второго порядка имеет вид
г = V<Ж),
где г - радиус-вектор звезды; V - векторный оператор набла. В декартовых координатах имеем
г =
V =
д а д
^дх' ду' дг
Значок Т предписывает транспонирование вектора. Начальные условия заданы для произвольного момента времени Ьо-
г(е£0)=г0, у(ег0)=^0,
где v = г - вектор скорости тела.
Время Ь (как и все переменные и функции) безразмерно. Уместно использовать в качестве безразмерной единицы времени период движения по круговой орбите с наибольшей скоростью ус (3). Размерная же единица времени £ следует из основных единиц (1):
2жгс - г
t — -- 1 — -
ис — , и --.
Ус у/ф
В случае принятых эквипотенциалей (9) и закона потенциала (10) находим
Уу>(0 =2<Р'(0(Х, У > .
Вводя компоненты скорости ух, уу, можно перейти к очевидной системе шести обыкновенных дифференциальных уравнений первого порядка, которые не выписываем.
Так как модель стационарна и обладает ротационной симметрией, то имеют место интегралы энергии Е и площадей I [8]
(12) 37
Здесь vr, vg, vz - компоненты скорости в цилиндрической системе координат, vm представляет собой меридиональную скорость, Е - энергия с обратным знаком. Используя интеграл площадей, уравнения движения можно записать как систему пяти уравнений первого порядка [9]
R = vr, vr = I2/R3 + Fr , ' d = I/R2, vz=Fz, . i = vz.
6. Исследование траекторий. Выявим связь траекторий с рельефом поля. Соответственно для задачи Коши будем выбирать начальные условия. Выбором начальных точки и скорости выделяются следующие группы траекторий.
Первая группа:
I - для звезд, падающих с контура гребня силы, когда полная начальная скорость v0 = 0. На диаграмме Линдблада [10] с осями I, Е начальные точки лежат на оси Е с координатой Е = Точки, изображающие орбиты, очевидно, сохраняют свое
положение во время движения.
Вторая группа:
Начальные точки берутся на изодинаме F — к • Лтах, где Fmax значение силы на верхней точке гребня, к = 0.1,0.2,0.3 в зависимости от вида рельефа сил для выбранного набора параметров:
Па - полная начальная скорость vo = 0;
lib - падающие с контура нулевой меридиональной скорости, когда начальная меридиональная скорость vm = 0.
На диаграмме Линдблада точки лежат на характеристической параболе с уравнением
= (13)
которое получается из (12) исключением vg и приравниванием нулю меридиональной скорости. Парабола имеет вершину на оси Е в точке Е = </?(£)• При этом vg выбирается пропорциональной локальной скорости освобождения (4): vg = 0.1г • vesc, г = 1,..., 10. Согласно (12), интеграл площадей имеет ту же сетку узлов, растянутую с коэффициентом R.
Вычисления проводятся в двух средах - Maple 9 и Borland Delphi 7. Система дифференциальных уравнений решается методом Мерсона либо методом Дормана-Принса [11], которые обеспечивают достаточно малую относительную погрешность интеграла энергии SE = АЕ/Ео ■ Здесь АЕ = Е-Е0, где Е0 ~ значение интеграла энергии в начальной точке, Е - в текущей точке, вычисляемое по формулам (12).
Прежде всего рассчитываются траектории звезд, падающих с различных точек гребня. Естественно сравнивать их с такими же траекториями в однородном шаре. Там, согласно выражениям для энергии (12) и потенциала (8), все падающие звезды колеблются по диаметрам. Скорость в центре равна круговой скорости на границе и составляет 1 /\/3 « 0.5774 в единицах (4). Если звезда начинает движение из центра с такой скоростью, то она достигает границы, где на мгновение замирает, а затем падает.
Ограничимся демонстрацией траекторий для одного набора параметров модели: ц = 1, е = 0.5, а = 2, к = 1. Рельеф сил в этом случае представляет собой кра-тероподобный гребень (см. рис. 1, а; 2, б). Две верхние точки гребня находятся на оси г, его нижняя точка лежит на окружности в плоскости z = 0. Соответствующие
координаты (R, z) в единицах (3) равны (0, ±0.27), (0.42,0). Отношение величины сил составляет 1.21.
Траектории свободного падения с контура гребня силы (группа I) изображены на рис. 3. Как и в однородном шаре, они лежат в меридиональной плоскости внутри эквипотенциали, проходящей через начальную точку. На рис. 4 также изображены изодинамы и контур гребня. В общем случае траектории не проходят через центр. Исключение составляют падения с верхних и нижних точек гребня, когда получаются колебания по отрезку оси z, либо по диаметру окружности в плоскости z = 0. В случае колебаний вдоль оси г погрешность 6Е не превышает 0.8 • Ю-8 для 150 колебаний (рис. 4, а).
На рис. 4, б в пространстве скоростей построена линия концов разнонаправленных векторов скоростей звезд, стартующих из центра и достигающих гребня либо с нулевой скоростью, либо по касательной. Для сравнения изображена линия скорости освобождения в центре (4).
0.2
.....' 0
50 100 150 0.2 0.6 1 1.4
Рис. 4• Относительная погрешность интеграла энергии (а) и линии концов векторов скоростей звезд (1) и скорости освобождения (2) (б).
4е-0.9
5Е
8е-0.9 1
В то время как траектории падения являются плоскими, траектории с отличным
от нуля интегралом площадей (группа IIb ) заполняют своими витками пространственную область внутри осесимметричного тора. Меридиональное сечение тора описывает уравнение (13). Это уравнение контура нулевой скорости и в то же время уравнение характеристической параболы на диаграмме Линдблада. В первом случае здесь фиксируются интегралы движения и изменяются со временем координаты, а во втором -наоборот: фиксируются координаты, а интегралы движения служат переменными.
В данной группе траекторий начальная точка лежит на контуре нулевой скорости. Как правило, орбиты получаются ящичными (рис. 5).
z а z б
Рис. 5. Траектории, в которых начальная точка лежит на контуре нулевой скорости.
Среди траекторий найдены периодические (рис. 6). При задании в начальных условиях малых возмущений координат в фазовом пространстве эти орбиты сохраняют форму, заполняя переменную по величине полосу, которая включает первоначальную орбиту.
Рис. 6. Периодические орбиты.
а - полная начальная скорость равна нулю (группа IIa); б - начальная меридиональная скорость равна нулю (группа IIb).
В случае е = 0, к ц2 — 1, когда в плоскости z = 0 присутствует диск, погружен-
ный в гало, поле становится бицентральным, и движение происходит вдоль множества плоскостей, наклоненных к плоскости диска [12].
Звезда перескакивает случайным образом с одной плоскости на другую. Как и у Хантера [13], получаются стохастические орбиты. Следуя работе [12], их можно описать гистограммой модуля момента импульса L = г х v , который становится случайной величиной. При этом интеграл площадей, являющийся z-компонентом момента импульса, по-прежнему сохраняется.
Заключение. Рассмотренное семейство моделей дает широкие возможности описания различных галактик и звездных скоплений и их динамики. Но намного более гибкий аппарат находим, пользуясь суперпозицией отдельных моделей этого семейства. Обозначим отдельный набор параметров
Рп = {£«, Рп, ап, кп} , n = l,N.
Тогда А^-компон'ентную модель гравитационного поля получим как взвешенную сумму потенциалов (9), (10):
N N
v - wnvn{z,n), Wn ~1'
n=l п—1
где wn - весовые нормированные коэффициенты. Точно так же складываются и плотности масс. Например, если взять N = 2, е\ = 0, £г > 0, а все остальные модельные параметры в обоих компонентах одинаковыми, то будем иметь пять свободных параметров, включая w\. В такой модели рельеф поля сил может включать кратер либо нет, либо кратер с центральным пиком.
Предполагается продолжить вычисление и анализ орбит в различных моделях, описывающих реальные галактики.
Авторы благодарят JI. П. Осипкова за ценные замечания. Summary
Kutuzov S. A., Raspopova N. V. Relief of force field and orbits in galaxy model.
A relief of a gravitational force field is described for a four-parametric family of galaxy models. The number of star trajectories with various initial conditions are computed. There are some closed stable orbits among them. Stochastic orbits occur in the particular case of a disk and halo model.
Литература
1. Binney ,/., Tremaine S. Galactic Dynamics. Princeton: Princeton Univ. Press, 1987. Vol. XV. 733 p.
2. Кутузов С. А., Осипков JI. П. Моделирование пространственного гравитационного потенциала звездных систем // Астрон. журн. 1980. Т. 57. С. 28-37.
3. Васильев Ф. П. Численные методы решения экстремальных задач. М.: Наука, 1980. 520 с.
4. Кутузов С. А., Осипков Л. П. Обобщенная модель крупномасштабного гравитационного поля галактик // Бюл. Абастуманск. астрофиз. обсерв. 1980. JV® 52. С. 93-108.
5. Кутузов С. А. Распределение масс в модели галактики, состоящей из диска и гало // Вестн. Ленингр. ун-та. Сер. 1: Математика, механика, астрономия. 1989. № 8. С. 79-85.
6. Satoh С. Dynamical models of axisymmetric galaxies and their application to the elliptical galaxy NGC 4697 // Pubis Astron. Soc. Japan. 1980. Vol. 32, N 1. P. 41-62.
7. Кузмин Г. Г., Маласидзе Г. А. Об одной форме гравитационного потенциала, допускающей решение задачи о плоских орбитах звезд в эллиптических интегралах // Публ. Таргуск. астрофиз. обсерв. 1969. Т. 38. С. 181-250.
8. Огородников К. Ф. Динамика звездных систем. М.: Физматгиз, 1958. 644 с.
9. Кутузов С. А., Осипков JI. П. Методы расчета галактических орбит звездных скоплений // Движение искусственных и естественных небесных тел. Свердловск: Изд-во Урал, ун-та, 1981. С. 46-62.
10. Линдблад Б. Динамика Галактики // Строение звездных систем / Пер. с нем.; Под ред. П. Н. Холопова. М.: Гос. изд-во иностр. лит-ры, 1962. С. 39-132.
11. Хайрер Э., Нерсеттп С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи / Пер. с англ. И. А. Кульчицкого, С. С. Филиппова; Под ред. С. С. Филиппова. М.: Мир, 1990. 512 с.
12. Kutuzov S. A. Orbits in a "disk+halo" galaxy model // ASP Conference Ser. 2004. Vol. 316. P. 37-42.
13. Hunter C. Disk-crossing orbits // Galaxies and Chaos / Eds: G. Contopoulos, N. Voglis. Lecture Notes in Physics. 2003. Vol. 626. P. 137-153'.
Статья рекомендована к печати проф. Л. А. Петросяном. Статья принята к печати 11 ноября 2007 г.
