Астрономический журнал, 2020, T. 97, № 6, стр. 443-475

Трехмерное численное моделирование структуры течения в асинхронном поляре CD Ind в приближении смещенного дипольного магнитного поля белого карлика

А. В. Соболев 1*, А. Г. Жилкин 1, Д. В. Бисикало 1, Д. А. Х. Бакли 2

1 Институт астрономии РАН
Москва, Россия

2 South African Astronomical Observatory
Cape, Town, South Africa

* E-mail: asobolev@inasan.ru

Поступила в редакцию 24.12.2019
После доработки 02.03.2020
Принята к публикации 02.03.2020

Полный текст (PDF)

Аннотация

В работе исследована структура течения в асинхронном поляре CD Ind в приближении смещенного диполя. Результаты расчетов позволили выделить такие особенности системы, как дрейф горячих пятен на поверхности белого карлика, зависимость структуры течения от фазы спин-орбитального периода биений, процесс переключения магнитных полюсов. Для исследования структуры течения мы использовали трехмерную численную МГД модель, основанную на приближении модифицированной магнитной гидродинамики. Численные расчеты проведены для десяти фаз периода биений при постоянном темпе массообмена. Кроме того, для более детального исследования процесса переключения полюсов проведены дополнительные серии расчетов в соответствующих диапазонах фазы спин-орбитального периода биений. Показано, что зоны энерговыделения в течение спин-орбитального периода смещаются по долготе в среднем на $20^\circ $, что соответствует величине 0.05 фазы орбитального периода. Установлено, что процесс переключения магнитных полюсов происходит за время, не превышающее 0.1 спин-орбитального периода биений. В середине этого процесса аккреция осуществляется на оба полюса с одинаковой интенсивностью, а струя вещества выглядит в форме арки. По результатам расчетов построены синтетические кривые блеска в оптическом диапазоне спектра.

1. ВВЕДЕНИЕ

Поляры [1] являются подклассом магнитных катаклизмических переменных звезд и представляют собой тесные двойные системы [2], состоящие, как правило, из белого карлика (аккретор, первичный компонент) с магнитным полем более 10 МГс и маломассивной звезды главной последовательности (донор, вторичный компонент). Оптическое (а также инфракрасное) излучение поляров характеризуется наличием значительной степени поляризации, что и отражено в их названии. Впервые этот эффект обнаружен Тапиа в 1976 г. у объекта AM Her [3]. Магнитное поле аккретора оказывает сильное влияние на обмен веществом между компонентами двойной системы. Донор, заполнив свою полость Роша, начинает терять вещество из оболочки через внутреннюю точку Лагранжа [4, 5], которое затем захватывается магнитным полем белого карлика. При этом вместо аккреционного диска формируется коллимированная струя, в которой вещество перетекает к магнитным полюсам аккретора. Другим эффектом, связанным с наличием сильного магнитного поля, является синхронизация периодов собственного вращения компонентов с их орбитальным периодом [6, 7].

Однако, когда белый карлик взрывается как новая, у него возникает некоторый асинхронизм (порядка 1–2%) собственного вращения относительно орбитального периода двойной системы [8]. В связи с этим принято разделять поляры на синхронные и асинхронные (системы типа BY Cam). В настоящее время известно около 100 поляров с орбитальным периодом от 80 мин до 8 ч, из них асинхронными являются всего 4 системы: V1500 Cyg, V1432 Aql, BY Cam и CD Ind. Для асинхронных поляров вводится такая характеристика, как период биений ${{P}_{{{\text{beat}}}}}$, определяемый из следующего соотношения:

(1)
$\frac{1}{{{{P}_{{{\text{beat}}}}}}} = \tfrac{1}{{{{P}_{{{\text{spin}}}}}}} - \tfrac{1}{{{{P}_{{{\text{orb}}}}}}},$
где ${{P}_{{{\text{spin}}}}}$ и ${{P}_{{{\text{orb}}}}}$ – периоды собственного вращения белого карлика относительно наблюдателя и орбитальный период двойной системы соответственно. В течение одного периода биений аккреционный поток от вторичной звезды вращается вокруг магнитного поля белого карлика. Таким образом, ориентация магнитного поля белого карлика относительно донора повторяется за один период биений. В зависимости от ориентации магнитного поля (и его геометрии) поток аккреции будет следовать за различными линиями магнитного поля и переключаться с одного полюса на другой [9].

Объект CD Ind (также известный как EUVE J2115–586 и RX J2115.7–5840) был обнаружен в 1996 г. [10] как вероятный кандидат в магнитные катаклизмические переменные звезды при обработке данных с ультрафиолетового телескопа EUVE, а также в 1997 г. [11] как рентгеновский источник. Кроме того, авторы работы [11] впервые идентифицировали систему как асинхронный поляр с самым коротким периодом биений, в котором напряженность магнитного поля белого карлика составляет 11 ± 2 МГс. Выполненные в дальнейшем оптическая поляриметрия и рентгеновские наблюдения [12, 13] подтвердили принадлежность CD Ind к классу асинхронных поляров, при этом период его биений оказался равным одной неделе. В 2017 г. после обработки данных, полученных за 10 лет наблюдений поляра, были уточнены значения орбитального периода (110.8 мин) и периода вращения белого карлика (109.6 мин) [14]. Определенный таким образом по формуле (1) период биения составил 7.03 дня, что делает на сегодня CD Ind асинхронным поляром с самым коротким периодом биения.

Возможно, что еще одной особенностью данного поляра является наличие смещенного дипольного магнитного поля. Предположение о такой его конфигурации было высказано в статье [15] на основе данных циклотронного картирования и численного моделирования в рамках приближения квазичастиц [1622]. В качестве другого объяснения наблюдаемой геометрии аккреции авторы статьи [15] допускают возможное наличие магнитного поля более высокого порядка. Например, можно рассматривать суперпозицию дипольного и квадрупольного полей. Однако в рамках данной работы мы исследуем только гипотезу смещенного диполя. В дальнейшем мы планируем выполнить численные расчеты с более сложной структурой магнитного поля.

Стоит отметить, что конфигурация магнитного поля со смещенным диполем свойственна некоторым другим объектам, например химически пекулярным звездам [23], а также планетам Солнечной системы – Урану [24] и, по последним данным, Юпитеру [25].

В данной работе на основе результатов трехмерного численного моделирования, а также с учетом интерпретации наблюдательных данных мы исследуем структуру течения в асинхронном поляре CD Ind в предположении смещенного дипольного магнитного поля белого карлика. Статья организована следующим образом. Во втором разделе приводится описание используемой нами численной модели. В третьем разделе представлены результаты численных расчетов. В Заключении кратко обсуждаются основные выводы по работе.

2. ОПИСАНИЕ МОДЕЛИ

Более детальное представление о структуре течения в полярах можно получить на основе совместного использования наблюдательных данных и результатов численного моделирования. В работах [4, 26, 27] (см. также монографию [5]) была разработана самосогласованная трехмерная численная модель для расчета структуры течений в тесных двойных системах с учетом магнитного поля. Для моделирования мы использовали модифицированную систему уравнений магнитной гидродинамики, которая позволяет описать все основные динамические эффекты, связанные с магнитным полем [28]. В численной модели были учтены процессы радиационного нагрева и охлаждения, нагрева за счет диссипации токов, а также диффузия магнитного поля. При этом формирование и последующая эволюция аккреционного потока происходят естественным образом в результате процесса массопереноса вещества через внутреннюю точку Лагранжа. Мы успешно применяли этот подход для моделирования структуры течения в полярах и промежуточных полярах [4, 9, 27, 2939] (см. также монографию [5]).

Для исследования структуры течения в системе CD Ind были выполнены трехмерные численные расчеты для десяти фаз спин-орбитального периода ${{P}_{{{\text{beat}}}}}$. В нашей модели использованы следующие параметры системы. Звезда-донор (красный карлик) имеет массу ${{M}_{{\text{d}}}} = 0.21{{M}_{ \odot }}$ и эффективную температуру ${{T}_{{\text{d}}}} = 3200$ К. Звезда-аккретор (белый карлик) имеет массу ${{M}_{{\text{a}}}} = 0.7{{M}_{ \odot }}$, радиус ${{R}_{{\text{a}}}} = 0.014{{R}_{ \odot }}$ и эффективную температуру ${{T}_{{\text{a}}}} = $ $ = 12{\kern 1pt} {\kern 1pt} 000$ К. Период обращения двойной системы ${{P}_{{{\text{orb}}}}} = 1.846$ ч, а большая полуось орбиты $A = $ $ = 0.735{{R}_{ \odot }}$. Индукция магнитного поля в районе северного магнитного полюса задавалась равной ${{B}_{{\text{a}}}} = 11$ МГс.

Численное моделирование структуры течения проводилось в неинерциальной системе отсчета, вращающейся вместе с двойной системой с угловой скоростью $\Omega = 2\pi {\text{/}}{{P}_{{{\text{orb}}}}}$ вокруг ее центра масс. Поле сил, действующих на вещество в такой системе отсчета, определяется потенциалом Роша

