Axion-Photon Conversion Signals from Neutron Stars with Spacetime Curvature Accounted for in the Magnetosphere Model
Jesse Satherley Chris Gordon Chris Stevens


APS/123-QED



 jesse.satherley@pg.canterbury.ac.nz

 chris.gordon@canterbury.ac.nz
 
 School of Physical and Chemical Sciences, University of Canterbury, Christchurch, New Zealand
 

 chris.stevens@canterbury.ac.nz
 
 School of Mathematics and Statistics, University of Canterbury, Christchurch, New Zealand
 







§ ABSTRACT 
 Axions are a well-motivated dark matter candidate. They may be detectable from radio line emission from their resonant conversion in neutron star magnetospheres. While radio data collection for this signal has begun, further efforts are required to solidify the theoretical predictions for the resulting radio lines. Usually, the flat spacetime Goldreich-Julian model of the neutron star magnetosphere is used, while a Schwarzschild geometry is assumed for the ray tracing. 
We assess the impact of incorporating the spacetime curvature into the magnetosphere model.
We examine a range of neutron star and axion masses and 
find an average difference of 26% in radiated power compared to the standard Goldreich-Julian magnetosphere model for a 10 mass axion and a 2.2M_⊙ mass neutron star. 
A much lesser difference is found for lower-mass neutron stars, as in that case, axion-photon conversion occurs further from the Schwarzschild radius. 


 
Axion-Photon Conversion Signals from Neutron Stars with Spacetime Curvature Accounted for in the Magnetosphere Model
    C. Stevens
    September 9, 2024
====================================================================================================================





§  INTRODUCTION





The axion was introduced 
to solve the strong CP problem in quantum chromodynamics
<cit.>. 
It was subsequently realized that it could account for a fraction, or all, of the
dark matter in the Universe <cit.>.

A recent and compelling proposal for indirectly searching for axions in astrophysical environments involves detecting radio photons produced by axion-photon mixing in neutron star magnetospheres. The strong magnetic fields and ambient plasma in these regions can resonantly amplify the mixing process, potentially revealing axions through distinctive radio signatures
<cit.>. 
Two main schemes exist to determine the signal received from the axion-photon conversions: (1) the emitter-to-observer scheme samples the conversion surface 
 and propagates photons forward through the plasma <cit.>; and
(2) the observer-to-emitter scheme sources photons at a distant detector and propagates them backwards onto the conversion surface <cit.>. 
Both methods have their own benefits and drawbacks. 

The early simulations of the axion-photon conversion process assumed a flat spacetime, did not consider the refractive effects of the plasma, and employed a simple flat-spacetime Goldreich-Julian (GJ) model <cit.> for the NS magnetosphere, e.g., <cit.>. 
More recently, work was done to improve this by including the dispersive effects of an isotropic (unmagnetized) plasma in Schwarzschild spacetime <cit.>. However, this still included the GJ model.
Other studies attempt to account for the random infall of axions onto the NS via a Monte Carlo style simulation <cit.>.
Although they include a magnetized plasma dispersion relationship for the photons, they only consider flat spacetime.
Most recently, work was done to include a magnetized plasma dispersion relationship for the photons in Schwarzschild spacetime <cit.>. 
This article also had the inclusion of using a recently derived axion-photon conversion probability incorporating 3D effects <cit.>.
However, all these models rely on the inclined GJ model of NS magnetospheres.

The GJ magnetosphere model assumes a dipole magnetic field 
in a flat spacetime. 
In the force-free NS model case, GJ shows that an NS must have a dense magnetosphere containing charged particles.
They derive an analytical form for the charge density and magnetic field vector.
This proves convenient for studying phenomena around NSs due to the simplicity of the model, which allows quick and approximate results to be shown.
More recently, strides have been made to create increasingly accurate models for NS magnetospheres. Numerical simulations can account for the complete Maxwell equations in the 3+1 formalism of a stationary background metric <cit.>.
Extensions have been done to include more multipolar components of the magnetic fields in a vacuum with analytical equations <cit.>.
A series of papers by Gralla et al. provide an analytical solution in the near field of an inclined dipole magnetic field in curved spacetime <cit.>. 
Their model[Hereafter referred to as the Gralla, Lupsasca, and Philippov (GLP) model.] provides corrections to the dipole magnetic field in GR, which alters the magnetic field strength and shape. 
Their article also includes a charge distribution constituting a plasma in the magnetosphere that depends on GR effects.

While our principal methods follow the work carried out in Refs. <cit.>, we wish to extend their results by using the GLP model and seeing what effect that has on the predicted power received from axion-photon conversions near an NS. 
We will do this by producing estimated signals in the emitter to observer scheme simulations.

In Sec. <ref>, we provide a review of the GJ model, the GLP model, and a brief aside on plasma.
In Sec. <ref>, we discuss our implementation of the conversion of axions to photons in our simulations.
In Sec. <ref>, the dispersion relationships are provided, which are then used in the ray tracing equations reviewed in Sec. <ref>.
In Sec. <ref>, we explain our implementation of axion-to-photon conversion simulations around neutron stars using the observer-to-emitter scheme.
Lastly, in Sec. <ref>, we present the results of our simulations and the conclusions that can be made. The article's main text considers the case where the magnetic field and rotation axes are aligned. We consider the misaligned case in the Appendices.

Throughout, we denote 4-dimensional abstract (where no particular coordinate system is specified) and coordinate tensor indices with Latin letters starting from a,b,c,… and i,j,k,… respectively, both in the range 0–3 and 3-dimensional coordinate indices with Greek letters μ,ν,… in the range 1–3. The 3-vectors with index range 1,2,3 will also be denoted by a boldface typeset where appropriate. We use the metric signature (-,+,+,+), along with the choice of natural units c=ħ=ϵ_0=1.



§  NEUTRON STAR MODELS


As already highlighted in the introduction, the leading NS magnetosphere model that is used in axion-photon ray tracing simulations is the GJ model, which assumes an inclined dipole magnetic field in a flatspace time. 
However, we wish to explore the effect on these simulations of including an NS model derived in a curved spacetime. 
This section will review and discuss the two models, focusing on their implementation into the simulations.





 §.§ Goldreich-Julian Neutron Star


The Goldreich-Julian (GJ) <cit.>
 assumes that the NS is surrounded by a dense plasma of charged particles in the star's magnetosphere. 
GJ derives forms for the magnetic and electric fields for an aligned rotator. The GJ model has since been extended to account for the misalignment angle <cit.>. 

