# Pneumatic bellows actuated parallel platform control with adjustable
stiffness using a hybrid feed-forward and variable gain I-controller

Martin Varga  Faculty of Mechanical Engineering, Technical University of
Košice, Slovakia  Ivan Virgala  Faculty of Mechanical Engineering, Technical
University of Košice, Slovakia  Michal Kelemen  Faculty of Mechanical
Engineering, Technical University of Košice, Slovakia  Ľubica Miková  Faculty
of Mechanical Engineering, Technical University of Košice, Slovakia  Zdenko
Bobovský  Faculty of Mechanical Engineering, Technical University of Ostrava,
Czech Republic  Peter Ján Sinčák  Faculty of Mechanical Engineering, Technical
University of Košice, Slovakia  Tomáš Merva  Faculty of Mechanical
Engineering, Technical University of Košice, Slovakia

###### Abstract

Redundant cascade manipulators actuated by pneumatic bellows actuators are
passively compliant, rugged and dexterous which are qualities making them
exceptionally well suited for applications in agriculture. Unfortunately
bellows actuators are notoriously difficult to precisely position. This paper
presents a novel control algorithm for the control of a parallel platform
actuated by pneumatic bellows actuators, which is serving as one module of a
cascade manipulator. The algorithm combines a feed-forward controller and a
variable gain I-controller. The feed-forward controller was designed using
experimental data and two regression steps to create a mathematical
representation of the data. The gain of the I-controller depends linearly on
the total reference error, which allows the I-controller to work in concert
with the feed-forward part of the controller. The presented algorithm was
experimentally verified and its performance was compared with two controllers,
an ANFIS controller and a constant gain PID controller, to satisfactory
results. The controller was also tested under dynamic loading conditions
showing promising results.

###### Index Terms:

pneumatic bellows, parallel platform, feed-forward controller, variable PID

††publicationid: pubid:

##  I Introduction

Industrial robots are an indispensable part of the manufacturing process in
many parts of the industry, where their traits, i.e., precision, speed, and
the ability to work basically nonstop, help increase productivity and decrease
cost. It is therefore understandable, that there is a strong incentive to use
industrial robots in other fields, like agriculture and medicine, to name a
few. The most common industrial robots are serial link 6R robots, 2R1T SCARA
robots or parallel Delta robots driven by, most commonly, electric actuators,
or in some cases hydraulic actuators [1]. These industrial robots were
developed for many decades, their design is standardized, and their
mathematical description and control design is fairly well researched.
Unfortunately, these robots lack some key features needed in the
aforementioned new fields of application. For example compliance, agility and
complex modes of motion [2]. These requirements fulfil new emerging classes of
robots i.e., redundant cascade and continuum robots [3]. The redundant robots
are all those that have more degrees of freedom than is necessary to perform a
certain task. Development of these robots took pace from the year 2000.
Redundant robots in general, but especially cascade and continuum robots, have
unique characteristics. For example, according to [4] and [5], compliance, a
good reach to weight ratio, modularity and other. These characteristics arise
from their specific design.

In general, continuum and cascade robots consist of several in series
connected parallel modules that, if underactuated, form a continuum robot and
if fully actuated form a cascade robot. To describe the motion of a redundant
robot and design a control algorithm, it is first necessary to focus on its
individual modules and their properties. The chosen actuator type influences
the achievable properties of a module. Electrical linear servo motors as used
by [6] within the structure of a module or outside as done by [7] provide high
force, stiffness and are usually equipped with position sensors simplifying
control. Hydraulic actuators as can be seen in [8] provide high forces and can
be precisely positioned, jet are slow. Nonstandard actuators for such
applications can also be used, for example SMA springs as seen in the work of
[9] or dielectric materials as described by [10]. Pneumatic actuators are a
popular class of actuators that are used in these applications. They provide
high power density, relatively low weight, can be easily manufactured to
custom specification, as described in [1] and [11], or bought off the shelve
in a variety of types and sizes [12], [13], [14] and the compressibility of
air give them a natural level of compliance. This ability makes them the
actuator of choice for medical applications, like rehabilitation equipment
[15], in flexible endoscopes [16], [17], in agriculture [18], as parts of
mobile robots [19], [20], as the actuator for a high precision positioning
system [21] or a stiffness regulating element in hybrid actuation schemes for
continuum tendon driven robots as presented by paper [22].

Control of a parallel platform module actuated by three or four pneumatic
actuators is a challenging task. One approach to control of these modules is
to use a feed-forward controller alone or in combination with other types.
[23] presented a modelling framework to design a model from which a feed-
forward controller with satisfactory performance was developed. The work [24]
developed a custom bellows type actuator and applied it in a parallel platform
with three degrees of freedom. The control algorithm was a feed-forward
controller based on experimental mapping between pose, external forces and
input pressure. In [25] authors used a feed-forward controller based on a
mathematical model in combination with a variable P gain PI controller to
combat the effect of hysteresis in the system. To positioning of a soft robot,
the paper [26] uses a model based controller using both feed-forward and
feedback components with a structure similar to a PD controller. A planar
platform actuated by pneumatic muscles is controlled by [27] using three fuzzy
controllers synchronized through an ANFIS (Adaptive neuro fuzzy inference
system) based controller. The paper [14] applied a simple constant gain PID
controller for positioning of a parallel platform actuated by four pneumatic
muscles. [28] have demonstrated a nonlinear SMC (sliding mode controller)
based on a PID type sliding surface combined with a lumped element model-based
controllers to control a soft pneumatically actuated robot. In the paper [29]
authors uses a fourier series-based adaptive sliding-mode controller with
$H\infty$ tracking performance to solve the high non-linearity and time-
varying problem for a parallel platform actuated by rod-less pneumatic
cylinders. While [30] dealing with a tendon driven redundant manipulator
proposes a population-based model-free control method that could be applied to
pneumatically actuated manipulators as well.

Based on previous works we have set up to develop the mechanical design and
control system for a rugged redundant cascade manipulator driven by pneumatic
bellows intended for both research and agricultural use. The most difficult
part of the development is the controller design for one separate module, as
noted in previously mentioned articles, this task is notoriously difficult due
to the inherent nonlinear hysteresis behaviour of the chosen actuator type and
MIMO system as a whole. A similar route to previous research in pneumatic
parallel platform module control was taken by relying on experimental data,
but this concept was expanded upon by applying two regression steps to the
data to get a mathematical module model, which is taking the place of a feed-
forward controller. This feed-forward controller was later supplemented by a
variable gain I-controller that facilitates disturbance rejection. The
controller design allows for a on demand change in stiffness of the system
during operation and lends itself well to be a part of complete control system
for the whole cascade manipulator.
Based on previous papers survey, the novelty of the paper can be defined as
follows:

  * •

Development of a novel hybrid FFvI controller

  * •

Establishment of controller design methodology for cascade redundant robots

  * •

Experimental positioning analyses under dynamic disturbance effects

This paper focuses on the development of a novel controller for a 2 DOF
pneumatic parallel platform that represents one module of the pneumatic
manipulator PneuTrunk (see Fig. 1), developed by our Cognitics Lab. The first
part of the paper presents the design and kinematic model of one module of
PneuTrunk. In the second part a feed-forward controller based on an
experimentally identified system with stiffness regulation capabilities
combined with a variable gain I controller for disturbance rejection is
presented. In the third part, the proposed controller is compared with a
simple PID controller and ANFIS controller.

Figure 1: CAD model of PneuTrunk

##  II Design of PneuTrunk module

As can be seen in Fig. 1, the redundant manipulator PneuTrunk is a cascade
type manipulator constructed out of parallel platform modules ordered in
series. The number of modules depends on the required degrees of freedom. One
module, shown in Fig. 2, consists of two duraluminium plates connected by an
universal joint and three evenly spaced pneumatic bellows. The tilt angles
between the top and the bottom plate are measured by two potentiometric
rotation sensors placed in such a way that the axis of the universal joint are
colinear with both axes of the sensors. The pneumatic actuators are of the
shelve Dunlop 2 3/4 x 3 bellows. The pressure in the pneumatic bellows is
controlled by three separate electropneumatic converters SMC ITV1050-31F20.
All tubing is of inner diameter 6 mm to eliminate the effects of tubing
diameter on the dynamic behavior of the bellows. The module is controlled by a
B&R PLC type 4PPC70-0702-20B.

Figure 2: One module of the manipulator PneuTrunk

The maximum operating pressure for one bellow is 7 bar, but to prevent damage
to the system and especially the universal joint, the allowable pressure range
is set to 0-5 bar. The tilt around the $x$-axis and $y$-axis at this pressure
range is written in Tab. I. The $x$-axis is oriented towards the center of one
bellows.

TABLE I: Tilt extremes | tilt | tilt
---|---|---
axis | min | max
x | -24.4 | 16.6
y | -21.2 | 20

It is expected that the module will be driven only by positive pressure. This
has important implications when designing the control algorithm for such a
device. While extension of one bellows is facilitated by simply supplying
pressure, compression is achieved by applying external forces coming
predominantly from the extension of one, or both remaining bellows. This fact
is also the reason why the minimum number of pneumatic bellows is three.
Coincidentally, because the module only has two degrees of freedom, this
design is inherently overactuated. This causes one posture of the module to be
reachable by an infinite number of bellows input pressure combinations and, in
theory, giving the system the ability to change its stiffness without changing
the posture. This was taken into account when designing the control algorithm
for one module.

Figure 3: Flow of information and energy

The flow of information and energy is visualized in Fig. 3. The whole system
is simple and contains only necessary components. It could be argued that
adding a center closed 2/2 valve between each bellows and its corresponding
electropneumatic converter could give the system the ability to pneumatically
lock the bellows extension, improving the systems positioning performance.
Unfortunately, this would also complicate the system and its regulation,
introducing other challenges and distracting from the aim of this paper.

##  III Mathematical modeling

###  III-A Mathematical model of one module

An important step before attempting to design a controller for a module is,
first to create an inverse kinematic model of the module. In other words,
finding a way to map the desired output parameters, here the tilting angles,
to the input parameters, in this case the extension/contraction of the
bellows, see Fig. 4. Inspiration is taken from the work of [31], where bellows
type actuators are represented by two elements connected by a translational
joints and connected to the bottom and top plate by universal joints. This
approach greatly simplifies kinematic modelling. For the purpose of modelling
the dynamics of an actuator, the model needs to be augmented by adding torque
on both universal joints that represents resistance to bending of the
actuator.

Figure 4: Schematic representation of the pneumatic module

For our design, there exists a closed form solution to inverse kinematics in
the form of

$l_{i}=\mid\textbf{a}_{i}-\textbf{H}_{t}\textbf{b}_{i(\alpha_{x},\alpha_{y})}\mid$
(1)

$\textbf{H}_{t}=\textbf{T}_{z34}\times\textbf{R}_{x23}\times\textbf{R}_{y12}\times\textbf{T}_{z01}$
(2)

where $\emph{i}\in\\{1,2,3\\}$ denotes the bellows,
$\textbf{a}_{i}\in\mathbb{R}^{3}$ is coordinates of the center of the bellows
on the bottom plate, $\textbf{b}_{i}\in\mathbb{R}^{3}$ is coordinates of the
center of the bellows on the top plate, $\emph{l}_{i}$ is distance between
point $\textbf{a}_{i}$ and $\textbf{b}_{i}$. Matrix
$\textbf{H}_{t}\in\mathbb{R}^{4\times 4}$ represents the transformation matrix
between fixed coordinate frame $x_{a}y_{a}z_{a}$ and top plate coordinate
frame $x_{b}y_{b}z_{b}$.
$\textbf{T}_{z01},\textbf{T}_{z34},\textbf{R}_{x23},\textbf{R}_{y12}\in\mathbb{R}^{4\times
4}$ where $\textbf{T}_{z01}$ is the translation matrix between the base frame
and a parallel but offset frame $x^{\prime}y^{\prime}z^{\prime}$,
$\textbf{R}_{y12}$ is the rotation matrix rotating frame
$x^{\prime}y^{\prime}z^{\prime}$ around its y axis by $\alpha_{y}$ into
$x^{\prime\prime}y^{\prime\prime}z^{\prime\prime}$, $\textbf{R}_{x23}$ is the
rotation matrix rotating frame
$x^{\prime\prime}y^{\prime\prime}z^{\prime\prime}$ around its x axis by
$\alpha_{x}$ into
$x^{\prime\prime\prime}y^{\prime\prime\prime}z^{\prime\prime\prime}$ and
$\textbf{T}_{z34}$ is the translation matrix between the frame
$x^{\prime\prime\prime}y^{\prime\prime\prime}z^{\prime\prime\prime}$ and a
parallel but offset top plate frame $x_{b}y_{b}z_{b}$, see Fig. 5. Angle
$\alpha_{x}$ is tilt angle of the top plate around axis _x_ and $\alpha_{y}$
is tilt angle of top plate around axis _y_.

Figure 5: Transformation steps between frame $x_{a}y_{a}z_{a}$ and
$x_{b}y_{b}z_{b}$

The role of the established kinematic model in relation to the tilt control of
the module described in later chapters is a central one. The end goal of the
control is to control the tilt, but this is achieved indirectly by controlling
the extension and total pressure in the respective bellows. The presented
inverse kinematic model converts the reference tilt and actual tilt sensed by
rotation sensors to the required extension/contraction and actual deformation
of the bellows.

###  III-B Model of pneumatic bellows

