Gravitational field of a homogeneous cube. Classical and relativistic case

Full Text

Введение

Вопрос исследования гравитационного поля тел сложной формы представляет большой интерес для геофизики, астрофизики, математической физики и других областей. В настоящей работе рассматривается гравитационное поле однородного куба для классического и релятивистского случаев. Форма куба выбрана не только из-за наличия нетривиальной симметрии. Ранее в литературе [1, 2] не встречались решения уравнений Эйнштейна для тел, имеющих форму однородного куба, в силу чего задача актуальна в рамках математической физики.

Кроме фундаментального значения, эта задача интересна своим практическим приложением. Несмотря на то, что астрономические тела кубической формы неизвестны, некоторые объекты могут иметь форму, обладающую схожими параметрами. Так, астероид Рюгу имеет форму, близкую к кубической, хотя ее можно описать как куполообразную [3]. Также в [4] рассматривалась модель столкновений частиц пыли с астероидами, приводящих к образованию больших относительно плоских областей и острых краев, что может свидетельствовать в пользу возможности образования астероидов кубической или подобной сложной формы.

В ряде случаев четкое представление о гравитационном поле объекта имеет важное значение для космических аппаратов, имеющих кубическую форму (например, при планировании высадки на небольшие астероиды и при выполнении других задач) [5]. Справедлива и актуальна также обратная задача по расчету орбит движения тел, находящихся в гравитационном поле тел с формой куба [6–8].

Также важно рассмотреть эту задачу с точки зрения ее значения в области геофизики. Для оценки гравитационного поля массивных объектов (возвышенностей, гор и горных хребтов) используется метод расчета их гравитационных полей, аналогичный методу конечных элементов, где в качестве единичного элемента выступает однородный куб или квадрат [9]. Аналогичным образом производится оценка потенциала гравитационного поля астероидов и галактик [8, 10].

Стоит отметить, что рассматриваемая задача имеет и научно-методологическое значение и была предложена для описания особенностей физики гравитационных полей для тел разной формы [11, 12].

К настоящему времени опубликовано достаточно много работ, посвященных расчетам гравитационного поля однородного куба в рамках классической ньютоновской механики, однако должной систематизации результатов этих работ и целостного представления обо всех методах решения до сих пор нет. Поэтому авторы посчитали важным привести литературный обзор подобных решений. Кроме того, эта тема поднимается в современных работах [13, 14]. Также сведения, приведенные в литературном обзоре, будут использоваться для рассмотрения релятивистского случая. Основное приложение этой задачи в рамках релятивистской физики относится к области математической физики (или, шире, математики).

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

1. Классический случай

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

1. получение аналитического решения;
2. как частный случай задачи нахождения гравитационного поля полиэдра;
3. применение метода конечных элементов;
4. использование мультипольного разложения.

Рассмотрим каждый способ решения.

1.1. Получение аналитического решения

Вероятно, впервые данная задача была сформулирована в [15], однако в [16, p. 522] решалась близкая по постановке задача. В работе [15] автор рассматривает формулу для вычисления проекции гравитационного потенциала параллелепипеда со сторонами $a$, $b$, $c$ вдоль оси $0z$, однако общую формулу для гравитационного потенциала куба (параллелепипеда) он не приводит, ссылаясь на излишне трудоемкие расчеты. Позже были получены общие выражения для гравитационного потенциала призмы по трем проекциям на оси [17, pp. 93–104].

В работе [18] автор исходя из теоремы Эйлера об однородных функциях доказывает, что первые производные потенциала и сам потенциал гравитационного поля можно легко вычислить из вторых производных, не интегрируя их. Данный результат позволяет вычислить гравитационный потенциал прямоугольной призмы. Подобный подход к решению данной задачи встречается в работе [19].

Практическое использование полученных результатов описано в работах [20, 21], в которых решались задачи вычисления гравитационного потенциала прямоугольной призмы в связи с моделированием гравитационного поля Земли и массивных горных образований. Потенциал (напряженность) гравитационного поля куба впервые обсуждался в работе [22], в том числе для анализа потенциала (напряженности) гравитационного поля массивных объектов сложной формы.

К основным недостаткам этих работ [15–23] можно отнести:

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

Этих недостатков лишена работа [24], где дается формула расчета потенциала, записанная с помощью сумм:
\[ \begin{multline}
U(x_0,y_0,z_0)=\sum\limits_{i=0}^1 \sum\limits_{j=0}^1
\sum\limits_{k=0}^1
\Bigl[ x_i y_j \cdot \operatorname{th}^{-1}
\Bigl( \frac{z_k}{r_{ijk}} \Bigr) + {} \\
{} +y_j z_k \cdot \operatorname{th}^{-1} \Bigl( \frac{x_i}{r_{ijk}} \Bigr)+ z_k x_i \cdot \operatorname{th}^{-1} \Bigl( \frac{y_j}{r_{ijk}} \Bigr) - \frac{x_i^2}{2} \operatorname{tg}^{-1} \Bigl( \frac{y_j z_k}{x_i r_{ijk}} \Bigr)- {}
\\
{} - \frac{y_j^2}{2} \operatorname{tg}^{-1} \Bigl( \frac{z_k x_i}{y_j r_{ijk}} \Bigr)
- z_k^2 \operatorname{tg}^{-1} \Bigl( \frac{x_i y_j}{z_k r_{ijk}} \Bigr)\Bigr],
\end{multline} \tag{1} \]
где $x_1=a-x_0$, $y_1=b-y_0$, $z_1=c-y_0$; $a$, $b$, $c$ — ребра прямоугольной призмы; $x_0$, $y_0$, $z_0$ — начало координат; $r_{ijk}=\sqrt{x_i^{2}+y_j^{2}+z_k^{2}}$, $i,j,k=0,1;$ здесь константы определены как $G=\rho=1$.

