PSRCHIVE user documentation: psrmodel

1.0 Purpose

The pulsar modeling program, psrmodel implements a novel algorithm to fit the Rotating Vector Model (RVM) to linear polarization data.

Important Note: The reference frame in which psrmodel fits the model to the data is consistent with the IAU convention for position angle, which may be inconsistent with other RVM-fitting programs (see Section 3 below).

2.0 Usage

Rotating Vector Model Fit

With fingers crossed, you can hope that psrmodel will find the global minimum without any intervention by simply typing 
psrmodel filename
Most likely, this will fail because the code that computes the first guess isn't perfect; therefore, psrmodel can also search a regular grid of first guesses. For example, to search an 18 by 18 grid in which alpha and zeta are varied from 5 to 175 degrees in steps of 10 degrees: 
> psrmodel -s 18x18 filename

psrmodel: search 18X18 grid
Pulsar::ComplexRVMFit::set_observation using psi0=-19.0705
Pulsar::ComplexRVMFit::set_observation using phi0=272.637
Pulsar::ComplexRVMFit::set_observation using delpsi_delphi=-5.49877
Pulsar::ComplexRVMFit::set_observation 29 bins
psrmodel: performing search of chi^2 map

chisq=16.6108/nfree=25 = 0.664431
PA_0=(55.8047+-11.9973) deg
zeta=(79.3135+-26.3904) deg
alpha=(49.7432+-10.6986) deg
phi_0=(189.267+-5.24661) deg
In this example, 29 on-pulse phase bins were chosen and included in the fit. Because the free parameters in the rotating vector model are highly covariant, it is also prudent to inspect the degeneracy between alpha and zeta. To output the reduced chi-squared surface add -x > chisq.map to the command line. The standard output (redirected to chisq.map) can be plotted using gnuplot (it is formatted as understood by the splot command).

The solution can be refined by reducing the threshold to 2 sigma and using the best-fit model from the above global search as a first guess; e.g. 

>  psrmodel -t 2 -d -p 56 -z 79 -a 50 -b 189 -d filename

peak phase=169.453 deg; PA=82.2966 deg
Pulsar::ComplexRVMFit::set_observation 45 bins
psrmodel: solving with initial guess: 
PA_0=56 deg
zeta=79 deg
alpha=50 deg
phi_0=189 deg

chisq=49.3536/nfree=41 = 1.20375
PA_0=(45.6649+-3.5955) deg
zeta=(51.1591+-15.3274) deg
alpha=(36.0193+-9.44647) deg
phi_0=(191.618+-1.41079) deg
The -d option is used to plot the best-fit model over the data.

3.0 Algorithms

Astronomical Convention for Position Angle

As described in Everett & Weisberg (2001), the typical formulation of the Rotating Vector Model is inconsistent with the IAU definition of position angle. To aid in comparison with position angles measured using standard astronomical convention, psrmodel uses an RVM formula that is consistent with the IAU defintion of position angle. Furthermore, to enable comparison of measurements made at different wavelengths, psrmodel first corrects the position angle to zero wavelength (or infinite frequency) using the rotation measure store in the header.

Complex-valued Rotating Vector Model

Rather than perform a one-dimensional fit to the real-valued position angle (P.A.) as a function of pulse, psrmodel performs a two-dimensional fit directly to the Stokes Q and Stokes U profiles by treating them as the real and imaginary components of a complex number. The phase of each complex number is given by the position angle predicted by the rotating vector model, and the magnitude (linearly polarized flux) is modeled as a free parameter.

Information is conserved: two numbers (Q and U) are used instead of one (P.A.) but the extra degree of freedom constrains the magnitude, which is orthogonal to the phase and therefore should not be covariant or cause degeneracy. Modeling the complex-valued L=Q+iU has a number of advantages over modeling the P.A.:

  • the errors on Stokes Q and U are normally distributed;
  • Stokes Q and U are not cyclic;
  • orthogonal mode transitions are trivially modeled by negative values of the complex magnitudes;
  • there is no need to avoid the noise, allowing every last bit of signal to be milked out of the data without biasing the result.
Aside from the magnitudes of the complex numbers, there are four free parameters in the fit:
  1. phi0 - pulse phase of magnetic meridian (steepest slope in P.A.)
  2. psi0 - position angle at phi0
  3. alpha - colatitude of magnetic axis
  4. zeta - colatitude of line of sight

When doing a global minimum search, the first guesses for all of the above parameters are computed using Equation 5 of Everett & Weisberg (2001). That is, the first guesses for psi0 and phi0 correspond to the pulse phase where dpsi/dphi is a maximum. The ratio, sin(alpha)/sin(beta) is then solved to yield the first guess for alpha by assuming beta=2*alpha-pi; the first guess for zeta is then simply alpha + beta.

When mapping out the chi-squared surface, the first guesses of both alpha and zeta are varied from 0 to 180 degress in steps of 180/N and 180/M, respectively where N and M specified with the -s NxM option; the poles are avoided and alpha=zeta is avoided. On each iteration, all parameters are allowed to vary freely. The guess that results in the best solution is solved one last time, and the best fit values are output as the solution.

4.0 Testing and examples

Please see section 2.0

5.0 Known bugs and features that require implementation

  • None at this time.