Угол отклонения найден из равенства $$\tan\alpha\approx\alpha=\frac{p_y}{p_x}.,$$ в предположении, что $\alpha\ll1$. Далее можно найти $$p_y=\int F_ydt,$$ и из закона гравитации $$F_y=\frac{GMm}{b^2}\cos^3\varphi.$$ Заменим переменную $$dt=\frac{dx}{v}=\frac{b}{v}\frac{d\varphi}{\cos^2\varphi}$$ и получим $$p_y=\frac{GMm}{bv}\int_{-\pi/2}^{\pi/2}\cos\varphi d\varphi=\frac{2GMm}{bv}.$$ Здесь мы предполагаем, что тело движется вдоль прямой, так как $\alpha\ll1$. Поэтому $$\alpha=\frac{p_y}{p}$$ и $$\alpha=\frac{2GM}{bv^2}=\frac{2b_1}{b}, \quad k=2$$

При пролете мимо массивного тела энергия звезды остается постоянной. Отсюда $$(p+\Delta p_x)^2+p_y^2=p^2.$$ Мы знаем, что $p_y\ll p$, поэтому $$\Delta p_x=-\frac{p_y^2}{2p}=-\frac{\alpha^2}{2}p=-\frac{2 G^2 M^2 m}{b^2 v^3}$$

Чтобы посчитать среднюю силу, мы можем проинтегрировать вклады звезд с разными прицельными параметрами. Количество звезд, пролетающих за время $\Delta t$ равно $\Delta N=2\pi bvn \,d b\Delta t$. Поэтому сила, тормозящая объект по оси $x$ равна $$F_{DF}=\frac{1}{\Delta t}\int\Delta p_x \,dN=-4\pi G^2M^2\frac{nm}{v^2}\int_{b_{min}}^{b_{max}}\frac{db}{b}=-4\pi G^2 M^2 \frac{\rho}{v^2} \ln \Lambda.$$ Формулы выше применимы только для $b>b_1$ поэтому $b_{min}=b_1$. Верхний предел определяется размером галактики $R=b_{max}$. Таким образом, мы имеем $$F_{D F}=-4 \pi G^2 M^2 \frac{\rho}{v^2} \log \Lambda,$$ где $\Lambda=R/b_1$.

Из второго закона Ньютона $$\frac{Mv^2}{a}=\frac{GM^2}{4a^2}.$$ Отсюда орбитальная скорость $$v_{bin}=\sqrt{\frac{GM}{4a}}.$$ Энергия системы равна $$E=E_\text{kin}+U=2\cdot\frac{Mv^2}{2}-\frac{GM^2}{2a}$$. Подставляя $v$, получаем $$E=-\frac{GM^2}{4a}$$

Из закона сохранения момента импульса $$b\sigma=r_mv_0.$$ Выразим $v_0$. Для этого запишем закон сохранения энергии $$\frac{\sigma^2}{2}=\frac{v_0^2}{2}-\frac{GM_2}{r_m}.$$ Подставим и получим $$
b=r_m \sqrt{1+\frac{2 G M_2}{\sigma^2 r_m}} .
$$

Чтобы оценить время между столкновениями воспользуемся аналогией с газом. Как известно из МКТ, $\pi r^2vn\Delta t=1$, где $r$ – радиус молекул, $v$ – скорость теплового движения, $n$ - молекулярная концентрация, а $\Delta t$ - время между столкновениями одной молекулы с другими. В нашей задаче $b_{max}$ эквивалентен радиусу молекулы. Поэтому для оценки можно записать $$(\Delta t)^{-1}=\pi \sigma b^2_{max}n.$$ Оценим максимальный прицельный параметр $b_{max}$, соответствующий столкновению звезды и бинарной системы. Звезда должна достичь расстояния $a$ до бинарной системы для столкновения. На больших расстояниях от бинарной системы звезда взаимодействует с ней как с точечным объектом массой $M_2=2M$. Из результата $B2$ предполагая $r_m=a$, получаем $b_{max}=a\sqrt{1+\frac{4GM}{\sigma^2a}}$. Принимая во внимание то, что $\sigma^2\ll \frac{GM}{a}$, упрощаем $$b_{max}=\frac{2}{\sigma}\sqrt{GMa}$$ и получаем $$\Delta t=\frac{m\sigma}{4\pi GM\rho a}$$

За одно соударение энергия звезды увеличивается в среднем на $$\Delta E_{star}=\frac{mv_{bin}^2}{2}-\frac{m\sigma^2}{2}.$$ Энергия бинарной системы уменьшается на ту же величину $\Delta E_{bin}=-\Delta E_{star}$. Учитывая, что $\sigma\ll v_{bin}$, получаем $$\Delta E_{bin}=-\frac{m}{2}v^2_{bin}=\frac{GmM}{8a}.$$ Средняя скорость изменения энергии бинарной системы равна $$\frac{dE}{dt}=\frac{\Delta E}{\Delta t}=-\frac{\pi G^2M^2\rho}{2\sigma}.$$ Из пункта $\text{B1}$ можно получить $$\frac{dE}{dt}=\frac{d}{dt}\left(-\frac{GM^2}{4a}\right)=\frac{GM^2}{4a^2}\frac{da}{dt}.$$ Таким образом, скорость изменения радиуса орбиты может быть оценена как $$\frac{da}{dt}=-\frac{2\pi G\rho a^2}{\sigma}$$

