
We present highprecision timing data over time spans of up to 11 years for
45 millisecond pulsars observed as part of the North American Nanohertz
Observatory for Gravitational Waves (NANOGrav) project, aimed at detecting and
characterizing lowfrequency gravitational waves. The pulsars were observed
with the Arecibo Observatory and/or the Green Bank Telescope at frequencies
ranging from 327 MHz to 2.3 GHz. Most pulsars were observed with approximately
monthly cadence, with six hightimingprecision pulsars observed weekly, and
all were observed at widely separated frequencies at each observing epoch in
order to fit for timevariable dispersion delays. We describe our methods for
data processing, timeofarrival (TOA) calculation, and the implementation of a
new, automated method for removing outlier TOAs. We fit a timing model for each
pulsar that includes spin, astrometric, and, if necessary, binary parameters,
in addition to timevariable dispersion delays and parameters that quantify
pulseprofile evolution with frequency. The new timing solutions provide three
new parallax measurements, two new Shapiro delay measurements, and two new
measurements of large orbitalperiod variations. We fit models that
characterize sources of noise for each pulsar. We find that 11 pulsars show
significant red noise, with generally smaller spectral indices than typically
measured for nonrecycled pulsars, possibly suggesting a different origin.
Future papers will use these data to constrain or detect the signatures of
gravitationalwave signals.

Pulsartiming datasets have been analyzed with great success using
probabilistic treatments based on Gaussian distributions, with applications
ranging from studies of neutronstar structure to tests of general relativity
and searches for nanosecond gravitational waves. As for other applications of
Gaussian distributions, outliers in timing measurements pose a significant
challenge to statistical inference, since they can bias the estimation of
timing and noise parameters, and affect reported parameter uncertainties. We
describe and demonstrate a practical endtoend approach to perform Bayesian
inference of timing and noise parameters robustly in the presence of outliers,
and to identify these probabilistically. The method is fully consistent (i.e.,
outlierness probabilities vary in tune with the posterior distributions of the
timing and noise parameters), and it relies on the efficient sampling of the
hierarchical form of the pulsartiming likelihood. Such sampling has recently
become possible with a "noUturn" Hamiltonian sampler coupled to a highly
customized reparametrization of the likelihood; this code is described
elsewhere, but it is already available online. We recommend our method as a
standard step in the preparation of pulsartimingarray datasets: even if
statistical inference is not affected, followup studies of outlier candidates
can reveal unseen problems in radio observations and timing measurements;
furthermore, confidence in the results of gravitationalwave searches will only
benefit from stringent statistical evidence that datasets are clean and
outlierfree.

We have searched for continuous gravitational wave (CGW) signals produced by
individually resolvable, circular supermassive black hole binaries (SMBHBs) in
the latest EPTA dataset, which consists of ultraprecise timing data on 41
millisecond pulsars. We develop frequentist and Bayesian detection algorithms
to search both for monochromatic and frequencyevolving systems. None of the
adopted algorithms show evidence for the presence of such a CGW signal,
indicating that the data are best described by pulsar and radiometer noise
only. Depending on the adopted detection algorithm, the 95\% upper limit on the
skyaveraged strain amplitude lies in the range $6\times
10^{15}<A<1.5\times10^{14}$ at $5{\rm nHz}<f<7{\rm nHz}$. This limit varies
by a factor of five, depending on the assumed source position, and the most
constraining limit is achieved towards the positions of the most sensitive
pulsars in the timing array. The most robust upper limit  obtained via a full
Bayesian analysis searching simultaneously over the signal and pulsar noise on
the subset of ours six best pulsars  is $A\approx10^{14}$. These limits, the
most stringent to date at $f<10{\rm nHz}$, exclude the presence of
subcentiparsec binaries with chirp mass $\cal{M}_c>10^9$M$_\odot$ out to a
distance of about 25Mpc, and with $\cal{M}_c>10^{10}$M$_\odot$ out to a
distance of about 1Gpc ($z\approx0.2$). We show that stateoftheart SMBHB
population models predict $<1\%$ probability of detecting a CGW with the
current EPTA dataset, consistent with the reported nondetection. We stress,
however, that PTA limits on individual CGW have improved by almost an order of
magnitude in the last five years. The continuing advances in pulsar timing data
acquisition and analysis techniques will allow for strong astrophysical
constraints on the population of nearby SMBHBs in the coming years.

