From 3420bbde034371a7c416858cc037daad9c2a9ef1 Mon Sep 17 00:00:00 2001 From: Aniket Khairnar Date: Wed, 14 Jan 2026 10:54:48 -0600 Subject: [PATCH 1/3] Adding script for computing boosted com-charge from analytical expressions --- scri/pn/boosted_comcharge.py | 106 +++++++++++++++++++++++++++++++++++ 1 file changed, 106 insertions(+) create mode 100644 scri/pn/boosted_comcharge.py diff --git a/scri/pn/boosted_comcharge.py b/scri/pn/boosted_comcharge.py new file mode 100644 index 00000000..be2d3167 --- /dev/null +++ b/scri/pn/boosted_comcharge.py @@ -0,0 +1,106 @@ +import numpy as np + +def PN_charges(m, ν): + """ + Determine a tuple of callables for passing to + com_transformation_map_to_superrest_frame function for a specific input of + mass (m) and symmetric mass ratio (ν). + + Parameters + ---------- + m: float, real + Total mass of the binary. + ν: float, real + Symmetric mass ratio of the binary. + + Returns + ------- + (energy_pn, angular_momentum_pn, com_charge_pn, phase): tuple of callables + energy_pn is the PN series of energy up to 2PN order. + angular_momentum_pn is the PN series of angular momentum up to 3PN order. + com_charge_pn is the PN series of CoM charge up to leading order. + orbital_phase is the phase obtained from the (2,1) mode of the strain. + + Our conventions for defining the orbital phase differ from the standard + conventions used in PN theory. This stems from the fact that SpEC uses + h_{ab} to define the metric perturbation while PN theory uses h^{ab} for + the metric perturbation, which results in h^{NR}_{lm} = − h^{PN}_{lm}. + Also, the π/2 phase difference is due to the leading order complex phase + of h_{21} mode from PN theory. + + All of the callables accept abd object as a parameter. + """ + + def compute_x(abd): + """Helper function to extract PN parameter x.""" + h = abd.h + ω_mag = np.linalg.norm(h.angular_velocity(), axis=1) + x = (m * ω_mag) ** (2. / 3.) + + return x + + def energy_pn(abd): + x = compute_x(abd) + E = m - (m * ν * x)/2 * (1 + (-3/4 - ν/12)* x + (-27/8 + 19/8 * ν - ν**2/24) * x**2) + return E + + def angular_momentum_pn(abd): + x = compute_x(abd) + J_mag = (m**2 * ν * x**(-1/2)) * (1 + (3/2 + ν/6)* x + (-27/8 + 19/8 * ν - ν**2/24) * x**2 + + (135/16 + (-6889/144 + 41/24 * np.pi**2)* ν + 31/24 * ν**2 + 7/1296 * ν**3) * x**3) + return J_mag + + def G_mag_pn(abd): + x = compute_x(abd) + G_mag = -(1142/105) * x**(5/2) * m**2 * np.sqrt(1 - 4*ν) * ν**2 + return G_mag + + def orbital_phase(abd): + h = abd.h + h_21 = h.data[:, h.index(2,1)] + ψ = (-1)* np.unwrap(np.angle(-h_21) - np.pi/2) + return ψ + + return energy_pn, angular_momentum_pn, G_mag_pn, orbital_phase + + +def analytical_CoM_func(θ, t, E, J_mag, G_mag, ψ): + """ + Computes the boosted center-of-mass-charge using PN expressions of other charges. + + Parameters + ---------- + θ: list[floats] + Parameters of the model. + θ[0:3] are the components of boost_velocity, θ[3:6] are components of + the spatial translation (l=1 supertranslation), and there can be any + extra fit parameters desired (but their fitted values are ignored). + t: ndarray, real + Time array corresponding to the size of the center-of-mass charge. + E: ndarray, real + Array corresponding to the energy. + J_mag: ndarray, real + Array corresponding to the magnitude of angular momentum. + G_mag: ndarray, real + Array corresponding to the magnitude of the center-of-mass charge. + ψ: ndarray, real + Orbital phase obtained from the (2,1) mode of the strain. + + Returns + ------- + G: ndarray, real, shape(..., 3) + Model timeseries of the boosted center-of-mass charge. + """ + β = θ[0:3][None] + X = θ[3:6][None] + α1, α2 = θ[6:] + + # Get the orbital direction vectors, and the angular momentum vector. + cosψ, sinψ = np.cos(ψ), np.sin(ψ) + nvec = np.stack((cosψ, sinψ, 0*ψ),axis=1) + λvec = np.stack((-sinψ, cosψ,0*ψ),axis=1) + J = np.array([0*t,0*t,J_mag]).T + + G = α1 * (G_mag/E)[:, None] * λvec + α2 * (G_mag/E)[:, None] * nvec - np.cross( β, J/E[:,None]) - t[:,None] @ β + X + + return G \ No newline at end of file From 1c4f35e3aa7baa97678bdb115b03876d6f394ef4 Mon Sep 17 00:00:00 2001 From: Aniket Khairnar Date: Wed, 4 Feb 2026 14:12:09 -0600 Subject: [PATCH 2/3] Update to PN_charges docstring Co-authored-by: Mike Boyle --- scri/pn/boosted_comcharge.py | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/scri/pn/boosted_comcharge.py b/scri/pn/boosted_comcharge.py index be2d3167..4d0fd1a2 100644 --- a/scri/pn/boosted_comcharge.py +++ b/scri/pn/boosted_comcharge.py @@ -15,20 +15,26 @@ def PN_charges(m, ν): Returns ------- - (energy_pn, angular_momentum_pn, com_charge_pn, phase): tuple of callables - energy_pn is the PN series of energy up to 2PN order. - angular_momentum_pn is the PN series of angular momentum up to 3PN order. - com_charge_pn is the PN series of CoM charge up to leading order. - orbital_phase is the phase obtained from the (2,1) mode of the strain. - - Our conventions for defining the orbital phase differ from the standard - conventions used in PN theory. This stems from the fact that SpEC uses - h_{ab} to define the metric perturbation while PN theory uses h^{ab} for - the metric perturbation, which results in h^{NR}_{lm} = − h^{PN}_{lm}. - Also, the π/2 phase difference is due to the leading order complex phase - of h_{21} mode from PN theory. + energy_pn : callable + Returns PN approximation for energy up to 2PN order + angular_momentum_pn : callable + Returns PN approximation for angular momentum up to 3PN order + com_charge_pn : callable + Returns PN approximation for CoM charge up to leading order + phase : callable + Returns the phase obtained from the (2,1) mode of the strain + + All of the returned callables accept an `ABD` object as a parameter. + + Notes + ------ + Our conventions for defining the orbital phase differ from the standard + conventions used in PN theory. This stems from the fact that SpEC uses + h_{ab} to define the metric perturbation while PN theory uses h^{ab} for + the metric perturbation, which results in h^{NR}_{l,m} = − h^{PN}_{l,m}. + Also, the π/2 phase difference is due to the leading-order complex phase + of h_{2,1} mode from PN theory. - All of the callables accept abd object as a parameter. """ def compute_x(abd): From 2a2bb913e18170d7495ad87f4a07bc154ffac662 Mon Sep 17 00:00:00 2001 From: Aniket Khairnar Date: Wed, 4 Feb 2026 14:49:38 -0600 Subject: [PATCH 3/3] Incorporating code feedback from Mike --- scri/pn/boosted_comcharge.py | 104 +++++++++++++++++++++++------------ 1 file changed, 68 insertions(+), 36 deletions(-) diff --git a/scri/pn/boosted_comcharge.py b/scri/pn/boosted_comcharge.py index 4d0fd1a2..79a5f30d 100644 --- a/scri/pn/boosted_comcharge.py +++ b/scri/pn/boosted_comcharge.py @@ -1,11 +1,11 @@ import numpy as np + def PN_charges(m, ν): - """ - Determine a tuple of callables for passing to - com_transformation_map_to_superrest_frame function for a specific input of - mass (m) and symmetric mass ratio (ν). - + """ + Determine a tuple of callables for a specific input of mass (m) and + symmetric mass ratio (ν). + Parameters ---------- m: float, real @@ -24,7 +24,17 @@ def PN_charges(m, ν): phase : callable Returns the phase obtained from the (2,1) mode of the strain - All of the returned callables accept an `ABD` object as a parameter. + All returned callables accept an `ABD` object as an input parameter and + return their respective quantities evaluated at the same time steps as the + `ABD` object. The PN expressions for energy and angular momentum are taken + from Eqs. (337) and (338) of Blanchet's Living Review (2014) + . The PN expression for CoM charge is + derived in Khairnar et al. . The tuple of these callables + can be passed as the `Gargsfun` keyword argument to + `map_to_superrest_frame`. These callables are then used to determine the + Gargs parameters to `transformation_from_CoM_charge`. They are called as + Gargs = [func(abd) for func in Gargsfun] if Gargsfun else None + within the `com_transformation_to_map_to_superrest_frame` function. Notes ------ @@ -41,30 +51,33 @@ def compute_x(abd): """Helper function to extract PN parameter x.""" h = abd.h ω_mag = np.linalg.norm(h.angular_velocity(), axis=1) - x = (m * ω_mag) ** (2. / 3.) - + x = (m * ω_mag) ** (2.0 / 3.0) return x def energy_pn(abd): x = compute_x(abd) - E = m - (m * ν * x)/2 * (1 + (-3/4 - ν/12)* x + (-27/8 + 19/8 * ν - ν**2/24) * x**2) + E = m - (m * ν * x) / 2 * (1 + (-3 / 4 - ν / 12) * x + (-27 / 8 + 19 / 8 * ν - ν**2 / 24) * x**2) return E def angular_momentum_pn(abd): x = compute_x(abd) - J_mag = (m**2 * ν * x**(-1/2)) * (1 + (3/2 + ν/6)* x + (-27/8 + 19/8 * ν - ν**2/24) * x**2 - + (135/16 + (-6889/144 + 41/24 * np.pi**2)* ν + 31/24 * ν**2 + 7/1296 * ν**3) * x**3) + J_mag = (m**2 * ν * x ** (-1 / 2)) * ( + 1 + + (3 / 2 + ν / 6) * x + + (-27 / 8 + 19 / 8 * ν - ν**2 / 24) * x**2 + + (135 / 16 + (-6889 / 144 + 41 / 24 * np.pi**2) * ν + 31 / 24 * ν**2 + 7 / 1296 * ν**3) * x**3 + ) return J_mag def G_mag_pn(abd): x = compute_x(abd) - G_mag = -(1142/105) * x**(5/2) * m**2 * np.sqrt(1 - 4*ν) * ν**2 + G_mag = -(1142 / 105) * x ** (5 / 2) * m**2 * np.sqrt(1 - 4 * ν) * ν**2 return G_mag def orbital_phase(abd): h = abd.h - h_21 = h.data[:, h.index(2,1)] - ψ = (-1)* np.unwrap(np.angle(-h_21) - np.pi/2) + h_21 = h.data[:, h.index(2, 1)] + ψ = (-1) * np.unwrap(np.angle(-h_21) - np.pi / 2) return ψ return energy_pn, angular_momentum_pn, G_mag_pn, orbital_phase @@ -72,41 +85,60 @@ def orbital_phase(abd): def analytical_CoM_func(θ, t, E, J_mag, G_mag, ψ): """ - Computes the boosted center-of-mass-charge using PN expressions of other charges. - + Compute a model time series for boosted center-of-mass (CoM) charge using PN + expressions of energy, angular momentum, and CoM charge, along with the + phase computed from the h_{21} mode. + + This model timeseries is derived in Khairnar et al. . + It serves as a fitting function that can be passed as the `Gfun` keyword + argument to the `map_to_superrest_frame` function. All the arguments are + computed over the window used for fixing the frame. + Parameters ---------- - θ: list[floats] + θ: ndarray, real, shape(8,) Parameters of the model. - θ[0:3] are the components of boost_velocity, θ[3:6] are components of - the spatial translation (l=1 supertranslation), and there can be any - extra fit parameters desired (but their fitted values are ignored). + - θ[0:3] : components of the boost velocity + - θ[3:6] : components of the spatial translation + - θ[6:] : two additional fit parameters referred to as the nuisance + parameters in Khairnar et al. ; their + fitted values are ignored. t: ndarray, real - Time array corresponding to the size of the center-of-mass charge. + Time array corresponding to the window over which the frame fixing is + performed. E: ndarray, real - Array corresponding to the energy. + PN approximation for the energy computed over the fitting window. J_mag: ndarray, real - Array corresponding to the magnitude of angular momentum. + PN approximation for the magnitude of angular momentum computed over the + fitting window. G_mag: ndarray, real - Array corresponding to the magnitude of the center-of-mass charge. + PN approximation for the magnitude of the CoM charge computed over the + fitting window. ψ: ndarray, real - Orbital phase obtained from the (2,1) mode of the strain. + Unwrapped orbital phase obtained from the (2,1) mode of the strain over + the fitting window. Returns ------- G: ndarray, real, shape(..., 3) - Model timeseries of the boosted center-of-mass charge. + Model time series of the boosted center-of-mass charge. """ - β = θ[0:3][None] - X = θ[3:6][None] + β = θ[0:3][np.newaxis] + X = θ[3:6][np.newaxis] α1, α2 = θ[6:] - + # Get the orbital direction vectors, and the angular momentum vector. cosψ, sinψ = np.cos(ψ), np.sin(ψ) - nvec = np.stack((cosψ, sinψ, 0*ψ),axis=1) - λvec = np.stack((-sinψ, cosψ,0*ψ),axis=1) - J = np.array([0*t,0*t,J_mag]).T - - G = α1 * (G_mag/E)[:, None] * λvec + α2 * (G_mag/E)[:, None] * nvec - np.cross( β, J/E[:,None]) - t[:,None] @ β + X - - return G \ No newline at end of file + nvec = np.array([cosψ, sinψ, 0 * ψ]).T + λvec = np.array([-sinψ, cosψ, 0 * ψ]).T + J = np.array([0 * t, 0 * t, J_mag]).T + + G = ( + α1 * (G_mag / E)[:, np.newaxis] * λvec + + α2 * (G_mag / E)[:, np.newaxis] * nvec + - np.cross(β, J / E[:, np.newaxis]) + - (t[:, np.newaxis] @ β) + + X + ) + + return G