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_{\plus}$), 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 1 : class CircularOrbit : public elliptic::analytic_data::Background,
163 : public elliptic::analytic_data::InitialGuess {
164 : public:
165 0 : struct BlackHoleMass {
166 0 : static constexpr Options::String help =
167 : "Kerr mass parameter 'M' of the black hole";
168 0 : using type = double;
169 : };
170 0 : struct BlackHoleSpin {
171 0 : static constexpr Options::String help =
172 : "Kerr dimensionless spin parameter 'chi' of the black hole";
173 0 : using type = double;
174 : };
175 0 : struct OrbitalRadius {
176 0 : static constexpr Options::String help =
177 : "Radius 'r_0' of the circular orbit";
178 0 : using type = double;
179 : };
180 0 : struct MModeNumber {
181 0 : static constexpr Options::String help =
182 : "Mode number 'm' of the scalar field";
183 0 : using type = int;
184 : };
185 0 : struct HyperboloidalSlicingTransitions {
186 0 : static constexpr Options::String help =
187 : "Enable hyperboloidal slicing by specifying the transition points for "
188 : "the boost function. The boost function transitions from -1 to zero "
189 : "between the first two points and from zero to 1 between the last "
190 : "two points. The effective source can only be evaluated where the "
191 : "boost function is zero, so the regularized region must be between "
192 : "the second and third points.";
193 0 : using type = Options::Auto<std::array<double, 4>, Options::AutoLabel::None>;
194 : };
195 0 : struct PenetratingHorizon {
196 0 : static constexpr Options::String help =
197 : "If 'False', use tortoise radial coordinate where the Kerr horizon is "
198 : "at negative infinity. If 'True', use Boyer-Lindquist radial "
199 : "coordinate where the Kerr horizon is at r_+.";
200 0 : using type = bool;
201 : };
202 0 : struct ImposeEquatorialSymmetry {
203 0 : static constexpr Options::String help =
204 : "Impose symmetry across the equatorial plane by using cos(theta)^2 "
205 : "as the angular coordinate instead of cos(theta). This means the "
206 : "domain should span [0, 1] instead of [-1, 1].";
207 0 : using type = bool;
208 : };
209 0 : using options = tmpl::list<BlackHoleMass, BlackHoleSpin, OrbitalRadius,
210 : MModeNumber, HyperboloidalSlicingTransitions,
211 : PenetratingHorizon, ImposeEquatorialSymmetry>;
212 0 : static constexpr Options::String help =
213 : "Quasicircular orbit of a scalar point charge in Kerr spacetime";
214 :
215 0 : CircularOrbit() = default;
216 0 : CircularOrbit(const CircularOrbit&) = default;
217 0 : CircularOrbit& operator=(const CircularOrbit&) = default;
218 0 : CircularOrbit(CircularOrbit&&) = default;
219 0 : CircularOrbit& operator=(CircularOrbit&&) = default;
220 0 : ~CircularOrbit() override = default;
221 :
222 0 : CircularOrbit(
223 : double black_hole_mass, double black_hole_spin, double orbital_radius,
224 : int m_mode_number,
225 : std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions,
226 : bool penetrating_horizon, bool impose_equatorial_symmetry);
227 :
228 0 : explicit CircularOrbit(CkMigrateMessage* m);
229 : using PUP::able::register_constructor;
230 0 : WRAPPED_PUPable_decl_template(CircularOrbit);
231 :
232 0 : tnsr::I<double, 2> puncture_position() const;
233 0 : double black_hole_mass() const { return black_hole_mass_; }
234 0 : double black_hole_spin() const { return black_hole_spin_; }
235 0 : double orbital_radius() const { return orbital_radius_; }
236 0 : int m_mode_number() const { return m_mode_number_; }
237 0 : double omega() const;
238 0 : std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions()
239 : const {
240 : return hyperboloidal_slicing_transitions_;
241 : }
242 0 : bool penetrating_horizon() const { return penetrating_horizon_; }
243 0 : bool impose_equatorial_symmetry() const {
244 : return impose_equatorial_symmetry_;
245 : }
246 :
247 0 : using background_tags = tmpl::list<Tags::Alpha, Tags::Beta, Tags::Gamma>;
248 0 : using source_tags = tmpl::list<
249 : ::Tags::FixedSource<Tags::MMode>, Tags::SingularField,
250 : ::Tags::deriv<Tags::SingularField, tmpl::size_t<2>, Frame::Inertial>,
251 : Tags::BoyerLindquistRadius>;
252 :
253 : // Background
254 0 : tuples::tagged_tuple_from_typelist<background_tags> variables(
255 : const tnsr::I<DataVector, 2>& x, background_tags /*meta*/) const;
256 :
257 : // Initial guess
258 0 : static tuples::TaggedTuple<Tags::MMode> variables(
259 : const tnsr::I<DataVector, 2>& x, tmpl::list<Tags::MMode> /*meta*/);
260 :
261 : // Fixed sources
262 0 : tuples::tagged_tuple_from_typelist<source_tags> variables(
263 : const tnsr::I<DataVector, 2>& x, source_tags /*meta*/) const;
264 :
265 : template <typename... RequestedTags>
266 0 : tuples::TaggedTuple<RequestedTags...> variables(
267 : const tnsr::I<DataVector, 2>& x, const Mesh<2>& /*mesh*/,
268 : const InverseJacobian<DataVector, 2, Frame::ElementLogical,
269 : Frame::Inertial>& /*inv_jacobian*/,
270 : tmpl::list<RequestedTags...> /*meta*/) const {
271 : return variables(x, tmpl::list<RequestedTags...>{});
272 : }
273 :
274 : // NOLINTNEXTLINE
275 0 : void pup(PUP::er& p) override;
276 :
277 : private:
278 0 : friend bool operator==(const CircularOrbit& lhs, const CircularOrbit& rhs);
279 :
280 0 : double black_hole_mass_{std::numeric_limits<double>::signaling_NaN()};
281 0 : double black_hole_spin_{std::numeric_limits<double>::signaling_NaN()};
282 0 : double orbital_radius_{std::numeric_limits<double>::signaling_NaN()};
283 0 : int m_mode_number_{};
284 0 : std::optional<std::array<double, 4>> hyperboloidal_slicing_transitions_{};
285 0 : bool penetrating_horizon_{false};
286 0 : bool impose_equatorial_symmetry_{false};
287 : };
288 :
289 0 : bool operator!=(const CircularOrbit& lhs, const CircularOrbit& rhs);
290 :
291 : } // namespace ScalarSelfForce::AnalyticData
|