We present new limits on an isotropic stochastic gravitationalwave
background (GWB) using a six pulsar dataset spanning 18 yr of observations from
the 2015 European Pulsar Timing Array data release. Performing a Bayesian
analysis, we fit simultaneously for the intrinsic noise parameters for each
pulsar, along with common correlated signals including clock, and Solar System
ephemeris errors, obtaining a robust 95$\%$ upper limit on the dimensionless
strain amplitude $A$ of the background of $A<3.0\times 10^{15}$ at a reference
frequency of $1\mathrm{yr^{1}}$ and a spectral index of $13/3$, corresponding
to a background from inspiralling supermassive black hole binaries,
constraining the GW energy density to $\Omega_\mathrm{gw}(f)h^2 <
1.1\times10^{9}$ at 2.8 nHz. We also present limits on the correlated power
spectrum at a series of discrete frequencies, and show that our sensitivity to
a fiducial isotropic GWB is highest at a frequency of $\sim 5\times10^{9}$~Hz.
Finally we discuss the implications of our analysis for the astrophysics of
supermassive black hole binaries, and present 95$\%$ upper limits on the string
tension, $G\mu/c^2$, characterising a background produced by a cosmic string
network for a set of possible scenarios, and for a stochastic relic GWB. For a
NambuGoto field theory cosmic string network, we set a limit
$G\mu/c^2<1.3\times10^{7}$, identical to that set by the {\it Planck}
Collaboration, when combining {\it Planck} and high$\ell$ Cosmic Microwave
Background data from other experiments. For a stochastic relic background we
set a limit of $\Omega^\mathrm{relic}_\mathrm{gw}(f)h^2<1.2 \times10^{9}$, a
factor of 9 improvement over the most stringent limits previously set by a
pulsar timing array.

We compute upper limits on the nanohertzfrequency isotropic stochastic
gravitational wave background (GWB) using the 9year data release from the
North American Nanohertz Observatory for Gravitational Waves (NANOGrav)
collaboration. We set upper limits for a GWB from supermassive black hole
binaries under power law, broken power law, and free spectral coefficient GW
spectrum models. We place a 95\% upper limit on the strain amplitude (at a
frequency of yr$^{1}$) in the power law model of $A_{\rm gw} < 1.5\times
10^{15}$. For a broken power law model, we place priors on the strain
amplitude derived from simulations of Sesana (2013) and McWilliams et al.
(2014). We find that the data favor a broken power law to a pure power law with
odds ratios of 22 and 2.2 to one for the McWilliams and Sesana prior models,
respectively. The McWilliams model is essentially ruled out by the data, and
the Sesana model is in tension with the data under the assumption of a pure
power law. Using the broken powerlaw analysis we construct posterior
distributions on environmental factors that drive the binary to the GWdriven
regime including the stellar mass density for stellarscattering, mass
accretion rate for circumbinary disk interaction, and orbital eccentricity for
eccentric binaries, marking the first time that the shape of the GWB spectrum
has been used to make astrophysical inferences. We then place the most
stringent limits so far on the energy density of relic GWs,
$\Omega_\mathrm{gw}(f)\,h^2 < 4.2 \times 10^{10}$, yielding a limit on the
Hubble parameter during inflation of $H_*=1.6\times10^{2}~m_{Pl}$, where
$m_{Pl}$ is the Planck mass. Our limit on the cosmic string GWB,
$\Omega_\mathrm{gw}(f)\, h^2 < 2.2 \times 10^{10}$, translates to a
conservative limit of $G\mu<3.3\times 10^{8}$  a factor of 4 better than the
joint Planck and high$l$ CMB data from other experiments.