Результат $\text{B4}$ может быть проинтегрирован $$\frac{da}{a^2}=-\frac{2\pi G\rho}{\sigma}dt.$$ Чтобы уменьшить радиус вдвое, требуется время $$T_{S S}=\frac{\sigma}{2 \pi G \rho a_1}=7.3 \times 10^{-4} \mathrm{~Gy}$$

Используя, что $$\omega=\frac{v_{bin}}{a}=\sqrt{\frac{GM}{4a^3}}$$ и формулы из условия, можно получить $$\frac{dE}{dt}=-\frac{1024 \cdot 4}{5}\cdot\frac{GM^2v_{bin}^6}{c^5a^2}=-\frac{64}{5}\cdot\frac{G^4M^5}{c^5a^5}.$$ Используя это и связь $\frac{dE}{dt}$ и $\frac{da}{at}$ можно получить $$\frac{da}{dt}=-\frac{256}{5}\cdot\frac{G^3M^3}{c^5a^3}$$

Интегрирование результата $\text{C1}$ дает $$a^3\,da=-\frac{256}{5}\cdot\frac{G^3M^3}{c^5} \Rightarrow\frac{a_2^4-r_g^4}{4}=\frac{256}{5}\cdot\frac{G^3M^3}{c^5}\cdot T_{GW}.$$ Принимая во внимание, что $a_2\gg r_g$ получаем конечный ответ $$T_{GW}=\frac{5}{1024}\cdot\frac{a^4_2c^5}{G^3M^3}$$

Из результата $\text{C2}$ с учетом того, что $T_{GW}=t_H$ получаем $$a_H=\sqrt[4]{\frac{1024}{5} \cdot \frac{G^3 M^3 t_H}{c^5}}=0.098 \mathrm{pc}$$

Галактика сферически симметрична, поэтому масса, расположенная внутри сферы радиусом $r$ равна $$m(r)=\int_0^r4\pi x^2\rho(x)\,dx=\frac{\sigma^2r}{G}.$$ Таким образом, ускорение свободного падения в гравитационном поле звезд равно $$g(r)=\frac{Gm(r)}{r^2}=\frac{\sigma^2}{r}.$$ Так как скорость тела определена уравнением $$\frac{v^2}{r}=g=\frac{\sigma^2}{r},$$ мы можем получить, что $$v=\sigma.$$ То есть скорость постоянна.

Энергия черной дыры в этом гравитационном поле равна $$E=\frac{M\sigma^2}{2}+U.$$ Так как кинетическая энергия постоянна и $$\frac{dE}{dt}=\frac{dU}{dt}=\frac{dU}{da}\frac{da}{dt},$$ из определения потенциальной энергии получаем $$\frac{dU}{da}=g(a)M=\frac{M\sigma^2}{a}.$$ Используя результат $\text{A3}$ получаем $$\frac{dE}{dt}=-F_{DF}v=-4\pi G^2M^2\frac{\rho(a)}{\sigma}\ln\Lambda=-\frac{GM^2\sigma \ln \Lambda}{a^2}.$$ Таким образом, получаем ответ $$\frac{da}{dt}=-\frac{GM\log\Lambda}{a\sigma}$$

$\textbf{Способ 1:}$ Для оценки радиуса $a_1$ можем предположить, что в момент образования бинарной системы масса звезд внутри сферы равна $M$ $$m(a)=\frac{\sigma^2a}{G}=M.$$ Отсюда $$a_1=\frac{GM}{\sigma^2}=10.8~\mathrm{pc}$$

$\textbf{Способ 2:}$ Сила взаимодействия черных дыр равна силе взаимодействия дыры со звездами $$\frac{Gm(a)}{a^2}=\frac{GM}{4a^2}.$$ Тогда ответ $$a_1=\frac{GM}{4\sigma^2}=2.7~\mathrm{pc}$$

Интегрируя результат $\text{D2}$ получаем $$\frac{a_0^2-a_1^2}{2}=\frac{GM\ln\Lambda}{\sigma}T_1.$$ Учитывая, что $a_1\ll a_0$, получаем $$T_1=\frac{a_0^2 \sigma}{2 G M \log \Lambda}=0.121 \mathrm{~Gy}$$

Потери энергии из-за гравитационных волн начинает доминировать, когда $$\frac{\pi G^2M^2\rho_1}{2\sigma}<\frac{64G^4M^5}{5c^5a^5}.$$ Тогда $$a_2=\frac{128}{5\pi}\cdot\frac{G^2M^3\sigma}{c^5\rho_1}=\frac{512}{5}\cdot\frac{G^3M^3a_1^2}{c^5\sigma}=0.018~\mathrm{pc} \quad (a_2=0.010~\mathrm{pc})$$

Для грубой оценки можно предположить, что на стадии рогатки ($a>a_2$) потери энергии вызваны только эффектом рогатки, поэтому $T_2$ вычисляется аналогично $\text{B5}$ $$T_2\approx\frac{\sigma}{2\pi G\rho_1a_2}=0.063~\mathrm{Gy} \quad (T_2\approx0.0068~\mathrm{Gy}).$$ На этапе гравитационных волн потери связаны только с ними, и $T_3$ вычисляется из $\text{C2}$ $$T_3\approx\frac{5}{1024}\cdot\frac{a_2^4c^5}{G^3M^3}=\frac{1}{8\pi}\cdot\frac{\sigma}{G\rho_1a_2}=0.016~\mathrm{Gy}\quad (T_3\approx0.0017~\mathrm{Gy})$$

Полное время слияния равно $$T_{ev}=T_1+T_2+T_{GW}=0.20~\mathrm{Gy}\quad(T_{ev}=0.13~\mathrm{Gy})$$