\section{Introduction}
\label{sec:intro}

Observations of \acp{GW} from \acp{CBC} by the LIGO~\cite{Harry:2010zz} and Virgo~\cite{TheVirgo:2014hva} detectors have brought tests of \ac{GR} in the strong-field regime to hitherto unachievable levels~\cite{TheLIGOScientific:2016src,Abbott:2018lct,LIGOScientific:2019fpa,Isi:2019asy}.
This is fundamentally dependent on the detailed knowledge about the structure of the \acp{GW} emitted from a binary of \acp{CO} (\acp{BH} or \acp{NS}) stemming from decades of analytical~\cite{Blanchet:1995ez, Blanchet:2004ek, Blanchet:2005tk, Blanchet:2013haa, Buonanno:1998gg, Buonanno:2000ef, Damour:2008qf, Damour:2009kr, Barausse:2009xi, Damour:2016bks} and numerical~\cite{Pretorius:2005gq, Campanelli:2005dd, Baker:2005vv, Boyle:2019kee} studies of \acp{GW} from binary systems of \acp{CO}.

So far, the majority of theories of gravity beyond \ac{GR} are unable to construct predictions for \acp{GW} emitted by coalescing binaries with generic \acp{CO} (but see~\cite{Yagi:2011xp,Berti:2013gfa,Berti:2015itd,Hirschmann:2017psw,Okounkova:2017yby,Witek:2018dmd,Okounkova:2018pql,Loutrel:2018ydv,Okounkova:2019dfo,Torsello:2019tgc,Okounkova:2019zjf,Okounkova:2020rqw,Julie:2020vov, Witek:2020uzz} for the status of current efforts), 
hence tests of \ac{GR} are generically formulated as consistency tests only, where the primary approach is to introduce \textit{ad-hoc} modifications of the \ac{GR} waveforms.
This can be constructed to test different regions and functional dependencies of the overall waveform~\cite{Agathos:2013upa, Li:2011cg, Li:2011vx, Cornish:2011ys, Meidam:2017dgf, Abbott:2018lct}, such as deviations from the analytical coefficients of the \ac{PN} expansion~\cite{Blanchet:1995ez, Blanchet:2004ek, Blanchet:2005tk, Blanchet:2013haa}, which has been successful when investigating constraints on each included \ac{PN}-order separately~\cite{TheLIGOScientific:2016src, TheLIGOScientific:2016pea, Abbott:2017vtc, Abbott:2018lct, LIGOScientific:2019fpa,Isi:2019asy}.
These constraints can later be mapped onto bounds on specific alternate theories of gravity~\cite{Chatziioannou:2012rf,Yunes:2013dva,Yunes:2016jcc,Nair:2019iur}, something which in turn highlights a potential flaw of this approach.
Since the \ac{PN}-coefficients themselves depend on the specific properties of the source's \acp{CO}, like their masses, it would be reasonable to also assume any deviations from the \ac{GR}-predicted values to also be source dependent. 
If a hypothetical theory modifies \ac{BH}-spin behaviour, but not any mass parameters, a general \ac{PN}-deviation would be different for two binaries with the same \ac{BH} spin magnitudes but different mass ratios.
This is not accounted for in most current analyses~\cite{TheLIGOScientific:2016src, TheLIGOScientific:2016pea, Abbott:2017vtc, Abbott:2018lct, LIGOScientific:2019fpa} (but see ~\cite{Isi:2019asy} for a more general approach) and could lead to misinterpreted inference if any deviation from \ac{GR} was observed~\cite{Zimmerman:2019wzo,Purrer:2019jcp}.
In addition, the strength of this class of tests is reduced when more than one \ac{PN}-term is simultaneously allowed to vary, where the addition of a large number of unconstraining degrees of freedom generates an overall null gain in information about any of the included terms (cf. Fig 7 of~\cite{TheLIGOScientific:2016src}).

