{'fraction_non_alphanumeric': 0.039964501652260495, 'fraction_numerical': 0.01722968585405273, 'mean_word_length': 4.1113364803068295, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 27, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 123, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A nested sampling code for targeted searches for continuous gravitational waves from pulsars

This document describes a code to perform parameter estimation and model selection in targeted searches for continuous gravitational waves from known pulsars using data from ground-based gravitational wave detectors. We describe the general workings of the code and characterise it on simulated data containing both noise and simulated signals. We also show how it performs compared to a previous MCMC and grid-based approach to signal parameter estimation. Details how to run the code in a variety of cases are provided in Appendix A.


OVERVIEW
There are a variety of sources of gravitational waves that are potentially observable using the current groundbased gravitational wave detectors (Aasi et al. 2015;Acernese et al. 2015). These include transient signals, such as those already observed from the merger of binary black hole systems (Abbott et al. 2016a,b), and also continuous quasi-monochromatic signals, such as expected from assymetrically deformed rotating neutron stars (see, e.g., Section III of Abbott et al. 2004). Here we will describe a code used in searches for continuous gravitational wave signals, in particular signals from sources that are seen as pulsars via electromagnetic observations.
One of the primary methods used in searches for gravitational waves from known pulsars is based on two stages: a heterodyne stage that pre-processes the calibrated gravitational wave detector data by removing a signal's expected intrinsic phase evolution (based on the observed solution from electromagnetic observations), whilst also low-pass filtering and down-sampling (via averaging) the heterodyned data (Dupuis and Woan 2005); and, a parameter estimation stage that takes the processed data and uses it to estimate posterior probability distributions of the unknown signal parameters (e.g., using a Markov chain Monte Carlo (MCMC) Abbott et al. 2010). This method has been variously called the Time Domain Method, the Heterodyne Method, the Bayesian Method, the Glasgow Method, or some combination of these names. Up to now this method has only been used to set upper limits on signal amplitudes, but has included no explicit criteria for deciding on signal detection/significance. A further limitation has been that the MCMC method used was inefficient 4 and required a lot of tuning. It also did not straightforwardly have the ability to perform a search on the data over small ranges in matthew.pitkin@glasgow.ac.uk the phase evolution parameters, which is required if the gravitational wave signal does not well match the known pulsar's phase evolution. A nested sampling algorithm (Skilling 2006), in particular a method based on that described in Veitch and Vecchio (2010), has been provided within the LALInference software library (Veitch et al. 2015) in LIGO Scientific Collaboration (2017). This has been used for parameter estimation and model selection for both modelled transient sources, such as the observed black hole mergers (Abbott et al. 2016c), and unmodelled transient signal searches (Lynch et al. 2015). Due to the above reasons, and to make use of the existing and well used library functions, the code has been re-written to use the available nested sampling algorithm. This allows us to evaluate the evidence or marginal likelihood for the signal model and compare it to other model evidences (i.e. the data is consistent with being just Gaussian noise), whilst also giving posterior probability distributions for the signal parameters.
This code is called lalapps pulsar parameter estimation nested (or lppen for short through the rest of this document). For more detailed descriptions of how the algorithm works we refer to Veitch and Vecchio (2010) and (Veitch et al. 2015), whilst here we will provide some information on the specific proposal distributions that can be used within the algorithm. This method has previously been briefly described in Pitkin et al. (2012). The previous code, called lalapps pulsar parameter estimation (or lppe for short), used in, e.g., Abbott et al. (2010), could perform posterior evaluation using two methods: an MCMC method, or, by gridding up the parameter space and explicitly evaluating the posterior at each grid point. In §4 we show comparisons of the various methods mainly to check lppen for consistency.
The code can also be used to perform parameter estimation and model selection for signals for any generic metric theory of gravity, however its use for this is discussed in more detail in a separate paper (Isi et al. 2017). When searching over parameters that vary the phase evolution there is potential to speed up the analysis through efficient likelihood evaluation via the reduced order quadrature method (see, e.g., Field et al. 2014;Canizares et al. 2015), but again that will be discussed in more detail in a future paper.

Background knowledge
The code calculates the Bayesian evidence, or marginal likelihood, for a particular signal model under a set of assumptions. The evidence, Z, for a given model, M , defined by a set of parameters, θ, is given by Z M = p(d|M, I) = θ p(d| θ, M, I)p( θ|M, I)d θ, (1) where p(d| θ, M, I) is the likelihood function for the data d given the model and its set of defining parameters, p( θ|M, I) is the prior probability distribution on θ, and I represents any additional prior assumptions. During nested sampling this integral is evaluated by transforming it into the one dimensional sum where w i = p( θ i |M, I)∆ θ i is the "prior weight" (i.e. the fraction of the prior occupied by point i).
By default the signal model evidence is compared to the evidence that the data consists of only Gaussian noise to form the odds for the two models O = p(d|M, I) p(d|noise, I) where on the right hand side we have explicitly set the prior odds for the two models to p(M |I)/p(noise|I) = 1. Other than this evidence value and odds, the code also returns the samples accumulated during the nested sampling process. These samples are not drawn from the posterior probability distribution as in an MCMC method, but they can be resampled to generate a subset of samples that are drawn from the posterior distribution. This resampling (performed using the lalapps nest2pos python script within lalapps LIGO Scientific Collaboration 2017) uses the value of L i w i for each point, normalised by (L i w i ) max , to give values proportional to the posterior probability, and then accepts a point with a probability given by its value, i.e. a point is accepted with the probability In our case the data d in equation 3 is the heterodyned, low-pass filtered and down-sampled data from a detector, or set of detectors, each giving a single data stream or two data streams depending on whether the heterodyne was performed at one or both potential signal frequencies near the rotation frequency and/or twice the rotation frequency. Here we will use B to represent a vector of these heterodyned data values (Dupuis and Woan 2005), for which a single time sample is often referred to using the jargon term B k ≡ B(t k ) ("B of k"), although we will not consistently use k as the time index for the data. Throughout this document if we refer to "heterodyned" data we mean data that has been processed, or simulated as if processed, by this method.

CORE FUNCTIONS
Here we will describe the core parts of the code defining the signal model and the probability functions used for the inference. These assume that the data comprises a number of calibrated narrow-band complex-heteroyned time series. These time series data streams can be from different detectors, and/or at different heterodyne frequencies. For example you could have a pair of times series from both the LIGO Hanford (H1) and LIGO Livingston (L1) detectors, with each pair containing one produced with a heterodyne at twice the rotation frequency of a known pulsar and the other produced with a heterodyne at the rotation frequency.

The signal model
Our code assumes that the electromagnetically observed rotational phase evolution of a pulsar is represented by its Taylor expansion where T is the time, corrected to an inertial reference frame (the solar system barycentre for isolated pulsars, or the binary system barycentre for pulsars in binaries or multiple component systems), and the f 's give the observed rotation frequency and its time derivatives. The value of T = (t + τ (t) − t 0 ) where t is the time at a detector, τ (t) is the time dependent correction to the inertial frame and t 0 is the epoch. The code can accept frequency derivatives of any order. We assume the calibrated detector data, d(t), has been heterodyned such that where Φ i,het (t) is the assumed phase evolution for a given data stream i, and we produce B(t) by low-pass filtering and resampling (via averaging) values of B (t).
Under the standard assumption that the general theory of relativity (GR) is correct, the code uses the most general form of the signal model defined in Equations 51-54 of Jones (2015), which, when heterodyned and assuming low-pass filtering, gives a signal at a pulsar's rotation frequency (where Φ 1,het (t) = φ(t)) of and at twice the pulsar's rotation frequency (where Φ 2,het (t) = 2φ(t)) of The F D + and F D × values are the detector (D) dependent antenna patterns (see, e.g., Forward (1978); Schutz and Tinto (1987); Finn (2009), or Jaranowski et al. (1998 for the convention used within this document in, e.g., equation 16), which are a function of the detector position, source sky position and source polarisation angle ψ (in equations 7 and 8 we only explicitly note the ψ dependence as for our sources the sky location is known). The C 21 , C 22 , Φ C 21 and Φ C 22 values are convenient ways of representing the waveform in terms of an amplitude and phase of the signal for the l = 2, m = 1 harmonic and l = m = 2 harmonic respectively. The ∆φ(t) values represent any time dependent phase deviation between the phase used in the heterodyne and the true signal phase (which does not necessarily have to precisely follow the electromagnetically observed rotational phase, see discussions in, e.g., Abbott et al. 2008), so ∆φ 1 (t) = (φ 1,true (t) − Φ 1,het (t)) and ∆φ 2 (t) = (φ 2,true (t) − Φ 2,het (t)). More generally, for emission at an arbitrary scaling of the rotation frequency, K (where Φ K,het (t) = Kφ(t)), an even more generic signal model (still assuming GR is correct) would be where H + and H × are the amplitude components for the '+' and '×' polarisations, and Φ 0 is the signal's initial phase.
To calculate the ∆φ K values using up to the (n − 1) th frequency derivative, and to avoid numerical overflow issues when dealing with large phases, we use where f (n) is the n th frequency derivative (for both the 'true' and 'het' the frequency/frequency dervitative values used here are those prior to scaling by K for a given signal model), and δt is the combination of any solar system barycentring and binary system barycentring time delays. By default the code assumes emission just from the l = m = 2 mode, i.e. there is only a signal at twice the rotation frequency. In this case C 22 and Φ C 22 can be related to the more familiar physical h 0 and φ 0 values via h 0 = −2C 22 (where the minus sign maintains consistency of equation 8 with the form given in Jaranowski et al. 1998) and pulsar rotational phase φ 0 = Φ C 22 /2. For the more general case the relations between the waveform amplitude and phase parameters and physical source parameters are given in Jones (2015) and . In general, for previous searches we have often assumed that we track the true signal phase perfectly with the heterodyne, and as such ∆φ K (t) = 0. In such cases the only time varying components of the signal are the antenna pattern functions, which allows great speed increases in the signal generation and likelihood calculations (see §2.2.4).

General tensor-vector-scalar signal model
A more generic heterodyned model, assuming gravity is described by an arbitrary metric theory, but not necessarily GR, is given by (see, e.g. Isi et al. 2015Isi et al. , 2017 h tvs where H p (for p ∈ {+, ×, x, y, b, l}) are the amplitudes for the tensor ('+' and '×'), vector ('x' and 'y'), and scalar ('b' and 'l') gravitational wave amplitude components, F p are their associated time dependent antenna responses (see, e.g. Nishizawa et al. 2009;B laut 2012;Isi et al. 2015), and Φ t,v,s and ψ t,v,s are angles related to the initial phase of, and rotation within, the tensor, vector and scalar components. For networks of quadrupolar antennas, like a standard gravitational wave interferometer, the two scalar components are entirely degenerate , and as such only one of them needs to be defined (e.g., only the 'b' mode is used in Equation 19 of Isi et al. 2017). When searching for a generic signal such as this the orientation of the wave-frame (given by ψ and ι) used to define the polarisation modes is arbitrary: one is free to pick any right-handed frame one wants and define the modes there, but we usually make a choice that we know simplifies things for a given sig-nal. Hence, for a generic search, we are free to arbitrarily set our frame such that ψ = 0 and the values F p given in equation 11 can be defined by, e.g., equations 24-29 of Isi et al. (2015). However, when we have a specific prediction of what a given waveform should look like, it is always in a particular frame. By convention, this is usually a frame in which the x-axis is aligned with the intersection of the equatorial plane of the source with the plane of the sky (though the convention is not universal). In other words, if a theory predicts certain amounts of polarisation with some specific phase evolution (say, that given by GR) then that is a frame dependent statement, and ψ in required to orient the model frame accordingly.
In such a case (e.g., GR or some other theory such as G4v in Mead 2015) the antenna patterns could be obtain from those with ψ = 0 using equations 30-35 (and ψ t and ψ v from equations 46-47, and 49-50, for GR and G4v respectively) of Isi et al. (2015) would be required.
2.2. The likelihood functions Our code can make use of two different likelihood functions. The default is a Student's t-likelihood function in which it is assumed that the standard deviation of the noise in the data is unknown, but can be marginalised over. A Gaussian likelihood function can also be used, for which the code can either take in estimates of the noise standard deviation at each data point, or calculate these internally based on stationary stretches of data. For the Student's t-likelihood function, and if calculating noise standard deviations for the Gaussian likelihood function internally, the code needs to break up the data into segments for which the noise can be assumed to have been drawn from independent Gaussian distributions. The method for doing this is given in §2.4.

Student's t-likelihood
A full derivation of the Student's t-likelihood function (see, e.g., Dupuis and Woan 2005) is given in Appendix B, but the final form of the joint likelihood (and its natural logarithm, which is actually used within the code to maintain precision) for multiple detectors and data streams is given by where N dets is the number of detectors used, N s is the number of data streams (e.g., heterodyned data from both the rotation frequency and twice the rotation frequency) per detector, M i,j is the total number of independent data chunks for detector i and data stream j with lengths m i,j,k , and n i,j,0 = k l=1 1 + m i,j,l−1 (with m i,j,0 = 0) being the index of the first data point in each chunk. The signal model y( θ) is that given by Eqns. 7 and/or 8 (or Eqn. 9) depending on which data streams are being analysed. For notational convenience we have made the substitution

Gaussian likelihood
The Gaussian likelihood, and its natural logarithm, are similarly given by where L i,j is the length of each dataset and B i,j,k = − ln (2πσ 2 i,j,k ). Note that normalisation factor reflects the fact that exponential is already the product of the real and imaginary data components.

The null likelihood
As we will most often want to perform model comparison for the signal model against one for which the data just contains noise (the null hypothesis in this case) we can define the null likelihoods for both of the above likelihoods by setting y( θ) = 0, giving for the Student's t-likelihood, and ln p( for and Gaussian likelihood. If we were only interested in comparing models calculated using equivalent likelihood functions we could in general ignore the factors that do not depend on the data or the model, as they would cancel in any odds ratio. But, in this code we keep them for cases when such a comparison is not performed, e.g., if we want to compare the joint multi-detector likelihood for a signal with the incoherent product of likelihoods from each detector then we would need these factors to still be present.

Fast likelihood evaluations
In cases when the only time varying components of the model are the antenna pattern functions (i.e. when ∆φ j (t) in equation 7 or 8 is zero) the likelihood evaluation can be greatly sped-up by pre-calculating the components in the internal summations. For a given sky position and detector, D, and using the standard reference frame, the antenna patterns can be written as (Jaranowski et al. 1998, following the form of) where ψ is the gravitational wave polarisation angle, ζ is the known angle between the detector arms (generally 90 • ), and a(t) and b(t) are the time dependent functions for a given detector position and source sky location that vary over a sidereal day. 5 These functions can be precomputed at a set of times over a sidereal day and used, via look-up table interpolation, to give the value at any other time (our code defaults to calculate a(t) and b(t) at 2880 points over a sidereal day). However, the vast majority of the speed-up comes from pre-computing summations over the combinations of the data and the a(t) and b(t) functions. This is described in more detail in Appendix C.
In the case where we want to search over phase evolution parameters, meaning that ∆φ K (t) = 0 in the model, then we cannot perform this pre-summing. However, reduced order quadrature methods (e.g. Field et al. 2014;Canizares et al. 2015) may help in these cases. Such a method is implemented in our code, but will be discussed in a separate paper.

The prior functions
An important part of any Bayesian inference method is the choice of parameter prior probability functions. Here we describe the prior functions allowed in the code for any of the required parameters.
Certain parameters are physically not allowed to take negative values, so it is hardcoded that the following parameters have zero probability below zero: gravitational wave amplitude, h 0 ; l = m = 2 mass quadrupole moment, Q 22 ; the pulsar distance; the pulsar parallax; the speed of gravitational waves; the projected semi-major axis of a binary orbit; the total binary system mass; and, the companion mass in a binary system.
There are also other hardcoded priors for waveform amplitudes in two particular cases when searching for 5 The functions a(t) and b(t) are just the antenna patterns in a frame defined with ψ = 0. The actual definition of the antenna pattern functions is given in, e.g., Equation 23 of Isi et al. (2015). emission at both once and twice the rotation frequency: if searching for a signal from a biaxial star then the two waveform amplitudes, C 21 and C 22 , must either both be positive or both be negative; if using the "source" model then the two amplitudes I 21 and I 31 , must have I 31 I 21 (see, e.g., .

Uniform prior
A parameter can be given a uniform (or flat, or tophat) prior, in which the probability is constant within a given range and zero outside that range, i.e.
where for parameter x the upper and lower bounds on the parameter values are x min and x max . This prior can be normalised by setting C = 1/(x max − x min ), although within nested sampling the normalisation is not specifically required. The prior requires a lower and upper bound to allow it to be normalisable and to enable an initial set of samples to be drawn from it. However, it should be noted that the choice of range will have an effect on evidences that are calculated as it affects the volume of parameter space allowed by the model. A histogram of a set of samples produced by the code when drawing from a uniform prior distribution are shown in Figure 1, along with the true prior function.
Some parameters for which the uniform prior is often considered appropriate are angle parameters, such as the Φ C 22 and ψ parameters in Equation 8, which can be restricted to a specific non-degenerate range (see, e.g., Table 1 in . In previous gravitational wave pulsar searches, such as Abbott et al. (2010); Aasi et al. (2014) a uniform prior was often used for the gravitational wave amplitude parameter h 0 . For those searches the prior upper bound was based on inspection of the data to find a value that was very large compared to the sensitivity suggested by the data. A uniform prior was used, as opposed to the more uninformative choice for such a scale parameter of the log-uniform prior (see §2.3.4), due to the fact that the main result was in producing upper limits, and wanting these to be trivially related to the likelihood, rather than being heavily influenced by the prior.

Gaussian prior
A parameter can be given a prior with a Gaussian distribution, defined by its mean, µ x , and standard deviation, σ x , i.e.
A histogram of a set of samples produced by the code when drawing from a Gaussian prior distribution are shown in Figure 2, along with the true prior function.
As well as single parameters being given Gaussian priors, we can define a set of k parameters, x, to have a multi-variate Gaussian prior where µ x are the means and C is the covariance matrix. In reality we input the correlation coefficient matrix and individual parameter standard deviations rather than the covariance matrix, but it is easy to convert between the two.

Gaussian Mixture Model prior
A parameter, or any set of parameters, can be given a prior distribution composed of a superposition of Gaussian distributions (a Gaussian Mixture Model), with each component defined by a set of means for each parameter, a covariance matrix, and a weight describing the relative probability concentrated in it. For an n component mixture on a set of k parameters, x, the prior function will be where, for component i, w i is the weight, ∆ x i = x − µ xi , µ xi are the means, and C i is the covariance matrix. Such a model can be used, for example, if there is a bimodal Gaussian prior on a particular parameter, as is the case for the inclination angle ι if it is calculated from fits to pulsar wind nebulae (see, e.g., Appendix B in Abbott et al. 2017). In such a case the Gaussian Mixture Model can define two Gaussian distributions of equal standard deviations and weights. A more complex use of this prior would be for reconstructing a smooth probability distribution using samples drawn from it, i.e. taking posterior samples from a search for a particular pulsar, and reconstructing a smooth model of those samples (for, say, the two parameters h 0 and cos ι) to then be input as a prior when looking at future data for the same source.
Two examples of histograms of samples drawn from one-dimensional Gaussian Mixture Models with two and three modes respectively are shown in Figure 3.

Log-Uniform prior
A parameter can be given a prior distribution that is uniform in the logarithm of the parameter, i.e.
This prior requires minimum and maximum values to be specified to make it normalisable and allow samples to be drawn from the distribution. It should be noted that the choice of range, and in particular the lower bound, will have some effect on the calculated evidence (see, e.g., Appendix B of Isi et al. 2017).
A histogram of a set of samples produced by the code when drawing from a log-uniform prior distribution, with a range between 10 −3 and 1 are shown in Figure 4, along with the true prior function.

Fermi-Dirac prior
A Fermi-Dirac distribution prior function was inspired by that used in Middleton et al. (2015). The prior has a sigmoid, or logistic, shape, although starts at large values and slopes downwards and is restricted to positive values. The form of the function, for parameter x, is and is defined by two parameters σ and µ. µ defines the point at which the distribution falls to 50% of its maximum value, whilst σ defines the range over which the attenuation happens. We also make use of a parameter r = µ/σ. The band over which the probability falls from 97.5% to 2.5% of its maximum value is given by µ ± Zµ/2r, where Z ≈ 7.33. Therefore, if we want, for example, the distribution to have this particular fall-off over a range that is 10% of µ, we would have Z/2r = 0.1 and σ = 0.027µ.
Three different examples of Fermi-Dirac priors, and samples drawn from them by the code, are shown in Figure 5, for a selection of µ and σ values.
The sampling from the prior uses inverse transform sampling. The cumulative distribution function (CDF) of Equation 22 is which can be inverted to give Thus, drawing CDF samples uniformly between 0 and 1, and inverting via the above equation gives values drawn from Equation 22. This prior can use useful for amplitude parameters, where a uniform prior may have previously been used. It has the property of being roughly flat at low values of the parameter and with an exponential fall-off at large values. This gives the advantage over a uniform prior with a sharp upper cut-off in that the distribution is continuous rather than (perhaps) arbitrarily truncated.
This form of prior was used for the gravitational wave amplitude parameter, h 0 , in Abbott et al. (2017) as opposed to the earlier use of a uniform prior with arbitrary upper cut-off used in, e.g., Aasi et al. (2014).

Splitting the data
The Student's t-likelihood function ( §2.2.1) requires us to have an idea about the timescales on which the data is stationary, i.e. periods when the noise in the data is best described as being drawn from a single Gaussian distribution. We therefore use a scheme similar to the BayesianBlocks change point algorithm (Scargle 1998), or BlockNormal gravitational wave burst finding algorithm (McNabb et al. 2004), to find points at which the statistics of the noise changes. This is done with a top-down iterative divide and conquer approach (Scargle 2000).
We start by taking a full complex time series data set for a detector and subtracting a running median from both the real and imaginary components. This subtraction is performed to try and remove the effect of any very strong signals in the data, as we want to assess the properties of the noise rather than any signal. 6 The running median is calculated over a running window of 30 data points (or a minimum of 15 points at the start or end of the data), but not accounting for gaps in the data. For the standard data sample rate of 1 per 60 seconds, this means that half an hour of data is used, over which time the signal (through the varying antenna patterns) will not change very much (unless there is a gap in the data). However, it should be noted that, if using a much more slowly sampled data set, then this hardcoded 30 sample window may not work well for very strong signals.
The next thing we do is calculate the evidence that the whole (running-median-removed) dataset is drawn from a single Gaussian distribution with unknown variance. To do this we just use Equation 14, with N dets = 1, N s = 1, M 1,1 = 1, n 1,1,0 = 1 and, given a dataset of length N , set m 1,1,1 = N . We will call the natural logarithm of this value ln Z single . We want to calculate the odds between this evidence and that for the data containing any single change point (i.e. point at which the data seems to be drawn from a different Gaussian distribution). So, for one single change point, at index i, we calculate ln Z cp,i , by using Equation 14 and setting N dets = 1, N s = 1, M 1,1 = 2, n 1,1 = {1, i} and m 1,1 = {i, N − i}. We also impose a minimum allowed data chunk size of n min (which defaults to a value of 5 data points in the code), so that i n min . As we want the evidence for any single change point, we must take the sum of all the single change point evidences ln Z any-cp = ln Finally, we are left with the odds (assuming equal prior probabilities for each hypothesis) ln O any-cp single = ln Z any-cp − ln Z single , for which values above zero favour there being a change point in the data. We could simply use this as the criterion on which to split the data into two chunks at the change point value i with the largest evidence. However, in practice this leads to data drawn from a single Gaussian distribution actually being split far too often, and also there being a dependence on the data length. So, instead, we empirically calculate an odds threshold above which the value of Equation 26 must be to impose a split in the data. 7 The empirical threshold we use is worked out by finding the 99% upper limit on the value of Equation 26 when data is purely drawn from a single Gaussian distribution (i.e. it gives a 1% false alarm probability) as a function of the data length N . The results of doing this for 30 000 instances of data for a range of n values, and three different false alarm probabilities (1%, 0.5% and 0.1%), are shown in Figure 6. Fitting a line to the 1% false alarm probability points gives an odds ratio threshold dependence of T = 4.07 + 1.33 log 10 N. The above process just splits the data once, so it must be applied iteratively to split the data further. When a change point is found then the process is repeated on the two split data chunks, and stops when either no change point is found within a chunk, or the split would give a segment smaller than n min .
2.5. Signal-to-noise ratio calculation Whether a signal is found or not, the code will generate a signal-to-noise ratio (SNR). The calculated SNR uses the set of signal parameters, θ ML , that give the maximum value for the likelihood, which are then used to create the best-fit complex signal template y( θ ML ). We then use calculations of the noise standard deviation for each data chunk after removal of a running median from the data (the splitting and running median are described in §2.4). The coherent SNR is then calculated as (28) using the notation given in §2.2. 8 A plot of the residual between the recovered singledetector and joint detector SNRs and their true values 9 versus the true signal SNRs, as calculated from Equation 28, for a set of simulated signals (described in §4.1.2) are shown in Figure 7. The standard deviation of the residuals is found to be one, with a mean around zero, and the distribution does not appear to change as a function of SNR (except at the lowest SNRs, where the noise in the data will generally give rise to a maximum likelihood model with a larger amplitude than the true signal, and thus an overestimate of the recovered SNR).

Odds values
This is not explicitly a core function of the code, but here we will define some model selections (odds between competing models) that can be produced from the output of the code. The code outputs two (natural logarithm) evidence values: the evidence for a coherent signal (Z S ) in all the datastreams (multi-detector and multi-frequency-band), and the evidence (Z noise ) for the data consisting of Gaussian noise (in sets of stationary segments as described in §2.4). These can be combined in various ways for different model selection situations, and we will define three 8 If using data input from the spectral interpolation code (Davies et al. 2017) the noise standard deviations are actually passed to the code (see §A.1), so in Equation 28 M i,j is replaced by L i,j and y( θ ML ) i,j,n is replace by y( θ ML ) i,j,k as in §2.2.2. 9 The true values are calculated by using the actual simulated signal model in equation 28, rather than the recovered maximum likelihood model. such situations here (see, e.g., Abbott et al. 2017, for two of these): 1. the odds for a (multi-detector and multi-frequencyband) coherent signal versus (non-stationary) Gaussian noise, O S/N ; 2. the odds for a (multi-detector and multi-frequencyband) coherent signal versus incoherent (between detector) signals, O S/I simple ; 3. the odds for a (multi-detector and multi-frequencyband) coherent signal versus combinations of incoherent (between detector) signals or noise, O S/I (cf. Equation A6 in Abbott et al. 2017).
The first of these (actually also given in Equation 3) just requires the two evidence values, Z S and Z noise The second requires that the code has been run for a multi-detector analysis, and also individually for each detector, i, to give Z Si and Z Ni , and is given by Such an odds is useful in vetoing strong incoherent signals that are therefore not astrophysical, e.g. instrumental lines (see, e.g., the similar line-robust statistic defined in Keitel et al. 2014). In both these cases it is assumed that the prior odds (p(M 1 |I)/p(M 2 |I)) between the two models are set to be equal and therefore are not explicitly given. The third odds makes use of a compound model in the denominator with all evidences for all permutations of Gaussian noise or incoherent signal being used. If we have the set of single detector signal and Gaussian noise models M j = {S j , N j }, then (see, e.g., Equation 50 of Isi et al. 2017) where p(S|I) is the prior on the coherent signal model, and p(S/N j |I) are the priors for the signal/noise hypotheses for each individual detector. As there are two possible models for data in each detector, signal or noise, the number of logical or combined sub-hypotheses (see, e.g., the four sub-hypotheses used in Equation A6 of Abbott et al. 2017, where the priors in the denominator are the combinations that would be given from expanding out equation 32) in the compound model is just given by 2 N dets . In Abbott et al. (2017), and later in this document, we have assigned the equal probabilities to the coherent signal prior and to each of the sub-hypothesis priors in the denominator, so that they cancel out. However, other priors could be used.
It is also worth noting that in all the above cases the models (and sub-hypotheses) are mutually exclusive. This may not at first seem obvious, as, for example, the Gaussian noise model is contained within the signal model provided the prior on the signal amplitude allowed a value of zero. But, as noted in Li et al. (2012); Isi et al. (2017), within the signal noise hypothesis having an amplitude of exactly zero (as in the Gaussian noise hypothesis) occupies an infinitesimally thin slice of the parameter space, and therefore has no weight. The same can be said for the situation of a coherent signal model within the incoherent signal model.

NESTED SAMPLING SETUP
Here we will briefly describe some of the core functions used for the Nested Sampling algorithm (Skilling 2006) within our code. We will not describe nested sampling itself, over-and-above the very brief summary given in §1.1, but instead refer the reader to Veitch and Vecchio (2010). We also refer the reader to Veitch et al. (2015) for a more thorough description of some of the concepts we will discuss below.
One main point about the nested sampling algorithm is that it starts by randomly drawing a certain number of samples (the number of which we call the number of live points, N live , which will be the constant number of "active", or "live", samples), where a sample is a vector, θ, containing a value for each parameter drawn from their respective prior distributions (see §2.3). The samples must be independent draws. For each of these samples the likelihood will be calculated and the sample θ min corresponding to the minimum likelihood, L min = p(B| θ min ), will be noted. As the algorithm progresses the number of live points stays the same, but at each iteration a new sample will be added in place of the one with the minimum likelihood, and therefore the minimum likelihood will continuously be updated to be the minimum value for the current set of points.

Proposal functions
The nested sampling algorithm requires a way to draw new samples from the parameter space that is being explored. The particular requirement for nested sampling is that any new sample is an independent draw from a constrained version of prior probabilty function for each parameter. The constraint is that the new sample's likelihood must be greater than the current value of L min , e.g. for a single parameter x, the probability of drawing a point at x i would be where p(x|I) is the prior on x.
There are various methods to sample from the constrained prior, the simplest of which is just to sample from the full prior and only accept a sample if it meets the likelihood constraint. However, such a simple method soon becomes very inefficient if the likelihood gets more and more tightly concentrated within the prior volume (which can happen rapidly for high-dimensional problems). One popular nested sampling method is Multi-Nest (Feroz et al. 2009), which performs the sampling by calculating a set of ellipsoids bounding the live points and drawing points uniformly from within them. Various other approaches exist such as Diffusive Nested Sampling (Brewer et al. 2011;Brewer and Foreman-Mackey 2016) and POLYCHORD (Handley et al. 2015).
The version of nested sampling used in the LALInference library, as used by our code, makes use of a Markov chain Monte Carlo (MCMC) approach to drawing new samples (Veitch and Vecchio 2010). The MCMC attempts to sample from the constrained prior distributions in an efficient manner, with the final sample of the chain being a statistically independent sample from the rest of the live points. To sample most efficiently the MCMC requires a proposal distribution that closely resembles the underlying constrained prior that it is trying to sample from. LALInference has a variety of potential proposal distributions (Veitch et al. 2015), but our code currently allows five different proposals, that can be used in different proportions: Ensemble walk: This is one of the affine invariant ensemble proposals described by Goodman and Weare (2010), which adopt the scale and shape of the current set of samples. The hard-coded setup in LALInference is that three samples randomly drawn from a cache of samples are used to define the "walker" samples. This is the main default proposal of our code, and is used as the proposal 75% of the time.
Ensemble stretch: This is another affine invariant ensemble proposal described by Goodman and Weare (2010), and makes use of two samples randomly drawn from the sample cache. By default this proposal is not used at all.
Differential evolution: This is similar to the stretch proposal in that it uses two samples randomly drawn from the sample cache. By default this proposal is not used at all.
Uniform proposal: This proposal just draws points from their full prior for parameters that have a uniform prior specified. Other parameters are not evolved. This proposal is used 25% of the time by default.
Frequency jump: This proposal tries jumping between Fourier frequency bins if frequency is one of the required parameters. This proposal has not been tested and, by default, is not used at all.
The cache of samples mentioned above, and used by the nested sampling algorithm, is the set of live points, but this cache is only updated on iterations of the algorithm that are a factor of 10% of the number of live points (e.g., for 1000 live points the cache is updated after every 100 new iterations). The length of each MCMC run can be chosen by the user, but by default it is automatically set within the code. To do this, at the same rate as the cache is updated, the autocorrelation length of the parameters is calculated. A random sample from the current live points in chosen, and evolved via the MCMC, and the autocorrelation length of each parameter, acl, is calculated.
The number of MCMC iterations, n MCMC , is then chosen as n MCMC = min[max( acl), 5000], where 5000 is a hardcoded maximum length.

Testing the proposals
We would like to test whether some of the above proposals lead to accurate calculations of the evidence integral as given in Equation 1. It is especially interesting to do this for a variety of prior volumes, starting of with volumes similar to the posterior volume, but then diverging to be much larger than the posterior volume (or in other words, as the information gain, H, or Kullback-Leibler divergence, increases). To do this our code has the option of setting a one-dimensional Gaussian distribution as the likelihood function, using a flat prior bounded between A and B, and working out the evidence. Analytically the evidence is given by As we are often interested in setting upper limits of, for example, the gravitational wave amplitude h 0 , we can also use the test to see if the upper limits we obtain from the posterior samples match the expected upper limit obtained (via root finding) from the CDF of the distribution where p(x|µ, σ, I) is the posterior, p(µ, σ|x, I) the likelihood, and p(x|I) is the prior (with p(x|I) = (B − A) −1 ), and Z is the evidence calculated via Equation 34. For the case of µ = 0 and A = 0, and provided B σ, this generally simplifies to A final thing that we can check, on top of the upper limit value, is whether the posterior samples appear to be drawn from the expected posterior distribution. To do this we use a Kolmogorov-Smirnov test comparing the CDF of the samples with the analytical CDF and calculating the two-sided p-value for the null hypothesis that the two distributions are the same.
In the tests below we have run the code on this Gaussian likelihood with a Gaussian of zero mean and σ = 10 −24 , truncated at zero by fixing the lower prior bound to A = 0. This is similar to the form of the likelihood seen in pulsar searches when no signal is present. We have then set a range of upper bounds, B, from 10 −23 to 10 −13 increasing by a factor of ten each time, and a range of numbers of live points from 2 9 to 2 13 increasing by powers of two each time. For each combination of upper bound and number of live points we have run the code on the test Gaussian likelihood 100 times to see the statistical spread of the output.
In the first test we check how using only the ensemble walk proposal (which is not the default configuration described above) fairs at estimating the evidence, upper limit, and posterior distribution. Figure 8 shows a variety of information, but it primarily shows the distribution of the ratio of evidences output from our code, Z, to the true evidence calculated via Equation 34 as a function of the value of the upper prior bound B (or equivalently on the top axis the information gain, in natural units of information (nats), calculated via Equation 36). This is shown for the range of different numbers of live points used. It is clear that as the information gain increases (i.e. the posterior is constrained within a successively smaller part of the prior range) we are systematically more biased towards overestimating the evidence. It can also be seen that the bias increases with increasing number of live points. Linear fits (in ln (Z/Z true ) and H, the information gain) for different numbers of live points are shown in Table 1. It is interesting to note that the statistical spread of values (if ignoring the systematic offset) are consistent with the expected range given by the error bars, and calculated using H /N live (Skilling 2006), where here H is the mean of the information gain (in nats) returned by the nested sampling algorithm for each run.     It can be seen that, other than for the narrowest prior range, there is a systematic bias towards the upper limit being underestimated by ∼ 2%. Figure 10 show the distribution of two-sided pvalues from a Kolmogorov-Smirnov test comparing the CDF of posterior samples against the analytical CDF in Equation 35. The p-value null hypothesis is that both distributions are the same, and it should be noted that we performed 100 trials for each prior range value and number of live points. For the lower numbers of live points the distributions appear to be consistent, however, for larger numbers of live points the distributions are rather inconsistent a large amount of the time. This is most likely due to differences in distributions showing up more obviously when there are more posterior samples to use.
The systematic biases in evidence values and upper limits may not be a major problem 10 , but does point towards there being some issue with the implementation of the ensemble walk proposal.
In the second test we have checked what happens with the code's current default proposal settings: using the ensemble walk proposal 75% of the time and the uniform proposal 25% of the time. It should be noted that this default is not arbitrary, but came about from tests similar to these when just using the ensemble walk proposal. Figure 11 shows the evidence distributions equivalent to those in Figure 8. However, what can be seen in this case is that there is no obvious systematic bias as a function of H, and the statistical uncertainty is consistent with expectations (as shown by the error bars). Similarly, the upper limits and Kolmogorov-Smirnov p-values shown in Figures 12 and 13 respectively are consistent. The reason why this combination of proposals seems to work whereas the ensemble walk proposal on its own does not is currently not understood, but we plan more work on figuring this out.
10 The evidence biases at the level seen are unlikely to sway model comparisons as they will most likely show up when one model is very highly favoured over another. The observed biases in the upper limits are at a level that is currently less than expected uncertainties due to instrumental calibration.   As discussed in §1.1 nested sampling codes, such as lppen, do not intrinsically output samples that can be histogrammed to give representations of marginalised posterior distributions for the model parameters. They instead output an ascending likelihood ordered list of samples, each of which occupies a known approximate amount prior volume. Samples must be drawn from these, with appropriate weighting (see Equation 4), to give a new set of samples that does represent the posterior distribution. During this process it is also possible to appropriately combine samples from parallel independent nested sampling runs (weighting each based on their calculated evidences) to recalculate the evidence and generate posterior samples.

Posterior sample generation
What is described above is not a core function of the lppen code, but is an integral part to generating results with it, and has been used to produce the posterior distributions discussed earlier in the section, and later in this document. The number of posterior samples generated by the code is dependent on the number of live points used when running nested sampling. The actual dependence of the number of posterior samples on N live can be seen in Figure 14 (this uses the simulations discussed in §4.1.3). It can be seen empirically, that the mean number of posterior samples goes roughly as N post ≈ 2.4N live , although there is a lower tail extending to N post ≈ 0.01N 1.14 live .
3.3. Code timing How long the code takes to runs depends on two main factors: the number of live points used and the information gain going from prior to posterior. Increases in both of these will generally mean that the code takes longer to run. Using the simulations from §3.1.1 we have assessed the code run time for the simple testing Gaussian likelihood function. This allows us to quantify the code run time as a function of both number of live points and information gain.
We find that the median calculation time for the simple Gaussian likelihood function used in §3.1.1 (on an Intel "Haswell" processor, specifically Intel CPU E5-2680 v3 @ 2.50GHz), including some standard overheads 11 is ∼ 5.9 × 10 −6 seconds. The median time for calculation of the likelihood (including model calculation time) for a standard search (Eqn. 12) over only non-phase evolution parameters, and for one chunk of data, is 8.1×10 −5 seconds. 12 The median run time per likelihood evaluation for the nested sampling part of the code for the standard likelihood function (using 1024 live points) is 8.5×10 −4 seconds, which is ∼ 10 times greater than the single likelihood evaluation time (this ratio will change as a function of the number of live points). The median run time per likelihood evaluation for the nested sampling part of the code for the test Gaussian likelihood function (using 1024 live points) is 1.8×10 −4 seconds, which is ∼ 30 times greater than the single likelihood evaluation time.
In Figure 15 we show the time taken for the runs in §3.1.1, using the default proposals, as a function of number of live points and information gain. The code run time is divided by T L = 5.9×10 −6 to provide a way to scale it for other values of T L , e.g. that found from the code's standard likelihood function. If we perform a 2D linear fit to the natural logarithm of the median run time, R/T L , we get the relation This relation should be considered as a rough lower limit to the run time for any particular likelihood function.
Above, we saw that for our standard likelihood function the nested sampling algorithm was roughly three times quicker per likelihood evaluation than the simple test function. In both cases the internal MCMC, which draws new samples within the algorithm, required very similar small lengths to give uncorrelated samples, so there must be other internal bottlenecks that make the simpler likelihood relatively slower. However, in general, longer MCMC chains will be required to give uncorrelated samples, so the number of likelihood evaluations will increase.

EVALUATING THE CODE
Here we will evaluate the code in a variety of ways. We compare a selection of posterior probability distributions generated using lppen, over a range of signal parameters, with those produced using a previous code. We also use a selection of many simulated signals to show how the odds comparing a signal model and a noise only model change with signal-to-noise ratio, along with comparing coherent verses incoherent signals (for coherent simulated signal and incoherent simulations). The simulations also allow us to evaluate the posterior probability distributions produced by the code using what is colloquially known as "P-P plots" (see §4.3.1)

Code comparison
In this section we will show how the posterior parameter distributions that we can generate from the output of our nested sampling code compare to evaluating them with a previous implementation of the analysis code (lalapps pulsar parameter estimation, or for the rest of this document shortened to lppe, used in, e.g., Aasi et al. 2014) using a (generally inefficient) MCMC algorithm and over a grid in parameters. 13 We will show comparisons of the 13 For the MCMC examples in this section, we have used "burn-in" periods of 100 000 iterations, followed by 500 000 pos-codes in two different regimes: no signal is present in the data, and an obvious signal is present in the data. We do this using simulated data both containing purely Gaussian noise and containing simulated signals added to Gaussian noise.
In each of the tests below we will be running our code (lppen) using the default proposal distributions discussed in §3.1 and with the number of live points fixed to be 2048. We will also assume a source at a fixed sky position and that the signal model is purely that for the l = m = 2 harmonic (Equation 8), and that we work in terms of the signal amplitude h 0 and signal phase Φ C

22
(this is for consistency with the old code which uses Φ C 22 rather than φ 0 , which we have earlier defined as the rotational phase). The simulated data in all cases will be generated as the complex heterodyned time series sampled once per minute. When comparisons between the current code and previous codes are made we run the current code in a way that splits datasets into 30 point chunks, rather than using the algorithm described in §2.4, as this makes it consistent with the old code (lppe). We also purely use the uniform prior ranges for the amplitude parameter to be consistent between codes.

Simulated noise
An initial test of the code is whether the output posterior probability distributions match those produced when evaluating the posterior over a fixed grid in the fourdimensional parameter space of θ = {h 0 , cos ι, ψ, φ C 22 }, when using purely Gaussian noise. 14 A comparison of the marginalised posteriors for individual parameters, and pairs of parameters, are shown when using simulated data lasting ten days from a single detector (in this case assumed to be the LIGO Hanford detector, H1) in Figure 16.
It can be seen that the output posteriors look qualitatively consistent. However, we can also quantify some aspects of consistency by comparing the evidence output from the nested sampling code with that estimated from the grid-based method from lppe, the upper limits on h 0 produced by the codes, and performing Kolmogorov-Smirnov consistency tests between the nested sampling posterior samples and MCMC posterior samples (giving a p-value for the null hypothesis that the samples are from the same distributions). These are shown in Table 2 and show very good consistency between the codes, although it should be noted that these are for one particular run and some statistical fluctuations in the exact values for different runs are to be expected (see, e.g., the variations in evidence values in §3.1.1). A jupyter notebook with this test can be found here.
We see a very similar situation, in terms of agreement between the codes, when running on simulated data assumed to be from two detectors (the LIGO H1 and L1 terior sample iterations, from which there have been roughly 1000 independent samples used to produce the posterior plots. For the grid-based examples the grid sampling used has been 14 It is worth noting that in many of these tests the polarisation angle prior covers −π/4 to π/4 rather than the 0 to π/2 range given in, e.g., , due to this being the range required by the older parameter estimation MCMC code. Although this range is degenerate to rotations over π/2, so spans the same range of waveform models.

TABLE 2
Consistency tests between outputs of the new code, lppen, and the old code, lppe, when running on simulated data and searching over the four parameters {h 0 , cos ι, ψ, Φ C 22 }.   Figure 17 and the second line of Table 2. A jupyter notebook with this test can be found here.

sites) as shown in
The above tests using simulated noise, and with searches covering the four standard unknown gravitational wave parameters, show that the code appears to be working as expected, and is in agreement with our previously used code. However, more generally we can perform parameter estimation over a larger number of signal parameters, and can again check for consistency with the previous code. We produce one day of Gaussian noise for a single detector and search over the four gravitational wave parameter with the same priors as before, but then also use a multi-variate Gaussian prior ( §2.3.2) over the phase parameters: rotational frequency and first frequency derivative, and binary parameters of binary period, time of periastron, angle of periastron and its first derivative, projected semi-major axis, and eccentricity ({f,ḟ , P b , T 0 , ω 0 ,ω 0 , a sin i, e}). This gives a total of 12 parameters in the search. The multi-variate prior is defined such that all parameters are uncorrelated except, in this case, the parameter pairs [T 0 , ω 0 ] and [P b ,ω 0 ] for which we use a very high correlation of 0.9999 (we do not Fig. 17.-Marginalised one and two-dimensional posterior probability distributions for the parameters {h 0 , cos ι, ψ, Φ C 22 } produced from the outputs of the nested sampling code lppen, and the previous code lppe, when running on simulated Gaussian noise assumed from the LIGO H1 and L1 detectors. The lppe code was run with both an MCMC sampler mode and a grid-based evaluation mode (shown as the black dashed line). set them to be fully correlated due to numerical issues inverting such matrices).
The one-and-two dimensional posteriors (expressed as offsets from the assumed heterodyne parameter values) output for this case when run with both lppe and lppen can be seen in Figure 18, which also overlays the marginalised priors on top. It can be seen qualitatively that the codes are consistent 15 and for the phase parameters the priors are recovered, except forḟ , which has some structure due to the specifics of the noise realisation. The Kolmogorov-Smirnov test p-values testing the null hypothesis that the posterior samples from each code are drawn from the same distribution are shown in Table 3, and suggest that the distributions are consistent.

Simulated signals
The next test is checking how the codes compare when there is a signal present in the data. When gener -TABLE 3 Kolmogorov-Smirnov test p-values testing the null hypothesis that the samples output by lppen and lppe, when running on simulated noise data and searching over the twelve parameters {h 0 , cos ι, ψ, Φ C 22 , f,ḟ , Pb, T 0 , ω 0 ,ω 0 , a sin i, e}, are drawn from the same distributions. ating the simulated signals 16 we initially assume that the signal's phase evolution is perfectly known and has been removed via the heterodyne described in §2.1, and therefore ∆φ 2 = 0 from Equation 10. As in §4.1.1 we search over the four-dimensional parameter space of θ = {h 0 , cos ι, ψ, Φ C 22 }. We create a simulated signal in the H1 detector spanning ten days of data with parameters (h 0 = 6.2×10 −23 , Φ C 22 = 2.4 rad, cos ι = 0.3 and ψ = 0.1 rad) that produce a signal-to-noise ratio of 8 when injected into Gaussian noise. The posterior probability distributions estimated for the parameters (expressed as offsets from the assumed heterodyne parameter values) are shown in Figure 19, which show consistency between the old and new codes and with the injected signal parameters. The quantitative consistency between the codes can be seen in Table 2, where is should be noted that the evidence ratio between the codes is larger than statistical uncertainty would expect, although the fractional difference between the values ( 1%) is still small enough to not greatly effect conclusions drawn from either values if used in model selection. Fig. 19.-Marginalised one and two-dimensional posterior probability distributions for the parameters {h 0 , cos ι, ψ, Φ C 22 } produced from the outputs of the nested sampling code lppen, and the previous code lppe, when running on simulated data containing Gaussian noise and a signal (of SNR 8) assumed from the LIGO H1 detector. The lppe code was run with both an MCMC sampler mode and a grid-based evaluation mode (shown as the black dashed line). The vertical dashed lines show the simulated signal parameter values.
Similarly, using a signal with a coherent multi-detector signal-to-noise ratio of 8 (for the same parameters as above except with h 0 = 3.5×10 −23 ) when simulated and injected into Gaussian noise for two detectors (H1 and L1) we see the posteriors shown in Figure 20, and con-sistency checks in Table 2. The notebooks for the above two tests can be found here and here. Fig. 20.-Marginalised one and two-dimensional posterior probability distributions for the parameters {h 0 , cos ι, ψ, Φ C 22 } produced from the outputs of the nested sampling code lppen, and the previous code lppe, when running on simulated data containing Gaussian noise and a signal (of coherent multi-detector SNR 8) assumed from the LIGO H1 and L1 detectors. The lppe code was run with both an MCMC sampler mode and a grid-based evaluation mode (shown as the black dashed line). The vertical dashed lines show the simulated signal parameter values.
A set of simulated signals with sky location and orbital parameters set to those of the Low-mass X-ray binary Sco X-1 were injected into Gaussian noise for the study in Messenger et al. (2015). As a test of our code when recovering a signal from a source in a binary system we have recovered one of these signals using approximately 1.4 days of simulated data. We have purposely performed the heterodyne data processing stage with values of several of the phase parameters (f 0 , α, T 0 , a sin i and P b ) set to not match the known simulated signal values. Therefore, to recover the signal, in addition to the four gravitational wave parameters we have had to allow the code to search over these extra parameters. For these additional parameters Gaussian priors were used, with the prior means set to be the parameter values used for heterodyning (not those for the actual signal), with standard deviations wide enough to encompass the true signal values. For h 0 a Fermi-Dirac prior was used, whilst flat prior were used for the other angle parameters, covering their minimal allowed ranges. The recovered signal posterior distributions can be seen in Figure 21, which show that the true parameters are correctly recovered (the plot shows recovered phase parameters as offsets from the heterodyne parameter values, i.e. the centres of the Gaussian priors). It can be seen that for f 0 and P b the posteriors contain the correct value, but are peaked well away from the value used for heterodyning (which would be at zero in the plots), showing the code's ability to explore the prior space in these parameters. The onedimensional marginal distribution for φ 0 spans the full prior range, however, this is due to it being highly correlated with f 0 , and, although not visible on the plot, this very strong correlation is present in the two-dimensional posterior for these parameters. Finally, for α and T 0 , it can be seen that the posteriors just fill the prior ranges as the data provides no additional information about these parameters. For completeness, we find that the signal was recovered with SNRs of ∼ 20 in both H1 and L1 individually, with a coherent SNR of ∼ 28, and odds values for the multi-detector analysis of log 10 O S/N = 152 and log 10 O S/I = 7.8.

Results evaluation
The above results have shown good qualitative, and quantitative, agreement between the old and new codes for individual cases, but it is also useful to see how they compare for a large number of noise realisations when assessing two important quantities: the 95% upper limit on gravitational wave amplitude h 0 , and the signal model evidence. It is also useful to see how these compare as a function of the number of live points used by the code. To this end we have chosen numbers of live points from 256 to 16384, increasing in a power of two for each step, and for each number generated 500 realisations of complex Gaussian noise with 1440 points over one day. Choosing a random source sky location for each realisation, and using uniform priors for the parameters h 0 , cos ι, Φ C 22 and ψ, we have performed parameter estimation using both lppen and lppe in its grid-based mode. For lppe the integrals required for calculating the signal evidence and 95% upper limit on h 0 are performed using the trapezium rule, whilst for the output of lppen the upper limit has been calculated using a greedy-binning approach bounded from zero. This assessment is very similar to what we did in §3.1.1, but has not just used a fixed Gaussian likelihood function in each case. We will make the assumption that the grid-based results are a representation of the true values for the evidence and upper limit, however, it should be noted that there will in fact be errors (and potentially biases) from the trapezium rule approximation.
In Figure 22 we see a comparison between the signal evidences evaluated using nested sampling and the trapezium rule. The left-hand axis shows the logarithm of the ratio of the two values, whilst the right hand axis shows the actual percentage difference between the evidences. These are plotted as a function of the number of live points used in the nested sampling evaluation by lppen. We find that the distribution of evidences very closely follows the theoretically expected distribution calculated as σ Z = H/N live (see §3.1.1), where given our prior range set up we found H ≈ 2.45. It can also be seen that there is a slight bias towards the nested sampling evidence being underestimated by a few percent, as was also observed for the simple Gaussian likelihood function case in Figure 11.
In Figure 23 we see a comparison between the 95% upper limits on h 0 as a function of the number of live points. The core distribution of upper limits produced via nested sampling are centred around the true (or gridbased) values, with the distribution decreasing in width with increasing numbers of live points. The standard deviation on the distribution roughly follows 2 6.8 N −1/2 live .
The figure shows some extreme outliers in the distributions, but it is found that these outliers occur in the cases when the number of posterior samples drawn from the nested samples is extremely low (the lower tails of the distributions in Figure 14).

Injections into real data
During all bar the first initial LIGO science run, simulated signals from a variety of sources have been physically "injected" into the detector datastreams. These are known as "hardware injections" (see Biwer et al. 2016, for a general discussion of hardware injections, in particular relating to their use, and extraction, in advanced LIGO's first observing run). These have included a range of continuous wave signals of the form given by, e.g., Equation 2 of Abbott et al. (2017), which in past have been searched for using a variety of analysis pipelines, including the one described in this document using lppe (see, e.g., Appendix B of Abbott et al. 2007).
Here we have used lppen to perform parameter estimation on the parameters θ = {h 0 , cos ι, ψ, φ 0 } for two hardware injections added to intial LIGO data from the sixth science run (S6). 1718 These are two out of the thirteen total continuous wave injections and were called "Pul-sar03" and "Pulsar05" respectively. The data have been heterodyned with the known phase evolution of the signals. The extracted posterior distributions can be seen in Figures 24 and 25 and are in very good agreement with the injection parameters. 19 For "Pulsar03" the signal is recovered with SNRs of 164 and 96 in H1 and L1 respectively, and a coherent SNR of 190, and with log 10 O S/N = 6729 and log 10 O S/I = 9.1. For "Pul-sar05" the signal is recovered with SNRs of 9.9 and 6.4 in H1 and L1 respectively, and a coherent SNR of 11.8, and with log 10 O S/N = 22.5 and log 10 O S/I = 4.6. The odds would imply a very high significance in distinguishing the coherent signal model from the noise models.
Several thousand simulated signals have been added to the LIGO sixth science run dataset (see Section IV of Walsh et al. 2016). For one particular signal (chosen for its moderate SNR over the run) we have used lppen to perform parameter estimation on the parameters θ = {h 0 , cos ι, ψ, φ 0 , f,ḟ , α, δ}. In this case, as for the binary signal discussed in §4.1.2, we have purposely heterodyned the signal using values of f andḟ offset from their known rotational values by ∆f = 5×10 −8 Hz and ∆ḟ = −4×10 −15 Hz s −1 . 20 We have used the whole of the S6 dataset, spanning almost 1.3 yr, with roughly 50% duty factors for both H1 and L1. The code was run using the ROQ mode for each individual detector and a joint detector analysis. Gaussian priors were used on 17 φ 0 is being used here as opposed to Φ C 22 as we are not having to compare to the previous code.
19 For these hardware injections exact agreement between the expected injection parameters and those recovered, or between detectors, is not entirely expected. This is due to the injections needing to be performed using actuation functions (converting the expected gravitational wave strain into a force required to be exerted in the interferometer test mass) that may not exactly match the actuations functions later required for calibrating the detectors. 20 The frequency offset is approximately two Fourier frequency bins away from the actual signal frequency. -Marginalised posterior plots for the recovered parameters of a simulated binary signal (with binary parameters similar to those of the low-mass X-ray binary system Sco X-1) using simulated LIGO H1 (red) and L1 (green) data, including a joint detector analysis (grey). The location of the simulated signal parameters are marked with a vertical black dashed line on the one-dimensional posteriors, or black cross on the two-dimensional posteriors. The priors used for each parameter are shown as the dark magenta dashed lines. the frequency and sky location parameters. The recovered marginalised posterior distributions can be seen in Figure 26, where all parameters are observed to be correctly recovered at their true values, even f andḟ which we offset (it is hard to see in the plots, but these values are constrained well within their allowed prior ranges). This signal was recovered with SNR of roughly 29 and 21 in H1 and L1 individually, and with a coherent SNR of 35. The odds values are log 10 O S/N = 251 and log 10 O S/I = 9.5 showing that the signal is very significantly distinguished from the noise model.

Monte-Carlo studies
In this section we will assess the outputs of lppen using a range of simulated signals. We have created two sets of 2000 signal parameters with which to inject simulated signals into Gaussian noise (with the same standard deviation of 1×10 −22 ) for the H1 and L1 detectors. Both sets have been generated with φ 0 , ψ and cos ι drawn from uniform distributions over their minimal allowed ranges (Table 1 of , but for the first set h 0 is calculated such that the coherent SNR (Equation 28) is drawn uniformly between 0 and 20; the second set draws the h 0 uniformly between 0 and 3.25×10 −22 , such that the maximum coherent SNR (given the known noise level) is ∼ 55. This latter set is used to evaluate the posterior distributions as described below in §4.3.1, whilst both sets are used to evaluate odds distributions. The sources' sky locations are drawn uniformly from over the sky-sphere. For these analyses each simulated data set is one solar day long, consisting of 1440 complex data points (generated as if heterodyned at precisely the known source on 500 simulated Gaussian noise realisations for each N live value. The ratio compares the evidence produced by the nested sampling algorithm in lppen to that produced from a grid over the parameter space with lppe (which we assume is a representation of the true value). A comparison of the estimated error on the evidence value using the information gain, H (which is ∼ 2.4 here), with the measured standard deviation of the distributions is also plotted and agree rather well. phase evolution) sampled once per minute. This data has been used to estimate the four parameters h 0 , cos ι, φ 0 and ψ, with each being given a flat prior: the angles and cosine of orientation angle over their minimal ranges, and h 0 between zero and 1 × 10 −20 (well above the extent of the bulk of the likelihood). When running lppen on these, two parallel runs with 1024 live points where used, along with the default sampler proposals described in §3.1, and with the samples from both runs being combined to produce posterior samples and evidence estimates.  -Posterior probability distributions for the recovered parameters of a hardware injection (named "Pulsar05") into LIGO sixth science run data. The black dashed lines and black crosses represent the expected signal parameters. Note that the probability contours for the joint detector analysis, in particular for the φ 0 versus ψ plot, are broader than the true distributions due to the smoothing kernel not handling the parameter degeneracy well.
Another set of 2000 simulations has been created that again create simulated signals and add them to Gaussian noise for the H1 and L1 detectors. In this case the total SNR shared between the detectors has be drawn from a uniform distribution between 0 and 30, however, the signals have been purposely created to be incoher- ent between detectors to simulate instrumental noise features. The parameters φ 0 , ψ and cos ι are drawn from their minimal allowed ranges, but are different for the signals input into the two detectors, whilst the sky positions are also independent for both detectors. The signal rotational frequency and frequency derivative are also offset between detectors, such that the offsets are drawn from Gaussian distributions of (1σ) width 1/2048 Hz and 1×10 −10 Hz s −1 respectively. The length of the data is the same as in the previous case, but it is assumed that the data for both detectors was heterodyned using the known phase evolution for the signal in only the H1 detector. The same prior and run settings as above were used. These simulations have been used for assessing the odds for a coherent versus incoherent signal (see §2.6).
Finally, two sets of simulations of 500 signals have been created to assess the code when estimating additional phase parameters: in this case the rotational frequency f 0 , frequency derivativeḟ 0 , and right ascension α. Gaussian priors on these three parameters have been used, with parameters given in Table 4. These Gaussian priors are used when generating signals and when estimating them using lppen, whilst the simulated data is created such that it was heterodyned using the mean values. As above, one set of these simulations has amplitudes calculated such that the coherent SNRs are drawn uniformly between 0 and 20, whilst the other has amplitude drawn uniformly between 0 and 3.25×10 −22 . The latter of these is used to evaluate the posterior distributions, whilst both are used for odds distributions assessment. Again, the angles are drawn uniformly from their minimal ranges and sky positions are drawn uniformly over the sky-sphere. Parameter mean standard deviation f 0 (Hz) 100 5×10 −5 f 0 −1×10 −9 2×10 −10 α (rads) π 0.0007272 (10 s )

Evaluating the posterior distributions
In Sidery et al. (2014) a method was developed for evaluating whether sky location credible intervals for compact binary coalescence gravitational wave signals, produced from posterior distribution output by LALInference codes (Veitch et al. 2015), behaved selfconsistently. The approach uses injected signals, with parameters drawn from the prior distributions used when reconstructing them 21 , and sees if the credible intervals effectively behave as frequentist confidence intervals, i.e. do 50% of the known parameter values fall within some pre-specified definition (like the minimal credible region, or credible region bounded by the lower extent of the prior) of the 50% credible region. This can test whether the underlying LALInference parameter sampling is working as expected, or if there are biases to wider, or narrower, credible regions. The method, and the "P-P plots" it produces through evaluation over a range of credible intervals, has now become commonly used to assess the LALInference parameter estimation codes (Veitch et al. 2015) as a way of highlighting any problems.
We have used this method here to evaluate our parameter posteriors produced from the simulations described above for which the amplitudes were drawn from a uniform distribution. For each simulation we have calculated the minimal credible region from the onedimensional marginalised posterior samples that the injected value falls within (using a greedy-binning approach). The cumulative histogram of these values, for each parameter when using the simulations that just estimate the four gravitational wave parameters, can be seen in Figure 27. Also, shown for each parameter is a Kolmogorov-Smirnov test p-value comparing our observed distribution to a uniform distribution, and a "cloud" of cumulative distributions that are random realisations of the expected cumulative distribution. We see that, overall, our posteriors give credible intervals that behave as expected, without any large deviations.
Similarly, we have done the same thing for the posteriors produced from the simulations that included searches over f 0 ,ḟ 0 , and α. In this case we know that the parameters are f 0 and φ 0 are completely correlated (i.e. there is a degeneracy), so we have held the value of φ 0 fixed at zero. 22 These "P-P plots" can be seen in Figure 28, 21 In this case the amplitudes of the simulations are not drawn from the full prior range, but are drawn from a flat distribution within that range, so are valid for this comparison. 22 In cases where there is a complete degeneracy between parameters the samplers used in the code can have problems, leading to observable effects in the "P-P plots". In certain cases, where the degeneracies are known, there are ways to deal with this: one parameter of a degenerate pair can be held fixed, a custom pro- where generally we see posterior credible intervals behaving as expected.

Evaluating the odds
It is useful to see how the odds values, described in §2.6, look when calculated for our simulations. As mentioned above, along with coherent signals between two detectors, we have purposely created incoherent signals to see how the odds look. Firstly, just considering the coherent simulations, we can see in Figure 29 how the odds O S/N , O S/I and O S/I simple (Equations 30, 32 and 31, respectively) vary as a function of coherent SNR. We see that O S/N rapidly increases, O S/I simple has a shallow growth, whilst O S/I transitions between following O S/N at very low SNR to matching O S/I simple at higher SNR. This behaviour is consistent with our expectations described in Appendix D with ln O S/N ∝ ρ 2 coh and O S/I simple ∝ ρ coh . We note that the locations of these distributions on the yaxis can be strongly effected by the prior volume, changes to which will move the distributions up or down, and can affect where they intersect. However, the general shape of the distributions should be very similar.
In Figure 30 there are two panels showing O S/I as a function of O S/N . In the upper panel (showing the coherent simulations) the colours of points show the signal's posal distribution can be hard coded for degenerate parameters that deals with the correlations, or, reparameterisations can be found that do not have the degeneracies. In the future the latter of these may be implemented in our code to deal with strong degeneracies between the initial phase and derivatives of it (frequency, and further frequency derivatives). Such a parameterisation could involve using, e.g. the start and end phase, rather than φ 0 and f 0 , with further derivatives being incorporated through additional phase parameters at fixed times throughout the model. So, as another example, rather than internally sampling in φ 0 , f 0 , andḟ 0 , one could use φstart, φ mid , and φ end , and convert between the two using linear algebra. coherent SNR. The plateauing of O S/I is just that seen in Figure 29. The bottom panel shows a zoomed out version, with the incoherent signals as dark circles (these are also seen in the lower left of the upper panel). Here we see that for incoherent signals the value of O S/I between the coherent and incoherent signals diverges quite quickly, with log 10 O S/I falling off roughly linearly with log 10 O S/N . This shows that, hopefully, for reasonable strength signals the value of O S/I will be useful for distinguishing a coherent from an incoherent (and therefore not astrophysical) signal. 23 It is useful to see how these distributions change when the parameter space searched over increases in complexity. So, in Figure 31 we see similar plots to those above, but for the simulations that included searches over f 0 , f , and α (but fixed φ 0 ). In the left panel we see O S/N 23 The caveats to this are that the way we created our incoherent signals may not well match true signals caused by instrumental (non-astrophysical) effects in gravitational wave detectors. However, we hope our simulations are extreme, in a conservative way, in that one detector is actually containing an unadulterated signal, and is therefore potentially producing a far larger O S/N value than might be expected whereas an instrumental feature would generally not match our signal template, even in a single detector. and O S/I as a function of ρ coh for these simulations, with the earlier simulations plotted underneath. It is obvious that, for these injections the odds follow the same form as before, however, at low SNR there is more scatter in O S/I than seen before. This can be interpreted as finding spikes in f 0 -ḟ -space, that are not related to the signal, in individual detectors, which are therefore not coherent. For higher SNRs ( 10) the O S/I value for this more complex search does tend to give higher values than for the simpler search (shown in Figure 31 as the lighter circles). A hand-wavy explanation for this could be that a truly coherent signal found (at large enough SNR) within this larger space is more convincing (more parameters have to match up) than in the simpler case, and this outweighs any prior Occam factor. The right panel of Figure 31 shows similar information.

Further parameter estimation
In Figures 32 and 33 we show the posterior distributions for two of the simulations performed for the searches over h 0 , cos ι, φ 0 , ψ, f 0 ,ḟ , and α above. 24 The simulation in Figure 32 had injected SNRs of 4.4 and 6.5 for H1 and L1 respectively, with a coherent SNR of 7.9, and the recovered values of log 10 O S/N and log 10 O S/I of 1.55 and 0.95 respectively. It can be seen that the search in H1 alone finds additional structure in the frequency space, whilst in L1 the signal dominates, and the joint 24 Earlier, during the discussing of the "P-P plots" in §4.3.1, it was mentioning that φ 0 was held fixed due to the degeneracies between f 0 and φ 0 . In this case, as we also haveḟ , the correlations become more complex, and whilst it would be expected that holding φ 0 fixed would give valid posteriors, it is in such cases that a re-parameterisation may be more appropriate. analysis also recovers the signal. The simulation in Figure 33 is a little stronger with injected SNRs of 7.3 and 8.5 for H1 and L1 respectively, and a coherent SNR of 11.2, and recovered values of log 10 O S/N and log 10 O S/I of 23.8 and 6.3 respectively. The signal is by far the strongest feature in the data, and no other structure is present in the posteriors. The first signal would probably only rate as marginal if determining whether it was a true signal or not, whereas the second signal would be a clear detection.

CONCLUSIONS
In this document we have described some of the workings of a code, lalapps pulsar parameter estimation nested, that is used for performing parameter estimation and Bayesian model selection in searches for gravitational wave signals from known pulsars. We have described the available signal models, likelihood functions and prior functions. We have compared the code to one that had been previously used for the same purpose, lalapps pulsar parameter estimation, and found that the results agree very well for a range of situations. The comparisons have not been completely exhaustive, and there are situations that have not been tested, however, it is worth noting that lppe is not necessarily expected to provide complete ground truth. The code has also been validated on additional software and hardware simulated signals and appears to return posterior distributions that behave as expected, and evidence values that show the expected behaviour.
As a test of the internal nested sampling algorithm, and proposal distributions, used by the code a large number of simulations were performed using a simple Gaussian likelihood function. These interestingly highlighted that if using, in particular, the ensemble walk proposal distribution the evidence values output became systematically more biased away from the true value as a function of information gain between the prior and posterior distributions. The test showed that this bias was alleviated by including a proposal distribution that drew values from the full prior a certain fraction of the time. The reason for this systematic bias is currently unknown, but in general, given the magnitude of the bias, it should have a relatively minor effect. It was also shown that with the ensemble walk proposal source amplitude upper limits could be systematically biased by about 2%, however, again this effect appears to be alleviated with the flat proposal distribution.
Some example use cases of lppen are documented in Appendix A.

A. CODE USAGE
Here we document the various command line options for the code and its use in a variety of situations. The majority of the command line options for lppen, as shown with the --help command, are given (in a slightly edited form) in Table 5.  Display these commands.

--verbose
Display progress information from the code, including details of the nested sampling progress.

--detectors
Input all interferometers whose data is to be analysed, e.g. H1,L1 (delimited by commas) for the LIGO Hanford and LIGO Livingston detectors (see the discussion in §A.1.4 for the allowed detector names). If generating fake data, described below, these should not be set, otherwise this is required.

--harmonics
The signal model frequency harmonics that you want to use (delimited by commas). Currently this can be either a single value (e.g., the default is to assume a 'triaxial' star emitting from the l = m = 2 mode and use a value of 2), or 1,2 for a model with emission at both the rotation frequency and twice the rotation frequency.

--input-files
The files containing the pre-processed (heterodyned, or spectrally interpolated) data for each detector and model harmonic specified by --detectors and --harmonics delimited by commas (these files must be listed in the same order as passed to --detectors and --harmonics, with the files for different harmonics for individual detectors listed together). If this is not set you can generate fake data (see --fake-data below), but it is otherwise this is required.

--outfile
The name for the output data file. This must be a HDF5 file with the extension '.hdf' or '.h5'. Required.

--output-chunks
Output lists of the stationary chunks into which the data has been split (see §2.4). Optional.

--chunk-min
The minimum stationary length of data to be used in the likelihood. This defaults to 5, e.g., 5 mins for data sampled at 1/60 Hz. Optional.

--chunk-max
The maximum stationary length of data to be used in the likelihood. This defaults to 0, which is the value for no maximum length to be applied, except if the --oldChunks flag is set, in which case it defaults to 30. If running with the --roq flag it can be worth setting a value here (e.g., 1440) to avoid large amounts of memory being required in training set generation. Optional.

--time-bins
Number of time bins over a sidereal day in the antenna response lookup table. This defaults to 2880 bins. Optional.

--prior-file
The file containing the prior ranges for parameters to search over (see §A.1.2). Required.

--ephem-earth
The Earth ephemeris file. If not supplied this will attempt to be found in the path.

--ephem-sun
The Sun ephemeris file. If not supplied this will attempt to be found in the path.

--ephem-timecorr
The Einstein delay time correction ephemeris file. If not supplied this will attempt to be found in the path.

--biaxial
Set this if using the waveform model parameters with two harmonics, and specifically for a biaxial star (see, e.g., ).

--gaussian-like
Set this if a Gaussian likelihood is to be used. If the input file(s) contains a column specifying the noise standard deviation of the data then that will be used in the Gaussian likelihood function, otherwise the noise variance will be calculated from the data.

--nonGR
Set this to allow non-GR polarisation modes and/or a variable speed of gravitational waves.

--randomise
Set this, with an integer seed, to randomise the data (through permutations of the time stamps) for use in Monte-Carlo studies. Note that this will not work if using the code to create injections.

--truncate-time
Only analyse data up to the given GPS time. Optional. --truncate-samples Only analyse data samples up to the number given here. Optional.
--truncate-fraction Only analyse the given fraction of data samples (this must be a value between 0 and 1). Optional.

Nested sampling parameters: --Nlive
Set the integer number of live points for nested sampling. Required.

--Nmcmc
Set the length of the MCMC used to find new live points (if not specified an adaptive number of points is used).

--Nmcmcinitial
Set the number of MCMC points to use in the initial resampling of the prior. This defaults to 5000, but for our analyses this can be set to zero.

--tolerance
The tolerance used as the stopping criterion for the nested sampling integrator. This defaults to 0.1.

--randomseed
A seed for the random number generator. By default, if not supplied, this is set by the system clock.
MCMC proposal parameters: (see §3.1) --diffev The (integer) relative weight of using differential evolution of the live points as the proposal. The default is 0, e.g., 0%.

--freqBinJump
The (integer) relative weight of using jumps to adjacent frequency bins as a proposal. The default is 0, and this is not required unless searching over frequency.

--ensembleStretch
The (integer) relative weight of the ensemble stretch proposal. The default is 0. Note that this proposal greatly increases the parameter autocorrelation lengths, so in general should be avoided.
Reduced order quadrature (ROQ) parameters: --roq Set this flag to use reduced order quadrature to compute the likelihood.

--ntraining
The (integer) number of training models used to generate an orthonormal basis of waveform models.

--roq-tolerance
The tolerance used during the basis generation. The default is 1e-11. --enrich-max The (integer) number of times to try and "enrich" (see, e.g., Smith et al. 2016) the basis set using new training data. The enrichment process stops before this value is reached if three consecutive enrichment steps produce no new bases. The default is 100.

--roq-uniform
Set this flag to cause training model parameters, for parameters with Gaussian prior distributions, to be drawn from a uniform distribution spanning µ ± 5σ. Otherwise, by default, parameters are drawn from their given prior distributions.

--output-weights
Give an output file name to which to output (a binary version of) the weights will be output. If this is set the programme will exit after outputting the weights. These files could be read in later instead of being regenerated. This allows the potential for the ROQ to be generated on a machine with a large amount of RAM, whilst the full parameter estimation can run on a machine with less RAM.

--input-weights
A binary file, created using the --output-weights command, containing all the weights in a defined format. If this is present then the ROQ will not be recalculated.

Signal injection parameters:
--inject-file A pulsar parameter (.par) file containing the parameters of a signal to be injected. If this is given a signal will be injected.

--inject-nonGR
A parameter (.par) file containing the values of specific non-GR signal parameters to be injected (this is an alternative to --inject-file, so both should not be passed at the same time).

--inject-output
A filename to which the injected signal will be output if specified.

--inject-only
If this is set, and an --inject-output file is set, then do not perform nested sampling on the created injection, but just exit the code after creation of, and writing out of, the injection.

--fake-data
A list of interferometers for which fake data will be generated, e.g., H1,L1 (delimited by commas). Unless the --fake-psd flag is set the power spectral density for the data will be generated from the noise models in described in §A.1.4. The noise will be white across the data bandwidth.

--fake-psd
If wanting to generate fake data with specific power spectral densities for each detector given by --fake-data then they should be specified here delimited by commas (e.g., for --fake-data H1,L1 then you could use --fake-psd 1e-48,1.5e-48) where values are single-sided PSDs in Hz −1 .

--fake-starts
The start times (in GPS seconds) of the fake data for each detector separated by commas (e.g., 910000000,910021000). If not specified these will all default to 900000000.

--fake-lengths
The length of each fake data set (in seconds) for each detector separated by commas. If not specified these will all default to 86400 (i.e., 1 solar day).

--fake-dt
The data sample rate (in seconds) for the fake data for each detector. If not specified these will default to 60 s.

--scale-snr
This gives a (multi-detector) SNR value to which you want to scale the injection. By default this is 1.
Legacy code flags: --oldChunks Set this if using fixed chunk sizes for dividing the data as in the old code, rather than the calculating chunks using the change point method.

--source-model
Set this if using both 1 and 2 multiples of the frequency and requiring the use of the original source model parameters from .
Benchmarking: --time-it Set this if wanting to time the various parts of the code. A output file with the "outfile" filename appended with " timings" will contain the timings --sampleprior Set this to output this (integer) number of samples generated from the prior. The nested sampling will not be performed.
A.1. Examples Here we show a few real examples of using the code for a variety of situations. In all cases there are a minimum of three files that must be specified to run the code: 1. One, or more, data files containing a complex time series. These will generally be (although do not have to be) the output of the heterodyne method described in (Dupuis and Woan 2005) (and implemented with lalapps heterodyne pulsar), or the Spectral Interpolation method described in (Davies et al. 2017). Files must be in ascii text format with three, or four, whitespace-separated columns. In either case the first three columns specify the data timestamp in GPS seconds, the real part of the data, and the imaginary part of the data, e.g., 900000000 1.2564374e-22 -4.5764347e-22. In the case where a fourth column is present this should give the noise standard deviation, and is for use if running with the Gaussian likelihood function. The code will automatically determine if the input contains three or four columns. The data does not have to be uniformly sampled and can contain gaps. However, if searching for a signal for which there could still be some frequency modulation it is useful for the data to have a large enough bandwidth to contain the modulation. The file can contain comments if the a line starts with a % or #. Gzipped versions of these files can also be input.

2.
A pulsar parameter file. This should be an ascii text tempo(2)-style .par file containing information on the source's parameters, e.g. sky location and frequency, as generated by tempo 25 or tempo2  from pulsar pulse time-of-arrival observations. An example .par file is given below.
3. A prior distribution file. This should be an ascii text file containing information on the prior distributions required for any variable parameters for which the code needs to sample posterior probability distributions. An example prior file is given below.
In the examples below we will assume that the LALSuite software package (LIGO Scientific Collaboration 2017) has been installed and that the executable binaries for lppen and lalapps nest2pos are in your path (i.e. the directory path containing those executables is set in the list of values in your PATH environment variable). The code can automatically find solar system ephemeris and timing files provided they are in a standard location, although these can also be specified explicitly on the command line. The .par file input should be the same file that was used when performing the heterodyning (or Spectral Interpolation) of the data used to produce the input data time series. A tempo(2)-style .par file could have the following form, e.g., for a pulsar in a binary system: The majority of the parameters that can be in a .par file are listed in Table A1 of . Note that our code does not compute effects of dispersion measure (i.e. it assumes observations at infinite frequency), or atmospheric effects, so it is preferable to use a .par file for which these parameters have been not used, or held fixed.
For pulsars in binary systems, the code can currently use the following binary system models (see, e.g., Taylor and Weisberg 1989;Pitkin and Woan 2007, for discussion of some of the models): BT, BTX, BT1P, BT2P, ELL1, DD, DDS, MSS, and T2 (although not all aspects of the T2 model are implemented). The code can use the TDB or TCB time systems, with the TCB system being the default if not otherwise specified (either by supplying the require timing files, or by setting "UNITS TDB" in the .par file). The code can accept .par files calculated using the JPL DE200, DE405, DE414, or DE421 solar system ephemerides, noting that it only takes into account the motion of the Earth-Moon system and the Sun, but no other solar system bodies.
The .par file can also contain a selection of gravitational wave parameters not present in the case of a standard electromagnetically observed pulsar. A full list of the allowed parameter names for a .par, and a prior file described below, is given in Table 6.
TABLE 6 A list of parameters allowed in a .par file or a prior file. The true-type font versions, in uppercase only, are what is allowed in the files. The specified units are for the .par file, whereas the prior file requires values in SI units.

Parameter Description
Gravitational wave parameters H0, h 0 gravitational wave amplitude for the l = m = 2 mode in the source model (see, e.g., Jaranowski et al. 1998) IOTA, ι pulsar inclination angle to the line-of-sight (rads) (0 • is gives a source whose rotation axis is pointed directly along the line-of-sight) COSIOTA, cos ι cosine of ι PSI, ψ polarisation angle of the gravitational waves (rads) (see, e.g., Jaranowski et al. 1998 Kopeikin term (arcsec −1 ) (see, e.g., Section 2.7.1 of ) A1 X projected semi-major axes for multiple orbital components (X = 2, 3 for BT1P and BT2P models) (lt-s) PB X periods for multiple orbital components (X = 2, 3 for BT1P and BT2P models) (days) OM X longitudes of periastron for multiple orbital components (X = 2, 3 for BT1P and BT2P models) (degs) T0 X times of periastron for multiple orbital components (X = 2, 3 for BT1P and BT2P models) (MJD) ECC X eccentricities for multiple orbital components (X = 2, 3 for BT1P and BT2P models) FBX orbital frequencies for multiple components (in BTX model) (Hz) A0, A first aberration parameter (see, e.g., Section 2.7.3 of ) B0, B second aberration parameter (se, e.g., Section 2.7.3 of  Glitch parameters (see Yu et al. 2013) GLEP X glitch epochs (MJD) GLPH X glitch rotational phase offsets (rads) GLF0 X glitch frequency offsets (Hz) GLF1 X glitch first frequency derivative offsets (Hz s −1 ) GLF2 X glitch second frequency derivative offsets (Hz s −2 ) GLF0D X glitch decaying frequency components (Hz) GLTD X glitch time constant for decaying components (s) Non-GR parameters CGW, cgw speed of gravitational waves as a fraction of the speed of light HPLUS, H + gravitational wave amplitude for the tensor +-polarisation (see, e.g., Isi et al. 2017, for definitions of this and the subsequent parameters) HCROSS, H × gravitational wave amplitude for the tensor ×-polarisation PHI0TENSOR, Φt tensor mode phase parameter PSITENSOR, ψt tensor mode phase parameter HVECTORX, Hx gravitational wave amplitude for the vector x-polarisation HVECTORY, Hy gravitational wave amplitude for the vector y-polarisation PHI0VECTOR, Φv vector mode phase parameter PSIVECTOR, ψv vector mode phase parameter HSCALARB, Hx gravitational wave amplitude for the scalar b-polarisation HSCALARL, Hy gravitational wave amplitude for the scalar l-polarisation PHI0SCALAR, Φs scalar mode phase parameter PSISCALAR, ψs scalar mode phase parameter

A.1.2. Example prior file
The prior functions used by the code are described in §2.3. To implement them in the code requires a prior file specifying the prior function require for each parameter that will be varied. A an example prior file containing all the allowed prior types (for descriptive purposed only) is: The units required for parameters in a prior file are all SI units, whereas for the .par file they often follow a different convention (see Table 6).
Going through each line we have: H0: this parameter has been assigned a Fermi-Dirac prior ( §2.3.5) using the fermidirac tag, which requires two values: σ and r.
PHI0: this parameter has been assigned a uniform prior ( §2.3.1) using the uniform tag, which requires two values: a lower limit and an upper limit.
PSI: this parameter has been assigned a Gaussian prior ( §2.3.2) using the gaussian tag, which requires two values: a mean and a standard deviation.
A1: this parameter has been assigned a prior uniform in the logarithm of the parameter ( §2.3.4) using the loguniform tag, which required two values: a lower limit and an upper limit on the parameter (not the logarithm of the parameter).
F0:F1: this pair of parameters has been assigned a Gaussian mixture model prior ( §2.3.3). This prior can be used for any number of parameters if given as a colon separated list in the prior file. Following the gmm assignment, a number gives the number of modes for the model. This is then followed by a square-bracketed list of lists, with each mode giving a comma separated main list entry containing a sub-list of comma separated means for each parameter (no additional whitespace is allowed in the lists). The list of means is then followed by a similar list of parameter covariance matrices for each mode (again no additional whitespace is allowed). This is then followed by a final bracketed list giving the relative weights for each mode (the weights do not have to be normalised to unity, as the code will automatically normalise them). If required this can then be followed by pairs of bracketed, comma separated values, for each parameter giving their minimum and maximum allowed ranges. If these are given then they have to be given for all parameters, not just some subset. These are not shown in this example.
COSIOTA: this single parameter has been assigned a Gaussian mixture model prior ( §2.3.3). This shows how to input a Gaussian mixture mode for a single parameter, and also shows at the end minimum and maximum ranges for the parameter have also been input.
If parameters in the prior file are specified with the gaussian tag then a multivariate Gaussian prior may also be used for them if a correlation matrix is passed to lppen (using the --cor-file flag) containing correlation coefficients for the required parameters. Only the lower diagonal of the correlation matrix is required, e.g. if we wanted a multivariate prior on f 0 andḟ we could have a correlation coefficient file containing: where the whitespace can be tabs or spaces, whilst the prior file contains the means and standard deviations of the parameters as noted above. It is worth noting that, for fully (anti)correlated parameters correlation coefficients of one can lead to numerical issues during matrix inversion. So, in such cases the (anti)correlation could be reduced to, say, 99.99%. However, better solutions to such a degeneracy would be to fix one of the correlated values (i.e. do not include it in the prior file) and only searching over the other, or finding a non-degenerate re-parameterisation.
A.1.3. Example 1: single detector, single harmonic, input data One of the simplest cases for running the code is on input data from a single detector and at a single frequency harmonic (in this case at twice the rotation frequency). Let us assume a search that is just interested in a signal model described by the four parameters h 0 , φ 0 , cos ι and ψ, and set a prior file, called prior.txt, containing and assuming the input data (from the LIGO Hanford detector, H1, in this case) is in a file called data H1.txt, then one would run: $ l a l a p p s _ p u l s a r _ p a r a m e t e r _ e s t i m a t i o n _ n e s t e d --detectors H1 --input -files data_H1 . txt --par -file pulsar . par --prior -file prior . txt --outfile nest_H1 . hdf --Nmcmcinitial 0 --Nlive 1024 --tolerance 0.1 where we have set an output HDF5-format file called nest H1.hdf in which the nested samples will be output. This setup has used 1024 nested sampling live points and a tolerance of 0.1 for stopping the code. The --Nmcmcinitial 0 flag tells the code that the initial live points drawn do not have to be resampled. 27 In this case the solar system ephemeris files will be found automatically. If the EPHEM or UNITS parameters are not explicitly set in the .par file then these will default to the DE405 ephemeris and the TCB time system. The progress of the above code can be monitored by including the --verbose flag, and watching how the reported dZ value approaches the required tolerance. Along with the file of nested samples two additional files called nest H1 SNR and nest H1 Znoise will be output. The former contains the recovered optimal SNR for the maximum likelihood sample. The later contains the null likelihood, or noise, evidence for the data, described in §2.2.3. The noise evidence is also contained in the nest H1.hdf file, but this additional file is useful in cases when multiple harmonics are used as it will contain noise evidences for individual harmonic datastreams, as well as that for the multiple combined streams.
The nested sample file can be converted to a file containing posterior samples, called, say, post H1.hdf, using: $ l a la p p s _ n e s t 2 p o s -p post_H1 . hdf nest_H1 . hdf where the -p flag give the name of the posterior output file. If the above use of lppen had been run more than once on the same data and prior file, with output names given by e.g. nest H1 01.hdf...nest H1 N .hdf, then all these could be listed as input to lalapps nest2pos to combined all the samples and evidence estimates. As well as the posterior samples, this file also contains (as does the nest H1.hdf), the signal model evidence, the noise model evidence, and the information gain (amongst other things). The posterior samples, and noise and signal model (natural logarithm) evidences, can be extracted using a python function within lalapps, via $ python >>> from lalapps import pulsarpputils as pppu >>> postsamples , sigevidence , noiseevidence = pppu . p u l s a r _ n e s t _ t o _ p o s t e r i o r ( ' post_H1 . hdf ') The postsamples variable is an object from which samples from different parameters can be extracted, and plotted as a marginalised posterior distribution, e.g., >>> from matplotlib import pyplot as pl >>> # plot the posterior samples for the H0 parameter >>> pl . hist ( postsamples [ ' H0 ']. samples , bins =20 , normed = True , histtype = ' stepfilled ') >>> pl . xlabel ( ' $h_0$ ') >>> pl . ylabel ( ' Posterior Probability Density ') >>> pl . show () where here, for the h 0 parameter it is contained in postsamples as the fully uppercase 'H0' value.
Alternatives to extract the (natural logarithm) evidence values are using the h5py python module, or HDF5 command line utilities. For the former one could use, e.g., >>>  However, note that in this latter example the evidence values are truncated to integer values.
A.1.4. Example 2: single detector, single harmonic, simulated data As well as taking real data as an input the code can simulate Gaussian noise, and/or create and add a simulated signal to the data. To create a signal to inject into the data requires a .par file for the injected signal parameters. This can be the same, or different, to the file input using the --par-file flag, but to produce a non-zero signal it must contain at least a gravitational wave amplitude parameter (any parameters not included will be assumed to be zero). Any differences in the phase evolution parameters between the two files will be incorporated into the signal via Equation 10 (which we see is the same as in §A.1.3, but with gravitational wave amplitude parameters), one would run the following (assuming various other names are the same as in §A.1.3) $ l a l a p p s _ p u l s a r _ p a r a m e t e r _ e s t i m a t i o n _ n e s t e d --fake -data H1 --inject -file injection . par --scale -snr 10 --par -file pulsar . par --prior -file prior . txt --outfile nest_H1 . hdf --Nmcmcinitial 0 --Nlive 1024 --tolerance 0.1 Here, rather than inputting a real data file, the --fake-data flag has been used to specify fake data from the LIGO Hanford detector, where the design sensitivity (Lazzarini and Weiss 1996) 28 power spectral density noise curve would be used to set the Gaussian noise level for an initial LIGO detector at the appropriate gravitational wave frequency (twice the value of F0 given in the .par file for the standard search). 29 Setting the noise level from the design sensitivity for a given detector can be overridden by also using the --fake-psd flag to input the required single sided power spectral density value. The simulated data will by default start at a GPS time of 900000000 (13 Jul 2088, 15:59:46 UTC), be sampled at a rate of once per minute, and last for one solar day. However, each of these can be set with the appropriate flags: --fake-starts, --fake-lengths, and --fake-dt, respectively. In this example a value is also set for --scale-snr, which means that after the creating of a signal with parameters as given in the injection.par file, it will be re-scaled so that it has the given SNR -in this case a value of 10. This enables signals to be injected with a known SNR without needing to know the noise level a priori. If the --scale-snr flag is not set, or has a value of zero, then no re-scaling takes place. In this case the output nest H1 SNR file would also contain the signal's injected SNR and the recovered SNR, e.g.,

28
https://labcit.ligo.caltech.edu/~jzweizig/ distribution/LSC_Data/ 29 Other named detectors that can be provided are L1 for the LIGO Livingston detector, G1 for the GEO600 detector (from Table IV of Damour et al. 2001), V1 for the initial Virgo detector, and T1 for the TAMA300 detector (from Table IV  where, after the detector name the harmonic frequency scaling factor is given, and for the injected signal the two SNR values are the pre-scaled value and the actually injected value after rescaling via --scale-snr value.
A simulated signal can also be injected into real data in the same way.
A.1.5. Example 3: multiple detectors, single harmonic Running with multiple detectors requires a simple change to the example given in §A.1.3. Given two input files, e.g. data H1.txt and data L1.txt, for the two LIGO detectors H1 and L1, one would simply use, e.g.
$ l a l a p p s _ p u l s a r _ p a r a m e t e r _ e s t i m a t i o n _ n e s t e d --detectors H1 , L1 --input -files data_H1 . txt , data_L1 . txt --par -file pulsar . par --prior -file prior . txt --outfile nest_H1L1 . hdf --Nmcmcinitial 0 --Nlive 1024 --tolerance 0.1 where the order of listing the detectors with the --detectors flag should match the order of listing the data files with the --input-files flag. No whitespace is allowed between detector values or input file values, and they must be separated by a comma. In this case the output nest H1L1 SNR file will contain the individual detector recovered SNRs and the coherent SNR, whilst the nest H1L1 Znoise file will contain both the individual detector noise model evidences and the combined noise model evidence.
A.1.6. Example 4: multiple detectors, two harmonics It is simple to extend the analysis to include the two harmonics from the l = m = 2 mode (with a gravitational wave frequency at twice the rotation rate) and the l = 2, m = 1 harmonic (with a gravitational wave frequency at the rotation rate). To search at these two harmonics requires that the parameters searched over are either the signal or waveform parameters given in . Using the (simpler and non-degenerate) waveform parameters C 21 , C 22 , Φ C 21 and Φ C 22 , from Equations 7 and 8, could lead to the following (prior.txt) prior file C22 uniform 0 1e -20 C21 uniform 0 1e -20 PHI22 uniform 0 6.28318530718 PHI21 uniform 0 6.28318530718 PSI uniform 0 1.57079632679 COSIOTA uniform -1 1 Assuming that for each detector (H1 and L1) the data has been heterodyned at both harmonics, leading to files data H1 1f.txt, data H1 2f.txt, data L1 1f.txt, data L1 2f.txt, where 1f and 2f note the harmonic rotation frequency scaling, then one would use $ l a l a p p s _ p u l s a r _ p a r a m e t e r _ e s t i m a t i o n _ n e s t e d --detectors H1 , L1 --harmonics 1 ,2 --input -files data_H1_1f . txt , data_H1_2f , data_L1_1f . txt , data_L1_2f . txt --par -file pulsar . par --prior -file prior . txt --outfile nest_H1L1 . hdf --Nmcmcinitial 0 --Nlive 1024 --tolerance 0.1 Here the ordering of the input files is important, with the files for each harmonic for a given detector listed together, and in the same order as given by the --harmonics flag. The noise evidence and SNRs for each harmonic and each detector will be output to the files nest H1L1 Znoise and nest H1L1 SNR.
A.1.7. Example 5: single detectors, additional parameters Finally, we show an example when searching over additional phase evolution parameters. We show two options for this: the first is using the standard model and likelihood method, and the second performs the same search, but sped-up by using the ROQ likelihood. If one were searching over frequency and frequency derivative, f 0 andḟ , then an example .par file (pulsar.par) could be Note that the frequency and frequency derivative values in the prior file are for the rotation frequency, so if searching for the standard l = m = 2 harmonic then the actual frequencies and derivatives (and in this case their standard deviation) covered will be twice these values. This would also be the case for any other form of the prior. As, in this example, the F0 and F1 values have be assigned Gaussian priors, if they have a known correlation then a correlation coefficient file (pulsar.cor) could also be set, e.g. However, if this is not provided, then the parameters will be assumed to a priori be uncorrelated.
If running without ROQ, for a data from a single detector, H1 (data H1.txt), the following commands could be used $ l a l a p p s _ p u l s a r _ p a r a m e t e r _ e s t i m a t i o n _ n e s t e d --detectors H1 --input -files data_H1 . txt --par -file pulsar . par --prior -file prior . txt --cor -file pulsar . cor --outfile nest_H1 . hdf --Nmcmcinitial 0 --Nlive 1024 --tolerance 0.1 To run with ROQ enabled, in the simplest case, requires the user to set a number of training waveforms and tolerance to use the produce a reduced order model. 30 In this example we will use 2500 training waveforms, and a tolerance of 5×10 −12 , and use data that starts at the same epoch as given in the .par file. The following commands could be used $ l a l a p p s _ p u l s a r _ p a r a m e t e r _ e s t i m a t i o n _ n e s t e d --detectors H1 --input -files data_H1 . txt --par -file pulsar . par --prior -file prior . txt --cor -file pulsar . cor --outfile nest_H1 . hdf --Nmcmcinitial 0 --Nlive 1024 --tolerance 0.1 --roq --roq -tolerance 5e -12 --ntraining 2500 In this example, for a single run on a random realisation of Gaussian noise lasting one day and sampled at a rate of one per minute, we find that the ROQ version runs about 1.8 times faster. Within the given tolerance it reduces the prior parameter space to 32 orthogonal waveform templates, thus reducing the the model and likelihood evaluations from using 1440 to 32 points. Due to various overheads the speed-up we see between the ROQ and regular run is only a factor of ∼ 1.8, rather than the 1440/32 = 45 factor that might be hoped for. However, more impressive speed-up can be observed for longer data sets. It should be noted that this observed speed-up is not a general rule, but is specific to the data and set up we used, but similar speed-ups for similar set ups would not be unexpected. We also find that the likelihoods produced by the ROQ method agree very well with that produced using the full likelihood evaluation, with the maximum-likelihood-template likelihoods agreeing to within 6×10 −10 %.
A.1.8. Example 6: single detector, non-GR model The standard analysis considers two hypotheses: GR signal and Gaussian noise. On top of these, we can also include non-GR signal hypotheses consisting of all possible combinations of tensor, vector and scalar modes: GR+s, GR+v, GR+sv, s, v, sv (see Isi et al. 2017, for definitions of these). 31 For any set of data, we can compute Bayes factors for each of these hypotheses using lppen and then combine them to produce odds for any-signal vs. noise and non-GR vs. GR.
A.2. The pipeline The process of performing searches for known pulsars does not just use lppen. It actually requires a whole pipeline of code to: a) gather the gravitational wave detector data and operational segments ("science mode"); b) pre-processing the raw data via heterodyning (or spectral interpolation) to give the complex down-sampled time series; c) set up the required prior distribution files; d) pass the data and priors to lppen; e) convert the output nested samples into posterior samples; and, finally, f) display those results (including odds values) in a useful fashion, e.g. on a webpage. This process can be performed manually in a step by step manner, but there also exists a pipeline script called lalapps knope that gathers together all these operations in a way that can be run on a computer cluster running the HTCondor 32 job management system. # Condor accounting group user a c c o u n t i n g _ g r o u p _ u s e r = user . name # the data find executable ( e . g . / usr / bin / gw_data_find ) , or a dictionary of pre -found cache files ( one entry for each detector ) if not wanting to use gw_data_find e . g . { ' H1 ': '/ home / username / analysis / H1cache . txt '} datafind = / usr / bin / gw_data_find ; inputs for running a data find job ( NOTE : this can be used to search for frame if using the heterodyne , or SFTs if using spectral interpolation ) [ datafind ] # a dictionary of frame types to be returned ; one for each detector type = { ' H1 ': ' H1_HOFT_C00 ' , 'L1 ': ' L1_HOFT_C00 '} # a string to match in the URL paths returned match = localhost ; inputs for running a science segment finding job [ segmentfind ] # path to segment database query script , or a dictionary of pre -found segment files ( one entry for each detector ) if not wanting to use segment finding script e . g . { ' H1 ': '/ home / username / analysis / H1segments . txt '} segfind = / usr / bin / l i g o l w _ s e g m e n t _ q u e r y _ d q s e g d b # path to ligolw_print ligolw_print = / usr / bin / ligolw_print # URL of segment database server server = https :// segments . ligo . org # the standard deviation threshold for removing outliers stddev_thresh = 3.5 # set to output the coarse heterodyne data in binary files ( if true then a binary input will be assumed for the fine heterodyne ) binary_output = True # gzip the coarse output files rather than outputting as binary g z i p _ c o a r s e _ o u t p u t = False # gzip the fine output files g z i p _ f i n e _ o u t p u t = True ; inputs for running the l a l a p p s _ S p l I n t e r code [ splinter ] # condor universe universe = vanilla # path to SplInter executable splinter_exec = / usr / bin / l a l a p p s _ S p l I n t e r # list with the start and end frequency ranges for the SFTs freq_range = [30. , 2000.] # the standard deviation threshold for removing outliers stddev_thresh = 3.5 # the bandwidth ( Hz ) around the signal frequency to use in interpolation bandwidth = 0.3 # this is the default value from the code # minimum length ( seconds ) of science segments to use min_ seg_leng th = 1800 # this is the default value ( half an hour ) from the code ; inputs for running the parameter estimation code l a l a p p s _ p u l s a r _ p a r a m e t e r _ e s t i m a t i o n _ n e s t e d [ pe ] # condor universe universe = vanilla # path to the parameter estimation executable pe_exec = / usr / bin / l a l a p p s _ p u l s a r _ p a r a m e t e r _ e s t i m a t i o n _ n e s t e d # the base output directory for the nested samples ( directories for each detector / joint analysis will be created ) pe_output_dir = / home / username / analyses / nest ed_sampl es # a pre -made prior file for the analysis ( if this is not set then the prior file will be generated for each source based on the par file and / or other information ) -this has priority over any other options in this configuration file p r e m a d e _ p r i o r _ f i l e = # a set of prior options in dictionary form e . g . # prior_options = { ' H0 ': { ' priortype ': ' uniform ' , ' ranges ': [0. , 1e -22]} , ' PHI0 ': { ' priortype ': ' uniform ' , ' ranges ': [0. , 3.14]}} prior_options = # if ' derive_amplitude_prior ' is true ( and the ' premade_prior_file ' is not set ) then amplitude priors will be derived # from the heterodyned / spectrally interpolated data based on the either previous upper limits or those # estimated from ASDs of previous runs . # The ' amplitude_prior_type ' can also be set to either ' fermidirac ' ( default ) or ' uniform '. All other parameters # will have their ranges taken from ' predefined_prior_options '. # go through the pulsar parameter file and use the errors to set Gaussian priors in the prior file ( also create a correlation coefficient file for these , setting everything to be uncorrelated except the extremely highly correlated binary parameters ) u s e _ p a r a m e t e r _ e r r o r s = False # the number of parallel l a l a p p s _ p u l s a r _ p a r a m e t e r _ e s t i m a t i o n _ n e s t e d runs for a given pulsar / detector combination n_runs = 5 # the number of live points for each run n_live = 2048 # the number of MCMC samples for each nested sample update ( if not set this will be automatically calculated ) n_mcmc = # the number of MCMC samples for initial shuffling of the prior points ( shouldn ' t be needed for pulsar code ) n_mc mc_initi al = # the tolerance ( stopping criterion ) for the runs tolerance = 0.1 # a random seed for the RNG ( if not set then this will default to it ' s standard method ) random_seed = # flag to set whether running with non -GR p a ra m e t e r i s a t i o n non_gr = False # flag to say whether to use the ' waveform ' or ' source ' p a ra m e t e r i s a t i o n ( defaults to ' waveform ') model_type = ' waveform ' # flag to set whether using a Gaussian likelihood , or the default Student ' s t -likelihood ( if using Splinter for the pre -processing engine then a Gaussian likelihood will automatically be used , and override anything set here ) gaussian_like = False # the number of of parallel runs for each n_ r u n s _ b a c k g r o u n d = 2 # the base output directory for the background nested samples ( directories for each detector / joint analysis will be created ) p e _ o u t p u t _ d i r _ b a c k g r o u n d = / home / username / analyses / background / n ested_sa mples # the base output directory for the background posteriors n 2 p _ o u t p u t _ d i r _ b a c k g r o u n d = / home / username / analyses / background / p o s t e r i o r _ s a m p l e s # flag to set whether to clean ( i . e . remove ) all the nested sample files ( keeping the posteriors ) c l e a n _ n e s t _ s a m p l e s = False # set to true if wanting the output to use the l = m =2 gravitational wave phase as the initial phase , rather than the default rotational phase use_gw_phase = False # flags for use of Reduced Order Quadrature ( ROQ ) # Any generated ROQ interpolant files will be placed in structure of the pe_output_dir use_roq = False # set the number of training sets for using the in reduced basis and interpolant generation roq_ntraining = 2500 # set the maximum data chunk length for when using ROQ roq_chunkmax = 1440 # set the tolerance for producing the ROQ bases roq_tolerance = 5e -12 # set if wanting the bases produced over a uniform parameter space even for parameters with Gaussian priors roq_uniform = False ; inputs for creating results pages [ results_page ] # condor universe universe = local # results page ( Condor ) log directory log_dir = / usr1 / username / logs # results page creation executable results_exec = / usr / bin / l a l a p p s _ k n o p e _ r e s u l t _ p a g e . py # results collation executable collate_exec = / usr / bin / l a l a p p s _ k n o p e _ c o l l a t e _ r e s u l t s . py # the output base web directory for the results web_dir = / home / username / public_html / results # the equivalent output base URL for the above path base_url = https :// myurl /~username / results  ' , ' f1rot ' , 'ra ' , 'dec ' , ' dist ' , ' sdlim '] # set whether to show posterior plots for all parameters s h o w _ a l l _ p o s t e r i o r s = False # set whether to subtract the true / heterodyned value from any phase parameters in a search for plotting s ub tr ac t _t ru t hs = False # set whether to show the priors on the 1 D posterior plots show_priors = True # set whether to copy par file , prior files , heterodyne files , and posterior files into results page directory copy _all_fil es = True Unless specified as otherwise within a .par file the pipeline assumes the use of the TCB timing convention (as is the default for TEMPO2), and the DE405 solar system ephemeris.

B. DERIVATION OF THE STUDENT'S T-LIKELIHOOD FUNCTION
Here we derive, in detail, the form of the Student's t-likelihood given in Section 2.2. This derivation can, in part, be found in Dupuis (2004) and Dupuis and Woan (2005), but we correct a slight error from those references.
If a stretch of data of length m is assumed to consist of a signal defined by a set of parameters θ and Gaussian noise with zero mean and standard deviation σ, then the standard deviation can be included as an unknown parameter in the likelihood function and marginalised over. The likelihood can therefore be given by where p(B| θ, σ, I) is the likelihood function for the data given the unknown signal parameters and the noise standard deviation, and p(σ|I) is the prior on σ.
Given that the noise is assumed Gaussian the likelihood within the integral is given by p(B| θ, σ, I) = 1 2πσ 2 m exp − m k=1 |B k − y( θ) k | 2 2σ 2 .
for our complex data B and model y, where for a complex number x we have used |x| 2 = (x) 2 + (x) 2 , and the noise in the real and imaginary parts is assumed to have the same distribution.
We chose a scale invariant prior on σ of p(σ|I) = 1/σ, for which the marginalisation of Eqn. B1 can be performed analytically as we will show. Substituting this prior and Eqn. B2 into Eqn. B1 we can begin the integration with the substitution giving du dσ = m k=1 |B k − y( θ) k | 2 2 σ −2 .
Putting everything into Eqn. B1 gives p(B| θ) =(2π) −m ∞ 0 2 (m+ 1 2 ) X −(m+ 1 2 ) u 2m+1 e −u 2 2 −1/2 X 1/2 u −2 du, The integral in Eqn. B7 can then be performed by making the substitutions x = u 2 and dx = 2u du, giving This integral is the definition of the Gamma function For positive integers, Stirling's formula shows that the Gamma function is given by Γ(n + 1) = n!, so our likelihood becomes Note that, again, this differs from Eqn. 2.31 in Dupuis (2004) by some constant factors. For our analysis the data is split into chunks for each of which the assumption of noise stationarity is thought to be valid, but between chunks is thought to be invalid (see discussion in §2.4). The joint likelihood of all data can be obtained from the product of the likelihood for each chunk given by where M is the total number of independent data chunks with lengths m j and k 0 = 1 + j i=1 m i−1 (with m 0 = 0) is the index of the first data point in each chunk.

C. FAST LIKELIHOOD EVALUATION
In §2.2.4 we stated that the likelihood can be evaluated quickly in cases where the ∆φ K values in Equations 7, 8 and 11 are zero, i.e., no phase evolution parameters are required. This can be done by pre-summing over data, B, and combinations of the time-varying antenna response function, after which only coefficients of these pre-summed terms need to be calculated. Here we will explicitly write out the derivation of these coefficients for the case of a signal in GR, for which the antenna pattern functions for the '+' and '×' components are given by Equation 16.
Assuming a signal in GR, we start with a model of the form h(t) =C + F + (ψ, t)e iΦ lm + iC × F × (ψ, t)e iΦ lm , =C + F + (ψ, t) cos Φ lm − C × F × (ψ, t) sin Φ lm + i (C + F + (ψ, t) sin Φ lm + C × F × (ψ, t) cos Φ lm ) , =C R,+ F + (ψ, t) + C R,× F × (ψ, t) + i (C I,+ F + (ψ, t) + C I,× F × (ψ, t)) , where, setting ζ = 90 • in Equation 16, the antenna pattern functions are F + (ψ, t) = a(t) cos 2ψ + b(t) sin 2ψ, F × (ψ, t) = b(t) cos 2ψ − a(t) sin 2ψ, and we have set the real and imaginary model amplitude coefficients for the two polarisations to be C R,+ = C + cos Φ lm , C R,× = −C × sin Φ lm , C I,+ = C + sin Φ lm , C I,× = C × cos Φ lm . 34 If we just take the likelihood (either the Student's 34 For our GR signal model emitting at twice the rotation frequency, given by Equation 8, we would have C + = − C 22 2 1 + cos 2 ι and C × = C 22 cos ι. t-likelihood of Equation 12, or the Gaussian likelihood of Equation 13) 35 for a single detector, data stream, and data chunk (see §2.2) then we see that it contains the summation over time samples We can expand this out to give where for convenience we have removed the summation subscripts and explicit time dependence, and used R and I subscripts to represent the real and imaginary components respectively. We are now interested in the summations that contain parts of the model function h. If we take the summation over the square of the real part of the model we get R,+ a 2 cos 2 2ψ + b 2 sin 2 2ψ + 2ab sin 2ψ cos 2ψ + C 2 R,× b 2 cos 2 2ψ + a 2 sin 2 2ψ − 2ab sin 2ψ cos 2ψ + 2C R,+ C R,× (b 2 − a 2 ) sin 2ψ cos 2ψ + ab(cos 2 2ψ − sin 2 2ψ) , where the summation is now implicit in the a and b terms, i.e., a 2 = a(t) 2 , b 2 = b(t) 2 , and ab = a(t)b(t), which are the only time varying components and can be pre-computed. If we now say where the K R s are the coefficients of a 2 , b 2 and 2ab, we find them to be K R a 2 = C 2 R,+ cos 2 2ψ + C 2 R,× sin 2 2ψ − 2C R,+ C R,× sin 2ψ cos 2ψ, K R b 2 = C 2 R,+ sin 2 2ψ + C 2 R,× cos 2 2ψ + 2C R,+ C R,× sin 2ψ cos 2ψ, K R 2ab = C R,+ C R,× cos 2 2ψ − sin 2 2ψ + sin 2ψ cos 2ψ C 2 R,+ − C 2 R,× .
Factorising these out we find that K R a = C R,+ cos 2ψ − C R,× sin 2ψ, K R b = C R,+ sin 2ψ + C R,× cos 2ψ, where K R a 2 = K R a 2 , K R b 2 = K R b 2 , and K R 2ab = K R a K R b . It is easy to see from above that things will be identical for the imaginary components and the equations hold just by swapping R for I. It is worth noting that the coefficients given in Equation C8 take the opposite form (covariant vs. contravariant rotation) than the antenna pattern functions given in Equation C2, and this difference is present in the model in the lppen code, but should not be unexpected.
We can also see that the coefficients given in Equation C8 are those required in the components of Equation C5 that sum over the product of the data and model. For the real part we have where we have set d a R = B R a and d b R = B R b. We can therefore see that the coefficients of d a R and d b R are indeed K R a and K R b respectively, i.e., Again, this holds for the imaginary parts by just swapping R for I. The same process can be applied for likelihood evaluations for the non-GR polarisation modes. For the 'x' and 'y' modes the antenna patterns can be defined as (Equations 32 and 33 of Isi et al. 2015) F x (ψ, t) = A x (t) cos ψ + A y sin ψ, F y (ψ, t) = A y (t) cos ψ − A x sin ψ.
In this case we find, in an identical fashion to above, that the required coefficients of A x and A y are K R/I Ax = C R/I,x cos ψ − C R/I,y sin ψ, K R/I Ay = C R/I,y cos ψ + C R/I,x sin ψ, where C R/I,x/y are the real/imaginary signal amplitude components for the two modes. Finally, for the breathing and longitudinal modes, where the antenna patterns are defined as (Equations 34 and 35 of Isi et al. 2015) and there is no ψ dependence, the coefficients of A b and A l are just the real/imaginary signal amplitude coefficients C R/I,b/l for the modes.

D. RELATION BETWEEN ODDS AND SNRS
It is interesting to look at how we might expect odds values (ratios of evidences) to vary with signal-to-noise ratio (SNR) (also see similar discussions in Isi et al. 2017). We will take the simple case of two Gaussian datasets both of the same length (N = 100) and standard deviation, σ. Both datasets contains a mean offset of µ h = 5, and we are concerned with a search over the mean parameter. We want to calculate two things: the odds for a mean offset versus the data being Gaussian noise with a known mean of zero; and, the odds for both datasets containing the same mean offset (a coherent signal) versus them containing two independent mean offsets (an incoherent signal). We can try and work out what to expect using approximations, and then see if that's actually the case when calculating the values analytically.
To set up the situation, we have the two datasets d 1 and d 2 , which have likelihoods given by p(d j |µ, M, I) = 1 2πσ 2 and priors on µ given by p(µ|M, I), which we will take as constant (i.e. a flat prior). We will also use the fact that d j,i = µ t + n j,i where n j,i is a value drawn from a Gaussian distribution of zero mean and standard deviation σ, and µ t is the true offset.
To re-cap from earlier in this paper, the evidence is given by where p(d|µ = µ t/ML , M, I) is the likelihood at the true/maximum likelihood value of µ. This value, for the joint datasets, can be seen to be which in turn (with d j,i = µ t + n j,i ) becomes where we have used N i n 2 i /σ 2 ≈ N ± √ N . Similarly, we have the noise evidence Z N = 1 2πσ 2 N exp − 1 2σ 2 N i 2µ 2 t + (n 2 1,i + n 2 2,i ) − 2µ t (n 1,i + n 2,i ) where we use the approximation that N i µ t n j,i ≈ 0. So, the odds becomes where we can also say that the single dataset SNR is defined as ρ 2 = N i µ 2 t /σ 2 = N (µ t /σ) 2 , and in this case ρ coh = √ 2ρ, and we can see that if we take the natural logarithm of this that ln O S/N = ln (p(µ|M, I)) + 1 2 ln 2πσ 2 N + 1 2 ρ 2 coh .
We see here that the log odds is proportional to ρ 2 , with a minor influence from the ln 2πσ 2 /N term which is roughly equivalent to ln 2πρ −2 . Empirically (and probably with some mathematical justification that can be derived) this approximation seems to give the median odds value on random realisations of noise, even down to low SNR (see Figure 34).

D.2. Coherent signal vs. incoherent signal
We can define a coherent signal versus incoherent signal odds ratio as where Z {d} , Z d1 , and Z d2 , are the evidences for a signal in the full dataset, and two individual datasets (this is equivalent to Equation 31). Using the same approximations as above, we find that O S/I simple ≈ 1 2 √ π 1 p(µ|M, I) In this we see that there is the Occam factor of 1/p(µ|M, I) accounting for the fact that the incoherent model has an extra parameter. We also see that this is proportional to √ N /σ which is roughly equivalent to ρ. So we see that for the coherent signal versus noise odds roughly scales with e ρ 2 , whilst the coherent signal versus incoherent signal odds roughly scales with ρ, i.e. its scaling is a lot shallower. Empirically (and probably with some mathematical justification that can be derived) this approximation seems to give the maximum odds value on random realisations of noise down to low SNR (see Figure 34).
In Figure 34 we see how these two approximations compare to analytical calculations of the odds values as the SNR changes. For a range of SNRs we created multiple realisations of two dataset of Gaussian random noise (with SNR altered by altering the standard deviation of the noise) with a given mean offset. In each case we have analytically calculated the evidences for a coherent µ for both datasets, and independent µs in both datasets. These are then compared to the approximations above. In this case we have purposely chosen a prior on µ such that the distributions of odds values for O S/I simple and O S/N cross, which is similar to the case we observed in our actual analysis shown in §4.1.2.
We are grateful for computational resources provided by Cardiff University, and funded by an STFC grant (ST/I006285/1) supporting UK Involvement in the Operation of Advanced LIGO. MP is funded by the STFC under grant number ST/N005422/1. JV is funded by the STFC under grant number ST/K005014/1. The versions of lalapps pulsar parameter estimation nested and lalapps pulsar parameter estimation used in the tests for this document can be found in the lalapps knope O2 branch of the LALSuite git repository with the git hash fe9b7d4. The various corner plots, showing posterior parameter distributions, in this document have been produced using https://github.com/mattpitkin/scotchcorner, and all plots have used Matplotlib Hunter (2007); Droettboom et al. (2017). This is expected to be a living document that will be updated with any changes to the code. The document has LIGO DCC number LIGO-P1700086.