Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <limits>
9 : #include <optional>
10 : #include <pup.h>
11 : #include <vector>
12 :
13 : #include "DataStructures/DataBox/Prefixes.hpp"
14 : #include "DataStructures/TaggedTuple.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "Elliptic/Systems/SelfForce/Scalar/Tags.hpp"
17 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
18 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
19 : #include "Options/Auto.hpp"
20 : #include "Options/String.hpp"
21 : #include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
22 : #include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
23 : #include "Utilities/Gsl.hpp"
24 : #include "Utilities/TMPL.hpp"
25 :
26 0 : namespace ScalarSelfForce::AnalyticData {
27 :
28 : /*!
29 : * \brief Scalar self force for a scalar charge on a circular equatorial orbit.
30 : *
31 : * This class implements Eq. (2.9) in \cite Osburn:2022bby . It does so by
32 : * defining the background fields $\alpha$, $\beta$, and $\gamma_i$ in the
33 : * general form of the equations
34 : * \begin{equation}
35 : * -\partial_i F^i + \beta \Psi_m + \gamma_i \partial_i \Psi_m = S_m
36 : * \text{.}
37 : * \end{equation}
38 : * with the flux [Eq. (2.39) in \cite Vu:2026ypc]
39 : * \begin{equation}
40 : * F^i = \{\frac{\Delta}{r^2+a^2}\partial_{r}, \frac{1-z^2}{r^2+a^2}
41 : \partial_{z}\} \Psi_m
42 : * \text{.}
43 : * \end{equation}
44 : * We make the following changes compared to \cite Osburn:2022bby :
45 : *
46 : * - Multiply by the factor $\Sigma^2 / (r^2 + a^2)^2$ so that we can easily
47 : * write the equations in first-order flux form.
48 : * - Use $z = \cos\theta$ as angular coordinate instead of $\theta$. This avoids
49 : * the $\cot\theta$ term by rewriting the angular derivatives as:
50 : * \f[
51 : * \partial_\theta^2+\cot\theta\partial_\theta =
52 : * \partial_{z}\sin^2\theta\partial_{z}
53 : * \f]
54 : * - Decompose $\Psi_m = \sin(\theta)^m u_m(r_\star, \theta)$. This avoids the
55 : * $m^2/\sin^2\theta$ term by factoring out the singular behavior at the
56 : * poles. The equations transform as:
57 : * \f[
58 : * -\partial_{z}\sin^2\theta\partial_{z} \Psi_m +
59 : * \frac{m^2}{\sin^2\theta}\Psi_m = \sin(\theta)^m \left( m(m+1)
60 : * + 2m z \partial_{z}
61 : * - \partial_{z}\sin^2\theta\partial_{z} \right) u_m
62 : * \f]
63 : * We divide by $\sin(\theta)^m$ to get the equations for $u_m$.
64 : * With the above three changes, Eq. (2.9) in \cite Osburn:2022bby becomes
65 : * Eq. (2.19) in \cite Vu:2026ypc .
66 : * - Use $r$ as the radial coordinate instead of $r_\star$. This has two
67 : * advantages: first, the horizon is placed at a finite radius rather than
68 : * $r_\star \rightarrow -\infty$; second, the flux component normal to the
69 : * boundary vanishes at the horizon ($r=r_+$), which reduces to
70 : * simple regularity conditions.
71 :
72 : * Written this way, the equations are regular at the poles and converge
73 : * exponentially. We also don't have to apply angular boundary conditions
74 : * because regularity at the poles is automatically enforced by the
75 : * $\sin^2\theta$ factor in the flux.
76 : *
77 : * The resulting factors in the equation are:
78 : * \begin{align}
79 : * &\beta = \left(\frac{1}{\Delta}\left(-m^2\Omega^2 \Sigma^2 + 4a m^2 \Omega M
80 : r\right) +
81 : * m (m + 1) + \frac{2M}{r}(1-\frac{a^2}{Mr}) + \frac{2iam}{r}
82 : * \right) \frac{1}{r^2 + a^2} \\
83 : * &\gamma_{r} = -\frac{2iam}{r^2+a^2} + \frac{2 a^2 \Delta}{r(r^2+a^2)^2} \\
84 : * &\gamma_{z} = \frac{2 m z}{r^2 + a^2}
85 : * \end{align}
86 : *
87 : * This class also provides the effective source $S_m^\mathrm{eff} = \Delta_m
88 : * \Psi_m^P$ and the singular field $\Psi_m^P$ in the regularized region (see
89 : * `ScalarSelfForce::FirstOrderSystem` and Sec. III in \cite Vu:2026ypc ).
90 : * The effective source is computed using the scalar EffectiveSource code by
91 : * Wardell et. al. (https://github.com/barrywardell/EffectiveSource and
92 : * \cite Wardell:2011gb ). It is then transformed to correspond to the m-mode
93 : * decomposition used in \cite Vu:2026ypc Eq. (2.16) as:
94 : * \begin{align}
95 : * \Psi_m^P &= \frac{r}{2 \pi \sin(\theta)^{|m|}}
96 : * e^{-i m \left( \varphi - \phi\right)} \Phi_m^\mathrm{Wardell} \\
97 : * S_m^\mathrm{eff} &= \frac{r}{2 \pi \sin(\theta)^{|m|}}
98 : * e^{-i m \left( \varphi - \phi\right)}
99 : * \frac{\Delta\,(r^2 + a^2\cos^2\theta)}{(r^2 + a^2)^2}
100 : * S_m^\mathrm{Wardell}
101 : * \text{,}
102 : * \end{align}
103 : * where $\varphi = \phi + \frac{a}{r_+ - r_-} \ln(\frac{r-r_+}{r-r_-})$
104 : * (Eq. (2.14) in \cite Vu:2026ypc ).
105 : * Compared to Eq. (2.8) in \cite Osburn:2022bby, we additionally divide by
106 : * $\sin(\theta)^{|m|}$ to account for the change of variable from $\Psi_m$ to
107 : * $u_m$ described above.
108 :
109 : *
110 : * \par Impose equatorial symmetry
111 : * Since this work is restricted to equatorial orbits, we can enforce equatorial
112 : * symmetry by reformulating the equations in terms of the coordinate
113 : * $z^2 = \cos^2\theta$ instead of $z = \cos\theta$. This transforms the factors
114 : * of first-order flux formulation as [Eqs. (2.44)-(2.45) in \cite Vu:2026ypc]:
115 : * \begin{equation}
116 : * F^{\cos^2\theta} = 4 \cos^2\theta F^{\cos\theta}
117 : * \text{.}
118 : * \end{equation}
119 : * \begin{equation}
120 : * \gamma_{\cos^2\theta} = 2 \cos\theta \gamma_{\cos\theta} + 2 \alpha
121 : * \text{.}
122 : * \end{equation}
123 : *
124 : * \par Hyperboloidal slicing
125 : * Transforming to a hyperboloidal time coordinate $s = t - h(r_*)$ can simplify
126 : * the solution a lot because the slice will only intersect a finite number of
127 : * wave fronts. In frequency domain this transformation means that partial
128 : * derivatives transform as [Eqs. (2.31)-(2.32) in \cite Vu:2026ypc]:
129 : * \begin{align}
130 : * &\partial_{r_*} \rightarrow \partial_{r_*} + i m\Omega H(r_*) \\
131 : * &\partial_{r_*}^2 \rightarrow \partial_{r_*}^2 +
132 : * 2 i m\Omega H(r_*) \partial_{r_*} + i m\Omega H'(r_*)
133 : * -m^2\Omega^2 H(r_*)^2
134 : * \end{align}
135 : * where $H(r_*) = h'(r_*)$ is the boost function that asymptotes to
136 : * $H(r_* \rightarrow \pm \infty) = \pm 1$.
137 : * This maps to the following additional terms in $\beta$ and $\gamma_{r_*}$:
138 : * \begin{align}
139 : * &\beta \rightarrow \beta + i m\Omega \gamma_{r_*} H(r_*)
140 : * - i m\Omega H'(r_*) + m^2\Omega^2 H(r_*)^2 \\
141 : * &\gamma_{r_*} \rightarrow \gamma_{r_*} - 2 i m\Omega H(r_*)
142 : * \end{align}
143 : * The complete expressions for these factors are given in Eqs. (2.41)-(2.43)
144 : * in \cite Vu:2026ypc .
145 : * For the boost function we use the piecewise-constant "vtu" slicing
146 : * [Eq. (2.34) in \cite Vu:2026ypc]:
147 : * $H = -1$ in the $v$ domain ($r_* \leq r_*^v$), $H = 0$ in the $t$ domain
148 : * ($r_*^v \leq r_* \leq r_*^u$), and $H = +1$ in the $u$ domain
149 : * ($r_* \geq r_*^u$), where $r_*^v < r_{*0} < r_*^u$ and $r_{*0}$ is
150 : * the particle location. The jump in $H$ at the domain interfaces produces a
151 : * jump in the radial derivative of $\Psi_m$, which we impose as junction
152 : * conditions [Eqs. (2.36)-(2.37) in \cite Vu:2026ypc]:
153 : * \begin{align}
154 : * &\left.\frac{\partial \Psi_m}{\partial r_*}\right|_{r_*=(r_*^v)^+} -
155 : * \left.\frac{\partial \Psi_m}{\partial r_*}\right|_{r_*=(r_*^v)^-} =
156 : * -im\Omega\Psi_m\big|_{r_*=r_*^v} \\
157 : * &\left.\frac{\partial \Psi_m}{\partial r_*}\right|_{r_*=(r_*^u)^+} -
158 : * \left.\frac{\partial \Psi_m}{\partial r_*}\right|_{r_*=(r_*^u)^-} =
159 : * -im\Omega\Psi_m\big|_{r_*=r_*^u}
160 : * \end{align}
161 : *
162 : * \par Compactification
163 : * To resolve the large wavezone without truncating the domain at a large
164 : * but finite $R$, we compactify the $u$ region to $r=\infty$.
165 : * With $r_u$ the radius at the $t$-$u$ transition, we redefine
166 : * the domain coordinate as
167 : * \begin{align}
168 : * \sigma = 2 r_u - \frac{r_u^2}{r} \text{,}
169 : * \end{align}
170 : * so $\sigma \rightarrow r_u$ as $r \rightarrow r_u$ and $\sigma
171 : * \rightarrow 2 r_u$ as $r \rightarrow \infty$. Both $\sigma$ and its
172 : * Jacobian $J \equiv d\sigma/dr = r_u^2/r^2$ are continuous ($C^1$) at the
173 : * transition ($J=1$ there), so no coordinate jump condition is needed.
174 : *
175 : * Since $\partial_r = J \partial_\sigma$, the flux form must be rewritten
176 : * in terms of $\sigma$-derivatives. Dividing the transformed equation by
177 : * $J$ to match the generic structure of first order flux form gives the
178 : * substitutions
179 : * \begin{align}
180 : * &\alpha_r \rightarrow \alpha_r J, &&
181 : * \alpha_\theta \rightarrow \alpha_\theta / J, \\
182 : * &\beta \rightarrow \beta / J, &&
183 : * \gamma_\theta \rightarrow \gamma_\theta / J
184 : * \end{align}
185 : * with $S_m = 0$ in the $u$ region. As $r\rightarrow\infty$, $J\rightarrow
186 : * 0$, so the radial flux vanishes, removing the need for the approximate
187 : * (Bayliss-Turkel) outer boundary condition. To avoid $\infty/\infty$
188 : * terms, all coefficients are written as functions of $1/r$ rather than
189 : * $r$.
190 : */
191 1 : class CircularOrbit : public elliptic::analytic_data::Background,
192 : public elliptic::analytic_data::InitialGuess {
193 : public:
194 0 : struct BlackHoleMass {
195 0 : static constexpr Options::String help =
196 : "Kerr mass parameter 'M' of the black hole";
197 0 : using type = double;
198 : };
199 0 : struct BlackHoleSpin {
200 0 : static constexpr Options::String help =
201 : "Kerr dimensionless spin parameter 'chi' of the black hole";
202 0 : using type = double;
203 : };
204 0 : struct OrbitalRadius {
205 0 : static constexpr Options::String help =
206 : "Radius 'r_0' of the circular orbit";
207 0 : using type = double;
208 : };
209 0 : struct MModeNumber {
210 0 : static constexpr Options::String help =
211 : "Mode number 'm' of the scalar field";
212 0 : using type = int;
213 : };
214 0 : struct HyperboloidalSlicingTransitions {
215 0 : static constexpr Options::String help =
216 : "Enable hyperboloidal slicing by specifying the transition points for "
217 : "the boost function. The boost function transitions from -1 to zero "
218 : "between the first two points and from zero to 1 between the last "
219 : "two points. The two points can be the same, in which case the "
220 : "transition "
221 : "is discontinuous (vtu slicing) and the jump must be handled by the DG "
222 : "scheme (see option 'NullSlicingBlocks'). "
223 : "The effective source can only be evaluated where the "
224 : "boost function is zero, so the regularized region must be between "
225 : "the second and third points.";
226 0 : using type = Options::Auto<std::array<double, 4>, Options::AutoLabel::None>;
227 : };
228 0 : struct PenetratingHorizon {
229 0 : static constexpr Options::String help =
230 : "If 'False', use tortoise radial coordinate where the Kerr horizon is "
231 : "at negative infinity. If 'True', use Boyer-Lindquist radial "
232 : "coordinate where the Kerr horizon is at r_+.";
233 0 : using type = bool;
234 : };
235 0 : struct ImposeEquatorialSymmetry {
236 0 : static constexpr Options::String help =
237 : "Impose symmetry across the equatorial plane by using cos(theta)^2 "
238 : "as the angular coordinate instead of cos(theta). This means the "
239 : "domain should span [0, 1] instead of [-1, 1].";
240 0 : using type = bool;
241 : };
242 0 : using options = tmpl::list<BlackHoleMass, BlackHoleSpin, OrbitalRadius,
243 : MModeNumber, HyperboloidalSlicingTransitions,
244 : PenetratingHorizon, ImposeEquatorialSymmetry>;
245 0 : static constexpr Options::String help =
246 : "Quasicircular orbit of a scalar point charge in Kerr spacetime";
247 :
248 0 : CircularOrbit() = default;
249 0 : CircularOrbit(const CircularOrbit&) = default;
250 0 : CircularOrbit& operator=(const CircularOrbit&) = default;
251 0 : CircularOrbit(CircularOrbit&&) = default;
252 0 : CircularOrbit& operator=(CircularOrbit&&) = default;
253 0 : ~CircularOrbit() override = default;
254 :
255 0 : CircularOrbit(
256 : double black_hole_mass, double black_hole_spin, double orbital_radius,
257 : int m_mode_number,
258 : std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions,
259 : bool penetrating_horizon, bool impose_equatorial_symmetry);
260 :
261 0 : explicit CircularOrbit(CkMigrateMessage* m);
262 : using PUP::able::register_constructor;
263 0 : WRAPPED_PUPable_decl_template(CircularOrbit);
264 :
265 0 : tnsr::I<double, 2> puncture_position() const;
266 0 : double black_hole_mass() const { return black_hole_mass_; }
267 0 : double black_hole_spin() const { return black_hole_spin_; }
268 0 : double orbital_radius() const { return orbital_radius_; }
269 0 : int m_mode_number() const { return m_mode_number_; }
270 0 : double omega() const;
271 0 : std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions()
272 : const {
273 : return hyperboloidal_slicing_transitions_;
274 : }
275 0 : bool penetrating_horizon() const { return penetrating_horizon_; }
276 0 : bool impose_equatorial_symmetry() const {
277 : return impose_equatorial_symmetry_;
278 : }
279 :
280 0 : using background_tags = tmpl::list<Tags::Alpha, Tags::Beta, Tags::Gamma>;
281 0 : using source_tags = tmpl::list<
282 : ::Tags::FixedSource<Tags::MMode>, Tags::SingularField,
283 : ::Tags::deriv<Tags::SingularField, tmpl::size_t<2>, Frame::Inertial>,
284 : Tags::BoyerLindquistRadius>;
285 :
286 : // Background
287 0 : tuples::tagged_tuple_from_typelist<background_tags> variables(
288 : const tnsr::I<DataVector, 2>& x, background_tags /*meta*/) const;
289 :
290 : // Initial guess
291 0 : static tuples::TaggedTuple<Tags::MMode> variables(
292 : const tnsr::I<DataVector, 2>& x, tmpl::list<Tags::MMode> /*meta*/);
293 :
294 : // Fixed sources
295 0 : tuples::tagged_tuple_from_typelist<source_tags> variables(
296 : const tnsr::I<DataVector, 2>& x, source_tags /*meta*/) const;
297 :
298 : template <typename... RequestedTags>
299 0 : tuples::TaggedTuple<RequestedTags...> variables(
300 : const tnsr::I<DataVector, 2>& x, const Mesh<2>& /*mesh*/,
301 : const InverseJacobian<DataVector, 2, Frame::ElementLogical,
302 : Frame::Inertial>& /*inv_jacobian*/,
303 : tmpl::list<RequestedTags...> /*meta*/) const {
304 : return variables(x, tmpl::list<RequestedTags...>{});
305 : }
306 :
307 : // NOLINTNEXTLINE
308 0 : void pup(PUP::er& p) override;
309 :
310 : private:
311 0 : friend bool operator==(const CircularOrbit& lhs, const CircularOrbit& rhs);
312 :
313 0 : double black_hole_mass_{std::numeric_limits<double>::signaling_NaN()};
314 0 : double black_hole_spin_{std::numeric_limits<double>::signaling_NaN()};
315 0 : double orbital_radius_{std::numeric_limits<double>::signaling_NaN()};
316 0 : int m_mode_number_{};
317 0 : std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions_{};
318 0 : bool penetrating_horizon_{false};
319 0 : bool impose_equatorial_symmetry_{false};
320 : };
321 :
322 0 : bool operator!=(const CircularOrbit& lhs, const CircularOrbit& rhs);
323 :
324 : } // namespace ScalarSelfForce::AnalyticData
|