In this letter we implement a null test of \ac{GR}, probing the validity of the current knowledge about \ac{GR}, and specifically its nonlinear behaviour originating from \ac{GW} tail effects~\cite{Thorne:1980ru, Blanchet:1987wq, Blanchet:1992br, Blanchet:1993ec, Tanaka:1993pu, Blanchet:1994ez, Blanchet:1997jj}, with the mathematical constant $\pi$ treated as a variable.
$\pi$ can here be considered as a universal parameter across all \ac{GW} observations of \acp{CBC}~\cite{Zimmerman:2019wzo}, and simultaneously tests 4 (out of the included 8) \ac{PN}-orders.
This enables an unprecedentedly powerful test, as it is both theory-agnostic and conceptually generalisable to probe a population of \acp{GW} through a quantity that is formally consistent across independent observations (while also being comparatively inexpensive computationally).
Throughout this letter we denote the true value of $\pi$ as $\ensuremath{\pi^\mathrm{T}}=3.141592653\ldots$~\cite{OEISpi}, a number which has been independently evaluated through several methods~\cite{Ramaley1969,Bailey1997,Galperin2003,Dumoulin:2014,Arndt2001,yCruncher}.
We assume $G = c = 1$.

\section{Method}
\label{sec:method}

The \ac{GW} signal from a \ac{CBC} can be generally expressed in the form
\begin{equation}
\tilde{h}(\theta,f) = A(\theta,f)e^{i\Psi(\theta,f)}\;,
\label{eq:htilde}
\end{equation}
where $\tilde{h}(f)$ is the emitted \ac{GW} strain in the frequency domain, with amplitude $A(\theta,f)$ and phase $\Psi(\theta,f)$ being functions of the source parameters $\theta$, e.g. \ac{CO} masses $m_{1,2}$, spin vectors $\vec{S}_{1,2}$ and tidal deformabilities $\Lambda_{1,2}$ (we fix $\Lambda_{1,2} = 0$ for \acp{BH}).
When the two \acp{CO} are sufficiently separated, for an orbital velocity $u \ll 1$ with $u = \ensuremath{\pi^\mathrm{T}} M f$ and $M$ being the binary's total mass, Eq.~\eqref{eq:htilde} can be described accurately through a \acl{PN} expansion in $u$.
Under the \ac{SPA}~\cite{Sathyaprakash:1991mt,Cutler:1994ys,Apostolatos:1994mx,Poisson:1995ef,Droz:1999qx,Buonanno:2009zt}, $\Psi(\theta,f)$ is given for the \ac{TF2} model~\cite{Sathyaprakash:1991mt,Damour:2000zb, Damour:2002kr, Arun:2004hn, Buonanno:2009zt} as
\allowdisplaybreaks
\begin{equation}
\begin{split}
\Psi_{\text{TF2}} &= \, 2 \ensuremath{\pi^\mathrm{T}} f t_c - \varphi_c - \ensuremath{\pi^\mathrm{T}}/4\\
                        & + \frac{3}{128 \, \eta}u^{-5/3} \sum_{i=0}^{7}
\left(\varphi_i + \log(u)\varphi^l_i\right) u^{i/3}\;,
\end{split}
\label{eq:TF2phase}
\end{equation}
where $t_c$ and $\varphi_c$ are the overall time and phase defined at coalescence and $\eta = m_1m_2/(m_1 + m_2)^2$ is the symmetric mass ratio.
The expansion coefficients $\varphi_i$ and $\varphi^l_i$ are then given as functions of $\theta$~\cite{Sathyaprakash:1991mt,Buonanno:2009zt, Blanchet:2013haa, Bohe:2013cla, Poisson:1997ha, Arun:2008kb, Mikoczi:2005dn, Khan:2015jqa}.
In this formalism, multiples of $\pi$ appear in $\varphi_3$, $\varphi_5$, $\varphi^l_5$, $\varphi_6$ and  $\varphi_7$.
To also capture the post-inspiral \ac{PN} description of a \ac{CBC} signal, we employ the \textsc{IMRPhenomPv2} model, where the analytical inspiral description from \ac{TF2} is smoothly extended with a phenomenological description of a \ac{BBH} merger-ringdown section~\cite{Husa:2015iqa,Khan:2015jqa} together with an effective-precession treatment~\cite{Hannam:2013oca}.
For the \ac{BNS} events, this is further extended with a description of \ac{NS} matter effects~\cite{Dietrich:2018uni,Dietrich:2017aum}.
Neither extension depends on the variable $\pi$, as the extensions are phenomenological rather than analytical in their nature.