We extend the formalisms developed in Gair et al. and Cornish and van
Haasteren to create maps of gravitationalwave backgrounds using a network of
groundbased laser interferometers. We show that in contrast to pulsar timing
arrays, which are insensitive to half of the gravitationalwave sky (the curl
modes), a network of groundbased interferometers is sensitive to both the
gradient and curl components of the background. The spatial separation of a
network of interferometers, or of a single interferometer at different times
during its rotational and orbital motion around the Sun, allows for recovery of
both components. We derive expressions for the response functions of a laser
interferometer in the smallantenna limit, and use these expressions to
calculate the overlap reduction function for a pair of interferometers. We also
construct maximumlikelihood estimates of the + and xpolarization modes of the
gravitationalwave sky in terms of the response matrix for a network of
groundbased interferometers, evaluated at discrete times during Earth's
rotational and orbital motion around the Sun. We demonstrate the feasibility of
this approach for some simple simulated backgrounds (a single point source and
spatiallyextended distributions having only grad or curl components),
calculating maximumlikelihood sky maps and uncertainty maps based on the
(pseudo)inverse of the response matrix. The distinction between this approach
and standard methods for mapping gravitationalwave power is also discussed.

Many dataanalysis problems involve large dense matrices that describe the
covariance of stationary noise processes; the computational cost of inverting
these matrices, or equivalently of solving linear systems that contain them, is
often a practical limit for the analysis. We describe two general, practical,
and accurate methods to approximate stationary covariance matrices as lowrank
matrix products featuring carefully chosen spectral components. These methods
can be used to greatly accelerate dataanalysis methods in many contexts, such
as the Bayesian and generalizedleastsquares analysis of pulsartiming
residuals.

In this work we review the application of the theory of Gaussian processes to
the modeling of noise in pulsartiming data analysis, and we derive various
useful and optimized representations for the likelihood expressions that are
needed in Bayesian inference on pulsartimingarray datasets. The resulting
viewpoint and formalism lead us to two improved parametersampling schemes
inspired by Gibbs sampling. The new schemes have vastly lower chain
autocorrelation lengths than the Markov Chain Monte Carlo methods currently
used in pulsartiming data analysis, potentially speeding up Bayesian inference
by orders of magnitude. The new schemes can be used for a fullnoisemodel
analysis of the large datasets assembled by the International Pulsar Timing
Array collaboration, which present a serious computational challenge to
existing methods.

We describe a new method for extracting gravitational wave signals from
pulsar timing data. We show that any gravitational wave signal can be
decomposed into an orthogonal set of sky maps, with the number of maps equal to
the number of pulsars in the timing array. These maps may be used as a basis to
construct gravitational wave templates for any type of source, including
collections of point sources. A variant of the standard HellingsDowns
correlation analysis is recovered for statistically isotropic signals. The
template based approach allows us to probe potential anisotropies in the signal
and produce maps of the gravitational wave sky.

A new Bayesian software package for the analysis of pulsar timing data is
presented in the form of TempoNest which allows for the robust determination of
the nonlinear pulsar timing solution simultaneously with a range of additional
stochastic parameters. This includes both red spin noise and dispersion measure
variations using either power law descriptions of the noise, or through a
modelindependent method that parameterises the power at individual frequencies
in the signal. We use TempoNest to show that at noise levels representative of
current datasets in the European Pulsar Timing Array (EPTA) and International
Pulsar Timing Array (IPTA) the linear timing model can underestimate the
uncertainties of the timing solution by up to an order of magnitude. We also
show how to perform Bayesian model selection between different sets of timing
model and stochastic parameters, for example, by demonstrating that in the
pulsar B1937+21 both the dispersion measure variations and spin noise in the
data are optimally modelled by simple power laws. Finally we show that not
including the stochastic parameters simultaneously with the timing model can
lead to unpredictable variation in the estimated uncertainties, compromising
the robustness of the scientific results extracted from such analysis.

A new model independent method is presented for the analysis of pulsar timing
data and the estimation of the spectral properties of an isotropic
gravitational wave background (GWB). We show that by rephrasing the likelihood
we are able to eliminate the most costly aspects of computation normally
associated with this type of data analysis. When applied to the International
Pulsar Timing Array Mock Data Challenge data sets this results in speedups of
approximately 2 to 3 orders of magnitude compared to established methods. We
present three applications of the new likelihood. In the low signal to noise
regime we sample directly from the power spectrum coefficients of the GWB
signal realization. In the high signal to noise regime, where the data can
support a large number of coefficients, we sample from the joint probability
density of the power spectrum coefficients for the individual pulsars and the
GWB signal realization. Critically in both these cases we need make no
assumptions about the form of the power spectrum of the GWB, or the individual
pulsars. Finally we present a method for characterizing the spatial correlation
between pulsars on the sky, making no assumptions about the form of that
correlation, and therefore providing the only truly general Bayesian method of
confirming a GWB detection from pulsar timing data.

