ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ, 2020, том 46, № 5, с. 366-380

ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ ПЛАНЕТНОЙ СИСТЕМЫ

НА КОСМОГОНИЧЕСКИХ ВРЕМЕНАХ

¹Санкт-Петербургский государственный университет, Санкт-Петербург, Россия

Поступила в редакцию 13.02.2020 г.

После доработки 19.04.2020 г.; принята к публикации 28.04.2020 г.

Рассматривается динамическая эволюция планетных систем, структура которых близка к круговой

и компланарной. Исследование выполняется с помощью метода осреднения Хори-Депри в рамках

теории первого порядка по планетным массам. Для получения уравнений движения используются

удобный набор канонических элементов и редко употребляемая разновидность астроцентрических

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

правых частей осредненных уравнений содержат относительно небольшое количество слагаемых.

По сравнению с другими распространенными системами координат, применяемые нами астроцен-

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

не требуют ее разложения в ряд по степеням малого параметра. На временах порядка 10⁵-10⁷ лет

с помощью численного интегрирования осредненных уравнений изучена долговременная эволюция

планетных систем HD 12661, υ Andromedae, а также некоторых модельных систем. В рассмотренных

системах выявлены возможные вековые резонансы.

Ключевые слова: планетная задача, метод осреднения, метод Хори-Депри, гелиоцентрические коор-

динаты, астроцентрические координаты, вековые резонансы, гамильтониан, возмущающая функция,

ряды Пуассона, канонические элементы Пуанкаре, коэффициенты Лапласа.

DOI: 10.31857/S0320010820050058

1. ВВЕДЕНИЕ

(предполагаем, что барицентр системы неподви-

жен), сформулированный тип устойчивости пред-

Первые внесолнечные планеты были открыты в

ставляет собой специальный случай устойчивости

конце прошлого века. В настоящее время количе-

по Лагранжу¹ . Время существования устойчивой

ство известных экзопланет измеряется тысячами.

планетной системы сравнимо со временем суще-

По своим физическим и динамическим характери-

ствования центральной звезды.

стикам экзопланетные системы существенно отли-

Устойчивость планетных орбит обычно изуча-

чаются друг от друга, поэтому в этих системах мо-

ют аналитическими средствами, в основу которых

гут наблюдаться разные эволюционные картины.

положены идеи метода осреднения. Для получения

уравнений движения сначала выбирают систему

Под устойчивым движением планетной системы

координат. Наибольшее распространение получи-

обычно (см., например, Кузнецов, Холшевников,

ли координаты Якоби (Холшевников и др., 2001,

2009; Холшевников, Кузнецов, 2011) понимают

2002; Ли, Пил, 2003; Балуев, 2011; Перминов,

такое поведение оскулирующих эллипсов, при ко-

Кузнецов, 2015) и введенные в употребление Пу-

тором исключены тесные сближения и выбросы

анкаре (1965) канонические астроцентрические ко-

планет из системы. Вместе с тем предполагается,

ординаты (Красинский, 1973; Ласкар, Робютель,

что планетная конфигурация все время остается

близкой к компланарной, а также, что в любой паре

¹ Напомним, что частное решение системы обыкновенных

соседних эллипсов афелийное расстояние внутрен-

дифференциальных уравнений называется устойчивым по

ней орбиты всегда меньше перигелийного рассто-

Лагранжу, если оно остается определенным и равномерно

ограниченным при всех t ≥ t0, где t0 — начальный момент

яния внешней орбиты. Поскольку движение пла-

времени. Примером системы, устойчивой по Лагранжу,

нетной системы на протяжении всей ее эволюции

могут служить лагранжевы треугольные решения в общей

происходит в ограниченной области пространства

задаче трех тел. Точное определение устойчивости по

Лагранжу можно найти в руководстве Гребеникова (1976)

^*Электронный адрес: d.mikryukov@spbu.ru

и Демидовича (1967).

366

ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ ПЛАНЕТНОЙ СИСТЕМЫ

367

1995; Родригес, Гальярдо, 2005; Либер, Дилсейт,

Терминологическое замечание. В исследова-

2012; Либер, Сансотера, 2013; Робютель, Нидер-

ниях, посвященных динамике планетных систем,

ман, Пусс, 2016; Сансотера, Либер, 2019).

можно встретить разные трактовки термина “ве-

ковая эволюция планетной системы”; употребля-

Менее распространена система астроцентри-

ются также термины-синонимы “вековая модель”,

ческих координат Уинтнера (1967, §340-§341).

“вековое поведение” и т.п. В настоящей работе мы

Использование координат Уинтнера в планетной

под всеми этими терминами будем подразумевать

задаче с некоторых точек зрения представляет-

динамическую эволюцию, которая определяется

ся более естественным. Например, в отличие от

исключительно теорией первого порядка по пла-

координат Якоби, астроцентрические координаты

нетным массам. При построении этой теории все

Уинтнера не требуют разложения возмущающего

резонансные слагаемые из правых частей осред-

потенциала системы в ряд по степеням малого

ненной системы отбрасываются (см. разделы 4 и 5).

параметра. Далее, в канонических астроцентриче-

В разделах 2, 3 и 4 мы приводим основной

ских координатах Пуанкаре кинетическая энергия

системы выражается через барицентрические ско-

аналитический аппарат, используемый в работе (в

работах Микрюкова, 2016, 2018, содержится его

рости. Представление же кинетической энергии в

более подробное описание). В разделе 5 мы рас-

системе координат Уинтнера записывается с помо-

сматриваем вековую эволюцию некоторых модель-

щью астроцентрических скоростей (ср. формулу (7)

ных и реальных внесолнечных планетных систем.

в работе Ласкара и Робютеля, 1995, с формулой

В разделе 6 обсуждаются результаты.

(6₂) в §341 руководства Уинтнера, 1967). Благо-

даря этому, оскулирующие эллипсы в координатах

Уинтнера определяются более удобным образом.

2. ИСПОЛЬЗУЕМЫЕ ОБОЗНАЧЕНИЯ

Целью настоящего исследования является изу-

Для кеплеровых элементов будем использовать

чение устойчивости нерезонансных и слаборезо-

стандартный набор букв: a, e, i, M, g, Ω обозначают

нансных планетных систем на космогонических

соответственно большую полуось, эксцентриситет,

временах. Наше исследование выполняется в аст-

наклонение, среднюю аномалию, аргумент пери-

роцентрических координатах Уинтнера, так что по-

центра и долготу восходящего узла. Точка сверху

путно мы показываем, что эти координаты явля-

и черта сверху будут всюду означать соответствен-

ются достойной альтернативой ныне широко ис-

но дифференцирование по времени и комплексное

пользуемым координатам Якоби и каноническим

сопряжение. Мнимую единицу обозначим через i и

координатам Пуанкаре. Основным инструментом

положим Exp ϕ = exp iϕ для любого комплексно-

является метод осреднения. Мы записываем урав-

го ϕ. Нам понадобятся два символа суммирования

нения движения в удобной системе оскулирую-

щих элементов, после чего составляем уравнения

∑ N∑ ∑

∑

движения в средних элементах. Осредненная си-

(1)

стема интегрируется численно. Аналитический ап-

k=1