(2)
$\Phi = - \frac{{G{{M}_{{\text{a}}}}}}{{\left| {\vec {r} - {{{\vec {r}}}_{{\text{a}}}}} \right|}} - \tfrac{{G{{M}_{{\text{d}}}}}}{{\left| {\vec {r} - {{{\vec {r}}}_{{\text{d}}}}} \right|}} - \tfrac{1}{2}{{[\vec {\Omega } \times (\vec {r} - {{\vec {r}}_{{\text{c}}}})]}^{2}},$
где $G$ – гравитационная постоянная, а радиусы-векторы ${{\vec {r}}_{{\text{a}}}}$, ${{\vec {r}}_{{\text{d}}}}$ и ${{\vec {r}}_{{\text{c}}}}$ определяют положение центра аккретора, центра донора и центра масс двойной системы. Первое и второе слагаемые в выражении (2) описывают гравитационные потенциалы аккретора и донора. Последнее слагаемое описывает центробежный потенциал относительно центра масс. В выбранной системе отсчета используется декартова система координат ($x$, $y$, $z$), начало которой расположено в центре аккретора (${{\vec {r}}_{{\text{a}}}} = 0$). Ось $x$ совпадает с линией, соединяющей центры донора и аккретора, при этом ее положительные значения направлены в сторону, противоположную донору. Отрицательное направление оси $x$ указывает на нулевую фазу орбитального периода, при этом возрастание орбитальной фазы происходит в направлении по часовой стрелке. Вращение системы в этом случае осуществляется против часовой стрелки. Центр донора расположен на расстоянии $A$ от начала координат вдоль оси $x$, поэтому вектор ${{\vec {r}}_{{\text{d}}}} = ( - A,0,0)$. Ось $y$ перпендикулярна оси $x$ и лежит в орбитальной плоскости системы. Ось $z$ направлена вдоль оси вращения двойной системы: $\vec {\Omega } = (0,0,\Omega )$.

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

(3)
${{\vec {B}}_{*}} = \frac{\mu }{{{{R}^{3}}}}[3(\vec {d} \cdot \vec {n})\vec {n} - \vec {d}],$
где $\mu $ – магнитный момент, вектор $\vec {R} = \vec {r} - {{\vec {r}}_{{{\text{dip}}}}}$, $\vec {n} = \vec {R}{\text{/}}R$ – единичный вектор, направленный из центра аккретора в точку наблюдения поля. Радиус-вектор ${{\vec {r}}_{{{\text{dip}}}}}$ соответствует центру диполя. В нашей модели центр диполя был смещен по $z$ на расстояние половины радиуса звезды ниже орбитальной плоскости поляра. Поэтому вектор ${{\vec {r}}_{{{\text{dip}}}}} = (0,0, - {{R}_{{\text{a}}}}{\text{/}}2)$. Единичный вектор $\vec {d}$ определяет ось симметрии диполя. Компоненты вектора $\vec {d}$ в декартовой системе координат могут быть записаны в виде
(4)
${{d}_{x}} = sin\theta cos\varphi ,\quad {{d}_{y}} = sin\theta sin\phi ,\quad {{d}_{z}} = cos\theta ,$
где углы $\theta $, $\varphi $ определяют ориентацию магнитной оси в пространстве. Угол $\theta $ отсчитывается от северного географического полюса аккретора, а $\phi $ – от положительного направления оси $x$ в сторону вращения системы. Вектор магнитного момента $\vec {\mu } = \mu \vec {d}$, а магнитный момент аккретора вычислялся по формуле $\mu = {{B}_{{\text{a}}}}R_{{\text{a}}}^{3}{\text{/}}2$. Поскольку мы рассматриваем систему с асинхронным вращением аккретора, азимутальный угол $\phi $ будет зависеть от времени. Ориентация магнитного диполя, соответствующая нулевой фазе периода биений, задана следующими значениями углов: ${{\theta }_{0}} = 70^\circ $, ${{\phi }_{0}} = 90^\circ $ [15]. Далее при переходе от одной фазы спин-орбитального периода к другой угол $\phi $ возрастает в направлении вращения двойной системы против часовой стрелки. Угол $\theta $ при этом остается неизменным.

Заметим, что магнитное поле ${{\vec {B}}_{*}}$, задаваемое формулой (3), является потенциальным, $\nabla \times {{\vec {B}}_{*}} = $ = 0. Кроме того, поскольку в данной работе структуру течения мы моделировали отдельно для разных фаз спин-орбитального периода ${{P}_{{{\text{beat}}}}}$ без учета собственного вращения аккретора, то магнитное поле ${{\vec {B}}_{*}}$ в каждом варианте расчета можно считать стационарным, $\partial {{\vec {B}}_{*}}{\text{/}}\partial t = 0$. Эти обстоятельства позволяют его частично исключить из соответствующих уравнений, описывающих структуру МГД течения [5, 4042]. Этот прием в численной модели удобно использовать, чтобы избежать в процессе расчета накопления ошибок при операциях с большими числами. Для этого полное магнитное поле $\vec {B}$ можно представить в виде суперпозиции поля аккретора ${{\vec {B}}_{*}}$ и поля $\vec {b}$, индуцированного электрическими токами в аккреционной струе и оболочке двойной системы: $\vec {B} = {{\vec {B}}_{*}} + \vec {b}$. При этом в разностной схеме вычисляется только возмущение магнитного поля $\vec {b}$. Такое расщепление магнитного поля часто используется при моделировании МГД аккреции (см., напр., [4346]).

Структура течения в магнитных тесных двойных системах может быть описана следующей системой уравнений [4, 5]:

(5)
$\frac{{\partial \rho }}{{\partial t}} + \nabla \cdot (\rho \vec {v}) = 0,$
(6)
$\begin{gathered} \frac{{\partial{ \vec {v}}}}{{\partial t}} + (\vec {v} \cdot \nabla )\vec {v} = - \frac{{\nabla P}}{\rho } - \frac{{\vec {b} \times (\nabla \times \vec {b})}}{{4\pi \rho }} - \\ \, - \nabla \Phi + 2(\vec {v} \times \vec {\Omega }) - \frac{{{{{\vec {v}}}_{ \bot }}}}{{{{t}_{w}}}}, \\ \end{gathered} $
(7)
$\frac{{\partial{ \vec {b}}}}{{\partial t}} = \nabla \times [\vec {v} \times \vec {b} + \vec {v} \times {{\vec {B}}_{*}} - {{\eta }_{w}}(\nabla \times \vec {b})],$
(8)
$\begin{gathered} \rho \left[ {\frac{{\partial \varepsilon }}{{\partial t}} + (\vec {v} \cdot \nabla )\varepsilon } \right] = \\ = \; - {\kern 1pt} P(\nabla \cdot \vec {v}) + {{n}^{2}}(\Gamma - \Lambda ) + \frac{{{{\eta }_{w}}}}{{4\pi }}{{(\nabla \times \vec {b})}^{2}}, \\ \end{gathered} $
где $\rho $ – плотность, $\vec {v}$ – скорость, $P$ – давление, $\varepsilon $ – удельная внутренняя энергия газа, $n$ – концентрация, ${{\eta }_{w}}$ – коэффициент магнитной вязкости. Слагаемое $2(\vec {v} \times \vec {\Omega })$ в уравнении движения (6) описывает силу Кориолиса. Плотность, внутренняя энергия и давление связаны между собой уравнением состояния идеального газа
(9)
$P = (\gamma - 1)\rho \varepsilon ,$
где $\gamma = 5{\text{/}}3$ – показатель адиабаты.

В уравнении энергии (8) учтены эффекты радиационного нагрева и охлаждения [4750], а также нагрев за счет диссипации токов. В нашей численной модели используются линейные аппроксимации функций нагрева $\Gamma $ и охлаждения $\Lambda $ от температуры $T$ в окрестности ее равновесного значения ${{T}_{*}} = 8278$ К, соответствующей эффективной температуре горячего компонента ${{T}_{{\text{a}}}} = 12\,000$ К:

(10)
$\Gamma = {{\Gamma }_{*}} + \Gamma _{*}^{'}(T - {{T}_{*}}),\quad \Lambda = {{\Lambda }_{*}} + \Lambda _{*}^{'}(T - {{T}_{*}}),$
где ${{\Gamma }_{*}} = {{\Lambda }_{*}} = 3.46 \times {{10}^{{ - 25}}}$ (эрг см3)/с, $\Gamma _{*}^{'} = - 8.39 \times $ $ \times \;{{10}^{{ - 29}}}$ (эрг см3)/(с К), $\Lambda _{*}^{'} = 3.05 \times {{10}^{{ - 29}}}$ (эрг см3)/(с К). При вычислении разности $\Gamma - \Lambda $ величины ${{\Gamma }_{ * }}$ и ${{\Lambda }_{ * }}$ сокращаются.