The number density of electrons and positrons in the magnetosphere of the GJ model is given by,

    n_GJ(r)=2Ω·B/|q_e|1/1-Ω^2r^2sin^2θ,

where Ω=(2π/P_ns)z is the NSs angular velocity vector with P_ns being the period of rotation, z the unit vector in line with the rotation axis of the star, q_e the charge of an electron, r  the distance from the center of the NS, and θ the polar angle from the rotation axis. 
The magnetic field is B, which will be defined shortly for the GJ magnetic field (<ref>). 
The regions that return a positive number density are dominated by positrons, whereas regions with a negative number density are dominated by electrons. 

The GJ model described above assumes that the star's interior magnetic field takes the form of a dipole. From this assumption, expressions for the fields can be derived. 
Since we wish to work with the near-zone fields, both of the electromagnetic fields simplify. 
This is done by taking the leading order terms as r→0 as long as r>R_NS so that the fields are external to the star, where R_NS is the radius of the NS's surface.

In the near-zone, the magnetic fields to a good approximation are given by that of an idealized inclined rotator (e.g. <cit.>),

    B_r    = 2μ/r^3(cosχcosθ + sinχsinθcosλ),
    
            B_θ   = μ/r^3(cosχsinθ - sinχcosθcosλ),
    
            B_ϕ   = μ/r^3sinχsinλ,

with μ being the magnetic dipole moment of the star, χ the misalignment angle between the rotation axis and the magnetic field axis, and λ=ϕ-Ω t, with ϕ the azimuthal angle around the NS where ϕ=0 is inline with x and increases in the anticlockwise direction.
These coordinates are shown in Fig. <ref>.



Because the GJ model is the simplest NS model which includes a plasma, with a well-defined charge density, it is the usual choice for axion-photon conversion papers, for example <cit.>. 
However, we wish to study the effect of including a dipole magnetic field derived in curved spacetime as compared to the GJ model.
In the following section we detail necessary equations from the GLP model.



 §.§ Gralla, Lupsasca and Philippov Magnetosphere


Recently, derivations have been published that extend the dipole magnetic field around an NS to a curved spacetime.
The GR pulsar model we explore in this work is taken from a series of papers by Gralla et al. <cit.>.
In their work, they use the Hartle-Thorne metric (<ref>) to describe the electromagnetic fields around a pulsar.
The Hartle-Thorne metric is used to describe a slowly rotating star in general relativity.
It is given as (e.g., <cit.>)

    g_ij = 
        [ -f^2+Ω^2r^2sin^2θ                 0                 0      -Ω r^2sin^2θ;                 0             1/f^2                 0                 0;                 0                 0               r^2                 0;      -Ω r^2sin^2θ                 0                 0         r^2sin^2θ;                   ],

 where f = 1 - r_s/r is the Schwarzschild function
and we are using Schwarzschild coordinates x^i = {t,r,θ,ϕ}. 
The rotation is contained with the terms including Ω, which is the angular velocity of the star. 
When the rotation is slow enough, frame-dragging terms that contain Ω become insignificant.
The Hartle-Thorne metric then simplifies to the Schwarzschild metric. 

In our work, we continue using the Schwarzschild metric for both the magnetic field and the ray tracing[We checked that the Hartle-Thorne metric does not make a significant difference to our results. However, it may matter for faster rotating NSs.].
Gralla et al. provide solutions for the near fields in the situation of a force-free axisymmetric field, after which they extend to an inclined pulsar where the misalignment angle is non-zero. 
In the following, we detail their work, and we provide insights on how to apply their results to the problem of axion photon conversion near NSs.
We begin by describing a set of electromagnetic relationships.

We can define the magnetic field for an arbitrary observer and metric by using (e.g. <cit.>),

    B^d = 1/2ϵ^abcdF_abU_c

where F_ab is the electromagnetic tensor, U_c is the 4-velocity of an observer and ϵ^abcd the covariant Levi-Civita tensor to account for the metric, which is related to the Levi-Civita symbol multiplied by the determinant of the metric tensor such that (e.g. <cit.>),

    ϵ_abcd   = (√(|g|))ε_abcd,
    ϵ^abcd   = sign(g)/(√(|g|))ε^abcd,

where ε_abcd is the fourth-rank Levi-Civita symbol.
The last relationship we require is the current density 4-vector, as the charge density is contained within the current density 4-vector. It is related to the 3-space quantities by the relationship,

    J^a=(ρ, J),

where ρ is the charge density and J=ρu is the current density  with velocity u^μ=dx^μ/dτ, where τ is the proper time. The current density 4-vector is related to the electromagnetic tensor by,

    J^a = ∇_b F^ab,

where ∇_b is the covariant derivative. 



  §.§.§ Aligned Rotator

The first paper in the series by Gralla et al. <cit.> begins by deriving an analytical method for studying the force-free magnetosphere of a slowly rotating aligned rotator (χ = 0), including the effects of GR. 
They find that a ∼60% correction to the dipole component of the surface magnetic field is introduced by accounting for GR. 

In deriving the equations for the magnetic field, they assume that the electromagnetic field is force-free (F_abJ^b=0).
This is the same assumption as the GJ  model.
The electromagnetic tensor can then be given by potentials ψ_i=ψ_i(r, θ, ϕ - Ω t), such that,

    F_ij =  ∂_iψ_1∂_jψ_2 - ∂_jψ_1∂_iψ_2.

As the field is axisymmetric, the time dependence does not alter the field configuration. 
However, we include it for completeness as it is required when considering the inclined magnetic field case. 
In the case of an aligned rotator, the magnetic flux function for a dipole in the near-field region is given as,

    
    ψ_1, near(r, θ)    = μ R^>_1(r)sin^2θ,
    ψ_2(ϕ-Ω t)    = ϕ-Ω t,


where μ=B_nsR_ ns^3/(2Δ_1) is the dipole moment with 
B_ns being the surface magnetic field strength at the pole,[This is related to the magnetic moment of (<ref>) via B_ns=μ.] R^>_1(r) is the first radial harmonic, and Δ_ℓ=R_ ns^ℓ R^>_ℓ(R_ ns), all with R_ ns the 
radius of the star. 
The function Δ_ℓ is dimensionless and depends only on the compactness of the star. 
It provides the GR correction to the dipole moment at the surface.
The first radial harmonic is (e.g. <cit.>),

    R^>_1(r)=-3/2r[3 - 4f + f^2 + 4log f/(1-f)^3],