Direct detection of gravitational waves by pulsar timing arrays will become
feasible over the next few years. In the low frequency regime ($10^{7}$ Hz 
$10^{9}$ Hz), we expect that a superposition of gravitational waves from many
sources will manifest itself as an isotropic stochastic gravitational wave
background. Currently, a number of techniques exist to detect such a signal;
however, many detection methods are computationally challenging. Here we
introduce an approximation to the full likelihood function for a pulsar timing
array that results in computational savings proportional to the square of the
number of pulsars in the array. Through a series of simulations we show that
the approximate likelihood function reproduces results obtained from the full
likelihood function. We further show, both analytically and through
simulations, that, on average, this approximate likelihood function gives
unbiased parameter estimates for astrophysically realistic stochastic
background amplitudes.

This is a summary of the methods we used to analyse the first IPTA Mock Data
Challenge (MDC), and the obtained results. We have used a Bayesian analysis in
the time domain, accelerated using the recently developed ABCmethod which
consists of a form of lossy linear data compression. The TOAs were first
processed with Tempo2, where the design matrix was extracted for use in a
subsequent Bayesian analysis. We used different noise models to analyse the
datasets: no red noise, red noise the same for all pulsars, and individual red
noise per pulsar. We sampled from the likelihood with four different samplers:
"emcee", "twalk", "MetropolisHastings", and "pyMultiNest". All but emcee
agreed on the final result, with emcee failing due to artefacts of the
highdimensionality of the problem. An interesting issue we ran into was that
the prior of all the 36 (red) noise amplitudes strongly affects the results. A
flat prior in the noise amplitude biases the inferred GWB amplitude, whereas a
flat prior in logamplitude seems to work well. This issue is only apparent
when using a noise model with individually modelled red noise for all pulsars.
Our results for the blind challenges are in good agreement with the injected
values. For the GWB amplitudes we found h_c = 1.03 +/ 0.11 [10^{14}], h_c =
5.70 +/ 0.35 [10^{14}], and h_c = 6.91 +/ 1.72 [10^{15}], and for the GWB
spectral index we found gamma = 4.28 +/ 0.20, gamma = 4.35 +/ 0.09, and gamma
= 3.75 +/ 0.40. We note that for closed challenge 3 there was quite some
covariance between the signal and the red noise: if we constrain the GWB
spectral index to the usual choice of gamma = 13/3, we obtain the estimates:
h_c = 10.0 +/ 0.64 [10^{15}], h_c = 56.3 +/ 2.42 [10^{15}], and h_c = 4.83
+/ 0.50 [10^{15}], with onesided 2 sigma upperlimits of: h_c <= 10.98
[10^{15}], h_c <= 60.29 [10^{15}], and h_c <= 5.65 [10^{15}].

The analysis of pulsar timing data, especially in pulsar timing array (PTA)
projects, has encountered practical difficulties: evaluating the likelihood
and/or correlationbased statistics can become prohibitively computationally
expensive for large datasets. In situations where a stochastic signal of
interest has a power spectral density that dominates the noise in a limited
bandwidth of the total frequency domain (e.g. the isotropic background of
gravitational waves), a linear transformation exists that transforms the timing
residuals to a basis in which virtually all the information about the
stochastic signal of interest is contained in a small fraction of basis
vectors. By only considering such a small subset of these "generalised
residuals", the dimensionality of the data analysis problem is greatly reduced,
which can cause a large speedup in the evaluation of the likelihood: the
ABCmethod (Acceleration By Compression). The compression fidelity, calculable
with crude estimates of the signal and noise, can be used to determine how far
a dataset can be compressed without significant loss of information. Both
direct tests on the likelihood, and Bayesian analysis of mock data, show that
the signal can be recovered as well as with an analysis of uncompressed data.
In the analysis of IPTA Mock Data Challenge datasets, speedups of a factor of
three orders of magnitude are demonstrated. For realistic PTA datasets the
acceleration may become greater than six orders of magnitude due to the low
signal to noise ratio.