1≤j<k≤N

парат, необходимый для получения используемых

Второй символ (двойной) суммы будем применять

нами здесь осредненных уравнений, разработан и

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

детально описан автором в работах (Микрюков,

и k. В первом символе индекс может обозначаться

2016, 2018).

различными буквами. В формулах (1) натуральное

Мы будем рассматривать теорию первого и вто-

число N ≥ 2.

рого порядка по планетным массам. Поскольку

Все вычисления в работе (см. раздел 5) выпол-

теория второго порядка — довольно специальная

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

тема, требующая рассмотрения сложной проблемы

массы мы принимаем массу Солнца, за единицу

малых знаменателей² , мы решили разделить наше

длины — астрономическую единицу, а за единицу

исследование на две части. В настоящей работе мы

времени — сидерический земной год. Численное

рассмотрим теорию первого порядка — так назы-

значение постоянной тяготения G в этих единицах

ваемую вековую эволюцию планетной системы (см.

равно 4π².

также работы Кордани, 2004; Чазов, 2013). Мы

Остальные обозначения нам будет удобно вво-

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

дить по ходу изложения.

ных переменных осредненной системы. Мы пока-

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

для описания эволюции планетных систем весьма

3. ГАМИЛЬТОНИАН СИСТЕМЫ В

широкого класса. Следующую работу мы посвятим

АСТРОЦЕНТРИЧЕСКИХ КООРДИНАТАХ

теории второго порядка.

В евклидовом пространстве R³ рассмотрим

² Строго говоря, малые знаменатели появляются уже в

материальные точки Q₀, Q₁, Q₂, . . . , Q_N с массами

первом порядке, см. первый абзац раздела 5.

соответственно m₀, μm₁m₀, μm₂m₀, . . . , μm_N m₀.

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

368

МИКРЮКОВ

(

μm_k )

Gm₀m_k

Здесь N ≥ 2, безразмерный малый параметр μ

β_k = m_k

κ2k =

(κ_k > 0).

равен отношению массы самой массивной из

β_k

“планет” Q₁, . . . , Q_N к массе “Солнца” Q₀, а со-

Функция h зависит только от 6N компонент аст-

ответствующие планетам безразмерные множители

роцентрических векторов r_k, r_k (1 ≤ k ≤ N), так

m_k ≤ 1.

что дальнейшее исследование системы происходит

в 6N-мерном фазовом пространстве. Далее под

3.1. Астроцентрическая система координат

гамильтонианом мы будем всюду подразумевать

величину h.

Введем декартову инерциальную систему ко-

ординат с началом в точке O и зададим в этой

системе положение точек Q₀, Q₁, . . . , Q_N соответ-

3.2. Система оскулирующих элементов

ственно радиусами-векторами ρ₀, ρ₁, . . . , ρ_N . Аст-

роцентрическая система координат r определяется

Для большинства планетных систем характерна

ситуация, когда оскулирующие орбиты на протя-

(Уинтнер, 1967; Микрюков, 2016) равенствами

жении длительного интервала времени сохраняют

∑

малые значения наклонов и эксцентриситетов. В

r₀ =

ρ₀ +

m_jρ_j,

(2)

таких случаях исследование динамической эво-

люции удобно (Ласкар, Робютель, 1995) вести в

r_k = ρ_k - ρ₀, k = 1,2,... ,N,

канонических комплексных элементах Пуанкаре

√

где m = 1 + μm₁ + . . . + μm_N . Определитель пре-

P =X

Λ/2, p = -

(6)

образования (2) между абсолютными ρ и аст-

√

роцентрическими r координатами равен едини-

Q=Y

2Λ, q = -iQ,

це. Бесстолкновительное конфигурационное про-

Λ=βκ

√a, λ = M + g,

странство A₁ системы координат r получается из

R3N+3 удалением N(N - 1)/2 плоскостей

где

√

r_j = r_k,

1≤j <k≤N,

(3)

X =

2(1 - η)Exp g,

(7)

√

и N плоскостей

Y =

η(1 - cos i)/2Exp Ω,

√

r_j = 0,

1≤j ≤N.

(4)

η=

1-e²,

g = g + Ω.

Каждая из N(N + 1)/2 плоскостей (3), (4) имеет

Мы предполагаем, что в формулах (6) и (7) все

размерность 3N. Фазовое пространство A₂ = A₁ ×

переменные снабжены одним и тем же индексом j,

× R3N+3. Легко видеть, что каждое из множеств

равным номеру планеты (1 ≤ j ≤ N). Далее во всех

A₁ и A₂ открыто и связно и поэтому является

случаях, когда не могут возникнуть недоразумения,

областью.

мы будем опускать планетную индексацию.

Движение точек Q₀, Q₁, . . . , Q_N полностью

В отличие от кеплеровой части

определяется функцией Гамильтона H, представ-

∑

ляющей в координатах r полную механическую

β3sκ4s

h₀ = -

энергию системы. При составлении канонических

2Λ2s

уравнений движения удобно вместо H использо-

вать величину

гамильтониана (5), возмущающая функция h₁ за-

висит сложным образом от всех 6N перемен-

ных (6). В результате уравнения движения прини-

μm₀

мают форму

имеющую размерность квадрата скорости (см. так-

∂h₁

же работу Робютель и др., 2016). После исклю-

= -μ

p=μ

(8)

∂p

∂P

чения центра инерции функция h представляется

(Микрюков, 2016) стандартной суммой невозму-

∂h₁

Q=-μ∂h1

q=μ

щенной и малой возмущающей частей:

∂q

∂Q

h = h₀ + μh₁,

(5)

Λ=-μ∂h1

λ=ω+μ∂h1

где

∂λ

∂Λ

)

∑

⁽ ˙r2k

κ2k

где ω = ∂h₀/∂Λ = κ⁴β³Λ-3 = κa-3/2. Порядок

h₀ =

β_k

автономной системы (8) равен 6N.

|r_k|

(

)

Ниже мы будем иногда пользоваться векторны-

∑

Gm₀

r_j r_k

ми обозначениями

h₁ = -

m_jm_k

|r_j - r_k|

P, Q, Λ, p, q, λ.

(9)

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ ПЛАНЕТНОЙ СИСТЕМЫ

369

Положим также Γ = (P, Q, Λ, p, q). Каждый из

рядов (11) вещественны и безразмерны, однако

векторов (9) следует представлять себе как матри-

в отличие от C^jkℓvn, коэффициенты I_ℓvn сохраня-

цу 1 × N, а вектор медленных переменных Γ — как

ют фиксированное значение для любой планет-

матрицу 1 × 5N.

ной пары, поскольку являются постоянными ра-

циональными числами. Коэффициенты рядов (11)

загружены нами в базу данных Mendeley Data и

3.3. Разложение возмущающей функции

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

http://dx.doi.org/10.17632/3cb75grjz4.1.

Занумеруем планеты так, чтобы было a₁ < a₂ <

< ... < a_N. Для любых натуральных j и k та-

ких, что 1 ≤ j < k ≤ N, введем обозначения α_jk =

4. УРАВНЕНИЯ ДВИЖЕНИЯ В СРЕДНИХ

= a_j/a_k, R_jk = a_k|r_j - r_k|-1,