В основе нашей модели лежит приближение модифицированной магнитной гидродинамики [27, 4, 5 ], которое подробно описано в работе [28]. Это приближение соответствует МГД в присутствии очень сильных внешних магнитных полей с учетом волновой альфвеновской турбулентности при малых магнитных числах Рейнольдса (${{R}_{{\text{m}}}} \ll 1$) [51]. Динамика плазмы в сильном внешнем магнитном поле характеризуется относительно медленным усредненным движением вдоль магнитных силовых линий, дрейфом частиц под действием внешних сил поперек магнитных силовых линий и распространением с очень большими скоростями альфвеновских и магнитозвуковых волн. Поскольку быстрые МГД волны за характерное динамическое время могут много раз пересекать область потока, взаимодействовать между собой и формировать турбулентный каскад [52], мы можем исследовать усредненную картину течения по аналогии с МГД турбулентностью [5355].

В уравнении движения (6) последнее слагаемое описывает эффективную электромагнитную силу, действующую на плазму со стороны магнитного поля аккретора, которая влияет на компонент скорости плазмы, перпендикулярный (символ $ \bot $) магнитным силовым линиям [36, 5658]. Выражение для этой силы является аналогом силы трения между компонентами плазмы, состоящей из нескольких видов частиц (см., напр., [59, 60]). Иными словами, сильное внешнее магнитное поле играет роль эффективной жидкости, с которой взаимодействует плазма.

Шкала времени релаксации для поперечного компонента скорости определяется выражением:

(11)
${{t}_{w}} = \frac{{4\pi \rho {{\eta }_{w}}}}{{B_{*}^{2}}}.$
Здесь коэффициент магнитной турбулентной вязкости
(12)
${{\eta }_{w}} = {{\alpha }_{w}}\frac{{{{l}_{w}}{{B}_{*}}}}{{\sqrt {4\pi \rho } }},$
где ${{\alpha }_{w}}$ – безразмерный коэффициент, характеризующий эффективность волновой турбулентности, а ${{l}_{w}} = {{B}_{*}}{\text{/}}\left| {\nabla {{B}_{*}}} \right|$ – характерная пространственная шкала волновых пульсаций. При расчетах значение ${{\alpha }_{w}}$ в модели принималось равным 1/3, что соответствует изотропной турбулентности [61]. Последний член в уравнении энергии (8) в рассматриваемом приближении волновой МГД турбулентности может быть представлен в виде следующего выражения:

(13)
$\frac{{{{\eta }_{w}}}}{{4\pi }}{{(\nabla \times \vec {b})}^{2}}\frac{{\rho \mathop {{\vec {v}}}\nolimits_ \bot ^2 }}{{{{t}_{w}}}}.$

В численной модели использованы следующие начальные и граничные условия. В оболочке звезды-донора нормальный компонент скорости по отношению к поверхности ${{v}_{n}}$ задавался равным локальной скорости звука ${{c}_{s}}$, соответствующей эффективной температуре донора ${{T}_{{\text{d}}}} = 3200$ К. Плотность газа в оболочке донора $\rho ({{L}_{1}})$ определяется из выражения для темпа массообмена через внутреннюю точку Лагранжа ${{L}_{1}}$:

(14)
$\dot {M} = \rho ({{L}_{1}}){{v}_{n}}S,$
где площадь сечения струи $S$ из донора вычисляется по формуле [5, 62]:
(15)
$S = \frac{{\pi c_{s}^{2}}}{{4{{\Omega }^{2}}}}{{g}_{y}}(q){{g}_{z}}(q),$
${{g}_{y}}(q) = 1.113$ и ${{g}_{z}}(q) = 1.036$ – безразмерные параметры, зависящие от отношения масс $q = {{M}_{{\text{d}}}}{\text{/}}{{M}_{{\text{a}}}} = 0.3$ компонентов двойной системы и определяющие соответственно большую и малую полуоси эллиптического сечения струи. Аккретор был определен сферой радиусом $0.014A$, на границе которой были заданы условия свободного втекания. На внешних границах вычислительной области были заданы постоянные граничные условия: плотность ${{\rho }_{b}} = {{10}^{{ - 8}}}\rho ({{L}_{1}})$, температура ${{T}_{b}} = {{T}_{*}}$, магнитное поле ${{\vec {b}}_{b}} = 0$. Для скорости $\mathop {\vec {v}}\nolimits_b $ были заданы условия свободного истечения: когда скорость направлена наружу, использовались симметричные граничные условия $\partial {{\vec {v}}_{b}}{\text{/}}\partial{ \vec {n}} = 0$, а когда скорость направлена внутрь, использовалось условие ${{\vec {v}}_{b}} = 0$. Начальные условия в вычислительной области: плотность ${{\rho }_{0}} = {{10}^{{ - 8}}}\rho ({{L}_{1}})$, температура ${{T}_{0}} = {{T}_{*}}$, скорость ${{\vec {v}}_{0}} = 0$ и магнитное поле ${{\vec {b}}_{0}} = 0$.

Для проведения расчетов использовался трехмерный численный код на основе разностной схемы годуновского типа повышенного порядка точности для уравнений магнитной гидродинамики. Соответствующая разностная схема подробно описана в нашей недавней работе [39]. Задача решалась в расчетной области $ - 2.0 \leqslant x{\text{/}}A \leqslant 1.0$, $ - 1.5 \leqslant y{\text{/}}A \leqslant 1.5$ и $ - 0.75 \leqslant z{\text{/}}A \leqslant 0.75$ с числом ячеек $256 \times 256 \times 128$ и неравномерным шагом сетки, экспоненциально уменьшающимся к центру аккретора. Такая расчетная область полностью включает в себя полости Роша аккретора и донора.

3. РЕЗУЛЬТАТЫ РАСЧЕТОВ

3.1. Структура течения

В данной работе мы представляем десять вариантов численного расчета структуры течения в системе CD Ind, соответствующих различным фазам спин-орбитального периода, с шагом фазе 0.1. Во всех моделях расчет проводился до выхода на квазистационарный режим, который определялся примерным (с точностью до 1%) постоянством полной массы вещества в расчетной области. Как правило, выход на квазистационарный режим достигался за время порядка 1–2 орбитальных периодов. При проведении вычислений темп массообмена между компонентами считался постоянной величиной, равной $\dot {M} = {{10}^{{ - 9}}}{{M}_{ \odot }}{\text{/год}}$ [15].

Результаты трехмерных численных расчетов представлены на рис. 1–5. Оттенками цвета показаны изоповерхности десятичного логарифма плотности в единицах $\rho ({{L}_{1}})$ для значений $ - 1$, $ - 2$, $ - 3$ и $ - 4$. Поверхность звезды-аккретора представлена в виде светлой сферы. Магнитным силовым линиям соответствуют линии со стрелками. Показаны также ось вращения (синяя вертикальная линия) и магнитная ось (зеленая наклонная линия) аккретора. Вид системы представлен с точки зрения наблюдателя с Земли с учетом того, что двойная система вращается вокруг общего центра масс, а угол наклона плоскости орбиты $i = 70^\circ $. Для иллюстрации структуры течения на каждой фазе периода биений было выделено такое орбитальное положение системы, при котором хорошо видны характерные детали течения. Так, для фаз 0.0–0.4 выбрана орбитальная фаза 0.7. В этом положении можно видеть разделение струи между полюсами аккретора на фазе 0.1 (режим переключения струи между полюсами будет описан ниже), а также проследить за изменением траектории движения вещества при асинхронном вращении аккретора. Помимо течения вещества из точки ${{L}_{1}}$, на рис. 1–5 видно скопление материи в плоскости магнитного экватора аккретора и дополнительный поток вещества из общей оболочки двойной системы, который создает вторичную, значительно более слабую по сравнению с основной, зону аккреции. Поэтому на всех фазах спин-орбитального периода аккрецию вещества всегда можно наблюдать в окрестности обоих магнитных полюсов, тогда как основная струя, за исключением моментов переключения, течет только на один полюс.

Рис. 1.

Результат трехмерного численного моделирования структуры течения вещества в системе CD Ind для фаз 0.0 (верхняя панель) и 0.1 (нижняя панель) периода биений. Показаны изоповерхности десятичного логарифма плотности (цвет) в единицах $\rho ({{L}_{1}})$. Поверхность аккретора представлена в виде белой сферы. Линии со стрелками соответствуют силовым линиям магнитного поля. Вертикальная линия, проходящая через аккретор, совпадает с его осью вращения, а наклонная линия – с его магнитной осью. Вид системы показан с точки зрения земного наблюдателя (система вращается вокруг общего центра масс, угол наклона плоскости орбиты равен $i = {{70}^{ \circ }}$).

Рис. 2.

То же, что и на рис. 1, но для фаз 0.2 (верхняя панель) и 0.3 (нижняя панель).

Рис. 3.

То же, что и на рис. 1, но для фаз 0.4 (верхняя панель) и 0.5 (нижняя панель).

Рис. 4.

То же, что и на рис. 1, но для фаз 0.6 (верхняя панель) и 0.7 (нижняя панель).