recalling that
f=1-r_s/r is the Schwarzschild function. 

For this model to be beneficial, it must have a magnetosphere containing charged particles.
The charge density around the star can be described by investigating the charge-current 4-vector.
Upon applying the covariant derivative to the electromagnetic tensor with the Hartle Thorne metric (<ref>), to leading order derivative terms, the time component of J^a is[We have added back in factors of G compared to GLP's derivation due to our choice of units.],

    J^t = 2(Ω - Ω_z)/r(r - 2 G M)[(r - 3 G M)∂_rψ_1 + θ∂_θψ_1],

which with,

    ρ_e=U_aJ^a=J^t√(1-r_s/r)=J^t̂,

gives the charge density around the NS for the aligned rotator case.
This relationship is used to find the charge density for the plasma frequency, which will be discussed shortly.
We can then compare the results of the fields and charge density described here to the GJ model (<ref>) and (<ref>)  to see the effect these corrections have on the axion-photon conversion signal.

In the simulations presented in this paper, we take the magnetic field and rotation axis to be aligned so that we may compare with Ref. <cit.> (henceforth referred to as mcdonaldGeneralizedRayTracing2023MW23).
So, for the results presented in this paper, it is only necessary to understand the GLP-aligned rotator reviewed in the main body.
However, our code implements the inclined rotator case of the GLP model.
In Appendix <ref>, we detail the equations and relationships necessary to include the inclined GLP model.



 §.§ Plasma


A fundamental parameter that characterizes a plasma is its plasma frequency. In the absence of a magnetic field, the plasma frequency is the oscillation frequency for the charge distribution about its equilibrium and is given as (e.g. <cit.>)

    ω_p=√(∑_i4πα n_i/m_i),

where m_i and n_i is the mass and number density of species i, and α is the fine-structure constant.

It is assumed that the plasma consists only of electrons and positrons.[If ions were considered for the positively charged regions instead, the plasma frequency would decrease due to the ion's greater mass.] 
This type of plasma is used, for example, by Refs. <cit.>. 
The number density of the charged particles and the charge density are just  related via,

    n_e = |ρ_e|/q_e,

where q_e is the charge of an electron.

In previous studies three cases are considered for the plasma:

    
  * The plasma consists of no charged particles in the magnetosphere, considering only a vacuum around the star. 
    
  * An isotropic plasma, where charged particles are present, but the magnetic field of the NS has no effect on the medium. In this case, the plasma frequency only depends on the distance from the star.
    
  * An anisotropic plasma, where the magnetic field induces new effects on the medium. 

The choice of plasma is important for the ray tracing of photons as it alters their dispersion.



§ AXION-PHOTON CONVERSION
 

As previously indicated, axions undergo resonant conversion to photons in the presence of plasma (e.g. <cit.>).
The benefit of this is an enhancement to the photon signal from axion-photon conversions, potentially leading to detectable signals from Earth.
The resonance occurs due to plasma effects generating an effective mass for the photon, allowing the axion and photon dispersion relations to intersect. 
This resonant condition is maximized when the plasma frequency is close to or equal to the mass of the axion[In SI units, this condition is ω_p=ħ m_a.] 

    ω_p ≈ m_a  .

mcdonaldGeneralizedRayTracing2023MW23 uses this condition for their isotropic plasma cases.
For an anisotropic plasma, they advocate using full kinematic matching of axion to photon conversion, which is given by k_a = k_γ. 





 §.§ Conversion Surface
 

Using the resonance condition in (<ref>),
we can define a surface in three-dimensional space surrounding a magnetic field source with a plasma.
This surface is referred to as the conversion surface.
The conversion surface is the region that will have the most predominant axion-photon conversion flux.
When considering the GJ dipole magnetic field, where B∝1/r^3, the charge number density (<ref>) can have its r dependency explicitly shown as,

    n(r, θ, ϕ, t) ≈1/r^3n(r=1 eV^-1, θ, ϕ, t),

where we have ignored the rightmost fraction in (<ref>), as it is near unity for small radii.
The expression above can then be combined with (<ref>) and (<ref>)
to give the radius of the conversion surface for a given θ, ϕ and t:

    r = √(ω_p^2(r=1 eV^-1, θ, ϕ, t)/m_a^2).

For the GLP model, where the dependency on r is more complicated, root-solving methods are required to find the conversion surface. 



 §.§ Probability of Conversion


As the process of an axion converting to a photon is based on classical field theory and is mediated by the interaction term in the axion's Lagrangian,  
an associated probability of conversion can be found (e.g. <cit.>).
This probability will affect the total radiated power predicted in simulations.
Hence, the choice of conversion probability method that is used has a significant impact on the results.
For this reason, we chose to use the conversion probability for an isotropic plasma that was given in Eq. (69) of mcdonaldGeneralizedRayTracing2023MW23 so that we may compare results and it conveniently already incorporates a curved spacetime.
This conversion probability is expressed as,

    P_aγγ = π g^2_aγγ|B|^2sin^2θ̃E_γ/|k^i∂_i(ω_p^2)|,

where g_aγγ is the axion-photon coupling constant, |B|^2=B_μ B^μ, θ̃ is the angle between the axion's momentum and the magnetic field, and E_γ is the energy of the photon at the point of conversion. 



§ CURVED SPACETIME DISPERSION RELATIONS


The warping of spacetime from the mass of an NS can be extreme near the star, significantly affecting the path that particles and light would follow compared to a flat spacetime.
By accounting for 
these 
influences on geodesics around the star, the results of axion-photon conversion simulations will be changed (see, for example, Refs. <cit.>). 

For the curvature of spacetime to be accounted for in the dispersion relationships above, the metric must be included in some form.
A simple method of converting the flat spacetime relationships to a curved spacetime is by simply taking the 3+1 approach where squared parameters can be converted to Einstein sums, as done in Ref. <cit.>. 
We can, however, extend this further by introducing covariant relationships for the photon 4-vector components.
In the following, we discuss how gravitational effects are accounted for in the photon dispersion relation inside different types of plasma. 

Firstly, the refractive index of the medium takes on a covariant form and becomes <cit.>,

    n^2= 1 + k_ak^a/(k_bU^b)^2,

where k^a is the photons 4-momentum[Remember that in natural units, the 4-momentum and 4-wavevector are equivalent. In SI, p^a=ħ k^a.] and U^a is a global unit timelike vector such that U^aU_a=-1, recalling that U^a is the 4-velocity of an observer.
When the Minkowski metric is used, (<ref>) reduces to the flat-space case of n=ω/k where ω is the photon frequency and k is the photon wavenumber.



 §.§ Vacuum

For the simple case when no plasma is present, only gravity will affect the path the photon travels. 
The dispersion relationship in a vacuum is just, 
    
    D(k)=k_ak^a=0,

where the sum over the indices contains the metric, which will account for the curvature of space.
With the Minkowski metric, this simply becomes k^2-ω^2. 
In this case, photons travel along geodesics around the star.



 §.§ Isotropic Plasma

With the inclusion of plasma, there is now a scalar function ω_p(x^a) present.
It can can essentially be considered a forcing term in the dispersion relationship, which alters the trajectory of the photons from the vacuum case. 
The function ω_p is a scalar independent of the metric. 
Hence, it remains the same as the flat-space case, as it is unaffected by the inclusion of GR. 
In an isotropic plasma, the dispersion relationship is,

    D(k)=k_ak^a+ω_p=0.

For the derivation of this relationship see Eqn. (10) in Ref. <cit.>, which also matches with  Eqn. (17) of Ref. <cit.>.



 §.§ Anisotropic  Plasma
 
The covariant form of the anisotropic dispersion relationship requires some covariant plasma expressions presented in Ref. <cit.>[They use a metric with the opposite signature (+, -, -, -). We modify their relationships to be compatible with our choice of metric signature, which is (-, +, +, +). This leads to a difference in signs on some terms.].
As in previous sections, let U^a be a global unit timelike vector such that U^aU_a=-1.
Then, the photon's effective energy measured by an observer with 4-velocity U^a is W=-k_aU^a.
Also, define the unit vector in the direction of the magnetic field b^a = B^a/√(B^cB_c) with b^ab_a=1 (where in Schwarzschild coordinates B^0=B^t=0). 
Then, the wave vector component parallel to the magnetic field can be represented by the sum K_=k_ab^a. 
In the non-relativistic plasma limit of Eqn. (12) of Ref. <cit.>, we have that the GR anisotropic dispersion relation is,

    D(k)=k_ak^a+ω_p^2(1 - K_^2/W^2)=0.

When the appropriate Minkowski limit is taken, where K_=kcosθ̃ and W=ω, the above equation simplifies to the flat-space case.



§ RAY TRACING
 
When a photon propagates through a plasma, it may undergo refraction and reflection due to the plasma's varying refractive index. To trace the path of the photons through the plasma, ray tracing is used. 
At a simple level, this involves a system of coupled ordinary differential equations (ODEs), which are constructed using one of the plasma dispersion relationships in Sec. <ref>. 
The ODEs `tell' the photon which direction to travel, how its momentum should change direction, and how its energy should evolve. 
The ODEs can then be solved analytically or integrated using a numerical solver, depending on the complexity of the dispersion relation that forms the ODEs. 
The solutions allow the photon to be followed through the magnetosphere of the NS and can be used to reproduce the expected photon signal from axion-photon conversions.
When the dispersion relation contains a plasma frequency term, the ray path is most affected when the photon's frequency is close to the plasma frequency. 




 §.§ Ray Tracing in General Relativity
 

Paths of geodesics in flat space-time are simply straight lines. In a curved spacetime additional terms involving the Christoffel symbols appear in the geodesic equation, altering the paths <cit.>. In a curved spacetime with the presence of a plasma, the paths will change further <cit.>. Below we give the equations governing null geodesics in a curved spacetime with plasma.



  §.§.§ Geometric Optics 

The geometric optics limit is formed by taking the WKB approximation with the eikonal form <cit.>. The general relativistic equations for ray propagation can then be found by first representing a wave packet in the form,

    A_b=∫A_b(k)e^ik_bx^a√(|g|)d^4k

having used the eikonal form with a Fourier transform. In the exponential, k_b should satisfy D(k,x)=0, and in particular, there should be a sharp maximum when k=k_0. So, we can make the substitution k=q+k_0. Then, along the ray x^a=x^a(λ), where λ is an affine parameter along the light path, the phase q_ax^a should be stationary. So, we will have q_a(dx^a/dλ)=0. One can then expand the dispersion relation about k_0 to obtain q_a(∂ D/∂ k_a)=0. This yields,

    dx^a/dλ = ∂ D/∂ k_a.

Then using dD/dλ=0,

    ∂ D/∂ k_adk_a/dλ + ∂ D/∂ x_adx^a/dλ = 0.

Hence, we will have a system of ODEs that describes the ray path. They are expressed as,

    
    dx^a/dλ   = ∂ D/∂ k_a, 
    dk_a/dλ   = -∂ D/∂ x^a.


The solution to the first equation (<ref>) prescribes the spacetime position, whilst the second equation (<ref>) controls the energy and direction of propagation.
The convenience of this form is once we have defined the dispersion relation in covariant form, we may then directly use it to find the ray paths. The only caveat is this method requires the dispersion relation to have the 4-wavevector expressed as a covariant vector due to the derivative present in (<ref>) while the position 4-vector remains in contravariant form.





§ METHODS
 

We now describe our procedure to produce the expected axion-photon conversion signals around an NS.
We begin by describing the required numerical methods, allowing us to compare the effects of changing the metric, dispersion relation, magnetic field, and conversion probability on the signal received from axion to photon conversion around NSs.



 §.§ Numerical Methods


Numerical methods are used to run simulations of the photons propagating through the plasma surrounding the star, source the axion-photon conversion points, and produce the estimated signal/flux from these phenomena. 
Our choice of programming language is Python[<https://www.python.org/>] due to its relative simplicity and an exhaustive library of modules capable of running the simulations carried out in this work.

Due to the complexity of the coupled ODE systems present from the dispersion relations and the ray tracing equations, we use numerical solvers to compute the paths of the photons through the plasma. 
We employed the Scipy library's solve_ivp function to numerically trace the photons' position and momentum through the dispersive plasma.
At each step, the magnetic field, plasma frequency, and angle between the magnetic field and photon momenta are computed.

To compute the derivatives of functions, the Python library Autograd[<https://github.com/HIPS/autograd>] was employed. Autograd is capable of automatically differentiating native Python and Numpy code. 
It works by using Automatic Differentiation to compute an approximation of the derivative of a function, with machine precision accuracy, without computing a symbolic expression of the derivative.
Hence, only the function must be known, but its related derivatives are unnecessary.
If the plasma function is to be modified, rather than finding its potentially difficult derivative, only the function itself needs to be known for Autograd to return the derivative/gradient of the plasma function in each coordinate direction.

We implement the module Multiprocessing[<https://docs.python.org/3/library/multiprocessing.html>] to leverage the multiple threads available to us during a simulation by propagating multiple photons through the plasma simultaneously. 
This allowed us to divide the time for a simulation to run by approximately the number of available processing threads. 

The code was developed so that changes to the dispersion relation, magnetic field, or probability calculations can be easily made by modifying the defining Python functions. 
This also allows for adding more complex magnetic fields, dispersion relations, or future probability calculations.



 §.§ Observer to Emitter Scheme


In this scheme, photons are sourced from an image plane and back-propagated to the axion-photon conversion surface <cit.>. 
This concept is based on the physical process of photons propagating from the conversion surface and falling incident on a detector found on Earth. 
Based on the fact that asymptotically far from the star, once the photons are no longer refracting due to the plasma, photon trajectories are perfectly straight and will be orthogonal to the detector on incidence and parallel to neighboring ray paths. 
Hence, the propagation of photons can be reversed, beginning orthogonal to and at the detector, and then integrating backwards in time to find the location of the photons from the source.
This means that only photons converging on the detector plane are considered in this simulation scheme.
Because of this, the observer-to-emitter scheme provides a less computationally heavy workload, as the number of photons traced is significantly smaller than that of the emitter-to-observer scheme. 
Another benefit of this method is that it produces the image a detector would `see' from the source. 
However, the viewing angle from the NS needs to be known as it can greatly affect the results.
The following method is adapted from Refs. <cit.>.

Due to the large separation between a detector and the NS, the rays are assumed to be parallel to each other and perpendicular to the detector. Hence, the rays are sourced in such a way at the detector. The validity of this assumption relies on the rays ceasing refraction due to the plasma once ω≫ω_p, which occurs within the light cylinder of the NS. 
Once refraction stops, the rays take on the free-space dispersion relation (<ref>), the solution of which is a ray traveling along geodesics. 
Hence, for a ray to reach the detector at some asymptotic distance, it must travel straight at the detector once the plasma's influence stops. Thus, the ray will be perpendicular to the detector plane and parallel to all other rays reaching the detector.

To run the Observer to Emitter Scheme, photons must first be sourced at an initial position. 
This is done by placing a detector at some substantial separation from the source.
Photons are then integrated backward in time (back-propagated) from the detector until the photons reach the source of axion-photon conversion.
More precisely, the photons are given initial positions along a detector plane (or image plane), which is divided into pixels of side length Δ b. 
A single photon is sourced at the center of each pixel before being propagated through the plasma.
They are initiated to be perpendicular to the detector plane and parallel to each other.
The propagation is described using one of the ray tracing methods with one of the dispersion relations. 
The ray path is then solved by integrating the ray tracing equations over time, or their affine parameter, to track the refraction of the rays through the plasma. Integration is finished when either: (1) the photon intercepts the axion-photon conversion surface; or (2) the photon misses the surface and reaches the end of the integration interval.

From here, the value of the radiated power received by each pixel can be calculated with the photon values at the conversion surface. The total differential power (dP/dΩ dω) received by the detector is given by Eqn. (11) of Ref. <cit.>,

    dP/dΩ dω = ∑_i,jΔ b^2 ω^3 f^i,j_γ,

where the sum is over i, j is the index of the pixels on the detector plane,  Δ b is the pixel side length (so that Δ b^2 is the pixel's area), ω is the frequency/energy of the photon, and f^i,j_γ is the phase space factor of the photons. 
The photon phase space factor can be related to the axion phase space factor via the axion to photon conversion probability, such that f_γ = P_aγγf_a.

We wish to evaluate the quantity (<ref>) at the conversion surface for each photon that is traced to the conversion surface. To do this we need the energy of the photon, conversion probability and the axion phase space all at the point of axion to photon conversion. The energy is found via ray tracing, the conversion probability is found using (<ref>), and the phase space factor of the axions is given by the expression (e.g. see mcdonaldGeneralizedRayTracing2023MW23 Eqns. (58) and (59)),

    f_a(x, k) = v_a n_a,∞k_c(|x|)/k_0δ(ω-ω_c)/4π |k|^2,

where v_a is the velocity of the axions at that point, n_a,∞=ρ_a,∞/m_a is the asymptotic number density of axions, k_c(|x|)^2 = k_0^2 + 2GM_nsm_a^2/|x| the square of the in-fall momentum of the axions, k_0 = m_a v_0 is the momentum dispersion, and ω_c is the energy of the photon at the conversion surface. 
In natural units, |k|^2=k_μ k^μ is the three-momentum magnitude of the axion/photon and is found via ray tracing.



  §.§.§ Implementation


In the following, we describe our implementation of the observer-to-emitter method in GR using the Schwarzschild metric. 
We will use the Schwarzschild coordinate system where x^i=(t, r, θ, ϕ). 
Hence, the 4-momentum will take on the form k^i=(ω, k^r, k^θ, k^ϕ).
Throughout this description, we will also highlight code variables in our algorithm using mygraythis formatting. 
We do this to clarify what each variable controls and what values have been chosen in this paper. 

The detector's observing angle (θ, ϕ) must be chosen for each simulation. 
This gives the center line of the detector relative to z of the NS (recall that z is in line with the rotation axis of the star). 
These angles are related to the code variables mygrayObs_theta and mygrayObs_phi. 
The detector plane distance from the NS also needs to be assigned to mygrayObs_r0. A balance must be struck with this initial distance. 
It should be far enough away to approximate what an extremely distant observer would see but not so far that computation time is significantly increased.
The dispersion relationship must also be picked via mygraydispersion_relation, where the method is chosen with the strings mygray[language=python]'vacuum', mygray[language=python]'isotropic', mygray[language=python]'anisotropic'. 
Lastly a metric must be chosen using mygraymetric_choice with the options mygray[language=python]'flat', mygray[language=python]'schwarzschild', or mygray[language=python]'hartle_thorne'.



Below, we outline the procedure used to produce an axion-photon simulation. 
We used the GR dispersion relations in Sec. <ref> and geometric optics in Subsec. <ref>.
A detailed explanation of the observer-to-emitter algorithm used in this article is as follows:

    
  * The center of a detector plane is initialized by setting a distance mygrayObs_r0 from the NS, and viewing angles mygrayObs_theta and mygrayObs_phi. 
    Pixels are then evenly spaced out on a rectangular detector plane, up to a max size of mygraymax_x and mygraymax_y with the number of pixels along each row and column given by mygraytotal_resolution.
    The center of the detector is aligned with the center of the NS.
    The distance between adjacent pixels is given by Δ b =mygraymax_x/mygraytotal_resolution.

    
  * A single photon is then sourced at the center of each pixel and given initial data:
    
        
  * The initial position of each photon is set as the center of each pixel and given a time mygrayt_0, the starting time for the simulation (typically t_0=0).
        Hence, each photon has an initial 4-position given by,
        
    x^i=(t_0, r_pixel, θ_pixel, ϕ_pixel),

        where pixel refers to the pixel at which the photon is sourced.
        
        
  * The initial 4-momentum of the photon is found using the photon energy and the refractive index.
        Energy conservation with the axion gives the energy of each photon during the conversion. 
        All the axion's energy is converted to the photon's energy. 
        One may then use E^2=m^2+p^2,
        
    ω=√(m_a^2[1+v_a^2(r_pixel)]).

        where ω is the photon's energy, v_a^2=v^2_min+v_∞^2 with v^2_min=r_s/r and v_∞ the dark matter velocity at asymptotic infinity. 
        The equation  (<ref>) must be translated via the effective energy equation ω=-k_iU^i to give the covector form k_i. 
        In the case that the plasma is static, such that U^t=√(-g^tt) and U^μ=0, we have that 
        
    k_t=-ω/√(-g^tt),

        where ω is given by (<ref>).
        
        
  * For the spatial components k_μ, we take the refractive index of the medium at distance mygrayObs_r0 to be n≃1.
        Hence, we will have the covariant relationship from the refractive index (<ref>) as -k_tk^t=k_μ k^μ.
        The centermost photon of the detector is taken to travel perfectly radial from/towards the star at the detector, such that k^α=(k^r, 0, 0).
        So we will have that,
        
    k^r = √(-g_ttk^tk^t/g_rr),
 
        for this photon. 
        Every other photon is assigned the same value to its spatial components but rotated such that the momentum is parallel to the centermost photon's momentum[If this rotation is not done, and instead the momentum of each photon is set identical to (<ref>), the photons will all have velocities towards the center of the star rather than generating rays that are initially parallel.].
    
    
    
  * Once the 
    photons'
     initial conditions are defined, they can then be backpropagated through the medium using the ray tracing equations (<ref>) and (<ref>) along with one of the dispersion relationships (<ref>), (<ref>) or (<ref>).
    This is done by using a numerical ODE solver for a finite time.
    An example of the paths rays will travel along in this during this algorithm is shown in Fig. <ref>.
    The figure also shows the origin from which the photons are sourced on the detector plane with red dots.
    
    
  * At the beginning of the simulation, a coarse pixel search over the entire detector plane is done. 
    This search uses the detector resolution[Resolution here means the number of pixels along the side length of the detector.] defined by mygray|coarse_resolution|.
    The intention of this is to find the regions of the image plane that likely receive photons from the conversion surface.
    This decreases the total number of photon being back tracked in the next step by removing detector regions that will not have photons intercepting the conversion surface. 
    Hence, this procedure decreases the run-time and increases the efficiency of the simulation by removing unnecessary ray tracing of photons.
    
    
  * After the coarse search has identified fine pixels that may have photons that will back-trace onto the conversion surface, a fine search is started over the smaller pixels. 
    This is done by selecting each coarse pixel identified before and completing a higher-resolution search of that coarse pixel. 
    The resolution is defined by mygray|fine_resolution|. 
    Each back-traced photon from a fine pixel that intercepts the conversion surface is recorded, while all photons that never reach the conversion surface are ignored.
    
    
  * Numerical integration of a photon's path ceases when either the photon intercepts the NS surface or the end of the integration interval is reached.
    The photon continues after intercepting the conversion surface, and the solver records the photon position and momentum 4-vectors to find the power from each intersection of the photon with the conversion surface.

    
  * Each photon from the fine search that intercepts the conversion surface can then have its probability of conversion calculated using the intercepts obtained from the ODE solver. 
    Then the photon values and the conversion probability can be used to find the radiated power received by each fine pixel using the terms inside the summation of (<ref>).
    

Once the algorithm is completed, the results can be plotted to form the image on the detector. 
Each radiated power found in the last step is assigned to the index of the fine pixel from which the photon originated.
These indices can be used to reproduce the image on the detector plane using appropriate plotting software.
Otherwise, the powers for a particular viewing angle can be summed together using (<ref>) to find the total estimated power received by a distant observer. 


To improve the effectiveness of the coarse search and avoid missing regions that may receive photons, a relative tolerance mygray|coarse_search_rel_tol| is defined in the search algorithm. 
This alters the search so that coarse pixel photons that approach near the conversion surface, but do not intercept it, are also included. 
This is to avoid the case that fine pixels may be missed if a coarse pixel is excluded by only considering coarse search intercepts. 




§ RESULTS AND DISCUSSION
 




We carry out the observer-to-emitter scheme as detailed above using the parameters in Table <ref>.
Because the effects of GR become greater closer to the Schwarzschild radius, we also include a higher mass NS of M_ns=2.2M_⊙ (compared to the typical M_ns=M_⊙).

We specify the resolution of the detector plane in Table <ref>.
To ensure that the total radiated power had converged, we ran simulations of increasing detector resolution a specified observing polar angles.
Hence, our choice to use a total resolution of 250 × 250 pixels provided a good balance of accuracy and simulation runtime, whilst having converged to a consistent radiated power. 
We found that using fewer pixels than our choice of 250 did not sufficiently alter the radiated power received in the simulation, until we reached a total resolution of 150. 
This total resolution of the detector in each dimension equals the number of coarse resolution pixels times the number of fine resolution pixels.

For the differing axion masses, the conversion surface changes size (e.g. see Fig. <ref>). To maintain a resolution of a constant number of pixels, we therefore have to use different pixel sizes for different axion masses.
In the case of m_a=10, the pixel size is Δ b=120.
For simulations using m_a=1, it has Δ b=480.





 §.§ GJ vs GLP Conversion Surface





Firstly, we explore the effect on the conversion surface by changing the magnetic field model.
Each magnetic field choice alters the charge density, hence the plasma frequency and ultimately the point at which the resonance condition ω_p=m_a is satisfied. 
A 2D cross-section of various conversion surfaces using the GJ and GLP models is shown in Fig. <ref>. 
This plot clearly shows that the conversion surface moves closer to the surface of the star when considering the GLP model. 
However, the general shape of the conversion surface remains the same, with a bulb at either pole and a central torus around the equator. 
Both of these regions of the conversion surface are separated by throats that extend below the NS's surface. 
The throats occur due to a change in sign of the charge density, meaning the plasma frequency tends to zero in the region by the throats.
Hence, from (<ref>), the conversion surface radius also tends to zero.

We also display the effect of altering the NS's mass in <ref>.
The effects of spacetime curvature are enhanced when the mass of the NS is increased, and hence, there is a greater difference between the two magnetosphere models.
This is because the GLP model depends on the mass of the NS[The GJ model has no dependence on the NS's mass as it is derived in flat space. However, even in that case,
the NS mass will still affect the ray tracing by affecting the spacetime geometry.], ultimately affecting the charge density and, hence, the plasma frequency that gives the conversion surface. 
The effect of increasing the NS mass is even clearer for the higher axion mass as the conversion surface is brought closer to the NS when choosing a heavier axion mass.
This is shown for m_a=10 in <ref>(b) as opposed to the lighter mass axion used for (a).
Most notably, the throats' shape, position, and size are altered for the higher mass NS and m_a=10.




 §.§ Isotropic Plasma




By using the simpler isotropic plasma case (<ref>), we can glimpse any initial difference in the total power a distant observer receives. 
The main effect of changing the magnetosphere model will be on the charge density, and hence, through (<ref>), the plasma frequency around the NS.
Ultimately, this will cause the back-traced photons to have altered trajectories through the plasma.
This will also result in a change to the values returned by (<ref>), (<ref>) and (<ref>) at the point of conversion for a photon.


The result of switching the magnetosphere model on the period averaged radiated power across different observing angles is shown in Figs. <ref> and <ref>.
The latter figure uses a higher mass NS in the simulations.
For the results presented in Fig. <ref>, we see only minor differences between the GJ and GLP models.
However, by increasing the mass of the NS, the changes introduced by spacetime curvature will become more important.
Most importantly, increasing the mass significantly alters the GLP magnetosphere model while leaving the GJ magnetosphere model unchanged. 
Changing the mass of the NS also has consequences on the ray tracing due to the changes in the metric. 
Also, the axions will have a greater momentum through the plasma as seen in (<ref>).
The result of all of this can be seen in Fig. <ref>, that by increasing the mass of the NS, the total radiated power is increased.
It also increases the difference between the two magnetosphere models, especially in the m_a=10 case.
This suggests that a GR magnetosphere may be important to consider in the results of searches around higher-mass NSs. 



As a test to ensure that we have implemented the physical processes occurring correctly in our numerical simulations, we compare the results we obtain with the GJ magnetosphere to that of the previously published work of mcdonaldGeneralizedRayTracing2023MW23. 
We see in Fig. <ref> that our results reproduce theirs 
reasonably well.
The slight deviation is likely due to differences in numerical solver tolerances and the resolution of the detector plane. 
We also only need to evaluate the polar angles from 0 to π/2 radians due to the symmetry of the magnetic field.




To study whether there is a reasonable difference between the total radiated power of the GJ and GLP models, we present the absolute percentage difference between our results using the GJ and GLP models in Figs. <ref> and <ref>.
To justify if the difference between the two models is significant, we use the difference between our simulations using the GJ magnetosphere and the results from mcdonaldGeneralizedRayTracing2023MW23.
This provides an estimate of `uncertainty' in implementing the two models.

In both Figs. <ref> and <ref>, the simulations using m_a=1 have large differences present at θ_ obs∼ 1 radian which coincides with the substantial flux received from the throat of the conversion surface.
When considering the conversion surfaces in Fig. <ref>, the throats are not as deep for the GLP case. 
Hence, back-traced photons will not `bounce' off the conversion surface as much down its throat, yielding less radiated power.
In the m_a=10 simulations, this issue is not prevalent due to the throat not being as deep and intercepting the NS surface much closer to the opening. 
Interestingly, the interception of the NS surface leads to a slight dip in radiated power near θ_obs∼ 1 rad in the bottom panel Fig. <ref>.




In Fig. <ref>, for the m_a=10 case, a reasonable difference is present between the GJ and GLP models across most viewing angles.
This difference is also greater than the difference between GJ and mcdonaldGeneralizedRayTracing2023MW23 data.
For the GLP model, from Table <ref>, there is an average absolute difference of 26% over all the viewing angles, which, when compared to the difference with mcdonaldGeneralizedRayTracing2023MW23 having an average of 9.8%, appears to be a significant change in power.
Looking closer at Fig. <ref>, we see that some viewing angles have larger changes in power, while a few have a minimal difference.

In Fig. <ref> for both axion masses, the difference between the GJ and GLP models is similar to the difference between GJ and mcdonaldGeneralizedRayTracing2023MW23.
Hence, the GLP model does not introduce a significant difference in the lower mass NS case.
Hence, for a higher mass NS and a conversion surface that will be close to the Schwarzschild radius, a GR magnetosphere induces a reasonable difference in power.








 §.§ Anisotropic Plasma


Earlier in this article, we discussed the anisotropic plasma relationships (<ref>). 
However, our implementation of this more complex plasma, discussed for example in Ref. <cit.>, did not produce results that reliably matched the results presented by mcdonaldGeneralizedRayTracing2023MW23 for the GJ model. 
Hence, we have excluded any of our results using the dispersion relation (<ref>).
In future work, we will try and determine the reason for the difference. 



§ CONCLUSION


In the work carried out here, we discussed the implementation of the GLP magnetosphere model that is derived in a curved spacetime using the Schwarzschild metric. 
This GR model was then employed in numerical simulations to study the effect on axion-photon conversion signals from NS.
This is in comparison to the recent numerical simulations of mcdonaldGeneralizedRayTracing2023MW23, which uses the flatspace-derived GJ model in their simulations.
Using an isotropic dispersion relation, we compared the difference in total radiated power across the observing angles of the neutron star between the two models that were simulated for this article and the simulation data supplied by mcdonaldGeneralizedRayTracing2023MW23.

For the M_ ns=2.2M_⊙ case, we found the average absolute percentage difference between the GLP and GJ models was 26% for the m_a=10μeV. This was about 2.7 times greater than the average difference between our GJ model and that of mcdonaldGeneralizedRayTracing2023MW23. The m_a=1μeV had an even larger 32% difference between the GLP and GJ models, but this was only about 2.1 times greater than the average difference between our GJ model and that of mcdonaldGeneralizedRayTracing2023MW23 for this case.

In the cases of M_ ns=1.0M_⊙, implementing the GLP model appears to have little effect on the radiated power for an isotropic plasma.

The conclusion that the higher mass case has the greatest difference intuitively makes sense.
This is due to the stronger dependence on spacetime curvature effects when the conversion surface is close to the Schwarzschild radius of the NS.





We thank Jamie McDonald, Harrison Ploeg, and Sam Witte for helpful comments. Sam Witte also for making available some of the simulation data from MW23 to include in our Figs. <ref>, <ref>, and <ref>.
We also thank Alexandru Lupsasca for helpful comments on the difference between the GLP model's choice of coordinates and our definition of the same coordinates.



apsrev4-2-author-truncate.bst
 





§ GLP - INCLINED ROTATOR
 

The third paper, in the series by Gralla et al. <cit.>, extends their work by including a misalignment between the rotation and magnetic field axis. 
They suggest that the results from their first paper can be modified by a spatial coordinate change, using a set of spatial coordinates about the rotation axis and another set about the magnetic field axis. 

The rotation axis Ω is chosen to be in line with the Cartesian axis z, after which the typical spherical polar coordinates (r, θ, ϕ) are defined about this axis.
This is the coordinate system a stationary observer will be using. 
The magnetic field symmetry axis e is inclined by a polar angle χ to the rotation axis and will have azimuthal angle Ω t due to the star's rotation. 
The polar coordinates (r, ϑ, φ) are defined around the axis e such that ϑ is the polar angle measured away from the axis and φ is the azimuthal angle around the axis measured from a line pointing at the rotation axis Ω, such that in this coordinate system the magnetic field is independent of time.
We can also introduce an azimuthal angle, which measures the angle from e to ϕ in the (r,θ,ϕ) frame. Explicitly, this can be expressed as λ=ϕ-Ω t.



Given axisymmetric field functions defined about the magnetic field symmetry axis, the coordinates (r, ϑ, φ) can be used to maintain the field 
functions' 
symmetry.
These coordinates are shown in Fig. <ref>. 
Hence, we may take the two functions from Gralla et al. <cit.> that describe an aligned dipole, (<ref>) and (<ref>), and perform the coordinate change θ→ϑ and ϕ→φ. This results in,

    
    ψ_1 near(r, ϑ)    = μ R^>_1(r)sin^2ϑ,
    ψ_2(φ)    = φ,


where the time dependence no longer exists in the rotated frame due to the axisymmetric field functions.

The coordinate angles ϑ and φ need to be related to the angles θ and ϕ so that we may compute the fields and charge density in an observer's frame. A method to derive these relations is using spherical triangles, and in particular, the relationships are,

    
    cosϑ   = cosθcosχ + sinθcosλsinχ, 
    tanφ   = sinθsinλ/cosθsinχ + sinθcosλcosχ,


which are found using the equations and figures in Appendix <ref>.[These relationships differ from Ref. <cit.>. 
This is due to different coordinate definitions. 
However, the results remain unaffected.] 
The relationships (<ref>), (<ref>), (<ref>) and (<ref>) can all be combined with (<ref>) to give the electromagnetic tensor around an inclined star with GR corrections due to the curved spacetime.
We can then use this electromagnetic tensor with (<ref>) to find the magnetic field strength throughout the star's magnetosphere.

Lastly, we need the charge density for this NS model to provide the plasma. From Ref. <cit.>, in the (r, θ, ϕ) frame, they give charge density as,


    J^t̂=Ω-Ω_z/r(r-2 G M){   ∂_θψ_1 ∂_θ∂_ϕψ_2 - ∂_θψ_2 ∂_θ∂_ϕψ_1 + r(r-2 G M)(∂_r ψ_1 ∂_r ∂_ϕψ_2-∂_r ψ_2 ∂_r ∂_ϕψ_1) 
       -∂_ϕψ_1[(1-2 G M/r) ∂_r(r^2 ∂_r ψ_2)+∂_θ(sinθ∂_θψ_2)/sinθ] 
       + ∂_ϕψ_2[(1-2 G M/r) ∂_r(r^2 ∂_r ψ_1)+∂_θ(sinθ∂_θψ_1)/sinθ]}.


Where (<ref>) was found via (<ref>) using (<ref>) with the coordinate dependencies induced by (<ref>) and (<ref>). For a stationary Schwarzschild observer,[Who will have a 4-velocity U^a=(√(-g^tt), 0, 0, 0).] the charge density ρ_e is given by,

    ρ_e=U_aJ^a=J^t√(1-r_s/r)=J^t̂.


This gives us a well-defined near-field approximation of a misaligned dipole magnetic field in covariant form which accounts for the curvature of spacetime.





§ SOLUTIONS OF SPHERICAL TRIANGLES



The spherical Law of Cosines is given as (e.g. <cit.>),

    cos c = cos a cos b + sin a sin b cos C,

and the spherical Law of Sines is given as,

    sin a/sin A = sin b/sin B,

where the lower case represents the side length and the upper case represents the corresponding angle (see Fig. <ref>).
Upon rearranging the spherical Law of Cosines for cos C, we get,

    cos C=cos c - cos acos b/sin a sin b.

We can replace a, b, and c with the coordinate angles around the rotation axis (θ, ϕ) and around the magnetic field axis (ϑ, φ) (see Fig. <ref>),

    
    cosϑ = cosχcosθ + sinχsinθcosλ,
    sinϑ/sinλ = sinθ/sin-φ,
    cos-φ = cosθ-cosχcosϑ/sinχsinϑ.


All three equations can be combined to remove the dependency on ϑ from (<ref>). By replacing sinϑ using (<ref>) and cosϑ using (<ref>)

    cosφ = -sinφ[cosθ - cosχ(cosχcosθ + sinχsinθcosλ)]/sinχsinλsinθ

where we have also simplified the negative arguments. This then becomes,

    
    cosφ/sinφ   = -cosθ + cos^2χcosθ - cosχsinχsinθcosλ/sinχsinλsinθ,
    tanφ   = sinχsinλsinθ/sin^2χcosθ - cosχsinχsinθcosλ,
    tanφ   = sinλsinθ/sinχcosθ - cosχsinθcosλ.