We note that $\pi$ is included in the orbital velocity $u$, originating in a conversion from angular to linear orbital frequencies.
From Eq.~\eqref{eq:TF2phase}, $u$ is already strongly constrained at the leading-order phase term and can be taken as known to sufficiently high precision in \ac{GR}.
Since we are here interested in specifically probing the \acl{PN} formalism, expressed through the $\varphi_i$ and $\varphi^l_i$ coefficients, we fix $\pi = \ensuremath{\pi^\mathrm{T}}$ in $u$ throughout. 
Similarly, the $\pi$ in $2 \ensuremath{\pi^\mathrm{T}} f t_c$, originating from a Fourier transform of the time-domain \ac{GW} signal, and the factor $\ensuremath{\pi^\mathrm{T}}/4$ appears in Eq.~\eqref{eq:TF2phase} out of convention.
As both describe an overall phase shift, perfectly degenerate with the variables $t_c$ and $\varphi_c$ respectively,  we fix those two $\pi = \ensuremath{\pi^\mathrm{T}}$ in this analysis.

It is important to note that the appearance of $\pi$ in the \ac{PN}-coefficients follows purely from definitions of mathematics itself~\cite{Blanchet:1997jj,Foffa:2011ub}, e.g. through the use of known identities to evaluate integrals of a specific form (cf. Eq. 5.4 of~\cite{Blanchet:1997jj}), and does not depend on the specific assumptions of \ac{GR} as a theory of gravity.
This formally justifies treating $\pi$ in all $\varphi_i$ and $\varphi^l_i$ coefficients as fundamentally the same quantity, and also treating it as a universal parameter across multiple independent \ac{CBC} observations~\cite{Zimmerman:2019wzo}.

We also note that the \ac{PN}-orders where $\pi$ appears are primarily describing so called \ac{GW} tail effects~\cite{Thorne:1980ru, Blanchet:1987wq, Blanchet:1992br, Blanchet:1993ec, Tanaka:1993pu, Blanchet:1994ez, Blanchet:1997jj}, where the outgoing \acp{GW} backscatter off the (approximately) static spacetime of the \ac{CBC} source.
Tail effects are an inherently nonlinear behaviour present in \ac{GR}, hence the use of a variable $\pi$ can directly probe the validity of the nonlinear terms expressed through the \ac{PN}-representation of \ac{GR} itself.

While the analysis with a variable $\pi$ formally is an extension of \ac{GR}, we do not argue that results presented here are direct suggestions for alternative theories of gravity.
Instead, we interpret this study primarily as a strong null test, validating the current understanding of \ac{GR} through a multi-order probe of the \ac{PN}-formalism.

Finally, we acknowledge that the \ac{GW} detectors, as well as the data they record, are constructed and calibrated for $\pi = \ensuremath{\pi^\mathrm{T}}$ only.

\subsection{Bayesian methods}
\label{sec:bayesian}

We explore the parameter space $\theta$ defined by the \ac{CBC} models using Bayes' theorem to infer the posterior \ac{PDF}:
\begin{equation}
\label{Eq:BayesTheorem}
p(\theta|d, H) = \frac{p(\theta|H)p(d|\theta,H)}{p(d|H)}\;,
\end{equation}
where $p(\theta|H)$ is the prior \ac{PDF} of $\theta$ given the model $H$, $p(d|\theta,H)$ is the likelihood of observing the data $d$ assuming $\theta$ and $p(d|H) = \int p(\theta|H)p(d|\theta,H)d\theta$ is the evidence for $H$.
We preform Bayesian inference using the \textsc{LALInference}\xspace package~\cite{Veitch:2014wba,LALInference-code, lalsuite}, following the analysis configuration from~\cite{LIGOScientific:2018mvr} which includes a fixed noise \ac{PSD} (defined for the analysed data $d$ and generated as a median \ac{PSD} using \textsc{BayesWave}~\cite{Littenberg:2014oda,Cornish:2014kda,Chatziioannou:2019zvs,GWTC1_PE_release,GW190425_PE_release,GW190412_PE_release}) and marginalisation over uncertainties in the calibration of $d$~\cite{SplineCalMarg-T1400682,GWTC1_PE_release,GW190425_PE_release,GW190412_PE_release,Cahillane:2017vkb}.
All GW events are analysed using publicly available data~\cite{GWOSC_O1, GWOSC_O2, GWOSC_GW190425, GWOSC_GW190412, Vallisneri:2014vxa, Abbott:2019ebz}.