Рис. 5.

То же, что и на рис. 1, но для фаз 0.8 (верхняя панель) и 0.9 (нижняя панель).

Формирование пояса вещества в плоскости магнитного экватора объясняется действием эффекта “магнитной ловушки”. Поэтому этот элемент структуры в какой-то мере можно считать аналогом радиационных поясов Земли. Электромагнитная сила со стороны магнитного поля аккретора (последнее слагаемое в правой части выражения (6)) действует только в поперечном к силовым линиям направлении. Поэтому в области магнитосферы, где поле сильное, плазма не может свободно течь вдоль плоскости магнитного экватора. Ускорение в продольном к магнитным линиям направлении мало, поскольку в плоскости магнитного экватора вектор гравитационной силы аккретора ортогонален магнитному полю. В результате вблизи плоскости магнитного экватора происходит постепенное накопление вещества. Численные эксперименты показывают, что этот эффект проявляется при использовании разностных схем высокого порядка точности (в нашем это схема Роу–Ошера–Эйнфельдта для уравнений магнитной гидродинамики [39]).

В силу значительного наклона магнитного диполя (70°) и его смещения ниже плоскости экватора белого карлика на половину радиуса струя вещества из оболочки донора незначительно отклоняется от орбитальной плоскости системы. На нулевой фазе периода биений направляющий вектор магнитного диполя $\vec {d}$ перпендикулярен прямой линии, соединяющей донор и аккретор. При этом северный магнитный полюс расположен слева и ближе к экватору белого карлика, если смотреть на систему со стороны донора. Меньшее расстояние до плоскости экватора от данного полюса по сравнению с южным магнитным полюсом означает, что он при асинхронном вращении аккретора должен подходить ближе к донору (к внутренней точке Лагранжа ${{L}_{1}}$), а значит, вызывать рост темпа аккреции. Соответствующая зона аккреции должна быть более активной. Это должно проявляться в увеличении потока излучения на кривой блеска при наблюдении северного магнитного полюса.

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

Плотность вещества струи вблизи внутренней точки Лагранжа ${{L}_{1}}$ для темпа массообмена ${{10}^{{ - 9}}}{{M}_{ \odot }}{\text{/год}}$ составляет $\rho ({{L}_{1}}) = 1.92 \times {{10}^{{ - 7}}}$ г/см3. Вблизи фаз периода биений 0.0 и 0.5, т.е. когда зона аккреции расположена под значительным углом к донору, под действием силы Кориолиса увеличивается длина баллистической части траектории струи по отношению к ее магнитной части. Напротив, в случае максимального сближения магнитных полюсов с точкой ${{L}_{1}}$ на фазах 0.3 и 0.7 вещество струи практически сразу же захватывается магнитным полем первичного компонента и устремляется на его поверхность. Поэтому баллистическая часть траектории оказывается наиболее короткой. Кроме того, из-за неоднородности течения его структура в поперечном сечении имеет слоистый характер с увеличением плотности к центру. Это означает, что внешние, менее плотные слои будут быстрее захватываться магнитным полем аккретора, тогда как внутренние будут успевать продвинуться дальше вдоль баллистической части траектории. Поскольку менее плотные части аккреционного потока останавливаются магнитным полем белого карлика раньше, чем более плотные, взаимодействие вещества аккреционного потока из оболочки донора с магнитным полем аккретора должно приводить к формированию иерархической структуры магнитосферы. Учет такой структуры магнитосферы представляется важным, поскольку это может повлиять на результаты анализа и интерпретации наблюдательных данных. Кроме того, такая неоднородность вещества струи ведет к ее незначительному расширению, что и наблюдается в представленных трехмерных расчетах. Особый интерес представляют фазы периода биений, на которых происходит переключение течения с одного магнитного полюса на другой. На центральной стадии этого процесса, когда струя движется одновременно на оба полюса, ее ширина значительно увеличивается за счет действия суперпозиции магнитных сил от соответствующих полюсов.

3.2. Горячие пятна на поверхности аккретора

Чтобы продемонстрировать эффект дрейфа горячих пятен в течение периода биений, мы рассчитали распределение температуры на поверхности аккретора. Для этого мы использовали методику, описанную в работе [45] и применявшуюся нами ранее при анализе структуры течения в асинхронном поляре BY Cam [9].

Аккретор в спокойном состоянии в отсутствие аккреции вещества имеет эффективную температуру ${{T}_{{\text{a}}}}$, так что поток излучения с его поверхности равен ${{\sigma }_{{{\text{SB}}}}}T_{{\text{a}}}^{4}$, где ${{\sigma }_{{{\text{SB}}}}}$ – постоянная Стефана–Больцмана. Если предположить, что при падении вещества на поверхность аккретора в излучение переходит его тепловая и кинетическая энергии, то плотность потока энергии аккрецирующего вещества в точке $\mathop {\vec {R}}\nolimits_{\text{a}} $ поверхности аккретора определяется выражением:

(16)
$q(\mathop {\vec {R}}\nolimits_{\text{a}} ) = - \rho \vec {n} \cdot \vec {v}\left( {\varepsilon + \tfrac{{{{{\vec {v}}}^{2}}}}{2} + \tfrac{P}{\rho }} \right),$
где $\vec {n}$ – вектор нормали к поверхности. Знак “минус” определяется тем, что на поверхности аккретора нормальный компонент скорости падающего вещества отрицателен, $\vec {n} \cdot \vec {v} < 0$. Поэтому в выражении (16) величина плотности потока энергии $q({{\vec {R}}_{{\text{a}}}})$ оказывается положительной.

В нашей модели считается, что излучение из зон аккреции, так же как и от остальной поверхности звезды, имеет чернотельный характер. Тогда можно утверждать, что локальная эффективная температура $T({{\vec {R}}_{{\text{a}}}})$ в данной точке поверхности должна удовлетворять следующему соотношению:

(17)
${{\sigma }_{{{\text{SB}}}}}T{{({{\vec {R}}_{{\text{a}}}})}^{4}} = {{\sigma }_{{{\text{SB}}}}}T_{{\text{a}}}^{4} + q({{\vec {R}}_{{\text{a}}}}).$
Заметим, что в этой формуле мы не учитываем наличие ударной волны в основании аккреционной колонки, за фронтом которой температура вещества может существенно возрасти. Учет ударной волны необходим для расчета потока рентгеновского излучения. Однако для наших целей этим эффектом можно пренебречь, поскольку на характеристики потока излучения в оптическом диапазоне ударная волна практически никак не влияет. Кроме того, мы хотим определить местоположения зон аккреции, а не исследовать подробно их физические характеристики.

Полученные распределения температуры на поверхности аккретора для всех фаз спин-орбитального периода показаны на рис. 6–10. Левая диаграмма каждого рисунка соответствует полушарию белого карлика, содержащего северный магнитный полюс. При этом считается, что магнитный полюс находится точно на срединном меридиане полушария. На правой диаграмме представлено противоположное полушарие белого карлика, содержащее южный магнитный полюс. Для наглядности указаны положения северного и южного магнитных полюсов в виде “шариков”: синий обозначает северный полюс, красный – южный. Светлые области соответствуют эффективной температуре аккретора ${{T}_{{\text{a}}}} = 12\,000$ К в спокойном состоянии, а темные области – зонам энерговыделения.

Рис. 6.

Распределение температуры на поверхности аккретора для фаз 0.0 (верхняя панель) и 0.1 (нижняя панель) периода биений. Темные области соответствуют зонам энерговыделения. “Шарики” указывают на положения северного (вверху) и южного (внизу) магнитных полюсов. Левое полушарие соответствует области, включающей северный магнитный полюс, правое полушарие – южный магнитный полюс.

Рис. 7.

То же, что и на рис. 6, но для фаз 0.2 (верхняя панель) и 0.3 (нижняя панель) периода биений.

Рис. 8.

То же, что и на рис. 6, но для фаз 0.4 (верхняя панель) и 0.5 (нижняя панель) периода биений.

Рис. 9.

То же, что и на рис. 6, но для фаз 0.6 (верхняя панель) и 0.7 (нижняя панель) периода биений.

Рис. 10.

То же, что и на рис. 6, но для фаз 0.8 (верхняя панель) и 0.9 (нижняя панель) периода биений.

Рисунки показывают, что зоны энерговыделения сосредоточены около магнитных полюсов. При этом в течение периода биений положение этих зон относительно каждого полюса меняется. Горячие пятна около северного и южного магнитных полюсов с увеличением фазы спин-орбитального периода движутся справа налево вдоль линии долготы белого карлика. Шаг сетки на рис. 6–10 по широте и долготе составляет 5°. Анализируя приведенные результаты, можно заключить, что горячие пятна смещаются по долготе относительно северного и южного магнитных полюсов в среднем на 20°, что соответствует величине 0.05 фазы орбитального периода. Смещение по широте при этом практически не происходит.

