Чёрные дыры: Ньютон, Эйнштейн

В этой статье моделируются ключевые эффекты общей теории относительности в окрестности чёрной дыры Шварцшильда. Мы последовательно разбираем:

Скорость убегания — как Ньютон и Эйнштейн по-разному приходят к идее чёрной дыры;
Замедление времени и красное смещение — наблюдаемые следствия сильного гравитационного поля;
Приливные силы («спагеттификация») — механизм разрушения объектов;
Релятивистское радиальное падение — решение уравнений движения в координатах собственного и координатного времени с анимацией трёх панелей;
Дополнительные орбиты — фотонная сфера и последняя стабильная круговая орбита (ISCO);
Таблица реальных объектов — сравнение чёрных дыр разной массы (звёздной, Стрелец A*, M87*, TON 618).

Словарь основных сокращений, используемых в статье:

Сокращение	Расшифровка	Пояснение
ЧД	Чёрная дыра	Объект с гравитационным полем, из которого ничто не может выйти
ОТО	Общая теория относительности	Теория гравитации Эйнштейна (1915)
r_s	Радиус Шварцшильда	Гравитационный радиус, горизонт событий ((r_s = 2GM/c^2))
ISCO	Innermost Stable Circular Orbit	Последняя устойчивая круговая орбита (для Шварцшильда (r = 3r_s))
EHT	Event Horizon Telescope	Телескоп горизонта событий, получивший первые изображения теней ЧД (M87* в 2019, Sgr A* в 2022)
G	Гравитационная постоянная	В коде нормирована (G = 1) (геометризованные единицы)
c	Скорость света	В коде (c = 1)
M	Масса чёрной дыры	В единицах, где (GM/c^2 = 1) → (r_s = 2)
τ	Собственное время	Время по часам, падающим вместе с частицей
t	Координатное время	Время удалённого неподвижного наблюдателя
dτ/dt	Отношение времён	Мера гравитационного замедления времени
z	Красное смещение	Относительное изменение длины волны ((z = \Delta\lambda/\lambda))
ΔL	Размер объекта	Характерная длина (например, рост человека 2 м) для расчёта приливных сил
GPS	Global Positioning System	Спутниковая система, учитывающая релятивистское замедление времени (упоминается как пример наблюдаемого эффекта)
BH	Black Hole	Английское сокращение «чёрная дыра» (используется в аннотациях)

В проекте используются следующие библиотеки Julia:

DifferentialEquations — мощный инструмент для численного решения дифференциальных уравнений. Здесь используется для интегрирования уравнения релятивистского радиального падения (ODEProblem + solve с солвером Tsit5).
Printf — стандартная библиотека для форматированного вывода строк (например, @sprintf в анимации для заголовка с текущими значениями времени).
DataFrames (используется в коде, хотя явно не указан в using в приведённой строке, но подключается далее) — для создания и красивой печати таблицы с параметрами реальных чёрных дыр (масса, радиус Шварцшильда, приливной радиус, замедление времени).

# Pkg.add(["Plots", "DifferentialEquations"])
using DifferentialEquations, Printf, DataFrames
gr()

Plots.GRBackend()

В коде используются геометризованные единицы, где (G = 1) и (c = 1). Это стандартный приём в ОТО: все формулы упрощаются, а расстояния и время выражаются в единицах массы. При (M = 1) радиус Шварцшильда становится (r_s = 2). Такая нормировка не меняет физики, но делает код компактным и наглядным. Для реальных объектов константы «возвращаются» отдельно (как в таблице в конце примера).

G = 1.0      # гравитационная постоянная (нормирована)
c = 1.0      # скорость света (нормирована)
M = 1.0      # масса чёрной дыры (в единицах, где GM/c² = 1 -> r_s = 2)
r_s = 2*G*M/c^2   # радиус Шварцшильда = 2
println("📊 Радиус Шварцшильда: r_s = $r_s")

📊 Радиус Шварцшильда: r_s = 2.0

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

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

Эйнштейн (для фотонов): — релятивистская формула, выведенная из условия, что свет не может покинуть горизонт.

