11institutetext:  Institute of Astrophysics, FORTH, N. Plastira 100, 70013,
Heraklion, Greece22institutetext: Max-Planck-Institut für Radioastronomie, Auf
dem Hügel 69, 53121 Bonn, Germany33institutetext: Department of Physics,
Indian Institute of Technology Roorkee, Roorkee-247667, India44institutetext:
Department of Electrical Engineering, IIT Hyderabad, Kandi, Telangana 502284,
India 55institutetext: Cosmology, Universe and Relativity at Louvain (CURL),
Institute of Mathematics and Physics, University of Louvain, 2 Chemin du
Cyclotron, 1348 Louvain-la-Neuve, Belgium 66institutetext: Université Paris
Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France77institutetext:
The Institute of Mathematical Sciences, C. I. T. Campus, Taramani, Chennai
600113, India 88institutetext: Homi Bhabha National Institute, Training School
Complex, Anushakti Nagar, Mumbai 400094, India 99institutetext: Fakultät für
Physik, Universität Bielefeld, Postfach 100131, 33501 Bielefeld,
Germany1010institutetext: Scuola Internazionale Superiore di Studi Avanzati:
Via Bonomea 265, I-34136 Trieste, Italy and INFN Sezione di Trieste
1111institutetext: ASTRON, Netherlands Institute for Radio Astronomy, Oude
Hoogeveensedijk 4, 7991 PD, Dwingeloo, The Netherlands1212institutetext:
Department of Physical Sciences, Indian Institute of Science Education and
Research, Mohali, Punjab 140306, India 1313institutetext: Laboratoire de
Physique et Chimie de l’Environnement et de l’Espace, Université d’Orléans /
CNRS, 45071 Orléans Cedex 02, France 1414institutetext: Observatoire
Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université
d’Orléans, CNRS, 18330 Nançay, France1515institutetext: Dipartimento di Fisica
“G. Occhialini”, Universitá degli Studi di Milano-Bicocca, Piazza della
Scienza 3, I-20126 Milano, Italy1616institutetext: INFN, Sezione di Milano-
Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy1717institutetext: INAF
- Osservatorio Astronomico di Brera, via Brera 20, I-20121 Milano,
Italy1818institutetext: Institute for Gravitational Wave Astronomy and School
of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15
2TT, UK1919institutetext: INAF - Osservatorio Astronomico di Cagliari, via
della Scienza 5, 09047 Selargius (CA), Italy2020institutetext: Hellenic Open
University, School of Science and Technology, 26335 Patras,
Greece2121institutetext: Université de Genève, Département de Physique
Théorique and Centre for Astroparticle Physics, 24 quai Ernest-Ansermet,
CH-1211 Genéve 4, Switzerland2222institutetext: CERN, Theoretical Physics
Department, 1 Esplanade des Particules, CH-1211 Genéve 23,
Switzerland2323institutetext: Kavli Institute for Astronomy and Astrophysics,
Peking University, Beijing 100871, P. R. China2424institutetext: Department of
Astronomy and Astrophysics, Tata Institute of Fundamental Research, Homi
Bhabha Road, Navy Nagar, Colaba, Mumbai 400005, India 2525institutetext:
Department of Physics, IIT Hyderabad, Kandi, Telangana 502284, India
2626institutetext: Department of Physics and Astrophysics, University of
Delhi, Delhi 110007, India 2727institutetext: Department of Earth and Space
Sciences, Indian Institute of Space Science and Technology, Valiamala,
Thiruvananthapuram, Kerala 695547,India 2828institutetext: School of
Mathematics and Physics, Faculty of Engineering and Physical Science,
University of Surrey, Guildford GU2 7XH, UK2929institutetext: School of
Physics, Faculty of Science, University of East Anglia, Norwich NR4 7TJ,
UK3030institutetext: Sternberg Astronomical Institute, Moscow State
University, Universitetsky pr., 13, Moscow 119234, Russia3131institutetext:
Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am
Mühlenberg 1, 14476 Potsdam, Germany3232institutetext: Gran Sasso Science
Institute (GSSI), I-67100 L’Aquila, Italy 3333institutetext: INFN, Laboratori
Nazionali del Gran Sasso, I-67100 Assergi, Italy 3434institutetext: National
Centre for Radio Astrophysics, Pune University Campus, Pune 411007, India
3535institutetext: Kumamoto University, Graduate School of Science and
Technology, Kumamoto, 860-8555, Japan 3636institutetext: Universitá di
Cagliari, Dipartimento di Fisica, S.P. Monserrato-Sestu Km 0,700 - 09042
Monserrato (CA), Italy3737institutetext: Department of Astrophysics/IMAPP,
Radboud University Nijmegen, P.O. Box 9010, 6500 GL Nijmegen, The
Netherlands3838institutetext: Department of Physical Sciences,Indian Institute
of Science Education and Research Kolkata, Mohanpur, 741246, India
3939institutetext: Center of Excellence in Space Sciences India, Indian
Institute of Science Education and Research Kolkata, 741246, India
4040institutetext: School of Physics, Trinity College Dublin, College Green,
Dublin 2, D02 PN40, Ireland4141institutetext: Jodrell Bank Centre for
Astrophysics, Department of Physics and Astronomy, University of Manchester,
Manchester M13 9PL, UK4242institutetext: Department of Physics, St. Xavier’s
College (Autonomous), Mumbai 400001, India 4343institutetext: Department of
Astronomy,School of Physics, Peking University, Beijing 100871, P. R.
China4444institutetext: National Astronomical Observatories, Chinese Academy
of Sciences, Beijing 100101, P. R. China4545institutetext: E.A. Milne Centre
for Astrophysics, University of Hull, Cottingham Road, Kingston-upon-Hull, HU6
7RX, UK4646institutetext: Centre of Excellence for Data Science, Artificial
Intelligence and Modelling (DAIM), University of Hull, Cottingham Road,
Kingston-upon-Hull, HU6 7RX, UK4747institutetext: Laboratory of Astrophysics,
École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne,
Switzerland4848institutetext: Department of Physics, BITS Pilani Hyderabad
Campus, Hyderabad 500078, Telangana, India 4949institutetext: Joint Astronomy
Programme, Indian Institute of Science, Bengaluru, Karnataka, 560012, India
5050institutetext: Arecibo Observatory, HC3 Box 53995, Arecibo, PR 00612,
USA5151institutetext: IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-
Yvette, France 5252institutetext: Raman Research Institute India, Bengaluru,
Karnataka, 560080, India 5353institutetext: Institut für Physik und
Astronomie, Universität Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, 14476,
Potsdam, Germany5454institutetext: Kazan Federal University, 18 Kremlyovskaya,
420008 Kazan, Russia5555institutetext: Department of Physics, IISER Bhopal,
Bhopal Bypass Road, Bhauri, Bhopal 462066, Madhya Pradesh, India
5656institutetext: Ollscoil na Gaillimhe — University of Galway, University
Road, Galway, H91 TK33, Ireland5757institutetext: Center for Gravitation,
Cosmology, and Astrophysics, University of Wisconsin-Milwaukee, Milwaukee, WI
53211, USA 5858institutetext: Division of Natural Science, Faculty of Advanced
Science and Technology, Kumamoto University, 2-39-1 Kurokami, Kumamoto
860-8555, Japan 5959institutetext: International Research Organization for
Advanced Science and Technology, Kumamoto University, 2-39-1 Kurokami,
Kumamoto 860-8555, Japan 6060institutetext: Laboratoire Univers et Théories
LUTh, Observatoire de Paris, Université PSL, CNRS, Université de Paris, 92190
Meudon, France6161institutetext: Max Planck Institute for Gravitational
Physics (Albert Einstein Institute), Leibniz Universität Hannover,
Callinstrasse 38, D-30167, Hannover, Germany6262institutetext: Florida Space
Institute, University of Central Florida, 12354 Research Parkway, Partnership
1 Building, Suite 214, Orlando, 32826-0650, FL, USA6363institutetext: Ruhr
University Bochum, Faculty of Physics and Astronomy, Astronomical Institute
(AIRUB), 44780 Bochum, Germany 6464institutetext: Advanced Institute of
Natural Sciences, Beijing Normal University, Zhuhai 519087, China

I : The European Pulsar Timing Array
II : The Indian Pulsar Timing Array

# The second data release from the European Pulsar Timing Array

IV. Implications for massive black holes, dark matter, and the early
Universe††thanks: Data availability. The EPTA+InPTA DR2 data used to perform
the analysis presented in this paper can be found at:
https://zenodo.org/record/8091568
https://gitlab.in2p3.fr/epta/epta-dr2

J. Antoniadis  0000-0003-4453-776 The second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing Array P. Arumugam  0000-0001-9264-8024 The second data release
from the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing Array S. Arumugam
0009-0001-3587-6622 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array P. Auclair  0000-0002-4814-1406 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array S. Babak  0000-0001-7469-4250 The second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing Array M.
Bagchi  0000-0001-8640-8186 The second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array A.-S. Bak Nielsen  0000-0002-1298-9392 The second data release
from the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array E. Barausse  0000-0001-6499-6263  The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array C. G. Bassa  0000-0002-1429-9010 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array A. Bathula  0000-0001-7947-6703 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array A. Berthereau The
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array M. Bonetti
0000-0001-7889-6810 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array E. Bortolas The second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array P. R. Brook  0000-0003-3053-6538 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array M. Burgay
0000-0002-8265-4344 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array R. N. Caballero  0000-0001-9084-9427 The
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array C. Caprini  0000-0001-5393-2205 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array A. Chalumeau
0000-0003-2111-1001 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array D. J. Champion  0000-0003-1361-7723 The
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array S. Chanlaridis  0000-0002-9323-9728 The second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing Array S. Chen
0000-0002-3118-5963 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array I. Cognard  0000-0002-1775-9692 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array M. Crisostomi
0000-0002-7622-4911 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing Array S. Dandapat
0000-0003-4965-9220 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array D. Deb  0000-0003-4067-5283 The second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing Array S.
Desai  0000-0002-0466-3288 The second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing Array G. Desvignes
0000-0003-3922-4055 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array N. Dhanda-Batra The second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing Array C. Dwivedi
0000-0002-8804-650X The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array M. Falxa The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing Array F. Fastidio The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array R. D. Ferdman The
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array A. Franchini  0000-0002-8400-0969 The second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array J. R. Gair  0000-0002-1671-3668 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array B. Goncharov A. Gopakumar  0000-0003-3189-5807 The second data release
from the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array 0000-0003-4274-4369 The second data release
from the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing Array E. Graikou
The second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array J.-M. Grießmeier  0000-0003-3362-7996 The second data release
from the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array A. Gualandris  0000-0002-9420-2679 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array L. Guillemot  0000-0002-9049-8716 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing Array Y. J. Guo The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array Y. Gupta
0000-0001-5765-0619 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array S. Hisano  0000-0002-7700-3379 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array H. Hu  0000-0002-3407-8071 The second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array F. Iraci The second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing Array D. Izquierdo-Villalba
0000-0002-6143-1491 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array J. Jang  0000-0003-4454-0204 The second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array J. Jawor
0000-0003-3391-0011 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array G. H. Janssen  0000-0003-3068-3677 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array A. Jessner
0000-0001-6152-9504 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array B. C. Joshi  0000-0002-0863-7781 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array F. Kareem
0000-0003-2444-838X The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array R. Karuppusamy  0000-0002-5307-2919 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array E. F. Keane
0000-0002-4553-655X The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array M. J. Keith  0000-0001-5567-5492 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array D. Kharbanda  0000-0001-8863-4152 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array T. Khizriev  The
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing Array T. Kikunaga  0000-0002-5016-3567
The second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array N. Kolhe  0000-0003-3528-9863 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array M. Kramer The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array M. A. Krishnakumar
0000-0003-4528-2745 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array K. Lackeos
0000-0002-6554-3722 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array K. J. Lee The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing Array K. Liu The second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array Y. Liu  0000-0001-9986-9360
The second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array A. G. Lyne The
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array J. W. McKee  0000-0002-2885-8485 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing Array Y. Maan The second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array R. A. Main The second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing Array M.
B. Mickaliger  0000-0001-6798-5682 The second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array H. Middleton
0000-0001-5532-3622 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing Array A. Neronov
The second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array I. C. Niţu  0000-0003-3611-3464 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array K. Nobleson
0000-0003-2715-4504 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array A. K. Paladi  0000-0002-8651-9510 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array A. Parthasarathy  0000-0002-4140-5616 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array B. B. P. Perera
0000-0002-8509-5947 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array D. Perrodin  0000-0002-1806-2483 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array A. Petiteau  0000-0002-7371-9695 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayeeThe
second data release from the European Pulsar Timing ArrayeeThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing Array N. K. Porayko<EMAIL_ADDRESS>second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array A. Possenti The
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array T. Prabu The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array K. Postnov  0000-0002-1705-617X The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array H. Quelquejay Leclere<EMAIL_ADDRESS>second data release
from the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing Array P. Rana
0000-0001-6184-5195 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array A. Roper Pol  0000-0003-4979-4430 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array A. Samajdar  0000-0002-0857-6018 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array S. A. Sanidas The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array D. Semikoz The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array A. Sesana<EMAIL_ADDRESS>second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array G. Shaifullah  0000-0002-8452-4834 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing Array J.
Singha  0000-0002-1636-9414 The second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing Array C. Smarra The second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing Array L. Speri  0000-0002-5442-7267 The second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing Array R.
Spiewak The second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing Array A. Srivastava  0000-0003-3531-7887 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array B. W. Stappers The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array D. A. SteerM. Surnis  0000-0002-8781-1273 The
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing Array 0000-0002-9507-6985 The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array S. C. Susarla  0000-0003-4332-8201 The second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array A. Susobhanan
0000-0002-2820-0931 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array K. Takahashi P. TarafdarG. Theureau
0000-0002-3034-5769 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
Array 0000-0001-6921-4195 The second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing Array 0000-0002-3649-276X The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing Array C. Tiburzi The second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array R. J. Truant The second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing Array E. van der Wateren  0000-0003-0382-8463 The
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing ArrayThe
second data release from the European Pulsar Timing Array S. Valtolina The
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing Array A. Vecchio  0000-0002-6254-1617
The second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array V. Venkatraman Krishnan  0000-0001-9518-9819 The second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing ArrayThe second data release from the European Pulsar Timing Array J.
P. W. Verbiest  0000-0002-4088-896X The second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing Array J. Wang
0000-0003-1933-6498 The second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing Array L. Wang and Z. Wu The
second data release from the European Pulsar Timing ArrayThe second data
release from the European Pulsar Timing ArrayThe second data release from the
European Pulsar Timing ArrayThe second data release from the European Pulsar
Timing Array 0000-0002-1381-7859 The second data release from the European
Pulsar Timing ArrayThe second data release from the European Pulsar Timing
ArrayThe second data release from the European Pulsar Timing ArrayThe second
data release from the European Pulsar Timing ArrayThe second data release from
the European Pulsar Timing ArrayThe second data release from the European
Pulsar Timing Array

