PsrPopPy: An open-source package for pulsar population simulations
MMon. Not. R. Astron. Soc. , 1–12 (2012) Printed 21 January 2014 (MN L A TEX style file v2.2)
PsrPopPy: An open-source package for pulsar populationsimulations
S. D. Bates , , D. R. Lorimer , , A. Rane and J. Swiggum Department of Physics and Astronomy, West Virginia University, Morgantown, WV, 26506 USA Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK National Radio Astronomy Observatory, PO Box 2, Green Bank, WV 24944, USA
Accepted; 21 January 2014
ABSTRACT
We have produced a new software package for the simulation of pulsar populations,
PsrPopPy , based on the
Psrpop package. The codebase has been re-written inPython (save for some external libraries, which remain in their native Fortran), utilis-ing the object-oriented features of the language, and improving the modularity of thecode. Pre-written scripts are provided for running the simulations in ‘standard’ modesof operation, but the code is flexible enough to support the writing of personalisedscripts. The modular structure also makes the addition of experimental features (suchas new models for period or luminosity distributions) more straightforward than withthe previous code. We also discuss potential additions to the modelling capabilities ofthe software. Finally, we demonstrate some potential applications of the code; first,using results of surveys at different observing frequencies, we find pulsar spectral in-dices are best fit by a normal distribution with mean − . .
0. Second, we model pulsar spin evolution to calculate the best-fit for a relationshipbetween a pulsar’s luminosity and spin parameters. We used the code to replicate theanalysis of Faucher-Gigu`ere & Kaspi, and have subsequently optimized their power-law dependence of radio luminosity, L , with period, P , and period derivative, ˙ P . Wefind that the underlying population is best described by L ∝ P − . ± . ˙ P . ± . and is very similar to that found for γ -ray pulsars by Perera et al. Using this rela-tionship, we generate a model population and examine the age-luminosity relation forthe entire pulsar population, which may be measurable after future large-scale surveyswith the Square Kilometer Array. Key words: pulsars: general - stars: neutron
The known Galactic population of pulsars stands at over2000 sources, but represents only a small fraction of thetotal detectable Galactic population of ∼ Psrpop , and as oneof the few publicly-available population simulation softwarepackages has been commonly used as a tool either for pre-dicting the results of future pulsar surveys (e.g. Keith et al.2010) or for study of some feature of the pulsar population(e.g. Ridley & Lorimer 2010; Bates et al. 2013).In this paper we present a new version of this software,rewritten in the Python programming language, providingincreased modularity, accessibility for new contributors, andmaking it easier to add new features than in the Fortran http://psrpop.sourceforge.net/c (cid:13) a r X i v : . [ a s t r o - ph . I M ] J a n S. D. Bates et al. codebase. In § § § § PsrPopPy uses a series of statistical models in order to gen-erate a simulated Galactic pulsar population. These modelsdescribe the pulse period, luminosity, and spatial distribu-tions, and the currently supported models are described be-low. Since the pulsar population has been so well studied atradio frequencies of ∼ . Normal distribution
A logical starting point in many studies may be to use anormal distribution to describe pulsar periods. Any valuescan be given for the mean and standard deviation of thisdistribution, and periods are drawn from it randomly. FK06used a normal distribution to describe the birth spin periodsof pulsars in their simulations of pulsar evolution.
Log-normal distribution
A more popular model for the current distribution of pulsarspin periods is to use the log-normal distribution. Lorimeret al. (2006) found the best-fit of this model to the knownpopulation of non-recycled pulsars has mean µ = 2 .
7, andstandard deviation σ = − .
34 (in logarithm to the base-10), used as default values in
PsrPopPy . The end point ofthe simulations mentioned above, by FK06, was also a log-normal period distribution; that is, the value of x ≡ log P is drawn from the distribution, f ( x ; µ, σ ) = 1 √ πσ e − ( x − µ ) / (2 σ ) . (1) Millisecond pulsars
Although outside the ‘normal usage’ of
PsrPopPy , the sim-ulation of Galactic millisecond pulsars (MSPs) is an impor-tant problem to address, especially as many modern pulsarsurveys are aimed at the discovery of MSPs with the po-tential for extremely high timing precision (see, e.g. Keithet al. 2010; Bates et al. 2011a; Boyles et al. 2012).Currently,
PsrPopPy provides the MSP period modeldevised by Cordes & Chernoff (1997) which can be used tomodel MSP periods. However, Lorimer (2012) suggests thatthis period distribution does not well describe the currentlyknown population of MSPs, which has grown substantially since that model was proposed. We intend to add more MSPperiod models based upon future results.
Log-normal distribution
The most simple pulse width model draws pulse widths froma log-normal distribution with standard deviation 0.3 inpulse phase, and a mean given by the user.
Using a model of pulsar beam widths
The following discussion describes the pulsar beam modelfirst described for use in pulsar simulations by Smits et al.(2009b), which assumes a simple geometry — pulsars in thismodel have circular beams of radius ρ , the angle between thepulsar’s spin and magnetic axes is given by χ and the anglebetween the magnetic axis and the line of sight to Earth is β .Kramer et al. (1998) found the following empirical relationbetween the value of ρ and spin period P ; ρ = (cid:40) . ◦ P − / , P > . ◦ , P (cid:54) ρ init . The valuesare then dithered by a value, p , drawn from a uniform dis-tribution between − .
15 and 0 . ρ = log ρ init + p. (3)The value of β is then chosen from a uniform distribu-tion as − ρ (cid:54) β (cid:54) ρ , and χ is calculated by χ = arccos q, (4)where q is selected, again, from a uniform distribution suchthat 0 < q <
1. Using these values of χ , β and ρ , the equationsin (cid:18) W (cid:19) = sin (cid:0) ρ (cid:1) − sin (cid:0) β (cid:1) sin α sin( α + β ) (5)is used to calculate the pulse width (e.g. Gil et al. 1984). Power law
Use of a power-law to describe the distribution of pulsar lu-minosities has been common since studies of the pulsar pop-ulation began. For example, Taylor & Manchester (1977) fit-ted a power-law with index d log N/ d log L = − . ∼ − . c (cid:13)000
Use of a power-law to describe the distribution of pulsar lu-minosities has been common since studies of the pulsar pop-ulation began. For example, Taylor & Manchester (1977) fit-ted a power-law with index d log N/ d log L = − . ∼ − . c (cid:13)000 , 1–12 srPopPy Log-normal distribution
To avoid using a luminosity cut-off (as discussed by Rid-ley & Lorimer 2010), pulsar luminosities are drawn froma log-normal distribution (see Equation 1) with a mean of (cid:104) log L (cid:105) = − . σ log L = 0 .
9, fol-lowing FK06.
Parameterised in P and ˙ P For evolutionary models where a pulsar is given both a ro-tation period, P , and a period derivative, ˙ P , the radio lu-minosity may be parameterised as see, e.g. FK06,log L = log ( P α ˙ P β γ ) + L corr , (6)where L corr is a dithering factor chosen from a normal distri-bution centred on zero with a variable standard deviation. Simple models
Three basic models are provided for distributing the pul-sars around the Galaxy. In the disk model, pulsars are dis-tributed along the Galactic plane ( z = 0), with randomGalactic x - y coordinates from −
15 to +15 kpc. In the slab model, Galactic x and y coordinates are calculated as withthe disk distribution, but z -coordinates run from − isotropic model, pulsarsare distributed randomly around the Earth at a distance of1 kpc. Radial distribution
More complex radial distribution models which attempt totake into account the structure of the Galactic disk havealso been incorporated into
PsrPopPy ; the first is the ra-dial model produced by Yusifov & K¨u¸c¨uk (2004). The sec-ond, and default, is an adaptation of this model based uponthe analysis of Lorimer et al. (2006). Thirdly, pulsars canbe distributed around the Galactic centre using a Gaussianradial density profile of variable width (e.g. Narayan 1987).
Galactic scale height
The distribution of pulsars in Galactic z coordinates is com-monly approximated by a two-sided exponential (e.g. Lyne1998; Lorimer et al. 2006). The default value used in Psr-PopPy is 330 pc, as obtained by Lorimer et al. (2006), butcan be set at any value — for instance, the scale height forMSPs is considerably higher than for non-recycled pulsars.Recent work by Levin et al. (2013) obtains a best-fit MSPscale height of 500 pc.
PsrPopPy also supports the use ofGaussian distributions with mean of zero and variable widthto model pulsar Galactic scale height.
Two popular models for the conversion between dispersionmeasure (DM; the column density of free electrons along a given line of sight) and true distance are included in
Psr-PopPy . The popular Cordes & Lazio (2002) model is in-cluded which is, to date, the best available model for ob-taining distance estimates for pulsars at a given DM. Thelibrary for performing this transformation has been kept inFortran; sometimes it can be prohibitively slow, so the modelfrom Lyne et al. (1985) is also included. The recent modelsby Schnitzeler (2012) are not yet included, but work to in-clude these and other electron density models is ongoing.
We model variations in pulse intensity, known as scintilla-tion, caused by the relative motions of the Earth, pulsarsand the interstellar medium. Our discussion of this effect isbased upon the outline given by Lorimer & Kramer (2005).If a pulsar of flux S is seen to vary with standard de-viation σ s , the intensity of these variations is described bythe modulation index, m , m = σ s / (cid:104) S (cid:105) , (7)whose value must be computed at observing frequency f from the scintillation strength, u = (cid:115) f ∆ f DISS , (8)where ∆ f DISS is the diffractive scintillation bandwidth,given by ∆ f DISS = 2 πτ s C , (9)where τ s is the scattering time, given by Equation 28, andwe assume a Kolmogorov spectrum, and hence the constant C = 1 . u <
1, known as weak scintillation,the modulation index is simply computed as m = √ u / , (10)but for u >
1, strong scintillation, the situation becomesmore complicated, and we must compute modulation in-dices for both refractive and diffractive scintillation, m RISS and m DISS respectively. The total modulation index is thencomputed as m = m + m + m RISS m DISS . (11)The modulation index of refractive scintillation is sim-ply related to the scintillation strength as m RISS = u − / , (12)while the modulation index for diffractive scintillation is re-lated to the number of scintles, N , sampled in the observingtime ∆ t obs and observing bandwidth ∆ f obs , m DISS = 1 √ N t N f , (13)where N t ≈ κ ∆ t obs ∆ t DISS (14) c (cid:13) , 1–12 S. D. Bates et al. and N f ≈ κ ∆ f obs ∆ f DISS (15)and we take an average value of κ = 0 .
15. In the strongregime, we rely on the NE2001 model to compute the valuesof ∆ f DISS and ∆ t DISS (Cordes & Lazio 2002).In the cases of both weak and strong scintillation, thevalue of m is then applied to the pulsar (for which flux S has already been calculated based upon the distance to thepulsar and its luminosity). Following the definition of themodulation index in Equation 7, the value of σ s can be com-puted. A modified flux value is then chosen at random from anormal distribution with mean S and standard deviation σ s ,and is applied only for the current realisation in the simula-tion. The original value of S is stored so that the modulationwill vary each time code is used to simulate a pulsar survey(see § PsrPopPy allows radio spectral indices to be normally dis-tributed with a given mean and standard deviation α and β , respectively. The default values are α = − . , σ = 0 . α = − . , σ = 0 .
2, and further work by Bateset al. (2013) finds an underlying spectral index distributiongiven by α = − . , σ = 1 . For evolutionary simulations of the pulsar population,
Psr-PopPy includes models for generating period derivatives foreach pulsar, based upon work by FK06 and Contopoulos &Spitkovsky (2006). The method follows that discussed byRidley & Lorimer (2010) and includes a variety of pulsarbeaming models and the option of including a decay in theangle between the spin and magnetic axes of each pulsar.Although the evolution code is contained in a separateexecutable, the output models are simply serialised versionsof the models stored in memory. This format is identicalto that used in the rest of
PsrPopPy , and so the modelsare completely compatible with the output from the snap-shot simulations. As well as the distributions outlined in §
2, additional parameters need to be modelled, using thedistributions described below.
By default, the pulsar magnetic fields are selected from alog-normal distribution with mean (cid:104) log B (cid:105) = 12 .
65 andstandard deviation σ log B = 0 .
55. These values are chosenper FK06, but may be altered.
Models of the pulsar spindown, discussed in § χ , between the rotational and magnetic axes ofthe pulsar.The most simple model treats all pulsars as orthogonalrotators; that is, χ = 90 ◦ for every simulated pulsar, while a more realistic model aligns the spin and magnetic axes atrandom. The value of χ is calculated in the same way as inEquation 4.Additionally, a beaming model is implemented basedupon the the work of Weltevrede & Johnston (2008) whofound evidence for an alignment of the magnetic and rota-tional axes over a timescale t d ∼ × yr. Therefore, weinclude the possibility of χ decaying from an initial value χ as sin χ = sin χ exp( − t/t d ) (16)after a time t . Following Ridley & Lorimer (2010), this for-mula does not exactly replicate the model as described byWeltevrede & Johnston, but describes a simplified version. Magnetic dipole model
The general expression for the pulse period, P as a functionof time in the magnetic dipole model (see Ridley & Lorimer2010, for details) for a pulsar with braking index n is givenby P ( t ) = (cid:20) P n − + n − t d kB sin χ (1 − exp( − t/t d )) (cid:21) n − (17)where the constant k = 8 π R Ic (18)in cgs units, where we assume the parameters of the canon-ical pulsar of radius R and moment of inertia I , assumed tobe 10 cm and 10 g cm respectively. Equation 17 is usedto calculate a value for P at time t , and then the periodderivative is calculated as P n − ˙ P = kB sin χ, (19)where, for models with time-varying χ , we use χ ( t ) instead.The so-called “death line” which demarcates the area of the P - ˙ P diagram in which few radio pulsars are observed. Bhat-tacharya et al. (1992) described this line as P death = (cid:114) B . × G , (20)and that pulsars with period P > P death will be radio-quiet.
CS06 Model
For the pulsar spindown model of Contopoulos & Spitkovsky(2006), the pulsar death line is represented by the equation P death = (cid:20) . × (cid:18) B G (cid:19) (cid:18) P (cid:19)(cid:21) n +1 (21)and, as with Equation 20, radio emission ceases when P >P death . However, P and ˙ P are now calculated by integrating c (cid:13)000
For the pulsar spindown model of Contopoulos & Spitkovsky(2006), the pulsar death line is represented by the equation P death = (cid:20) . × (cid:18) B G (cid:19) (cid:18) P (cid:19)(cid:21) n +1 (21)and, as with Equation 20, radio emission ceases when P >P death . However, P and ˙ P are now calculated by integrating c (cid:13)000 , 1–12 srPopPy ˙ P = 3 . × − (cid:18) PP (cid:19) − n (cid:18) B G (cid:19) (cid:18) P (cid:19) − (cid:18) − PP death cos χ (cid:19) (22)with respect to time, and then solving for the period. Once a pulsar’s time-evolved period and period derivativehave been calculated, it must be assigned a position in theGalaxy. Again, following FK06, initial positions are chosenin the following way:(i) A radial position is chosen using the radial model ofYusifov & K¨u¸c¨uk (2004);(ii) The pulsar is positioned along one of the Galactic spiralarms, at the radius given in the previous step. The positionis represented by x - y coordinates in the plane of the Galaxy;(iii) The pulsar is assigned a third coordinate, z , perpendic-ular from the Galactic plane, using the same method as in § x , y and z di-rections are typically assigned from a Gaussian distributioncentered on 0 km s − with a width of 180 km s − . The meanand standard deviation of this distribution may be varied,and additional distributions are simple to implement. Thepulsar is then evolved from its initial position and velocityfor a time equal to the age of the pulsar, using the modelby Carlberg & Innanen (1987) as modified by Kuijken &Gilmore (1989), giving a final position of the pulsar in theGalaxy. Models of the Galactic electron distribution maythen be applied as discussed in § All of the algorithms discussed in § Psrpop executables. These scripts form apipeline for creating a model population and applying dif-ferent pulsar survey parameters to it.The processes for generating synthetic populations us-ing both the ‘snapshot’ and evolutionary methods are dis-cussed in § § § Pulsars in the model population can be run through a seriesof survey parameters to see if they would, in theory, be de-tected in such a survey. As we saw in the previous section,this can be used to constrain the population based uponknown detections, but equally this could be used to make
Figure 1.
Screenshot of the GUI for inspecting population mod-els. Shown are the whole model population generated by
Psr-PopPy , and the subset detected in a simulated Parkes Multibeamsurvey.
Table 1.
Observational parameters required for simulating a pul-sar survey using
PsrPopPy . The examples are the Parkes south-ern pulsar survey (PKS70, Manchester et al. 1996), the Parkesmulti-beam pulsar survey (PMPS, Manchester et al. 2001) andthe Parkes 6.5 GHz multi-beam pulsar survey (MMB, Bates et al.2011b). PKS70 PMPS MMBDegradation factor, β G (K Jy − ) 0.64 ∼ . t obs (s) 157.3 2100 1055Sampling interval, t samp ( µ s) 300 250 125System Temperature, T sys (K) 35 25 40Centre frequency, f (MHz) 436 1352 6591Bandwidth, BW (MHz) 32 288 576Channel width, ∆ f (MHz) 0.125 3 3Number of polarizations, n p ◦ ) 0 0 0Max. RA ( ◦ ) 360 360 360Min. Dec ( ◦ ) − − − ◦ ) 0 +90 +90Min. Galactic longitude ( ◦ ) − − − ◦ ) +180 +50 +30Min.Galactic latitude ( ◦ ) − − − . ◦ ) +90 +6 +0 . projections about future surveys with hypothetical surveyparameters. Describing survey parameters
Surveys parameters are defined in plain text files, making iteasy for users to add or edit their own surveys. Examples c (cid:13) , 1–12 S. D. Bates et al.
Table 2.
Parameters used to simulate the pulsar population in § z -scale height 50 pcLuminosity distribution Log-normal (cid:104) log L (mJy kpc ) (cid:105) − . L (mJy kpc ))) 0 . (cid:104) α (cid:105) − . α ) 0.96Initial period distribution Gaussian (cid:104) P (ms) (cid:105) B field ditribution Log-normal (cid:104) log B (G) (cid:105) . B (G)))) 0 . of the parameters used in these ‘survey files’ are shown inTable 1, and can then be used to predict whether the surveywould be able to detect each of the pulsars.For added precision, instead of describing the surveyregion by the bounds in either Galactic or equatorial coor-dinates, it is possible to provide a list of pointing coordinates(in either coordinate system), with associated gain and ob-servation length values, for the survey. This may be useful,for example, for drift-scan surveys which are not easily de-scribed by a bounding region or for on-going surveys whichcould also be described, though less precisely, using a verylow completed fraction, or for multibeam surveys which havehighly variable gain values away from the central beam. Detection thresholds
The first test to be done is whether the pulsar is inside theregion described by the RA/Dec or l / b ranges in the textfile. Any pulsars which are inside this region then have theirtheoretical signal-to-noise ratioS / N = S ν G (cid:112) n p t obs ∆ fβT tot × (cid:114) − δδ (23)where total temperature is the sum of the system, sky andcosmic microwave background (CMB) temperatures, T tot = T sys + T sky + T CMB , (24)and many of the other terms are defined in Table 1. Thepulse duty cycle, δ = W eff /P, (25) where W eff is the effective pulse width and P the pulse pe-riod.The effective pulse width is not the intrinsic width ofpulses from the pulsar. It is given by W eff = (cid:113) W + t + ∆ t + τ (26)where W int is the intrinsic pulse width, t samp is the samplingtime of the hardware used to record survey data, and∆ t = 8 . × ms × DM × ∆ f MHz f . (27)is the dispersive smearing across a single frequency channelof bandwidth ∆ f MHz at frequency f MHz , both in units ofMHz.The final term in Equation 26 is the pulse smearing dueto scattering by free electrons in the interstellar medium, τ s .For this, we use the empirical scattering fit of Bhat et al.(2004),log τ s = − .
46 + 0 .
154 log DM + 1 . DM) − .
86 log f GHz , (28)for a pulsar with dispersion measure DM, and where fre-quency f GHz is now given in GHz. Since there is a largescatter about this relationship, we pick the final value of τ s from a Gaussian distribution centred on the value computedfrom Equation 28. To enable investigations of this scatteringrelationship, users are also able to use custom values for thefrequency coefficient in this equation.During a pulsar survey where the sky is ‘tiled’ with ob-servations, pulsars are not commonly found at the centre ofthe survey beam, rather, they are offset by some angle whichcauses them to be detected with a lower S/N than if theywere positioned at the beam centre. There are two methodsincorporated in PsrPopPy to reproduce this behaviour. Inthe most simple model, following Lorimer et al. (1993), aGaussian telescope beam is assumed, where the gain, G atany offset, r , from the beam centre is given by G = G exp (cid:18) − . r w (cid:19) (29)for a telescope with FWHM w and where the gain at thecentre of the beam is G . The square of the offset is chosenfrom a pseudorandom uniform distribution, 0 (cid:54) r (cid:54) w / G = G (cid:18) ( ka sin r ) ka sin r (cid:19) (30)where J is a Bessel function of the first kind with index 1, k = 2 π/λ is the wavenumber for an observing wavelength λ ,and a is the effective aperture radius. This is slightly slower,but may be useful in certain situations — for example thePALFA survey at Arecibo, where the first sidelobe is roughlyas sensitive at the main beam of the Parkes radio telescope(see Swiggum et al. in press).More accurate values of the modified gain can be ob-tained by using a list of survey beam positions and, in thecase of multibeam surveys, the value of G for each obser-vation. From this list, the offset from each pulsar to the c (cid:13)000
86 log f GHz , (28)for a pulsar with dispersion measure DM, and where fre-quency f GHz is now given in GHz. Since there is a largescatter about this relationship, we pick the final value of τ s from a Gaussian distribution centred on the value computedfrom Equation 28. To enable investigations of this scatteringrelationship, users are also able to use custom values for thefrequency coefficient in this equation.During a pulsar survey where the sky is ‘tiled’ with ob-servations, pulsars are not commonly found at the centre ofthe survey beam, rather, they are offset by some angle whichcauses them to be detected with a lower S/N than if theywere positioned at the beam centre. There are two methodsincorporated in PsrPopPy to reproduce this behaviour. Inthe most simple model, following Lorimer et al. (1993), aGaussian telescope beam is assumed, where the gain, G atany offset, r , from the beam centre is given by G = G exp (cid:18) − . r w (cid:19) (29)for a telescope with FWHM w and where the gain at thecentre of the beam is G . The square of the offset is chosenfrom a pseudorandom uniform distribution, 0 (cid:54) r (cid:54) w / G = G (cid:18) ( ka sin r ) ka sin r (cid:19) (30)where J is a Bessel function of the first kind with index 1, k = 2 π/λ is the wavenumber for an observing wavelength λ ,and a is the effective aperture radius. This is slightly slower,but may be useful in certain situations — for example thePALFA survey at Arecibo, where the first sidelobe is roughlyas sensitive at the main beam of the Parkes radio telescope(see Swiggum et al. in press).More accurate values of the modified gain can be ob-tained by using a list of survey beam positions and, in thecase of multibeam surveys, the value of G for each obser-vation. From this list, the offset from each pulsar to the c (cid:13)000 , 1–12 srPopPy Figure 2.
Grey-scale plots of likelihoods overlaid with confidence-level contour lines (as defined in Bates et al. 2013) for populationsimulations with varying spectral index parameters α and σ (see Section 2.7). nearest survey beam is calculated and then used in one ofEquations 29 or 30.Using these relationships, the value of S/N can then becomputed using Equation 23, and if it is greater than thedetection threshold (see Table 1), then the pulsar is countedas detected by the survey. PsrPopPy also keeps track ofhow many pulsars are outside the survey region, are smearedout completely (that is, W eff > P ) or are simply too faintto be detected. Results are then reported and can be storedin a text file. A new population model (in the same formatas that generated in § PsrPopPy also records the number of discoveries (notonly detections) in each survey. This allows the potential offuture surveys to be more carefully calculated.
The most basic method for creating a model population isas follows(i) the user selects a number of pulsars to be generated;(ii) for each simulated pulsar in turn, values for each of thepulsar parameters are drawn from the user-specified (or de-fault) distributions.This method might be suitable for situations where a Galac-tic population of X pulsars is hypothesised, which could thenbe tested from simulating survey results, or other means.Commonly, users wish to generate a population basedupon constraints provided by the large-scale pulsar surveyswhich have been performed to date. In this case, the usercan provide a list of surveys they wish to use to constrain themodel, and the total number of pulsars, n , that should bedetected in these surveys. Then, the method differs slightly.(i) the code continually generates new synthetic pulsars;(ii) for each simulated pulsar in turn, values for each of thepulsar parameters are drawn from the user-specified (or de-fault) distributions; Figure 3.
Marginalized probability density functions obtainedfor α (left) and β (right). Best fits are shown with a solid line,giving the results ¯ α = − .
4, ¯ β = 0 . (iii) each pulsar is run through each of the simulated surveysin turn (see § n , the loop terminates. For evolutionary simulations of the pulsar population,
Psr-PopPy includes models for generating period derivatives foreach pulsar, based upon work by FK06 and Contopoulos &Spitkovsky (2006). The method follows that discussed byRidley & Lorimer (2010) and includes a variety of pulsarbeaming models and the option of including a decay in theangle between the spin and magnetic axes of each pulsar.Although the evolution code is contained in a separateexecutable, the output file format is identical to that usedin the rest of
PsrPopPy , and so is completely compatible c (cid:13) , 1–12 S. D. Bates et al.
Figure 4. P - ˙ P diagram of pulsars detected in model PMSURV and SWIN surveys produced using evolve and the luminosity distributionparameters α = − . β = 0 .
48 (bold points). P - ˙ P values that were taken from the pulsar catalogue are shown for comparison (smallerpoints). Histograms of the P and ˙ P distributions are also shown, with grey bars for catalogue sources, and with solid lines for the simulatedpulsars. For comparison, step histograms (dotted) are shown for a simulation using the values α = − . , β = 0 . with the output from the snapshot simulations. The Evolve code is outlined below;(i) the code continually generates new synthetic pulsars asin § τ init , of thepulsar is chosen from a flat ditribution between 0 and a max-imum age, t max (corresponding to the age of the Galaxy),and a magnetic field is chosen from a log-normal distribu-tion;(iii) an alignment angle between the spin and magnetic axesis chosen, and a decay timescale for this alignment can bechosen, if required;(iv) the pulsar is assigned a braking index, and then oneof the two spindown models is applied. If the user wants toapply a ‘deathline’, the pulsar is evaluated to see if it hascrossed into the ‘dead’ region;(v) each ‘live’ pulsar is randomly assigned a birth positionin Galactic plane and a birth velocity;(vi) the pulsar is then evolved through the Galactic poten-tial for a time equal to the age of the pulsar;(vii) the pulsar is assigned a radio luminosity and spectral index, and run through any selected model surveys, as dis-cussed in § Once a population model has been created, with the evo-lutionary or snapshot method, a common use of the modelis to run alternate pulsar surveys over the model — eitherin order to test the input distributions in some way (forexample, using surveys at multiple frequencies to constrainspectral evolution), or to make predictions of how many pul-sars you might expect future surveys to detect or discover.The script dosurvey is provided for exactly this purpose,which performs the following steps:(i) a specified population model is read in; c (cid:13)000
48 (bold points). P - ˙ P values that were taken from the pulsar catalogue are shown for comparison (smallerpoints). Histograms of the P and ˙ P distributions are also shown, with grey bars for catalogue sources, and with solid lines for the simulatedpulsars. For comparison, step histograms (dotted) are shown for a simulation using the values α = − . , β = 0 . with the output from the snapshot simulations. The Evolve code is outlined below;(i) the code continually generates new synthetic pulsars asin § τ init , of thepulsar is chosen from a flat ditribution between 0 and a max-imum age, t max (corresponding to the age of the Galaxy),and a magnetic field is chosen from a log-normal distribu-tion;(iii) an alignment angle between the spin and magnetic axesis chosen, and a decay timescale for this alignment can bechosen, if required;(iv) the pulsar is assigned a braking index, and then oneof the two spindown models is applied. If the user wants toapply a ‘deathline’, the pulsar is evaluated to see if it hascrossed into the ‘dead’ region;(v) each ‘live’ pulsar is randomly assigned a birth positionin Galactic plane and a birth velocity;(vi) the pulsar is then evolved through the Galactic poten-tial for a time equal to the age of the pulsar;(vii) the pulsar is assigned a radio luminosity and spectral index, and run through any selected model surveys, as dis-cussed in § Once a population model has been created, with the evo-lutionary or snapshot method, a common use of the modelis to run alternate pulsar surveys over the model — eitherin order to test the input distributions in some way (forexample, using surveys at multiple frequencies to constrainspectral evolution), or to make predictions of how many pul-sars you might expect future surveys to detect or discover.The script dosurvey is provided for exactly this purpose,which performs the following steps:(i) a specified population model is read in; c (cid:13)000 , 1–12 srPopPy (ii) survey models are created from text files outlining thesurvey parameters. Multiple surveys may be selected;(iii) for each survey, in turn, the pulsars in the populationmodel are run through the equations listed in § Two scripts have been developed as part of
PsrPopPy toallow users to visually inspect the population models gener-ated by the processes discussed in § § § PsrPopPy provides an interactive graphicaluser interface for plotting pulsar parameters against one an-other (for a screenshot see Figure 1). Users can select modelsin the current working directory to be plotted, select whichparameters to plot, and choose to plot using linear or loga-rithmic x and y axes. Plots can also be saved in various fileformats, including vector formats such as postscript. PsrPopPy has been designed to allow user-generatedscripts to be simple to implement. All population modelsare written in a uniform format and, in fact, are simply se-rialised objects from the code. This means that no code isrequired to parse models back into memory, and that high-level scripts can be used to pass a model directly betweencodes, without having to write to disk and then read it outagain. This also removes the need to modify the codebasewhenever additional parameters are added to the models,making the code more robust.
In order to test the ability of
PsrPopPy to reproduce re-sults previously obtained with
PsrPop , we first ran simula-tions following Bates et al. (2013, see § α = − . β = 1 . As discussed in § P , and period derivative, ˙ P , using the formulation show in Equa-tion 6. In this section, we use the evolutionary simulationcode evolve to estimate the value of the parameters α and β in this equation. We kept value of γ in this equation fixedat 0 .
18 mJy kpc , the optimal value according to FK06.For values of α and β covering a grid covering the range − . < α < . . < β < .
7, model pulsar popula-tions were generated using evolve as per the specificationsoutlined in Table 2. The populations were grown until 1206pulsars per population were detected in simulations of theParkes multibeam pulsar survey (PMSURV, Manchesteret al. 1996) and the two Swinburne pulsar surveys at higherGalactic latitudes (Edwards et al. 2001; Jacoby et al. 2009).This is the number of normal pulsars (defined as
P >
30 ms,˙
P < × − ) listed in the pulsar catalogue (Manchesteret al. 2005) as detected by any one or more of these threesurveys.For each grid point ( α, β ), 25 realisations of the sim-ulation were performed. For each, the P and ˙ P values ofthe resulting population model were compared to those ofthe 1206 pulsars listed in the pulsar catalogue using the 2-dimensional Kolmogorov-Smirnov test (Press et al. 1986).While the probabilities generated by this algorithm are notstatistically rigorous, this test does provide a helpful way ofquantifying how well two 2-d distributions match one an-other. The average value of this probability over the 25 re-alisations was stored, as were the models generated by do-survey .Marginalising over α and then β in turn, by averag-ing the probabilities for each α or β value, the distributionsshown in Figure 3 were obtained. Fitting Gaussian func-tions to these probability distributions, we obtained best-fitvalues of ¯ α = − . ± .
09 and ¯ β = 0 . ± .
04. An exam-ple P - ˙ P diagram generated using these values for α and β is then shown in Figure 4. Histograms are also shown fora simulation using the values α = − . , β = 0 .
7, and areclearly skewed to longer periods and higher values of theperiod derivative.It is worth noting that while we have derived values of α and β by comparing our simulations to the known pulsarpopulation, we have obtained a very similar result to Per-era et al. (2013). Their analysis considered the relationshipfor gamma-ray luminosities, P and ˙ P , and derived values of α = − . ± .
03 and β = 0 . ± .
02. This one-to-one cor-respondence between empirical luminosities is remarkablegiven that the fraction of the spin-down energy loss goinginto gamma-ray emission is substantially greater than in theradio. The X-ray emission does not scale in a similar man-ner for rotation-powered pulsars (e.g. Possenti et al. 2002;Kargaltsev et al. 2012). Further theoretical studies to ex-plain these results are definitely warranted, and we plan toexplore this further in future work.
Having, in the previous section, constrained the values of α and β in Equation 6, we then should be able to see howpulsar luminosities vary with age (both the ‘real’ age, τ real inour simulations, and the characteristic age, given by τ char = P/ P ).First, a model was created using the best values α = − . β = 0 .
48. Using only pulsars in the population which c (cid:13) , 1–12 S. D. Bates et al.
Figure 5.
Histograms of pulsar luminosities for each ‘real age’bin. Fitted Gaussians are shown with a solid black line.
Figure 6.
Normalised histograms (grey bars) of pulsar luminosi-ties for each characteristic age bin. Fitted Gaussians are shownwith a solid black line. The dotted lines show normalised his-tograms of the pulsars detected in a model all-sky survey withthe SKA.
Figure 7.
Mean (upper) and standard deviation (lower) values,with associated errors, from the fits performed in Figures 5 (left)and 6 (right). Upper panels are fit by an exponential decay, thelower panels with a linear decay. lie above the deathline (and are, therefore, considered to beemitting radiation), the population was then divided intologarithmic bins in both real and characteristic ages (withcentre values log τ = 4 . , . ..., . τ in years). For eachof the age bins, histograms were made of pulsar luminosity(shown in Figures 5 and 6). The histograms were fitted withGaussians of mean, µ and standard deviation σ , and thefitted values of µ and σ , with the errors on these valuestaken from the resultant covariance matrix, are then plottedin Figure 7.The results from Figure 7 are that we can characterisethe change in pulsar luminosity, L , with the logarithm ofage (both real age and characteristic age) as an exponentialdecay, log L ∝ exp (cid:18) − log τ log τ decay (cid:19) + a. (31)We obtained best fits of log τ decay = 2 . ± . a = − . ± . τ decay = 2 . ± . a = − . ± . σ L = b log τ + c. (32)We obtained best fitting values of b = − . ± . , c =1 . ± .
05 for pulsar real ages, and b = − . ± . , c =1 . ± .
05 in the case of characteristic ages.Our simulations show, however, that the currently-known pulsar population (simulated using a sub-sample de-tected by PMSURV and the Swinburne pulsar surveys) doesnot contain enough low-luminosity pulsars to be able to mea-sure this age-luminosity relationship. However, as shown inFigure 6, the large number of pulsars expected to be discov-ered in all-sky surveys using the SKA (see, e.g., Smits et al.2009a), may begin to probe this relationship. c (cid:13)000
05 in the case of characteristic ages.Our simulations show, however, that the currently-known pulsar population (simulated using a sub-sample de-tected by PMSURV and the Swinburne pulsar surveys) doesnot contain enough low-luminosity pulsars to be able to mea-sure this age-luminosity relationship. However, as shown inFigure 6, the large number of pulsars expected to be discov-ered in all-sky surveys using the SKA (see, e.g., Smits et al.2009a), may begin to probe this relationship. c (cid:13)000 , 1–12 srPopPy We have developed a software packge for the simulation ofthe Galactic pulsar population. The package has been writ-ten mainly in Python, with reliance on some older Fortrancode for complex calculations. Using this software, we havedemonstrated possible applications of the code: • in § − . . • in § P , and period deriva-tive, ˙ P . We obtained best-fit values for the power-law indices α and β of α = − . ± .
09 and β = 0 . ± . • in § L falls exponentially with log τ with adecay constant of 2.1–2.6.Further functionality will be added to the PsrPopPy code in the future, including: • an improved treatment of millisecond pulsars, includingmodels of binary parameters; • using these binary parameters to model the effect of bi-nary motion on signal strength in pulsar surveys (e.g. Bagchiet al. 2013); • the inclusion of pulse intensity distributions, and theextension of these to include rotating radio transients(RRATs); • model the effects of scintillation on low-DM pulsars (e.g.Cordes & Chernoff 1997).Once some of these features have been implemented,it may become possible to repeat the simulations we haveperformed here, or others similar to them, instead focussingon the millisecond pulsars, or RRATs. It may also be use-ful to include these pulsars in predictions for future surveys.Of course, as the number of such known sources increases,there will be better statistical models to be included in Psr-PopPy , in turn increasing the accuracy of the software.
ACKNOWLEDGEMENTS
The authors acknowledge support from WVEPSCoR in theform of a Research Challenge Grant. DRL is also supportedby the Research Corporation for Scientific Advancement asa Cottrell Scholar.We thank the referee for helpful comments on themanuscript.
REFERENCES
Bagchi M., Lorimer D. R., Wolfe S., 2013, MNRAS, 432,1303Bates S. D. et al., 2011a, MNRAS, 416, 2455Bates S. D. et al., 2011b, MNRAS, 411, 1575Bates S. D., Lorimer D. R., Verbiest J. P. W., 2013, MN-RAS, 431, 1352Bhat N. D. R., Cordes J. M., Camilo F., Nice D. J., LorimerD. R., 2004, ApJ, 605, 759Bhattacharya D., Wijers R. A. M. J., Hartman J. W., Ver-bunt F., 1992, A&A, 254, 198Boyles J. et al., 2012, ArXiv e-printsBurgay M. et al., 2006, MNRAS, 368, 283Carlberg R. G., Innanen K. A., 1987, AJ, 94, 666Contopoulos I., Spitkovsky A., 2006, ApJ, 643, 1139Cordes J. M., Chernoff D. F., 1997, ApJ, 482, 971Cordes J. M., Lazio T. J. W., 2002, ArXiv Astrophysicse-prints (astro-ph/0207156)Edwards R. T., Bailes M., van Straten W., Britton M. C.,2001, MNRAS, 326, 358Faucher-Gigu`ere C. A., Kaspi V. M., 2006, ApJ, 643, 332,in pressGil J., Gronkowski P., Rudnicki W., 1984, A&A, 132, 312Jacoby B. A., Bailes M., Ord S. M., Edwards R. T., Kulka-rni S. R., 2009, ApJ, 699, 2009Kargaltsev O., Durant M., Pavlov G. G., Garmire G., 2012,ApJS, 201, 37Keith M. J. et al., 2010, MNRAS, 409, 619Kramer M., Xilouris K. M., Lorimer D., Doroshenko O.,Jessner A., Wielebinski R., Wolszczan A., Camilo F.,1998, ApJ, 501, 270Kuijken K., Gilmore G., 1989, MNRAS, 239, 651Levin L. et al., 2013, MNRAS, 434, 1387Lorimer D. R., 2012, ArXiv e-printsLorimer D. R., Bailes M., Dewey R. J., Harrison P. A.,1993, MNRAS, 263, 403Lorimer D. R. et al., 2006, MNRAS, 372, 777Lorimer D. R., Kramer M., 2005, Handbook of Pulsar As-tronomy. Cambridge University PressLorimer D. R., Yates J. A., Lyne A. G., Gould D. M., 1995,MNRAS, 273, 411Lyne A. G., 1998, Advances in Space Research, 21, 149Lyne A. G., Manchester R. N., Taylor J. H., 1985, MNRAS,213, 613Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005,AJ, 129, 1993Manchester R. N. et al., 2001, MNRAS, 328, 17Manchester R. N. et al., 1996, MNRAS, 279, 1235Maron O., Kijak J., Kramer M., Wielebinski R., 2000,A&A, 147, 195Narayan R., 1987, ApJ, 319, 162Perera B. B. P., McLaughlin M. A., Cordes J. M., Kerr M.,Burnett T. H., Harding A. K., 2013, ApJ, 776, 61Possenti A., Cerutti R., Colpi M., Mereghetti S., 2002,A&A, 387, 993Press W. H., Flannery B. P., Teukolsky S. A., VetterlingW. T., 1986, Numerical Recipes: The Art of ScientificComputing. Cambridge University Press, CambridgeRidley J. P., Lorimer D. R., 2010, MNRAS, 404, 1081Schnitzeler D. H. F. M., 2012, MNRAS, 427, 664Smits R., Kramer M., Stappers B., Lorimer D. R., Cordes c (cid:13) , 1–12 S. D. Bates et al.
J., Faulkner A., 2009a, A&A, 493, 1161Smits R., Lorimer D. R., Kramer M., Manchester R., Stap-pers B., Jin C. J., Nan R. D., Li D., 2009b, A&A, 505,919Taylor J. H., Manchester R. N., 1977, ApJ, 215, 885Weltevrede P., Johnston S., 2008, MNRAS, 387, 1755Yusifov I., K¨u¸c¨uk I., 2004, A&A, 422, 545 c (cid:13)000