We assume prior choices consistent with those used in~\cite{Abbott:2018exr,LIGOScientific:2018mvr,Abbott:2020uma,LIGOScientific:2020stg}.
For the two \acp{BNS}, we perform only analyses with $|\vec{S}_{1,2}| \leq 0.05$, and parametrize the \ac{NS} tidal deformability following the \acl{EoS} independent relations from~\cite{Chatziioannou:2018vzf}.
We assume a prior distribution for $\pi$ that is uniform between $-20 \leq \pi \leq 20$.

As $\pi$ can be considered a formally universal parameter, as defined by~\cite{Zimmerman:2019wzo}, it is trivial to evaluate joint constraints on $\pi$ from a set of $N$ individual observations by multiplying the $1D$ likelihood distributions (marginalised over all other parameters), dividing by one instance of the common prior and normalising the resulting posterior \ac{PDF}.

Finally, as \ac{GR} is nested inside the model which allows for a variable $\pi$ it is possible to compute a \ac{BF} in favour of \ac{GR}, more directly where $\pi = \ensuremath{\pi^\mathrm{T}}$, using the Savage-Dickey density ratio ~\cite{dickey1971, BF_SavageDickey} as
\begin{equation}
\label{Eq:SavageDickey}
\mathrm{BF} = \frac{p(\pi = \ensuremath{\pi^\mathrm{T}}|d, H) }{p(\pi = \ensuremath{\pi^\mathrm{T}}|H)}\; ,
\end{equation}
i.e. the ratio of the posterior and prior \acp{PDF} evaluated at $\pi = \ensuremath{\pi^\mathrm{T}}$.

\section{Astrophysical measurement of $\pi$}
\label{sec:realpi}

The \ac{LVC} has so far, from its \ac{O1}, \ac{O2} and \ac{O3}~\cite{Aasi:2013wya}, confirmed 13 \ac{GW} observations, 2 \acp{BNS} and 11 \acp{BBH}~\cite{LIGOScientific:2018mvr,Abbott:2020uma,LIGOScientific:2020stg}.
Whereas other studies restrict themselves to high-significance events only \cite{LIGOScientific:2019fpa,Isi:2019asy}, primarily due to computational restrictions, the analysis presented here is easily extendable to and informed by all available \ac{GW} observations.
The individual-event posterior \acp{PDF} for $\pi$ are shown in Fig.~\ref{fig:IndPi}, visualised through \acp{KDE}. 
All \ac{GW} events support the region near $\ensuremath{\pi^\mathrm{T}}$, with the strongest constraints coming from the two \acp{BNS} and the lowest-mass \acp{BBH} (GW151226 and GW170806).
This agrees with prior expectations as lower-mass \ac{CBC} signals are dominated by the binary inspiral, described by the \ac{PN}-series, in turn constrained by this analysis.
Apart from a general broadening of the recovered posterior \acp{PDF} in other source parameters, consistent with the addition of a new degree of freedom, we note no general degeneracies between $\pi$ and other parameters.
This is especially noticeable as $\varphi_3$ contains the leading-order terms for both $\pi$ and the effects from \ac{CO}-spins.
As both $\pi$ and spin-parameters however appear jointly at higher \ac{PN}-orders, with different interdependences than in $\varphi_3$, the potentially strong degeneracy is thus broken in this analysis.