(Received XXX; accepted YYY)

The European Pulsar Timing Array (EPTA) and Indian Pulsar Timing Array (InPTA)
collaborations have measured a low-frequency common signal in the combination
of their second and first data releases, respectively, with the correlation
properties of a gravitational wave background (GWB). Such a signal may have
its origin in a number of physical processes including a cosmic population of
inspiralling supermassive black hole binaries (SMBHBs); inflation, phase
transitions, cosmic strings, and tensor mode generation by the non-linear
evolution of scalar perturbations in the early Universe; and oscillations of
the Galactic potential in the presence of ultra-light dark matter (ULDM). At
the current stage of emerging evidence, it is impossible to discriminate among
the different origins. Therefore, for this paper, we consider each process
separately, and investigated the implications of the signal under the
hypothesis that it is generated by that specific process. We find that the
signal is consistent with a cosmic population of inspiralling SMBHBs, and its
relatively high amplitude can be used to place constraints on binary merger
timescales and the SMBH-host galaxy scaling relations. If this origin is
confirmed, this would be the first direct evidence that SMBHBs merge in
nature, adding an important observational piece to the puzzle of structure
formation and galaxy evolution. As for early Universe processes, the
measurement would place tight constraints on the cosmic string tension and on
the level of turbulence developed by first-order phase transitions. Other
processes would require non-standard scenarios, such as a blue-tilted
inflationary spectrum or an excess in the primordial spectrum of scalar
perturbations at large wavenumbers. Finally, a ULDM origin of the detected
signal is disfavoured, which leads to direct constraints on the abundance of
ULDM in our Galaxy.

###### Key Words.:

gravitational waves – black holes: physics – early universe – dark matter –
methods:data analysis – pulsars:general

##  1 Introduction

The recent observation of a common signal with excess power in the nanohertz
frequency ranges (i.e. a common red signal, as defined in Arzoumanian et al.,
2020; Goncharov et al., 2021; Chen et al., 2021) in pulsar timing array (PTA)
datasets, with emerging evidence for quadrupolar correlations111Readers can
refer to EPTA and InPTA Collaborations (2023b) for more details on inter-
pulsar correlations. opens a new era in the exploration of the Universe. This
important milestone has been achieved thanks to the efforts of the European
Pulsar Timing Array (EPTA, Desvignes et al., 2016), the Indian PTA (InPTA
Joshi et al., 2022), the North American Nanohertz Observatory for
Gravitational Waves (NANOGrav, McLaughlin, 2013), the Parkes PTA (PPTA,
Manchester et al., 2013), and the Chinese PTA (CPTA, Lee, 2016). Although the
significance of the signal does not yet reach the 5$\sigma$ mark, which is
generally accepted as the ‘golden rule’ for a firm detection claim, the
evidence reported by the different collaborations ranges between $2\sigma$ and
$4\sigma$ (EPTA and InPTA Collaborations, 2023b; The Nanograv Collaboration et
al., 2023; Reardon et al., 2023; Xu et al., 2023), strongly suggesting a
genuine gravitational wave (GW) origin. Awaiting decisive confirmation within
the International PTA (IPTA, Verbiest et al., 2016; Perera et al., 2019)
framework, with the additional contribution of the MeerKAT PTA (Miles et al.,
2023), we are hearing, for the first time, the faint murmur of the GW Universe
at frequencies of 1-to-30 nano-Hz, which is ten orders of magnitude lower than
the frequencies currently probed by ground-based interferometers (Abbott et
al., 2016). This opens a completely new window onto the Universe, allowing us
to look at different phenomena, probe new astrophysical and cosmological
sources, and, potentially, new physics as well.

By monitoring an array of millisecond pulsars (MSPs) for decades with a weekly
cadence, PTAs are sensitive to GWs in the $10^{-9}$–$10^{-7}$Hz range (Foster
& Backer, 1990). At those frequencies, the most anticipated signal to be
detected is a stochastic GW background (GWB) produced by the incoherent
superposition of waves emitted by adiabatically inspiralling supermassive
black hole binaries (SMBHBs, Rajagopal & Romani, 1995; Jaffe & Backer, 2003).
The signal manifests as a stochastic Gaussian process characterised by a
power-law Fourier spectrum of delays-advances to pulse arrival times, with
characteristic inter-pulsar correlations of general relativity identified by
Hellings & Downs (1983). The statistical properties of the signal are expected
to significantly deviate from the typical isotropy, Gaussianity and perhaps
even stationarity that is typical of many stochastic signals from the early
Universe (e.g. Sesana et al., 2008; Ravi et al., 2012). In fact, due the shape
of the SMBHB mass function and their redshift distribution, the overall signal
is often dominated by a few sources, and particularly massive, nearby SMBHBs
might result in loud enough signals to be individually resolved as continuous
GWs (CGWs, Sesana et al., 2009; Babak & Sesana, 2012; Kelley et al., 2018)
emerging from the GWB. The exact amplitude and spectral shape of the spectrum
are intimately related to the cosmological galaxy merger rate and to the
dynamical properties of the emitting binaries forming in the aftermath of the
merger process (Enoki & Nagashima, 2007; Kocsis & Sesana, 2011; Sesana, 2013a,
b; Ravi et al., 2014). Therefore, the demonstration of an SMBHB origin of the
signal observed by PTAs provides invaluable insights into the formation,
evolution, and dynamics of these objects. Moreover, it brings decisive
evidence that SMBHBs merge in nature, thus overcoming the ‘final parsec
problem’ (Milosavljević & Merritt, 2003), which is still an open question in
our understanding of galaxy and structure formation.

Beyond SMBHBs, a number of processes (potentially) occurring in the early
Universe can also produce a stochastic GWB at nanohertz frequencies. Tensor
modes can be produced as early as during the first tiny fraction of a second
after the Big Bang through quantum fluctuations of the gravitational field
stretched by the accelerated expansion of inflation (Grishchuk, 1975; Rubakov
et al., 1982; Starobinskii, 1985; Fabbri & Pollock, 1983; Abbott & Wise,
1984). In the literature, these GWs are referred to as ‘primordial’. In this
case, the shape of the power spectrum is defined by the specific model of
inflation. Classical tensor mode production invoking the presence of a source
term in the GW equation of motion can also take place in the early Universe.
There are a plethora of physical processes that can lead to such a source
term, and trigger the production of GWs. Topological defects, for example
decaying cosmic string loops (Vilenkin & Shellard, 2000; Damour & Vilenkin,
2001; Jones et al., 2003; Dvali & Vilenkin, 2004; Damour & Vilenkin, 2005),
particle production during inflation (Sorbo, 2011; Barnaby et al., 2012; Cook
& Sorbo, 2013; Anber & Sorbo, 2012), (magneto-)hydrodynamic turbulence ((M)HD,
Kamionkowski et al., 1994; Kosowsky et al., 2002; Dolgov et al., 2002; Caprini
& Durrer, 2006; Gogoberidze et al., 2007; Caprini et al., 2009), the collision
of bubble walls during a first-order primordial phase transition (Kosowsky et
al., 1992; Kosowsky & Turner, 1993; Caprini et al., 2008; Huber & Konstandin,
2008; Jinno & Takimoto, 2017; Cutting et al., 2018), sound waves in the
aftermath of a first-order phase transition (Hindmarsh et al., 2014, 2015,
2017), as well as scalar perturbations at second order in cosmological
perturbation theory (Baumann et al., 2007; Ananda et al., 2007), are among the
most commonly considered scenarios. GWs decouple from the rest of the matter
and radiation immediately after their generation at essentially any
temperature in the Universe, so that in the case of clear observational
evidence for these types of signals, we can infer nearly unaltered information
on the physical processes occurring during or just after the birth of the
Universe (Caprini & Figueroa, 2018).

Contrary to the incoherent superposition of GWs from SMBHBs, the stochastic
GWB from sources in the early Universe is usually assumed to be statistically
homogeneous and isotropic, unpolarised, and Gaussian (Allen, 1996; Maggiore,
2000; Caprini & Figueroa, 2018). Statistical homogeneity and isotropy are
inherited from the same property of the FLRW Universe. The absence of
polarisation holds provided no macroscopic source of parity violation is
present in the Universe. Gaussianity follows by the central limit theorem in
most cases of GWBs generated by processes operating independently in many
uncorrelated, sub-horizon regions. This also applies to the irreducible GWB
generated during inflation in the simplest scenarios, when the tensor metric
perturbation can be quantised as a free field, and hence with Gaussian
probability distribution for the amplitude. There are, however, notable
exceptions, among which for example rare GWB bursts from cosmic strings cusps
and kinks (Damour & Vilenkin, 2000, 2001), or the GWB from particle production
during inflation (Cook & Sorbo, 2013; Sorbo, 2011; Anber & Sorbo, 2012).
Therefore, although statistical properties might be useful for discriminating
between SMBHBs and several processes in the early Universe, a full assessment
of the nature of the GW signal will not be trivial.