r_newton = range(0.1, 10, length=300)
v_escape_newton = @. sqrt(2*G*M/r_newton)          # классическая формула
r_einstein = range(r_s, 10, length=300)
v_escape_einstein = @. c * sqrt(r_s / r_einstein)  # релятивистская для фотонов

p1 = plot(r_newton, v_escape_newton,
     label="Ньютон",
     linewidth=3, color=:blue,
     xlabel="Расстояние от центра r",
     ylabel="Скорость убегания v / c",
     title="Предсказание чёрных дыр: Ньютон vs Эйнштейн",
     legend=:topright, grid=true, ylim=(0, 2.5))
plot!(p1, r_einstein, v_escape_einstein,
      label="Фотоны",
      linewidth=3, color=:red, ls=:dash)
hline!(p1, [c], ls=:dot, color=:black, label="Скорость света", linewidth=2)
vline!(p1, [r_s], ls=:dash, color=:green, label="Горизонт", linewidth=2)
annotate!(p1, r_s+0.5, 2.2, text("Горизонт событий", :green, :bold, 11))
annotate!(p1, 3.0, 1.8, text("При r < r_s : v_esc > c → ничто не может уйти!", :red, 10))
display(p1)

Что показывает график:

При r > r_s обе скорости различаются, но при r = r_s релятивистская достигает c, а ньютоновская превышает c уже на бóльших расстояниях.
Горизонт событий r_s отмечен вертикальной линией. Для r < r_s ньютоновская кривая даёт — это намёк на невозможность выхода, но точное предсказание чёрной дыры даёт только ОТО.
Горизонтальная линия c — предел скорости света.

Далее описано замедления времени: — отношение собственного времени падающего наблюдателя к координатному времени t удалённого наблюдателя. Формула красного смещения: — относительное увеличение длины волны света, идущего из окрестности чёрной дыры.

r_time = range(r_s*1.01, 10, length=300)
time_dilation = @. sqrt(1 - r_s / r_time)           # dτ/dt
redshift = @. 1 / sqrt(1 - r_s / r_time) - 1;       # z = Δλ/λ

p2 = plot(r_time, time_dilation,
     linewidth=3, color=:purple,
     xlabel="Расстояние r", ylabel="Отношение dτ/dt",
     title="Гравитационное замедление времени",
     legend=:bottomright, grid=true, ylim=(0, 1.1))
hline!(p2, [1.0], ls=:dot, color=:green, label="Вдали от ЧД (время нормально)", linewidth=1)
vline!(p2, [r_s], ls=:dash, color=:red, label="Горизонт (время останавливается)", linewidth=2)
annotate!(p2, r_s+0.3, 0.1, text("На горизонте: время останавливается!", :red, :bold, 9))
annotate!(p2, 5.0, 0.9, text("Вдали от ЧД:время идётобычно", :green, 9))
annotate!(p2, 2.0, 0.7, text("При r = 1.5r_s: 1 час near BH = 1.7 часа на Земле", :purple, 8))
display(p2)

По вертикали — (от 0 до 1). Вдали от ЧД значение стремится к 1 (время идёт одинаково). При приближении к горизонту знаменатель , поэтому — время для падающего наблюдателя останавливается с точки зрения удалённого. Аннотации поясняют: на горизонте время «останавливается», вдали — нормально, а на один час вблизи ЧД соответствует 1.7 часа на Земле. Горизонт отмечен вертикальной красной линией, на нём время растягивается до бесконечности, отсюда и название горизонт событий.

p3 = plot(r_time, redshift,
     linewidth=3, color=:orange,
     xlabel="Расстояние r", ylabel="Красное смещение z",
     title="Гравитационное красное смещение",
     legend=:topright, grid=true, ylim=(0, 10))
vline!(p3, [r_s], ls=:dash, color=:red, label="Горизонт (z → ∞)", linewidth=2)
display(p3)