ЭЛЕМЕНТАХ

Gm₀m_jm_k

m_jm_k

Система уравнений (8) принадлежит к классу

R_jk = -

R_jk, V_jk = -

r_j r_k,

многочастотных систем с медленными и быстрыми

a_k

переменными (Гребеников, 1986). Поведение реше-

σ_jk = -Gm₀m_jm_kκ2kβ2k,

ний таких систем на длительных интервалах време-

m_jm_kκ2jκ2kβ_jβ_k

ни обычно исследуют с помощью метода осредне-

τ_jk = -

ния. Благодаря тому, что правые части осредненной

системы зависят лишь от медленных переменных,

ее решение изучать гораздо проще, чем решение

(XY )^ℓvjk = Xℓ1jXk2

Yℓ⁴k

Y v⁴k^,

точной системы. Общая теория метода осреднения

(nλ)_jk = n₁λ_j + n₂λ_k.

излагается в стандартных руководствах (Моисеев,

1969; Волосов, Моргунов, 1971; Митропольский,

В определениях величин (XY )^ℓvjk и (nλ)_jk буквами

1971; Боголюбов, Митропольский, 1974; Холшев-

ℓ, v, n обозначены целочисленные векторы. При

ников, 1985; Гребеников, 1986).

этом компоненты n₁ и n₂ вектора n могут при-

Важным на практике является случай, когда

нимать произвольные значения, а компоненты ℓ_s

исходные уравнения имеют вид автономной кано-

и v_s соответственно векторов ℓ и v — произволь-

нической системы, гамильтониан h которой, по-

ные неотрицательные значения. Теперь справедли-

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

во представление

аналитически от малого параметра μ:

∑

h₁ = R_jk + V_jk,

(10)

h = h₀ + μh₁ + μ²h₂ + ...

(12)

Обычным требованием осреднения гамильтоновых

в котором первая и вторая сумма называются

систем является сохранение канонической формы

соответственно главной и дополнительной частью

осредненных уравнений. Для случаев слабовоз-

возмущающей функции. Величины R_jk и V_jk раз-

мущенных консервативных систем, описываемых

лагаются в ряды Пуассона

гамильтонианом типа (12), новый (осредненный)

∑

гамильтониан H строится в виде

σ_jk

R_jk =

C^jkℓvn(XY )^ℓvjkExp (nλ)_jk,

(11)

Λ²

H = H₀ + μH₁ + μ²H₂ + ...

(13)

k K

τ_jk

∑

В нерезонансном случае коэффициенты H₀, H₁,

V_jk =

I_ℓvn(XY )^ℓvjkExp (nλ)_jk

Λ_jΛ_k

H₂,... ряда (13) должны зависеть только от мед-

K₁

ленных переменных.

с множествами суммирования K = {ℓ, v ∈ Z⁴, n ∈

Уравнения большинства планетных теорий за-

∈Z² :

даются гамильтонианом вида (12). В частности, для

нашего гамильтониана (5) справедливо

ℓ₁,ℓ₂,ℓ₃,ℓ₄,v₁,v₂,v₃,v₄ ≥ 0,

h(Γ, λ) = h₀(Λ) + μh₁(Γ, λ).

ℓ₁ + ℓ₂ + ℓ₃ + ℓ₄ - v₁ - v₂ - v₃ - v₄ + n₁ + n₂ = 0,

ℓ₃ + ℓ₄ + v₃ + v₄ = четное}

Мы будем рассматривать системы первого и второ-

го приближения, которые определяются соответ-

и K₁ = K\{ℓ,v ∈ Z⁴,n ∈ Z² : |n₁| + |n₂| = 0}. Ко-

ственно гамильтонианами

эффициенты

H(Γ) = H₀(Γ) + μH₁(Γ)

(14)

C^jkℓvn = C^jkℓ

(α_jk),

1ℓ2ℓ3ℓ4v1v2v3v4n1n2

I_ℓvn = Iℓ1ℓ₂ℓ₃ℓ₄v₁v₂v₃v₄n₁n₂

H(Γ) = H₀(Γ) + μH₁(Γ) + μ²H₂(Γ).

(15)

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

370

МИКРЮКОВ

Если планеты движутся вне областей острого ре-

Запишем уравнения движения в средних эле-

зонанса, то решения систем первого и второго

ментах

приближения дают общее представление о пове-

= -∂H/∂p,

p = ∂H/∂P,

(18)

дении решения исходной системы (8) на временах

Q=-∂H/∂q,

˙q = ∂H/∂Q,

порядка t ∼ μ-1 и t ∼ μ-2 соответственно.

Λ=0,

λ=∂H/∂Λ.

Метод Хори-Депри (Найфэ, 1976; Маркеев,

1978; Джакалья, 1979; Холшевников, 1985; Мор-

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

биделли, 2014; Шевченко, 2017) указывает эффек-

медленных переменных, поэтому из 6N уравнений

тивный путь вычисления коэффициентов осреднен-

(18) для нас представляют интерес только 4N

ного гамильтониана (13). Основная задача состоит

уравнений относительно P , Q, p, q. Далее, для

в построении производящей функции

определения P , Q, p, q достаточно решать лишь 2N

уравнений

T = μT₁ + μ²T₂ + ...,

(16)

˙P

= -∂H/∂p,

(19)

которая полностью определяет осредненный га-

мильтониан (13) и формулы перехода между стары-

Q=-∂H/∂q,

ми и новыми переменными. Для получения функций

поскольку координаты p = -

P и q = -i ^Q явля-

(14), (15) достаточно вычислить только первый

ются простейшими функциями импульсов P и Q

коэффициент T₁ = T₁(Γ, λ) асимптотического ряда

соответственно.

(16). Введем обозначения, с помощью которых

На практике уравнениям (19) удобно придать

можно будет выписать явный вид коэффициентов

несколько иной вид. Согласно основным разло-

H₀, H₁, H₂ отрезков разложений (14), (15).

жениям (11), исходный гамильтониан (5) зависит

Пусть некоторая функция f = f(Γ, λ) в некото-

от медленных переменных P , Q, p, q только через

рой области изменения своих аргументов разлага-

посредство

ется по компонентам вектора λ в N-кратный ряд

√

X =P

2/Λ, Y = Q/

2Λ,

(20)

Фурье

√

∑

X= ip

2/Λ,

Y = iq/

2Λ.

f (Γ, λ) =

f_k(Γ)Exp (kλ).

|k|≥0

Таким образом, можно написать

∑

h(P, Q, Λ, p, q, λ) = h(X,X, Y

Y ,Λ,λ).

(21)

Здесь |k| =

|k_s|, (kλ) =

k_sλ_s. Положим для

В результате для осредненного гамильтониана в

такой функции

(19) справедливо

∫2π

∫

2π

H(P, Q, Λ, p, q) = H(Λ, X,X, Y

Y ).

(22)

〈f〉 = (2π)^-N . . . f(Γ, λ)dλ₁ . . . dλ_N = f₀(Γ),

Мы предполагаем, что в равенствах (21), (22) все

аргументы являются N-мерными векторами.

∑

f_k(Γ)

Замечание 1. Равенства (21), (22) следует по-

f =