Spatially correlated delays of the time of arrivals (TOAs) in an array of MSPs
are not a unique imprint of GWs. For example, it is well known (e.g. Tiburzi
et al., 2016) that such delays can emerge due to the imperfect fitting of the
solar system ephemerides (dipolar correlated noise), or due to a
miscalibration of the time standard to which the measured TOAs are referred
(monopolar correlated noise). Furthermore, individual Fourier harmonics of a
common signal in PTA data may include contributions from the oscillations of
the gravitational potential associated with the presence of ultralight dark
matter (ULDM, Smarra et al., 2023)222In the pulsar timing community, it is
customary to refer to the dispersion measure with the acronym DM. Please note
that in this paper, DM will exclusively stand for Dark Matter., also known as
fuzzy dark matter (FDM), in the Galactic halo (Khmelnitsky & Rubakov, 2014).
The existence of ultralight scalars, generally motivated by string-theoretical
frameworks (Green et al., 1988; Svrcek & Witten, 2006; Arvanitaki et al.,
2010), is also particularly appealing from the astrophysical and cosmological
point of view. In fact, several potential issues in the small-scale structure
of the Universe, such as the cusp vs core (Flores & Primack, 1994; Moore,
1994; Karukes, E. V. et al., 2015) or the missing satellite (Klypin et al.,
1999; Moore et al., 1999) problems, could be disposed of or, at least,
mitigated assuming that dark matter is made of ultralight particles. As
predicated by Khmelnitsky & Rubakov (2014), the presence of ULDM induces
harmonic delays in the arrival times, with a frequency proportional to the
ultralight boson mass.

In this paper, we provide a broad overview of the implications of the signal
observed in the second data release of the EPTA+InPTA (EPTA and InPTA
Collaboration, 2023) for the different physical processes mentioned above.
More in-depth analysis of several of these scenarios will be the subject of
separate future publications. Unless otherwise stated, we consider each
process separately, and we discuss the implications of the detected signal
under the hypothesis that it is generated by that specific process. We do not
attempt any Bayesian model selection on the signal origin, although a general
framework for that is being developed (Moore & Vecchio, 2021). The main reason
for this choice is that, at this stage of data taking and analysis, the
information carried by the signal is not particularly constraining; the
evidence of the measurement is still at $\approx 3\sigma$, and the amplitude
and spectral shape of the signal are not very well measured. With these
premises, the result of any model selection is bound to be severely influenced
by the priors employed for each of the models under examination. This exercise
becomes more meaningful as data get more informative, which we strive to
achieve with the analysis of the third release of the combined IPTA data,
which is now being assembled.

The paper is organised as follows. In Sec. 2, we describe the signal observed
by EPTA+InPTA and its main features, including its free spectrum and best-fit
parameters. We then proceed with detailing the main implications of the
detected signals under the assumption that it is generated by a cosmic
population of SMBHBs (Sec. 3) or by a number of processes occurring in the
early Universe (Sec. 4). In Sec. 5, we investigate the compatibility of the
observed signal with a DM origin and place constraints on ULDM candidates.
Finally, in Sec. 6, we summarise our main results and discuss future
prospects.

##  2 The observed signal in the EPTA DR2 dataset

Figure 1: Properties of the common correlated signal detected in DR2new. Left
panel: free spectrum of the RMS induced by the excess correlated signal in
each frequency resolution bin (with width defined by the inverse of the data
span, $\Delta{f}=T^{-1}$). The straight line is the best power-law fit to the
data. Right panel: joint posterior distribution in the $A-\gamma$ plane. Note
that we normalise $A$ to a pivotal frequency $f_{0}=10{\rm yr}^{-1}$.

Our investigation is based on the results reported in EPTA and InPTA
Collaborations (2023b, hereinafter PaperIII), which analyses the data of 25
MSPs collected by the EPTA using five of the largest radio telescopes in
Europe: the Lovell telescope at the Jodrell Bank Observatory, the Nançay
decimetric radio telescope, the Westerbork synthesis radio telescope, the
Effelsberg 100 m radio telescope, and the Sardinia radio telescope. The
dataset also includes the Large European Array for Pulsars (LEAP) data, in
which individual telescope observations are coherently phased to form an
equivalent dish with a diameter of up to 194 m (Bassa et al., 2016). These
data are complemented by low-frequency observations of a subset of 10 MSPs
performed by the InPTA using the upgraded Giant Metrewave Radio Telescope
(uGMRT) and covering about 3.5 years.

The data of each individual pulsar are combined as described in EPTA and InPTA
Collaboration (2023) and the noise properties of each pulsar are then
extracted according to the optimised custom noise models presented in EPTA and
InPTA Collaborations (2023a). The final result is a dataset of unprecedented
sensitivity spanning up to 24.7 years. Four versions of the dataset were
analyzed:

  1. 1.

DR2full. 24.7 years of data taken by the EPTA;

  2. 2.

DR2new. 10.3 years of data collected by the EPTA using new-generation wide-
band backends;

  3. 3.

DR2full+. The same as DR2full, but with the addition of InPTA data;

  4. 4.

DR2new+. The same as DR2new, but with the addition of InPTA data.

The analysis presented in this paper refers to the DR2new dataset only. We do
not consider DR2full and DR2full+ because evidence of quadrupolar correlation
(usually referred to as HD correlation, from Hellings & Downs, 1983) of the
common process is weaker in those datasets, potentially due to the lower
quality of early data that were collected with narrowband backends (see
discussion in PaperIII). On the other hand, although the analysis of DR2new+
produced results in broad agreement with DR2new, that dataset was assembled
relatively recently and has not been analysed as thoroughly. For example, the
binned free-spectra that we will use in some of the following analyses have
only been produced after this work was completed.

Before proceeding with the description of the signal detected in DR2new, here
we summarise some notations used in PTA analysis for the benefit of the
reader. The perturbation affecting the TOAs, whether produced by GWs or DM, is
described in terms of its dimensionless strain $h$. A broad-band stochastic
perturbation is defined by its characteristic dimensionless strain $h_{c}(f)$,
often modelled as a power law

$h_{c}(f)=A\left(\frac{f}{f_{0}}\right)^{\alpha}.$ (1)

For example, a population of GW-driven circular SMBHBs produces a spectrum
with $\alpha=-2/3$ and amplitude $A\approx 10^{-15}$, assuming $f_{0}=1{\rm
yr}^{-1}$. $h_{c}(f)$ is connected to the differential energy content of the
signal per logarithmic frequency through the equation:

$\Omega(f)=\frac{2\pi^{2}}{3H_{0}^{2}}f^{2}h_{c}^{2}(f),$ (2)

where $H_{0}$ is today’s Hubble expansion parameter. We note that $h_{c}(f)$
and $\Omega(f)$ provide equivalent parametrizations of the spectrum. The
former is more popular in the astrophysics domain, whereas the latter is the
preferred choice for early Universe and cosmology.

Given $h_{c}(f)$, the one-sided power spectral density induced by the GW
signal in the timing residuals is given by (Lentati et al., 2015):

$S(f)=\frac{h_{c}^{2}(f)}{12\pi^{2}f^{3}}=\frac{A^{2}}{12\pi^{2}f_{0}^{2\alpha}}f^{-\gamma},$
(3)

where $\gamma=3-2\alpha$. PTAs search for HD correlated time delays with such
a power spectrum in the data, and measure the parameters $A$ and $\gamma$. For
an observation timespan $T$, measurements are discretised in frequency bins
$\Delta{f_{i}}=f_{i+1}-f_{i}$, where $f_{i}=i/T$. It is then customary to
convert $S(f)$ in RMS residual induced in the TOAs in each frequency bin:

${\rm RMS}_{i}=\left(\int_{\Delta{f_{i}}}S(f){\rm
d}f\right)^{1/2}\approx\left(S(f_{i})\Delta{f_{i}}\right)^{1/2}=\left(\frac{S(f_{i})}{T}\right)^{1/2}.$
(4)

The main properties of the GWB signal observed in DR2new and examined in this
paper are shown in Fig. 1. The length of the dataset is $T=$10.3 years, and
excess common correlated power is detected in several frequency bins up to
$\approx 30\,$nHz (Fig. 1 left panel). Conversely, some bins are
unconstrained, which results in a relatively loose determination of the
spectral properties of the observed signal. In the literature, $h_{c}$ and $S$
in Eqs. (1) and (3) are usually anchored to the pivotal frequency $f_{0}=1{\rm
yr}^{-1}$. The data are, however, most informative at the lowest frequencies,
while the common power at $1{\rm yr}^{-1}$ is essentially unconstrained. This
naturally leads to a strong degeneracy of the $A-\gamma$ 2D posterior, as
shown for example in Figure 1 of PaperIII. Therefore, unless otherwise stated,
we change the reference frequency to $f_{0}=10{\rm yr}^{-1}$, where the data
are actually constraining, which results in a weaker dependence of $A$ upon
$\gamma$, as shown in the right panel of Fig. 1.

In the following three sections, we discuss three possible contributions to
the signal, probing completely different epochs and scales of our Universe,
and the implications for the associated physical processes. Namely, the cosmic
population of SMBHBs (at redshifts $z\lesssim 1$), the early Universe
($z>1000$), and DM (within our Galaxy).

##  3 Implications I: supermassive black hole binaries

A cosmic population of SMBHBs is the primary astrophysical candidate to
produce a signal in the nanohertz band detectable by PTAs. If we define ${\rm
d}^{5}N/({\rm d}z{\rm d}m_{1}{\rm d}q{\rm d}e{\rm d}t_{r})$ as the cosmic
merger rate of SMBHBs as a function of redshift, primary black hole mass, mass
ratio, and eccentricity, the general form of the generated GWB as a function
of observed frequency $f$ can be written as (Sesana, 2013a)

$h_{c}^{2}(f)=\int_{0}^{\infty}{\rm d}z\int_{0}^{\infty}{\rm
d}m_{1}\int_{0}^{1}{\rm d}q\frac{{\rm d}^{5}N}{{\rm d}z{\rm d}m_{1}{\rm
d}q{\rm d}e{\rm d}t_{r}}\frac{{\rm d}t_{r}}{{\rm d}{\rm ln}f_{{\rm
K},r}}\times\\\ h^{2}(f_{{\rm K},r})\sum_{n=1}^{\infty}\frac{g[n,e(f_{{\rm
K},r})]}{(n/2)^{2}}\bigg{|}_{f_{K,r}=f(1+z)/n}.$ (5)

Here, ${\rm d}t_{r}/{{\rm d}{\rm ln}f_{{\rm K},r}}$ is the time spent by the
shrinking binary within a given logarithmic orbital frequency bin, which
converts the merger rate into the distribution of rest-frame orbital
frequencies of the emitting population. This quantity depends on the physical
processes driving the evolution of the SMBHBs including, besides GW emission,
interaction with the stellar and gaseous environment surrounding them. As
such, it is generally a function of the binary parameters $m_{1},q,e$, and
extra parameters describing the environment, such as the stellar density in
the nucleus of the galaxy host (for more details, see Sesana, 2013a). The
second line of Eq. (5) is the strain amplitude produced by each individual
eccentric SMBHB binary, cast as the sum of harmonics fulfilling the condition
${f_{K,r}=f(1+z)/n}$. In that expression, $h(f_{K,r})$ is the strain of the
equivalent circular binary emitted at twice the orbital frequency of the
system, as given in Eq. 4-7 of Rosado et al. (2015), and $g(n,e)$ is a
combination of Bessel functions (see, e.g. Bonetti & Sesana, 2020, for
details). For a distribution of circular, GW-driven binaries, the only
relevant mass parameter is the chirp mass
$\mathcal{M}=(m_{1}m_{2})^{3/5}(m_{1}+m_{2})^{1/5}$, and Eq. (5) takes the
familiar form (Sesana et al., 2008)

$h_{c}^{2}(f)=\int_{0}^{\infty}{\rm d}z\int_{0}^{\infty}{\rm
d}\mathcal{M}\frac{{\rm d}^{3}N}{{\rm d}z{\rm d}{\cal M}{\rm d}{\rm
ln}f_{r}}h^{2}(f_{r}).$ (6)

This can be recast in terms of the comoving number density of merging binaries
${{\rm d}^{2}n}/({{\rm d}z{\rm d}\mathcal{M}})$ (Phinney, 2001)

$h_{\rm c}^{2}(f)=\frac{4G^{5/3}}{3\pi^{1/3}c^{2}}f^{-4/3}\int{\rm
d}\mathcal{M}\int{\rm d}z\left(1+z\right)^{-1/3}\mathcal{M}^{5/3}\frac{{\rm
d}^{2}n}{{\rm d}z{\rm d}\mathcal{M}}\,,$ (7)