Although it is widely understood that pulsar timing observations generally
contain timecorrelated stochastic signals (TCSSs; red timing noise is of this
type), most data analysis techniques that have been developed make an
assumption that the stochastic uncertainties in the data are uncorrelated, i.e.
"white". Recent work has pointed out that this can introduce severe bias in
determination of timingmodel parameters, and that better analysis methods
should be used. This paper presents a detailed investigation of timingmodel
fitting in the presence of TCSSs, and gives closed expressions for the postfit
signals in the data. This results in a Bayesian technique to obtain
timingmodel parameter estimates in the presence of TCSSs, as well as
computationally more efficient expressions of their marginalised posterior
distribution. A new method to analyse hundreds of mock dataset realisations
simultaneously without significant computational overhead is presented, as well
as a statistically rigorous method to check the internal consistency of the
results. As a byproduct of the analysis, closed expressions of the rms
introduced by a stochastic background of gravitationalwaves in
timingresiduals are obtained. Using $T$ as the length of the dataset, and
$h_c(1\rm{yr}^{1})$ as the characteristic strain, this is: $\sigma_{\rm GWB}^2
= h_{c}(1\rm{yr}^{1})^2 (9\sqrt[3]{2\pi^4}\Gamma(10/3) / 8008) \rm{yr}^{4/3}
T^{10/3}$.

Markov Chain Monte Carlo (MCMC) methods have revolutionised Bayesian data
analysis over the years by making the direct computation of posterior
probability densities feasible on modern workstations. However, the calculation
of the prior predictive, the Bayesian evidence, has proved to be notoriously
difficult with standard techniques. In this work a method is presented that
lets one calculate the Bayesian evidence using nothing but the results from
standard MCMC algorithms, like MetropolisHastings. This new method is compared
to other methods like MultiNest, and greatly outperforms the latter in several
cases. One of the toy problems considered in this work is the analysis of mock
pulsar timing data, as encountered in pulsar timing array projects. This method
is expected to be useful as well in other problems in astrophysics, cosmology
and particle physics.

Pulsar timing arrays (PTAs) are designed to detect gravitational waves with
periods from several months to several years, e.g. those produced by by wide
supermassive blackhole binaries in the centers of distant galaxies. Here we
show that PTAs are also sensitive to mergers of supermassive black holes. While
these mergers occur on a timescale too short to be resolvable by a PTA, they
generate a change of metric due to nonlinear gravitationalwave memory which
persists for the duration of the experiment and could be detected. We develop
the theory of the singlesource detection by PTAs, and derive the sensitivity
of PTAs to the gravitationalwave memory jumps. We show that mergers of
$10^8M_{\odot}$ black holes are $2\sigma$detectable (in a direction,
polarization, and timedependent way) out to comoving distances of $\sim 1$
billion light years. Modern prediction for blackhole merger rates imply
marginal to modest chance of an individual jump detection by currently
developed PTAs. The sensitivity is expected to be somewhat higher for
futuristic PTA experiments with SKA.

Longterm precise timing of Galactic millisecond pulsars holds great promise
for measuring the longperiod (monthstoyears) astrophysical gravitational
waves. Several gravitationalwave observational programs, called Pulsar Timing
Arrays (PTA), are being pursued around the world. Here we develop a Bayesian
algorithm for measuring the stochastic gravitationalwave background (GWB) from
the PTA data. Our algorithm has several strengths: (1) It analyses the data
without any loss of information, (2) It trivially removes systematic errors of
known functional form, including quadratic pulsar spindown, annual modulations
and jumps due to a change of equipment, (3) It measures simultaneously both the
amplitude and the slope of the GWB spectrum, (4) It can deal with unevenly
sampled data and coloured pulsar noise spectra. We sample the likelihood
function using Markov Chain Monte Carlo (MCMC) simulations. We extensively test
our approach on mock PTA datasets, and find that the algorithm has significant
benefits over currently proposed counterparts. We show the importance of
characterising all red noise components in pulsar timing noise by demonstrating
that the presence of a red component would significantly hinder a detection of
the GWB. Lastly, we explore the dependence of the signaltonoise ratio on the
duration of the experiment, number of monitored pulsars, and the magnitude of
the pulsar timing noise. These parameter studies will help formulate observing
strategies for the PTA experiments.