A pneumatic bellows is a linear pneumatic actuator consisting of a bellows
type body and mounting flanges whose free length is dependent on the
difference between the ambient pressure and the pressure inside the bellows.
From a physical point of view, the pneumatic bellows is a pneumatic spring
with variable equilibrium length, dependent on the passive properties of the
bellows and the internal pressure within the bellows Fig. 6.

Figure 6: Bellow length in dependence on internal pressure

The bellows can be modelled using the standard mass, spring, damper model
represented by eq. 3, see Fig. 7. The dominant effect on the system has spring
force and it in turn is dependent on the pneumatic spring stiffness and the
equilibrium height of the bellows at the given internal pressure. The eq. 4
shows this relationship. The bellow without being pressurized behaves like a
spring whose stiffness depends on the shape of the bellow, the current
material properties that are also dependent on other factors like ambient
temperature. Therefore, if the bellows is deformed a spring force appears in
the direction of free height. The equilibrium length is the length of the
bellows at which the deformation force from the internal pressure is at
equilibrium with the spring force. It can be seen that, the equilibrium height
is a nonlinear parameter that depends on multiple other coupled parameters.
Therefore, instead of a physical modelling approach, the model for the
equilibrium height was derived from experimental data by measuring the
equilibrium height at different internal pressures, see Fig. 8. The data was
then interpolated by a third order polynomial function resulting in eq. 5.

Figure 7: Mass-spring-damper model of pneumatic bellow

$F_{M}+F_{b}+F_{k}(k_{p},P_{b})+F_{o}=0$ (3)
$F_{k}=k_{p}\Bigl{(}z_{e}\bigl{(}F_{k},P_{b}\bigr{)}-z_{r}\Bigr{)}$ (4)
$z_{e}=0.45P_{b}^{3}+5.6P_{b}^{2}+23P_{b}+1200$ (5)

where $F_{M}$ is inertial force, $F_{b}$ is damping force, $F_{k}$ is
pneumatic spring force, $F_{o}$ is outside disturbance force, $k_{p}$ is
pneumatic spring stiffness, $P_{b}$ is internal pressure, $z_{e}$ is
equilibrium length, $z_{r}$ is actual length and $F_{m}$ is material spring
force.

Figure 8: Relationship between internal pressure and Equilibrium height

To create a simulation model of a pneumatic spring the stiffness of the spring
is necessary. This parameter can be derived from the eq. 6, 7 and 8.

$k_{p}=\frac{dF_{k}}{dz}$ (6)

where

$F_{k}=(P_{0}-P_{A})A$ (7)

Assuming that the change in bellows internal volume is polytropic we get

$P_{0}V^{2}=constant$ (8)

Combining the above equations

$k_{p}=\frac{P_{0}nA^{2}}{V}+P_{B}\frac{dA}{dz}$ (9)

where $V$ is total volume of air within the bellows and corresponding
pneumatic tube,$P_{0}$ is the absolute pressure inside the bellows,$A$ is an
effective surface of the bellow, $n$ is polytropic constant (for this process
$n=1$). According to eq. 3 \- 9 a simulation model in MATLAB was developed.

Figure 9: Experiment vs. simulation (1st and 3rd row) and the difference
between simulation and experiment (2nd and 4th row)

The results of this model were compared with experimental data, where the
bellow was pressurized to different pressures, a positive extension force was
applied to the bellow and the total extension was measured. The results are in
Fig. 9.

The model of pneumatic bellow gives satisfactory results. The maximum
deviation for the pressures 1 bar to the pressure 6 bar does not exceed 1 mm
while for pressure 0 bar the deviation is nearly 4 mm, which points to either
to a measurement error or to some unknown effect that is much less pronounced
in higher pressures.

This model represents the static behavior of an air bellow performing linear
deformation. It does not capture its bending behavior or its dynamics. To be
able to design a controller for one module of the manipulator PneuTrunk, it is
necessary to also have a basic understanding of the dynamic behavior of one
bellow. This can be seen in the step response of one bellow to an input
pressure step of 5 bar, shown in Fig. 10.

Figure 10: Step response of one bellows to step input 5 bar

There is no discernible overshoot and the rise time from 0 s to maximum value
is about 0.4 s. This means that the system is overdamped and 0.4 s represents
the maximum possible regulation speed. It is also important to note the
behavior of the bellows when going back from 5 bar to 0 bar, where the
actuator is passively returning to the original length. Here, not even after 8
s does the actuator reach the original length.

One important property of a pneumatic bellow, as noted by [32], is its
hysteresis behavior, where inflating and deflating a bellow results in a
different free length at zero internal pressure. In our experiments, this
behavior resulted in a deviation of ± 2mm. To combat this effect, the bellow
was forced by an external stop to always be extended at zero internal pressure
securing a stable free length. Creating a comprehensive bellow model falls
outside the scope of this paper and will be a topic of further research. Non
the less, it gives important insights into the behavior of one bellow
regarding controller development.

##  IV Controller design

Controlling the posture of one module requires the combined effort of all
three of its bellows actuators. The presented controller is designed to deal
both with the non-linearity of the actuators and the over-actuation of the
system. We define the controller consisting of two parts, a feed-forward
controller and a variable gain I-controller (FFvI).

The feed-forward control is widely used in research for these applications,
for example [24] and [23]. It uses an inverse model of a controlled system
without a feedback loop. For this application, it will provide the rough
estimate input. This leverages the lack of overshoot of the actuators even at
large input pressure steps, as can be seen in Fig. 10, it maximizes the
controller speed, and it is generally easy to design and implement. On the
other hand, because of its lack of a feedback loop, as seen in [23], it is
unable to compensate for disturbance forces and system-model deviation. These
are the reasons why the feed-forward controller is supplemented by a variable
gain I-controller designed to complement the feed-forward controller and
dynamically react to any differences between the reference values and actual
values of the controlled variables.

###  IV-A Feed-forward controller design

The feed-forward controller developed in this paper was designed using
experimental data. Various pressure combinations were supplied to each bellows
and the resulting tilt was measured. The supplied pressures ranged from 0 bar
to 5 bar with a 0.2 bar increment. This results in 18275 different pressure
combinations and their corresponding tilting angles. The module workspace can
be seen on Fig. 11.

Figure 11: Module workspace

The pointcloud matrix structure is organized as seen in eq. 10

$\textbf{P}_{PC}=[\alpha_{x},\alpha_{y},P_{1C},P_{2C},P_{3C}]$ (10)

where $\alpha_{x}$ and $\alpha_{y}$ are the measured stable tilt angles which
are the result of corresponding input pressures for the respective bellows
$P_{1C}$, $P_{2C}$ and $P_{3C}$. Because of over actuation and the parallel
nature of the module design, one orientation of the module is achievable by an
infinite combination of input pressures. This can be seen on Fig. 12. Here the
$x$-axis and $y$-axis are the tilt around the respective axis in degrees and
the $z$-axis is the aggregate pressure, which is the sum of all bellows input
pressures in bars. A higher aggregate pressure corresponds to a higher
mechanical stiffness of the system. The control algorithm needs to take this
into account.

Figure 12: Aggregate pressure augmented workspace

The experimentally measured data represent points in aggregate pressure
augmented workspace. To find the inlet pressures from the measured point-
cloud, Algorithm 1 was applied.

Algorithm 1 Extract required bellows input pressures from point cloud

Input: Reference tilt angles $\alpha_{xref}$, $\alpha_{yref}$, required
aggregate pressure level $P_{agr}$ and
${P}_{PC}=[\alpha_{x},\alpha_{y},P_{1C},P_{2C},P_{3C}]$

Output: Required bellows input pressures $P_{1}$, $P_{2}$ and $P_{3}$ to reach
$\alpha_{xref}$ and $\alpha_{yref}$

while isEmpty(Region) do Add all points to Region that pass the criterion:
$\sqrt{(\alpha_{xcloud}-\alpha_{xref})^{2}}\leq$ anT &
$|P_{1cloud}+P_{2cloud}+P_{3cloud}-P_{agr}|\leq$ aggrT

if isEmpty(Region) then

$anT=anT+$ incrementAngle

$aggrT=aggrT+$ incrementPressure

end if

end whileFind point $\textbf{Q}\in Region$ with minimal expression

$\sqrt{(\alpha_{xregion}-\alpha_{xref})^{2}-(\alpha_{yregion}-\alpha_{yref})^{2}}$

if
$\sqrt{(\alpha_{xQ}-\alpha_{xref})^{2}-(\alpha_{yQ}-\alpha_{yref})^{2}}\leq$
acceptableTol then $[P_{1},P_{2},P_{3}]=[P_{1Q},P_{2Q},P_{3Q}]$

else$[P_{1},P_{2},P_{3}]=[mean(P_{1region}),mean(P_{2region}),mean(P_{regionQ})]$

end if

$Region$\- matrix of measured input pressures and corresponding tilts that
will be used to calculate the ended input pressures to reach desired tilt;
$anT$\- maximum Euclidean distance of a measured point from the reference
point in the augmented workspace in the $\alpha_{x}$ $\alpha_{y}$ plane to be
eligible for inclusion in $Region$; $aggrT$\- maximum Euclidean distance of a
measured point from the reference point in the augmented workspace along the
aggregate pressure axis to be eligible for inclusion in $Region$;
$incrementAngle$\- increment to expand $anT$ in case the previous search
yealdet empty $Region$; $incrementPressure$\- increment to expand $aggrT$ in
case the previous search yealdet empty $Region$

Alg. 1 will already supply a set of usable input pressures. Unfortunately, the
results are influenced by errors in measurement and effects of hysteresis. In
a smooth trajectory tracking task this can produce erratic, non-smooth input
pressures. To solve these issues, the above Alg. 1 was supplied with a set of
reference angles ranging from -10° to 10° with an increment of 0.05° for both
$\alpha_{x}$ and $\alpha_{y}$ and a constant aggregate pressure $P_{agr}$. The
result are three 3D meshes representing the relationship between the reference
angles and the three input pressures separately. These meshes were then
separately interpolated as a surface using a second order $x$ and second order
$y$ surface plot. The result for input pressure 1 can be seen in Fig. 13 and
the equation describing this surface is eq. 11

$\begin{split}&P_{1}=2.964-0.1113\alpha_{x}+0.000344\alpha_{y}+0.000726\alpha_{x}^{2}+\\\
&0.00407\alpha_{x}\alpha_{y}-0.00123\alpha_{y}^{2}\end{split}$ (11)

Figure 13: Calculated input pressures and interpolated surface for bellow 1
Figure 14: Surface interpolation for $P_{agr}$ = 3.6:0.6:15 for bellows 1
Figure 15: Output from algorithm 1 vs. Interpolation

This process was repeated for aggregate pressures between 3.6 bar and 15 bar
with increment of 0.6 bar for all three bellows. The result is a set of smooth
surfaces representing the complete augmented workspace, see Fig. 14 for
bellows 1.

The coefficients describing all the surfaces can be further interpolated to
get six equations approximating the complete aggregate pressure augmented
workspace. The coefficients were interpolated by a 7-th order polynomial. Fig.
15 compares the output of Alg. 1 and the interpolated feed-forward controller.
The output of Alg. 1 follows the smooth reference signal, but it shows non-
smooth, erratic step behavior, while the output of the interpolated feed-
forward controller is smooth.

###  IV-B Variable gain I-controller design

The goal of adding a I-part to the controller design is to facilitate
disturbance rejection by integrating the reference error over time and scaling
it by using a gain. In a constant gain I-controller the gain is tuned to and
fixed at a value dependent on the controlled plant. It is simple, easy to
implement and does not necessarily require the plant model for correct design
and unlike a proportional controller it allows for complete error compensation
in a step response. The problem with using a constant gain I-controller is as
follows. The feed-forward controller can immediately supply a rough estimate
for input pressures, but it is never clear before the movement ends how good
this estimation is. Therefore, if the estimate is optimal, a constant gain
I-controller would cause overshoot, requiring the controller to be slow. On
the other hand, if the estimation is sub-optimal, an aggressive constant gain
I-controller is needed to quickly compensate the error. This contradiction in
requirements can be solved by applying a variable gain I-controller, with the
gain dependent on the error, see Eq. 12

$u(t)=K_{i}\Bigl{(}e(t)\Bigr{)}\int_{0}^{t}e(t)\,dt$ (12)

where $K_{i}\Bigl{(}e(t)\Bigr{)}$ is controller gain, $e(t)$ is tilt error,
$t$ is time and $u(t)$ is controller output. The relationship between
I-controller gain and tilt error can be described by different types of smooth
monotonic functions, like linear, exponential etc. For this controller, as a
proof of concept a linear relationship was chosen, see Eq. 13.

$K_{i}=ae_{t}(t)+b$ (13)

$e_{t}(t)=\sqrt{e_{x}^{2}(t)+e_{y}^{2}(t)}$ (14)

where $e_{t}(t)$ is total tilt error, $e_{x}(t)$ is tilt error around $x$-axis
and $e_{y}(t)$ is tilt error around $y$-axis. Parameters $a$ and $b$ were
calculated from experimental data, where $K_{i}=350$ was found to work well
for small total error values below 1° and $K_{i}=75$ was found to not cause
significant overshoot at error values bellow 5°. This relationship is
described by eq. 15 and can be seen on Fig. 16

$K_{i}=-68.75e_{t}+418.75$ (15)

Figure 16: Tilt error in dependence on I-controller gain

