This is the central location where all variables should be defined, and any relationships between them should be given. Having all definitions collected in one file is useful because other files can reference this one, so there is no need for duplication, and less room for mistakes. In particular, the relationships between variables are defined only here, so we don't need to reimplement those relationships. And if you do ever find a mistake, you only have to fix it in one place, then just re-run the other notebooks.
There are two main classes of variables:
The distinction is only required for code output, to ensure that everything gets calculated correctly. The PN equations you write down and manipulate can be in terms of any of these variables.
The fundamental variables that go into PN equations are things like the mass, spins $\chi_1$, and $\chi_2$, orbital angular velocity $\hat{\ell}$, and unit separation vector $\hat{n}$. We also include the tidal-coupling parameters in this list. Also, note that only $M_1$ is included. This is because the total mass is always assumed to be 1, so $M_2 = 1-M_1$.
The derived variables are further distinguished by whether they will need to be recalculated at each time step or not. For example, though we define the spins fundamentally as $\chi_1$ and $\chi_2$, we can also define derived spins $S$ and $\Sigma$, which need to be recalculated if the system is precessing. On the other hand, the masses are constant and fundamentally defined by $M_1$, so $M_2$ and $\nu$ only need to be calculated from that information once.
In [10]:
# Make sure division of integers does not round to the nearest integer
from __future__ import division
import sys
sys.path.insert(0, '..') # Look for modules in directory above this one
# Make everything in python's symbolic math package available
from sympy import * # Make sure sympy functions are used in preference to numpy
import sympy # Make sympy. constructions available
from sympy import Rational as frac # Rename for similarity to latex
from sympy import log as ln
# Print symbolic expressions nicely
init_printing()
# We'll use the numpy `array` object for vectors
from numpy import array, cross, dot
# We'll use a custom object to keep track of variables
from Utilities.PNObjects import PNCollection
PNVariables = PNCollection()
Only the most basic variables should be defined here.
Note that we will be using (quaternion) logarithmic rotors to describe the orientations of the spins, and the orientation and velocity of the binary itself. This allows us to reduce the number of constraints in the system, and only evolve the minimal number of equations. For example, the spins are constant, so only two degrees of freedom are needed. These can be expressed without ambiguities or singularities in the form of logarithmic rotors: $\mathfrak{r}_1 = \mathfrak{r}_{\chi_1 x} \hat{x} + \mathfrak{r}_{\chi_1 y} \hat{y}$, so that $\vec{\chi}_1 = S_{\chi_1}\, e^{\mathfrak{r}_1}\, \hat{z}\, e^{-\mathfrak{r}_1}\, \bar{S}_{\chi_1}$. Here, $S_{\chi_1}$ is a constant spinor with magnitude $\sqrt{\lvert \chi_1 \rvert}$ encoding the initial direction of the spin. This may look complicated, but it performs very well numerically.
We will still be able to write and manipulate the PN equations directly in terms of familiar quantities like $\vec{S}_1 \cdot \hat{\ell}$, etc., but the fundamental objects will be the rotors, which means that the substitutions made for code output will automatically be in terms of the rotors.
In [11]:
# Unit basis vectors
PNVariables.AddBasicConstants('xHat, yHat, zHat', datatype='Quaternions::Quaternion', commutative=False)
# Dimensionful quantities, just in case anybody uses them...
PNVariables.AddBasicConstants('G, c')
# Masses of objects 1 and 2.
PNVariables.AddBasicConstants('M1')
PNVariables.AddBasicConstants('M2')
# Angular speed of separation vector
PNVariables.AddBasicVariables('v')
# Initial spins expressed as spinors taking zHat onto those spins (assumed to have constant magnitudes)
PNVariables.AddBasicConstants('S_chi1', datatype='Quaternions::Quaternion', commutative=False)
PNVariables.AddBasicConstants('S_chi2', datatype='Quaternions::Quaternion', commutative=False)
# Dynamic spin directions
PNVariables.AddBasicVariables('rfrak_chi1_x, rfrak_chi1_y')
PNVariables.AddBasicVariables('rfrak_chi2_x, rfrak_chi2_y')
# Tidal deformabilities, in units where the total mass is 1
PNVariables.AddBasicConstants('lambda1, lambda2')
# Frame aligned to Orbital angular velocity vector and magnitude ("Newtonian" angular momentum)
PNVariables.AddBasicVariables('rfrak_frame_x, rfrak_frame_y, rfrak_frame_z')
Any variable that can be derived from the variables above should be put in this section.
These variables should probably be left in arbitrary form, unless a particular simplification is desired. The substitutions dictionary should map from the general names and their definitions in terms of basic variables. In numerical codes, their values can be calculated once per time step and then stored, so that the values do not have to be re-calculated every time they appear in an expression.
Various common combinations of the two masses:
In [12]:
PNVariables.AddDerivedConstant('M', M1+M2)
PNVariables.AddDerivedConstant('delta', (M1-M2)/M)
PNVariables.AddDerivedConstant('nu', M1*M2/M**2)
PNVariables.AddDerivedConstant('nu__2', (M1*M2/M**2)**2)
PNVariables.AddDerivedConstant('nu__3', (M1*M2/M**2)**3)
PNVariables.AddDerivedConstant('q', M1/M2)
The system's vector basis is given by $(\hat{\ell}, \hat{n}, \hat{\lambda})$, and will be computed by the code in terms of the fundamental logarithmic rotors defined above. Here, we give all the substitutions that will be needed in the code.
In [13]:
# This rotor encodes all information about the frame
PNVariables.AddDerivedVariable('R', exp(rfrak_frame_x*xHat + rfrak_frame_y*yHat + rfrak_frame_z*zHat),
datatype='Quaternions::Quaternion', commutative=False)
# Unit separation vector between the compact objects
PNVariables.AddDerivedVariable('nHat', R*xHat*conjugate(R), datatype='Quaternions::Quaternion')
# Unit vector orthogonal to the other two; in the direction of velocity
PNVariables.AddDerivedVariable('lambdaHat', R*yHat*conjugate(R), datatype='Quaternions::Quaternion')
# Unit vector in direction of angular velocity
PNVariables.AddDerivedVariable('ellHat', R*zHat*conjugate(R), datatype='Quaternions::Quaternion')
# Components of the above
for i,d in zip(['1','2','3'],['x','y','z']):
PNVariables.AddDerivedVariable('nHat_'+d, substitution_atoms=[nHat], substitution='nHat['+i+']')
for i,d in zip(['1','2','3'],['x','y','z']):
PNVariables.AddDerivedVariable('lambdaHat_'+d, substitution_atoms=[lambdaHat], substitution='lambdaHat['+i+']')
for i,d in zip(['1','2','3'],['x','y','z']):
PNVariables.AddDerivedVariable('ellHat_'+d, substitution_atoms=[ellHat], substitution='ellHat['+i+']')
Various spin components and combinations:
In [15]:
# These rotors encode all information about the spin directions
PNVariables.AddDerivedVariable('R_S1', exp(rfrak_chi1_x*xHat + rfrak_chi1_y*yHat),
datatype='Quaternions::Quaternion', commutative=False)
PNVariables.AddDerivedVariable('R_S2', exp(rfrak_chi2_x*xHat + rfrak_chi2_y*yHat),
datatype='Quaternions::Quaternion', commutative=False)
# The spins are derived from rfrak_chi1_x, etc.
PNVariables.AddDerivedVariable('chiVec1', S_chi1*R_S1*zHat*conjugate(R_S1)*conjugate(S_chi1), datatype='Quaternions::Quaternion')
PNVariables.AddDerivedVariable('chiVec2', S_chi2*R_S2*zHat*conjugate(R_S2)*conjugate(S_chi2), datatype='Quaternions::Quaternion')
for i,d in zip(['1','2','3'],['x','y','z']):
PNVariables.AddDerivedVariable('chi1_'+d, substitution_atoms=[chiVec1], substitution='chiVec1['+i+']')
for i,d in zip(['1','2','3'],['x','y','z']):
PNVariables.AddDerivedVariable('chi2_'+d, substitution_atoms=[chiVec2], substitution='chiVec2['+i+']')
PNVariables.AddDerivedConstant('chi1chi1', substitution_atoms=[chiVec1], substitution='chiVec1.normsquared()')
PNVariables.AddDerivedConstant('chi1chi2', substitution_atoms=[chiVec1,chiVec2], substitution='chiVec1.dot(chiVec2)')
PNVariables.AddDerivedConstant('chi2chi2', substitution_atoms=[chiVec2], substitution='chiVec2.normsquared()')
PNVariables.AddDerivedVariable('chi1_n', substitution_atoms=[chiVec1,nHat], substitution='chiVec1.dot(nHat)')
PNVariables.AddDerivedVariable('chi1_lambda', substitution_atoms=[chiVec1,lambdaHat], substitution='chiVec1.dot(lambdaHat)')
PNVariables.AddDerivedVariable('chi1_ell', substitution_atoms=[chiVec1,ellHat], substitution='chiVec1.dot(ellHat)')
PNVariables.AddDerivedVariable('chi2_n', substitution_atoms=[chiVec2,nHat], substitution='chiVec2.dot(nHat)')
PNVariables.AddDerivedVariable('chi2_lambda', substitution_atoms=[chiVec2,lambdaHat], substitution='chiVec2.dot(lambdaHat)')
PNVariables.AddDerivedVariable('chi2_ell', substitution_atoms=[chiVec2,ellHat], substitution='chiVec2.dot(ellHat)')
PNVariables.AddDerivedConstant('sqrt1Mchi1chi1', sqrt(1-chi1chi1))
PNVariables.AddDerivedConstant('sqrt1Mchi2chi2', sqrt(1-chi2chi2))
PNVariables.AddDerivedVariable('S', chiVec1*M1**2 + chiVec2*M2**2, datatype=chiVec1.datatype)
PNVariables.AddDerivedVariable('S_ell', chi1_ell*M1**2 + chi2_ell*M2**2)
PNVariables.AddDerivedVariable('S_n', chi1_n*M1**2 + chi2_n*M2**2)
PNVariables.AddDerivedVariable('S_lambda', chi1_lambda*M1**2 + chi2_lambda*M2**2)
PNVariables.AddDerivedVariable('Sigma', M*(chiVec2*M2 - chiVec1*M1), datatype=chiVec1.datatype)
PNVariables.AddDerivedVariable('Sigma_ell', M*(chi2_ell*M2 - chi1_ell*M1))
PNVariables.AddDerivedVariable('Sigma_n', M*(chi2_n*M2 - chi1_n*M1))
PNVariables.AddDerivedVariable('Sigma_lambda', M*(chi2_lambda*M2 - chi1_lambda*M1))
PNVariables.AddDerivedVariable('S1', chiVec1*M1**2, datatype=chiVec1.datatype)
PNVariables.AddDerivedVariable('S1_ell', chi1_ell*M1**2)
PNVariables.AddDerivedVariable('S1_n', chi1_n*M1**2)
PNVariables.AddDerivedVariable('S1_lambda', chi1_lambda*M1**2)
PNVariables.AddDerivedVariable('S2', chiVec2*M2**2, datatype=chiVec1.datatype)
PNVariables.AddDerivedVariable('S2_ell', chi2_ell*M2**2)
PNVariables.AddDerivedVariable('S2_n', chi2_n*M2**2)
PNVariables.AddDerivedVariable('S2_lambda', chi2_lambda*M2**2)
PNVariables.AddDerivedVariable('chi_s', (chiVec1 + chiVec2)/2, datatype=chiVec1.datatype)
PNVariables.AddDerivedVariable('chi_s_ell', (chi1_ell+chi2_ell)/2)
PNVariables.AddDerivedVariable('chi_s_n', (chi1_n+chi2_n)/2)
PNVariables.AddDerivedVariable('chi_s_lambda', (chi1_lambda+chi2_lambda)/2)
PNVariables.AddDerivedVariable('chi_a', (chiVec1 - chiVec2)/2, datatype=chiVec1.datatype)
PNVariables.AddDerivedVariable('chi_a_ell', (chi1_ell-chi2_ell)/2)
PNVariables.AddDerivedVariable('chi_a_n', (chi1_n-chi2_n)/2)
PNVariables.AddDerivedVariable('chi_a_lambda', (chi1_lambda-chi2_lambda)/2)
Other functions of the angular velocity that find frequent use:
In [6]:
PNVariables.AddDerivedVariable('x', v**2)
PNVariables.AddDerivedVariable('Omega_orb', v**3/M)
PNVariables.AddDerivedVariable('logv', log(v))
Note the simple expressions $M_{1,2} = M \frac{1 \pm \delta}{2}$, which lead to the useful translations \begin{align} \frac{\vec{S}}{M^2} &= (1-2\nu)\vec{\chi}_s + \delta \vec{\chi}_a, \\ \frac{\vec{\Sigma}}{M^2} &= -\delta\vec{\chi}_s - \vec{\chi}_a, \end{align} and equivalently \begin{align} \chi_s &= \frac{\vec{S} + \delta \vec{\Sigma}}{2M^2\nu}, \\ \chi_a &= - \frac{\delta \vec{S} + (1-2\nu) \vec{\Sigma}} {2M^2\nu}. \end{align}