which highlights that, in this case, the expected spectrum follows a power law
$h_{c}\propto f^{-2/3}$, and the only free parameter is its overall amplitude.
The latter is set by the function ${{\rm d}^{2}n}/({{\rm d}z{\rm
d}\mathcal{M}})$, which contains all the astrophysical knowledge of the cosmic
population of merging SMBHBs, and is determined by the cosmological
hierarchical assembly of galaxies and their central SMBHs. Conversely, in the
most general case described by Eq. (5) there is also information in the
spectral shape of the signal, since coupling with the environment as well as
eccentricity affect the function ${\rm d}t_{r}/{{\rm d}{\rm ln}f_{{\rm
K},r}}$, suppressing signal at the lowest frequencies. Moreover, eccentricity
distributes the GW power at several higher harmonics of the orbital frequency,
slightly modifying the power-law behaviour at high frequencies. In general,
the GWB cannot be cast in term of a simple analytical form, although a broken
power-law is a sufficient approximation for most situations (see, e.g. Ravi et
al., 2014; Sesana, 2015; Kelley et al., 2017; Chen et al., 2017b).

The literature investigating the GWB produced by a population of SMBHBs is
vast, dating back to the mid-nineties and early 2000s (Rajagopal & Romani,
1995; Jaffe & Backer, 2003; Wyithe & Loeb, 2003; Sesana et al., 2004), and
predictions have been made by employing different models and techniques.
Models can be broadly classified into two categories: self-consistent
theoretical models for SMBH evolution within their galaxies (Sesana et al.,
2008, 2009; Ravi et al., 2012; Kulier et al., 2015; Kelley et al., 2017;
Bonetti et al., 2018a; Siwek et al., 2020; Izquierdo-Villalba et al., 2022),
and empirical models based on observed properties of galaxy pairs coupled to
SMBH-host galaxy relations (Sesana, 2013b; Rosado et al., 2015; Ravi et al.,
2015; Simon, 2023), or on the evolution of the SMBH mass function inferred
from observations (McWilliams et al., 2014). Note that we group both
semianalytic models (SAMs) and large cosmological simulations in the first
class. The main difference between these two classes is that self-consistent
models are constructed to reproduce a large array of observations related to
galaxies and the SMBH they host, such as the redshift-dependent galaxy mass
function, quasar luminosity function, and so on. Conversely, empirical models
are, by construction, consistent with the observations upon which they are
based, but are usually not tested against independent constraints. As a
consequence, they can generally produce a wider distribution of GWB
amplitudes, but consistency with other observations is not necessarily
guaranteed.

In this Section, we investigate the implications of the signal observed in the
DR2new dataset for representative models of the two classes. In Sec. 3.1, we
perform a semi-quantitative comparison between the measured signal and
predictions of an extended version of the Rosado et al. (2015) models
(hereinafter RSG15) including binary eccentricity and environmental coupling.
In Sec. 3.2, we exploit the framework developed in Middleton et al. (2016);
Chen et al. (2017a, 2019) to draw inference on SMBHB astrophysics from the
data, either by assuming astrophysical priors from independent observations,
or by using a completely generic model for the SMBHB mass function with
minimal assumptions. In Sec. 3.3, we demonstrate how the measured signal can
inform galaxy and SMBH formation models by examining its constraining power on
two state-of-the-art SAMs, namely L-galaxies (Henriques et al., 2015) and the
model developed by Barausse and collaborators (Barausse, 2012; Bonetti et al.,
2018b; Barausse et al., 2020). We discuss caveats and future directions of
improvement in Sec. 3.4.

###  3.1 Qualitative analysis of empirical SMBHB population models

To carry out a semi-quantitative comparison between observations and empirical
models, we use an extended set of SMBHB populations based on the work of
Sesana (2013b) (S13 hereinafter) and RSG15.

####  3.1.1 Description of the models

As described in S13, the models are constructed around the number density of
merging SMBHs per unit primary mass, mass ratio, and redshift, ${\rm
d}^{3}n/({\rm d}m_{1}{\rm d}q{\rm d}z)$ obtained by combining different
observations of the galaxy mass function and pair fraction, estimated galaxy
pair merger timescales, SMBH-host galaxy relations, and prescription for SMBH
accretion during mergers (see Section 2 of S13 for full details). Given the
large uncertainties in all of the ingredients, the models predict a broad
distribution of expected GW amplitudes, as shown in Fig. 2.

Figure 2: GWB amplitude distributions predicted by the RSG15 models. The thin-
dashed yellow line is for the full set of models in RSG15, whereas the thick-
dashed orange line is for the subset considered here. The solid blue line is
the distribution predicted by the 108 down-selected sample used in the
analysis. The vertical line marks the median value of $A$ at $f_{0}=1{\rm
yr}^{-1}$ reported in PaperIII when fixing $\gamma=13/3$ in the search. Note
that the lower $x$-axis scale is for A at $f_{0}=1{\rm yr}^{-1}$, whereas the
upper $x$-axis is for A at $f_{0}=10{\rm yr}^{-1}$ (the normalization used in
this paper). Since $\alpha=-2/3$ for circular GW-driven binaries, there is a
shift of 0.666 dex between the two.

Guided by the relatively large amplitude of the detected signal and by
theoretical and observational advancements in the last decade, we select a
sub-sample of those models, as we now justify. First, hydrodynamical
simulations of merging galaxies at different scales as well as deep X-ray
observations of merging galaxies support accretion activation onto the
individual SMBHs prior to merger (e.g. Capelo et al., 2017; Koss et al.,
2018). Moreover, hydrodynamical simulations of sub-pc scale binaries, have
found most of the accretion to occur on the secondary (i.e. less massive) SMBH
(Farris et al., 2014). We therefore restrict the analysis to models where
SMBHs accrete prior to the merger, either with an equal Eddington ratio or
with preferential accretion onto the secondary.333This is implemented
according to a simple scheme whereby accretion is quenched on the primary SMBH
until the binary becomes equal mass. If/when this occurs, further accretion is
equally distributed among the two (now equal mass) binary components. Second,
observations of overmassive black holes in the centre of large ellipticals
(McConnell et al., 2011) has led to an upward revision of the SMBH-galaxy
relations. Contrary to S13 and RSG15, here we consider only those revised
realations, namely the ones reported by Kormendy & Ho (2013); McConnell & Ma
(2013); Graham & Scott (2013). Finally, given the large number of models, to
save computing power, we perform an ad hoc down-selection of 108 models that
preserves the overall distribution of the expected GWB amplitudes, as shown in
Fig. 2.

As opposed to RSG15, we go beyond the circular-GW driven binary approximation
and consider the self-consistent evolution of SMBHBs within their stellar
environment.444Most of the signal comes from SMBHBs hosted in massive, low-
redshift galaxies (but see Simon, 2023), which are relatively gas-poor. This
is done by employing the hardening models of Sesana (2010) that self-
consistently evolve the SMBHB semimajor axis and eccentricity under the
combined effect of stellar hardening and GW emission, once a given initial
eccentricity $e_{0}$ at binary formation is given. Those evolutionary tracks
allow us to evaluate the term ${\rm d}t/{{\rm d}{\rm ln}f_{{\rm K},r}}$ in Eq.
(5), and thus to reconstruct from the density distribution of merging
binaries, ${\rm d}^{3}n/({\rm d}m_{1}{\rm d}q{\rm d}z)$, the numerical
distribution of systems emitting at any time in the whole sky as a function of
mass, mass ratio, redshift, orbital frequency and eccentricity, ${\rm
d}^{5}N/({\rm d}m_{1}{\rm d}q{\rm d}z{\rm d}f{\rm d}e)$. For each model, we
consider 10 values of $e_{0}=0,0.1,...,0.9$ and three different normalizations
of the stellar density profile, described as
$\rho=C\times\rho_{0}(r/r_{0})^{-1.5}$, with $C=0.1,1,10$ (details in Sesana,
2010).

In total, we explore 108$\times$10$\times$3$=3240$ models, spanning different
eccentricities and densities of the stellar environment. For each model, we
draw 100 Monte Carlo samplings of the distribution of the emitting binaries,
we discretise the frequency domain in bins of $\Delta{f}=10.3{\rm yr}^{-1}$,
and add the signals of binaries falling in the same bin in quadrature. This
leads to the binned characteristic strain spectrum $h_{c}(f)$ that we then
convert in $S(f)$ and RMS residuals using Eqs. (3) and (4). The full procedure
is described in Amaro-Seoane et al. (2010); Bonetti & Sesana (2020). In this
way, we generate a grand total of 324k Monte Carlo realizations of the
predicted GW spectrum in the PTA band.

####  3.1.2 Comparison with the observed signal

Figure 3: Free spectrum violin plot comparing measured (orange) and expected
(green) signals. Overlaid to the violins are the 100 Monte Carlo realizations
of one specific model; among those, the thick one represents an example of a
SMBHB signal consistent with the excess power measured in the data at all
frequencies.

The binned spectrum shown in Fig. 3 contrasts expectations from the 324k
models (green) to the measured correlated signal in DR2new (orange). The two
sets of violin plots are in good agreement in the few lowest frequency bins,
where measurements are the most constraining. Note that the model prediction
distributions are highly non-Gaussian and asymmetric, with long tails
extending upwards. This is due to the fact that sparse very massive/nearby
binaries can sometimes produce exceptionally loud signals, as illustrated by
the 100 individual GWBs overplotted to the violins. In fact, this might
explain the extra power measured in the 4th and, most strikingly, in the 9th
lowest bins compared to the bulk of the model predictions. We caution that the
9th bin is close to the 1yr-1 mark, where PTAs are blind due to fitting for
the Earth orbital motion, and leakage from imperfect fitting might affect that
measurement. In any case, if this extra power is indeed due to GWs, it can be
easily accommodated by theoretical models, as demonstrated by the realization
highlighted by the tick grey line.

Figure 4: Expected properties of CGWs as a function of frequency. Top panel:
free spectrum violin plot comparing the measured signal (orange) to the power
distribution of CGWs (green). Empty violins show the full GWB produced by the
models for comparison. Bottom panel: the probability of detecting a CGW with
S/N$>3$ as a function of frequency (green circles, left $y-$axis scale). The
average S/N of CGWs is also shown as red crosses (right $y-$axis scale).

Our Monte Carlo approach to generate the SMBHB population and its associated
GW signal also allows us to investigate the occurrence of CGWs in the data,
for which evidence in DR2new is found to be inconclusive (EPTA and InPTA
Collaborations, 2023c). Since the search performed in that paper was limited
to circular binaries, we only carry out this analysis for the 32.4k models
with $e_{0}=0$.555For binaries starting with $e_{0}=0$, eccentricities remain
well below $0.1$ in the course of stellar hardening-driven evolution, and the
GW signals can be approximated as monochromatic. A full assessment of the
detectability of CGWs requires the evaluation of the detection probability of
each individual binary for a given false alarm rate, as detailed in RSG15. For
the sake of simplicity, and given the qualitative nature of this analysis, we
just compute the signal-to-noise ratio (S/N) of each individual binary
according to Eq. (46) of RSG15 (thus also restricting to the Earth term only).
When computing the S/N of a source, we model each pulsar noise by using the
maximum likelihood values of the single pulsar noise analysis presented in
EPTA and InPTA Collaborations (2023a), and add the GWB produced by all of the
other binaries to the noise spectral density. We arbitrarily set the
detectability threshold at S/N$=3$ in the following.

Results are shown in Fig. 4, which compares the power distribution of
resolvable CGWs to the binned spectra of the overall predicted GW signal and
of the DR2new measurements. In line with RSG15, the probability of detecting a
CGW is maximum at the lowest frequency, rapidly decaying to less than 0.01
past the 6th bin. Although this seems to suggests that the feature observed at
the 9th frequency bin is unlikely to be due to a CGW, it should be noticed
that we considered here a threshold of $S/N=3$. The overall GW signal in our
data has a total $S/N\approx 3.5-4$ (EPTA and InPTA Collaborations, 2023b),
mostly accumulated at the lowest frequency bins. It might still be the case
that the feature at the 9th bin is due to an unresolved CGW with S/N$<3$,
which would be a more common occurrence in the data. Note that the average S/N
of CGWs slightly increases at higher frequencies, which is primarily due to
the frequency dependence of the CGW characteristic strain, $h_{c}\propto
f^{7/6}$. If $p_{i}$ is the probability of having a CGW of S/N $>3$ in the
$i$-th bin, we can compute the probability of detecting at least one CGW with
S/N $>3$ in DR2new according to these models as $p=1-\prod(1-p_{i})$, which
gives $p=0.49$. It is therefore quite possible that the observed signal is
dominated by few bright sources, which might be resolvable in the future with
more sensitive datasets. Thus far, searches for CGWs on the current dataset
yielded inconclusive evidence (EPTA and InPTA Collaborations, 2023c). This
probability is obviously S/N threshold dependent. For example, by increasing
this threshold to S/N$>5$, we get $p=0.13$. This is comparable to the 6%
chance found by Bécsy et al. (2022). and the slightly larger probability in
our models is likely due to the louder overall amplitudes of the signals
considered here. We stress, however, that these findings apply to models where
binaries remain essentially circular. The number of resolvable CGWs, tends, in
fact, to slightly decrease when the eccentricity increases (Truant et al. in
preparation).