The comparison between the performance of a constant gain I-controller and our
variable gain controller in combination with our feed-forward controller is
depicted in Fig. 17. One can see, that in cases a) and b) a gain of 350
results in a significant overshoot, but for small changes in tilt and a large
residual error after feed-forward controller action, like in case c), the
controller is fast and has acceptable overshoot. On the other hand, a gain of
75 has no overshoot in any case but is slow and has best performance if the
residual error from feed-forward controller action is small, like in case b).
The performance of both fixed gain controllers in case c) is not satisfactory.
Our variable gain controller performs satisfactory in all cases.

Figure 17: Comparison between constant and variable gain I-controller under
different conditions

As mentioned above, for our controller, we have chosen a linear relationship
to govern the variable gain. A linear relationship is simple, is easy to tune
and should result in predictable behaviour of the module. Nevertheless, other
types of functions like exponential, polynomial or even logarithmic can also
be applied and could result in better controller performance than with a
linear equation. The purpose of this article is to show the viability of the
by us presented controller and using a linear relationship for this purpose is
sufficient. The complexity of comparing different types of governing equations
and the required tuning methods warrant its own article.

##  V Experimental verification of FFvI controller

This section will focus on experimental comparison between FFvI controller and
controllers designed according to established algorithms. First of all, the
performance of the feed-forward part of the FFvI algorithm will be compared
with a ANFIS controller designed using the same dataset as our feed-forward
controller. In the second part the complete FFvI controller will be compared
to a constant gain PID controller.

The ANFIS controller was designed using the MATLAB neuro-fuzzy designer. Three
controllers for each bellows separately were created. The input data are tilt
angles $\alpha_{x}$ and $\alpha_{y}$ and the output is the corresponding
pressure. The teaching data are picked from the same data, that is used to
design the feed-forward controller, but is limited to having an aggregate
pressure of $9\pm{1.5}$ bar. The feed-forward controller is also set to the
same level of aggregate pressure. This will decrease the teaching time and
ensure more reliable results. The neural network used is a Sugeno type network
[33] with 10 linear generalized bell-shaped membership functions for each
input. The minimum achieved teaching error is in Tab. II.

TABLE II: Minimum teaching error Bellows | Teaching error
---|---
1 | 0.0205
2 | 0.0227
3 | 0.0272

Both controllers were supplied with the sets of reference tilts and their
outputs were compared to input pressures creating the reference tilts. The
resulting comparison between the outputs of the feed-forward controller and
the ANFIS controller are shown in Fig. 18. The minimum error and mean error
are in Tab. III.

Figure 18: Comparison between feed-forward and ANFIS controller TABLE III: Maximum error Bellows | Controller | maximum error | mean error
---|---|---|---
1 | feed-forward | 0.0959 | 0.0217
ANFIS | 0.0231 | -0.0028
2 | feed-forward | 0.1463 | 0.0219
ANFIS | 0.0174 | -0.0011
2 | feed-forward | 0.1286 | 0.025
ANFIS | 0.0294 | 0.0029

The ANFIS controller has a consistently lower error and the mean error is also
lower. Newer the less, as can be seen from the above results, both controllers
perform well, with the maximum error of the feed-forward controller not
exceeding 0.15°. The magnitude of this error is still well within what would
be considerate acceptable for this application. Considering that the feed-
forward controller can approximate the whole pressure augmented workspace
while the ANFIS controller is specifically designed to work in the tested
aggregate pressure region, these results are encouraging.

When assessing the complete FFvI controller, one must keep in mind that the
physical module is meant to be part of a robot. Therefore, we have specified
that the overshoot, when performing a movement, should not exceed 2° on
$\alpha_{x}$ and $\alpha_{y}$ separately. This number, although arbitrary,
should be a good controller benchmark for this case and the controller is
expected to be further tuned after being installed in the complete robot
control system.

As a comparison to our FFvI controller a PID-based controller was used. The
controller consists of three identical constant gain PID controller
controlling all the bellows separately. The PID controllers were tuned using
the P-I-D tuning approach described in [34]. The relevant constants can be
seen in Tab. IV. Although the values of the tuned constants are not optimal,
the authors believe that they are close enough to the optimal values combining
adequate speed while not exceeding our criterion of 2° overshoot on
$\alpha_{x}$ and $\alpha_{y}$.

TABLE IV: PID controller gains | Gain
---|---
$K_{P}$ | 40
$K_{I}$ | 240
$K_{D}$ | 5

The performance of both controllers can be seen on Fig. 19. The reference
signal is sequentially alternating between $\alpha_{xref}=8^{\circ}$,
$\alpha_{yref}=-10^{\circ}$, $\alpha_{xref}=-8^{\circ}$ and
$\alpha_{yref}=10^{\circ}$. These values were chosen because all bellows need
to be engaged simultaneously to different degrees, they represent values in
the middle part of the module workspace and it is expected that most movement
will be in this area and they also demonstrate the asymmetric behavior of the
platform.

Figure 19: Comparison PID vs. FFvI controller - input step
$\alpha_{xref}=-8^{\circ}$ to $8^{\circ}$; $\alpha_{yref}=-10^{\circ}$ to
$10^{\circ}$

The performance of both controllers types is shown in Tab. V. It needs to be
noted that, this comparison deviates from a standard comparison of step
responses by not having the platform at zero tilt and testing a combined
rotation around $\alpha_{x}$ and $\alpha_{y}$. Nevertheless, this comparison
is much closer to comparing the controller performance under more realistic
conditions.

TABLE V: Performance comparison

$\alpha_{x}-8\rightarrow 8$ | FFvI | PID
---|---|---
rise time [s] | 0.42 | 0.4
setling time 5% [s] | 2.03 | 5.42
first maximum [s] | 0.54 | 2.15
overshoot [%] | 21.1 | 10.3

$\alpha_{y}10\rightarrow-10$ | FFvI | PID
---|---|---
rise time [s] | 2.74 | 4.13
setling time 5% [s] | 3.43 | 4.5
first maximum [s] | 0.65 | 7.45
overshoot [%] | 3 | 3

$\alpha_{x}8\rightarrow-8$ | FFvI | PID
---|---|---
rise time [s] | 0.5 | 1.23
setling time 5% [s] | 4.46 | 5¡
first maximum [s] | 0.57 | 3.4
overshoot [%] | 8 | 21.3

$\alpha_{y}-10\rightarrow 10$ | FFvI | PID
---|---|---
rise time [s] | 0.53 | 3
setling time 5% [s] | 2.73 | 8.04
first maximum [s] | 0.64 | 7.5
overshoot [%] | 5 | 9.6

TABLE VI: Sinus reference signal parameters sinus signal | $\alpha_{xref}$ | $\alpha_{yref}$
---|---|---
amplitude [°] | 10 | 10
period [s] | $2\pi$ | $2\pi$
phase [rad] | 0 | $\pi/2$
bias [°] | 0 | 0

It can be seen, that both controllers fulfill the condition of not having more
than $2^{\circ}$ overshoot. From Tab. V it can be seen, that the FFvI
controller is faster in all categories and has more overshoot only on when
going from $-8^{\circ}\rightarrow 8^{\circ}$. The second comparison between
both controllers is in following a sine reference signal with parameters
written in Tab. VI, as can be seen in Fig. 20. Here the FFvI controller can
tightly and smoothly follow the sinus reference signal while overshooting at
the maximum and minimum of the reference signal. The overshoot, again, does
not exceed 2° for either tilt angles. The PID controller, while not
overshooting nearly as much, lags constantly behind the reference signal by
about 0.7 s and the plot is in some parts jittery. The overall tracking error
is much smaller for FFvI than for PID controller. It can be said that both
controllers perform acceptably, while the FFvI controller is faster while
still passing the overshoot criteria.

Figure 20: Comparison PID controller vs. FFvI controller - input sinus wave
Figure 21: Module response to stiffness change

One ability that the FFvI controller has is to change the stiffness of the
system while maintaining reference tilt. This is shown in Fig. 21. The module
is first tilted to $\alpha_{xref}=10^{\circ}$, $\alpha_{yref}=5^{\circ}$ at a
requested aggregate pressure of 3 bar. Then, the aggregate pressure is changed
stepwise 6 bar, 9 bar, 12 bar and 15 bar, respectively. A change in aggregate
pressure corresponds to a change in stiffness of the system. As can be seen
from the Fig. 21 pressure and therefore stiffness can be changed online. This
change results in a momentous destabilization of the system resulting in
slight position loss. The maximum error for our test was at the transition
between aggregate pressure 12 bar and 15 bar with maximum error for
$\alpha_{xref}=1.3^{\circ}$, $\alpha_{yref}=1.65^{\circ}$. This can be again,
attributed to nonsynchronous pressure change between the bellows and possible
measurement errors in the original feed-forward controller input data.
Therefore, instead of a sharp aggregate pressure step a smooth aggregate
pressure transition should decrease this issue.

##  VI Dynamic test of FFvI tilt platform controller

To be able to apply the proposed controller as part of the control system of
the whole manipulator it is necessary to test the controller under dynamic
conditions.

For this purpose a two axis loading mechanism was developed, see Fig. 22. It
consists of two linear motion axis stacked perpendicular on top of each other
and a attached to the top axis. This mechanism is mounted on top of a tilt
module with the axis axis of the tilt platform aligned with the linear motion
axis. It can be seen that the potential loading momentum will be different for
both axis because the bottom axis is loaded not only by the loading weight but
also the top axis itself. The possible motion of one linear axis is 0$mm$ to
330$mm$ and is centered on the central axis of the tilt module. The weight of
the axis and loading weight is in Tab. VII .This mechanism is used to generate
dynamic loading forces to study the capability of the FFvI controller to
reject dynamic disturbances.

Figure 22: Dynamic two axis loading mechanism TABLE VII: Loading weigths | weight [g]
---|---
$linear~{}axis$ | 2467
$loading~{}weigth$ | 1802

The experiment seen in Fig. 24 is done for reference tilts $\alpha_{x}$=0° and
$\alpha_{y}$=0°. Error rejection is facilitated by the I controller part of
the algorithm, hence the test will be performed with the I controller active
and, for comparison, with the I controller inactive. Also, the experiment will
be done for different values of aggregate pressure and for different speeds of
movement of the loading weight. The movement of the load was converted to
loading momentum around axis $x$ and $y$ (see Fig. 23). Fig. 23 and Fig. 24
show only the the extremes of the test, different combinations of load speed
and aggregate pressure were also tested.

Figure 23: Disturbance momentum at different loading speeds

The disturbance momentum was indirectly established from the measured movement
of the load and the geometry of the mechanism. The shown time window is the
same for all movement speeds. The maximum total error for all measurements is
in Tab. VIII.

Figure 24: Comparison between active and inactive I regulator at tilt holding at extreme aggregate pressures and dynamic effects TABLE VIII: Maximum total error under dynamic load velocity; pressure | I on [°] | I off [°]
---|---|---
10$mm/s$; 6$bar$ | 0.93 | 3.57
30$mm/s$; 6$bar$ | 0.66 | 4.29
50$mm/s$; 6$bar$ | 1.52 | 3.48
70$mm/s$; 6$bar$ | 1.51 | 6.18
10$mm/s$; 9$bar$ | 0.31 | 3.84
30$mm/s$; 9$bar$ | 0.62 | 3.75
50$mm/s$; 9$bar$ | 1.02 | 3.66
70$mm/s$; 9$bar$ | 1.05 | 6.82
10$mm/s$; 12$bar$ | 0.28 | 2.78
30$mm/s$; 12$bar$ | 0.60 | 2.58
50$mm/s$; 12$bar$ | 0.88 | 3.17
70$mm/s$; 12$bar$ | 1.15 | 5.42
10$mm/s$; 15$bar$ | 0.52 | 2.74
30$mm/s$; 15$bar$ | 0.63 | 2.68
50$mm/s$; 15$bar$ | 0.87 | 3.24
70$mm/s$; 15$bar$ | 1.14 | 3.44
mean | 0.86 | 3.85

From the above figures one can deduce important findings regarding the
behaviour of the regulator under dynamic changing load. The controller greatly
benefits the precision of the control by rejecting the residual error created
by feed-forward error imperfection where the maximum total error with the I
part active was 1.52° and the mean error was 0.86° and the maximum total error
with the I part inactive was 6.85° and the mean error was 3.85°. Apart from
that a well tuned variable gain I-controller has a significant stabilisation
effect at higher dynamic loads. It can also be seen, especially at higher
speeds of the load, that a higher aggregate pressure has a positive effect on
the rejection of dynamic loads.

It can be concluded that his controller in its current state is able to
control the prototype manipulator, especially at lower speeds and higher
aggregate pressure settings.

##  VII Conclusion

In this paper a new type of controller for the control of the tilt of a
pneumatic bellows actuated module of the cascade robot PneuTrunk was
presented. The controller consists of a feed-forward controller designed using
experimental data and a variable gain I-controller. The feed-forward
controller is created by fitting the data at a certain pressure level using a
polynomial function and subsequently, again fitting the resulting set of
polynomial function constants by another set of polynomials. This allows for a
simple and fast controller, that not only allows to control the tilt of the
module, but also its stiffness on demand. The variable gain I-controller
supplements the feed-forward controller by adding a feedback loop, hence
facilitating disturbance rejection and correcting for feed-forward controller
imperfections. The variable gain allows for fast error correction while
limiting overshoot and windup at the same time. This hybridisation approach
allows for a simple controller design for a complex MIMO systems that can be
easily adjusted and updated.