На фазах 0.0, 0.7, 0.8 и 0.9, когда на северный полюс аккрецирует вещество не из струи, а из околозвездной оболочки, мы видим зону симметрично распределенной аккреции вокруг полюса. Как только полюс становится активным и принимает вещество из струи, картина аккреции содержит явно выраженную область повышенного энерговыделения. На фазе 0.1, соответствующей середине процесса первого переключения полюсов, горячее пятно является довольно размытым и растянутым по широте белого карлика. Однако уже на фазах 0.2–0.4, когда северный полюс расположен близко к точке ${{L}_{1}}$, зона основного энерговыделения сужается и становится явно выраженной. Размер этой зоны в указанном положении аккретора составляет порядка $10^\circ \times 5^\circ $, увеличиваясь до $10^\circ \times 10^\circ $ на фазе 0.4. Однако здесь стоит отметить, что несмотря на рост площади зоны энерговыделения на фазе 0.4, максимум потока излучения приходится на фазу 0.3 (см. болометрические кривые блеска ниже). Как видно из рис. 6–10, длительность активного состояния аккреции северного магнитного полюса составляет около $30\% $ периода биений.

Аналогичная картина наблюдается для южного магнитного полюса. Здесь сформированная зона энерговыделения соответствует фазам 0.7–0.0 периода биений. Ее размер, а также его изменение полностью повторяют соответствующие параметры для северного полюса. Максимум потока излучения при этом приходится на фазу 0.9.

Отметим еще одну особенность распределения горячих пятен. Северное пятно формируется ниже, а южное – выше соответствующего магнитного полюса по широте. Такую конфигурацию зон аккреции можно объяснить геометрией магнитных силовых линий и их ориентацией относительно экваториальной плоскости двойной системы. На баллистическом участке движение аккреционной струи вещества из точки ${{L}_{1}}$ происходит в экваториальной плоскости. Вблизи аккретора струя под действием магнитного поля начинает отклоняться от орбитальной плоскости вниз, двигаясь вдоль магнитных силовых линий. В результате формирование горячего пятна происходит выше (южнее) северного магнитного полюса и ниже (севернее) южного магнитного полюса.

3.3. Синтетические кривые блеска

По результатам трехмерных расчетов нами были построены болометрические и оптические кривые блеска, которые приведены на рис. 11–15. Используемый метод синтеза этих кривых подробно описан в работе [63]. В данной работе приведем лишь основные положения метода.

Рис. 11.

Синтетические болометрические кривые блеска, построенные без учета поглощения потока в расчетной области. По оси $x$ отложена орбитальная фаза, по оси $y$ – абсолютное значение потока в единицах 1033 эрг/с. Вверху над каждой панелью указана фаза периода биений.

Рис. 12.

То же, что и рис. 11, но с учетом поглощения потока в расчетной области.

Рис. 13.

Синтетические кривые блеска видимого диапазона, построенные без учета поглощения потока в расчетной области. По оси $x$ отложена орбитальная фаза, по оси $y$ – абсолютное значение потока в единицах 1029 эрг/с. Вверху над каждой панелью указана фаза периода биений.

Рис. 14.

То же, что и рис. 13, но с учетом эффекта прогрева донора рентгеновским излучением горячих пятен.

Рис. 15.

То же, что и рис. 13, но с учетом поглощения потока в расчетной области и эффекта прогрева донора рентгеновским излучением горячих пятен.

Перед расчетом потока излучения от компонентов двойной системы необходимо поверхности донора и аккретора разбить на элементарные площадки с постоянными шагами сферической сетки по углам. При этом форма поверхности донора определяется соответствующей полостью Роша (см. выражение (2)), а аккретора – сферой радиуса ${{R}_{a}}$. Площадь каждой такой площадки рассчитывается по формулам сферической геометрии. Принимая во внимание сложную форму поверхности донора, его элементарные площадки рассматривались как проекции на сферу, помещенную внутрь звезды, с соответствующим пересчетом площади.

Интегральный поток ${{F}_{{\nu ,0}}}$ излучения от каждого компонента двойной системы в заданном диапазоне частот от ${{\nu }_{{min}}}$ до ${{\nu }_{{max}}}$ можно определить из следующего выражения:

(18)
${{F}_{{\nu ,0}}} = \int\limits_{{{\nu }_{{min}}}}^{{{\nu }_{{max}}}} \,d\nu \int\limits_S \,dS{{B}_{\nu }}(T)cos\beta ,$
где ${{B}_{\nu }}(T)$ – функция Планка, $\beta $ – угол между единичными векторами нормали к площадке на поверхности звезды $\vec {n}$ и луча зрения ${{\vec {a}}_{0}}$. Во внутреннем интеграле интегрирование проводится по всей поверхности $S$ каждой звезды. Шаг интегрирования по частоте выбирается исходя из требуемой точности расчета кривой блеска. В данной модели мы использовали логарифмическое разбиение интервала частот на 100 отрезков.

При расчете потока излучения температура площадки на поверхности каждой звезды с учетом эффекта прогрева донора рентгеновским излучением горячих пятен и наличия зон энерговыделения на аккреторе рассчитывалась по формулам (23) и (17) соответственно.

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

Расчет поглощения проводится на основе формулы для коэффициента поглощения, обусловленного свободно-свободными переходами электрона в поле протона [64]:

(19)
${{\alpha }_{\nu }} = \frac{{n_{{\text{e}}}^{2}16{{\pi }^{2}}{{e}^{6}}kT}}{{{{\nu }^{3}}ch{{{(6\pi {{m}_{{\text{e}}}}kT)}}^{{3/2}}}}},$
где ${{n}_{{\text{e}}}}$ – концентрация электронов, $e$ – элементарный заряд, ${{m}_{{\text{e}}}}$ – масса электрона, $c$ – скорость света, $h$ – постоянная Планка, $k$ – постоянная Больцмана.

Коэффициент рассеяния излучения на свободных электронах вычислялся по формуле Томсона:

(20)
$\sigma = {{\sigma }_{0}}{{n}_{e}},$
где ${{\sigma }_{0}}$ – постоянная Томсона.

Интегральное действие эффектов поглощения и рассеяния оценивается с помощью оптической толщины ${{\tau }_{\nu }}$ слоя вещества на пути прохождения луча из следующего выражения:

(21)
${{\tau }_{\nu }} = \int\limits_{{{l}_{0}}}^{{{l}_{1}}} \,({{\alpha }_{\nu }} + \sigma )dl.$
Пределы интегрирования ${{l}_{0}}$ и ${{l}_{1}}$ являются точками на луче зрения, при этом ${{l}_{0}}$ – это точка на этой прямой, отстоящая от поверхностной площадки звезды на расстоянии одного шага интегрирования $\Delta l$, а ${{l}_{1}}$ – точка прямой, соответствующая границе расчетной области. Шаг интегрирования $\Delta l$ задается исходя из требуемой точности построения кривой блеска. Приведенные в данной работе кривые построены с параметром $\Delta l = 2{{R}_{a}}$.

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

(22)
${{F}_{\nu }} = {{F}_{{\nu ,0}}}exp( - {{\tau }_{\nu }}).$
При построении кривых блеска производилось суммирование вычисленных по формуле (22) потоков излучения донора и аккретора.

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

Среди болометрических кривых блеска мы построили графики, иллюстрирующие поток излучения от системы без учета поглощения и с учетом поглощения. Чистые болометрические кривые (рис. 11) позволяют исследовать поток излучения от горячих пятен аккретора без учета светимости донора, поскольку величина потока от последнего на несколько порядков ниже, чем излучение пятен, а также дрейф этих пятен в течение периода биений. На приведенных графиках по оси $x$ отложены фазы орбитального периода, а по оси $y$ – абсолютное значение потока от поляра в эрг/с. Выбор абсолютного значения потока обоснован возможностью унифицированного сравнения всех приведенных ниже кривых блеска. Вверху каждой панели указано значение фазы периода биений.

Анализ болометрических кривых позволяет прийти к следующим выводам. На фазах 0.9 и 0.0 периода биений (см. рис. 8) поглощение излучения в струе вещества отсутствует. В этом положении активным полюсом является южный. По оценке соотношения потоков от двух горячих пятен можно заключить, что поток от северного пятна в среднем в 8–9 раз меньше потока от южного пятна. При этом температуры зон энерговыделения отличаются в 2 раза: южное пятно разогрето до 240 000 K, а северное – до 120 000 K.

На фазе 0.1 процесс переключения полюсов почти завершился, и теперь активным становится северное пятно. Поскольку оно ближе к географическому экватору аккретора, а значит и приближается на меньшее расстояние к точке ${{L}_{1}}$, то поток энергии аккрецирующего вещества в окрестности северного магнитного полюса будет несколько выше, чем для южного полюса. Начавшийся процесс роста потока излучения от северного пятна достигает своего максимума на фазе 0.3, когда его значение составляет $4 \times {{10}^{{33}}}$ эрг/с.