Figure 5: $A-\gamma$ distribution of the measured signal (orange) compared to
model predictions (green). 1$\sigma$ and 2$\sigma$ contours are displayed.
Shown are also the marginalised $A$ (left) and $\gamma$ (right) distributions,
with their 1$\sigma$ credible intervals highlighted as shaded areas.

Finally, we once again propose the comparison first shown by Middleton et al.
(2021), who contrasted the measured 2D $A-\gamma$ posterior to model
expectations. For the latter, we just fit the 9 lowest frequency bins of the
GWB spectrum of each Monte Carlo realization of the Universe with a straight
line in the ${\rm log}A-{\rm log}f$ plane. As described in Sec. 2, we
normalise the measurement to $f_{0}=10$yr-1, where our data are informative,
which alleviates the $A-\gamma$ degeneracy in the posterior. Results are shown
in Fig. 5. Although the measured spectrum tends to be shallower than the
theoretical one (see also Sec. 3.4), the contours overlap at 2$\sigma$ and the
marginalised amplitudes are broadly consistent.

###  3.2 Inference on the SMBHB population.

After checking the general compatibility of the observed signal with
expectations from empirical SMBHB population models, we take a more
quantitative approach and exploit Bayesian inference to constrain the
properties of SMBHBs from the data. We repeat the analysis carried out by
Middleton et al. (2021) on the common signal detected in the NANOGrav 12.5
year dataset (Arzoumanian et al., 2020), exploiting the framework laid out in
Middleton et al. (2016); Chen et al. (2017a, 2019). The SMBHB population of a
given model M is described by a set of parameters $\\{\theta\\}$, we then use
Bayesian inference to find the posterior distribution $p(\theta|d,{\bf M})$
for the parameters $\theta$ of model ${\bf M}$ given the observation data $d$:

$p(\theta|d,{\bf M})=\frac{p(\theta|{\bf M})p(d|\theta,{\bf M})}{p(d|{\bf
M})}.$ (8)

Here, $p(\theta|{\bf M})$ is the prior distribution on the model parameters,
$p(d|\theta,{\bf M})$ is the likelihood of model ${\bf M}$ with parameters
$\theta$ of producing the data, and $p(d|{\bf M})$ is the evidence. For each
set of $\theta$, the likelihood is computed by comparing the signal amplitude
produced by ${\bf M}$ at frequencies $f_{i}=i/(10.3{\rm yr}),i=1,...,9$, to
the posterior distribution of the correlated signal measured in DR2new at the
same frequencies. We treat each bin as independent, therefore the likelihood
takes the form

${\rm log}_{10}p(d|\theta,{\bf M})=\sum_{i=1}^{9}{\rm
log}_{10}p(A,f_{i})|_{A=A_{\bf M}}$ (9)

where $p(A,f_{i})$ is the probability density of the amplitude $A$ of the
correlated signal measured in the $i$-th bin, and it is evaluated at the value
$A_{\bf M}$ predicted by the model. We estimate the likelihood in each bin
using a KDE estimate of the DR2new posteriors, similar to the method used in
Moore & Vecchio (2021).

Note that the models we use are deterministic in the sense that they have a
1:1 correspondence between $\theta$ and the predicted $h_{c}(f)$. In reality,
given $\theta$, the ensemble of SMBHBs generating the signal is statistically
drawn from the deterministic distribution function, which results in a
significant scatter of $h_{c}(f)$, as demonstrated by the individual spectra
shown in Fig. 3. We caution that this variance is not captured by these
models, and its inclusion in the Bayesian inference pipeline is the subject of
ongoing work.

As in Middleton et al. (2021), we use two models to describe the SMBHB
population, which are described in turn below.

####  3.2.1 Agnostic SMBHB population model

The agnostic model, developed in Middleton et al. (2016), makes minimal
assumptions about the underlying population of SMBHBs. Binaries are assumed to
be circular, GW-driven and the characteristic strain is computed according to
Eq. (7) where the source distribution is described by a Schechter function
(Schechter, 1976) in $z$ and $\mathcal{M}$,

$\frac{{\rm d^{2}}n}{{\rm d}z{\rm d}\log_{10}\mathcal{M}}={\dot{n}}_{\rm
0}\left[\left(\frac{\mathcal{M}}{10^{7}{\rm
M_{\odot}}}\right)^{-\alpha_{\rm\mathcal{M}}}e^{-\mathcal{M}/\mathcal{M}_{\rm\star}}\right]\left[\left(1+z\right)^{\beta_{\rm
z}}e^{-z/z_{\rm 0}}\right]\frac{{\rm d}t_{\rm R}}{{\rm d}z}\,,$ (10)

where $t_{R}$ is the time in the source frame and we assume cosmological
parameters from Planck18 (Planck Collaboration et al., 2020a). The five model
parameters are $\theta=\\{{\dot{n}}_{\rm
0},\alpha_{\rm\mathcal{M}},\mathcal{M}_{\rm\star},\beta_{\rm z},z_{\rm 0}\\}$,
where ${\dot{n}}_{\rm 0}$ is the merger rate per unit rest-frame time, co-
moving volume, and logarithmic $\mathcal{M}$ interval, and the parameter pairs
$\\{\alpha_{M},\mathcal{M}_{\rm\star}\\}$ and $\\{\beta_{\rm z},z_{\rm 0}\\}$
control the shape of the $\mathcal{M}$ and $z$ distributions, respectively.
The integration limits in $\mathcal{M}$ and $z$ are
$10^{6}\leq\mathcal{M}/{\rm M_{\odot}}\leq 10^{11}$ and $0\leq z\leq 5$,
respectively. The prior ranges of the five model parameters are identical to
those used in (Middleton et al., 2021).

####  3.2.2 Astrophysically-informed SMBHB population model

The astrophysically-informed model was developed in Chen et al. (2019). This
model captures the interaction between the SMBHBs and their environment and
allows for eccentric orbits, both of which lead to a characteristic amplitude
that does not follow a simple single power law, as in Eq. (5). The model has
$18$ parameters, $16$ of which describe astrophysical observables linking the
number of SMBHB mergers to the number of galaxy mergers. The galaxy stellar
mass function is modelled as a redshift-dependent Schechter function defined
by five parameters: $\\{\Phi_{\rm 0},\Phi_{I},M_{0},\alpha_{\rm 0},\alpha_{\rm
I}\\}$. Both the galaxy pair fraction and merger timescales have power law
dependencies on the primary galaxy stellar mass $M$, mass ratio $q$ and
redshift function $(1+z)$ and each of them is defined by a set of four
parameters: $\\{f_{\rm 0},\alpha_{\rm f},\beta_{\rm f},\gamma_{\rm f}\\}$ for
the pair fraction, and
$\\{\tau_{0},\alpha_{\rm\tau},\beta_{\rm\tau},\gamma_{\rm\tau}\\}$ for the
merger timescale:

$f_{\rm pair}=f_{0}\left(\frac{M}{10^{11}{\rm
M_{\odot}}}\right)^{\alpha_{f}}(1+z)^{\beta_{f}},\,\,\,{\rm
with}\,\,\,\,\,f_{0}=f^{\prime}_{0}\int q^{\gamma_{f}}{\rm d}q;$ (11)
$\tau=\tau_{0}\left(\frac{M}{10^{11}{\rm
M_{\odot}}}\right)^{\alpha_{\tau}}(1+z)^{\beta_{\tau}}q^{\gamma_{\tau}}.$ (12)

Galaxy pairs are then populated with SMBHs of mass $m$ following a standard
black hole to stellar bulge mass relation of the form

$m={\cal N}\Biggl{\\{}M_{*}\left(\frac{M_{b}}{10^{11}{\rm
M_{\odot}}}\right)^{\alpha_{*}},\epsilon\Biggl{\\}}$ (13)

where ${\cal N}\\{x,y\\}$ is a log normal distribution with mean value $x$ and
standard deviation $y$, which adds three further parameters,
$\\{M_{*},\alpha_{*},\epsilon\\}$, to the model. The final two parameters
describe the eccentricity at SMBHB pairing $\\{e_{\rm 0}\\}$ and the density
of the stellar environment $\\{\zeta_{\rm 0}\\}$. For the 18 parameters listed
above, in the analysis presented here, we use the extended prior intervals
listed in Table I of Chen et al. (2019).

####  3.2.3 Results of the inference

Figure 6:  Marginalised posterior distributions for ${\dot{n}}_{\rm 0}$ using
two SMBHB population models. The orange and green open histograms show
marginalised posteriors for the agnostic and astrophysically-informed models,
respectively. The filled-green histogram shows the prior for the
astrophysically-informed model (the prior for the agnostic model is uniform in
the range $[-20,3]$). The vertical dotted lines show the $5{\rm th}$ and
$95{\rm th}$ percentiles of the posteriors.  Figure 7: Posterior distribution
of the chirp mass function of merging SMBHBs for both the agnostic (orange)
and astrophysically informed (green) models. For both models, shaded areas are
the central 50% and 90% credible regions and the dashed lines show the medias.
The black-dotted lines show the central 99% region for the astrophysical
prior. Figure 8:  Posterior distribution of selected parameters for the
astrophysically-informed model using nine frequency bins of the free spectrum
for the inference. Parameters are described in Sec. 3.2.2.

The main results of the inference are shown in Figs. 6, 7 and 8. Fig. 6 shows
the marginalised posterior distribution for the normalization of the merger
rate density ${\dot{n}}_{\rm 0}$ from the agnostic model compared to an
equivalent parameter derived from the astrophysically-informed model. The
constraint on the amplitude of the signal imposes an informative constraint on
the normalization of the SMBHB merger density. Using nine frequency bins, the
median and central $90\%$ credible regions for $\log_{10}{\dot{n}}_{\rm
0}/[{\rm Mpc}^{3}{\rm Gyr}]$ are $-1.95^{+2.91}_{-5.87}$ and
$-3.51^{+0.59}_{-0.62}$ for the agnostic and astrophysically-informed models,
respectively. The measurement essentially constrains the amplitude of the
signal, which imposes an informative constraint on ${\dot{n}}_{\rm 0}$. The
astrophysically-informed model clearly shows that the signal favours an
${\dot{n}}_{\rm 0}$ at the upper edge of the astrophysical prior. All other
parameters of the agnostic model are unconstrained and the posterior is very
similar to the prior (see Appendix A for full posterior distributions for both
models).

Fig. 7 displays the posterior on the SMBHB chirp mass function for the two
models integrated over the redshift range $0<z<1.5$. Although the agnostic
model results in a loosely constrained mass function, the measured PTA signal
alone places interesting lower limits on the SMBHB binary merger rate in the
Universe. For example, we can say at 95% confidence that for each comoving
Gpc3, there have been at least 104 SMBHB mergers with ${\cal M}\approx
10^{7}$M⊙ since cosmic noon. When astrophysical priors weights in, the mass
function of merging SMBHBs is well constrained by the PTA signal and, as
expected from Fig. 6, it lies in the upper range of the astrophysically
informed prior. Within this model, the DR2new measured signal implies there
have been about 106 SMBHB mergers for each comoving Gpc3, with ${\cal
M}\approx 10^{9}$M⊙ since $z=1.5$. This points towards a very active merger
history of massive galaxies and a very efficient dynamical evolution of the
SMBHBs forming in the merger process.