Формула описывает, как гравитация растягивает длину волны света, идущего из окрестности чёрной дыры. При знаменатель стремится к нулю, поэтому (z \to \infty) — свет бесконечно краснеет, а частота падает до нуля. На графике по вертикали — z, по горизонтали — расстояние r. Вертикальная линия на r_s показывает горизонт, где смещение уходит в бесконечность. Вдали от ЧД (смещение исчезает).

Ниже представленны расчёты приливного ускорения: , где м — характерный размер объекта (рост человека). Она показывает разницу ускорения свободного падения между головой и ногами (или между двумя концами объекта). Чем меньше r, тем быстрее растёт . Условный порог разрушения задан как (в относительных единицах). Из этого уравнения находится радиус спагеттификации: . В коде при объект разрывается.

ΔL = 2.0   # рост человека ~2 м
r_tidal = range(r_s*1.01, 20, length=300)
tidal_accel = @. (2 * G * M / r_tidal^3) * ΔL;

p4 = plot(r_tidal, tidal_accel,
     linewidth=3, color=:darkred,
     xlabel="Расстояние r", ylabel="Приливное ускорение (отн. ед.)",
     title="Приливные силы — «спагеттификация»",
     legend=:topright, grid=true)

r_spaghetti = (2*G*M*ΔL / 100)^(1/3) # Условный порог разрушения
vline!(p4, [r_spaghetti], ls=:dot, color=:blue, label="Порог разрушения")
annotate!(p4, r_spaghetti+0.5, maximum(tidal_accel)*0.8,
          text("При r < $(round(r_spaghetti,digits=2)) объект разрывается", :blue, 8))
display(p4)

По вертикали — приливное ускорение (в относительных единицах), по горизонтали — расстояние r от центра чёрной дыры. Кривая резко растёт при приближении к горизонту. Вертикальной пунктирной линией отмечен порог разрушения. Область левее этой линии подписана: «объект разрывается». График наглядно демонстрирует, почему спагеттификация наиболее опасна вблизи малых чёрных дыр (для сверхмассивных порог находится внутри горизонта, как показано в таблице в конце примера).

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

где — собственное время падающего наблюдателя. Знак минус указывает на движение к центру.

Функция radial_fall_full! вычисляет эту производную с защитой от отрицательного подкоренного выражения (которое может возникнуть при из-за численных погрешностей): если подкоренное выражение становится отрицательным, производная принудительно полагается равной нулю. Интегрирование выполняется методом Tsit5 (адаптивный Рунге–Кутта 5-го порядка) в интервале собственного времени с начальным условием . Решение сохраняется с шагом saveat=0.05, а максимальный шаг интеграции ограничен dtmax=0.1.

После получения дискретной последовательности вычисляется координатное время из соотношения:

Численно это делается методом средних значений (midpoint rule): для каждого шага вычисляется среднее радиальное положение . Пока , приращение координатного времени добавляется к накопленному значению. Как только частица достигает горизонта событий (), накопление прекращается, и далее координатное время остаётся постоянным.

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

function radial_fall_full!(du, u, p, τ)
    r = u[1]
    E = p[1]
    r_s = p[2]
    radicand = E^2 - (1 - r_s/r)
    if radicand < 0
        du[1] = 0.0
    else
        du[1] = -sqrt(radicand)   
    end
end
E = 1.0                    
u0 = [10.0]                
τspan = (0.0, 40.0)        
prob_rad = ODEProblem(radial_fall_full!, u0, τspan, (E, r_s))
sol_rad = solve(prob_rad, Tsit5(), saveat=0.05, dtmax=0.1)
r_vals = sol_rad[1,:]
τ_vals = sol_rad.t
t_vals = zeros(length(τ_vals))
last_t = 0.0
for i in 2:length(τ_vals)
    r_mid = (r_vals[i-1] + r_vals[i]) / 2
    if r_mid > r_s
        dt_dτ = E / (1 - r_s / r_mid)
        last_t += dt_dτ * (τ_vals[i] - τ_vals[i-1])
        t_vals[i] = last_t
    else
        t_vals[i] = last_t
    end