Exp (kλ),

нимать в физическом, а не математическом смысле:

i(kω)

|k|≥1

функциональные зависимости в правой и левой ча-

∑

сти (21) отличаются. Ситуация с (22) та же самая.

где (kω) = k_sω_s. Величина

f определяется в

Замечание 2. При разложении возмущающей

предположении отсутствия резонансов в системе.

функции (10) методом Ласкара и Робютеля (1995)

соотношения (20) рассматриваются как определе-

Коэффициенты H₀, H₁, H₂ гамильтонианов (14)

ния удобных безразмерных нормировок фазовых

и (15) находятся (Холшевников, 1985; Микрюков,

переменных P , Q, p, q. Поэтому при нахожде-

2018) по формулам

нии производных ∂H/∂X, ∂H/∂¯ гамильтониана,

H₀ = h₀, H₁ = 〈h₁〉,

(17)

определяемого правой частью (22), аргументы X_s,

H₂ = 〈{T₁,h₁ + H₁}〉/2.

X_s, Y_s

Y_s (1 ≤ s ≤ N) следует считать независи-

мыми друг от друга переменными (несмотря на то,

В равенствах (17) все функции зависят от средних

что, например, X_s и

X_s являются элементарными

элементов, причем T₁ =h1, а фигурные скобки

функциями друг друга).

означают скобку Пуассона по основной системе

Из последних двух равенств в (20) с учетом

элементов (6). Далее все встречающиеся функции

замечания 2 получаем

будут зависеть от средних элементов и, начиная с

√

∂H

2 ∂H

∂H

формул (17), мы будем опускать указание на эту

√

(23)

зависимость.

∂p

∂q

Λ∂X,

2Λ

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ ПЛАНЕТНОЙ СИСТЕМЫ

371

Далее, поскольку (средние) аналоги больших полу-

Отсюда

осей Λ неподвижны, справедливо

∑

√

∂H₁

∂〈R_js〉

∂〈R_sj 〉

Q₌Y˙

(28)

=X˙

Λ/2,

2Λ.

(24)

∂X

∂ X_s

j=1

j=s+1

Подстановка (23), (24) в (19) дает окончательно

Первая и вторая сумма в правой части (28) со-

-2i ∂H

-i ∂H

держат соответственно s - 1 и N - s слагаемых,

Λ ∂X,

2Λ

так что имеем N - 1 слагаемых вида ∂〈R_jk〉/∂^¯s,

Здесь осредненный гамильтониан определяется

причем в каждом из этих слагаемых либо j =

= s, либо k = s. Мы считаем, что первая и вторая

правой частью (22), при этом, согласно замеча-

сумма в правой части (28) равны нулю для случаев

нию 2, производные ∂H/∂X и ∂H/

Y находятся

соответственно s = 1 и s = N. Аналогично

по обычным правилам дифференцирования функ-

ций многих переменных.

∑

∂H₁

∂〈R_js〉

∂〈R_sj 〉

(29)

∂Y

Y_s

j=1

j=s+1

5. ИНТЕГРИРОВАНИЕ СИСТЕМЫ

ПЕРВОГО ПРИБЛИЖЕНИЯ

Комбинированием равенства R_jk = σ_jkR_jkΛ-2k с

соотношениями (25), (26), (28), (29) получаем на-

Запишем систему первого приближения

конец рабочее представление системы (25):

)

X_s = μ-2i∂H1 ,

(25)

(s-1∑

∑

Λ_s ∂Xs

X_s = μ

ξ_sju2sj +

ξ_sju¹

(30)

j=1

j=s+1

Y_s = μ-i∂H1 , 1≤s≤N.

2Λ_s

Y_s

)