For the astrophysically informed model, DR2new also provides interesting
information on several model parameters. This is because the astrophysical
prior considerably narrows down the range of signal amplitudes allowed by the
model, and the measured signal pushes towards the upper bound of this range.
This results in informative constraints on several key parameters, related in
particular to the SMBHB merger timescale and the SMBH-bulge mass relation. As
shown in Fig. 8, the SMBH merger timescale $\tau_{0}$ following galaxy pairing
must be shorter than $\approx 1\,{\rm Gyr}$ (90% confidence), with the data
mildly favouring shorter merger times for massive galaxies at low redshifts
(i.e. $\alpha_{\tau},\beta_{\tau}<0$). Moreover, the data favour a high
normalization of the SMBH-bulge mass relation ${\rm log}_{10}M_{*}\approx
8.4_{-0.32}^{+0.18}$, compared to a much wider prior range extending all the
way down to ${\rm log}_{10}M_{*}=7.8$. This is in line with the qualitative
analysis of Sec. 3.1, which showed that the signal is consistent with recent,
upward-revised, SMBH-galaxy relations. There is also a slight preference for a
high normalization of the pair fraction $f_{0}$ with a positive $z$
dependence, $\beta_{f}>1$. Despite the low $\gamma$ value favoured by the
data, indicative of a flatter spectrum compared to the canonical value
predicted by circular GW-driven binaries, SMBHB dynamics is largely
unconstrained, perhaps with a marginal preference for eccentric binaries
evolving in dense environments ($e_{0}$ and $\zeta_{0}$ posteriors slightly
rising towards the right bound of the prior).

Altogether, these results point towards efficient orbital decay of SMBHBs in
the aftermath of galaxy mergers, providing direct evidence that the ‘final
parsec problem’ is solved in nature and that compact sub-parsec SMBHBs must be
common in the centre of massive galaxies.

###  3.3 Implications for SAMs

We now explore the implication of this signal for the joint modelling of the
galaxy and SMBH formation and evolution by taking a close look at two state-
of-the-art SAMs: the model constructed by Barausse and collaborators
(Barausse, 2012; Sesana et al., 2014; Antonini et al., 2015; Klein et al.,
2016; Bonetti et al., 2018b; Barausse et al., 2020) and L-Galaxies (Henriques
et al., 2015; Izquierdo-Villalba et al., 2022). In this preliminary study, we
do not model the dynamical evolution of the binaries and we assume them to be
circular and GW driven, thus resulting in a characteristic strain spectrum
with $\alpha=-2/3$.

####  3.3.1 SAMs and SMBHB delays

Figure 9: Predictions for the GWB characteristic strain amplitude at
$f=1/10$yr from a range of SAMs published in the literature, assuming
quasicircular orbits and no environmental interactions (i.e. $\gamma=13/3$),
but different physical prescriptions for the delays (increasing from left to
right) between galaxy mergers and black hole mergers. The “no delays”, “medium
delays” and “long delays” models correspond respectively to the classes of
models (i), (ii) and (iii) defined in the text. The ranges account for the
finite resolution of the models. The shaded area is the DR2new 95% confidence
bound. More details about the models are provided in the text.

In Fig. 9, we show this comparison for the model of Barausse (2012) in its
original version (B12) and its subsequent evolutions, which were used to
produce the results of Klein et al. (2016) (K+16), Bonetti et al. (2018b)
(B+18) and Barausse et al. (2020) (B+20). Besides the specific SAM
implementation and (astro)physics, these models mainly differ for the physical
prescriptions describing the delays between galaxy and MBH mergers, with (i)
models “LS-nod (B12)”, “HS-nod (B12)”, “Q3-nod (K+16)”, “LS-nod-noSN (B+20)”,
“HS-nod-noSN (B+20)” and “HS-nod-SN-high-accr (B+20)” assuming no such delays
(except for the delays between the mergers of the halos and those of the
corresponding baryonic components)666Models “LS-nod-noSN (B+20)”, “HS-nod-noSN
(B+20)” and “HS-nod-SN-high-accr (B+20)” were not presented in B+20, but are
produced using the model of that paper, setting to zero the delays between
galaxy and MBH mergers (except for the dynamical friction timescale –
including tidal effects – between dark matter halos).; (ii) models “popIII-d
(K+16)”, “Q3-d (K+16)”, “LS-d (B+18)”, “HS-d (B+18)” additionally introducing
the effect of stellar hardening, triple MBH interactions and gas-driven
migration; and (iii) models “LS-noSN-d (B+20)”, “LS-SN-d (B+20)”, “HS-noSN-d
(B+20)” and “HS-SN-d (B+20)” accounting for even longer delays (including
large contributions from SMBHB separations of hundreds of pc). Note that the
labels “SN” (and “noSN”) refer respectively to a putative effect of SN
feedback on black hole accretion (and the absence thereof), while
“LS”/“popIII” and “HS”/“Q3” denote respectively light and heavy high-redshift
seeds for the black hole population.

The predictions are computed by summing the gravitational energy spectra of
all the SMBHBs in each model’s theoretical population, assuming quasi-circular
orbits. As a result, the spectrum has a slope of $\gamma=13/3$ and has no
cosmic variance (i.e. we do not account at this stage for the scatter from one
realization of the SMBHB population to another). The range shown for each
model represents the uncertainty due to the correction for the finite
resolution of the SAM’s merger tree. In more detail, the lower end of the
range represents a model’s prediction at the finite resolution, while the
upper end is the extrapolation – performed as described in Figure 4 of Klein
et al. (2016) – to infinite resolution. The lower arrow (pointing up) should
therefore be interpreted as a lower limit, while the upper arrow (pointing
down) should be understood as an upper bound (due to the uncertainty of the
extrapolation procedure). The extrapolation has not been performed for the
model HS-nod-SN-high-accr (B+20), for which we report only the (more robust)
prediction at finite resolution. The latter already agrees with the measured
amplitude, as a result of a higher accretion rate (by a factor $\sim 4$) for
SMBHs.

For two of the models in better qualitative agreement with the data (i.e. “HS-
nod-SN-high-accr (B+20)” and “HS-nod-noSN (B+20)”), we compare the predicted
signal with the measured DR2new free spectrum in Fig. 10. Unlike in the case
of Fig. 9, these predictions were obtained for multiple specific realizations
of the SMBHB population, following Sesana et al. (2008)777Note however that
following Sesana et al. (2008), and different from RSG15, we average over sky
position and binary inclination.. The probability distribution function
plotted in each bin represents the scatter of these multiple realizations, and
should therefore be interpreted as a “cosmic variance”. Similarly, in Fig. 11
we show the theoretical forecasts for $A(f=1/10{\rm yr})$ from a subset of the
models presented above (namely those in qualitative agreement with the data in
Fig. 9). These forecasts are obtained by fitting the model predictions (from
multiple realizations of the SMBHB population) in the first 9 frequency bins
with a power law with $\gamma=13/3$. The error bars represent the 95%
confidence regions (accounting for the scatter due to cosmic variance), while
the shaded area indicates the 95% confidence region of the posterior for
$A(f=1/10{\rm yr})$ (assuming again $\gamma=13/3$).

Figure 10: Binned spectrum of the predicted GWB amplitude for models “HS-nod-
SN-high-accr (B+20)” and “HS-nod-noSN (B+20)”. The distribution of the
predictions represents the scatter among different realizations of the SMBHB
population (“cosmic variance”). Also shown are power-law fits to the
predictions.

Overall, this qualitative comparison, while somewhat dependent on the details
of the SAM implementation, suggests that (i) large delays arising from the
dynamics of black hole pairs at large $\sim 100$ pc separations are
disfavoured, (ii) SMBHB mergers proceed efficiently after their host galaxies
have coalesced. Moreover, these results seem to suggest that (iii) accretion
onto SMBHs proceeds efficiently, possibly resulting in a larger local SMBH
mass function at high masses.

Figure 11: Predictions for $A(f=1/10{\rm yr})$ in various SAMs, obtained by
fitting the spectrum in the first 9 frequency bins with $\gamma=13/3$ for
multiple realizations of the SMBHB population. The error bars represent the
95% confidence interval for the predictions, and account for the scatter due
to cosmic variance. For each model (except for the boosted accretion model HS-
nod-SN-high-accr (B+20)), the higher prediction is the extrapolation to
infinite SAM resolution, while the lower one is the finite-resolution
prediction. The shaded area is the 95% confidence interval for the measurement
of $A(f=1/10{\rm yr})$, fixing $\gamma=13/3$. For HS-nod-SN-high-accr (B+20)
we only show the result uncorrected for resolution.

####  3.3.2 L-Galaxies

Figure 12: Predictions for the GWB characteristic strain amplitude at
$f=10/$yr-1 from a range of L-Galaxies semi-analytical model versions,
assuming that $h_{c}(f)\,{\propto}\,f^{-2/3}$. The error bars are computed
taking into account the cosmic variance. To this end, we divided the
Millennium box into 125 sub-boxes and we compute the GWB in each one. The
standard deviation provided by the 125 GWBs corresponds to the extension of
our error bars.

Next, we explore the implications that the EPTA results have for L-Galaxies
(Henriques et al., 2015; Izquierdo-Villalba et al., 2022), a sophisticated SAM
constructed on the dark matter merger trees extracted from the Millennium
simulation suite (Springel et al., 2005; Boylan-Kolchin et al., 2009). On top
of the galaxy physics, L-Galaxies features a comprehensive modelling for the
assembly of SMBHs, including gas accretion triggered by galactic mergers and
disc instabilities and dynamical evolution of SMBHBs within the galaxy merger
process. The latter accounts for dynamical friction (DF), stellar and gas
hardening and, eventually, GW emission. All of these processes are governed by
a set of parameters that are tuned to reproduce a vast array of observations
including, among others, the galaxy mass function and morphological
distribution, the quasar luminosity function and the SMBH-host galaxy
relations.

Izquierdo-Villalba et al. (2022) found that the standard L-Galaxies tuning
results in a GWB with ${\rm log_{10}}A=-14.9$ at $f_{0}=1{\rm yr}^{-1}$, lower
than that measured in DR2new. Here we perform a systematic investigation of
how the stochastic GWB at nanohertz frequencies depends upon the parameters
governing the physics of SMBHs and SMBHBs in the SAM. To this aim, we run
L-Galaxies in the following configurations:

  * •

std: standard configuration (Izquierdo-Villalba et al., 2022);

  * •

t_DF_x0.1: SMBH dynamical friction (DF) time reduced by a factor of ten;

  * •

NO_GAS_HARD: only stellar hardening;

  * •

growthDF_x10: accretion boosted by ten in the DF phase;

  * •

boostBH1: gas accretion doubled after galaxy mergers and disc instabilities;

  * •

boostBH2 gas accretion doubled after galaxy mergers and tripled after disc
instabilities.

  * •

boostBH1_growthDF_x10: adding accretion boost in the DF phase to model
boostBH1;

  * •

boostBH2_growthDF_x10: adding accretion boost in the DF phase to model
boostBH2;

Results are shown in Fig. 12. Changes to the dynamics of SMBHs appear to have
a minor effect on the amplitude of the GWB. While shortening the DF time
(t_DF_x0.1) allows more SMBHBs to merge within the Hubble time, the most
massive ones, which are responsible for the bulk of the GW signal, already
have short DF timescales, and the overall GW signal is only mildly increased.
Turning off gas hardening results in longer-lived SMBHBs that tend to merge at
lower redshifts, resulting in louder GW signals. The effect is, however,
negligible. Conversely, the tuning of gas accretion can significantly change
the masses of the SMBHBs, resulting in a larger impact on the GWB. Model
growthDF_x10 leaves the general population of SMBHs untouched, only boosting
the growth of those in the dynamical friction phase. This alone increases the
level of the GWB by a factor $\approx 1.5$ compared to model std. Finally, the
models boostBH1 and boostBH2 increase the gas accretion onto the whole
population of SMBHs, making the GWB a factor of 2.5 louder with respect to the
baseline model. Boosting accretion onto the whole population and in the
hardening phase further amplifies the expected GWB, reaching the upper bound
of the measured value (models boostBH1_growthDF_x10 and boostBH2_growthDF_x10
in the figure). Although these models can accommodate the GWB signal measured
in PTAs, the boosted accretion and subsequently larger SMBH masses can make it
harder to reproduce the observed SMBH mass and quasar luminosity functions
(Izquierdo-Villalba et al., 2022). Additionally, more work is required to find
a model that can reproduce all observational constraints in the light of the
PTA GW signal (Izquierdo-Villalba et al. in preparation).

###  3.4 Further considerations on the measured spectrum: eccentricity and
statistical biases.

The analyses presented so far give strong indications that the signal is
compatible with a cosmic population of SMBHBs swiftly coalescing in the
aftermath of galaxy mergers. The relatively flat slope of the measured
spectrum might be indicative of strong environmental coupling and high
eccentricities, although inference from the data is inconclusive in this
respect (see Fig. 8).