Однако формула (1), представленная в работе [24], имеет неточность (или опечатку), а именно вместо $z_k^{2}$ необходимо записать ${z_k^{2}}/{2}$. Таким образом, корректная формула выглядит так:
\[ \begin{multline}
U(x_0,y_0,z_0)=\sum\limits_{i=0}^1 \sum\limits_{j=0}^1
\sum\limits_{k=0}^1 \Bigl[ x_i y_j \cdot \operatorname{th}^{-1}
\Bigl( \frac{z_k}{r_{ijk}} \Bigr) + {}
\\
{} +y_j z_k \cdot \operatorname{th}^{-1} \Bigl( \frac{x_i}{r_{ijk}} \Bigr)+ z_k x_i \cdot \operatorname{th}^{-1} \Bigl( \frac{y_j}{r_{ijk}} \Bigr) -
\frac{x_i^2}{2} \operatorname{tg}^{-1} \Bigl( \frac{y_j z_k}{x_i r_{ijk}} \Bigr)- {}
\\
{} - \frac{y_j^2}{2} \operatorname{tg}^{-1} \Bigl( \frac{z_k x_i}{y_j r_{ijk}} \Bigr)
- \frac{z_k^{2}}{2} \operatorname{tg}^{-1} \Bigl( \frac{x_i y_j}{z_k r_{ijk}} \Bigr)\Bigr].
\end{multline} \tag{2} \]

Кроме того, уравнение (1) несимметрично, хотя от куба ожидается симметричность. В дальнейшем подобные задачи решались в [13, 14], но они не лишены проблем, которые были описаны ранее.

1.2. Как частный случай задачи нахождения гравитационного поля полиэдра