This controller was compared to other established controllers. On the feed-
forward level, the controller was compared with an ANFIS controller,
delivering comparable results. Comparing the complete controller to a tuned
PID controller showed that our FFvI controller is faster, can reliably follow
a harmonic reference signal while maintaining required performance parameters.
This controller was also testet under dynamic load to satisfactory results.
The maximum FFvI controller error during the positioning of the module with
consideration of dynamic disturbance was only 1.52°.

To further improve the performance of the controller, it is necessary to
create a comprehensive mathematical model of the module, mainly to combat the
detrimental effects of the hysteretic behavior of the bellows. It is also
appropriate to compare different types of functions driving the variable gain
of the I-controller. In the future, this controller will be applied as a part
of a larger control system controlling the pneumatic cascade robot PneuTrunk.

## Acknowledgments

This research was funded by Slovak Grant Agency VEGA 1/0436/22 Research on
modeling methods and control algorithms of kinematically redundant mechanisms
and VEGA 1/0201/21 Mobile mechatronic assistant. This research has been also
elaborated under support of the project Research Centre of Advanced
Mechatronic Systems, reg. No. CZ.02.1.01$/0.0/0.0/16_{0}19/$0000867 in the
frame of the Operational Program Research, Development and Education.

## References

  * [1] G. Gregov, T. Ploh, and E. Kamenar, “Design, development and experimental assessment of a cost-effective bellow pneumatic actuator,” _Actuators_ , vol. 11, no. 6, p. 170, Jun 2022.
  * [2] M. Kelemen, I. Virgala, T. Liptak, L. Mikova, F. Filakovsky, and V. Bulej, “A novel approach for a inverse kinematics solution of a redundant manipulator,” _Applied Sciences_ , vol. 8, no. 11, p. 2229, Nov 2018. [Online]. Available: http://dx.doi.org/10.3390/app8112229
  * [3] I. Virgala, M. Kelemen, M. Varga, and P. Kuryło, “Analyzing, modeling and simulation of humanoid robot hand motion,” _Procedia Engineering_ , vol. 96, pp. 489–499, 2014, modelling of Mechanical and Mechatronic Systems. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1877705814031798
  * [4] G. Chirikjian and J. Burdick, “Hyper-redundant robot mechanisms and their applications,” in _Proceedings IROS ’91:IEEE/RSJ International Workshop on Intelligent Robots and Systems ’91_ , 1991, pp. 185–190 vol.1.
  * [5] ——, “A hyper-redundant manipulator,” _IEEE Robotics & Automation Magazine_, vol. 1, no. 4, pp. 22–29, 1994.
  * [6] Y. Zhao, X. Song, X. Zhang, and X. Lu, “A hyper-redundant elephant’s trunk robot with an open structure: Design, kinematics, control and prototype,” _Chinese Journal of Mechanical Engineering_ , vol. 33, no. 1, p. 96, Dec 2020\.
  * [7] D. Grace, J. Lee-Ortiz, M. Garcia, A. Contreras-Esquen, A. Tekes, and A. A. Amiri Moghadam, “Development of a novel six dof soft parallel robot,” in _SoutheastCon 2022_ , 2022, pp. 81–86.
  * [8] Y. Pi and X. Wang, “Trajectory tracking control of a 6-dof hydraulic parallel robot manipulator with uncertain load disturbances,” _Control Engineering Practice_ , vol. 19, no. 2, pp. 185–193, 2011.
  * [9] C. Cheng, J. Cheng, and W. Huang, “Design and development of a novel sma actuated multi-dof soft robot,” _IEEE Access_ , vol. 7, pp. 75 073–75 080, 2019.
  * [10] A. Wingert, M. D. Lichter, S. Dubowsky, and M. Hafez, “Hyper-redundant robot manipulators actuated by optimized binary-dielectric polymers,” in _SPIE Smart Structures and Materials + Nondestructive Evaluation and Health Monitoring_ , 2002.
  * [11] C. L. Nall and P. A. Bhounsule, “A miniature 3d printed on-off linear pneumatic actuator and its demonstration into a cartoon character of a hopping lamp,” _Actuators_ , vol. 8, no. 4, p. 72, Oct 2019.
  * [12] G. Tao, C. Shang, and D. Meng, “Adaptive robust posture control of a 3-rps pneumatic parallel platform with unknown deadzone,” _Mathematical Problems in Engineering_ , vol. 2016, p. 2034923, Jun 2016.
  * [13] F. Trebua, M. Kelemen, M. Pastor, and I. Virgala, “Trajectory tracking controller of air bellow,” _Journal of Automation and Control_ , vol. 4, no. 2, pp. 51–55, 2016.
  * [14] D. Gryparis, G. Andrikopoulos, and S. Manesis, “Parallel robotic manipulation via pneumatic artificial muscles,” in _2014 11th International Conference on Informatics in Control, Automation and Robotics (ICINCO)_ , vol. 02, 2014, pp. 29–36.
  * [15] Y. Wang and Q. Xu, “Design and testing of a soft parallel robot based on pneumatic artificial muscles for wrist rehabilitation,” _Scientific Reports_ , vol. 11, no. 1, p. 1273, Jan 2021.
  * [16] M. Cianchetti, T. Ranzani, G. Gerboni, I. De Falco, C. Laschi, and A. Menciassi, “Stiff-flop surgical manipulator: Mechanical design and experimental characterization of the single module,” in _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2013, pp. 3576–3581.
  * [17] N. Garbin, L. Wang, J. H. Chandler, K. L. Obstein, N. Simaan, and P. Valdastri, “A disposable continuum endoscope using piston-driven parallel bellow actuator,” in _2018 International Symposium on Medical Robotics (ISMR)_ , 2018, pp. 1–6.
  * [18] G. Z., C. B., L. J., F. X., L. Z., W. T., and W. L., “An opposite-bending-and-extension soft robotic manipulator for delicate grasping in shallow water,” _Front. Robot. AI_ , 2019.
  * [19] Z. Liu and K. Karydis, “Position control and variable-height trajectory tracking of a soft pneumatic legged robot,” _CoRR_ , vol. abs/2108.01338, 2021.
  * [20] G. Granosik and M. Kaczmarski, “Bellows driven, muscle steered caterpillar robot,” in _Climbing and Walking Robots_ , M. O. Tokhi, G. S. Virk, and M. A. Hossain, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2006, pp. 743–750.
  * [21] T. Fujita, K. Sakai, Y. Takagi, ¡br/¿Kenji Kawashima, , and T. Kagawa, “Ultra precise positioning of a stage driven by pneumatic bellows,” _International Journal of Automation Technology_ , vol. 5, no. 4, pp. 508–515, 2011.
  * [22] E. Harsono, J. Yang, S. Bhattacharya, and H. Yu, “Design and analysis of a novel hybrid-driven continuum robot with variable stiffness,” _Mechanism and Machine Theory_ , vol. 177, p. 105067, 2022.
  * [23] X. Huang, X. Zhu, and G. Gu, “Kinematic modeling and characterization of soft parallel robots,” _IEEE Transactions on Robotics_ , pp. 1–15, 2022.
  * [24] J. Zhang, H. Wei, Y. Shan, P. Li, Y. Zhao, L. Qi, and H. Yu, “Modeling and experimental study of a novel multi-dof parallel soft robot,” _IEEE Access_ , vol. 8, pp. 62 932–62 942, 2020.
  * [25] M. Calabria, F. Schreiber, Y. Sklyarenko, D. Inkermann, A. Raatz, T. Vietor, and W. Schumacher, “Redundancy resolution and control of manipulators driven by antagonistic pneumatic muscles,” in _2012 17th International Conference on Methods & Models in Automation & Robotics (MMAR)_, 2012, pp. 119–124.
  * [26] R. K. Katzschmann, C. D. Santina, Y. Toshimitsu, A. Bicchi, and D. Rus, “Dynamic motion control of multi-segment soft robots using piecewise constant curvature matched with an augmented rigid body model,” in _2019 2nd IEEE International Conference on Soft Robotics (RoboSoft)_ , 2019, pp. 454–461.
  * [27] L. D. Khoa, D. Q. Truong, and K. K. Ahn, “Synchronization controller for a 3-r planar parallel pneumatic artificial muscle (pam) robot using modified anfis algorithm,” _Mechatronics_ , vol. 23, no. 4, pp. 462–479, 2013.
  * [28] A. H. Khan and S. Li, “Sliding mode control with pid sliding surface for active vibration damping of pneumatically actuated soft robots,” _IEEE Access_ , vol. 8, pp. 88 793–88 800, 2020.
  * [29] M.-H. Chiang and H.-T. Lin, “Development of a 3d parallel mechanism robot arm with three vertical-axial pneumatic actuators combined with a stereo vision system,” _Sensors_ , vol. 11, no. 12, p. 11476–11494, Dec 2011.
  * [30] H. Mo, X. Li, B. Ouyang, G. Fang, and Y. Jia, “Task autonomy of a flexible endoscopic system for laser-assisted surgery,” _Cyborg and Bionic Systems_ , vol. 2022, p. 9759504, Aug 2022.
  * [31] O. Lakhal, A. Melingui, A. Chibani, C. Escande, and R. Merzouki, “Inverse kinematic modeling of a class of continuum bionic handling arm,” in _2014 IEEE/ASME International Conference on Advanced Intelligent Mechatronics_ , 2014, pp. 1337–1342.
  * [32] M. K. Mishra, A. K. Samantaray, and G. Chakraborty, “Fractional-order bouc-wen hysteresis model for pneumatically actuated continuum manipulator,” _Mechanism and Machine Theory_ , vol. 173, p. 104841, 2022.
  * [33] T. Takagi and M. Sugeno, “Fuzzy identification of systems and its applications to modeling and control,” _IEEE Transactions on Systems, Man, and Cybernetics_ , vol. SMC-15, pp. 116–132, 1985.
  * [34] F. Haugen, _PID Control_. Tapir Academic Press, 2004. [Online]. Available: https://books.google.sk/books?id=8CIlAAAACAAJ


# Gravitational-Wave Phasing of Compact Binary Systems
to the Fourth-and-a-Half post-Newtonian Order

Luc Blanchet<EMAIL_ADDRESS>$\mathcal{G}\mathbb{R}\varepsilon{\mathbb{C}}\mathcal{O}$, Institut
d’Astrophysique de Paris,
UMR 7095, CNRS, Sorbonne Université, 98bis boulevard Arago, 75014 Paris,
France  Guillaume Faye<EMAIL_ADDRESS>$\mathcal{G}\mathbb{R}\varepsilon{\mathbb{C}}\mathcal{O}$, Institut
d’Astrophysique de Paris,
UMR 7095, CNRS, Sorbonne Université, 98bis boulevard Arago, 75014 Paris,
France  Centre for Strings, Gravitation and Cosmology, Department of Physics,
Indian Institute of Technology Madras, Chennai 600036, India  Quentin Henry
<EMAIL_ADDRESS>Max Planck Institute for Gravitational Physics
(Albert Einstein Institute), D-14476 Potsdam, Germany  François Larrouturou
<EMAIL_ADDRESS>Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany  David Trestini<EMAIL_ADDRESS>Laboratoire Univers et Théories, Observatoire de Paris, Université PSL,
Université Paris Cité, CNRS, F-92190 Meudon, France
$\mathcal{G}\mathbb{R}\varepsilon{\mathbb{C}}\mathcal{O}$, Institut
d’Astrophysique de Paris,
UMR 7095, CNRS, Sorbonne Université, 98bis boulevard Arago, 75014 Paris,
France

###### Abstract

The inspiral phase of gravitational waves emitted by spinless compact binary
systems is derived through the fourth-and-a-half post-Newtonian (4.5PN) order
beyond quadrupole radiation, and the leading amplitude mode
$(\ell,\mathrm{m})=(2,2)$ is obtained at 4PN order. We also provide the
radiated flux, as well as the phase in the stationary phase approximation.
Rough numerical estimates for the contribution of each PN order are provided
for typical systems observed by current and future gravitational wave
detectors.

###### pacs:

04.25.Nx, 04.30.-w, 97.60.Jd, 97.60.Lf

††preprint: DESY-23-043