(s-1∑

∑

Уравнения (25) задают вековую модель движения

Y_s = μ

η_sju4sj +

η_sju³

планетной системы. В нерезонансном случае эта

j=1

j=s+1

модель воспроизводит основные динамические ха-

Здесь

рактеристики эволюции орбит. При составлении

⎧

уравнений (25) информация о резонансах стирает-

⎨-2iσjs , если 1 ≤ j < s ≤ N,

Λ3s

ся, так что если последние присутствуют в системе,

ξ_sj =

⎩-2iσsj

если 1 ≤ s < j ≤ N,

решения систем (8) и (25) могут значительно от-

ΛsΛ²

личаться друг от друга. Информация о резонансах

восстанавливается после перехода к оскулирую-

η_sj = ξ_sj/4. Нетрудно проверить, что ξ_sj и η_sj име-

щим элементам по формулам замены перемен-

ют размерность обратного времени, а все осталь-

ных, см. формулы (41)-(45) в работе Микрюкова

ные величины в правых частях (30) безразмерны.

(2018).

Каждое слагаемое разложения величины H₁

имеет четную степень по совокупности переменных

Преобразуем правые части (25) к такому виду,

к которому уже удобно непосредственно применять

X,^X, Y

Y , поэтому согласно (26) правые части

машинные алгоритмы интегрирования. Для любых

(30) расположены лишь по нечетным степеням

натуральных j и k таких, что 1 ≤ j < k ≤ N, введем

этих величин (при дифференцировании H₁ степень

обозначения

каждого одночлена по X,

^X, Y,

Y понижается

на единицу). В нашей работе мы использовали

u1jk = ∂〈R_jk〉/∂Xj, u2kj = ∂〈R_jk〉/∂Xk,

(26)

разложение H₁ до десятой степени включительно.

u3jk = ∂〈R_jk〉/

Y_j, u4kj = ∂〈R_jk〉/

Y_k.

На практике построение уравнений (30) легко

выполняется даже в случае нескольких планет в

Здесь

системе. По существу, дело сводится лишь к одно-

∑

〈R_jk〉 =

C^jkℓv(XY )^ℓvjk,

(27)

типному формированию отрезков разложений (26)

и их подстановке в (30).

K₀

Ниже мы рассматриваем результаты численного

причем C^jkℓv = Cjkℓv00, а множество K₀ получается

интегрирования системы (30) на примерах двупла-

в результате пересечения K с плоскостью {ℓ,v ∈

нетных и трехпланетных систем. Интегрирование

∈ Z⁴,n ∈ Z² : |n₁| + |n₂| = 0}. Далее, среднее зна-

мы выполняли реализованным на языке C++ яв-

чение дополнительной части возмущающей функ-

ным одношаговым семистадийным методом Рунге-

ции равно нулю, так что

Кутты шестого порядка точности. Использованная

нами таблица Бутчера для коэффициентов метода

∑

H₁ = 〈h₁〉 =

R_jk + V_jk

〈R_jk〉.

приведена в руководстве Хайрера, Нерсетта и Ван-

нера (1990, глава II.6, таблица 6.1).

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

372

МИКРЮКОВ

Таблица 1. Начальные значения средних элементов для модельной двупланетной системы

Планета

a (а.е.)

Внутренняя планета (“Венера”)

0.723315

0.00676

3.39448^◦

54.8978^◦

76.6243^◦

Внешняя планета (“Земля”)

1.000027

0.01672

0.00262^◦

287.9199^◦

175.0383^◦

Таблица 2. Границы изменения средних эксцентриситетов и средних наклонений в модельной двупланетной системе.

Интервал интегрирования составляет 2 × 10⁶ лет

Планета

e_min

e_max

i_min

i_max

Внутренняя планета (“Венера”)

0.00564

0.0188

0.6367^◦

3.3942^◦

Внешняя планета (“Земля”)

0.00836

0.01701

0.005^◦

2.7597^◦

5.1. Модельная двупланетная система

ключевых факторов, поскольку ею обеспечивается

длительное поддержание на поверхности планеты

Эффективность методов теории возмущений в

более менее постоянных условий.

планетной задаче обусловлена тем, что в уравне-

ния движения можно простым и удобным образом

Изучение движения землеподобных внесолнеч-

ных планет разумно начинать на модельных приме-

ввести малый параметр μ. В результате во всех

рах. Мы рассмотрим сейчас динамическую эволю-

планетных теориях получаются уравнения типа (8)

цию модельной двупланетной системы, которая со-

и (18), к которым можно применять асимптотиче-

стоит из звезды солнечной массы и двух обращаю-

скую теорию и делать определенные выводы отно-

щихся около нее планет с массами и орбитальными

сительно эволюции планетных орбит. Большинство

исследований сосредоточено на орбитальной дина-

элементами Венеры и Земли. Массовые параметры

мике газовых гигантов — наиболее массивных пла-

μ = 3.0404 × 10-6,

(31)

нет в системе. В таких работах параметр μ обычно

m₁ = 0.8051, m₂ = 1

имеет порядок 10-3-10-2. Например, в работе

Кузнецова и Холшевникова (2006) рассматрива-

и начальные значения оскулирующих элементов

ется эволюция двупланетной системы Солнце-

(эпоха JD 2459000.5) мы взяли из ежегодника (Ко-

Юпитер-Сатурн³ , в которой принято μ = 10-3.

четова и др., 2019). Начальные значения средних

Интерес к изучению поведения систем, состо-

элементов требуют расчета (Волосов, Моргунов,

ящих из планет-гигантов, представляется есте-

1971; Микрюков, 2018); мы приводим их в табл. 1.

ственным. Во-первых, количество подтвержден-

Интегрирование системы (30) выполнялось на

ных экзопланет земного типа на порядки мень-

2 × 10⁶ лет вперед. На всем временном интервале

ше числа открытых внесолнечных планет-гигантов.

наблюдается близкое к периодическому поведение

Во-вторых, динамика именно больших планет ре-

(средних) эксцентриситетов и наклонов в системе.

шает вопрос об устойчивости системы в целом.

На рис. 1 изображена вековая эволюция этих эле-

Экзопланеты земного типа сложнее обнаружить

ментов на интервале 5 × 10⁵ лет. Мы видим, что ко-

из-за их относительно небольшого размера и ма-

лебания совершаются в противофазе, что являет-

лой массы. Между тем исследование орбиталь-

ся классическим следствием сохранения интеграла

ной динамики этих объектов также заслуживает

площадей. Период P_e колебаний эксцентрисите-

внимания. Согласно современным представлени-

тов составляет примерно 131 × 10³ лет, а период

ям такие планеты являются наиболее благоприят-

ным местом для возникновения жизни. Устойчи-

P_i колебаний наклонов — около 106 × 10³ лет. В

вость планетной орбиты является здесь одним из

табл. 2 приведены предельные значения, принима-

емые средними наклонами и эксцентриситетами в

³ Под двупланетной системой “Солнце-Юпитер-Сатурн”

данной системе.

Кузнецов и Холшевников (2006) понимают такую мо-

Одной из важных динамических характеристик

дельную двупланетную систему, в которой около звезды

солнечной массы обращаются две планеты с массами и

планетной системы является вековой резонанс

орбитальными параметрами Юпитера и Сатурна.

(Мюррей, Дермотт,

2009). Вековой резонанс

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ ПЛАНЕТНОЙ СИСТЕМЫ

373

0.020

3.5

0.018

3.0

0.016

2.5

0.014

2.0

0.012

1.5

0.010

1.0

0.008

0.006

0.5

0.004

10⁵

t (годы)

Рис. 1. Эволюция средних эксцентриситетов (слева) и наклонений (справа) в модельной двупланетной системе согласно

вековой модели движения (30). Сплошная линия соответствует внешней (более тяжелой) планете.

500

10⁵

t (годы)

Рис. 2. Иллюстрация векового резонанса в модельной двупланетной системе. Разность Δg = g1 - g2 является близкой

к периодической функцией времени.

может отсутствовать, однако это довольно рас-

кривых повторяется. На основании рис. 4 за-

пространенное явление. Чаще всего под вековым

ключаем, что эволюция лагранжевых элементов

резонансом понимают такое поведение планет,

sin

cos Ω и sin

sin Ω имеет также близкий к

= g₁ - g₂ оказывается

периодическому характер. Период изменения этих

близкой к периодической функцией времени (см.,

элементов совпадает с P_i, поскольку на плоскостях

например, Ли, Пил, 2003). Обычно Δg колеблется

(

)

i_s

либо около нуля, либо около π. На основании

sin

cos Ω_s, sin

sin Ω_s

, s = 1,

2, движение

рис. 2 заключаем, что в отсутствие всех остальных

тел в Солнечной системе, движение Венеры и

изображающей точки совершает один цикл-оборот

Земли происходит в вековом резонансе. Амплитуда

(см. рис. 4) за P_i ≈ 106 × 10³ лет.

колебаний Δg равна 48◦, причем период этих

колебаний совпадает с P_e.

5.2. Двупланетные системы HD 12661

На рис. 3 мы изображаем вековую эволю-

и v Andromedae

цию лагранжевых элементов e cos g и e sin g

В разобранном модельном примере малый па-

на интервале времени

2 × 10⁶

лет. На обеих

раметр μ имеет порядок 10-6. Рассмотрим теперь

плоскостях (e_s cos g_s, e_s sin g_s), s = 1,

2, движе-

примеры реальных планетных систем, в которых

ние изображающей точки оказывается близким

μ ∼ 10-3.

к периодическому. Продолжительность одного

Солнцеподобная звезда HD 12661 имеет спек-

цикла-периода составляет 10P_e ≈ 1.31 × 10⁶ лет,

тральный класс G6 и находится в созвездии Ов-

после которого качественное поведение обеих

на на расстоянии примерно 120 световых лет от

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

374

МИКРЮКОВ

0.02

0.01

-0.01

-0.02

-0.01

0.01

0.02

-0.02

-0.01

0.01

0.02

e₁ cos g₁

e₂ cos g₂

Рис. 3. Эволюция лагранжевых элементов es cos gs, es sin gs (s = 1, 2) на интервале времени 2 × 10⁶ лет. Модельная

двупланетная система.

0.03

0.02

0.01

0.02

0.010

0.005

0.010

0.015

Рис. 4. Эволюция лагранжевых элементов sin

cos Ωs, sin

sin Ωs (s = 1, 2) на интервале 5 × 10⁵ лет. Модельная

двупланетная система.

Солнца. Вокруг звезды обращаются, по крайней

В табл. 3 и 4 приводятся массовые параметры и

мере, две планеты (газовые гиганты HD 12661 b

начальные значения средних элементов для обеих

и HD 12661 c). В отличие от HD 12661, звезду

систем. На основании графиков, представленных

v Andromedae можно наблюдать невооруженным

на рис. 5, мы делаем вывод, что период колебаний

глазом, поскольку, имея примерно такие же фи-

средних эксцентриситетов в системах HD 12661 и

зические характеристики, она находится ближе к

v Andromedae составляет около 29000 и 9570 лет

нам почти в три раза. В настоящее время из-

соответственно. В табл. 5 приведены максималь-

вестно о существовании четырех планет в системе

ные и минимальные значения.

v Andromedae. Мы ограничимся изучением дви-

Согласно графикам, изображенным на рис. 6,

жения только двух внутренних планет v And c

в обеих системах наблюдается вековой резонанс.

иv Andd,посколькумассапланет v Andbиv Ande

существенно меньше.

Таблица

3. Значения массовых параметров μ,

Системы HD 12661 и v Andromedae — доволь-

m₀, m₁ и m₂ для двупланетных систем HD 12661

но известные экзопланетные объекты. Система

и v Andromedae

v Andromedae была первой системой, в которой

была обнаружена третья планета. На примере

m₀

m₁

m₂

системы HD 12661 было впервые установлено

(Ли, Пил, 2003), что две реальные внесолнечные

HD 12661

2.079 × 10-3

1.07

0.7854

планеты могут двигаться таким образом, при ко-

тором разность долгот перицентров Δg совершает

v And

3 × 10-3

1.28

0.4875

колебания около значения π.

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ ПЛАНЕТНОЙ СИСТЕМЫ

375

0.35

0.30

0.25

0.20

0.15

0.10

0.05

10⁴

10⁵

10⁴

t (годы)

Рис. 5. Эволюция средних эксцентриситетов в двупланетных системах HD 12661 (слева) и v Andromedae (справа).

Сплошная линия соответствует наиболее массивной планете (см. табл. 3). Результаты интегрирования системы первого

приближения.

Амплитуда колебаний Δg в системах HD 12661

мкнутой кривую. Если Δg колеблется около π, то

и v Andromedae равна 56^◦ и 43^◦ соответственно.

на плоскостях (32) кривая располагается слева от

Для иллюстрации векового резонанса часто (см.,

оси ординат. Если Δg совершает колебания около

например, Родригес, Гальярдо, 2005) используют

нуля, то на плоскостях (32) получается фигура,

переменные e_s cos Δg, e_s sin Δg (s = 1, 2). Если

расположенная справа от оси ординат. (Кривая

планетная система испытывает вековой резонанс,

находится целиком по одну сторону от оси ординат

то периоды изменения e₁, e₂, Δg совпадают. В

тогда и только тогда, когда амплитуда Δg меньше

результате на обеих плоскостях

90^◦.) На рис. 7 мы изображаем вековую эволюцию

систем HD 12661 и v Andromedae на плоскостях

(e_s cos Δg, e_s sin Δg), s = 1, 2,

(32)

(e₁ cos Δg, e₁ sin Δg) и (e₂ cos Δg, e₂ sin Δg) соот-

изображающая точка описывает близкую к за-

ветственно.

Таблица 4. Начальные значения средних элементов для

5.3. Модельная трехпланетная система

систем HD 12661 и v Andromedae

Добавление третьей планеты в систему может

значительно изменить характер движения первых

Системы

a (а.е.)

двух планет. В подразделе 5.1 мы исследовали ве-

HD 12661 b

0.821

0.34

0.5^◦

290.6^◦

6^◦

ковую эволюцию модельной двупланетной систе-

мы, которая определяет движение Венеры и Земли

HD 12661 c

2.855

0.066

1.5^◦

94.1^◦

8^◦

вокруг Солнца в предположении отсутствия возму-

щений со стороны всех остальных тел Солнечной

v And c

0.832

0.224

2^◦

229.4^◦

3^◦

системы. Рассмотрим теперь поведение этой си-

v And d

2.53

0.267

5^◦

123.5^◦

135^◦

стемы с учетом влияния Марса. Именно, образуем

такую модельную трехпланетную систему, в кото-

рой около звезды солнечной массы обращаются

Таблица 5. Границы изменения средних эксцентрисите-

три планеты с массами и орбитальными парамет-

тов в системах HD 12661 и v Andromedae. Интервал

рами, соответствующими Венере, Земле и Марсу.

интегрирования для систем HD 12661 и v Andromedae

Планеты в данной системе будем далее называть

составляет соответственно 10⁵ и 4 × 10⁴ лет

“Венерой” (внутрення планета, соответствует Ве-

нере), “Землей” (промежуточная планета, соответ-

Системы

e_min

e_max

ствует Земле) и “Марсом” (внешняя планета, со-

ответствует Марсу). Использованием кавычек мы

HD 12661 b

0.1506

0.3405

хотим подчеркнуть, что рассматриваем поведение

Венеры, Земли и Марса в предположении отсут-

HD 12661 c

0.0636

0.2624

ствия всех остальных тел в Солнечной системе.

v And c

0.053

0.256

В табл. 1 и 6 мы приводим начальные значения

средних элементов для всех трех планет. Наряду с

v And d

0.259

0.29

(31) примем m₃ = 0.1061.

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

376

МИКРЮКОВ

240

220

200

180

160

140

120

10⁴

10⁵

500

10⁴

t (годы)

Рис. 6. Вековой резонанс в системах HD 12661 (слева) и v Andromedae (справа).

0.3

0.2

0.1

-0.1

-0.2

-0.3

-0.2

-0.36

-0.32

-0.28

-0.24

-0.20

-0.16

0.20

0.22

0.24

0.26

0.28

0.30

e₁ cos Δg

e₂ cos Δg

Рис. 7. Вековое поведение систем HD 12661 (слева) и v Andromedae (справа) на плоскостях (e1 cos Δg, e1 sin Δg) и

(e2 cos Δg, e2 sin Δg) соответственно. См. подробности в тексте.

Интегрирование выполнялось на 3 × 10⁶ лет

ситетов приобретает еще более сложный харак-

вперед. На основании рис. 8 заключаем, что влия-

тер (см. рис. 9). В разобранном выше вариан-

ние более легкого “Марса” существенно изменяет

те (при a₃ = 1.5238) прослеживается систематич-

поведение эксцентриситетов “Венеры” и “Земли”

ность в изменении всех эксцентриситетов, причем

(ср. с рис. 1). Эволюция эксцентриситетов e₁, e₂, e₃

период этой систематичности, согласно графику на

в системе имеет сложный характер, однако “Марс”

рис. 8, можно оценить примерно в 735 × 10³ лет.

сохраняет самую вытянутую орбиту в течение всего

Теперь же (см. рис. 9) эксцентриситеты “Венеры” и

интервала интегрирования. Эксцентриситет e₃ из-

“Земли” изменяются практически непредсказуемо.

меняется в пределах 0.082-0.104.

На всем интервале интегрирования продолжает

Рассмотрим теперь случай, когда орбита “Мар-

сохраняться верным неравенство e₃ ≫ e₁, e₂. Заме-

са” располагается ближе к орбите “Земли”. Иссле-

тим, что для изучения вековой эволюции планетных

дуем вариант α₁₂ = α₂₃. Для этого в табл. 1 и 6

систем с близкими орбитами более эффективны-

сохраним все значения прежними, кроме одно-

ми могут оказаться специальные методы разло-

го: положим a₃ = 1.382. Поведение эксцентри-

жения возмущающей функции (см. работу Вашко-

вьяк и др., 2015).

больше

Рассмотрим, наконец, случай, когда a₃

Таблица 6. Начальные значения средних элементов

значения a₃ = 1.5238. (Значение a₃ = 1.5238 близ-

третьей планеты (“Марса”) в модельной трехпланетной

ко к текущему оскулирующему значению большой

системе

полуоси Марса, см. табл. 6.) Положим a₃ = 1.84,

а все остальные числовые значения в табл. 1 и 6

a (а.е.)

вновь оставим прежними. При a₃ = 1.84 получаем

α₂₃ ≈ 0.5434, что примерно соответствует отноше-

1.5238

0.09345

1.8479^◦

286.59^◦

49.5^◦

нию больших полуосей Юпитера и Сатурна. Дан-

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

378

МИКРЮКОВ

0.10

0.08

0.06

0.04

0.02

10⁶

t (годы)

Рис. 10. Эволюция средних эксцентриситетов в модифицированной модельной трехпланетной системе. Точечная линия

соответствует “Венере”, сплошная — “Земле”, а пунктирная — “Марсу”. Для среднего значения большой полуоси

орбиты “Марса” положено a3 = 1.84 а.е.

ный вариант интегрировался вперед на 8 × 10⁶ лет,

большими, чтобы в уравнениях (33) производные

см. рис. 10. На этот раз все эксцентриситеты ведут

∂H₂/∂X и ∂H₂/

Y давали заметный вклад в эво-

себя более понятным образом. Изменение e₁, e₂,

люцию X и Y . Если большие полуоси в системе

e₃ становится близким к периодическому. Согласно

далеки от резонансных значений, то при столь ма-

рис. 10 период изменения e₁, e₂, e₃ равен примерно

лом μ результаты интегрирования систем первого и

3.57 × 10⁶ лет. Интересно, что теперь наблюдаются

второго приближения будут слабо отличаться друг

промежутки времени, в течение которых e₁, e₂ > e₃.

от друга.

По сравнению с двумя предыдущими модельными

В настоящей работе мы рассматривали эволю-

вариантами, эксцентриситет “Марса” изменяется

цию нерезонансных модельных и реальных планет-

в более широких пределах 0.017-0.098. Эксцен-

ных конфигураций на временных шкалах порядка

триситеты “Венеры” и “Земли” могут превышать

10⁵-10⁷ лет. Как правило, в каждом случае для

значение 0.034.

интегрирования мы выбирали такой промежуток

времени, на котором может быть выявлено ка-

6. ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

чественное поведение основных орбитальных ха-

Чем меньше значение параметра μ в основных

рактеристик исследуемой системы (период коле-

уравнениях (18), тем выше эффективность асимп-

баний эксцентриситетов, наклонов, разности Δg,

тотической теории. Напишем систему второго при-

вековые резонансы и т.п.). Стоит отметить, что в

ближения

нерезонансных системах согласно общей теории

(

)

-2i

∂H

∂H₂

метода осреднения (Боголюбов, Митропольский,

(33)

1974; Волосов, Моргунов, 1971; Митропольский,

∂ X+μ2

∂X¯

(

)

1971) длина временного интервала, на котором

-i

∂H₁

∂H₂

осредненное решение сохраняет главные черты ис-

+μ²

2Λ

тинного решения, существенно зависит от малого

параметра μ. Чем меньше значение μ, тем длиннее

Пусть, например, при изучении взаимодействия

должен быть этот интервал.

экзопланет земного типа мы имеем μ ∼ 10-6. В

результате μ² ∼ 10-12 и мы заключаем, что коэф-

Значения μ ∼ 10-6 и μ ∼ 10-3 примерно со-

фициенты разложения H₂ должны быть достаточно ответствуют планетам земного типа со скалистой

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ ПЛАНЕТНОЙ СИСТЕМЫ

379

Таблица 7. Число N слагаемых в разложении 〈R_jk〉

группы Солнечной системы может существенно за-

в зависимости от d. На практике всегда берется d > 0

висеть от устойчивости движения планет-гигантов.

Таким образом, устойчивости движения планет

земной группы благоприятствует или устойчивость

системы планет-гигантов, или их отсутствие.

N 1 9 61 261

878

2446

5982

13182

26807

В заключение отметим, что употребление ком-

плексных переменных (6) приводит к относительно

небольшому количеству слагаемых в разложении

Таблица 8. Число N^′ слагаемых в разложении каждой

H₁ (см. работы Ласкар, Робютель, 1995; Микрю-

из четырех величин u1jk, u2kj, u3jk, u4kj в зависимости

ков, Холшевников, 2016). На практике бесконеч-

от d^′ = d - 1

ный восьмикратный степенной ряд (27) усекается

условием

d^′

ℓ₁ + ℓ₂ + ℓ₃ + ℓ₄ + v₁ + v₂ + v₃ + v₄ ≤ d,

N^′

122

472

1452

3804

8844

18744

где d — четное натуральное. В результате макси-

мальная степень d^′ разложений величин (26) ока-

зывается равной d^′ = d - 1. В табл. 7 мы приводим

поверхностью и газовым гигантам. Не стоит за-

число N слагаемых в разложении (27) в зависи-

бывать о существовании планет с промежуточ-

мости от d, а в табл. 8 — число N^′ слагаемых в

ными физическими характеристиками. Например,

каждом из четырех разложений (26) в зависимости

масса планеты Kepler-11 f (система Kepler-11) по

от d^′. В случае, например, N = 2 и d = 4 система

текущим оценкам превышает земную лишь в два

(30) принимает вид

раза, хотя по плотности Kepler-11 f сопоставима с

X₁ = μξ₁₂u¹¹²,

Y₁ = μη₁₂u³¹²,

(34)

Сатурном. Обратим также внимание на то, что в

нашем исследовании мы были сосредоточены в ос-

X₂ = μξ₂₁u²²¹,

Y₂ = μη₂₁u⁴²¹,

новном на поведении эксцентриситетов и движении

перицентров в планетной системе.

при этом H₁ содержит 61 слагаемое (ср. с раз-

ложением 〈R₁₂〉, приведенным на стр. 215 работы

В подразделах 5.1 и 5.3 мы изучали динами-

Ласкара и Робютеля, 1995) и генерирует в правой

ческую эволюцию систем, образованных только

части каждого из уравнений (34) по 22 слагаемых.

планетами земного типа. В подразделе 5.2 ис-

следовалась эволюция конфигураций, состоящих

только из планет-гигантов. На основании резуль-

татов подразделов 5.1 и 5.3 можно сделать вывод,

что если система образована только планетами

Автор выражает благодарность К.В. Холшев-

земного типа, то в течение длительного интервала

никову за руководство работой, В.Ш. Шайдулину

времени она может сохранять устойчивость. Для

за помощь в оформлении рисунков, а также ре-

устойчивости такой системы большие планеты мо-

цензенту за полезные замечания, способствовав-

гут оказаться не нужными.

шие улучшению работы. Все вычисления в работе

производились с помощью оборудования вычис-

Более близким к реальному будет случай, когда

лительного центра научного парка СПбГУ. Работа

в систему входят планеты обоих типов. Лучше

выполнена при финансовой поддержке РНФ (грант

всего в этом убеждает пример нашей Солнечной

19-72-10023).

системы, которая содержит одинаковое количество

планет земного типа и газовых гигантов. Пример

Солнечной системы показывает, что формирование

СПИСОК ЛИТЕРАТУРЫ

землеподобных экзопланет около звезд должно

1. Балуев (R.V. Baluev), Celest. Mech. Dynam. Astron.

быть не менее распространенным явлением, чем

111, 235 (2011).

образование больших планет. Движение систем,

2. Боголюбов Н.Н., Митропольский Ю.А., Асимпто-

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

тические методы в теории нелинейных коле-

массами, является предметом отдельного и особого

баний (М.: Наука, 1974).

исследования. Например, Ласкаром (1994) было

3. Вашковьяк М.А., Вашковьяк С.Н., Емельянов

выполнено интегрирование осредненных уравне-

Н.В., Астрон. вестник

49,

208

(2015)

ний для модели Солнечной системы, состоящей

[M.A. Vashkov’yak et al., Solar System. Res.

из Солнца и восьми основных планет. Интегриро-

49, 191 (2015)].

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

4. Волосов В.М., Моргунов Б.И., Метод осредне-

вперед и назад. В результате этого исследования

ния в теории нелинейных колебательных си-

был сделан вывод, что устойчивость планет земной

стем (М.: МГУ, 1971).

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020

380

МИКРЮКОВ

Гребеников Е.А., Качественная небесная меха-

23.

Митропольский Ю.А., Метод усреднения в нели-

ника (М.: Наука, сб. Справочное руководство

нейной механике (К.: Наукова думка, 1971).

по небесной механике и астродинамике, под

24.

Моисеев Н.Н., Асимптотические методы нели-

ред. Г.Н. Дубошина, 1976), с. 788.

нейной механики (М.: Наука, 1969).

Гребеников Е.А., Метод усреднения в приклад-

25.

Морбиделли А., Современная небесная механи-

ных задачах (М.: Наука, 1986).

ка. Аспекты динамики Солнечной системы (М.:

Демидович Б.П., Лекции по математической

ИКИ, 2014).

теории устойчивости (М.: Наука, 1967).

26.

Мюррей К., Дермотт С., Динамика Солнечной

Джакалья Г.Е.О., Методы теории возмущений

системы (М.: Физматлит, 2009).

для нелинейных систем (М.: Наука, 1979).

27.

Найфэ А.Х., Методы возмущений (М.: Мир,

Кордани (B. Cordani), Celest. Mech. Dynam. Astron.

1976).

89, 165 (2004).

28.

Перминов А.С., Кузнецов Э.Д., Астрон. вестник 49,

10.

Кочетова О.М., Кузнецов В.Б., Медведев Ю.Д.,

469 (2015)

[A.S. Perminov, E.D. Kuznetsov, Solar

Чернетенко Ю.А., Шор В.А., Эфемериды малых

System. Res. 49, 430 (2015)].

планет на 2020 год (СПб.: ИПА РАН, 2019).

29.

Пуанкаре А., Лекции по небесной механике (М.:

11.

Красинский Г.А., Основные уравнения планет-

Наука, 1965).

ной теории (М.: Наука, сб. Малые Планеты, под

30.

Робютель, Нидерман, Пусс (P. Robutel,

ред. Н.С. Самойловой-Яхонтовой, 1973), с. 81.

L. Niederman, and A. Pousse), Comp. Appl. Math.

12.

Кузнецов Э.Д., Холшевников К.В., Астрон. вестник

35, 675 (2016).

40, 263 (2006) [E.D. Kuznetsov, K.V. Kholshevnikov,

31.

Родригес, Гальярдо (A. Rodr´ıguez and T. Gallardo),

Solar System. Res. 40, 239 (2006)].

Astrophys. J. 628, 1006 (2005).

13.

Кузнецов Э.Д., Холшевников К.В., Астрон. вестник

32.

Сансотера, Либер (M. Sansottera and A.-S. Libert),

43, 230 (2009) [E.D. Kuznetsov, K.V. Kholshevnikov,

Celest. Mech. Dynam. Astron. 131, 38 (2019).

Solar System. Res. 43, 220 (2009)].

33.

Уинтнер А., Аналитические основы небесной ме-

14.

Ласкар (J. Laskar), Astron. Astrophys., 287, L9-L12

ханики (М.: Наука, 1967).

(1994).

34.

Хайрер Э., Нерсетт С., Ваннер Г., Решение

15.

Ласкар, Робютель (J. Laskar and P. Robutel), Celest.

обыкновенных дифференциальных уравнений.

Mech. Dynam. Astron. 62, 193 (1995).

Нежесткие задачи (М.: Мир, 1990).

16.

Ли, Пил (M.H. Lee and S.J. Peale), Astrophys. J.

35.

Холшевников К.В., Асимптотические методы

592, 1201 (2003).

небесной механики (Л.: ЛГУ, 1985).

17.

Либер, Дилсейт (A.-S. Libert and N. Delsate),

36.

Холшевников К.В., Греб А.В., Кузнецов Э.Д., Аст-

MNRAS 422, 2725 (2012).

рон. вестник 35, 267 (2001) [K.V. Kholshevnikov et

18.

Либер, Сансотера (A.-S. Libert and M. Sansottera),

al., Solar System. Res. 35, 243 (2001)].

Celest. Mech. Dynam. Astron. 117, 149 (2013).

37.

Холшевников К.В., Греб А.В., Кузнецов Э.Д., Аст-

19.

Маркеев А.П., Точки либрации в небесной меха-

рон. вестник 36, 75 (2002) [K.V. Kholshevnikov et al.,

нике и космодинамике (М.: Наука, 1978).

Solar System. Res. 36, 68 (2002)].

20.

Микрюков Д.В., Письма в Астрон. журн. 42, 611

38.

Холшевников, Кузнецов (K.V. Kholshevnikov and

(2016)

[D.V. Mikryukov, Astron. Lett. 42,

555

E. D. Kuznetsov), Celest. Mech.Dynam. Astron. 109,

(2016)].

201 (2011).

21.

Микрюков Д.В., Письма в Астрон. журн. 44, 361

39.

Чазов В.В., Астрон. вестник

47,

112

(2013)

(2018)

[D.V. Mikryukov, Astron. Lett. 44,

337

[V.V. Chazov, Solar System. Res. 47, 99 (2013)].

(2018)].

22.

Микрюков Д.В., Холшевников К.В., Письма в

40.

Шевченко И.И., The Lidov-Kozai Effect —

Астрон. журн. 42, 302 (2016)

[D.V. Mikryukov,

Applications in Exoplanet Research and

K.V. Kholshevnikov, Astron. Lett. 42, 268 (2016)].

Dynamical Astronomy (Springer, 2017).

ПИСЬМА В АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 46

№5

2020