The eccentricity of SMBHBs is of particular relevance as it might carry
important information on the dynamical processes driving binary evolution at
sub-pc scales (see e.g. Roedig & Sesana, 2012). While gas-driven dynamics is
expected to saturate the binary eccentricity at a value $e\approx 0.4-0.6$
(Roedig et al., 2011; D’Orazio & Duffell, 2021), stellar hardening is known to
statistically increase eccentricity without any saturation point (Quinlan,
1996), potentially leading to extremely eccentric systems (Sesana, 2010). A
large binary eccentricity has two important implications for the
interpretation of the current data: it flattens the low-frequency spectrum and
it speeds up the SMBHB merger process, as inferred by the small $\tau_{0}$
derived in Sec. 3.2.

Figure 13: Orbital parameters (distance between the SMBHs, semi-major axis and
eccentricity) of a SMBHB formed in a representative $N$-body simulation of a
galactic merger with parameters drawn from progenitors of likely PTA sources
in the IllustrisTNG100-1 cosmological simulation. Mergers are selected from
the merger trees of the 100 most massive galaxies at $z=0$, based on galaxy
mass ratio (major mergers with mass ratio $1:4$ or higher) and redshift
($z\leq 2$). The dashed lines indicate the critical separation $a_{f}$ and the
corresponding eccentricity $e_{f}$ at the time in the evolution marking
approximately the end of the SMBH inspiral due to DF and the beginning of the
hardening phase. Dots represent $a$ and $e$ computed from the apoastron-
periastron separation of the two SMBHs before pairing in a bound binary.

Figure 14: Posterior distributions of the recovered GWB from injections on
synthetic data mimicking DR2new. Top panel: statistical offset in an ideal
dataset. The square marks the injected value and the blue contours are
1$\sigma$ and 2$\sigma$ of the recovered posterior. Bottom panel: effect of
high-frequency noise mismodeling on the recovery. The orange, blue and green
contours are respectively obtained when EFAC$=0.8,1,1.2$ are used for the
recovery (injected EFAC$=1$).

Low-redshift massive galaxies are generally gas-poor, and stellar-driven
hardening represents the main mechanism driving the evolution of the binaries
comprising the bulk of the PTA GW signal. Modelling the whole dynamical
evolution, from the first galaxy encounter to black hole pairing and final
coalescence, is theoretically and numerically challenging and has been the
subject of many studies (e.g. Preto et al., 2011; Khan et al., 2012, 2016;
Nasim et al., 2021; Gualandris et al., 2022). In particular, the binary
eccentricity is very sensitive to the initial orbits of the merging galaxies
(Gualandris et al., 2022). In ongoing work (Fastidio et al. in prep), we are
connecting the sub-pc dynamics of SMBHBs to the large-scale parameters of the
galactic encounters extracted from the IllustrisTNG100 simulation (Pillepich
et al., 2018). Preliminary results show that mergers of massive galaxies occur
preferentially on nearly radial orbits, potentially resulting in highly
eccentric binaries. Fig. 13 shows the orbital parameters of a SMBHB formed in
a high-accuracy $N$-body simulation of a representative galactic merger with
properties taken from a merger tree in IllustrisTNG100. Merger trees are
selected to represent likely PTA sources at low redshifts. The galactic merger
is followed from early times through the inspiral, pairing and hardening of
the SMBHs via a hybrid numerical scheme able to model the evolution self-
consistently from kpc to mpc scales (Dehnen, 2014). Despite the intrinsic
stochasticity of the processes driving binary formation and hardening (Nasim
et al., 2020), binaries tend to form with a large eccentricity, which then
further grows due to encounters with background stars, as also found by
scattering experiments (Sesana, 2010). Although more work is needed to
determine the distribution of expected binary eccentricities and current EPTA
data are not yet strongly informative, this pilot study shows the potential of
using these measurements in the near future to constrain the physics of galaxy
and SMBH mergers.

When drawing conclusions on the physical properties of the sources of the GW
signal, it is useful to bear in mind not only that the constraints on the
spectrum are relatively loose (see Fig. 1), but also that the measured
parameters can be subject to statistical and systematic biases. To address
this, we are conducting an extensive campaign of simulations by injecting and
recovering different types of signals in synthetic PTAs mimicking the
properties of the EPTA DR2new dataset (Valtolina et al., 2024). We generated
individual noises for 25 pulsars using the maximum likelihood values of the
measured white noise and drew the red noise and dispersion measure parameters
from the posterior distribution of the customised noise models (EPTA and InPTA
Collaborations, 2023a). We simulated TOAs from multi-frequency observations
and added a GWB spectrum from an astrophysical population of circular SMBHBs
producing a nominal GWB with $A_{1{\rm yr}^{-1}}=2.5\times 10^{-15}$,
consistent with the DR2new measurement at $\gamma=13/3$. We performed 100
experiments by changing the specific noise realization and the sampling of the
injected SMBHB population. The analysis was carried using the ENTERPRISE
software package (Ellis et al., 2020).

Two illustrative examples of injection-recovery mismatch are shown in Fig. 14.
The top panel shows one of the GWB recoveries. Although the injection did not
present particular features (e.g. loud CGWs), for this specific noise
realization, the recovered signal has a shallow slope with a median value of
$\gamma=3.10$. Similar cases have been observed when injecting a pure
$\gamma=13/3$ power law with the createGWB function of libstempo (Vallisneri,
2020). This shows that even with an ideal setup when simultaneously fitting
multiple parameters (102 in this case) in a complex problem, the stochastic
properties of the noise can easily bias the recovered signal, especially if
the S/N is low (S/N$\approx 3.5$ for DR2new). Multiple injections with the
realistic GWB model and createGWB show systematic biases of recovery of the
realistic GWB signal, when modelled with an ideal power law. This is explored
in detail in a follow-up work. In the bottom panel of Fig. 14, we show how the
presence of some extra high-frequency noise unaccounted for in the MSP noise
model can also influence the recovery of the parameters. The setup is the same
as in the left panel, but we simulate high-frequency noise mismatch by setting
different values of EFAC$=0.8,1,1.2$ in the recovery. Although the posterior
of the signal amplitude is hardly affected, $\gamma$ can shift significantly
depending on whether the high-frequency noise is slightly over- or under-
estimated. While these are only two examples, they highlight the complexity of
PTA measurements and invite us to be cautious when drawing conclusions that
might strongly be influenced by potential biases in the recovered signal.

##  4 Implications II: physics of the early Universe

Although a GWB generated by an ensemble of the putative SMBHBs is the most
plausible source of the observed common-red noise process in pulsar data, more
exotic explanations are possible, such as signals originated in the early
Universe. The various possible types of cosmological backgrounds of GWs
associated with early Universe physics are reviewed in Caprini & Figueroa
(2018) and are found to be stochastic. Similarly to the traditional case
invoking SMBHBs, the angular spatial correlation for those scenarios follows
the HD curve. However, the spectral shapes of the predicted GW spectra are
generally different, which can help to disentangle between different types of
backgrounds. In this work, we focus on four possible scenarios:

  1. 1.

an inflationary GWB from the amplification of quantum fluctuations of the
gravitational field,

  2. 2.

a GWB from a network of cosmic string loops,

  3. 3.

a GWB from vortical (M)HD turbulence at the QCD energy scale,

  4. 4.

a scalar-induced GWB arising from inflationary scalar perturbations at the 2nd
order in perturbation theory.

Given the low significance of the detected signal and the limited number of
probed frequency bins due to the short timespan of the data, one cannot
currently perform a reliable model selection. Therefore, throughout the
section, we consider these scenarios separately and assume that each of them
can fully explain the detected signal independently. Analysis invoking more
complex models with simultaneous fits for multiple scenarios as well as
opportunities to disentangle between those (e.g. Goncharov et al., 2022;
Kaiser et al., 2022) will be considered in a number of future works.

###  4.1 Implications on a stochastic background of primordial (inflationary)
gravitational waves

Here we address the GWB possibly generated during inflation (Grishchuk, 1975;
Rubakov et al., 1982; Starobinskii, 1985; Fabbri & Pollock, 1983; Abbott &
Wise, 1984). In the standard inflationary scenario, tensor quantum vacuum
fluctuations of the metric are amplified by the accelerated expansion, leading
to a GWB as they subsequently re-enter the horizon during the radiation (or
matter) era. The cosmic microwave background (CMB) and Big Bang
Nucleosynthesis (BBN) provide precise measurements of the radiation energy
density, from which one can derive weak upper bounds on the amplitude of such
a GWB (see e.g. Caprini & Figueroa, 2018, and references therein).
Furthermore, tensor metric perturbations lead to CMB temperature anisotropies
and polarisation at large angular scales (Sachs & Wolfe, 1967; Starobinskii,
1985; Kosowsky, 1996; Allen & Koranda, 1994). Since the anisotropies and
polarisation detected so far are due to scalar perturbations, it is possible
to constrain the energy density of a GW background by placing an upper limit
on the tensor-to-scalar ratio $r$ at CMB scales: recent upper bounds are given
e.g. in Tristram et al. (2022); Galloni et al. (2023). Another parameter to
consider is the tensor spectral index $n_{T}$ of the tensor perturbations. In
the context of slow-roll single-field inflation, these two parameters are
linked via the consistency relation $r=-8n_{T}$ . By fixing the consistency
relation, Tristram et al. (2022) finds $r<0.032$ at 95% CL, while by relaxing
it, Planck Collaboration et al. (2020b) finds $r<0.076$ and $-0.55<n_{T}<2.54$
at 95% CL.

Within the slow-roll consistency relation, the GWB spectral slope is therefore
slightly red-tilted, causing this signal to be out of reach of most current
and planned GW detectors such as PTAs, LISA, aLIGO, aVirgo or the Einstein
Telescope. On the other hand, it is fair to consider that the spectral slope
could vary over the large frequency span ranging from CMB scales to those
probed by GW detectors (Lasky et al., 2016). Inflationary scenarios breaking
the consistency relation at small scales might therefore produce a detectable
GWB, if they lead to a blue-tilted spectrum. In this case, PTAs, LISA and
ground-based devices can place upper limits on $n_{T}$ (see e.g. Abbott et
al., 2017). It is interesting to investigate which kinds of processes could
give rise to a blue-tilted GW background while still obeying CMB constraints
at large scales. One possibility is the presence, after inflation, of a stiff
component in the Universe, with an equation of state $w>1/3$ (Boyle &
Buonanno, 2008; Giovannini, 1998; Boyle & Steinhardt, 2008). The enhancement
of the tensor spectra can also be produced during inflation thanks to
processes such as, for example, (i) particle production during inflation (see
e.g. Sorbo (2011); Barnaby et al. (2012); Cook & Sorbo (2013); Anber & Sorbo
(2012); Bartolo et al. (2016)) (ii) enhancement of tensor perturbations for
example by spectator fields, or space-dependent inflation (see e.g. Bartolo et
al. (2007); Biagetti et al. (2013, 2015); Fujita et al. (2015)) (iii) modified
gravity theories such as $f(R)$ or Horndeski gravity (Horndeski (1974);
Sotiriou & Faraoni (2010)) and iv) enhanced scalar perturbations at small
scales and/or primordial black holes, which are treated in Sec. 4.4.

####  4.1.1 Analysis

Similarly to what was done in Lasky et al. (2016), we constrain the key
parameters defining the GWB, $r$ and $n_{T}$, while being agnostic on the
underlying process generating the blue-tilted spectrum. If we assume that the
common quadrupolar red noise signal present in EPTA data is of inflationary
origin, these two parameters can be estimated using the DR2new dataset. Note
that the spectral index $n_{T}$ is expected to vary with the frequency scale
considered (see e.g. Giarè & Melchiorri, 2021; Giarè et al., 2023; Auclair &
Ringeval, 2022). However, for simplicity, here we consider a constant $n_{T}$
all the way from CMB scales to those corresponding to the (narrow) EPTA
frequency band. It is then possible to approximate the fractional
characteristic GW energy density using (Lasky et al., 2016; Caprini &
Figueroa, 2018)

$\displaystyle\Omega_{\rm GW}(f)=\frac{3}{128}\Omega_{\rm
rad}\,r\,\mathcal{P}^{*}_{\mathcal{R}}\left(\frac{f}{f_{*}}\right)^{n_{T}}\left[\frac{1}{2}\left(\frac{f_{\rm
eq}}{f}\right)^{2}+\frac{16}{9}\right]$ (14) $\displaystyle\approx 1.5\times
10^{-16}\left(\frac{r}{0.032}\right)\left(\frac{f}{f_{*}}\right)^{n_{T}},$