end
println("Релятивистское падение:")
println("   Собственное время до горизонта: $(round(τ_vals[end], digits=2))")
println("   Координатное время до горизонта: $(round(t_vals[end], digits=2)) (расходится)")

Релятивистское падение:
   Собственное время до горизонта: 40.0
   Координатное время до горизонта: 27.95 (расходится)

Далее представлена интерактивная визуализация приближения тела к ЧД.

Первый график (Радиальное падение) – показывает, как расстояние от падающего объекта до центра чёрной дыры уменьшается с течением собственного времени (времени по часам самого падающего). Горизонтальная красная пунктирная линия – это горизонт событий. Видно, что объект пересекает горизонт за конечное собственное время и продолжает падать к центру.

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

Третий график (Связь времён) – сравнивает координатное время (по часам удалённого наблюдателя) и собственное время падающего. Пунктирная серая линия показывает, как они совпадали бы в отсутствие чёрной дыры. Реальная зелёная кривая демонстрирует, что по мере приближения к горизонту собственное время всё сильнее отстаёт. После пересечения горизонта координатное время практически перестаёт расти (стремится к конечному пределу), тогда как собственное время продолжает увеличиваться. Аннотация в левом нижнем углу указывает, что горизонт пересекается при собственном времени примерно 13,55.

dilation_vals = [r > r_s ? sqrt(1 - r_s/r) : 0.0 for r in r_vals]

anim = @animate for i in 1:5:length(τ_vals)
    r_cur = r_vals[i]
    τ_cur = τ_vals[i]
    t_cur = t_vals[i]
    dil_cur = dilation_vals[i]
    title_str = @sprintf "Падение в ЧД: t = %.2f, τ = %.2f, dτ/dt = %.3f" t_cur τ_cur dil_cur
    
    p1 = plot(τ_vals[1:i], r_vals[1:i],
              color=:blue, linewidth=2,
              label="r(τ)",
              xlabel="Собственное время τ",
              ylabel="Радиус r",
              legend=:bottomleft,
              title=title_str,
              ylim=(0, maximum(r_vals)+0.5))
    hline!(p1, [r_s], ls=:dash, color=:red, label="Горизонт r_s")
    scatter!(p1, [τ_cur], [r_cur], color=:orange, ms=6, label="Частица")

    p2 = plot(τ_vals[1:i], dilation_vals[1:i],
              color=:purple, linewidth=2,
              label="dτ/dt",
              xlabel="Собственное время τ",
              ylabel="dτ/dt",
              legend=:bottomleft,
              ylim=(0, 1.1))
    hline!(p2, [0], ls=:dot, color=:red, label="Горизонт и внутри")
    scatter!(p2, [τ_cur], [dil_cur], color=:red, ms=6, label="Текущее значение")

    p3 = plot(t_vals[1:i], τ_vals[1:i],
              color=:green, linewidth=2,
              label="τ(t)",
              xlabel="Координатное время t",
              ylabel="Собственное время τ",
              legend=:bottomleft)
    plot!(p3, t_vals[1:i], t_vals[1:i], ls=:dot, color=:gray, label="τ = t (без ЧД)")
    scatter!(p3, [t_cur], [τ_cur], color=:orange, ms=6, label="Частица")
    x_min = minimum(t_vals[1:i])
    y_min = minimum(τ_vals[1:i])
    annotate!(p3, x_min + 0.05*(maximum(t_vals[1:i])-x_min),
                   y_min + 0.05*(maximum(τ_vals[1:i])-y_min),
              text("Горизонт пересечён при τ ≈ 13.55", :black, :left, 8))

    plot(p1, p2, p3, layout=(3,1), size=(800, 900))
end

gif(anim, "relativistic_fall_full.gif", fps=20)

Ниже представлен график эффективного потенциала он показывает, как «энергетический барьер» зависит от расстояния до чёрной дыры для частицы, движущейся по орбите. По вертикали — эффективная потенциальная энергия, по горизонтали — расстояние. Минимумы на кривой соответствуют устойчивым круговым орбитам, максимумы — неустойчивым.