\begin{figure
\centering
\includegraphics[width=\columnwidth]{WhatIsPi_Individual.pdf}
\caption{Individual posterior \acp{PDF} of $\pi$ for the thirteen \ac{LVC} \ac{GW} observations~\cite{LIGOScientific:2018mvr, Abbott:2020uma,LIGOScientific:2020stg}.
The dashed line indicates $\ensuremath{\pi^\mathrm{T}}$, the true value of $\pi$.
The shaded gray region indicates the prior \ac{PDF}.
}
\label{fig:IndPi}
\end{figure}

The chronological progression of the joint posterior \ac{PDF} of $\pi$, from the population of \acp{CBC} reported by the \ac{LVC} is shown in Fig.~\ref{fig:JoinPi}, again highlighting the significant contribution of the four lowest-mass events.
Together, these 13 events give a maximum a posteriori value, with associated $90\%$ \ac{CI}, of $\pi=$\mypiMaxPthirteen\ensuremath{^{+0.049}_{-0.091}}\xspace.
For this set of events, the \acl{BF} in favour of \ac{GR} being an accurate description of strong-field gravity is \ensuremath{301}\xspace.

\begin{figure
\centering
\includegraphics[width=\columnwidth]{WhatIsPi_Joint.pdf}
\caption{Joint posterior \acp{PDF} of $\pi$ for the thirteen \ac{LVC} \ac{GW} observations~\cite{LIGOScientific:2018mvr, Abbott:2020uma,LIGOScientific:2020stg}.
The observations are added chronologically, corrresponding to the light to dark blue transition.
The dashed line indicates $\ensuremath{\pi^\mathrm{T}}$, the true value of $\pi$.
\textit{Inset}: The joint posterior \acp{PDF} for the eleven \acp{BBH} (green), the two \acp{BNS} (orange) and all 13 \ac{LVC} \acp{CBC}.
The shaded blue region in the inset corresponds to the $90\%$ \ac{CI}.
}
\label{fig:JoinPi}
\end{figure}

In addition to the eleven \acp{BBH} reported by the \ac{LVC} in~\cite{LIGOScientific:2018mvr}, independent analyses (\citet{Zackay:2019tzo,Venumadhav:2019lyq,Nitz:2019hdf,Zackay:2019btq} hereafter collectively labelled ZVNZ) have claimed an additional nine \ac{BBH} observations, whose posterior \acp{PDF} of $\pi$ are shown in Fig.~\ref{fig:ExtraPi}.
It should be noted that out of all 22 included \acp{CBC}, only GW151216~\cite{Zackay:2019tzo} and GW170304~\cite{Venumadhav:2019lyq} recover $\pi$ disfavouring $\ensuremath{\pi^\mathrm{T}}$ with single-event \acp{BF} in support of \ac{GR} of $1/2$ and $1/5$ respectively.
These two events have previously been identified as especially sensitive to overall prior and analysis choices~\cite{Galaudage:2019jdx,Huang:2020ysn}.
The population of 22 \ac{CBC} observations gives a measurement of $\pi=$\mypiMaxPtwentyTwo\ensuremath{^{+0.048}_{-0.088}}\xspace, and a \ac{BF} in favour of \ac{GR} of \ensuremath{321}\xspace.
This constitutes the strongest constraints on the validity of the positive \ac{PN}-order coefficients to date~\cite{MinusTwoPN}, with a fractional width of the joint $\pi$ $90\%$ \ac{CI} $< 0.04$, more than a factor of 2 improvement over previous single-\ac{PN}-order variability results~\cite{Abbott:2018lct,LIGOScientific:2019fpa,GWTC-1_TGR_dataRelease}.
We also note a more significant improvement when comparing against the constraints on $\varphi_3$, the lowest \ac{PN}-order directly probed by this analysis.
This can be attributed to a combination of the inherent multi-order nature of the variable $\pi$ analysis and the inclusion of a larger population of \ac{CBC} observations than previous studies, thus together enabling a stronger constraint on the validity of the tested theory.

\begin{figure
\centering
\includegraphics[width=\columnwidth]{WhatIsPi_Extra.pdf}
\caption{Individual posterior \acp{PDF} of $\pi$ for the nine ZVNZ \acp{BBH} observations~\cite{Zackay:2019tzo,Venumadhav:2019lyq,Nitz:2019hdf,Zackay:2019btq}.
The dashed line indicates $\ensuremath{\pi^\mathrm{T}}$, the true value of $\pi$.
The shaded gray region indicates the prior \ac{PDF}.
\textit{Inset}: The joint posterior \acp{PDF} for the nine ZVNZ \acp{BBH} (purple), the eleven \ac{LVC} \acp{BBH} (green, cf. Fig.~\ref{fig:JoinPi}), all 20 \acp{BBH} (brown) and all 22 \acp{CBC} (red).
The shaded red region in the inset corresponds to the $90\%$ \ac{CI}.
}
\label{fig:ExtraPi}
\end{figure}

\section{\ac{BBH}-like noise transients -- $\pi$ estimation}
\label{sec:background}

In order to test the reliability of this analysis against spurious false-positives we analyse a set of background triggers, where sections of real data from LIGO~\cite{GWOSC_O1, Abbott:2019ebz} have been offset in time by longer than the light-travel time between sites.
This time-shifted data is thus guaranteed to not contain any real coincident \ac{GW} events, and primarily represent noise-transients from the LIGO instruments.
We select 10 high-significance BBH background triggers produced by the \textsc{PyCBC} search pipeline~\cite{Nitz:2017svb, Usman:2015kfa, Canton:2014ena, PyCBC_zenodo} and used by~\cite{Isi:2018vst}.
Following the same procedure as above, we recover their individual and joint posterior \acp{PDF} on $\pi$, shown in Fig.~\ref{fig:background}.
When the analysed signals do not correspond to \ac{GR}, as is the case for these noise-transients, there is significant scatter of the recovered posterior \acp{PDF} in $\pi$.
The apparent strong constraint from the joint analysis of the 10 triggers stems primarily from the narrow region hosting the black curve in Fig.~\ref{fig:background} being the only range where $p(\pi|d, H) \neq 0$ for all triggers, and that the joint posterior \ac{PDF} has a unit area.
The corresponding \ac{BF} $\lesssim 1/10^{300}$ for $\pi = \ensuremath{\pi^\mathrm{T}}$, further highlighting the variable $\pi$ analysis correctly identifying the non-\ac{GR} features present in these noise triggers.
We conclude that the analysis appears stable against generating false-positive results when exposed to even a small population of 10 known non-\ac{GR} (background noise) signals.

\begin{figure
\centering
\includegraphics[width=\columnwidth]{WhatIsPi_Background.pdf}
\caption{Individual posterior \acp{PDF} (coloured) of $\pi$ for ten \ac{BBH}-like noise triggers from a background distribution.
The dashed line indicates $\ensuremath{\pi^\mathrm{T}}$, the true value of $\pi$.
The shaded gray region indicates the prior\ac{PDF}.
The black line corresponds to the joint posterior \acp{PDF} for the background triggers.
}
\label{fig:background}
\end{figure}

\section{Presence of massive graviton}
\label{sec:lambdag}

To show that this analysis can reveal the presence of realistic beyond-\ac{GR} effects, we simulate three \ac{BBH} systems with parameters consistent with the three \acp{BBH} detected by the \ac{LVC} during \ac{O1}~\cite{TheLIGOScientific:2016pea, LIGOScientific:2018mvr, GWTC1_PE_release}.
We modify the $\varphi_2$ \ac{PN} coefficients of the simulated signals to mimic a massive graviton with a Compton wavelength $\lambda_G$~\cite{Will:1997bb,Zimmerman:2019wzo}.
We choose $\lambda_G$ to be in the range between $2.48\times10^{13}\:\mathrm{km}$, consistent with the current \ac{GW} observational lower bound from~\cite{LIGOScientific:2019fpa}, and $10^{12}\:\mathrm{km}$, a value already ruled-out observationally.
We also simulate ``pure'' \ac{GR}, with $\lambda_G = \infty$.
It should be noted that $\lambda_G$ enters at a \ac{PN}-order where $\pi$ is not present.

In Fig.~\ref{fig:lambdag} we show the joint posterior \acp{PDF} on $\pi$ from the three \ac{BBH} signals for each value of $\lambda_G$.
Given that $\lambda_G \geq 2.48\times10^{13}\:\mathrm{km}$ is not ruled out by current \ac{GW} observations~\cite{LIGOScientific:2019fpa}, it is not surprising that an analysis using this bound yields a posterior \ac{PDF} in agreement with \ac{GR}.
For $\lambda_G = 10^{12}\:\mathrm{km}$, more than an order of magnitude below the current lower bound, the recovered $\pi$ posterior \ac{PDF} is biased away from $\ensuremath{\pi^\mathrm{T}}$ but only marginally informative over the assumed prior. 
This indicates that a strong beyond-\ac{GR} effect, acting partially orthogonal to the changes to the signal from a varying $\pi$, can be sufficient to saturate the constraining power of this test.
Namely, if no allowed value of $\pi$ is able to sufficiently ``correct'' for the beyond-\ac{GR} modification present in the signal, the variable $\pi$ degree of freedom becomes uninformative.
It is instead the case in between these extremes that is the most illustrative, where a presence of a marginal beyond-\ac{GR} effect induces a clear bias in the recovered $\pi$ and a \ac{BF}$\sim 1/10^{15}$ for $\pi = \ensuremath{\pi^\mathrm{T}}$.
Hence, a detection of $\pi \neq \ensuremath{\pi^\mathrm{T}}$ in a population of real observations can be interpreted as first indication of the presence of beyond-\ac{GR} behaviour, with the variable $\pi$ test being especially powerful from its generalisable and multi-\ac{PN}-order nature.
The identification of $\pi \neq \ensuremath{\pi^\mathrm{T}}$ does itself not guide what beyond-\ac{GR} effect is present.
Such questions can only be answered by performing theory-specific model comparison analyses~\cite{Yunes:2013dva,Yunes:2016jcc, Abbott:2018lct,LIGOScientific:2019fpa} over the population of observations for which $\pi \neq \ensuremath{\pi^\mathrm{T}}$.

\begin{figure
\centering
\includegraphics[width=\columnwidth]{WhatIsPi_LambdaG.pdf}
\caption{Joint distributions of $\pi$ shown for four different instances of $\lambda_G$, for \ac{GW} signals consistent with the three \ac{O1} events}
\label{fig:lambdag}
\end{figure}

\section{Discussion}
\label{sec:conclusions}

In the \acl{PN} formalism of \acl{GR}, the mathematical constant $\pi$ presents a powerful null test of our currently preferred theory of gravity.
With $\pi$ simultaneously probing four \ac{PN}-orders, fundamentally describing the same conceptual quantity in all instances, doing so in a theory-agnostic way that is also generalisable and universal across independent \ac{GW} observations, it provides an unmatched capability for validating our understanding of \ac{GR}. 
Using the current set of 22 \ac{CBC} observations in data from LIGO and Virgo, identified by both the \ac{LVC} and independent researchers~\cite{LIGOScientific:2018mvr,Abbott:2020uma,LIGOScientific:2020stg,Zackay:2019tzo,Venumadhav:2019lyq,Nitz:2019hdf,Zackay:2019btq}, we achieve an astrophysical measurement of $\pi=$\mypiMaxPtwentyTwo\ensuremath{^{+0.048}_{-0.088}}\xspace, consistent with the accepted $\ensuremath{\pi^\mathrm{T}}$ value.
This is the most stringent constraint on the positive-order \ac{PN}-series to date~\cite{MinusTwoPN}, and the first viable multi-\ac{PN}-order constraint from \ac{GW} observations.
The analysis also allows, through the construction of a \acl{BF}, direct validation of \ac{GR} with \ac{BF}=\ensuremath{321}\xspace in support of it as the currently favoured theory of gravity.
We have shown the analysis to be robust when exposed to non-signal, but high-significance, \ac{GW} triggers as well as being able to indicate the presence of beyond-\ac{GR} effects in the case where such signals were to exist.
The method presented in this letter is easily extended to future \ac{GW} observations, of both \acp{CBC} and other modelled sources, such as the quasi-monochromatic \acp{GW} emitted by spinning \acp{NS}~\cite{Pisarski:2019vxw, Abbott:2019bed, Authors:2019ztc}, and capable of accommodating observations from across the \ac{GW} spectrum~\cite{Berti:2015itd, Sesana:2016ljz, Vitale:2016rfr, Chamberlain:2017fjl, Carson:2019kkh, Carson:2019rda, Marsat:2020rtl, Toubiana:2020vtf,Kelley:2017vox, Burke-Spolaor:2018bvk, Aggarwal:2018mgp,Becsy:2019dim,Gnocchi:2019jzp}.

\begin{acknowledgments}
The author thanks Katerina Chatziioannou, Riccardo Sturani, Salvatore Vitale and Aaron Zimmerman for helpful suggestions and discussion.
I also thank Maximiliano Isi, and the authors of~\cite{Isi:2018vst}, for providing access to the background BBH-like trigger information.
The author acknowledges support of the National Science Foundation, and the LIGO Laboratory. 
LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the National Science Foundation and operates under cooperative agreement PHY-1764464. 
The author for computational resources provided by the LIGO Laboratory and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459.
This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (\url{https://www.gw-openscience.org}), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. 
LIGO is funded by the U.S. National Science Foundation. 
Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes.
This analysis was made possible by the {\tt LALSuite}~\cite{lalsuite}, {\tt numpy}~\cite{numpy}, {\tt SciPy}~\cite{Virtanen:2019joe} and {\tt matplotlib}~\cite{Hunter:2007ouj} software packages. 
This is LIGO Document Number DCC-P2000159.
\end{acknowledgments}