where the second line is valid in the PTA frequency band, and has been
obtained by setting $h^{2}\Omega_{\rm rad}=2.47\cdot 10^{-5}$ with $h=0.67$,
the amplitude of the scalar spectrum $\mathcal{P}^{*}_{\mathcal{R}}=2\cdot
10^{-9}$, and $f_{*}\approx 7.7\times 10^{-17}$ Hz related to the CMB pivot
scale $k_{*}=0.05/{\rm Mpc}$ (Planck Collaboration et al., 2014). $f_{\rm eq}$
denotes the frequency entering the horizon at matter-radiation equality.

We then use the nine lowest frequency posteriors of the RMS free spectrum
shown in Fig. 1 (see Moore & Vecchio, 2021; Lamb et al., 2023; Leclere et al.,
2023, for details on the method) to fit the inflationary spectrum of Eq. 14
and obtain posteriors on $\log_{10}r$ and $n_{T}$. Results are reported in
Fig. 15. Note that, since $\gamma=5-n_{T}$, the correlation between the
amplitude and spectral index of the signal is compatible with Fig. 5. The 90%
credible (symmetric) intervals are $\log_{10}r=-12.18^{+8.81}_{-7.00}$ and
$n_{T}=2.29^{+0.87}_{-1.11}$. The obtained value of $n_{T}$ corresponds to a
PSD spectral index of $\gamma\simeq 2.7$, as in Fig. 1. The excessively small
value of $r$ is a consequence of the simplistic parameterisation of Eq. 14,
which assumes a constant $n_{T}$ at all scales. The fractional energy density
spectrum obtained from the maximum a posteriori parameter values is plotted in
Fig. 17.

We have so far considered a primordial background to be the only source of GWs
in our data. We now recall that the most plausible and loud source of a GW
background at these frequencies remains that of a SMBHB background. It is
therefore likely that any signature for a cosmological background needs to be
considered in parallel with a SMBHB background, or in this case more
accurately termed ‘foreground’. Kaiser et al. (2022) have explored the
likelihood of detecting a cosmological background in the presence of a SMBHB
foreground using simulations, and found that the shallower the slope of the
cosmological background (for example $\gamma=4$ as opposed to $\gamma=5$), the
harder it is to detect (and the longer it takes, possibly more than 20 years).
According to these simulations, this does not bode well for an even shallower
slope, similar to the one detected in DR2new with a possible $\gamma<3$.

Here we explore a superposition of these two backgrounds in the DR2new
dataset. Considering a two-component GWB for the common red noise model, we
place constraints on $\log_{10}r$ for given values of $n_{T}$ spanning the
range $[-1,3]$. In this case, our null hypothesis is a GWB from a population
of GW-driven circular SMBHBs parameterised only by the PSD amplitude
$\log_{10}A$ of Eq. (3) (we fix $\gamma=13/3$). We run several analyses with a
fixed $n_{T}$ for the inflationary background, sampling over
$(\log_{10}r,\log_{10}A)$. For each of the $n_{T}$ values, we obtain a
distribution for $\log_{10}r$ and take the 95% quantile as an upper bound. As
found in Lasky et al. (2016), $n_{T}$ and the $\log_{10}r$ upper bounds are
related with good precision by a linear relation:

$n_{T}=a\log_{10}\left(\frac{r}{0.032}\right)+b.$ (15)

Our analysis gives $a=-0.16$ and $b=0.70$, which is comparable to the forecast
values given in Lasky et al. (2016) (note that they normalise $r$ to 0.11).

####  4.1.2 Discussion

From the analysis of the DR2new dataset above, we have obtained credible
intervals for the tensor-to-scalar ratio $r$ and the spectral index $n_{T}$.
This was performed assuming that reheating is instantaneous, and that
inflation is followed directly by the radiation-dominated era, for which the
equation of state parameter of the Universe is $w=1/3$. Under this assumption,
one finds that the best-fit value for the tensor spectral index is
$n_{T}=5-\gamma\simeq 2.3$, which is directly linked to the best-fit PSD
spectral index $\gamma\simeq 2.7$. This high value of $n_{T}$ is not
consistent with slow roll inflation. However, if inflation is followed by a
stage in which $w\neq 1/3$, the relation between the PSD spectral index
$\gamma$ and the primordial tensor spectral index $n_{T}$ changes to
(Arzoumanian et al., 2016; Caprini & Figueroa, 2018)

$\gamma=5-n_{T}+\frac{2(1-3w)}{3w+1},$ (16)

again with $\gamma\simeq 2.7$. If a stiff fluid component ($w>1/3$) were to
dominate the Universe for a finite amount of time after inflation, the last
term in Eq. 16 would be bounded between 0 and $-2$. Hence, $n_{T}\gtrsim 0.3$,
meaning that even allowing for the presence of a stiff component after
inflation, it does not seem possible to explain the common red noise in the
context of slow roll inflation ($n_{T}\simeq 0$) for the best fit value
$\gamma\simeq 2.7$. However, by broadening the range of possible values to
$\gamma\geq 3$, $n_{T}\simeq 0$ does become compatible with the common red
noise.

Figure 15: 2D posteriors of the tensor-to-scalar ratio (in $\log_{10}$) and
the fractional energy density spectral index $n_{T}$ in the PTA frequency
range. The 68% and 95% credible regions are displayed. The black dashed line
represents the tensor-to-scalar ratio upper bound found in Tristram et al.
(2022) assuming single-field slow-roll inflation.

###  4.2 Implications on a background of cosmic strings

Cosmic strings are line-like topological defects that may form after a
symmetry-breaking phase transition in the early Universe (Kibble, 1976;
Hindmarsh & Kibble, 1995); they are generic predictions of most Grand
Unification Theories scenarios (Jeannerot et al., 2003). These 1D objects are
characterised by the string tension $G\mu$ (or equivalently their energy per
unit length) which is related to the energy scale of the phase transition.

Overall, cosmic strings combine relativistic velocities and large energy
densities, making them natural sources of GWs. These GWs may take the form of
bursts from cusps, kinks and kink-kink collisions on the loops (Damour &
Vilenkin, 2001), and have been searched for in the LIGO/Virgo/KAGRA (LVK)
detectors (Abbott et al., 2018, 2021). Cusps are points on the string which,
in the Nambu-Goto approximation, propagate at the speed of light, and the
string doubles back on itself. On the other hand, kinks are discontinuities in
the tangent vector of a string and are formed at each intercommutation of
strings. The future space-based detector LISA will be sensitive to cosmic
string bursts for tensions as low as $G\mu\sim 10^{-11}$ (Auclair et al.,
2023a). Most notably, in the event that LISA detects GW bursts from cosmic
strings, they will likely repeat multiple times a year due to the periodic
nature of the cosmic string loops (Auclair et al., 2023b).

The incoherent sum of these GW bursts also produces a stochastic GWB (SGWB),
which has been looked for by LVK (Abbott et al., 2018, 2021). LISA is expected
to be able to detect a SGWB from cosmic strings with tensions as low as
$G\mu\sim 10^{-17}$ (Auclair et al., 2020). The strings background has already
been looked for in EPTA (Sanidas et al., 2012; Leclere et al., 2023), in
NANOGrav (Blasi et al., 2021; Ellis & Lewicki, 2021) and PPTA (Bian et al.,
2022; Chen et al., 2022). These earlier analyses consistently found a good
agreement between the common uncorrelated red signal (CURN) present in the
data and a SGWB from cosmic strings with tensions $G\mu\sim
10^{-10}-10^{-11}$.

####  4.2.1 Description of the models

Since the SGWB is sourced by sub-horizon cosmic string loops, the central
quantity in our analysis is the loop density distribution,
$\textbf{n}(\ell,t)$, where $\ell$ is the invariant length of the loop
(defined by its total energy divided by $\mu$) and $t$ is cosmic time. We
focus on the most up-to-date loop distributions, calibrated on large scales
with Nambu-Goto simulations (Lorenz et al., 2010; Blanco-Pillado et al.,
2014). Loops are formed from the intercommutation of super-Horizon infinite
strings, with a maximal size $\ell\approx 0.1t$ at time $t$. (Note that the
simulations assume an intercommutation probability of $1$, as is the case for
field theory strings.) The two models — which we refer to as BOS (Blanco-
Pillado et al., 2015) and LRS (Lorenz et al., 2010), after the names of their
authors — differ in their treatment of small loops, particularly on scales at
which gravitational radiation (not included in numerical simulations) can have
important effects. Indeed, compared to the BOS model, the LRS model has an
additional population of small loops, which leaves an imprint in the SGWB
spectrum (Auclair, 2020). The explicit expressions for the two loop
distributions are given in Abbott et al. (2018), and both have been considered
by the LVK (Abbott et al., 2021) and LISA (Auclair et al., 2020)
collaborations. The third observing run of LVK placed constraints on $G\mu$,
based on the non-detection of a SGWB. These are $G\mu\lesssim 10^{-8}$ for BOS
and $G\mu\lesssim 10^{-15}$ for LRS. At the frequency of ground-based
detectors, the SGWB signal is produced by loops formed during the radiation
era. At low PTA frequencies, the SGWB signal is dominated by larger loops,
namely those formed at recent times, in transition from the radiation to
matter era and also in the matter era (Ringeval & Suyama, 2017; Auclair, 2020;
Leclere et al., 2023).

Other than $G\mu$, a further parameter appears in $\textbf{n}(\ell,t)$, namely
$\Gamma$, which describes the rate at which loops lose energy through
gravitational radiation: $\dot{\ell}=-\Gamma G\mu$. If $\Gamma$ is large,
fewer loops are present in the distribution since loops decay more rapidly. On
the other hand, GWs are emitted more intensively: the final effect is a
combination of these two, which impacts the SGWB of the BOS and LRS models in
subtly different ways (Auclair, 2020). The value of $\Gamma$ depends on the
properties of loops, and in particular on how many cusps, kinks, and kink-kink
they contain as (Damour & Vilenkin, 2001; Siemens et al., 2007; Ringeval &
Suyama, 2017)

$\Gamma=N_{c}\Gamma^{c}+N_{k}\Gamma^{k}+\frac{N_{k}^{2}}{4}\Gamma^{kk},$ (17)

where $N_{c,k}$ is the number of cusps/kink events per oscillation period of
the loop, $N_{k}^{2}/4$ the number of kink-kink collisions, and

$\Gamma^{c}=\frac{3(\pi
g_{1}^{c})^{2}}{2^{1/3}g_{2}^{2/3}},\>\Gamma^{k}=\frac{3(\pi
g_{1}^{k})^{2}}{2^{2/3}g_{2}^{1/3}},\>\Gamma^{kk}=2(\pi g_{1}^{kk})^{2}\;,$
(18)

where $g_{2}=\sqrt{3}/4,g_{1}^{c}\approx 0.85,g_{1}^{k}\approx
0.29,g_{1}^{kk}\approx 0.1$.

In this paper, we consider 2 cases. For the first model (i) $N_{c}=2$ and
$N_{k}=0$, for which $\Gamma=57$ (Leclere et al., 2023), a value close to that
observed in numerical simulations of individual loops (Vachaspati & Vilenkin,
1985; Blanco-Pillado & Olum, 2017). Therefore, the only free parameter for
this model is $G\mu$. As for the second model (ii), we fix $N_{c}=1$ and allow
$N_{k}$ to vary between 1 and 200 (with a uniform prior) so as to account for
theoretical uncertainties on the initial number of kinks at loop creation, and
on the efficiency of the gravitational backreaction that should smooth out the
loops (see the LVK analysis, Abbott et al., 2021, for a similar approach).
Therefore, this is a two-parameter model $(G\mu,N_{k})$.

The fractional energy density of the SGWB per logarithmic interval of
frequency is given by

$\Omega_{\rm{GW}}(t_{0},f)=\frac{16\pi(G\mu)^{2}}{3H_{0}^{2}}\sum_{b}\frac{N_{b}\Gamma^{(b)}}{\zeta(q_{b})}\\\
\times\sum_{n=1}^{+\infty}\int\frac{n^{1-q_{b}}\mathrm{d}\,{z}}{(1+z)^{5}H(z)}\textbf{n}\left[\frac{2n}{(1+z)f},t(z)\right],$
(19)