Фотонная сфера (голубая вертикальная линия) — расстояние, на котором свет может двигаться по круговой орбите вокруг чёрной дыры. Эта орбита неустойчива: малейшее возмущение — и фотон либо упадёт на дыру, либо улетит прочь.
ISCO (последняя устойчивая круговая орбита) (розовая вертикальная линия) — самый внутренний радиус, на котором возможно стабильное круговое движение массивной частицы. Ближе к дыре любая круговая орбита становится неустойчивой, и частица быстро падает внутрь. ISCO играет ключевую роль в астрофизике: именно на этом радиусе заканчивается аккреционный диск, а вещество начинает безвозвратно падать в чёрную дыру, выделяя огромную энергию.
Горизонтальная серая линия — уровень энергии частицы, покоящейся на бесконечности. Где эта линия пересекается с кривой потенциала, там возможны точки поворота движения.

r_photon = 1.5 * r_s
r_ISCO = 3.0 * r_s

println("Дополнительные орбиты:")
println("   Фотонная сфера: r = $r_photon")
println("   Последняя стабильная орбита (ISCO): r = $r_ISCO")

Дополнительные орбиты:
   Фотонная сфера: r = 3.0
   Последняя стабильная орбита (ISCO): r = 6.0

L_isco = 2.0 * sqrt(3) * M
r_pot = range(r_s+0.01, 20, length=500)
V_eff = @. (1 - r_s/r_pot) * (1 + L_isco^2 / r_pot^2)

p5 = plot(r_pot, V_eff,
     label="V_eff(r) = (1 - r_s/r) * (1 + L^2 / r^2)",
     linewidth=2, color=:darkgreen,
     xlabel="Расстояние r", ylabel="Эффективный потенциал V_eff",
     title="Эффективный потенциал: ISCO и фотонная сфера",
     legend=:topright, grid=true)
vline!(p5, [r_photon], ls=:dash, color=:cyan, label="Фотонная сфера (1.5 r_s)")
vline!(p5, [r_ISCO], ls=:dash, color=:magenta, label="ISCO (3 r_s)")
hline!(p5, [1.0], ls=:dot, color=:gray, label="E = 1 (покой на бесконечности)")
display(p5)
savefig("effective_potential.png")
println("✅ График 5 сохранён и отображён: effective_potential.png")

✅ График 5 сохранён и отображён: effective_potential.png

Ниже представлена таблица реальных объектов, на которую мы ссылались выше, в таблице представлены четыре типа чёрных дыр разной массы: звездной массы (10 масс Солнца), Стрелец A* (сверхмассивная дыра в центре нашей Галактики, 4,3 миллиона солнечных масс), M87* (6,5 миллиарда солнечных масс) и TON 618 (одна из самых массивных известных дыр — 66 миллиардов солнечных масс).

Для каждого объекта рассчитаны:

Радиус Шварцшильда r_s (горизонт событий) в километрах.
Приливной радиус — расстояние, на котором приливные силы разрывают объект размером 2 метра (рост человека) с ускорением, равным земному .
Минимальное безопасное расстояние — максимум из приливного радиуса и (чтобы не заходить за горизонт).
Замедление времени на этом минимальном расстоянии — во сколько раз 1 час вблизи дыры длиннее для удалённого наблюдателя.

G_si = 6.67430e-11
c_si = 299792458
M_sun = 1.989e30
ΔL = 2.0          # рост человека, м
g_earth = 9.8     # м/с²

r_s_km(mass_solar) = 2 * G_si * (mass_solar * M_sun) / c_si^2 / 1000 # Радиус Шварцшильда в км
function r_tidal_km(mass_solar)
    M_kg = mass_solar * M_sun
    r_m = (2 * G_si * M_kg * ΔL / g_earth)^(1/3)
    return r_m / 1000
end
function dilation_factor(r_km, r_s_km)
    if r_km <= r_s_km
        return Inf
    else
        return 1 / sqrt(1 - r_s_km / r_km)
    end
end

objects = [
    ("ЧД звёздной массы", 10),
    ("Стрелец A*", 4.3e6),
    ("M87*", 6.5e9),
    ("TON 618", 6.6e10)
]