Если проследить изменение максимума потока от системы в целом по орбитальной фазе, то при смене магнитных полюсов происходит его смещение на величину несколько меньшую, чем 0.5 периода. Как указано в предыдущем разделе, это связано с тем, что область аккреции вещества никогда не совпадает точно с магнитным полюсом вследствие конфигурации магнитных линий. Поэтому фазовый угол между горячими пятнами может быть либо больше, либо меньше 0.5, в зависимости от структуры магнитного поля и расположения диполя. В данном случае имеет место первый вариант. Аналогичную картину можно наблюдать при обратной смене магнитных полюсов на фазах периода биений 0.5–0.6. На последней хорошо виден процесс переключения.

Как показывают кривые блеска, максимум потока излучения от южного горячего пятна достигается на фазе 0.9, он составляет $1.9 \times {{10}^{{33}}}$ эрг/с, что почти в 2 раза меньше максимального потока северного пятна. Смещение максимума к фазе 0.9 при теоретически ожидаемом на фазе 0.7–0.8 может быть связано с траекторией струи и меньшей активностью южного полюса. Баллистическая часть этой траектории имеет большую длину, чем магнитная, в окрестности указанных фаз, поэтому энергетически выгодное взаимоположение аккретора и струи, при котором аккреция максимальна, приходится как раз на фазу 0.9.

При поглощении в веществе струи (см. рис. 8) излучение аккретора ослабляется максимально в 2 раза по сравнению с чистым потоком. Такая величина ослабления свойственна фазам 0.1 и 0.5 периода биений. Соответствующая оптическая толщина слоя струи составляет 0.7. Это означает, что его можно считать оптически тонким.