At the time when the LIGO and Virgo gravitational-wave detectors were
approved, there was no theoretical prediction available for gravitational
waves (GWs) generated by compact binary systems, apart from that of the famous
Einstein quadrupole formula [1, 2, 3]. However, it was soon realized that,
given the frequency band and the expected sensitivity of these ground-based
detectors, the waveform modeling was to be drastically improved in order to
extract all the potential information from the signal, at least in the case of
the inspiral of two neutron stars [4, 5]. The breakthrough came with the
merging of the post-Newtonian (PN) and the multipolar post-Minkowskian (MPM)
expansions into a single formalism [6, 7, 8, 9, 10], that was applied with
success to derive step by step the waveform of compact binary systems up to
3.5PN order [11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. Since then, many works
(outlined below) have aimed at extending the precision of this result to the
next level, namely 4PN or even 4.5PN order beyond the Einstein quadrupole
formula.

This Letter provides the final results of these efforts, _i.e._ , the GW
phasing of non-spinning compact binary systems on quasi-circular orbits up to
4.5PN order, as well as the dominant GW mode, given by
$(\ell,\mathrm{m})=(2,2)$, at 4PN order. Ready to be used for building
accurate PN template banks for the detection and analysis of the inspiral
phase of compact binaries, they should be important for third generation
ground-based detectors (Einstein Telescope and Cosmic Explorer), future space-
borne detectors (LISA and TianQin) and of course the current second-generation
detectors (LIGO, Virgo and KAGRA). All results presented in this Letter are to
be found in the ancillary file [21] associated with the companion paper [22].

Besides improving the detectors’ data analysis, the motivation for computing
high PN orders is also to perform high-accuracy tests of general relativity
(GR), since the PN coefficients directly probe the non-linear structure of the
theory. By confronting results from the PN expansion against data, one can put
constraints on potential deviations from GR [23, 24]. This has already allowed
for the confirmation of the signature of GW tails [25, 26], and is promising
for tests with future multi-band detections between LISA and ground-based
detectors [27].

We denote by $f(t)$ the frequency of the dominant $(2,2)$ mode of the GW as
measured by an observer in the asymptotically flat region far from the source
(recall that this is twice the orbital frequency), and by
$\psi(t)=\pi\int\\!\mathrm{d}tf(t)$ the corresponding phase. As usual, it is
convenient to define the directly-measurable PN parameter
$x=\mathcal{O}(c^{-2})$ by

$x\equiv\left(\frac{\pi Gmf}{c^{3}}\right)^{2/3}\,,$ (1)

where $m=m_{1}+m_{2}$ is the binary’s total mass, $m_{1}$ and $m_{2}$ being
the constant masses of the progenitors. For circular orbits, $x$ may be
defined invariantly from the Killing vector of the helical symmetry in the
asymptotically flat space-time. Since compact binaries tend to have
circularized by the time they enter the detector’s frequency band [28], we
only consider the case of quasi-circular orbits, for which the time-evolution
of the frequency and phase (or “chirp”) is entirely driven by the energy flux-
balance equation,

$\frac{\mathrm{d}E}{\mathrm{d}t}=-\mathcal{F}\,,$ (2)

where $E$ denotes the invariant energy of the compact binary and $\mathcal{F}$
the total energy flux (or GW luminosity). Both $E$ and $\mathcal{F}$ in the
balance equation are unique functions of the PN parameter $x$ and the two
masses. They have to be evaluated with the same relative PN precision, in the
present case 4.5PN ($\sim x^{9/2}$). From Eq. (2), we derive a simple ordinary
differential equation for the frequency as a function of time, and, once it is
solved, a further integration yields the phase as a function of frequency.

The invariant energy $E$ follows from the conservative dynamics of the compact
binary at 4PN order, which have been obtained by various groups using
different methods: (i) the Arnowitt-Deser-Misner (ADM) Hamiltonian formalism
[29, 30, 31, 32] led to complete results except for one “ambiguity” parameter,
which was fixed by resorting to a comparison with gravitational self-force
(GSF) results [33]; (ii) the Fokker Lagrangian formalism in harmonic
coordinates [34, 35, 36, 37] yielded for the first time a complete result
without any ambiguity parameter; (iii) the effective field theory (EFT)
approach [38, 39, 40, 41, 42, 43, 44, 45] also led to the complete and
unambiguous result. From this series of works, the binary’s invariant energy
was obtained as the Noetherian quantity associated with temporal dilatation,
and reads at 4PN order

$\displaystyle E$ $\displaystyle=-\frac{m\nu
c^{2}x}{2}\Biggl{\\{}1+\biggl{(}-\frac{3}{4}-\frac{\nu}{12}\biggr{)}x+\biggl{(}-\frac{27}{8}+\frac{19}{8}\nu-\frac{\nu^{2}}{24}\biggr{)}x^{2}$
$\displaystyle\qquad\quad+\biggl{[}-\frac{675}{64}+\biggl{(}\frac{34445}{576}-\frac{205}{96}\pi^{2}\biggr{)}\nu-\frac{155}{96}\nu^{2}-\frac{35}{5184}\nu^{3}\biggr{]}x^{3}$
$\displaystyle\qquad\quad+\Biggl{[}-\frac{3969}{128}+\biggl{(}-\frac{123671}{5760}+\frac{9037}{1536}\pi^{2}+\frac{896}{15}\gamma_{\text{E}}+\frac{448}{15}\ln(16x)\biggr{)}\nu$
$\displaystyle\qquad\quad\qquad+\biggl{(}-\frac{498449}{3456}+\frac{3157}{576}\pi^{2}\biggr{)}\nu^{2}+\frac{301}{1728}\nu^{3}+\frac{77}{31104}\nu^{4}\Biggr{]}x^{4}+\mathcal{O}\bigl{(}x^{5}\bigr{)}\Biggr{\\}}\,.$
(3)

We denote by $\nu\equiv m_{1}m_{2}/m^{2}$ the symmetric mass ratio
($\gamma_{\text{E}}$ is the Euler constant). Since there are no terms of half-
integer PN order for circular orbits, this expression is actually valid up to
4.5PN order.

The second input is the energy flux, which we have computed using the PN-MPM
formalism applied to compact binaries at 4.5PN beyond the leading quadrupole
formula. Crucial to this computation was the recently-completed source mass
quadrupole moment at 4PN order [46, 47, 48, 49], the source current quadrupole
moment at 3PN order [50] and the non-linear tail-of-memory effect [51, 52]. We
provide the technical details of the derivation in the companion paper [22],
and report here only the final result:

$\displaystyle\mathcal{F}=\frac{32c^{5}}{5G}\nu^{2}x^{5}\Biggl{\\{}$
$\displaystyle 1+\biggl{(}-\frac{1247}{336}-\frac{35}{12}\nu\biggr{)}x+4\pi
x^{3/2}+\biggl{(}-\frac{44711}{9072}+\frac{9271}{504}\nu+\frac{65}{18}\nu^{2}\biggr{)}x^{2}+\biggl{(}-\frac{8191}{672}-\frac{583}{24}\nu\biggr{)}\pi
x^{5/2}$
$\displaystyle+\Biggl{[}\frac{6643739519}{69854400}+\frac{16}{3}\pi^{2}-\frac{1712}{105}\gamma_{\text{E}}-\frac{856}{105}\ln(16\,x)+\biggl{(}-\frac{134543}{7776}+\frac{41}{48}\pi^{2}\biggr{)}\nu-\frac{94403}{3024}\nu^{2}-\frac{775}{324}\nu^{3}\Biggr{]}x^{3}$
$\displaystyle+\biggl{(}-\frac{16285}{504}+\frac{214745}{1728}\nu+\frac{193385}{3024}\nu^{2}\biggr{)}\pi
x^{7/2}$
$\displaystyle+\Biggl{[}-\frac{323105549467}{3178375200}+\frac{232597}{4410}\gamma_{\text{E}}-\frac{1369}{126}\pi^{2}+\frac{39931}{294}\ln
2-\frac{47385}{1568}\ln 3+\frac{232597}{8820}\ln x$
$\displaystyle\qquad+\biggl{(}-\frac{1452202403629}{1466942400}+\frac{41478}{245}\gamma_{\text{E}}-\frac{267127}{4608}\pi^{2}+\frac{479062}{2205}\ln
2+\frac{47385}{392}\ln 3+\frac{20739}{245}\ln x\biggr{)}\nu$
$\displaystyle\qquad+\biggl{(}\frac{1607125}{6804}-\frac{3157}{384}\pi^{2}\biggr{)}\nu^{2}+\frac{6875}{504}\nu^{3}+\frac{5}{6}\nu^{4}\Biggr{]}x^{4}$
$\displaystyle+\Biggl{[}\frac{265978667519}{745113600}-\frac{6848}{105}\gamma_{\text{E}}-\frac{3424}{105}\ln(16\,x)+\biggl{(}\frac{2062241}{22176}+\frac{41}{12}\pi^{2}\biggr{)}\nu$
$\displaystyle\qquad-\frac{133112905}{290304}\nu^{2}-\frac{3719141}{38016}\nu^{3}\Biggr{]}\pi
x^{9/2}+\mathcal{O}\bigl{(}x^{5}\bigr{)}\Biggr{\\}}\,.$ (4)

In the test-mass limit $\nu\to 0$, we exactly retrieve the result of linear
black-hole perturbation theory [53, 54, 55, 56, 57]. Since BH perturbations
have recently been extended numerically to second order in the mass ratio
$\nu$ [58, 59, 60], it would be interesting to verify the consistency of this
numerical result with the PN prediction (Gravitational-Wave Phasing of Compact
Binary Systems  to the Fourth-and-a-Half post-Newtonian Order). Note also that
in the case of black holes, the contributions due to the absorption by the BH
horizons are not included in the PN calculation, and should be added
separately (see [61] for Schwarzschild black holes, and [62, 63, 64, 65, 66]
for spinning ones).

With both (Gravitational-Wave Phasing of Compact Binary Systems  to the
Fourth-and-a-Half post-Newtonian Order) and (Gravitational-Wave Phasing of
Compact Binary Systems  to the Fourth-and-a-Half post-Newtonian Order) in
hand, we apply the flux-balance equation (2) and readily obtain the time
evolution of the GW frequency. Sophisticated techniques exist to increase the
precision of the PN results and the overlap with numerical relativity [67,
68], but we do not discuss them here and simply present the results in the
form of a fully expanded Taylor PN series. We employ the time variable

$\tau\equiv\frac{\nu c^{3}}{5Gm}\bigl{(}t_{0}-t\bigr{)}\,,$ (5)

where $t$ is the coordinate time in the asymptotic radiative coordinate
system, and $t_{0}$ an integration constant. We have the freedom to redefine
it as $t_{0}\longrightarrow t_{0}+\alpha\,\frac{Gm}{c^{3}}$ where $\alpha$ is
any constant, which amounts to the replacement
$\tau\longrightarrow\tau[1+\alpha\,\nu/(5\tau)]$. Although $t_{0}$ is not
uniquely defined, it might be formally interpreted as the instant of
coalescence, when $x\to+\infty$, and then it satisfies
$t_{0}-t=\mathcal{O}(c^{5})$ in the PN regime. At the 4PN order, using the
fact that $\tau^{-1}=\mathcal{O}(c^{-8})$ is a small 4PN quantity, we
conveniently adjust $\alpha$ so as to simplify as much as possible the result:

$\displaystyle x$
$\displaystyle=\frac{\tau^{-1/4}}{4}\Biggl{\\{}1+\biggl{(}\frac{743}{4032}+\frac{11}{48}\nu\biggr{)}\tau^{-1/4}-\frac{1}{5}\pi\,\tau^{-3/8}$
$\displaystyle\quad\quad+\biggl{(}\frac{19583}{254016}+\frac{24401}{193536}\nu+\frac{31}{288}\nu^{2}\biggr{)}\tau^{-1/2}+\biggl{(}-\frac{11891}{53760}+\frac{109}{1920}\nu\biggr{)}\pi\,\tau^{-5/8}$
$\displaystyle\quad\quad+\Biggl{[}-\frac{10052469856691}{6008596070400}+\frac{1}{6}\pi^{2}+\frac{107}{420}\gamma_{\text{E}}-\frac{107}{3360}\ln\biggl{(}\frac{\tau}{256}\biggr{)}$
$\displaystyle\quad\quad\qquad+\biggl{(}\frac{3147553127}{780337152}-\frac{451}{3072}\pi^{2}\biggr{)}\nu-\frac{15211}{442368}\nu^{2}+\frac{25565}{331776}\nu^{3}\Biggr{]}\tau^{-3/4}$
$\displaystyle\quad\quad+\biggl{(}-\frac{113868647}{433520640}-\frac{31821}{143360}\nu+\frac{294941}{3870720}\nu^{2}\biggr{)}\pi\tau^{-7/8}$
$\displaystyle\quad\quad+\Biggl{[}-\frac{2518977598355703073}{3779358859513036800}+\frac{9203}{215040}\gamma_{\text{E}}+\frac{9049}{258048}\pi^{2}+\frac{14873}{1128960}\ln
2+\frac{47385}{1605632}\ln 3-\frac{9203}{3440640}\ln\tau$
$\displaystyle\quad\quad\qquad+\Biggl{(}\frac{718143266031997}{576825222758400}+\frac{244493}{1128960}\gamma_{\text{E}}-\frac{65577}{1835008}\pi^{2}+\frac{15761}{47040}\ln
2-\frac{47385}{401408}\ln 3-\frac{244493}{18063360}\ln\tau\Biggr{)}\nu$
$\displaystyle\quad\quad\qquad+\Biggl{(}-\frac{1502014727}{8323596288}+\frac{2255}{393216}\pi^{2}\Biggr{)}\nu^{2}-\frac{258479}{33030144}\nu^{3}+\frac{1195}{262144}\nu^{4}\Biggr{]}\tau^{-1}\ln\tau$
$\displaystyle\quad\quad+\Biggl{[}-\frac{9965202491753717}{5768252227584000}+\frac{107}{600}\gamma_{\text{E}}+\frac{23}{600}\pi^{2}-\frac{107}{4800}\ln\biggl{(}\frac{\tau}{256}\biggr{)}$
$\displaystyle\quad\quad\qquad+\biggl{(}\frac{8248609881163}{2746786775040}-\frac{3157}{30720}\pi^{2}\biggr{)}\nu-\frac{3590973803}{20808990720}\nu^{2}-\frac{520159}{1634992128}\nu^{3}\Biggr{]}\pi\,\tau^{-9/8}+\mathcal{O}\bigl{(}\tau^{-5/4}\bigr{)}\Biggr{\\}}\,.$
(6)

Then, the GW phase $\psi$ of the dominant harmonics is related to the binary’s
orbital phase $\phi$ by

$\psi=\phi-\frac{2\pi
G\mathrm{M}f}{c^{3}}\ln\biggl{(}\frac{f}{f_{0}}\biggr{)}\,,$ (7)

where $\mathrm{M}$ is the ADM mass of the binary, and $f_{0}$ is an arbitrary
unphysical scale, reflecting the different origins of time between the local
coordinates covering the source and the radiative coordinates. The logarithmic
phase modulation was determined in [69, 70] and is physically due to the
scattering of GWs on the Schwarzschild background associated with
$\mathrm{M}$. While the GW phase $\psi$ and the corresponding GW frequency
$f=\dot{\psi}/\pi$ are directly measurable, the orbital phase $\phi$ can only
be inferred _via_ the theoretical prediction (7). The explicit expression of
the GW phase is

$\displaystyle\psi$
$\displaystyle=\psi_{0}-\frac{x^{-5/2}}{32\nu}\Biggl{\\{}1+\biggl{(}\frac{3715}{1008}+\frac{55}{12}\nu\biggr{)}x-10\pi
x^{3/2}$
$\displaystyle\quad+\biggl{(}\frac{15293365}{1016064}+\frac{27145}{1008}\nu+\frac{3085}{144}\nu^{2}\biggr{)}x^{2}+\biggl{(}\frac{38645}{1344}-\frac{65}{16}\nu\biggr{)}\pi
x^{5/2}\ln x$
$\displaystyle\quad+\Biggl{[}\frac{12348611926451}{18776862720}-\frac{160}{3}\pi^{2}-\frac{1712}{21}\gamma_{\text{E}}-\frac{856}{21}\ln(16\,x)$
$\displaystyle\quad\qquad+\biggl{(}-\frac{15737765635}{12192768}+\frac{2255}{48}\pi^{2}\biggr{)}\nu+\frac{76055}{6912}\nu^{2}-\frac{127825}{5184}\nu^{3}\Biggr{]}x^{3}$
$\displaystyle\quad+\biggl{(}\frac{77096675}{2032128}+\frac{378515}{12096}\nu-\frac{74045}{6048}\nu^{2}\biggr{)}\pi
x^{7/2}$
$\displaystyle\quad+\Biggl{[}\frac{2550713843998885153}{2214468081745920}-\frac{9203}{126}\gamma_{\text{E}}-\frac{45245}{756}\pi^{2}-\frac{252755}{2646}\ln
2-\frac{78975}{1568}\ln 3-\frac{9203}{252}\ln x$
$\displaystyle\quad\qquad+\biggl{(}-\frac{680712846248317}{337983528960}-\frac{488986}{1323}\gamma_{\text{E}}+\frac{109295}{1792}\pi^{2}-\frac{1245514}{1323}\ln
2+\frac{78975}{392}\ln 3-\frac{244493}{1323}\ln x\biggr{)}\nu$
$\displaystyle\quad\qquad\qquad+\biggl{(}\frac{7510073635}{24385536}-\frac{11275}{1152}\pi^{2}\biggr{)}\nu^{2}+\frac{1292395}{96768}\nu^{3}-\frac{5975}{768}\nu^{4}\Biggr{]}x^{4}$
$\displaystyle\quad+\Biggl{[}-\frac{93098188434443}{150214901760}+\frac{1712}{21}\gamma_{\text{E}}+\frac{80}{3}\pi^{2}+\frac{856}{21}\ln(16x)$
$\displaystyle\quad\qquad+\biggl{(}\frac{1492917260735}{1072963584}-\frac{2255}{48}\pi^{2}\biggr{)}\nu-\frac{45293335}{1016064}\nu^{2}-\frac{10323755}{1596672}\nu^{3}\Biggr{]}\pi
x^{9/2}+\mathcal{O}\bigl{(}x^{5}\bigr{)}\Biggr{\\}}\,,$ (8)

where the integration constant $\psi_{0}$ is determined by initial conditions,
_e.g._ , when the wave frequency enters the detector’s band. The results
hereabove (Gravitational-Wave Phasing of Compact Binary Systems  to the
Fourth-and-a-Half post-Newtonian Order)–(Gravitational-Wave Phasing of Compact
Binary Systems  to the Fourth-and-a-Half post-Newtonian Order) give the
prediction of Einstein’s general relativity for the GW frequency and phase
chirp of non-spinning compact binaries up to 4.5PN precision. On the other
hand, most of the frameworks for data analysis rely on the stationary phase
approximation (SPA) [71], for which the phase of the dominant quadrupole mode
reads

$\displaystyle\psi_{\text{SPA}}$ $\displaystyle=2\pi F\,T_{0}+\Psi_{0}$
$\displaystyle+\frac{3\,v^{-5}}{128\nu}\Biggl{\\{}1+\biggl{(}\frac{3715}{756}+\frac{55}{9}\nu\biggr{)}v^{2}-16\pi
v^{3}$
$\displaystyle\quad\qquad\quad+\biggl{(}\frac{15293365}{508032}+\frac{27145}{504}\nu+\frac{3085}{72}\nu^{2}\biggr{)}v^{4}+\biggl{(}\frac{38645}{252}-\frac{65}{3}\nu\biggr{)}\pi
v^{5}\ln v$
$\displaystyle\quad\qquad\quad+\Biggl{[}\frac{11583231236531}{4694215680}-\frac{640}{3}\pi^{2}-\frac{6848}{21}\gamma_{\text{E}}-\frac{6848}{21}\ln(4v)+\biggl{(}-\frac{15737765635}{3048192}+\frac{2255}{12}\pi^{2}\biggr{)}\nu$
$\displaystyle\quad\qquad\qquad\quad+\frac{76055}{12}\nu^{2}-\frac{127825}{1296}\nu^{3}\Biggr{]}v^{6}$
$\displaystyle\quad\qquad\quad+\Biggl{[}\frac{77096675}{254016}+\frac{378515}{1512}\nu-\frac{74045}{756}\nu^{2}\Biggr{]}\pi
v^{7}$
$\displaystyle\quad\qquad\quad+\Biggl{[}-\frac{2550713843998885153}{276808510218240}+\frac{90490}{189}\pi^{2}+\frac{36812}{63}\gamma_{\text{E}}+\frac{1011020}{1323}\ln
2+\frac{78975}{196}\ln 3+\frac{18406}{63}\ln v$
$\displaystyle\quad\qquad\qquad\quad+\biggl{(}\frac{680712846248317}{42247941120}-\frac{109295}{224}\pi^{2}+\frac{3911888}{1323}\gamma_{\text{E}}+\frac{9964112}{1323}\ln
2-\frac{78975}{49}\ln 3+\frac{1955944}{1323}\ln v\biggr{)}\nu$
$\displaystyle\quad\qquad\qquad\quad+\biggl{(}-\frac{7510073635}{3048192}+\frac{11275}{144}\pi^{2}\biggr{)}\nu^{2}-\frac{1292395}{12096}\nu^{3}+\frac{5975}{96}\nu^{4}\Biggr{]}v^{8}\ln
v$
$\displaystyle\quad\qquad\quad+\Biggl{[}\frac{105344279473163}{18776862720}-\frac{640}{3}\pi^{2}-\frac{13696}{21}\gamma_{\text{E}}-\frac{13696}{21}\ln(4v)$
$\displaystyle\quad\qquad\qquad\quad+\biggl{(}-\frac{1492917260735}{134120448}+\frac{2255}{6}\pi^{2}\biggr{)}\nu+\frac{45293335}{127008}\nu^{2}+\frac{10323755}{199584}\nu^{3}\Biggr{]}\pi
v^{9}+\mathcal{O}\bigl{(}v^{10}\bigr{)}\Biggr{\\}}\,,$ (9)

where $v\equiv\left(\frac{\pi GmF}{c^{3}}\right)^{1/3}$ with $F$ being the
Fourier frequency, and where $T_{0}$ and $\Psi_{0}$ are two integration
constants. Again we have adjusted $T_{0}$ in order to simplify the result (and
we have absorbed the usual $-\frac{\pi}{4}$ into $\Psi_{0}$). The coefficients
up to 3.5PN, as well as the 4.5PN piece, are already in use; see _e.g._ App. A
of [72].

In order to get intuition on the relative contribution of each PN order to the
signal, we provide in Table 1 rough numerical estimates for the number of
accumulated GW cycles in the frequency band of current and future detectors.
Our naive estimation does not take the various detector noises into account,
and a more realistic estimation should be performed [73]. Nevertheless, it can
be useful to gain insight on the behavior of the PN expansion, which seems to
converge well, as we see from Table 1. For all the typical compact binaries in
Table 1, we find that the 4PN and 4.5PN orders amount to about a tenth of a
cycle (less than 1 radian). This suggests that systematic errors due to the PN
modeling may be dominated by statistical errors and negligible for LISA.
However, this should be confirmed by detailed investigations along the lines
of [74].

  Detector LIGO/Virgo ET LISA Masses ($M_{\odot}$) $1.4\times 1.4$
$10\times 10$ $1.4\times 1.4$ $500\times 500$ $10^{5}\times 10^{5}$
$10^{7}\times 10^{7}$ PN order cumulative number of cycles Newtonian
$2\,562.599$ $95.502$ $744\,401.36$ $37.90$ $28\,095.39$ $9.534$ 1PN $143.453$
$17.879$ $4\,433.85$ $9.60$ $618.31$ $3.386$ 1.5PN $-94.817$ $-20.797$
$-1\,005.78$ $-12.63$ $-265.70$ $-5.181$ 2PN $5.811$ $2.124$ $23.94$ $1.44$
$11.35$ $0.677$ 2.5PN $-8.105$ $-4.604$ $-17.01$ $-3.42$ $-12.47$ $-1.821$ 3PN
$1.858$ $1.731$ $2.69$ $1.43$ $2.59$ $0.876$ 3.5PN $-0.627$ $-0.689$ $-0.93$
$-0.59$ $-0.91$ $-0.383$ 4PN $-0.107$ $-0.064$ $-0.12$ $-0.04$ $-0.12$
$-0.013$ 4.5PN $0.098$ $0.118$ $0.14$ $0.10$ $0.14$ $0.065$

Table 1: Contribution of each PN order to the total number of accumulated
cycles inside the detector’s frequency band, for typical (but non-spinning)
quasi-circular compact binaries observed by current and future detectors. We
have approximated the frequency bands of LIGO/Virgo, Einstein Telescope (ET)
and LISA with step functions, respectively between
$\bigl{[}30\,\text{Hz},10^{3}\,\text{Hz}\bigr{]}$,
$\bigl{[}1\,\text{Hz},10^{4}\,\text{Hz}\bigr{]}$ and
$\bigl{[}10^{-4}\,\text{Hz},10^{-1}\,\text{Hz}\bigr{]}$. When the merger
occurs within the frequency band of the detector, the exit frequency is taken
to be the Schwarzschild ISCO, $f_{\text{ISCO}}=c^{3}/(6^{3/2}\pi Gm)$. The
contributions due to the non-linearities of GR (_e.g._ , tails) increase with
the PN order and are detailed in [22].

Besides the chirp described by the results (Gravitational-Wave Phasing of
Compact Binary Systems  to the Fourth-and-a-Half post-Newtonian
Order)–(Gravitational-Wave Phasing of Compact Binary Systems  to the Fourth-
and-a-Half post-Newtonian Order), it is also important to compute the wave
amplitude, in view of the data analysis of LISA [75, 76, 77] and high-accuracy
comparisons with numerical relativity (see _e.g._ [78, 79, 80, 81]). We
decompose the waveform, at leading order in the distance $R$ to the source,
onto a basis of spin-weighted spherical harmonics (following the conventions
of [82, 83])

$h_{+}-\mathrm{i}h_{\times}=\frac{8Gm\nu
x}{Rc^{2}}\,\sqrt{\frac{\pi}{5}}\,\sum_{\ell=2}^{+\infty}\sum_{\mathrm{m}=-\ell}^{\ell}H_{\ell\mathrm{m}}e^{-\mathrm{i}\mathrm{m}\psi}Y_{-2}^{\ell\mathrm{m}}\,,$
(10)

where the phase variable is given by (Gravitational-Wave Phasing of Compact
Binary Systems  to the Fourth-and-a-Half post-Newtonian Order). All
$H_{\ell\mathrm{m}}$ modes are currently known at 3.5PN order for spinning,
non-precessing, quasi-circular orbits [84, 82, 83, 85, 86]. Although we were
able to derive the phase with 4.5PN accuracy, the same precision for the modes
is yet out of reach, since, even though the 4.5PN radiation-reaction terms in
the equations of motion are known [87, 88], neither the source quadrupole
moment nor the non-linear contributions to the GW propagation are fully
controlled at 4.5PN order (only the contributions that enter the 4.5PN flux
for circular orbits are known). We thus report the extension of the dominant
quadrupole mode $(\ell,\mathrm{m})=(2,2)$ for non-spinning, quasi-circular
orbits up to 4PN order:

$\displaystyle H_{22}$
$\displaystyle=1+\biggl{(}-\frac{107}{42}+\frac{55}{42}\nu\biggr{)}x+2\pi
x^{3/2}+\biggl{(}-\frac{2173}{1512}-\frac{1069}{216}\nu+\frac{2047}{1512}\nu^{2}\biggr{)}x^{2}+\left[-\frac{107\pi}{21}+\left(\frac{34\pi}{21}-24\,\mathrm{i}\right)\nu\right]x^{5/2}$
$\displaystyle\qquad+\Biggl{[}\frac{27027409}{646800}-\frac{856}{105}\,\gamma_{\text{E}}+\frac{428\,\mathrm{i}\,\pi}{105}+\frac{2\pi^{2}}{3}+\biggl{(}-\frac{278185}{33264}+\frac{41\pi^{2}}{96}\biggr{)}\nu-\frac{20261}{2772}\nu^{2}+\frac{114635}{99792}\nu^{3}-\frac{428}{105}\ln(16x)\Biggr{]}x^{3}$
$\displaystyle\qquad+\Biggl{[}-\frac{2173\pi}{756}+\biggl{(}-\frac{2495\pi}{378}+\frac{14333\,\mathrm{i}}{162}\biggr{)}\nu+\biggl{(}\frac{40\pi}{27}-\frac{4066\,\mathrm{i}}{945}\biggr{)}\nu^{2}\Biggr{]}x^{7/2}$
$\displaystyle\qquad+\Biggl{[}-\frac{846557506853}{12713500800}+\frac{45796}{2205}\gamma_{\text{E}}-\frac{22898}{2205}\mathrm{i}\pi-\frac{107}{63}\pi^{2}+\frac{22898}{2205}\ln(16x)$
$\displaystyle\qquad\qquad+\biggl{(}-\frac{336005827477}{4237833600}+\frac{15284}{441}\gamma_{\text{E}}-\frac{219314}{2205}\mathrm{i}\pi-\frac{9755}{32256}\pi^{2}+\frac{7642}{441}\ln(16x)\biggr{)}\nu$
$\displaystyle\qquad\qquad+\biggl{(}\frac{256450291}{7413120}-\frac{1025}{1008}\pi^{2}\biggr{)}\nu^{2}-\frac{81579187}{15567552}\nu^{3}+\frac{26251249}{31135104}\nu^{4}\Biggr{]}x^{4}+\mathcal{O}\bigl{(}x^{9/2}\bigr{)}\,.$
(11)

Satisfyingly, this result is in perfect agreement with linear black-hole
perturbation theory in the limit when $\nu\to 0$; see App. B of [54]. Again,
it would be interesting to compare the PN prediction (Gravitational-Wave
Phasing of Compact Binary Systems  to the Fourth-and-a-Half post-Newtonian
Order) with second-order (numerical or analytical) BH perturbation theory.

###### Acknowledgements.

We thank Adam Pound for discussions on BH perturbation and GSF results, Bala
Iyer for his longstanding support and interest in this project, and Alessandra
Buonanno for suggestions. F.L. is grateful to the Institut d’Astrophysique de
Paris for its hospitality during this project. He received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (grant agreement No 817791). G.F. thanks IIT
Madras for a Visiting Faculty Fellow position under the IoE program during the
completion of this work.

## References

  * Einstein [1918] A. Einstein, Über gravitationswellen, Sitzber. Preuss. Akad. Wiss. Berlin 1, 154 (1918).
  * Landau and Lifshitz [1971] L. Landau and E. Lifshitz,  _The classical theory of fields_ (Pergamon, Oxford, 1971).
  * Peters and Mathews [1963] P. Peters and J. Mathews, Gravitational radiation from point masses in a keplerian orbit, Phys. Rev. 131, 435 (1963).
  * Cutler  _et al._ [1993] C. Cutler, T. Apostolatos, L. Bildsten, L. Finn, E. Flanagan, D. Kennefick, D. Markovic, A. Ori, E. Poisson, G. Sussman, and K. Thorne, The last three minutes: Issues in gravitational-wave measurements of coalescing compact binaries, Phys. Rev. Lett. 70, 2984 (1993).
  * Cutler and Flanagan [1994] C. Cutler and E. Flanagan, Gravitational waves from merging compact binaries: How accurately can one extract the binary’s parameters from the inspiral waveform?, Phys. Rev. D 49, 2658 (1994).
  * Blanchet and Damour [1986] L. Blanchet and T. Damour, Radiative gravitational fields in general relativity. i. general structure of the field outside the source, Phil. Trans. Roy. Soc. Lond. A 320, 379 (1986).
  * Blanchet [1987] L. Blanchet, Radiative gravitational fields in general relativity. II. Asymptotic behaviour at future null infinity, Proc. Roy. Soc. Lond. A 409, 383 (1987).
  * Blanchet and Damour [1988] L. Blanchet and T. Damour, Tail-transported temporal correlations in the dynamics of a gravitating system, Phys. Rev. D 37, 1410 (1988).
  * Blanchet and Damour [1992] L. Blanchet and T. Damour, Hereditary effects in gravitational radiation, Phys. Rev. D 46, 4304 (1992).
  * Blanchet [1998a] L. Blanchet, On the multipole expansion of the gravitational field, Class. Quant. Grav. 15, 1971 (1998a), gr-qc/9801101 .
  * Blanchet  _et al._ [1995a] L. Blanchet, T. Damour, and B. R. Iyer, Gravitational waves from inspiralling compact binaries: Energy loss and wave form to second post-Newtonian order, Phys. Rev. D 51, 5360 (1995a), gr-qc/9501029 .
  * Blanchet  _et al._ [1995b] L. Blanchet, T. Damour, B. R. Iyer, C. M. Will, and A. G. Wiseman, Gravitational radiation damping of compact binary systems to second post-Newtonian order, Phys. Rev. Lett. 74, 3515 (1995b), gr-qc/9501027 .
  * Blanchet  _et al._ [1996] L. Blanchet, B. R. Iyer, C. M. Will, and A. G. Wiseman, Gravitational wave forms from inspiralling compact binaries to second-post-Newtonian order, Class. Quant. Grav. 13, 575 (1996), gr-qc/9602024 .
  * Blanchet [1996] L. Blanchet, Energy losses by gravitational radiation in inspiralling compact binaries to five halves post-Newtonian order, Phys. Rev. D 54, 1417 (1996), erratum Phys. Rev. D 71, 129904(E) (2005), gr-qc/9603048 .
  * Blanchet [1998b] L. Blanchet, Gravitational-wave tails of tails, Class. Quant. Grav. 15, 113 (1998b), gr-qc/9710038 .
  * Blanchet  _et al._ [2002a] L. Blanchet, B. R. Iyer, and B. Joguet, Gravitational waves from inspiralling compact binaries: Energy flux to third post-Newtonian order, Phys. Rev. D 65, 064005 (2002a), erratum Phys. Rev. D, 71:129903(E), 2005, gr-qc/0105098 .
  * Blanchet  _et al._ [2002b] L. Blanchet, G. Faye, B. R. Iyer, and B. Joguet, Gravitational-wave inspiral of compact binary systems to 7/2 post-Newtonian order, Phys. Rev. D 65, 061501(R) (2002b), erratum Phys. Rev. D, 71:129902(E), 2005, gr-qc/0105099 .
  * Blanchet and Iyer [2004] L. Blanchet and B. R. Iyer, Hadamard regularization of the third post-Newtonian gravitational wave generation of two point masses, Phys. Rev. D 71, 024004 (2004), gr-qc/0409094 .
  * Blanchet  _et al._ [2004] L. Blanchet, T. Damour, G. Esposito-Farèse, and B. R. Iyer, Gravitational radiation from inspiralling compact binaries completed at the third post-Newtonian order, Phys. Rev. Lett. 93, 091101 (2004), gr-qc/0406012 .
  * Blanchet  _et al._ [2005] L. Blanchet, T. Damour, G. Esposito-Farèse, and B. R. Iyer, Dimensional regularization of the third post-Newtonian gravitational wave generation of two point masses, Phys. Rev. D 71, 124004 (2005), gr-qc/0503044 .
  * [21] The ancillary file GW_4PN_supplementary_material.m collects the results associated with the gravitational-wave phase at 4.5PN and the $(\ell,\text{m})=(2,2)$ mode at 4PN.
  * Blanchet  _et al._ [2023] L. Blanchet, G. Faye, Q. Henry, F. Larrouturou, and D. Trestini, Gravitational Wave Flux and Quadrupole Modes from Quasi-Circular Non-Spinning Compact Binaries to the Fourth Post-Newtonian Order (2023), arXiv:2304.11186 [gr-qc] .
  * Blanchet and Sathyaprakash [1995] L. Blanchet and B. S. Sathyaprakash, Detecting the tail effect in gravitational wave experiments, Phys. Rev. Lett. 74, 1067 (1995).
  * Arun  _et al._ [2006] K. Arun, B. Iyer, M. Qusailah, and B. Sathyaprakash, Testing post-Newtonian theory with gravitational wave observations, Class. Quant. Gravity 23, L37 (2006), gr-qc/0604018 .
  * Abbott  _et al._ [2016] B. Abbott _et al._ (LIGO Scientific Collaboration and VIRGO Collaboration), Tests of general relativity with GW150914, Phys. Rev. Lett. 116, 221101 (2016), arXiv:1602.03841 [gr-qc] .
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Tests of General Relativity with GW170817, Phys. Rev. Lett. 123, 011102 (2019), arXiv:1811.00364 [gr-qc] .
  * Sesana [2016] A. Sesana, Prospects for Multiband Gravitational-Wave Astronomy after GW150914, Phys. Rev. Lett. 116, 231102 (2016), arXiv:1602.06951 [gr-qc] .
  * Peters [1964] P. Peters, Gravitational radiation and the motion of two point masses, Phys. Rev. 136, B1224 (1964).
  * Jaranowski and Schäfer [2012] P. Jaranowski and G. Schäfer, Towards the fourth post-Newtonian Hamiltonian for two-point-mass systems, Phys. Rev. D 86, 061503(R) (2012), arXiv:1207.5448 [gr-qc] .
  * Jaranowski and Schäfer [2013] P. Jaranowski and G. Schäfer, Dimensional regularization of local singularities in the 4th post-Newtonian two-point-mass Hamiltonian, Phys. Rev. D 87, 081503(R) (2013), arXiv:1303.3225 [gr-qc] .
  * Damour  _et al._ [2014] T. Damour, P. Jaranowski, and G. Schäfer, Non-local-in-time action for the fourth post-Newtonian conservative dynamics of two-body systems, Phys. Rev. D 89, 064058 (2014), arXiv:1401.4548 [gr-qc] .
  * Damour  _et al._ [2016] T. Damour, P. Jaranowski, and G. Schäfer, On the conservative dynamics of two-body systems at the fourth post-Newtonian approximation of general relativity, Phys. Rev. D 93, 084014 (2016), arXiv:1601.01283 [gr-qc] .
  * Bini and Damour [2013] D. Bini and T. Damour, Analytical determination of the two-body gravitational interaction potential at the fourth post-Newtonian approximation, Phys. Rev. D 87, 121501(R) (2013), arXiv:1305.4884 [gr-qc] .
  * Bernard  _et al._ [2016] L. Bernard, L. Blanchet, A. Bohé, G. Faye, and S. Marsat, Fokker action of non-spinning compact binaries at the fourth post-Newtonian approximation, Phys. Rev. D 93, 084037 (2016), arXiv:1512.02876 [gr-qc] .
  * Bernard  _et al._ [2017a] L. Bernard, L. Blanchet, A. Bohé, G. Faye, and S. Marsat, Energy and periastron advance of compact binaries on circular orbits at the fourth post-Newtonian order, Phys. Rev. D 95, 044026 (2017a), arXiv:1610.07934 [gr-qc] .
  * Bernard  _et al._ [2017b] L. Bernard, L. Blanchet, A. Bohé, G. Faye, and S. Marsat, Dimensional regularization of the IR divergences in the Fokker action of point-particle binaries at the fourth post-Newtonian order, Phys. Rev. D 96, 104043 (2017b), arXiv:1706.08480 [gr-qc] .
  * Bernard  _et al._ [2018] L. Bernard, L. Blanchet, G. Faye, and T. Marchand, Center-of-mass equations of motion and conserved integrals of compact binary systems at the fourth post-Newtonian order, Phys. Rev. D 97, 044037 (2018), arXiv:1711.00283 [gr-qc] .
  * Foffa and Sturani [2013] S. Foffa and R. Sturani, The dynamics of the gravitational two-body problem in the post-Newtonian approximation at quadratic order in the Newton’s constant, Phys. Rev. D 87, 064011 (2013), arXiv:1206.7087 [gr-qc] .
  * Foffa and Sturani [2012] S. Foffa and R. Sturani, Tail terms in gravitational radiation reaction via effective field theory, Phys. Rev. D 87, 044056 (2012), arXiv:1111.5488 [gr-qc] .
  * Galley  _et al._ [2016] C. R. Galley, A. K. Leibovich, R. A. Porto, and A. Ross, Tail effect in gravitational radiation reaction: Time nonlocality and renormalization group evolution, Phys. Rev. D 93, 124010 (2016), arXiv:1511.07379 [gr-qc] .
  * Foffa  _et al._ [2017] S. Foffa, P. Mastrolia, R. Sturani, and C. Sturm, Effective field theory approach to the gravitational two-body dynamics at fourth post-Newtonian order and quintic in the Newton constant, Phys. Rev. D 95, 104009 (2017), arXiv:1612.00482 [gr-qc] .
  * Porto and Rothstein [2017] R. Porto and I. Rothstein, On the apparent ambiguities in the post-Newtonian expansion for binary systems, Phys. Rev. D 96, 024061 (2017), arXiv:1703.06463 [gr-qc] .
  * Foffa and Sturani [2019] S. Foffa and R. Sturani, Conservative dynamics of binary systems to fourth post-Newtonian order in the EFT approach I: Regularized Lagrangian, Phys. Rev. D 100, 024047 (2019), arXiv:1903.05113 [gr-qc] .
  * Foffa  _et al._ [2019] S. Foffa, R. Porto, I. Rothstein, and R. Sturani, Conservative dynamics of binary systems to fourth post-Newtonian order in the EFT approach II: Renormalized Lagrangian, Phys. Rev. D 100, 024048 (2019), arXiv:1903.05118 [gr-qc] .
  * Blümlein  _et al._ [2020] J. Blümlein, A. Maier, P. Marquard, and G. Schäfer, Fourth post-Newtonian Hamiltonian dynamics of two-body systems from an effective field theory approach, Nucl. Phys. B 955, 115041 (2020), arXiv:2003.01692 [gr-qc] .
  * Marchand  _et al._ [2020] T. Marchand, Q. Henry, F. Larrouturou, S. Marsat, G. Faye, and L. Blanchet, The mass quadrupole moment of compact binary systems at the fourth post-Newtonian order, Class. Quant. Grav. 37, 215006 (2020), arXiv:2003.13672 [gr-qc] .
  * Larrouturou  _et al._ [2022a] F. Larrouturou, Q. Henry, L. Blanchet, and G. Faye, The quadrupole moment of compact binaries to the fourth post-Newtonian order: I. Non-locality in time and infra-red divergencies, Class. Quant. Grav. 39, 115007 (2022a), arXiv:2110.02240 [gr-qc] .
  * Larrouturou  _et al._ [2022b] F. Larrouturou, L. Blanchet, Q. Henry, and G. Faye, The quadrupole moment of compact binaries to the fourth post-Newtonian order: II. Dimensional regularization and renormalization, Class. Quant. Grav. 39, 115008 (2022b), arXiv:2110.02243 [gr-qc] .
  * Blanchet  _et al._ [2022] L. Blanchet, G. Faye, and F. Larrouturou, The quadrupole moment of compact binaries to the fourth post-Newtonian order: from source to canonical moment, Class. Quant. Grav. 39, 195003 (2022), arXiv:2204.11293 [gr-qc] .
  * Henry  _et al._ [2021] Q. Henry, G. Faye, and L. Blanchet, The current-type quadrupole moment and gravitational-wave mode ($\ell$, m) = (2, 1) of compact binary systems at the third post-Newtonian order, Class. Quant. Grav. 38, 185004 (2021), arXiv:2105.10876 [gr-qc] .
  * Trestini  _et al._ [2023] D. Trestini, F. Larrouturou, and L. Blanchet, The quadrupole moment of compact binaries to the fourth post-Newtonian order: relating the harmonic and radiative metrics, Class. Quant. Grav. 40, 055006 (2023), arXiv:2209.02719 [gr-qc] .
  * Trestini and Blanchet [2023] D. Trestini and L. Blanchet, Gravitational-wave tails of memory (2023), arXiv:2301.09395 [gr-qc] .
  * Sasaki [1994] M. Sasaki, Post-Newtonian expansion of the ingoing-wave regge-wheeler function, Prog. Theor. Phys. 92, 17 (1994).
  * Tagoshi and Sasaki [1994] H. Tagoshi and M. Sasaki, Post-Newtonian expansion of gravitational-waves from a particle in circular orbit around a schwarzschild black-hole, Prog. Theor. Phys. 92, 745 (1994), gr-qc/9405062 .
  * Tanaka  _et al._ [1996] T. Tanaka, H. Tagoshi, and M. Sasaki, Gravitational waves by a particle in circular orbit around a schwarzschild black hole: 5.5 Post-Newtonian formula, Prog. Theor. Phys. 96, 1087 (1996), gr-qc/9701050 .
  * Fujita [2012a] R. Fujita, Gravitational radiation for extreme mass ratio inspirals to the 14th post-Newtonian order, Prog. Theor. Phys. 127, 583 (2012a), arXiv:1104.5615 [gr-qc] .
  * Fujita [2012b] R. Fujita, Gravitational waves from a particle in circular orbits around a schwarzschild black hole to the 22nd post-Newtonian order, Prog. Theor. Phys. 128, 971 (2012b), arXiv:1211.5535 [gr-qc] .
  * Wardell  _et al._ [2021] B. Wardell, A. Pound, N. Warburton, J. Miller, L. Durkan, and A. Le Tiec, Gravitational waveforms for compact binaries from second-order self-force theory (2021), arXiv:2112.12265 [gr-qc] .
  * Albertini  _et al._ [2022a] A. Albertini, A. Nagar, A. Pound, N. Warburton, B. Wardell, L. Durkan, and J. Miller, Comparing second-order gravitational self-force, numerical relativity, and effective one body waveforms from inspiralling, quasicircular, and nonspinning black hole binaries, Phys. Rev. D 106, 084061 (2022a), arXiv:2208.01049 [gr-qc] .
  * Albertini  _et al._ [2022b] A. Albertini, A. Nagar, A. Pound, N. Warburton, B. Wardell, L. Durkan, and J. Miller, Comparing second-order gravitational self-force and effective one body waveforms from inspiralling, quasicircular and nonspinning black hole binaries. II. The large-mass-ratio case, Phys. Rev. D 106, 084062 (2022b), arXiv:2208.02055 [gr-qc] .
  * Poisson and Sasaki [1995] E. Poisson and M. Sasaki, Gravitational radiation from a particle in circular orbit around a black hole. v. black-hole absorption and tail corrections, Phys. Rev. D 51, 5753 (1995), gr-qc/9412027 .
  * Tagoshi  _et al._ [1997] H. Tagoshi, S. Mano, and E. Takasugi, Post-Newtonian expansion of gravitational waves from a particle in circular orbits around a rotating black hole, Prog. Theor. Phys. 98, 829 (1997), gr-qc/9711072 .
  * Alvi [2001] K. Alvi, Energy and angular momentum flow into a black hole in a binary, Phys. Rev. D 64, 104020 (2001), arXiv:gr-qc/0107080 .
  * Porto [2008] R. A. Porto, Absorption effects due to spin in the worldline approach to black hole dynamics, Phys. Rev. D 77, 064026 (2008), arXiv:0710.5150 [hep-th] .
  * Chatziioannou  _et al._ [2013] K. Chatziioannou, E. Poisson, N. Yunes, K. Chatziioannou, E. Poisson, and N. Yunes, Tidal heating and torquing of a kerr black hole to next-to-leading order in the tidal coupling, Phys. Rev. D 87, 044022 (2013), arXiv:1211.1686 [gr-qc] .
  * Saketh  _et al._ [2022] M. V. S. Saketh, J. Steinhoff, J. Vines, and A. Buonanno, Modeling horizon absorption in spinning binary black holes using effective worldline theory (2022), arXiv:2212.13095 [gr-qc] .
  * Damour  _et al._ [1998] T. Damour, B. R. Iyer, and B. S. Sathyaprakash, Improved filters for gravitational waves from inspiralling compact binaries, Phys. Rev. D 57, 885 (1998), arXiv:gr-qc/9708034 .
  * Buonanno  _et al._ [2009] A. Buonanno, B. Iyer, Y. Pan, E. Ochsner, and B. Sathyaprakash, Comparison of post-Newtonian templates for compact binary inspiral signals in gravitational-wave detectors, Phys. Rev. D 80, 084043 (2009), arXiv:0907.0700 [gr-qc] .
  * Wiseman [1993] A. Wiseman, Coalescing binary-systems of compact objects to (post)5/2-Newtonian order. iv. the gravitational-wave tail, Phys. Rev. D 48, 4757 (1993).
  * Blanchet and Schäfer [1993] L. Blanchet and G. Schäfer, Gravitational wave tails and binary star systems, Class. Quant. Grav. 10, 2699 (1993).
  * Tichy  _et al._ [2000] W. Tichy, E. E. Flanagan, and E. Poisson, Can the postNewtonian gravitational wave form of an inspiraling binary be improved by solving the energy balance equation numerically?, Phys. Rev. D 61, 104015 (2000), arXiv:gr-qc/9912075 .
  * Pratten  _et al._ [2020] G. Pratten, S. Husa, C. Garcia-Quiros, M. Colleoni, A. Ramos-Buades, H. Estelles, and R. Jaume, Setting the cornerstone for a family of models for gravitational waves from compact binaries: The dominant harmonic for nonprecessing quasicircular black holes, Phys. Rev. D 102, 064001 (2020), arXiv:2001.11412 [gr-qc] .
  * Arun  _et al._ [2005] K. Arun, B. Iyer, B. Sathyaprakash, and P. Sundararajan, Parameter estimation of inspiralling compact binaries using 3.5 post-Newtonian gravitational wave phasing: The non-spinning case, Phys. Rev. D 71, 084008 (2005), gr-qc/0411146 .
  * Owen  _et al._ [2023] C. B. Owen, C.-J. Haster, S. Perkins, N. J. Cornish, and N. Yunes, Waveform accuracy and systematic uncertainties in current gravitational wave observations, arXiv e-prints (2023), arXiv:2301.11941 [gr-qc] .
  * Hellings and Moore [2003] R. W. Hellings and T. A. Moore, The Information content of gravitational wave harmonics in compact binary inspiral, Class. Quant. Grav. 20, S181 (2003), arXiv:gr-qc/0207102 .
  * Moore and Hellings [2000] T. A. Moore and R. W. Hellings, The Angular resolution of space based gravitational wave detectors, AIP Conf. Proc. 523, 255 (2000), arXiv:gr-qc/9910116 .
  * Arun  _et al._ [2007] K. G. Arun, B. R. Iyer, B. S. Sathyaprakash, S. Sinha, and C. Van Den Broeck, Higher signal harmonics, LISA’s angular resolution and dark energy, Phys. Rev. D 76, 104016 (2007), [Erratum: Phys.Rev.D 76, 129903 (2007)], arXiv:0707.3920 [astro-ph] .
  * Grandclement and Novak [2009] P. Grandclement and J. Novak, Spectral methods for numerical relativity, Living Rev. Rel. 12, 1 (2009), arXiv:0706.2286 [gr-qc] .
  * Boyle  _et al._ [2007] M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mroue, H. P. Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky, High-accuracy comparison of numerical relativity simulations with post-Newtonian expansions, Phys. Rev. D 76, 124038 (2007), arXiv:0710.0158 [gr-qc] .
  * Bernuzzi  _et al._ [2012] S. Bernuzzi, M. Thierfelder, and B. Bruegmann, Accuracy of numerical relativity waveforms from binary neutron star mergers and their comparison with post-Newtonian waveforms, Phys. Rev. D 85, 104030 (2012), arXiv:1109.3611 [gr-qc] .
  * Borhanian  _et al._ [2020] S. Borhanian, K. G. Arun, H. P. Pfeiffer, and B. S. Sathyaprakash, Comparison of post-Newtonian mode amplitudes with numerical relativity simulations of binary black holes, Class. Quant. Grav. 37, 065006 (2020), arXiv:1901.08516 [gr-qc] .
  * Blanchet  _et al._ [2008] L. Blanchet, G. Faye, B. R. Iyer, and S. Sinha, The third post-Newtonian gravitational wave polarisations and associated spherical harmonic modes for inspiralling compact binaries in quasi-circular orbits, Class. Quant. Grav. 25, 165003 (2008), arXiv:0802.1249 [gr-qc] .
  * Faye  _et al._ [2012] G. Faye, S. Marsat, L. Blanchet, and B. R. Iyer, The third and a half post-Newtonian gravitational wave quadrupole mode for quasi-circular inspiralling compact binaries, Class. Quant. Grav. 29, 175004 (2012), arXiv:1204.1043 [gr-qc] .
  * Arun  _et al._ [2004] K. Arun, L. Blanchet, B. R. Iyer, and M. S. Qusailah, The 2.5PN gravitational wave polarisations from inspiralling compact binaries in circular orbits, Class. Quant. Grav. 21, 3771 (2004), erratum Class. Quant. Grav., 22:3115, 2005, gr-qc/0404085 .
  * Henry  _et al._ [2022] Q. Henry, S. Marsat, and M. Khalil, Spin contributions to the gravitational-waveform modes for spin-aligned binaries at the 3.5PN order, Phys. Rev. D 106, 124018 (2022), arXiv:2209.00374 [gr-qc] .
  * Henry [2023] Q. Henry, Complete gravitational-waveform amplitude modes for quasicircular compact binaries to the 3.5PN order, Phys. Rev. D 107, 044057 (2023), arXiv:2210.15602 [gr-qc] .
  * Gopakumar  _et al._ [1997] A. Gopakumar, B. R. Iyer, and S. Iyer, Second post-Newtonian gravitational radiation reaction for two-body systems: Nonspinning bodies, Phys. Rev. D 55, 6030 (1997), gr-qc/9703075 .
  * Leibovich  _et al._ [2023] A. K. Leibovich, B. A. Pardo, and Z. Yang, Radiation Reaction for Non-Spinning Bodies at 4.5PN in the Effective Field Theory Approach (2023), arXiv:2302.11016 [gr-qc] .