Так как космические тела могут иметь сложную форму многогранника, нахождение потенциала (напряженности) гравитационного поля куба, в том числе, можно рассматривать как частный случай нахождения потенциала (напряженности) полиэдра. В настоящей работе подробный литературный обзор становления данного метода не проводится, так как он хорошо описан в [25]. Также в работе [25] представлен оригинальный подход к решению этой задачи с помощью элементарных теорем векторного исчисления. Значение потенциала гравитационного поля однородного полиэдра плотностью $\rho$ (в соответствии с работой [25, 26]) найдем как
\[ \begin{equation}
U(\vec{r})=\frac{\rho}{2} \sum_{i} \vec{n}_i \cdot \bigl( \vec{r}-\vec{r}\,' \bigr)
\oint_{C_i} \vec{Q_i} \bigl( \vec{r}-\vec{r}\,' \bigr)\,d\vec{r}\,',
\end{equation} \tag{3} \]
где $\vec{Q}(\vec{r})$ — дополнительный векторный потенциал, причем
\[ \begin{equation*}
\iint_{S} U(\vec{r})\,dS = \oint_{C} \vec{Q}(\vec{r}) \,d\vec{r},
\end{equation*} \]
а $\vec{n}_i$ — единичные нормали.

Суть метода состоит в оценке величины поверхностного интеграла однородного многогранника:
\[ \begin{equation}
I_r(\vec{r}) = \iint_{S_i} \frac{dS'}{|\vec{r}-\vec{r}\,'|},
\end{equation} \tag{4} \]
где $I_r(\vec{r})$ — интегральное выражение по $i$-той грани полиэдра; $S_i$ — любая из плоских поверхностей полиэдра.

Запишем связь уравнений (4) и (3):
\[ \begin{equation*}
\oint_{C_i} \vec{Q_i} \bigl( \vec{r}-\vec{r}\,' \bigr)\,d\vec{r}\,' = \sum_{j} I_{i,j}.
\end{equation*} \]

Основной трудностью использования данного метода являются вычисления интегрального выражения (4). Несмотря на то, что в результате вычисления разными методами получаются отличные друг от друга значения потенциала однородного куба, все они отличаются в 7–8 знаке после десятичной запятой [25].

В общем виде интегральное выражение имеет вид
\[ \begin{multline}
I_{i,j}=d_{i,j} \left[ \frac{c_{i,j}}{\sqrt{a_{i,j}^2-b_{i,j}^2-c_{i,j}^2}} \operatorname{tg}^{-1}
\biggl( \frac{c_{i,j} ( 1-b_{i,j} )}{\sqrt{a_{i,j}^2-b_{i,j}^2-c_{i,j}^2} \sqrt{1+a_{i,j}^2-2b_{i,j}^2}} \biggr) \right. + {}
\\
{} + \frac{c_{i,j}}{\sqrt{a_{i,j}^2-b_{i,j}^2-c_{i,j}^2}} \operatorname{tg}^{-1} \biggl( \frac{c_{i,j}b_{i,j}}{a_{i,j}\sqrt{a_{i,j}^2-b_{i,j}^2-c_{i,j}^2}} \biggr) +{}
\\
{}\left.+ \ln\biggl( \frac{1-b_{i,j}+\sqrt{1+a_{i,j}^2-2b_{i,j}^2}}{a_{i,j}-b_{i,j}} \biggr)\right],
\end{multline} \tag{5} \]
где
\[ \begin{equation*}
a_{i,j}=\frac{|\vec{r}-\vec{r}_{i,j}|}{|\vec{r}_{i,j+1}-\vec{r}_{i,j}|},
\quad
b_{i,j}=\frac{( \vec{r}-\vec{r}_{i,j} )\cdot ( \vec{r}_{i,j+1}-\vec{r}_{i,j} ) }{|\vec{r}_{i,j+1}-\vec{r}_{i,j}|^2},
\end{equation*} \]
\[ \begin{equation*}
c_{i,j}=\frac{\vec{n}_i \cdot ( \vec{r}-\vec{r}_{i,j})}{|\vec{r}_{i,j+1}-\vec{r}_{i,j}|},
\quad
d_{i,j}=\frac{ ( \vec{n}_i \times ( \vec{r}-\vec{r}_{i,j} ) ) \cdot ( \vec{r}_{i,j+1}-\vec{r}_{i,j})} {|\vec{r}_{i,j+1}-\vec{r}_{i,j}|}.
\end{equation*} \]

В уравнении (5) аргументы различных членов квадратного корня и натурального логарифма могут быть нулевыми только тогда, когда точка поля $r$ лежит на краевой линии или продолжении краевой линии полиэдра. Можно заметить явное сходство формул (2) и (5). Несмотря на высокую ресурсоемкость с вычислительной точки зрения, этот метод расчета потенциала гравитационного поля однородного куба используется в работах [6, 7].

1.3. Метод конечных элементов

Задачи расчета гравитационного потенциала системы $N$ тел анализируется на основе наложения на группу точечных тел квадратной сетки, например, как в [27]. Таким образом, условно такие задачи можно отнести к расчету гравитационного потенциала (напряженности) квадрата.

1.4. Мультипольное разложение

В силу своей простоты метод расчета потенциала (напряженности) гравитационного поля мультипольным разложением (рядом Лапласа) в последнее время стал востребован для исследования орбит, создаваемых вокруг тел кубической формы [8]. Впервые задача на вычисление потенциала за пределами куба до четвертого порядка мультипольного разложения была решена в [28], а позже обсуждалась в работе [21]:
\[ \begin{equation}
U_{ext}(x,y,z) = \frac{M}{r}-\frac{7Ma^4}{18r^9} ( x^4 + y^4 + z^4 + 3 x^2 y^2 + 3 y^2 z^2 + 3 z^2 x^2 ),
\end{equation} \tag{6} \]
где $U_{ext}(x,y,z)$ — потенциал гравитационного поля вне куба, $M$ — масса куба, $a$ — половина стороны ребра куба.

Эквипотенциальные поверхности, полученные с помощью формулы (6), изображены на рис. 1.

Рис. 1. Визуализация эквипотенциальных поверхностей, полученная с помощью формулы (6)
[Figure 1. Visualization of equipotential surfaces obtained using the formula (6)]

Рис. 2. Визуализация эквипотенциальных поверхностей, полученная с помощью формулы (7)
[Figure 2. Visualization of equipotential surfaces obtained using the formula (7)]

Рис. 3. Потенциалы гравитационного поля однородного куба от координаты $x$, полученные с помощью формул (6), (7) и (2), причем $G=c=\rho=1$; вертикальные линии представляют собой границы куба: $a=0.5$
[Figure 3. Potentials of the gravitational field of a homogeneous cube as functions of the coordinate $x$ obtained using the formulas (6), (7) и (2), with $G=c=\rho=1$; the vertical lines represent the boundaries of the cube: $a=0.5$]

Формула более высоких порядков мультипольного разложения для области внутри куба приводится в работе [14]:
\[ \begin{multline}
U_{int}(x,y,z) = \Bigl( 24 \ln{\frac{1+\sqrt{3}}{\sqrt{2}}} - 2 \pi \Bigr) -
\Bigl[ \frac{2}{3} \pi ( x^2 + y^2 + z^2 ) \Bigr] + {}
\\
{} +
\Bigl[ -\frac{4}{9 \sqrt{3}} ( x^4 + y^4 + z^4 ) + \frac{4}{3 \sqrt{3}} ( x^2 y^2 + x^2 z^2 + y^2 z^2 )\Bigr] + {}
\\
{}+ \Bigl[ - \frac{1}{162 \sqrt{3}} ( x^6 + y^6 + z^6 ) +
\frac{5}{108 \sqrt{3}} ( x^2 y^4 + x^2 z^4 + x^4 y^2 + x^4 z^2 + y^2 z^4 + y^4 z^2 ) -{} \\
{}- \frac{5}{9 \sqrt{3}} x^2 y^2 z^2 \Bigr],
\end{multline} \tag{7} \]
где $U_{int}(x,y,z)$ — потенциал гравитационного поля внутри куба.

Эквипотенциальные поверхности, полученные с помощью формулы (7), изображены на рис. 2.

Рассмотрим потенциал гравитационного поля однородного куба, полученный с помощью формул (6), (7) и (2) и изображенный на рис. 3.

По рис. 3 видно, что для анализа гравитационного поля за пределами куба в релятивистском случае вполне обосновано использование формулы (6) из-за ее простоты по сравнению с формулой (2) и способности обеспечить плоское пространство на $r\gg 2a$. Аналогичное утверждение справедливо и для формулы (7) для области внутри куба.

2. Релятивистский случай

2.1. Общая постановка задачи о гравитационном поле однородного куба

Современная теория тяготения в своей основе использует систему уравнений Эйнштейна, по которым и определяется гравитационное поле тела:
\[ \begin{equation}
\begin{split}
R^{\alpha\beta}-\dfrac{1}{2}R\,g^{\alpha\beta}=\dfrac{8\pi G}{c^{4}}T^{\alpha\beta}.
\end{split}
\end{equation} \tag{8} \]
Однако поскольку уравнения тяготения являются системой нелинейных дифференциальных уравнений в частных производных, решение такой системы в общем виде является очень тяжелой задачей как с математической, так и физической точки зрения. Поэтому в настоящее время получен лишь небольшой класс точных решений [1, 2]. В данной задаче при описании гравитационного поля куба в рамках ОТО следует решать именно уравнения (8), и авторы предполагают рассмотреть этот вопрос в будущих публикациях, но ввиду его чрезвычайной сложности рассмотрим другой метод решения.

В [29] описан способ приближенных решений уравнений тяготения. В современной литературе данный способ получил название постньютоновского формализма, или PNN [30, 31]. Применим данный подход в нашем исследовании и, как следствие, получим приближенные выражения для компонент метрического тензора однородного куба. Несмотря на общеизвестность данного метода, ранее к однородному кубу он не применялся. Постньютоновский формализм для решения этой задачи выбран по причине наличия известного нерелятивистского решения и эффективного использования в практических задачах. Также постньютоновский формализм хорошо подходит для расчета гравитационных полей малой напряженности, в особенности в применении к телам несферической формы, например, с формой куба.

В рамках рассматриваемой задачи определим куб следующим образом. Пусть в плоском пустом пространстве, которое задается обычной галилеевой метрикой, находится куб со стороной $2a$ и метрикой
\[ \begin{equation}
ds_{0}^{2}=c^{2}dt^{2}-dx^{2}-dy^{2}-dz^{2},
\end{equation} \tag{9} \]
центр куба совмещен с началом координат; $V=8a^{3}$ — объем куба.

В данной метрике все компоненты тензора энергии-импульса куба равны нулю, никакого воздействия на геометрию пространства он не оказывает. Когда куб заполнен веществом, данная метрика должна быть справедлива на бесконечности, поэтому на больших расстояниях компоненты метрического тензора следующие:
\[ \begin{equation*}
g_{00 \infty} =c^{2}, \quad g_{ik \infty}= -\delta_{ik}.
\end{equation*} \]
Теория Эйнштейна, основанная на уравнениях (8), в первом приближении должна давать то же, что и ньютоновская теория тяготения. Поскольку теория Ньютона нерелятивистская, а перед нами стоит задача решать релятивистские уравнения, необходимо выделить некий параметр, характеризующий переход. Таким параметром формально будет являться скорость света $c$.

Для решений уравнений тяготения с достаточной точностью мы будем использовать приближенный метод [29]. Этот метод основан на разложении всех искомых функций по обратным степеням скорости света $U/c^{2}$. Наша задача облегчается тем, что метрика везде мало отличается от евклидовой. Таким образом, в результате решения задачи мы получим приближенные выражения для метрического тензора.

2.2. Физическая модель

Пусть координатный куб в пространстве (9) заполнен несжимаемой жидкостью с нулевой скоростью и плотностью $\rho$. Равновесие обеспечивается изотропным давлением жидкости, которое можно найти как
\[ \begin{equation*}
\rho \frac{\partial U}{\partial x_{i}}=-\frac{\partial p_{ik}}{\partial x_{i}},
\end{equation*} \]
где $x_{i}=x, y, z$, а $U$ — гравитационный потенциал куба.

Когда $\rho=\mathrm{const}$, получим, что
\[ \begin{equation*}
p=\rho U.
\end{equation*} \]

Для тензора энергии-импульса в правой части (8) мы используем приближенные выражения, которые соответствуют евклидовой метрике:
\[ \begin{equation}
T^{00}= {\rho}/{c^{2}}, \quad
T^{0\alpha}= 0, \quad
T^{ik}= {p} / {c^{2}}.
\end{equation} \tag{10} \]

При решении внешних уравнений будем считать, что тензор энергии-импульса во всем пространстве равен нулю, кроме некоторого отдельного объема $8a^{3}$, поэтому масса куба
\[ \begin{equation}
M_{cube} = \int_{-a}^{a}\int_{-a}^{a}\int_{-a}^{a} T dxdydz=8\rho a^{3},
\end{equation} \tag{11} \]
где $T=T^{00}g_{00}$.

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

2.3. Первое приближение

После заполнения координатного куба веществом по заданному распределению масс можно найти распределение гравитационного потенциала. В первом приближении единственной компонентой, вносящей вклад, будет компонента
\[ \begin{equation*}
T^{00}= {\rho}/{c^{2}}.
\end{equation*} \]

Перепишем уравнение (8) в виде
\[ \begin{equation*}
R^{\mu\nu}=- \varkappa\Bigl(T^{\mu\nu}-\dfrac{1}{2}g^{\mu\nu}T\Bigr),
\end{equation*} \]
оставляя в первом приближении компоненты тензора Риччи:
\[ \begin{equation}
R^{00}=\dfrac{1}{2}\Delta g^{00},
\quad
R^{ik}=\dfrac{1}{2}\Delta g^{ik}\delta^{ik}.
\end{equation} \tag{12} \]
Таким образом, метрика внутри куба удовлетворяет уравнениям
\[ \begin{equation}
\Delta g^{00}=-\dfrac{8\pi G}{c^{4}}\rho, \quad
\Delta g^{ik}=-\dfrac{8\pi G}{c^{2}}\rho \delta_{ik}.
\end{equation} \tag{13} \]

Решение уравнений тяготения (13) определит приближенные выражения для компонент метрического тензора в виде
\[ \begin{equation}
g^{00}=\dfrac{1}{c^{2}} \Bigl(1+\dfrac{2U}{c^{2}}\Bigr), \quad
g^{ik}=- \left(1-\dfrac{2U}{c^{2}}\right)\delta_{ik} .
\end{equation} \tag{14} \]
Переход к ковариантным компонентам дает выражения
\[ \begin{equation*}
g_{00}=c^{2}-2U , \quad
g_{ik}=- \Bigl(1+\dfrac{2U}{c^{2}}\Bigr)\delta_{ik} .
\end{equation*} \]

Данные выражения для компонент метрического тензора удовлетворяют в первом приближении уравнениям тяготения. Соответствующее им выражение для квадрата интервала имеет вид
\[ \begin{equation*}
ds^{2}_{I}=(c^{2}-2U)dt^{2}-\Bigl(1+\dfrac{2U}{c^{2}}\Bigr) (dx^{2}+dy^{2}+dz^{2} )
\end{equation*} \]
с соответствующим ньютоновым потенциалом (7) и удовлетворяет выражению
\[ \begin{equation*}
\Delta U_{int} =-4\pi G\rho.
\end{equation*} \]

Приближение накладывает ограничение на потенциал в виде
\[ \begin{equation*}
U_{\max}\ll {c^{2}} /{2},
\end{equation*} \]
что в свою очередь ограничивает плотность вещества куба
\[ \begin{equation*}
\rho<\frac{c^{2}}{a^{2}G}.
\end{equation*} \]

2.4. Второе приближение

Во втором приближении необходимо еще раз решить уравнения тяготения, но на этот раз учесть следующий порядок точности по $1/c$. Поэтому точность, с которой необходимо вычислить временную компоненту метрического тензора, равна
\[ \begin{equation*}
g^{00}=\frac{1}{c^{2}}+\frac{f(x,y,z)}{c^{4}}+\frac{f(x,y,z)}{c^{6}},
\end{equation*} \]
а точность, с которой будут вычислены пространственные компоненты, —
\[ \begin{equation*}
g^{ik}=-\delta^{ik}+\frac{g(x,y,z)}{c^{4}}.
\end{equation*} \]
Правая часть уравнений тяготения останется без изменений, однако учтем пространственные компоненты. Тогда во втором приближении имеем
\[ \begin{equation}
T^{00}= {\rho} / {c^{2}} , \quad
T^{ik}= {\rho U} /{c^{2}}.
\end{equation} \tag{15} \]
Контравариантный тензор Риччи представим в виде
\[ \begin{equation}
R^{\mu\nu}=-\dfrac{1}{2} g^{\alpha\beta}\frac{\partial^{2} g^{\mu\nu}}{\partial x_{\alpha}\partial x_{\beta}}+\Gamma^{\mu\,\alpha\beta}\Gamma^{\nu}_{\alpha\beta}-\Gamma^{\mu\nu}.
\end{equation} \tag{16} \]
Первый член (16) в ньютоновом приближении дает (12); последний член в (16) является сложной квадратичной функцией, содержащей первые производные от $g^{\mu\nu}$. Поэтому во втором приближении должно учитываться попарное произведение скобок Кристоффеля с той точностью, которого требует приближение. Во втором приближении мы ограничимся членами порядка
\[ \begin{equation*}
\Gamma^{0}_{\mu \nu}\sim 1/c^{2}, \quad \Gamma^{0}_{ik}\sim O(1/c^{4}),\quad
\Gamma^{i}_{k l}\sim 1/c^{2}.
\end{equation*} \]

При решении уравнений тяготения мы свободны в выборе системы координат. Наиболее простой вид приближенные уравнения тяготения имеют в гармонических координатах, таких, что
\[ \begin{equation}
\frac{\partial}{\partial x_{\mu}} (\sqrt{-g}\:g^{\mu \nu} )=0, \quad \Gamma^{\mu\nu}=0.
\end{equation} \tag{17} \]

Для удобства дальнейшего представления введем величину
\[ \begin{equation*}
\mathcal{G}^{\mu \nu} = \sqrt{-g}g^{\mu \nu}, \quad \sqrt{-g}\approx c.
\end{equation*} \]
Компоненты $\mathcal{G}^{\mu \nu}$ в первом приближении равны $ds_{I}$:
\[ \begin{equation*}
\mathcal{G}^{00} =\frac{1}{c}+\dfrac{4}{c^{3}}U, \quad
\mathcal{G}^{ik} =- \delta^{ik} c.
\end{equation*} \]
Компоненты $\mathcal{G}^{\mu \nu}$ во втором приближении равны $ds_{II}$:
\[ \begin{equation}
\begin{split}
\mathcal{G}^{00} &=\frac{1}{c}+\dfrac{4}{c^{3}}V(x,y,z)+\dfrac{7}{c^{5}}V(x,y,z)^{2}, \\
\mathcal{G}^{ik} &=- \delta^{ik} c+\dfrac{4}{c^{3}}S_{ik}(x,y,z),
\end{split}
\end{equation} \tag{18} \]
где $V$ и $S$ — некоторые функции разложения, которые удовлетворяют приближенным уравнениям Эйнштейна (в соответствии с обозначениями, принятыми в [29]).

Используя условие (17), подставим разложение (18) в тензор Риччи (16). Приравнивая полученное выражение компонентам тензора энергии-импульса (15), получим два уравнения для релятивистских поправок. Поправка для координаты $t$:
\[ \begin{equation}
\Delta V=4\pi G \rho\Bigl(c^{2}+\dfrac{1}{2}U\Bigr).
\end{equation} \tag{19} \]
Поправка пространственной координаты, которая позволит вычислить давление $p$:
\[ \begin{equation}
\Delta S=-4\pi G\rho U-\frac{1}{2}\Bigl(\frac{\partial U}{\partial x_{i}}\Bigr)^{2}.
\end{equation} \tag{20} \]

Уравнения (20) и (19) для релятивистских поправок представляют собой уравнения Пуассона, которые решаются, например, с использованием функции Грина [32, 33]. Но так как физическая модель представляет собой координатный куб в пространстве (9), решать уравнения (20) и (19) стоит именно в декартовых координатах $(x,y,z)$.

Получение аналитических решений для уравнений (20) и (19) для области внутри куба (т.е. области на некотором расстоянии $R<a$) с использованием внутреннего потенциала куба (7) представляет собой крайне непростую задачу. Поэтому решение для уравнений (20) и (19) для области внутри куба получено с помощью численных методов.

Граничные условия при решении уравнений (20) и (19): $S(a)=0$; ${V(a)=0}$. Результат представлен на рис. 4 и 5.

Рис. 4. Результат расчета релятивистских поправок $V(x,y)$ и $V(x,y,z)$ в системе Wolfram Mathematica 13.0 с использованием численных методов ($c=G=\rho=1$)
[Figure 4. The result of the calculation of relativistic corrections $V(x,y)$ and $V(x,y,z)$ in the Wolfram Mathematica 13.0 system using numerical methods ($c=G=\rho=1$)]

Рис. 5. Результат расчета релятивистских поправок $S(x,y)$ и $S(x,y,z)$ в системе Wolfram Mathematica 13.0 с использованием численных методов ($c=G=\rho=1$)
[Figure 5. The result of the calculation of relativistic corrections $S(x,y)$ and $S(x,y,z)$ in the Wolfram Mathematica 13.0 system using numerical methods ($c=G=\rho=1$)]

Из рис. 4 видно, что значение $V$ в центре куба будет минимальным.

Из рис. 5 видно, что значение $S$ в центре куба будет максимальным, что согласуется с физической моделью, т.к. ожидается, что давление $p$ в центре куба будет также максимальным. Полученные распределения $S(x,y)$ и $V(x,y,z)$ соответствуют визуализации, представленной на рис. 2.

Рассмотрим аналитические решения уравнений (20) и (19) для области на некотором расстоянии от куба $R>a$. Уравнения (19) для временной поправки запишем так:
\[ \begin{equation}
V(R)=Gc^{2}\int_{V'}\frac{\rho'}{|R-r'|}dV'+\frac{G}{2}\int_{V'}\frac{\rho' U_{ext}}{|R-r'|}dV',
\end{equation} \tag{21} \]
где $r'=x'_{i}$ — компоненты радиус-вектора, принадлежащие элементу объема куба, $R=x_{i}$ — компоненты радиус-вектора до точки вычисления потенциала.

Для случая, когда $R\gg a$, разложим $|R-r'|$ в ряд по мультполям [28]. Для куба из-за наличия симметрии квадроупольный момент равен $Q_{\alpha\beta}=0$, $Q_{\alpha\beta\gamma}=0$. Вычислим следующий полином — октоуполь. Введем координаты радиус-вектора $r'$ и $R$, причем
\[ \begin{equation*}
r'=n_{i}x'_{i},\quad R=n_{j}x_{j}.
\end{equation*} \]
Октоуполь является четвертым полиномом Лежандра:
\[ \begin{equation*}
O=\dfrac{r^{4}}{8R^{5}}(35\cos^{4}\alpha-30\cos^{2}\alpha+3),
\end{equation*} \]
тогда получим
\[ \begin{equation*}
O=\dfrac{x_{i}x_{l}x_{\alpha}x_{\gamma}}{8R^{9}}(x'_{k}x'_{m}x'_{\mu}x'_{\beta})
\bigl[35\gamma_{il\alpha\gamma km\mu\beta}-30n_{iklm}\delta_{\beta\mu}\delta_{\alpha\gamma}+3\delta_{\beta\mu}\delta_{km}\delta_{il}\delta_{\alpha\gamma}\bigr],
\end{equation*} \]
где $\gamma_{il\alpha\gamma km\mu\beta} = n_{iklm}(\delta_{\alpha\beta}\delta_{\gamma\mu}+\delta_{\alpha\gamma}\delta_{\beta\mu}+\delta_{\alpha\mu}\delta_{\beta\gamma})/225$, $ n_{iklm}= (\delta_{ik}\delta_{
lm}+\delta_{il}\delta_{km} +\delta_{im}\delta_{kl})/15$ — усредненные значения единичного вектора.

Для куба разложение $1/|R-r'|$ в ряд по мультполям имеет следующий вид:
\[ \begin{equation}
\frac{1}{|R-r'|}=\dfrac{1}{R}+\dfrac{(x_{i})^{2}(x_{j})^{2}}{8R^{9}}(x'_{\alpha})^{2}(x'_{\beta})^{2}\gamma,
\end{equation} \tag{22} \]
где $\gamma=0.16(2\delta_{i\alpha}+1)(2\delta_{j\beta}+1)-2(2\delta_{i\alpha}+1)+3$.

Подставляя разложение (22) в (21), получим
\[ \begin{multline*}
V(R)= \frac{Gc^{2}}{R} \int_{V'} \rho'dV' + \frac{Gc^{2} \gamma (x_{i})^{2}(x_{j})^{2}}{8R^{9}} \int _{V'} \rho'(x'_{\alpha})^{2}(x'_{\beta})^{2}dV' + {}\\
{} +\frac{G}{2 R} \int_{V'} \rho'U_{ext}dV' +O(1/R^{9+n}).
\end{multline*} \]

Рассмотрим аналитическое решение уравнения для пространственных поправок (20). Поскольку первый член уравнения (20) зависит от $\rho$, вне куба он обращается в ноль, поэтому первый интеграл конечен. Второй член в (20) не обращается в ноль во всем пространстве, однако на бесконечности обращается в ноль, что также говорит о конечности интеграла. Перепишем уравнение (20) в виде
\[ \begin{equation}
\Delta S=-4\pi G\rho U_{ext}+ \frac{G^{2}}{2} \int_{V'}\int_{V''}\rho'\rho''\nabla'\Bigl(\frac{1}{|R-r'|}\Bigr)\nabla''\Bigl(\dfrac{1}{|R-r''|}\Bigr)dV'dV'' ,
\end{equation} \tag{23} \]
где $r'$ и $r''$ — радиус-векторы, принадлежащие элементам объема $V'$ и $V''$ куба.

Решение (23) принимает вид
\[ \begin{equation*}
S(R)= \frac G R \int_ V \rho U_{ext} dV +O(1/R^{9+n})+
\frac{G^{2}}{2} \int_{V'}\int_{V''}\rho'\rho''\frac{\partial^{2} (\ln{z})}{\partial x'_{i}\partial x''_{i}} dV' dV'',
\end{equation*} \]
где $\ln{z}=\ln(|R-r'|+|R-r''|-|r'-r''|)$ и удовлетворяет уравнению
\[ \begin{equation*}
\Delta \ln{z}=\dfrac{1}{|R-r'||R-r''|}.
\end{equation*} \]
Для нахождения внешнего решения на больших расстояниях от куба ${R\gg a}$ разложим функцию $z$ по обратным степеням $R$. В разложении отбросим члены, убывающие быстрее $1/R^{2}$, и учтем симметрию куба, тогда
\[ \begin{equation*}
\frac{\partial^{2} (\ln{z})}{\partial x'_{i}\partial x''_{i}}=\dfrac{1}{2R^{2}}-\frac{1}{R}\frac{1}{|r'-r''|}.
\end{equation*} \]

Запишем решение для функции $S$ на больших расстояниях от куба $R\gg a$:
\[ \begin{multline}
S(R)= \frac G R \int _V \rho U_{ext} dV +O(1/R^{9+n})+ {}
\\
{} +
\frac{16G^{2}\rho^2 a^{6}}{R^{2}}-
\frac{G^{2}}{2R} \int_{V'}\int_{V''}\frac{\rho'\rho''}{|r'-r''|}dV'dV''.
\end{multline} \tag{24} \]
Интегралы, входящие в $S$ и $V$, являются вполне определенными физическими величинами:
$\qquad \displaystyle W=\frac{1}{2}\int_{V}\rho U dV$ — гравитационная энергия тела [14, 34],
$\qquad \displaystyle M=\int_{V}\rho dV=8\rho a^{3}$ — масса куба (11).

В (24) 6-кратный интеграл можно упростить, используя определение ньютонового потенциала:
\[ \begin{equation*}
\frac{G^{2}}{2R}\int_{V'}\int_{V''}\frac{\rho'\rho''}{|r'-r''|}dV'dV''=\frac{G}{2R}\int_{V'}\rho'U(r')dV'=\frac{GW}{R}.
\end{equation*} \]
Таким образом, решение $S$ и $V$ снаружи куба имеет вид
\[ \begin{multline}
V(R) =\dfrac{GM}{R}c^{2}-\dfrac{7Ma^{4}}{18R^{9}}c^{2} (x^{4}+y^{4}+z^{4}+3x^{2}y^{2}+3y^{2}z^{2}+3z^{2}x^{2})+ {}
\\
{}+ \frac{GW}{R}+O(1/R^{9+n});
\end{multline} \tag{25} \]
\[ \begin{equation}
S(R)=\frac{G^{2}M^{2}}{4R^{2}}+O(1/R^{9+n}).
\end{equation} \tag{26} \]
Компоненты метрического тензора во внешнем пространстве будут иметь вид
\[ \begin{multline}
g^{00} \approx\mathcal{G}^{00}c=1+\dfrac{4GM}{R}-\dfrac{3Ma^{4}}{2R^{9}}A(x,y,z)+\frac{4GW}{Rc^{2}} +{}
\\
{}+\frac{7}{c^{4}}\Bigl(\dfrac{GM}{R}c^{2}-\dfrac{7Ma^{4}}{18R^{9}}A(x,y,z)+\frac{GW}{R}\Bigr)^{2},
\end{multline} \tag{27} \]
\[ \begin{equation}
g^{ik} \approx\mathcal{G}^{ik}c=-c^{2}\delta_{ik}\Bigl(1-\frac{G^{2}M^{2}}{R^{2}c^{4}}\Bigr),
\end{equation} \tag{28} \]
где $A(x,y,z)=\left(x^{4}+y^{4}+z^{4}+3x^{2}y^{2}+3y^{2}z^{2}+3z^{2}x^{2}\right)$.

Рассмотрим параметры гравитационного поля однородного куба, полученные с помощью формул (6) и (27) и отраженные на рис. 6.

Рис. 6. Сравнение параметров гравитационного поля во внешнем пространстве без релятивистских поправок и с их наличием от координаты $x$, причем $c=G=\rho=1$; вертикальные линии представляют собой границы куба: $a=1$
[Figure 6. Comparison of the gravitational field parameters in outer space without relativistic corrections and with their presence as functions of the $x$ coordinate, with $c=G=\rho=1$; the vertical lines represent the boundaries of the cube: $a=1$]

На рис. 6 очевиден вклад релятивистских поправок. Стоит обратить внимание, что функция $g^{00}(x)$ убывает с большей скоростью, чем функция $U_{ext}(x)$, а также имеет большее значение функции при равных значениях аргумента $x$. Черные линии, представляющие границы куба, являются также асимптотой к функции $g^{00}(x)$. Полученные результаты с точностью достаточной для нашей задачи, определяют параметры гравитационного поля для однородного куба, рассмотренного в рамках релятивистского подхода.

Выводы

В работе проведен литературный обзор расчета потенциала гравитационного поля однородного куба в рамках ньютоновской механики. Для релятивистского случая описана физическая модель, в которой координатный куб в пространстве (9) заполнен несжимаемой жидкостью с нулевой скоростью и плотностью $\rho$. Методом постньютоновского формализма (в первом и во втором приближении) аналитически описана временная и пространственная поправка $S$ и $V$ для области вне куба (см. формулы (25) и (26)), а для области внутри куба решение получено с помощью численных методов. Получены компоненты метрического тензора во внешнем пространстве в первом (см. формулы (14)) и во втором (см. формулы (27) и (28)) приближении.

Конкурирующие интересы. Конкурирующих интересов не имеется.
Авторский вклад и ответственность. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.
Финансирование. Исследование выполнялось без финансирования.
Благодарности. Авторы статьи выражают благодарность Авдееву Вячеславу Юрьевичу (Физический институт им. П.Н. Лебедева РАН, лаборатория математических методов обработки, Москва, Россия), Неясову Петру Петровичу (Оренбургский государственный университет, кафедра физики и методики преподавания физики, Оренбург, Россия) за бесценные советы по написанию рукописи статьи; Иванчику Александру Владимировичу и Барсукову Дмитрию Петровичу (ФТИ им. А.Ф. Иоффе РАН, сектор теоретической астрофизики, Санкт-Петербург, Россия) за рекомендации в решении задачи.

About the authors

Valery N. Makarov

Orenburg State University

Author for correspondence.
Email: makarsvet13@gmail.com
ORCID iD: 0000-0001-5749-1427
https://www.mathnet.ru/person190858

Cand. Phys. & Math. Sci.; Associate Professor; Dept. of Physics and Methodology of Teaching Physics

Russian Federation, 460018, Orenburg, pr. Pobedy, 13

Leonid A. Shleiger

Ioffe Physical-Technical Institute, Russian Academy of Sciences

Email: lslejger@gmail.com
ORCID iD: 0009-0007-9648-3172
https://www.mathnet.ru/person191490

Junior Researcher; Theoretical Astrophysics Section

Russian Federation, 194021, St. Petersburg, Polytekhnicheskaya str., 26

Alexander A. Karasev

Institute of Industrial Ecology, Ural Branch of RAS; Ural Federal University named after the First President of Russia B.N. Yeltsin

Email: KarasevAL@ecko.uran.ru
ORCID iD: 0000-0001-7394-7375
https://www.mathnet.ru/person130772

младший научный сотрудник; лаб. эколого-климатических проблем Арктики; ассистент; каф. прикладной математики и механики

Russian Federation, 620108, Ekaterinburg, Sofya Kovalevskaya str., 20; 620002, Ekaterinburg, Mira str., 19

References

Supplementary files