На рис. 13–15 представлены кривые блеска в видимом диапазоне (фильтр $V$). При синтезе этих кривых учитывался эффект прогрева поверхности донора за счет переработки рентгеновского излучения горячих пятен аккретора. Результирующая температура донора $T_{{\text{d}}}^{'}$ вычислялась по следующей формуле [2]:

(23)
$T_{{\text{d}}}^{'} = {{\left( {T_{{\text{d}}}^{4} + \frac{{{{L}_{{\text{X}}}}\kappa }}{{4\pi {{\sigma }_{{{\text{SB}}}}}{{d}^{2}}}}} \right)}^{{1/4}}},$
где ${{T}_{{\text{d}}}}$ – эффективная температура донора, ${{L}_{{\text{X}}}}$ – рентгеновский поток от горячего пятна, $\kappa $ – коэффициент переработки донором рентгеновского излучения аккретора, $d$ – расстояние от горячего пятна до данной площадки на поверхности донора. Коэффициент $\kappa $ может принимать значения от 0 до 1, однако на практике конкретная его величина выбирается, исходя из соответствия получаемой синтетической и наблюдаемой кривых блеска. В нашей модели мы использовали значение $\kappa = 0.1$.

Анализ приведенных оптических кривых блеска выявляет следующие особенности. С учетом чернотельного характера излучения системы в видимом диапазоне изменилось соотношение потоков излучения от горячих пятен, поскольку здесь используется лишь часть функции Планка, соответствующая данному диапазону длин волн. Так, общий поток видимого излучения поляра уменьшился в среднем на 3–4 порядка по сравнению с болометрическим. Динамика прогрева донора рентгеновским излучением горячих пятен показывает максимумы на фазах 0.3 для северного пятна и 0.7 для южного пятна. При этом излучение донора в окрестности орбитальной фазы 0.5, где оно имеет наибольшую величину с учетом прогрева, остается на порядок меньше, чем излучение зон энерговыделения аккретора.

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

Рассмотрим более подробно процесс переключения струи между магнитными полюсами. В данной двойной системе он происходит дважды за период биений: между фазами 0.0–0.1 и 0.5–0.6. На рис. 16–19 представлен первый процесс переключения с южного полюса на северный. Выбор шага для последовательности фаз периода биений, на которых происходит полный цикл переключения, обусловлен значительными временными затратами на расчет каждой фазы. Мы исходили из начального разбиения всего периода биений на 10 фаз. Было определено, что начало первого процесса переключения происходит между фазами 0.0 и 0.1. Далее этот фазовый отрезок каждый раз делился пополам в поисках момента начала течения струи на северный полюс.

Рис. 16.

Первый процесс переключения магнитных полюсов – с южного на северный. Показаны фазы периода биений 0.0 (верхняя панель) и 0.025 (нижняя панель). Орбитальная фаза системы равна 0.7.

Рис. 17.

То же, что и рис. 16, но для фаз периода биений 0.050 (верхняя панель) и 0.075 (нижняя панель).

Рис. 18.

То же, что и рис. 16, но для фаз периода биений 0.1 (верхняя панель) и 0.103125 (нижняя панель).

Рис. 19.

То же, что и рис. 16, но для фаз периода биений 0.10625 (верхняя панель) и 0.1125 (нижняя панель).

В результате получились следующие фазы периода биений: 0.025, 0.050 и 0.075. Поскольку фаза 0.1 находится примерно во второй половине процесса, был продолжен поиск конца переключения за данной точкой по приведенному алгоритму. Оказалось, что первое переключение полюсов заканчивается около фазы 0.12. На фазовом отрезке 0.1–0.12 были установлены метки 0.103125, 0.10625 и 0.1125. Подчеркнем еще раз, что для каждой такой фазовой метки проводился полноценный трехмерный МГД расчет структуры течения. Таким образом, определились 6 промежуточных фаз для демонстрации процесса переключения полюсов. На рис. 16–19 также повторно показаны фазы 0.0 и 0.1, чтобы представить полную картину переключения. Далее поэтапно рассмотрим особенности первого процесса переключения.

На нулевой фазе периода биений на северный магнитный полюс аккрецирует вещество из общей оболочки системы (рис. 16, верхняя панель). Однако уже на фазе 0.025 (рис. 16, нижняя панель) данный полюс начинает захватывать часть струи, в результате чего она становится шире. При этом поток вещества из струи на северный полюс имеет меньшую плотность, примерно в 8–10 раз, по сравнению с частью потока на южный полюс. Кроме того, аккреция вещества из общей оболочки на северный полюс по-прежнему сохраняется, хотя ее темп несколько уменьшается, поскольку часть вещества, следуя вдоль магнитных линий, прижимается к точке ${{L}_{1}}$.

Дальнейшее продвижение по фазе биений приводит к росту плотности части струи, аккрецирующей на северный полюс и к образованию двойного потока из точки ${{L}_{1}}$. Так, на фазе 0.050 (рис. 17, верхняя панель) можно видеть центральную точку процесса переключения, в которой плотность потока на оба полюса имеет равную величину. Также стоит отметить слоистую структуру течения: плотность струи к ее краю уменьшается, тогда как центральная часть является более плотной. На всех последующих фазах процесса переключения наблюдается описанная выше дополнительная аккреция вещества из атмосферы донора (из общей оболочки системы). На фазе 0.075 (рис. 17, нижняя панель), когда активность аккреции на южный полюс начинает спадать, становится заметным удлинение баллистической части траектории южного компонента струи. Северный полюс, напротив, демонстрирует рост магнитной части струи, доходящей практически до точки ${{L}_{1}}$. В окрестности фазы 0.1 поток вещества на южный полюс приобретает тонкую структуру и на фазе 0.10625 происходит отрыв струи от данного полюса, при этом еще некоторое время сохраняется баллистическая часть его южного компонента, что ведет к утолщению потока на границе магнитосферы. По нашим оценкам, радиус магнитосферы аккретора составляет 7.5${{R}_{{\text{a}}}}$.

Иллюстрация второго процесса переключения полюсов приведена на рис. 20–23. Так же, как и в первом случае, длительность процесса составляет около 0.1 фазы периода биений. В данном процессе плоскость течения разделенной струи повернута на 40° относительно предыдущей стадии переключения, поэтому для лучшей его детализации была выбрана орбитальная фаза 0.3. На начальном этапе, на фазе биений 0.525 (рис. 20, верхняя панель) захват вещества донора происходит на середине расстояния между аккретором и точкой ${{L}_{1}}$. При этом на южный полюс продолжает поступать значительное количество материи из общей оболочки системы, размер арки струи несколько превышает радиус магнитосферы белого карлика. Плотность вещества южного компонента течения меньше, чем в аккрецирующем потоке на северный полюс. Однако можно заметить несколько более плотное течение из атмосферы донора, прижатое к точке Лагранжа. Этот плотный поток сохраняется на протяжении всего процесса переключения, несмотря на то, что южный полюс впоследствии захватывает часть вещества из струи. Уже на фазе 0.550 (рис. 20, нижняя панель) значение плотности южного компонента струи увеличивается и приближается к северному потоку, равно как и кольцо струи выравнивается с границей магнитосферы.

Рис. 20.

Второй процесс переключения магнитных полюсов – с северного на южный. Показаны фазы периода биений 0.525 (верхняя панель) и 0.550 (нижняя панель). Орбитальная фаза системы равна 0.3.

Рис. 21.

То же, что и рис. 20, но для фаз периода биений 0.575 (верхняя панель) и 0.6 (нижняя панель).

Рис. 22.

То же, что и рис. 20, но для фаз периода биений 0.60625 (верхняя панель) и 0.609375 (нижняя панель).

Рис. 23.

То же, что и рис. 16, но для фаз периода биений 0.6125 (верхняя панель) и 0.625 (нижняя панель).

Фазы 0.575 и 0.6 (рис. 21) иллюстрируют формирование двойного потока на оба полюса и увеличение плотности его южного компонента. Стоит отметить значительную длину баллистической части траектории струи на данных фазах. Далее в течение 0.01 фазы биений процесс протекает довольно быстро и уже на фазе 0.609 происходит резкое уменьшение северного потока и его разрыв (рис. 22). По аналогии с первым переключением последние стадии данного процесса характеризуются утолщением струи в сторону северного магнитного полюса (рис. 23) на границе магнитосферы.

4. ЗАКЛЮЧЕНИЕ

В работе исследована структура течения в асинхронном поляре CD Ind в предположении, что собственное магнитное поле белого карлика является дипольным, но его центр смещен относительно центра звезды. Такое предположение было выдвинуто в работе [15] на основе интерпретации наблюдательных данных, полученных космической обсерваторией TESS (Transiting Exoplanet Survey Satellite [65]).

Для численного моделирования мы использовали трехмерный численный МГД код, развитый нами ранее для поляров [39]. В основе численной модели лежит приближение модифицированной магнитной гидродинамики, описывающее динамику плазмы в очень сильном внешнем магнитном поле с учетом волновой альфвеновской турбулентности при малых магнитных числах Рейнольдса. Использованная в работе расчетная область полностью включает в себя как полость Роша аккретора, так и полость Роша донора. Это означает, что в рамках данной модели в отличие от наших предыдущих работ (см., напр., [4]) формирование истечения из оболочки донора в окрестности внутренней точки Лагранжа происходит естественным путем, а не за счет искусственно заданных граничных условий. Кроме того, это позволяет по результатам трехмерных расчетов синтезировать кривые блеска в различных диапазонах спектра. Для численного решения уравнений магнитной гидродинамики мы использовали разностную схему Роу–Ошера–Эйнфельдта, подробно описанную в работе [39]. Эта схема обладает низкой численной вязкостью, что позволяет получать решение с более высоким пространственным разрешением.

Мы провели десять расчетов структуры течения в CD Ind, которые соответствуют десяти фазам спин-орбитального периода биений ${{P}_{{{\text{beat}}}}}$. Для каждого варианта расчет производился при условии синхронного вращения аккретора до выхода решения на квазистационарный режим, который определялся примерным постоянством полной массы вещества в расчетной области. Как правило, установление квазистационарного течения в численном решении происходит за время порядка одного орбитального периода. Поскольку для системы CD Ind ${{P}_{{{\text{beat}}}}} = 91.3{{P}_{{{\text{orb}}}}}$, то такой подход представляется корректным.

По результатам расчетов проведен анализ структуры течения в зависимости от фазы спин-орбитального периода биений ${{P}_{{{\text{beat}}}}}$, а также изучено соответствующее распределение зон аккреции на поверхности белого карлика. Кроме того, нами были синтезированы болометрические и оптические кривые блеска с целью продемонстрировать особенности изменения потока излучения от системы в течение как орбитального, так и спин-орбитального периодов. Полученные синтетические кривые блеска в оптическом диапазоне с учетом всех необходимых факторов ($V$-фильтр, прогрев донора рентгеновским излучением из зон аккреции, поглощение излучения веществом аккреционной струи) вполне согласуются с наблюдаемыми кривыми, приведенными в работе [15]. Однако более детальное сравнение этих кривых требует отдельного исследования.

Численное моделирование показало, что зоны энерговыделения сосредоточены около магнитных полюсов. Однако в течение спин-орбитального периода биений происходит заметный дрейф горячих пятен. Анализ результатов моделирования позволяет заключить, что за время ${{P}_{{{\text{beat}}}}}$ горячие пятна смещаются по долготе относительно северного и южного магнитных полюсов в среднем на 20°, что соответствует величине 0.05 фазы орбитального периода. При этом смещения по широте практически не происходит. Кроме того, из-за особенностей геометрии магнитных силовых линий и их ориентацией относительно экваториальной плоскости двойной системы северное пятно формируется ниже (южнее), а южное пятно формируется выше (севернее) соответствующего магнитного полюса по широте.

В ходе расчетов было обнаружено, что в системе CD Ind дважды за период биений происходит переключение магнитных полюсов, в ходе которого резко изменяется конфигурация аккреционной струи. В начале этого процесса активным является один магнитный полюс, на который идет аккреция, а в его конце активным становится другой магнитный полюс. При этом из-за смещенного центра диполя относительно центра звезды, а также из-за наклона диполя относительно оси вращения этот процесс происходит несимметричным по времени образом. Первое переключение полюсов осуществляется в интервале фаз периода биений от 0.0 до 0.1, а второй между фазами 0.5 и 0.6. Для более детального исследования процесса переключения полюсов были проведены дополнительные серии расчетов в соответствующих диапазонах фазы спин-орбитального периода биений.

Темп аккреции в течение спин-орбитального периода биений испытывает незначительные вариации. Зависимость темпа аккреции от фазы спин-орбитального периода представлена сплошной линией (полученным численным значениям соответствуют кружки) на рис. 24. Показаны также соответствующие зависимости для отдельных значений темпа аккреции на северный (пунктирная линия, квадратики) и южный (штрих-пунктирная линия, треугольники) магнитные полюса. Наиболее значимые вариации полного темпа аккреции (порядка 12%) можно наблюдать в процессе переключения магнитных полюсов. В начале этого процесса темп аккреции немного уменьшается, что говорит о накоплении вещества в магнитосфере аккретора. Затем это вещество выпадает на звезды, что приводит к увеличению темпа аккреции. После этого величина темпа аккреции выходит на постоянное значение. Величины темпов аккреции на отдельные магнитные полюса изменяются со временем в противофазе. Одну половину спин-орбитального периода активным является один полюс, а вторую половину – другой. В процессе переключения полюсов величины отдельных темпов аккреции меняются местами. Один резко уменьшается примерно в 10 раз, а другой, наоборот, резко увеличивается.

Рис. 24.

Зависимость темпа аккреции от фазы спин-орбитального периода биений. Показаны полный темп аккреции (сплошная линия, кружки), темп аккреции на северный магнитный полюс (пунктирная линия, квадратики) и темп аккреции на южный магнитный полюс (штрих-пунктирная линия, треугольники). Значения темпа аккреции приведены в единицах ${{10}^{{ - 9}}}{{M}_{ \odot }}{\text{/год}}$.

Анализ результатов расчетов позволяет сделать вывод о том, что процесс переключения происходит довольно быстро, за время, не превышающее 0.1 ${{P}_{{{\text{beat}}}}}$. В начале этого процесса аккреция идет только на один магнитный полюс. Затем на границе магнитосферы начинает формироваться второй аккреционный поток. В середине этого процесса аккреция осуществляется на оба полюса с одинаковой интенсивностью, а струя выглядит в форме арки. В какой-то момент арка достигает максимального размера и в ней накапливается максимальная масса вещества. В случае смещенного диполя магнитное поле является несимметричным, поскольку индукция поля на одном полюсе больше, чем на другом. Это обстоятельство должно выражаться в некотором отличии в форме арки в первом и втором процессах переключения полюсов. Далее первый поток начинает ослабляться и, наконец, диссипирует. В конце процесса переключения полюсов остается только вторая аккреционная струя, а практически все вещество аккрецирует на второй магнитный полюс. Мы моделировали процесс переключения полюсов при условии синхронного вращения аккретора до выхода решения на квазистационарный режим, также как и остальные варианты. Однако динамическое время формирования арки оказывается порядка одного–двух орбитальных периодов. Поэтому, строго говоря, такой подход не дает полной картины. Более корректным подходом при исследовании процесса переключения полюсов является моделирование нестационарной структуры течения с учетом собственного вращения аккретора. Отметим, что сложная структура струи (например, в форме арки) в процессе переключения полюсов должна особенным образом проявляться на кривой блеска на орбитальных фазах, на которых происходит ослабление потока излучения от горячих пятен за счет поглощения аккрецирующим веществом.

Использованная нами методика построения синтетических болометрических кривых блеска и кривых блеска в видимом диапазоне по результатам трехмерного численного моделирования позволяет проводить более детальное сравнение результатов расчетов с наблюдениями. С учетом наличия ударной волны в основании аккреционной колонки можно синтезировать кривую блеска в рентгеновском диапазоне. Кроме того, важная информация о физических процессах, происходящих в полярах, содержится и в ультрафиолетовом диапазоне спектра [66]. Наш подход позволяет получать соответствующие синтетические ультрафиолетовые кривые блеска. Поэтому можно надеяться, что использование наблюдательных возможностей космической обсерватории Спектр-УФ, запуск которой запланирован на 2025 г., предоставит уникальные возможности для более детального изучения процессов аккреции как в CD Ind, так и в других полярах.

Список литературы

  1. B. Warner, Cataclysmic Variable Stars (Cambridge: Cambridge Univ. Press, 2003).

  2. А. М. Черепащук, Тесные двойные звезды, т. I, II (М.: Физматлит, 2013).

  3. S. Tapia, Bull. Amer. Astron. Soc. 8, 511 (1976).

  4. А. Г. Жилкин, Д. В. Бисикало, А. А. Боярчук, Успехи физ. наук 182, 121 (2012).

  5. Д. В. Бисикало, А. Г. Жилкин, А. А. Боярчук, Газодинамика тесных двойных звезд (М.: Физматлит, 2013).

  6. C. G. Campbell, Magnetohydrodynamics in binary stars (Dordrecht: Kluwer Acad. Publishers, 1997).

  7. C. G. Campbell and A. D. Schwope, Astron. and Astrophys. 343, 132 (1999).

  8. J. Patterson, Publ. Astron. Soc. Pacific 106, 209 (1994).

  9. А. Г. Жилкин, Д. В. Бисикало, П. А. Масон, Астрон. журн. 89(4), 291 (2012).

  10. S. Vennes, D. T. Wickramasinghe, J. R. Thorstensen, D. J. Christian, and M. J. Bessell, Astron. J. 112, 2254 (1996).

  11. A. Schwope, D. H. Buckley, D. O’Donoghue, G. Hasinger, J. Trümper, and W. Voges, Astron. and Astrophys. 326, 195 (1997).

  12. G. Ramsay, D. H. Buckley, M. Cropper, and M. K. Harrop-Allin, Monthly Not. Roy. Astron. Soc. 303, 96 (1999).

  13. G. Ramsay, D. H. Buckley, M. Cropper, M. K. Harrop-Allin, and S. Potter, Monthly Not. Roy. Astron. Soc. 316, 225 (2000).

  14. G. Myers, J. Patterson, E. de Miguel, F.-J. Hambsch, et al., Publ. Astron. Soc. Pacific 129, 4204 (2017).

  15. P. Hakala, G. Ramsay, S. B. Potter, A. Beardmore, D. H. Buckley, and G. Wynn, Monthly Not. Roy. Astron. Soc. 486, 2549 (2019).

  16. A. R. King, Monthly Not. Roy. Astron. Soc. 261, 144 (1993).

  17. G. A. Wynn and A. R. King, Monthly Not. Roy. Astron. Soc. 275, 9 (1995).

  18. G. A. Wynn, A. R. King, and K. Horne, Monthly Not. Roy. Astron. Soc. 286, 436 (1997).

  19. A. R. King and G. A. Wynn, Monthly Not. Roy. Astron. Soc. 310, 203 (1999).

  20. A. J. Norton, J. A. Wynn, and R. V. Somerscales, Astrophys. J. 614, 349 (2004).

  21. N. R. Ikhsanov, V. V. Neustroev, and N. G. Beskrovnaya, Astron. and Astrophys. 421, 1131 (2004).

  22. A. J. Norton, O. W. Butters, T. L. Parker, and G. A. Wynn, Astrophys. J. 672, 524 (2008).

  23. Yu. V. Glagolevskij, Astrophys. Bull. 66(2), 144, (2011).

  24. Е. С. Беленькая, Успехи физ. наук 179(8), 809 (2009).

  25. K. M. Moore, R. K. Yadav, L. Kulowski, H. Cao, et al., Nature, 561, 76 (2018).

  26. А. Г. Жилкин, Матем. моделирование 22(1), 110 (2010).

  27. А. Г. Жилкин, Д. В. Бисикало, Астрон. журн. 87(12), 1155 (2010).

  28. Е. П. Курбатов, А. Г. Жилкин, Д. В. Бисикало, Успехи физ. наук 187, 857 (2017).

  29. А. Г. Жилкин, Д. В. Бисикало, Астрон. журн. 86(5), 475 (2009).

  30. A. G. Zhilkin and D. V. Bisikalo, Adv. Space Research 45, 437 (2010).

  31. А. Г. Жилкин, Д. В. Бисикало, Астрон. журн. 87(9), 913 (2010).

  32. Д. В. Бисикало, А. Г. Жилкин, П. В. Кайгородов, В. А. Устюгов, М. М. Монтгомери, Астрон. журн. 90, 366 (2013).

  33. В. А. Устюгов, А. Г. Жилкин, Д. В. Бисикало, Астрон. журн. 90(11), 885 (2013).

  34. А. М. Фатеева, А. Г. Жилкин, Д. В. Бисикало, Астрон. журн. 92, 977 (2015).

  35. П. Б. Исакова, А. Г. Жилкин, Д. В. Бисикало, Астрон. журн. 92(9), 720 (2015).

  36. П. Б. Исакова, Н. Р. Ихсанов, А. Г. Жилкин, Д. В. Бисикало, Н. Г. Бескровная, Астрон. журн. 93, 474 (2016).

  37. П. Б. Исакова, А. Г. Жилкин, Д. В. Бисикало, А. Н. Семена, М. Г. Ревнивцев, Астрон. журн. 94, 566 (2017).

  38. Е. П. Курбатов, А. Г. Жилкин, Д. В. Бисикало, Астрон. журн. 96, 27 (2019).

  39. А. Г. Жилкин, А. В. Соболев, Д. В. Бисикало, М. М. Габдеев, Астрон. журн. 96(9), 748 (2019).

  40. T. Tanaka, J. Comp. Phys. 111, 381 (1994).

  41. K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, and D. L. de Zeeuw, J. Comp. Phys. 154, 284 (1999).

  42. А. Г. Куликовский, Н. В. Погорелов, А. Ю. Семенов, Математические вопросы численного решения гиперболических систем уравнений (М.: Физматлит, 2001).

  43. A. V. Koldoba, M. M. Romanova, G. V. Ustyugova, and R. V. E. Lovelace, Astrophys. J. 576, L53 (2002).

  44. M. M. Romanova, G. V. Ustyugova, A. V. Koldoba, J. V. Wick, and R. V. E. Lovelace, Astrophys. J. 595, 1009 (2003).

  45. M. M. Romanova, G. V. Ustyugova, A. V. Koldoba, J. V. Wick, and R. V. E. Lovelace, Astrophys. J. 610, 920 (2004).

  46. M. M. Romanova, G. V. Ustyugova, A. V. Koldoba, J. V. Wick, and R. V. E. Lovelace, Astrophys. J. 616, L151 (2004).

  47. D. P. Cox and E. Daltabuit, Astrophys. J. 167, 113 (1971).

  48. A. Dalgarno and R. A. McCray, Ann. Rev. Astron. Astrophys. 10, 375 (1972).

  49. J. C. Raymond, D. P. Cox, and B. W. Smith, Astrophys. J. 204, 290 (1976).

  50. Л. Спитцер, Физические процессы в межзвездной среде (М.: Мир, 1981).

  51. С. И. Брагинский, ЖЭТФ 37, 1417 (1960).

  52. S. Galtier, S. V. Nazarenko, A. C. Newell, and A. Pouquet, J. Plasma Phys. 63, 447 (2000).

  53. В. Е. Захаров, Журн. приклад. мех. тех. физ. 1, 14 (1965).

  54. П. С. Ирошников, Астрон. журн. 40, 742 (1963).

  55. R. H. Kraichnan, Phys. Fluids 8, 575 (1965).

  56. S. D. Drell, H. M. Foley, and M. A. Ruderman, J. Geophys. Res. 70, 3131 (1965).

  57. А. В. Гуревич, А. Л. Крылов, Е. Н. Федоров, ЖЭТФ 75, 2132 (1978).

  58. Р. Р. Рафиков, А. В. Гуревич, К. П. Зыбин, ЖЭТФ 115, 542 (1999).

  59. Д. А. Франк-Каменецкий, Лекции по физике плазмы (М.: Атомиздат, 1968).

  60. Ф. Чен, Введение в физику плазмы (М.: Мир, 1987).

  61. A. G. Zhilkin, D. V. Bisikalo, and V. A. Ustyugov, AIP Conf. Proc. 1551, 22 (2013).

  62. S. H. Lubow and F. H. Shu, Astrophys. J. 198, 383 (1975).

  63. А. В. Соболев, А. Г. Жилкин, Д. В. Бисикало, Научные труды ИНАСАН 3, 231 (М.: Янус-К, 2019).

  64. В. В. Соболев, Курс теоретической астрофизики (М.: Наука, 1985).

  65. G. R. Ricker, J. N. Winn, R. Vanderspek, D. W. Latham, et al., J. Astron. Tel. Instr. and Systems 1, id. 014003 (2015).

  66. А. А. Боярчук, Б. М. Шустов, И. С. Саванов, М. Е. Сачков, и др., Астрон. журн. 93(1), 3 (2016).

Дополнительные материалы отсутствуют.