rows = []
for (name, mass) in objects
    rs = r_s_km(mass)
    r_tidal = r_tidal_km(mass)
    r_min = max(r_tidal, rs * 1.0001)
    dilation = dilation_factor(r_min, rs)
    push!(rows, (name, mass, rs, r_tidal, r_min, dilation))
end

df = DataFrame(
    Symbol("Объект") => [row[1] for row in rows],
    Symbol("M (M☉)") => [row[2] for row in rows],
    Symbol("r_s (км)") => [round(row[3], sigdigits=4) for row in rows],
    Symbol("r_tidal (км)") => [round(row[4], sigdigits=4) for row in rows],
    Symbol("r_min (км)") => [round(row[5], sigdigits=4) for row in rows],
    Symbol("Замедление (1 час → часов)") => [row[6] for row in rows]
)

println("Таблица: для каждого объекта указано минимальное расстояние, на котором можно находиться без разрыва приливными силами (но не ближе 1.0001 r_s), и соответствующее замедление времени.\n")
println(df)

Таблица: для каждого объекта указано минимальное расстояние, на котором можно находиться без разрыва приливными силами (но не ближе 1.0001 r_s), и соответствующее замедление времени.

4×6 DataFrame
 Row │ Объект             M (M☉)   r_s (км)  r_tidal (км)  r_min (км)  Замедление (1 час → часов)
     │ String             Real     Float64   Float64       Float64     Float64
─────┼────────────────────────────────────────────────────────────────────────────────────────────
   1 │ ЧД звёздной массы  10       29.54       8153.0      8153.0                         1.00182
   2 │ Стрелец A*          4.3e6    1.27e7   615300.0         1.27e7                    100.005
   3 │ M87*                6.5e9    1.92e10       7.062e6     1.92e10                   100.005
   4 │ TON 618             6.6e10   1.95e11       1.529e7     1.95e11                   100.005

Ключевые выводы таблицы:

Для чёрной дыры звёздной массы приливной радиус (около 8150 км) значительно больше горизонта (около 30 км). Поэтому астронавт будет разорван приливными силами задолго до подлёта к горизонту. Замедление времени на безопасном минимуме невелико (всего в 1,002 раза).
Для сверхмассивных дыр (Стрелец A*, M87*, TON 618) картина обратная: приливной радиус лежит глубоко внутри горизонта. Это означает, что горизонт можно пересечь, не ощутив спагеттификации — приливные силы на границе чёрной дыры очень слабы. Однако замедление времени на минимальном расстоянии (которое в данном случае ограничено самим горизонтом) колоссально: 1 час вблизи горизонта растягивается до 100 часов для внешнего наблюдателя (коэффициент 100). Это следует из того, что минимальное расстояние практически равно r_s, и знаменатель становится очень малым.

Вывод

Эта статья представляет собой практическое введение в моделирование чёрных дыр в рамках общей теории относительности, в ней продеманстрированно, что ньютоновская теория тяготения, используя формулу скорости убегания , предсказывает: если объект достаточно массивен и компактен, то при скорость убегания превышает скорость света, и свет не может покинуть такую «тёмную звезду» (идея Мичелла, 1783). Однако это лишь одна сторона вопроса, так как в ньютоновской физике свет не имеет массы и не должен подчиняться гравитации в том же смысле. Общая теория относительности Эйнштейна (1915) даёт строгое описание: гравитация есть искривление пространства-времени, а метрика Шварцшильда (1916) предсказывает горизонт событий , за которым ничто, даже свет, не может выйти. Кроме того, ОТО предсказывает качественно новые эффекты: гравитационное замедление времени (время на горизонте останавливается), красное смещение, приливные силы («спагеттификация»), а также существование фотонной сферы и последней устойчивой орбиты ISCO. В коде мы наглядно сравниваем ньютоновскую и релятивистскую скорости убегания, демонстрируя, что только ОТО даёт правильное значение радиуса горизонта и объясняет наблюдаемые явления.