# Flexural wave modulation and mitigation in airfoils using acoustic black
holes

Kaushik Sampath<EMAIL_ADDRESS>Caleb F Sieck  Matthew D Guild
Alec K Ikei  U.S. Naval Research Laboratory, Code 7165, Washington DC 20375,
USA  Charles A Rohde  U.S. Naval Research Laboratory, Code 6364, Washington DC
20375, USA

###### Abstract

This study introduces a framework for the design and implementation of
acoustic black holes (ABHs) in airfoils. A generalized multi-parameter damped-
ABH generation function is mapped onto NACA series airfoils. Representative
geometries and a uniformly distributed baseline, all with the same mass of
structure and damping are fabricated using multi-material PolyJet 3D printing.
Laser Doppler vibrometer measurements along the airfoil chord in response to a
broadband 0.1 - 12 kHz excitation show a decrease in trailing edge vibrations
by as much as 10 dB, a broadband 5 dB reduction across the entire chord as
well as substantial spatial and temporal modulation of flexural waves by ABH-
embedded foils. Finite element analysis (FEA) models are developed and
validated based on the measured data. Furthermore, a parametric FEA study is
performed on a set of comparable designs to elucidate the scope of modulation
achievable. These findings are applicable to trailing-edge noise reduction,
flow control, structural enhancement and energy harvesting for airfoils.

##  I  Introduction

Fluid-loaded structures such as turbomachine blades and aircraft wings are
often designed slender due to constraints on weight, making them susceptible
to vibrational excitation by flow [1]. This leads to an increased wear of
structures, affecting longevity and performance. Decades of past and ongoing
research has been aimed at finding better ways to mitigate such undesirable
consequences [2]. Substantial efforts have also gone into the redistribution
and harvesting of flow-induced vibrations for controlling turbulent flow and
overall energy efficiency. Trailing edge noise, which is in fact, a subset of
the above, still remains an active topic of research due to its relevance to
airframes, propellers and rotors [3, 4, 5].

New approaches of structural geometry modifications, such as, applying the so-
called acoustic black hole (ABH) effect have become increasingly popular.
Mironov [6] theorized that flexural waves can be ‘trapped’ in a beam with an
ideal power law-shaped tapering end. This is because the group velocity goes
to zero as it scales with square root of the edge thickness. In practice, this
is leveraged by adding viscoelastic damping wherever the taper truncates [7,
8, 9, 10, 11, 12]. A detailed review of ABH theory and applications has
recently been carried out by Pelat et al [13]. Evidently, despite their
popularity, applications of ABHs have been largely restricted to beams, plates
and more recently, cylindrical shells [14].

As far as aerodynamic applications are concerned, 1D power-law tapers have
been successfully incorporated into the trailing edge of turbo-fan blades
[15]. The study found that the measured acceleration, especially around
resonances, was substantially reduced in airfoils with power law tapers (ABHs)
when compared to those without. As this study notes, fabricating a trailing
edge with a power law taper is not trivial, but still achievable in an
otherwise seamless manner. However, for other use cases, where flow-structure
interactions may need to be modulated in sections of the chord besides the
trailing edge, structural modifications would need to be concealed internally
without affecting the aerodynamic external shape of the airfoils.

In fact, only recently have numerical studies even looked at embedding ABHs in
higher dimensional closed geometries such as cylindrical shells and beams [16,
17, 14]. Deng et al. [14] compare the performance of a uniformly damped
cylindrical steel shell with that of a ten-element ABH-embedded shell with the
same damping layer thickness. They also evaluate the effect of the truncation
thickness and found that even for their thickest (least favorable) truncation
case there is a 10 dB reduction in transmissibility of flexural vibrations
when compared to the uniformly damped case. It is important to reiterate that
unlike the previously mentioned turbo-fan blade edges [15], cylindrical shells
and beams do not have an obvious site to embed ABH tapers and there is a large
impact on the overall rigidity of the shells. A work-around proposed by Deng
et al. [14] is the addition of specially directed ‘stiffeners’ that
effectively support the structure without compromising the ABH effect. These
numerical studies, to the best of the authors’ knowledge, remain the only
examples of embedded (or closed-geometry) implementations of ABHs in higher-
dimensions. It still remains to be seen whether ABH-embedded designs can be
fabricated and demonstrated as such.

Recent advances in additive manufacturing, such as multi-material PolyJet
printing enable rapid fabrication of complex designs with hard and soft
materials in a single build. Subsequently, power law tapers, including
functionally graded ABHs, have been 3D-printed recently [18] for beams.
Comparisons of the measured reflection coefficients between a traditional (or
single material) ABH beam and those spatially distributed with softer and
higher loss materials towards the tapering end show an order-of-magnitude
reduction in the latter.

Motivated by the above-mentioned works, the present study aims to provide a
framework for the design and implementation of ABH-embedded airfoils for the
modulation and mitigation of flexural waves. The organization of the remaining
sections in this paper is as follows. In Section II, the methods used in this
paper are presented, where, a wide range of ABH parameters are mapped inside
an airfoil profile (Sections II.1 and II.2) followed by the materials and
fabrication of representative geometries with the same total mass of
structural and damping material in Section II.3. A description of the
experimental setup where Laser Doppler vibrometry (LDV) is used to
characterize chord-wise vibrations when the airfoils are subjected to a
leading-edge point excitation in the 0.1-12 kHz range is provided in Section
II.4. Subsequently, FEA-based simulation and modeling is introduced in Section
II.5. The results of this work are presented and discussed in Section III,
beginning with an examination of the wavenumber-frequency characteristics in
Section III.1, followed by chord-frequency characteristics in Section III.2, a
detailed discussion of the FEA results in Section III.3 and the ABH airfoil
parametric study in Section III.4. The conclusions of this work are then
summarized in Section IV.

##  II Methods

###  II.1 ABH generating function

An ABH generating function is formulated with several parameter inputs (Figure
1). These have been chosen after compiling recent works on optimization of the
ABH shape or ABH-embedded shells and beams [19, 17, 16]. Horizontal ($x$) and
vertical ($y$) axes are defined along the length and thickness of the sample
respectively. The thickness, $h$, of the taper changes from its starting
maximum value, $h_{\text{s}}$, to its minimum truncated value, $h_{\text{t}}$,
following a taper power, $n$. The length over which $h$ = $h_{\text{s}}$
(constant) is denoted $L_{\text{s}}$, and the taper length over which
$h_{\text{s}}\geq h\geq h_{\text{t}}$ is denoted $L_{\text{t}}$. The total
length is set to a constant $L_{c}$, prescribing
$L_{\text{s}}$+$L_{\text{t}}\leq L_{c}$, where the equality and inequality
represent ‘continuous’ and ‘truncated’ ABHs, respectively. The damping layer
is distributed in the direction of increasing thickness from the end. The
total length of the damping is denoted $L_{\text{d}}$ and its height from the
taper is denoted $h_{\text{d}}$. When the ABH is truncated, there exists a
damping layer with a total height of $h_{\text{t}}+h_{\text{d}}$ in the
truncated region to ensure continuity of the exterior. The number of ABHs,
$N$, can also be varied by adopting a unit cell approach on the entire length.
A single taper is designated by $N$ = $1/2$, and whole numbers represent even
number of tapers. The range of values adopted for the different parameters is
listed in Table 1 with representative cases illustrated in Figure 2. Length
and thickness are normalized by $L_{c}$.

Figure 1: Schematic defining input parameters for the ABH generating function Table 1: Range of input parameters for the ABH generating function | $L_{\text{s}}$ | $L_{\text{t}}$ | $h_{t}/h_{\text{s}}$ | $n$ | $L_{\text{d}}$ | $h_{\text{d}}/h_{\text{s}}$ | $N$
---|---|---|---|---|---|---|---
Min | 0.0 | 0.8 | 0.1 | 2 | 0.1 | 0.0 | $1/2$
Step | 0.1 | 0.1 | 0.1 | 1 | 0.1 | 0.1 | $1/2$
Max | 0.2 | 1-$L_{s}$ | 0.3 | 5 | 0.5 | 0.5 | 5.0
Figure 2: Sample geometries from the ABH generating function for (a)
$L_{\text{s}}$=0, $L_{\text{t}}$=1, $h_{\text{t}}/h_{\text{s}}$=0.3, $n$=2,
$N$=$1/2$, $L_{\text{d}}$=0.3, $h_{\text{d}}/h_{\text{s}}$=0.3, (b)
$L_{\text{s}}$=0.2, $L_{\text{t}}$=0.6, $h_{\text{t}}/h_{\text{s}}$=0.1,
$n$=3, $N$=$1/2$, $L_{\text{d}}$=0.4, $h_{\text{d}}/h_{\text{s}}$ =0.4 and (c)
$L_{\text{s}}$=0.1, $L_{\text{t}}$=0.8, $h_{\text{t}}/h_{\text{s}}$=0.2,
$n$=4, $N$=2, $L_{\text{d}}$=0.5, $h_{\text{d}}/h_{\text{s}}$=0.4.

###  II.2 ABH-embedded airfoil geometry

Due to its prevalence in the aerospace community, the National Advisory
Committee for Aeronautics (NACA) system of 4-digit airfoil geometries is used.
Specifically, the symmetric NACA0012 foil is chosen in this study due to its
ubiquity, although the framework can be extended to any shape. Following
Ladson et al. [20], Equation (1) is used to generate the ordinates of a
NACA00tt foil, where $tt$ is the chord to thickness ratio.

$\frac{y(tt,x)}{0.05tt}=0.30\sqrt{x}-0.13x-0.35x^{2}+0.28x^{3}-0.10x^{4}$ (1)

The ABH profile obtained from the generating function is applied as a normal
offset curve [21] on the interior of the foil starting from the leading edge
(LE) to the trailing edge (TE) of the foil. The foil shape adds constraints on
the ABH-embedded foil in some cases. For instance, the foil chord length,
$L_{c}$ prescribes the maximum thickness of the foil, which is 0.12$L_{c}$ for
the NACA0012, thereby constraining $h_{d}$. ABH profiles that fall below the
$x$ axis after the offset curve calculation are bumped back to the axis. The
masses (i.e. areas under the curves) of the structural and damping materials
are then calculated and used to generate a baseline geometry where the
structural and damping profiles are at a constant (uniform) offset from that
of the foil external shape as illustrated in Figure 3.

Figure 3: (a) ABH generating function geometry, (b) ABH-embedded NACA0012
profile and (c) corresponding uniformly distributed baseline profile for
$L_{\text{s}}$=0.1, $L_{\text{t}}$=0.8, $h_{\text{t}}/h_{\text{s}}$=0.3,
$n$=2, $N$=$2$, $L_{\text{d}}$=0.5, $h_{\text{d}}/h_{\text{s}}$=0.5

Based on the parameter space (Table 1), a look-up-table (LUT) of ABH-embedded
foil designs is created. The target mass of structural and damping material
(arbitrarily chosen) is used to shortlist designs that have the same mass (or
within a small percentage). Note that the LUT shortlist sizes can be made
arbitrarily large by refining parameter increments (Table 1).

###  II.3 Materials and fabrication

The Stratasys J750 PolyJet printer is used for fabrication. It is capable of
building prints with multiple hard plastics, soft rubber-like materials as
well as a large number of composite materials by combining hard and soft
materials. Recently, Huang et al. [18] used a predecessor of this printer to
successfully fabricate and demonstrate the ABH effect in beams. Following
their selection, the hard plastic - VeroGray (RGD850) is chosen as the
structural material and the soft rubber-like TangoPlus (FLX930) is chosen as
the damping material for fabricating ABH-embedded foils. These are printed in
a single build using the ‘high mix’ build mode characterized by a build layer
resolution of 27 $\mu$m.

Although these materials have widely been used for various applications,
including the fabrication of ABH beams, their stiffness and loss properties
are not well established close to the frequency range of the current interest.
Huang et al. [18] perform tests around the 20 kHz range, very close to the
present range (0.1-12 kHz), however, due to the lack of available data and to
focus primarily on the influence of a Young’s modulus gradient on the wave
attenuation, they assume a constant loss factor of 0.1 for all their
materials. To better characterize the complex moduli of the 3D printed
materials for this study, an ad-hoc non-destructive testing (NDT) technique
using commercial grade compressional and shear wave transducers was developed
[22]. The resulting complex moduli, $E$ the Poisson’s ratio, $\nu$, as well as
the density, $\rho$ are presented in Table 2.

Table 2: Mechanical properties of 3D printed materials Material |  $E$ [GPa] |  $\rho$[g/cm3] | $\nu$
---|---|---|---
VeroGray | 2.5 | 1.16 | 0.35
TangoPlus | 0.65 + 0.4i | 1.18 | 0.47

To demonstrate the potential of ABH-embedded foils in modulating and
mitigating chordwise vibrations, four geometries are fabricated in the same
multi-material print job, with varying ABH-generating functions as shown in
Table 3. The selection was made to ensure a wide range in measured performance
while specifically evaluating the effect of truncation
($L_{\text{s}}$+$L_{\text{t}}\leq L_{\text{c}}$) and number of ABH elements
($N$). All cases have $h_{\text{t}}/h_{\text{s}}$ = 0.2 and $n$ = 2, i.e.
values that have been considered in other studies as well [14, 19]. It must
however be noted that this selection is not otherwise optimized. It serves as
a proof of concept that also validates subsequent FEA models (Sections II.5
and III.3) and measured material properties (Table 2), based on which a
parametric study is discussed at the end in Section III.4.

Table 3: 3D printed ABH-foil parameters # | $L_{\text{s}}$ | $L_{\text{t}}$ | $h_{\text{t}}/h_{\text{s}}$ | $n$ | $L_{\text{d}}$ | $h_{\text{d}}/h_{\text{s}}$ | $N$
---|---|---|---|---|---|---|---
1 | 0.10 | 0.90 | 0.20 | 2 | 0.50 | 0.50 | 1
2 | 0.08 | 0.87 | 0.20 | 2 | 0.50 | 0.50 | 3
3 | 0.12 | 0.87 | 0.20 | 2 | 0.50 | 0.50 | 1
4 | 0.03 | 0.97 | 0.20 | 2 | 0.50 | 0.54 | 3

All fabricated samples have the same total mass of structural and damping
material. A fifth baseline sample is also fabricated with the same masses
uniformly distributed following the outer shape of the airfoil, as if one were
to make a hollow damped version without any ABH-inspired tapers. This is a
crucial aspect that has often been overlooked in prior studies. For instance,
most studies compare the vibrational response of a beam or plate to that with
the so-called ABH version where substantial mass has been removed in machining
out the taper [13]. While Deng et al. [14] partly address this by applying the
same damping layer thickness for their baseline and ABH-embedded shells, as
also noted by them, the difference in structural masses evokes a different
response, complicating comparative performance assessment.

The airfoils have a chord, $L_{c}$ = 203.2 mm, resulting in a maximum
thickness of 0.12$L_{c}$ = 24.4 mm. To allow for clamping, a $L_{c}$/4 = 50.8
mm long section of constant thickness, $L_{c}/16$ = 12.7 mm is added upstream
of the leading edge. Foils are extruded to a depth of $L_{c}/8$ = 25.4 mm.
Images of the fabricated designs are shown in Figure 4. The five fabricated
samples weighed 86.2 g on average with a standard deviation of 0.2 g (0.2%).

Figure 4: Fabricated ABH-embedded foil designs. VeroGray and TangoPlus are
shaded in blue and pink respectively.

###  II.4 Experimental setup

The ABH-embedded foils are fixed as shown in Figure 5.

Figure 5: Schematic of the experimental setup

Vibrational excitation (B&K Type 4809) over a frequency range of 0.1 - 12 kHz
is provided through a stinger to the foil around its leading edge. A
piezoelectric force sensor (Model 208A11, 112 mV/N, PCB Piezotronics) is
connected between the stinger and foil, allowing an accurate measurement of
the force input by the exciter. The vibrational response, i.e. velocity along
the chord is measured using a single-point LDV (Polytec CLV-2534). The LDV is
mounted on a motorized translation stage (Velmex XSlide) and acquires data
over vertical increments of 0.1 mm. Velocity and signal strength from the LDV,
as well as force are sampled at a rate of 200 kHz using a National Instruments
compactRIO 9035 chassis equipped with analog (NI-9223) and digital (NI-9402)
input modules. More details of this setup, including its remote operation have
been described by Ikei [23].

A 200 ms long chirp excitation signal spanning a frequency range 0.1-12 kHz is
sent using a function generator (Agilent 33500B) to a power amplifier (B&K
Type 2718) and serves as the input to the exciter. At each sampling location,
LDV data from 32 time sequences is averaged. Given the large input excitation
frequency range, a quadratic convex chirp is amplitude-weighted towards higher
frequencies as shown in Figure 6.

Figure 6: Spectrogram of the excitation input waveform

The time-averaged force response at the point of excitation is shown in Figure
7 for all five samples.

Figure 7: Frequency spectrum of the measured force input.

The horizontal and vertical axes are plotted in logarithmic scale denoting
frequency, $f$ in kHz and force, $F$ in mN respectively. For convenience,
frequency grid lines are chosen near local extrema. Evidently, the same modal
signatures, and profiles are found across the entire range for all the cases.
The frequency spectra of the LDV data (not shown) calculated at excitation
($x/L_{c}$=0) are also consistent with the trends in Figure 7 for all cases.
The forces are as high as 2-3 N around 100-200, 250-400 and 2830 Hz. Minima
around 220, 485, and 10445 Hz have forces in the 0.1-0.2 N range, that are
presumably anti-resonances associated with structural modes common to all the
samples. Consequently, the measured force and velocity (from LDV) at these
frequencies are extremely small, comparable to the noise level. However, on
either sides of these minima, the force recovers to a value above 1 N.
Therefore, it can be concluded that the amplitude weighting is for the most
part, effective at preventing any significant decay in the forcing with an
increase in frequency, allowing maximum utilization of the LDV’s dynamic
range. For subsequent analysis, frequency spectra are normalized by the value
at excitation, allowing direct inter-frequency and inter-sample comparisons to
be drawn.

As evident from Figure 5, the foil surface is curved in the direction of the
LDV laser beam leading to variations in the optical path length ($\Delta$OPL)
of 0.09$L_{c}$ = 18 mm. The LDV is positioned such that its nearest visibility
maximum is centered within $\Delta$OPL leading to a 9% variation around the
peaks that are 204 mm apart. Hence, this is not expected to have a substantial
impact on the measurements. Furthermore, the LDV signal strength is also
acquired for every measurement, based on which outliers in the data are
identified and flagged for subsequent processing steps. The raw LDV data is
analyzed in the frequency-wavenumber ($f$-$k$) space to identify bounds for
most of the energy content. Subsequently, a Tukey window is applied on the
data to remove noise associated with high wavenumbers. A threshold condition
of $|k|/2\pi<20/L_{c}$ is used for all the cases. For convenience, Table 4
enlists relevant spectral parameters for space and time.

Table 4: Relevant space and time spectrum parameters | Samples | Resolution | Max range | Range of interest
---|---|---|---|---
$t$ | 40,000 | 5 $\mu$s | 0 - 100 kHz | 0.1 - 12 kHz
$x$ | 3,001 | 0.1 mm | 0 - 5 km-1 | 0 - 618 m-1

###  II.5 Modeling

Three-dimensional finite element analysis using COMSOL Multiphysics is
performed on the five measured cases. In practice, most airfoil wings are thin
and hollow, lending themselves to a 2D plate approximation, making FEA
substantially less computationally intensive. The present samples also are
thin and plate-like in the tapering part of the geometry ($x/L_{c}\geq$ 0.03 -
0.12, refer Table 3). However, clamping and excitation requirements force a
beam-like structure upstream ($x/L_{c}\leq$ 0) making the present structures
complex and requiring a 3D model for accurate results. The CAD model used for
fabricating the foils is directly imported into COMSOL for analysis. A fixed
boundary condition is applied at the upstream edge (Figure 5). A force of 1 N,
based on the force measurements (Figure 7) is distributed where the force
sensor mounts to the foil. The computational domain is halved by leveraging
symmetry in the direction of extrusion to reduce solver time. Material
properties as per Table 2 are applied to the model. A frequency domain
simulation spanning the 0.1 - 12 kHz range (Table 4) with a resolution of 100
frequencies per decade is executed.

The mesh resolution is a critical component affecting the accuracy of the
model. Based on the sample geometry, an upper bound of $L_{c}$/32 = 6.35 mm is
obtained by requiring at least two mesh elements to span the $h$ =
$h_{\text{s}}$ region (Figure 5). A lower bound can obtained by requiring at
least two elements near the truncated edge, resulting in $h_{\text{t}}/2$ =
0.5 mm, which also corresponds to approximately 20 elements per wavelength at
the largest wavenumber of interest (Table 4). These constraints are imposed on
the mesh followed by refinement of the element growth rate and curvature until
there is less than a 0.5% change in the computed velocity over all frequencies
in the domain. Figure 8 shows a 2D slice of the mesh elements (a) over the
entire domain, as well as a sub-region of the mesh (b) focusing near the
truncated edge for Case 1. The tetrahedral volumetric mesh has roughly 200,000
elements with variations between cases dictated by the local geometry.

Figure 8: FEA mesh elements for (a) entire domain and (b) sub region near
truncation for Case 1. All dimensions in m. Figure 9: Contour plots of 20
log($\hat{V}(k,f)|_{x>0}$) for all the measured cases overlaid (black solid
curves) with Mindlin corrected RKU composite beam model prediction for the
uniformly distributed case, $k_{\text{uni}}(f)$.

##  III Results and Discussion

###  III.1 Wavenumber-frequency characteristics

The time Fourier spectrum of the wavenumber-filtered velocity data is computed
and as noted earlier, normalized by the corresponding value at $x/L_{c}$ = 0
over the entire frequency range, denoted as $\hat{V}(x,f)=V(x,f)/V(0,f)$.
Distributions of the wavenumber spectrum, $\hat{V}(k,f)$ restricting to
$x/L_{c}\geq$ 0 are shown for all cases in Figure 9. The part of the
excitation signal that goes in the other direction ($x/L_{c}<$ 0) couples to
the mounting structure (Figure 5). Evidently, this region has the same effect
on all the samples (Figure 7) and is thereby excluded from subsequent analysis
that aims to characterize relative effects of ABH tapers on the foil.

The Ross-Kerwin-Ungar (RKU) method is used to model the effective bending
stiffness of the uniformly distributed baseline foil, assuming the TangoPlus
to be a constant and thin absorbing layer on the VeroGray [24, 12]. The
Mindlin plate correction [25] is applied to this stiffness to derive the
composite wavenumber for the uniform case, denoted by $k_{\text{uni}}(f)$,
also shown in Figure 9. Evidently, $k_{\text{uni}}(f)$ provides a very good
prediction for the measured data over all the cases. Following convention,
$k>0$ refers to waves traveling from the exciter (LE) towards the TE (Figure
5), while $k<0$ contains waves reflected back from TE to LE.

The uniformly distributed case illustrated in Figure 9(a) has high amplitude
densely concentrated around $k_{\text{uni}}(f)$ for the entire $k>0$ range
over all frequencies. This is presumably because it has a constant thickness
throughout the foil. Conversely, in the four ABH foils, tapers and damping
layers originate and terminate at various chordwise locations, resulting in a
variable local bending stiffness. Hence, as expected, the amplitude of
$\hat{V}(k,f)$ is distributed over a larger region about $k_{\text{uni}}(f)$.
In several regions, alternative wave-paths become available. This ‘smearing’
of $\hat{V}(k,f)$ is also in agreement with other works [24, 26].

Another important conclusion from Figure 9 is that across all the cases,
reflected waves ($k<0$) symmetric to the incident waves exist only below 1.8
kHz. In other words, there is a standing wave below 1.8 kHz, corresponding to
$kL_{c}/2\pi$ = 2.75 on the curve, or a wavelength of $L_{c}/17$, which is
equivalent to the starting thickness, $L_{c}/16$ of the samples (Figure 5).
This indicates that past 1.8 kHz, waves enter all the samples, and most of
them do not make it back, based on the substantially higher magnitude of the
$k>0$ waves compared to those that are in the opposite direction. Therefore,
subsequent discussions, where the objective remains to study the modulation of
flexural characteristics in the foil regions by ABHs, are restricted to
frequencies above 1.8 kHz.

###  III.2 Chord-frequency characteristics

Contour plots of $\hat{V}(x,f)$ in dB for all the cases on top of their
corresponding geometries (not to scale) are shown in Figure 10.

Figure 10: LDV frequency-chordwise distributions of $\hat{V}(x,f)$ in dB for
all cases. Corresponding foil geometries (not to scale) below. Region
0-0.2$L_{c}$ close to the trailing edge (TE) across all frequencies is
outlined in yellow for emphasis.

Chordwise locations with geometric junctions, i.e. starts and ends of tapers
and damping are shown as grid lines in the bottom row to facilitate direct
correlations with spatial velocity distributions. To emphasize the differences
in the trailing edge vibrations between cases, the region within 0.2$L_{c}$ of
the trailing edge is outlined in yellow. Evidently, $\hat{V}(x,f)$ varies
substantially with $x/L_{c}$ across samples. Magnitudes are elevated near the
leading and trailing edges for all cases, presumably due to the absence of
damping. Lower amplitudes prevail in the mid-chord, where a series of trapped
waves present as vertical bands. Some examples of this inter-junction trapping
can be seen in (i) all cases above 8 kHz about their first junction, and (ii)
Cases 1 and 3 in the 6-9 kHz region between the 3rd and 4th junctions.

To quantify the extent of spatial modulation across frequencies,
$\hat{V}_{\text{rms}}(x)$ is computed by integrating across $f$. To facilitate
comparisons, profiles of $\hat{V}_{\text{rms}}(x)$ corresponding to the
baseline case are then subtracted from all the others, denoted as
$\hat{V}_{\text{rms}}^{\Delta}(x)$ and shown in Figure 11 expressed in dB for
the (a) 1.8 - 10 kHz and (b) 0.1 - 1.8 kHz ranges.

Figure 11: Profiles of $\hat{V}_{\text{rms}}^{\Delta}(x)$ in dB integrated
over (a) 1.8 - 10 kHz and (b) 0.1 - 1.8 kHz for all cases. Figure 12: Baseline
subtracted profiles of $\hat{V}(x,f)$ in dB averaged over the (a) first
(0$<x/L_{c}<$0.5) and (b) second (0.5$<x/L_{c}<$1) halves of the chord.

Although frequencies below 1.8 kHz may not interact with the ABH structures
(as discussed in Section III.1), it is still interesting to note that despite
everything else being similar, there is still an effect of the ABHs seen in
Figure 11(b), where the $N$=3 and truncated cases have substantially lower
magnitudes compared to those of the baseline, $N$=1 or continuous cases
respectively. Subsequent discussion is restricted to the 1.8 - 10 kHz range.

Figure 13: FEA frequency-chordwise distributions of $\hat{V}(x,f)$ in dB for
all cases. Corresponding foil geometries (not to scale) below. Region
0-0.2$L_{c}$ close to the trailing edge (TE) across all frequencies is
outlined in yellow for emphasis.

For $N$=1 geometries, i.e. Cases 1 and 3, damping only starts at
$x/L_{c}$=0.25, thereby extending the extent of the elevated amplitude region
near the leading edge. However, it is interesting to note that although
damping starts at $x/L_{c}$=0.02 for the baseline, compared to $x/L_{c}$=0.08
for the $N$=3 geometries (Cases 2 and 4), elevated amplitude regions extend
till $x/L_{c}$ = 0.16 in the 1.8-6 kHz range for all of them. Therefore,
despite the absence of damping, in this frequency range, ABH tapers perform
better in the 0.08 $<x/L_{c}<$ 0.16 region for Cases 2 and 4 when compared to
the baseline. This trend reverses above 6 kHz, and also when integrated over
the entire frequency range. The baseline case performs better than all the
current cases (Figure 11) near the leading edge, where the earlier onset of
damping outweighs all other factors.

In the mid-chord region (0.2 $<x/L_{c}<$ 0.8), there is a clear distinction
between different cases. $N$=3 geometries substantially outperform the
baseline case throughout, in fact the reduction in amplitude is above 6-7 dB
at $x/L_{c}$ = 0.25, 0.35, 0.65 and 0.7 as seen in Figure 10 and 11). $N$=1
geometries perform consistently worse than the baseline by as much as 5-7 dB
at $x/L_{c}$ = 0.2-0.35. The only location where $N$=1 geometries fair better
than the baseline, by 2-3 dB, is around $x/L_{c}$=0.6 and 0.7.

Near the trailing edge (0.95 $<x/L_{c}<$ 1), $N$=1 geometries perform very
similar to the baseline. $\hat{V}_{\text{rms}}(x)$ as well as the spatial
distributions along the entire frequency range (Figure 10) are very similar.
However, $N$=3 geometries showcase substantial reduction in magnitude, as high
as 10 dB for the truncated Case 2, when integrated over the entire frequency
range (Figure 11). This result can have significant implications for airfoil
design where trail-edge noise control is of interest. The spatial velocity
distributions (Figure 10) reveal that there is not only a broadband reduction
in amplitude near the trailing edge, but also a down-shift in the 4 kHz
baseline peak to around 3.3 kHz.

To further quantify the frequency modulation by the ABHs, Figure 12 shows the
baseline subtracted $\hat{V}(x,f)$ in dB averaged in the (a) first
(0$<x/L_{c}<$0.5) and (b) second (0.5$<x/L_{c}<$1) halves of the chord. As
discussed previously, there is a front-loading effect prevalent in the $N$=1
cases, that is also highlighted by Figure 12(a), with a 5 dB increase in
velocity above 6 kHz. However, for the second half of the foil (Figure 12(b))
in this same range, Cases 1 and 3 exhibit a 3 dB reduction compared to the
baseline. Therefore, the same frequency range is modulated 5 dB above and 3 dB
below the baseline in the first and second halves respectively by the $N$=1
geometries. The $N$=3 cases 2 and 4 show 3-5 dB reduction in the first half of
the foil for some frequency bands, i.e. 4-5 and 8-10 kHz while performing
similar to the baseline in other regions of the leading edge half. However, in
the second half (Figure 12), they exhibit a 5 dB reduction in amplitude across
4-10 kHz range. Such a large broadband reduction in amplitude across the
entire second half of the airfoil can have significant implications to
applications involving control of flow separation, stall and other
transitional and unsteady effects.

###  III.3 FEA results

Distributions of $\hat{V}(x,f)$ in dB obtained from the FEA simulations for
all the cases along with their corresponding geometries are shown in Figure
13. The plot extents, color map range, grid and placement are identical to
those of Figure 10 to facilitate a comparison of the simulations with the LDV
measurements. There is excellent agreement between the FEA model and the LDV
measurements. They capture the spatial extents of the elevated amplitude
regions near the leading and trailing edges, as well as the differences
between cases across the entire frequency range. For e.g. the 4 kHz peak and
neighboring features near the leading edge for the baseline and Cases 1 and 3,
albeit shifted slightly, match very well in the FEA results. The trailing edge
distributions here also show the close similarities between the baseline and
$N$=1 geometries, whereas the $N$=3 cases also show a significant broadband
reduction in amplitude, in good agreement with the LDV data. Furthermore, the
inter-junction trapping examples (i)-(ii) highlighted for the LDV data are
also captured extremely well by the FEA models. The frequency down-shift of
the trailing edge 4 kHz peak from the baseline to the $N$=3 cases is also
evident.

There are some minor discrepancies that are worth mentioning. First, as eluded
to earlier, the entire data, across all samples appears slightly downshifted
in frequency when compared to the measurements. This suggests that the printed
samples had slightly different material properties than those used in the FEA
model. This difference is small enough that there doesn’t seem to be any
significant change in vibrational modes. Second, not all the features in the
mid-chord, especially above 5 kHz for the baseline case are captured by the
FEA. Third, Cases 2 and 4 appear to have elevated amplitudes in the FEA
results (Figure 13(c) and (e)) above 8 kHz. This might be linked to the
frequency down-shift described earlier, which carries through the entire
range. Since this is a wavenumber-effect, the shift is expected to increase
with increasing frequency. Thereby, the elevated amplitudes around 9 kHz would
be expected around 10-11 kHz in the LDV data. As evident from Figure 7, the
measured value drops significantly in this range until at least 12 kHz, where
the measured LDV signal has a low signal-to-noise ratio, affecting the
frequency normalization and precluding any conclusions to be drawn in this
range. Despite these issues, the FEA model does a really good job at capturing
the key features of the velocity distributions and most importantly, the
differences between geometries. Therefore, for a given desired objective
function, this can help reduce the fabrication requirement substantially.

Furthermore, the FEA model can be used to visualize the ABH effect by looking
at frequency-specific displacement fields for all cases. Figure 14 shows a
snapshot of the displacement, amplified by 30,000 times, at 4 kHz for all the
samples obtained from the FEA models.

Figure 14: FEA displacements at 4 kHz for all samples amplified by 30,000
times.

The leading and trailing edges of the $N$=1 cases deflect even higher than the
baseline. However, $N$=3 geometries, as seen in all the data, have
substantially lower amplitudes throughout the foils, especially at the
trailing edge, elucidating the overall impact of the ABHs.

###  III.4 ABH-foil parametric study

Although a comprehensive optimization for various cost functions and
parametric spreads remains out of the present scope, a glimpse of the spatial
vibrational modulation that can be achieved by implementing ABHs in airfoils
is provided here. The validation of the FEA models also confirms the material
properties moving forward. As discussed earlier, FEA models were required to
be 3D to match better to the fabricated samples. However, given that most
applications of foils are expected to have a large span-to-chord ratio, a 2D
plate approximation may be sufficient. This also makes FEA computation for the
entire shortlisted LUT feasible, which is automated using COMSOL Livelink with
MATLAB.

For the present purpose, 24 geometries with masses within 2% of each other are
shortlisted from the LUT (Table 1). For brevity, we restrict only to whole-
numbered $N$ values, i.e. the geometries terminate with VeroGray. As noted
earlier, the LUT and shortlist thresholds can be made arbitrarily dense to
accommodate the objective at hand. Figure 15 contains profiles of
$\hat{V}_{\text{rms}}^{\Delta}(x)$ in dB for all 24 geometries.

Figure 15: ABH foil parametric study: 2D FEA-generated profiles of
$\hat{V}_{\text{rms}}^{\Delta}(x)$ in dB for shortlisted cases.

They are colored based on their $N$ values, while continuous and truncated
geometries are represented as solid and dotted lines respectively. Curves also
represent variations in $n$, $L_{\text{s}}$, $L_{\text{t}}$, $L_{\text{d}}$,
$h_{\text{s}}$ and $h_{\text{d}}$ while satisfying the mass constraints, and
these have not been explicitly identified to restrict to the present scope.
Evidently, the profiles vary as high as 20 dB compared to the baseline case.
As noted in the fabricated cases earlier, even here, it appears that
increasing $N$ reduces the overall vibration level, especially in the second
half of the chord, closer to the trailing edge. This is presumably due to a
compounding effect as the waves encounter more damped ABHs before reaching the
trailing edge. Conversely, $N$=1 cases, have the most elevated vibration
levels, especially in the first half of the chord. Truncation also reduces the
amplitude, again, more so, closer to the trailing edge.

##  IV  Conclusions

This study introduces a framework for the design and implementation of ABHs in
airfoils. A multi-parameter damped-ABH generation function is mapped onto a
NACA series airfoil. Four ABH geometries and a uniformly distributed baseline,
all with the same mass of structure and damping are fabricated using multi-
material PolyJet 3D printing. Laser Doppler vibrometer measurements of
velocity along the airfoil chord in response to a broadband 0.1 - 12 kHz
excitation are performed for all the cases. 3D FEA is also performed on the
fabricated geometries, to enable for model and material property validation.
Furthermore, a parametric 2D FEA study is performed on shortlisted geometries
using the validated material properties, to showcase the mitigation and
modulation that is achievable by implementing ABH design. Key findings of the
study are described below,

  * •

Wavenumber-frequency characteristics of the measured data follow the Mindlin-
corrected RKU model. The uniform baseline is densely concentrated about the
curves, whereas all the ABH cases show $k$-space smearing of energy in
agreement with findings from other ABH studies [26].

  * •

In general, spatial distributions of velocity as a function of frequency,
normalized by that at excitation reveal substantial variations between
samples. Magnitudes are elevated near the leading and trailing edges of the
foils, while lower amplitudes prevail in the middle. In the ABH cases, a
series of standing waves are trapped between local junctions where tapers and
damping transition.

  * •

In comparison to the uniform baseline, there is an reduction of 5 dB in the
magnitude across the entire frequency range for foils with $N$=3 embedded ABHs
on average over the chord length, with up to a 10 dB reduction near the
trailing edge for the truncated case. On the other hand, ABH foils with $N$=1
ABHs are associated with an increase in magnitude by as much as 5-7 dB in the
first half of the chord, while remaining comparable to the baseline towards in
the second half.

  * •

Baseline-subtracted velocity profiles averaged in the first half of the chord
show that the front-loading effect for $N$=1 ABH cases exists above 6 kHz.
Profiles in the second half elucidate a broadband (4-10 kHz) reduction in
amplitude by 5 dB with the $N$=3 cases.

  * •

3D finite element models of the five samples are in good agreement with LDV
measurements. They capture the spatial extents of the leading and trailing
edges as well as inter-sample variations very well. These validate the model
as well as the material properties.

  * •

Two-dimensional parametric FEA results indicate that a modulation in the
velocity amplitude of as much as 20 dB with ABH-embedded foil designs. The
effects of the number of ABHs and truncation are in agreement with the trends
measured in the experiments.

In conclusion, this study provides an insight into the design, fabrication,
performance, modeling and optimization of ABH-embedded airfoils. Given a
constant mass structure, airfoils can be designed to mitigate, focus or
modulate vibrations for any chordwise region by adapting the process presented
in this study. In applications where minimizing the noise radiation from
trailing edges is of concern, findings from the present study can be used to
achieve upwards of a 10 dB reduction in vibration. For cases where minimizing
broadband vibrations for structural integrity and wear of foils is important,
multiple truncated ABHs can be distributed in a spanwise orientation.
Alternatively, for achieving flow control at specific chordwise locations or
frequency bands, using the present framework, energy can be added or
subtracted as desired from the boundary layer close to the foil leading to
superior performance and efficiency. Furthermore, other applications involving
energy harvesting or restructuring can benefit from the front-loading effect
introduced by the $N$=1 ABH cases.

###### Acknowledgements.

This work was supported by the Office of Naval Research.

## References

  * [1] Michel Roger and Stéphane Moreau.  Broadband self noise from loaded fan blades.  AIAA journal, 42(3):536–544, 2004.  ISBN: 0001-1452.
  * [2] Thomas F. Brooks, D. Stuart Pope, and Michael A. Marcolini.  Airfoil self-noise and prediction, volume 1218.  National Aeronautics and Space Administration, Office of Management …, 1989.
  * [3] T.F. Brooks and T.H. Hodgson.  Trailing edge noise prediction from measured surface pressures.  Journal of Sound and Vibration, 78(1):69–117, September 1981.
  * [4] M.S. Howe.  A review of the theory of trailing edge noise.  Journal of Sound and Vibration, 61(3):437–465, December 1978.
  * [5] Mathieu Gruber.  Airfoil noise reduction by edge treatments.  PhD thesis, University of Southampton, 2012.
  * [6] M. A. Mironov.  Propagation of a flexural wave in a plate whose thickness decreases smoothly to zero in a finite interval.  Soviet Physics Acoustics-USSR, 34(3):318–319, 1988.  Issue: 3 Pages: 318-319 Publication Title:.
  * [7] V.V. Krylov and R.E.T.B. Winward.  Experimental investigation of the acoustic black hole effect for flexural waves in tapered plates.  Journal of Sound and Vibration, 300(1-2):43–49, February 2007.
  * [8] D.J. O’Boy, V.V. Krylov, and V. Kralovic.  Damping of flexural vibrations in rectangular plates using the acoustic black hole effect.  Journal of Sound and Vibration, 329(22):4672–4688, October 2010\.
  * [9] Victor V. Krylov.  Acoustic black holes: recent developments in the theory and applications.  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 61(8):1296–1306, August 2014.
  * [10] Philip A. Feurtado and Stephen C. Conlon.  An Experimental Investigation of Acoustic Black Hole Dynamics at Low, Mid, and High Frequencies.  Journal of Vibration and Acoustics, 138(6):061002, December 2016\.
  * [11] V.V. Krylov and F.J.B.S. Tilman.  Acoustic ‘black holes’ for flexural waves as effective vibration dampers.  Journal of Sound and Vibration, 274(3):605–619, July 2004.
  * [12] Alfonso Climente, Daniel Torrent, and José Sánchez-Dehesa.  Omnidirectional broadband insulating device for flexural waves in thin plates.  Journal of Applied Physics, 114(21):214903, December 2013.
  * [13] Adrien Pelat, François Gautier, Stephen C. Conlon, and Fabio Semperlotti.  The acoustic black hole: A review of theory and applications.  Journal of Sound and Vibration, page 115316, March 2020.
  * [14] Jie Deng, Oriol Guasch, Laurent Maxit, and Ling Zheng.  Reduction of Bloch-Floquet bending waves via annular acoustic black holes in periodically supported cylindrical shell structures.  Applied Acoustics, 169:107424, December 2020.
  * [15] E.P. Bowyer and V.V. Krylov.  Damping of flexural vibrations in turbofan blades using the acoustic black hole effect.  Applied Acoustics, 76:359–365, February 2014.
  * [16] Jie Deng, Oriol Guasch, Laurent Maxit, and Ling Zheng.  Transmission loss of plates with multiple embedded acoustic black holes using statistical modal energy distribution analysis.  Mechanical Systems and Signal Processing, 150:107262, March 2021\.
  * [17] Jie Deng, Oriol Guasch, and Ling Zheng.  A semi-analytical method for characterizing vibrations in circular beams with embedded acoustic black holes.  Journal of Sound and Vibration, 476:115307, June 2020.
  * [18] Wei Huang, Hui Zhang, Daniel J. Inman, Jinhao Qiu, Carlos E.S. Cesnik, and Hongli Ji.  Low reflection effect by 3D printed functionally graded acoustic black holes.  Journal of Sound and Vibration, 450:96–108, June 2019.
  * [19] Cameron A. McCormick and Micah R. Shepherd.  Design optimization and performance comparison of three styles of one-dimensional acoustic black hole vibration absorbers.  Journal of Sound and Vibration, 470:115164, March 2020.
  * [20] Charles L. Ladson, Cuyler W. Brooks Jr, Acquilla S. Hill, and Darrell W. Sproles.  Computer program to obtain ordinates for NACA airfoils.  page 27, 1996.
  * [21] Joss Duchateau.  Offset Curve, 2020.  Library Catalog: www.mathworks.com.
  * [22] Matthew D. Guild, Kaushik Sampath, and Caleb F Sieck.  Broadband ultrasonic characterization of the complex-valued elastic properties for 3D printed polymers.  In Preparation, 2021.
  * [23] Alec K. Ikei and Kaushik Sampath.  Remote Operation of a Single-Point LDV System to Acquire 2D Measurements.  arXiv:2102.13031 [cond-mat], February 2021.  arXiv: 2102.13031.
  * [24] Julien Leng.  Controlling flexural waves using subwavelength perfect absorbers : application to Acoustic Black Holes.  phdthesis, Université du Maine, November 2019.
  * [25] Andrew N. Norris.  Flexural waves on narrow plates.  The Journal of the Acoustical Society of America, 113(5):2647–2658, May 2003.  Publisher: Acoustical Society of America.
  * [26] Philip Andrew Feurtado.  Quiet Structure Design Using Acoustic Black Holes.  Ph.D., The Pennsylvania State University, United States – Pennsylvania, 2017.  ISBN: 9781369991451 Publication Title: ProQuest Dissertations and Theses Section: 0176.


††thanks: ∗Corresponding author<EMAIL_ADDRESS>
# Breaking up with the continuous exoplanet mass-radius relation

Kathryn Edmondson  Jordan Norris  Eamonn Kerins∗ Department of Physics and
Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK.

###### Abstract

We use a carefully selected sub-sample of 1053 confirmed exoplanets from the
NASA Exoplanet Archive to construct empirical power-law exoplanet mass-radius-
temperature ($M$-$R$-$T$) relations. Using orthogonal distance regression to
account for errors in both mass and radius, we allow the data to decide: 1)
the number of distinct planetary regimes; 2) whether the boundaries of these
regimes are best described by broken power laws joined at mass break points,
or by discontinuous power laws motivated by changes in equations of state and
temperature. We find strong support from the data for three distinct planetary
$M$-$R$ regimes and for those regimes to be discontinuous. Our most successful
model involves an $M$-$R$-$T$ relation in which ice/rock (rocky) and ice-giant
(neptunian) planets are segregated by a pure-ice equation of state, whilst
neptunes and gas giant (jovian) planets are segregated by a mass break at
$M_{\rm br}=115\pm 19~{}M_{\oplus}$. The rocky planet regime is shown to
follow $M\propto R^{0.34\pm 0.01}$, whilst neptunes have $M\propto R^{0.55\pm
0.02}$. Planets in both regimes are seen to extend to similar maximum masses.
In the jovian regime we find that $M\propto R^{0.00\pm 0.01}T^{0.35\pm 0.02}$,
where $T$ is the planet equilibrium temperature. This implies that, for jovian
planets detected so far, equilibrium temperature alone provides a robust
estimator of mass.

##  1 Introduction

Understanding the relationship between exoplanet mass and radius is crucial to
constraining internal composition and testing planet formation simulations, as
well as predicting the detectability of a planet to aid in the design of
future surveys. Some studies (e.g. Seager et al., 2007; Swift et al., 2012)
take the physically-motivated approach of modelling planetary composition
using equations of state to infer a mass-radius relation. Others take an
empirical approach, involving applying analytic, probabilistic, or machine-
learning models to data (e.g. Otegi et al., 2020; Chen & Kipping, 2017;
Mousavi-Sadr et al., 2023). Weiss et al. (2013) defined two planetary regimes
either side of $150~{}M_{\oplus}$, and fit a mass-radius-incident-flux
($M$-$R$-$F$) plane to each population. The small and large planet regimes
were found to follow $R\propto M^{0.53}F^{-0.03}$ and $R\propto
M^{0.039}F^{0.094}$ respectively. Wolfgang et al. (2016) focused on small
planets less than $4R_{\oplus}$ and used a Hierarchical Bayesian model to
obtain a probabilistic model with $M\propto R^{1.3}$. A later relation
developed by Bashi et al. (2017) again split exoplanets into two regimes, this
time with a floating break-point. This yielded $R\propto M^{0.55}$ and
$R\propto M^{0.01}$ for the small and large planet regimes, and a break-point
at $127M_{\oplus}$, which was attributed to the mass at which electrons in
hydrogen become degenerate.

Chen & Kipping (2017) used a similar methodology to Wolfgang et al. (2016),
extending their analysis to develop a probablistic model, Forecaster, that
also included brown dwarfs, small stars, and solar system dwarf planets and
moons. They identified four regimes: terran worlds, for which $R\propto
M^{0.28}$; neptunian worlds, for which $R\propto M^{0.59}$; jovian worlds, for
which $R\propto M^{-0.04}$ and stellar worlds, for which $R\propto M^{0.88}$.
The transition between terran and neptunian worlds occurred at $2M_{\oplus}$,
and between neptunian and jovian worlds at $130M_{\oplus}$.

Power laws are often used for exoplanet mass-radius ($M$-$R$) relations.
However, Ning et al. (2018) argued that such a model is not adequately
flexible to describe more complex features in a $M$-$R$ diagram and, hence,
they presented a non-parametric model through the use of Bernstein
polynomials.

Ulmer-Moll et al. (2019) investigated the dependence of radius on other
parameters, such as equilibrium temperature, semi-major axis and properties of
the host star through a machine-learning approach. This method circumvents the
need to classify planets, or use differing $M$-$R$ relations in different
regimes, and was better able to characterise the spread of hot-jupiter radii
than Chen & Kipping (2017).

All of the models discussed so far imposed the condition that a $M$-$R$
relation should be continuous, however Otegi et al. (2020) took a different
approach by categorising exoplanets below $120M_{\oplus}$ as rocky or
volatile-rich according to an equation of state of pure water. Their results
were unaffected by the exact equation of state and temperature assumption
used, yielding $R\propto M^{0.29}$ for rocky planets and $R\propto M^{0.63}$
for those rich in volatiles. The main downside of a discontinuous $M$-$R$
relation is that it permits relations that overlap in mass or radius,
resulting in non-uniqueness. But, this may allow a more accurate
characterisation of underlying discontinuities that arise as a result of
planet compositional transitions.

The large scatter in the radii of massive planets cannot be explained simply
by a deterministic power law and is generally attributed to atmospheric
bloating due to stellar insolation. Enoch et al. (2012) investigated which
parameters contribute to this scatter for Saturn-, Jupiter- and high-mass
planets and found that the radius of a Jupiter-mass planet could be predicted
from equilibrium temperature and semi-major axis, with no mass dependence at
all. Saturn- and high-mass planet radii were found to depend also on stellar
metallicity and planet mass, however the division of exoplanets into these
three populations is somewhat arbitrary. Thorngren & Fortney (2018) and
Sestovic et al. (2018) both used hierarchical Bayesian modelling to
investigate the cause of bloating and concluded that insolation is key in
characterising inflated radii, favouring a $M$-$R$-$F$ relation. Thorngren &
Fortney (2018) also found that radius inflation decreased for high enough
temperatures, deducing that the mechanism by which bloating occurs is Ohmic
dissipation, in which a planet’s magnetic field interacts with atmospheric
flows, thereby transferring heat to the planet interior. Ohmic dissipation
initially has a greater impact with increasing equilibrium temperature and
thus ionisation, however for very high temperatures the atmospheric winds are
subject to magnetic drag and the process becomes inhibited (Batygin et al.,
2011).

There have been several studies to investigate whether the composition of
rocky planets can be described by the abundances of rock-building elements in
the host star. Schulze et al. (2021) found that in general, the mass and
radius of a rocky exoplanet is consistent with a model of the planet derived
from the $\frac{\mathrm{Fe}}{\mathrm{Mg}}$ and
$\frac{\mathrm{Si}}{\mathrm{Mg}}$ ratios found in the host star’s photosphere.
However, there are mechanisms after planet formation that can alter planet
composition. For example, Mercury has been vastly enriched in iron,
potentially due to mantle-stripping collisions.

In this paper we seek to take advantage of the rapid expansion of exoplanet
data to revisit the exoplanet $M$-$R$ relation. By using a carefully selected
data subsample, we allow the data alone to decide on the required number and
location of the different planetary regimes. We also test for the support
between continuous and discontinuous $M$-$R$ relations, including extension to
a $M$-$R$-$T$ relation for massive planets. The paper is organized as follows.
In Section 2 we define our data subsample and in section 3 we use the data to
construct piece-wise power laws, allowing for floating planet mass break
points and accounting for errors in both mass and radius. In section 4, we
explore massive-planet $M$-$R$-$T$ relations as well as discontinuous $M$-$R$
relations for lower-mass planets. In Section 5 we compare the performance of
our fits to each other and to the widely-used probabilistic Forecaster code
(Chen & Kipping, 2017). In Section 6 we present our conclusions.

##  2 Data

Exoplanet data were retrieved from the NASA Exoplanet
Archive111https://exoplanetarchive.ipac.caltech.edu/, accessed during
September 2022, consisting at the time of 5171 confirmed exoplanets. Of these,
1057 have both mass and radius measurements quoted with uncertainties. In
cases where a planet had multiple measurements of the same parameter, we chose
the entry that optimized the fractional uncertainty in density. In cases where
equilibrium temperature measurements were available, the most precise
measurement was taken. While previous works have implemented significance cuts
that require the mass and radius to exceed a given multiple of their errors,
we do not consider such a selection cut here due to its potential for bias in
the small planet regime (Burt et al., 2018).

We impose two plausibility criteria: firstly, a planet must have a mass less
than the minimum mass for the onset of deuterium burning; secondly, a planet
must have a density less than that of pure iron. In the first criterion, there
is some ambiguity as to the exact mass above which deuterium burning occurs.
Spiegel et al. (2011) have demonstrated that this mass is not well defined in
general, due to the influence of several factors including helium abundance,
initial deuterium abundance and metallicity; leading to deuterium burning
limits between around $11M_{J}$ and $16M_{J}$, depending on model assumptions.
For our sample selection we adopt the canonical value of $13M_{J}$.

Some planets in the sample give a bulk density larger than that of pure iron.
While measurement error may be a cause, it has been suggested by Mocquet et
al. (2014) that there could be a physical explanation for these anomalous
cases. Mocquet et al. (2014) demonstrated that the bare core of a gas giant
which has been stripped of its atmosphere during migration could have been
irreversibly compressed, leading to a new regime of very high density
exoplanets. However, there has been no follow-up work to date and we consider
the current sample of very high density exoplanets to be too small to warrant
their treatment as a new sub-population of planets. Consequently, super-iron
density planets are excluded by our plausibility criteria.

The radius of a pure-iron planet, $R_{\rm Fe}$, has been obtained from Fortney
et al. (2007) by setting the rock mass fraction for a rock/iron planet to be
zero, yielding

$\displaystyle R_{\rm Fe}=0.0975(\log M)^{2}+0.4938\log M+0.7932,$ (1)

where both mass and radius are in Earth units. An equivalent expression may be
found for the mass of a pure iron planet, $M_{\rm Fe}$, by rearranging
Equation 1 and recasting $M\rightarrow M_{\rm Fe}$ and $R_{\rm Fe}\rightarrow
R$ to give

$\displaystyle\log M_{\rm Fe}=-2.532+5.128\sqrt{0.3900R-0.0655},$ (2)

where, again, quantities are in Earth units.

To quantify our plausibility criteria, we define a weighting factor for a
measurement to be the probability that an exoplanet satisfies our plausibility
criteria, assuming that the true parameter is distributed as a Gaussian such
that $R_{\rm t}\sim N(R,\sigma_{R}^{2})$, where $R_{\rm t}$ is the true
radius, $R$ is the measured radius with uncertainty $\sigma_{R}$, and $N$ is
the normal distribution. Similarly, the true mass $M_{\rm t}$ is assumed to
follow $M_{\rm t}\sim N(M,\sigma_{M}^{2})$, where $M$ is the measured mass
with uncertainty $\sigma_{M}$. The weighting factor for radius, $W_{R}$, is
therefore

$\displaystyle W_{R}=P(R_{\rm t}\geq R_{\rm Fe}),$ (3)

where $R_{\rm Fe}$ is obtained by substituting $M$ into Equation 1. To account
for both a high-density and high-mass planet, the total mass weighting factor,
$W_{M}$, is given by

$\displaystyle W_{M}=P(M_{\rm t}\leq M_{\rm Fe})\cdot P(M_{\rm t}\leq
13M_{J}),$ (4)

where $M_{\rm Fe}$ is obtained by substituting $R$ into Equation 2. To embed
this information in a single variable, we define a combined weighting factor
to be the product of $W_{R}$ and $W_{M}$. Four planets received a combined
weighting factor of effectively zero and so were excluded from the sample,
leaving a total of 1053 exoplanets to be considered in this analysis. The
weighting factors are carried forward and considered in the mass-radius
relations developed in this paper.

##  3 Piece-wise mass-radius relations

One of the simplest deterministic models for a $M$-$R$ relation is a power
law, which can be written as

$\displaystyle\frac{R}{R_{\oplus}}=k\left(\frac{M}{M_{\oplus}}\right)^{\beta},$
(5)

where $k$ is a constant and $\beta$ is the index of the power law. In this
section we consider piece-wise power-law models. If we impose that each piece
$n$ joins to form a continuous relation, then such models reduce to a simple
$y=m_{n}x+c$ form in log space, where

$\displaystyle y=m_{1}x+c+\sum_{i=2}^{n}(m_{i-1}-m_{i})b_{i-1},$ (6)

where $m_{n}$ is the gradient of the $n^{th}$ piece, $c$ is the intercept of
the first piece and $b_{n}$ is the $n^{th}$ break-point, provided that
$b_{n-1}\leq x<b_{n}$ for $n\geq 1$ with $b_{0}=0$. We treat $m_{n}$, $b_{n}$
and $c$ as free parameters of the model, and we use orthogonal distance
regression (ODR) to fit $y=\log R$ as a function of $x=\log M$ to Equation 6,
accounting for errors in both quantities. We also incorporate the weighting
factors derived in section 2 into our analysis by combining them with the
statistical weights from the measurement errors, such that

$\displaystyle W_{{\rm tot},X}=W_{X}\cdot\frac{1}{\sigma_{X}^{2}},$ (7)

where $X$ represents mass or radius, and $W_{X}$ is the result from Equations
3 or 4. $W_{{\rm tot},X}$ in Equation 7 is the weight used in the ODR fitting
routine, hence data points may have small weights stemming from large errors,
or a small probability of satisfying our plausibility criteria, or both. These
data points therefore carry low importance to the fit.

From visual inspection of a $M$-$R$ diagram (e.g. Figure 1), it is clear that
a single power law is not optimal to describe the data. The cases in which
$n=2$ and $n=3$ both yield reasonable fits, however for the case of $n=4$, the
last break-point was found to be at a mass greater than the largest mass in
the data set. We interpret this result as an indication that no more than
three distinct regimes are supported by current data.

We use planet radius to measure the success of a model, by defining the metric

${\cal M}\equiv\left\langle\frac{|R_{o}-R_{e}|}{\sigma}\right\rangle,$ (8)

where $R_{o}$ is the observed radius, $R_{e}$ is the radius expected by the
model given the measured mass, and $\sigma$ is the radius measurement
uncertainty. The choice to use the prediction of radius, rather than mass, as
the basis for the metric is a pragmatic one, stemming from a clear
insensitivity of radius to mass for the most massive planets that is driven by
underlying physics rather than by measurement uncertainty. We find that a two-
piece model gives a value of ${\cal M}=7.81$, whereas a three-piece model
gives ${\cal M}=7.64$. For a continuous piece-wise $M$-$R$ relation current
data prefer, though not strongly prefer, three rather than two distinct
exoplanetary regimes, which we label here onwards as rocky, neptunian and
jovian. The best-fit three-piece continuous model is presented in Figure 1.

Figure 1: A log-log plot of exoplanet radius against mass, in Earth units. The
colours represent the combined weighting factor such that blue planets satisfy
our plausibility criteria completely, whereas red planets do not. The red line
shows the three-piece function fitted by ODR, and the residuals from the model
are plotted in the lower panel. The average error-normalised absolute
difference between model and measured radius is 7.64. Note the significant
residual tail that curves downward at around $50~{}M_{\oplus}$.

The transitions between regimes are found to be at $4.95\pm 0.81M_{\oplus}$
and $115\pm 19M_{\oplus}$, and the power law parameters from Equation 5 are as
follows: $k=1.01\pm 0.03$ and $\beta=0.28\pm 0.03$ in the rocky regime;
$k=0.53\pm 0.05$ and $\beta=0.68\pm 0.02$ in the neptunian regime; and
$k=13.0\pm 1.2$ and $\beta=0.012\pm 0.003$ in the jovian regime. From the
residuals, planets in the rocky regime appear to be well constrained by a
power law, however in the neptunian regime there is a noticeable systematic
downturn indicative of a population of planets that coherently deviates from
the rest of the neptunian regime. This is evidence that a continuous $M$-$R$
may not be appropriate, and that the rocky and neptunian planetary regimes
overlap (e.g. Otegi et al., 2020). We consider a discontinuous model in
Section 4.2.

##  4 Beyond a continuous mass-radius relation

###  4.1 Temperature Considerations

Many hot jupiters show evidence of radius bloating due to high levels of
stellar insolation. This motivates a need to include temperature by
considering mass-radius-temperature $M$-$R$-$T$ relations for the jovian
regime. In order to include a temperature dependence, we extend the use of ODR
to fit a plane in mass-radius-temperature log space, accounting for errors in
all measured quantities. We carry forward the same analysis of weighting
factors that was performed in section 3, resulting in a multiplicative power
law of the form

$\displaystyle\frac{R}{R_{\oplus}}=CT_{\rm
eq}^{\beta_{1}}\left(\frac{M}{M_{\oplus}}\right)^{\beta_{2}},$ (9)

where $C$, $\beta_{1}$ and $\beta_{2}$ are constants and $T_{\rm eq}$ is the
equilibrium temperature of the planet. For the rocky and neptunian regimes,
and for jovian planets without an equilibrium temperature measurement, a
simple $M$-$R$ power law is used to predict their radii, relaxing the
constraint that the global $M$-$R$ relation must be continuous.

Figure 2: A $M$-$R$ diagram in which planets have been classified into rocky,
neptunian and jovian regimes according to the mass break-points
$4.95M_{\oplus}$ and $115M_{\oplus}$. The black line indicates the three-piece
$M$-$R$ relation found in section 3, and the contours of the $M$-$R$-$T$
relation in the jovian regime are plotted in purple. The size of the markers
is proportional to the weighting factor. Rocky planets are coloured brown,
neptunian planets are coloured cyan, and jovian planets follow a colour map
according to their equilibrium temperature. jovian planets without an
equilibrium temperature measurement are plotted in grey. The residuals are
plotted in the lower panel, as the measured radius subtracted by the radius in
the model. The average absolute difference between the model and measurements
is 5.88.

Combining the continuous three-piece model from section 3 with the $M$-$R$-$T$
model for jovian planets from Equation 9 provides the semi-continuous
$M$-$R$-$T$ model shown in Figure 2. The systematic downturn in the residuals
of the neptunian regime due to misclassified rocky planets is once again
evident as before in Figure 1. However, the scatter in jovian radii is
reduced, with ${\cal M}=5.88$, illustrating the importance of modeling
equilibrium temperature.

###  4.2 Discontinuous mass-radius-temperature relations

The coherent nature of the residual excess seen in the neptunian regime in
Figures 1 and 2 provides clear support for discontinuity in the transition
from the rocky to neptunian regime. Whilst different planetary regimes can be
segregated by specific mass or radius break points, we need a way to define
distinct regions of the $M$-$R$ plane in order to consider discontinuous
models (c.f. Otegi et al., 2020). To separate rocky from neptunian planets, we
consider the ice/rock equation of state of Fortney et al. (2007). The radius
is given by

$\displaystyle R$ $\displaystyle=$ $\displaystyle(0.0592f_{\rm
ice}+0.0975)(\log M)^{2}$ (10) $\displaystyle+(0.2337f_{\rm ice}+0.4938)\log
M$ $\displaystyle+(0.3102f_{\rm ice}+0.7932),$

where $f_{\rm ice}$ is the ice mass fraction (1 for pure ice and 0 for pure
rock), and mass and radius are in Earth units. We choose first to classify
rocky planets as those which, for some fixed value of $f_{\rm ice}$, have a
radius less than that calculated by substituting their mass measurement into
Equation 10. For the remaining planets, neptunian planets are defined as
having a mass less than a mass break-point $M_{\rm br}=115~{}M_{\oplus}$,
while jovian planets have a mass larger than $M_{\rm br}$. After
classification, discontinuous power law $M$-$R$ relations are fitted to the
rocky and neptunian regimes, and a $M$-$R$-$T$ relation is fitted to the
jovian regime as outlined in section 4.1. We investigate the impact of our
choice of $f_{\rm ice}$ using metric ${\cal M}$ defined by Equation 8. The
dependence of ${\cal M}$ on $f_{\rm ice}$ is shown in Figure 3.

Figure 3: A plot of the error-scaled average absolute difference between
radius measurements and predictions [metric ${\cal M}$ in Equation (8)],
against the assumed ice mass fraction $f_{\rm ice}$ used in Equation 10 to
separate rocky from neptunian planets. A mass break-point of $M_{\rm
br}=115M_{\oplus}$ is adopted to separate neptunian and jovian planets.

In Figure 3 we see a general downwards trend in ${\cal M}$, suggesting that a
larger ice mass fraction is favourable for separating the super-Earth and
mini-neptune regimes. Indeed, current data supports convergence to the
approach of Otegi et al. (2020), in which the composition line of water was
used to separate planets. The observational discontinuity of the two regimes
is consistent with a physical interpretation of a relatively sharp transition
from a rock/ice to an ice giant regime. For our discontinuous $M$-$R$-$T$
relation we therefore adopt $f_{\rm ice}=1$ to distinguish rocky planets from
neptunes, and $M_{\rm br}=115~{}M_{\oplus}$ to segregate neptunian and jovian
planets. The resulting $M$-$R$-$T$ relation is presented in Figure 4.

Figure 4: A $M$-$R$ diagram in which planets have been classified into rocky,
neptunian and jovian regimes according to the equation of state of a pure ice
planet and a mass break-point of $115M_{\oplus}$. The black lines indicate the
power law $M$-$R$ relations fitted to the data, and fixed example contours of
the $M$-$R$-$T$ relation in the jovian regime are plotted in purple. The size
of the markers is proportional to the combined weighting factor. Rocky planets
are coloured brown, neptunian planets are coloured cyan, and jovian planets
follow a colour map according to their equilibrium temperature. jovian planets
without an equilibrium temperature measurement are plotted in grey. The error-
normalised residuals are plotted in the lower panel, with an average absolute
difference of 3.81.

From this model we find that when compared to Equation 5, the rocky regime
yields $k=0.99\pm 0.02$ and $\beta=0.34\pm 0.01$; the neptunian regime gives
$k=0.97\pm 0.07$ and $\beta=0.55\pm 0.02$; and the jovian regime in the
absence of equilibrium temperature data gives $k=8.01\pm 0.48$ and
$\beta=0.087\pm 0.001$. Comparing to Equation 9, the radii in the jovian
regime are best described by $C=1.10\pm 0.15$ with temperature index
$\beta_{1}=0.35\pm 0.02$ and mass index $\beta_{2}=0.00\pm 0.01$. It is
interesting that even the weak dependence of jovian radius on mass seen in
Figure 1 can apparently be explained away as pure temperature dependence.

In the residuals in Figure 4, the radii of rocky planets appear to be well
modelled, supported by the small uncertainty in the index of the power law in
this regime. There remains some scatter in the neptunian regime, although the
systematic down-turn from Figure 1 is no longer present. Some scatter also
remains in the jovian regime, and there is a slight downward trend in the
residuals, which is made more apparent in Figure 5.

Figure 5: A plot of predicted radius against measured radius for jovian
planets. The colour map represents the equilibrium temperature and the line
for which the predicted radius is equal to the measured radius is plotted in
black. On the left, the parameters for the model in Equation 9 are as plotted
in Figure 4. The planets shown have been selected using the mass break point
$M_{\rm br}=115~{}M_{\oplus}$. The distribution of hotter planets is clearly
skewed with respect to the solid line demarcating perfect agreement between
measured and predicted radii. On the right, we illustrate how one could
correct for this by refitting only to the subset of planets selected using the
friends-of-friends algorithm of Huchra & Geller (1982), with a clustering
parameter, $b=0.3$. However, in this case we recover worsened predictions for
the radii of cooler jovian worlds.

From the left panel in Figure 5, it appears that the model is skewed away from
the visible trend by some cooler planets that have radii larger than expected.
This is likely an effect of implementing a break-point, as the transition
between neptunian and jovian planets may itself be temperature sensitive. The
right panel of Figure 5 uses the friends-of-friends algorithm of Huchra &
Geller (1982) with a clustering parameter, $b=0.3$, to isolate the main group
of hotter jovian planets in $M$-$R$-$T$ space. We then use this data to refit
parameters in the $M$-$R$-$T$ model. This new model is used to predict the
radii of the same planets in the left panel for direct comparison. As this
approach essentially removes outlying planets, the predicted radii correlate
better with the measured radii for the bulk of the jovian planets. However, it
does worsen the predictions for cooler planets, though some of these may well
be misclassified neptunian planets. A fixed mass break-point between neptunes
and jovian planets is well-motivated by the general trend, but may not be
optimal.

##  5 Comparison of relations

Forecaster is a widely-used, publicly-available program for predicting
exoplanet masses and radii on a probabilistic basis, developed by Chen &
Kipping (2017). For a given mass or radius with optional uncertainties,
Forecaster returns a prediction of the other measurement based on how likely
the planet is to fall into each regime and the uncertainty regions surrounding
the model. As this model is probabilistic, a different prediction value will
be returned each time the program is run. The terran, neptunian, jovian and
stellar regimes are split by mass break-points, and a continuous broken power
law relates mass and radius, where each segment of the power law has a
different uncertainty region. Figure 6 shows the $M$-$R$ diagram coloured
according the the probability Forecaster has assigned to each planet of being
neptunian. The residuals for one run of Forecaster are also shown and coloured
according to their bulk density.

Figure 6: A plot of radius against mass in the upper panel, where the colour
corresponds to the probability Forecaster has assigned to each planet of being
neptunian. The residuals plotted in the lower panel were calculated as the
measured radii subtracted by the predicted radii generated from a single run
of Forecaster, with the colour representing the bulk density of each planet.

Figure 6 demonstrates that Forecaster suffers from the same systematic
downturn in the residuals of the neptunian regime as our three-piece model in
Figure 1. From the residuals panel, it is also apparent that these planets
with over-predicted radii are among the densest planets, supporting the
presumption that these are misclassified rocky planets. This is once again a
feature of using a mass break-point to distinguish between rocky and neptunian
planets. Additionally, a large scatter in radii is evident in the residuals
for large masses due to an absence of temperature dependence in the model.

The slope of a log-log $M$-$R$ plot for our three-piece model is consistent
with those found both by Chen & Kipping (2017) and by Otegi et al. (2020) in
the rocky regime; all around 0.28. However, our discontinuous model gives a
steeper slope of $0.34\pm 0.01$ in this regime, consistent with $R\propto
M^{\frac{1}{3}}$ as expected for a solid body with constant density across all
rocky planets. This contrasts with Schulze et al. (2021) in that it argues for
a rocky planet bulk density that is, on average, insensitive to the
composition of its host. In the neptunian regime, the three-piece model gives
a slope larger than that used by Forecaster, potentially due to the smaller
break-point found by Chen & Kipping (2017) to separate between rocky and
neptunian planets. Our discontinuous model, on the other hand, gives a slope
smaller than that found by both Otegi et al. (2020), and Chen & Kipping
(2017). The difference from Forecaster can be attributed to the use of a
density-based classification scheme rather than a mass break-point, however
this was also the approach taken by Otegi et al. (2020). This difference in
slope may arise from the subtly different cutoff for higher mass planets of
$120M_{\oplus}$ compared to $115M_{\oplus}$ used in our model, as well as from
differences in the data used and in fitting approaches.

The slopes of the jovian $M$-$R$ relation found by Chen & Kipping (2017) and
our three-piece and discontinuous models are not consistent with one another.
We expect that this is due to the large scatter in radii and weak mass-
dependence, hence subtle changes to the exact data included or the location of
the break-point can have a significant impact on the results of the fit.
Furthermore, Chen & Kipping (2017) included brown dwarfs in their sample,
which anchor the fit behaviour within the jovian regime. Nonetheless, all of
the $M$-$R$ slopes are very shallow and therefore near degenerate. In our
discontinuous $M$-$R$-$T$ model, we find that $M$-$T$ dependence has a slope
of 0.35 $\pm$ 0.02, but that $M$ is essentially uncorrelated to $R$, with a
slope of $0.00\pm 0.01$.

The use of a density-based classification scheme has the interesting
consequence that planets we classify to be rocky have masses up to
$79~{}M_{\oplus}$, comparable to relatively high mass neptunian bodies. Whilst
Chen & Kipping (2017) found that data at the time indicated that rocky super-
Earths were not nearly as large as expected, current data indicates rocky
planets can be as large as $4.3~{}R_{\oplus}$.

Figure 7: Plots of the radii predicted by a two-piece function, three piece
function, discontinuous $M$-$R$-$T$ model and continuous $M$-$R$-$T$ model
against measured radius. The size of the markers is proportional to the
combined weighting factors used to down-weight points in the ODR fit which did
not satisfy the plausibility criteria. The line for which the predicted radius
is equal to the measured radius is plotted in black.

A plot of predicted radii against measured radii for the four models we have
developed are displayed in Figure 7. The 2-piece relation in Figure 7 shows
significant deviation from the unity line for the smallest planets, while for
the other three models agreement is much stronger, confirming that the data
provides strong support for three distinct planetary regimes. Both $M$-$R$
models in Figure 7 show a systematic degeneracy between measured and predicted
radii for the largest planets. There is much better correlation of predictions
with measurements for the $M$-$R$-$T$ models. All continuous models show large
scatter around the unity line for intermediate mass planets due to the
superposition of the super-Earth and mini-neptune regimes. This scatter is
greatly reduced by the discontinuous $M$-$R$-$T$ model.

Figure 8: A plot of the radius predicted by one run of Forecaster against
measured radius. The line for which the predicted radius is equal to the
measured radius is plotted in black.

The large scatter for intermediate mass planets is similarly seen in Figure 8,
which shows the results of predicted versus measured radii for a single run of
Forecaster. Additionally, whilst a probabilistic model goes some way towards
reproducing the radius scatter of large planets, it does not perform as well
as a model that includes temperature to account for bloating.

We can use Equation 8 to quantitatively assess the success of each model via
${\cal M}$. Since Forecaster is a probabilistic model, ${\cal M}$ will differ
for each run of the program. We therefore opt to calculate it for 10,000 runs
and to fit a Gaussian to the resulting distributions. The results are
displayed in Figure 9.

Figure 9: Metric ${\cal M}$ calculated for our two- and three-piece $M$-$R$
models as well as our continuous and discontinuous $M$-$R$-$T$ models applied
to the current data set. The distribution of ${\cal M}$ for 10,000 runs of
Forecaster has been fitted to a Gaussian distribution, with $\langle{\cal
M}\rangle=10.1$ and $\sigma({\cal M})=0.2$.

Figure 9 clearly demonstrates that when considering the current data
available, there is no case in which Forecaster’s average performance is
better than that of our deterministic models. $M$-$R$-$T$ models outperform
$M$-$R$ models and a discontinuous model is favoured over a continuous one.

One reason why Forecaster may not perform well compared to our models may
simply be because it is conditioned on much older data. At the time of
writing, the public version of
Forecaster222https://github.com/chenjj2/forecaster is calibrated upon data
available at the time of its initial development, which is prior to 2017. To
fairly compare models, we consider how our models predict data available to
Chen & Kipping (2017), though with some caveats.

As we are interested in an exoplanet mass-radius relation, we neglect any
measurements for objects that are classified as stars or brown dwarfs.
Furthermore, like Chen & Kipping (2017) we include solar system bodies in our
sample, but we use current measurements and uncertainties taken from the JPL
Planetary Satellite Physical
Parameters333https://ssd.jpl.nasa.gov/sats/phys_par/ database, Williams et al.
(2014) and the JPL Planetary Physical
Parameters444https://ssd.jpl.nasa.gov/planets/phys_par.html database. The full
table of solar system parameters is presented in Table 1 of the Appendix,
along with our equilibrium temperature calculation for Jupiter. We use these
updated values both within Forecaster and for our own relations for
consistency. The uncertainties for these bodies are small, hence we do not
expect their revision to significantly alter the performance of Forecaster.

We compare ${\cal M}$ for the dataset subset that existed at the time of Chen
& Kipping (2017) in Figure 10, with solar system objects excluded. This leaves
only one exoplanet in the old data set that Forecaster classifies as rocky,
and this lies within the uncertainty region of the $2M_{\oplus}$ break-point.
As a result, Forecaster has to rely on the inclusion of solar system values in
this regime. In Figure 11 we show the result of ${\cal M}$ for the model
predictions of the size of solar system bodies.

Figure 10: Metric ${\cal M}$ for our two- and three-piece $M$-$R$ models as
well as our continuous and discontinuous $M$-$R$-$T$ models applied to the old
data set without the inclusion of solar system bodies. The distribution of
${\cal M}$ for 10,000 runs of Forecaster has been fitted to a Gaussian
distribution with $\langle{\cal M}\rangle=4.68$ and $\sigma({\cal M})=0.06$.

In Figure 10 we see that while our $M$-$R$-$T$ models still significantly
outperform Forecaster, our $M$-$R$ models are now comparable. Using the mean
and width of the Gaussian fitted to the Forecaster ${\cal M}$ distribution, we
find that Forecaster can outperform the two-piece $M$-$R$ model 99.9% of the
time, but can out-perform our three-piece $M$-$R$ model only 8% of the time.

Figure 11: Metric ${\cal M}$ for our two- and three-piece $M$-$R$ models as
well as our continuous and discontinuous $M$-$R$-$T$ models applied to solar
system bodies. The ${\cal M}$ axis is four orders of magnitude larger than in
previous similar plots and so the metrics calculated for the $M$-$R$-$T$
models cannot be resolved in this figure. The distribution of ${\cal M}$ for
10,000 runs of Forecaster has been fitted to a Gaussian with ${\cal M}=9000$
and $\sigma({\cal M})=44000$. None of the models is able to predict the radii
of solar systems bodies to within their very small current uncertainties,
though the $M$-$R$-$T$ models perform much better, in relative terms, than the
other models.

We see in Figure 11 that our two-piece model performs comparatively poorly for
predicting the size of solar system bodies. Our other models are able to
extrapolate more reliably down to this low-mass regime though, unsurprisingly,
nowhere near the precision of current measurement uncertainty for these
bodies. We find that Forecaster outperforms a three-piece $M$-$R$ relation 67%
of the time, a continuous $M$-$R$-$T$ model 6% of the time and a discontinuous
$M$-$R$-$T$ model 3% of the time. We conclude that a $M$-$R$-$T$ model
calibrated on current data nearly always outperforms Forecaster calibrated on
data prior to 2017. We expect that the performance of Forecaster could be
significantly improved by updating the hyper-parameters of the Chen & Kipping
(2017) model to include current measurements, though the means to do so have
not been made publicly available. Nonetheless, the weaknesses of Forecaster
are also those inherent to a continuous $M$-$R$ model with mass break-points,
and as such are unlikely to be fully mitigated by updating the input dataset.

One aspect that has not been accounted for in any of the models discussed in
this report is the effect of detection and selection bias on the mass-radius
dataset. With the exception of Burt et al. (2018), prioritisation schemes for
the follow-up of exoplanet detections are not generally made available and
there are very few investigations into how these schemes may introduce bias
into the population of planets for which we have both a mass and a radius
measurement. This therefore makes it very difficult to de-bias $M$-$R$
relations calibrated on observations.

As for detection bias, the densities of planets calculated using transit
timing variations (TTV) tend be smaller than those of planets with mass
measurements obtained from radial velocities. Leleu et al. (2022) has found
evidence to suggest that some of this discrepancy is due to differing
detection biases, and they correct these measurements accordingly. A similar
attempt to account for biased TTV planet measurements is presented in Jontof-
Hutter et al. (2020). The number of TTV planets in our sample is however
vastly out-numbered by the number of planets meaaured using radial velocities,
so we expect any TTV biasing effect to be small.

##  6 Conclusions

We have compiled a catalogue of 1053 confirmed exoplanets with which to
calibrate power-law exoplanet mass-radius ($M$-$R$) and mass-radius-
temperature ($M$-$R$-$T$) relationships. We have strived to let the data
itself inform us as to the piece-wise structure of these relationships,
including whether continuous or discontinuous power-law forms are preferred.

Using orthogonal distance regression fits that account for errors in both mass
and radius, we find that current data is best explained by three distinct
planetary regimes that, under a continuous $M$-$R$ relation, transition from a
rocky to an ice giant (neptunian) regime at $4.95\pm 0.81~{}M_{\oplus}$ and
from ice giant to gas giant (jovian) at $115~{}M_{\oplus}$.

We find that the modeling of the jovian regime is improved through inclusion
of the effect of bloating via extension to an $M$-$R$-$T$ relation. In fact,
when doing this, we find that $M\propto R^{0.00\pm 0.01}T^{0.35\pm 0.02}$, so
that jovian mass planets can be well modeled with no radius dependence at all.

Our analysis also finds strong support from the data for a discontinuous
$M$-$R$-$T$ relation between rocky and neptunian planets, as has been
previously argued by Otegi et al. (2020). Modeling the boundary with analytic
ice-rock equations of state from Fortney et al. (2007) we find that the data
prefers a boundary corresponding to a pure-ice world, giving support for the
physical interpretation of the discontinuity as separating rocky from ice-
giant (neptunian) planet populations. Interestingly, we find that the
resulting upper mass of planets categorized within the rocky planet regime can
extend almost up to the upper mass limit of the neptunian population.

Given the significant increase in the amount of exoplanet data since the
publication of the widely-used Forecaster code (Chen & Kipping, 2017) we find
most of our models outperform Forecaster in the accuracy of radius
predictions. While this can to an extent be attributed to the hyper-parameters
of the Forecaster model being conditioned on an older and smaller dataset, the
models which perform the best against it are those that allow for
discontinuities arising from variations in temperature and equation of state
that are not included in the underlying $M$-$R$ model used in Forecaster.

Looking ahead, the current exoplanet dataset will see a massive expansion over
the coming decade, thanks largely to astrometric detection by ESA Gaia
(Perryman et al., 2014) and by the transit and microlensing samples of the
NASA Nancy Grace Roman Space Telescope (Penny et al., 2019; Wilson et al.,
2023). Roman alone is expected to expand the exoplanet catalogue from the
current size of under 6,000 planets to at least 60,000 and possibly up to
200,000 hot transiting planets, as well around 1,400 cool microlensing
planets. As these datasets will involve combinations of mass (astrometry and
microlensing) and radius (transit) measurements, a coherent analysis of
exoplanet demography will require an increasingly precise modelling of the
exoplanet $M$-$R$-$T$ relation.

## Acknowledgements

This research has made use of the NASA Exoplanet Archive, which is operated by
the California Institute of Technology, under contract with the National
Aeronautics and Space Administration under the Exoplanet Exploration Program.

## References

  * Bashi et al. (2017) Bashi D., Helled R., Zucker S., Mordasini C., 2017, A&A, 604, A83
  * Batygin et al. (2011) Batygin K., Stevenson D. J., Bodenheimer P. H., 2011, ApJ, 738, 1
  * Burt et al. (2018) Burt J., Holden B., Wolfgang A., Bouma L. G., 2018, AJ, 156, 255
  * Chen & Kipping (2017) Chen J., Kipping D., 2017, ApJ, 834, 17
  * Enoch et al. (2012) Enoch B., Collier Cameron A., Horne K., 2012, A&A, 540, A99
  * Fortney et al. (2007) Fortney J. J., Marley M. S., Barnes J. W., 2007, ApJ, 659, 1661
  * Huchra & Geller (1982) Huchra J. P., Geller M. J., 1982, ApJ, 257, 423
  * Jontof-Hutter et al. (2020) Jontof-Hutter D., Ford E., Lissauer J., Wolfgang A., Rowe J., Fabrycky D., 2020, in AAS/Division for Planetary Sciences Meeting Abstracts. p. 303.02
  * Leleu et al. (2022) Leleu A., et al., 2022, in European Planetary Science Congress. pp EPSC2022–61, doi:10.5194/epsc2022-61
  * Mocquet et al. (2014) Mocquet A., Grasset O., Sotin C., 2014, Philosophical Transactions of the Royal Society of London Series A, 372, 20130164
  * Mousavi-Sadr et al. (2023) Mousavi-Sadr M., Jassur D. M., Gozaliasl G., 2023, MNRAS,
  * Ning et al. (2018) Ning B., Wolfgang A., Ghosh S., 2018, ApJ, 869, 5
  * Otegi et al. (2020) Otegi J. F., Bouchy F., Helled R., 2020, A&A, 634, A43
  * Penny et al. (2019) Penny M. T., Gaudi B. S., Kerins E., Rattenbury N. J., Mao S., Robin A. C., Calchi Novati S., 2019, ApJS, 241, 3
  * Perryman et al. (2014) Perryman M., Hartman J., Bakos G. Á., Lindegren L., 2014, ApJ, 797, 14
  * Schulze et al. (2021) Schulze J. G., Wang J., Johnson J. A., Gaudi B. S., Unterborn C. T., Panero W. R., 2021, The Planetary Society Journal, 2, 113
  * Seager et al. (2007) Seager S., Kuchner M., Hier-Majumder C. A., Militzer B., 2007, ApJ, 669, 1279
  * Sestovic et al. (2018) Sestovic M., Demory B.-O., Queloz D., 2018, A&A, 616, A76
  * Spiegel et al. (2011) Spiegel D. S., Burrows A., Milsom J. A., 2011, ApJ, 727, 57
  * Swift et al. (2012) Swift D. C., et al., 2012, ApJ, 744, 59
  * Thorngren & Fortney (2018) Thorngren D. P., Fortney J. J., 2018, AJ, 155, 214
  * Ulmer-Moll et al. (2019) Ulmer-Moll S., Santos N. C., Figueira P., Brinchmann J., Faria J. P., 2019, A&A, 630, A135
  * Weiss et al. (2013) Weiss L. M., et al., 2013, ApJ, 768, 14
  * Williams et al. (2014) Williams J. G., et al., 2014, Journal of Geophysical Research (Planets), 119, 1546
  * Wilson et al. (2023) Wilson R. F., et al., 2023, arXiv e-prints, p. arXiv:2305.16204
  * Wolfgang et al. (2016) Wolfgang A., Rogers L. A., Ford E. B., 2016, ApJ, 825, 19

## Appendix

The only solar system object with a mass greater than $115M_{\oplus}$ and thus
requiring an equilibrium temperature measurement is Jupiter. Assuming that
both the Sun and Jupiter behave as blackbodies, it can be shown that the
equilibrium temperature of Jupiter, $T_{\rm eq,J}$, is

$\displaystyle T_{\rm
eq,J}=T_{\odot}\left(\frac{R_{\odot}}{2a}\right)^{\frac{1}{2}}(1-A)^{\frac{1}{4}},$
(11)

where $T_{\odot}=5772$ K is the solar effective temperature, $R_{\odot}$ is
the solar radius, $a=5.2$ au is the separation between Jupiter and the Sun,
and $A=0.5$ is the Bond albedo of Jupiter. The substitution of these values
into Eqn (11) and propagating their uncertainties yields an equilibrium
temperature for Jupiter of $102.3\pm 0.6$ K.

The masses and radii of solar system objects are presented in Table 1.

Table 1: The solar system objects considered in Forecaster with updated values for their masses and mean radii including uncertainties. Parameters for moons are taken from the JPL Planetary Satellite Physical Parameters database, excepting the mass of the Moon, which was taken from Williams et al. (2014). All planet parameters are taken from the JPL Planetary Physical Parameters database. Object | Mass (kg) | Mean Radius (km)
---|---|---
Moons |  |
Moon | (7.3463 $\pm$ 0.0088) $\cdot 10^{22}$ | 1737.4 $\pm$ 0.1
Io | (8.931938 $\pm$ 0.000018) $\cdot 10^{22}$ | 1821.49 $\pm$ 0.5
Europa | (4.799844 $\pm$ 0.000013) $\cdot 10^{22}$ | 1560.8 $\pm$ 0.3
Ganymede | (1.4819 $\pm$ 0.0001) $\cdot 10^{23}$ | 2631.2 $\pm$ 1.7
Callisto | (1.07594 $\pm$ 0.00014) $\cdot 10^{23}$ | 2410.3 $\pm$ 1.5
Rhea | (2.306520 $\pm$ 0.000035) $\cdot 10^{21}$ | 2410.3 $\pm$ 1.5
Titan | (1.3452 $\pm$ 0.0002) $\cdot 10^{23}$ | 2574.76 $\pm$ 0.02
Titania | (3.400 $\pm$ 0.061) $\cdot 10^{21}$ | 788.9 $\pm$ 1.8
Oberon | (3.076 $\pm$ 0.087) $\cdot 10^{21}$ | 761.4 $\pm$ 2.6
Triton | (2.1390 $\pm$ 0.0028) $\cdot 10^{22}$ | 1352.6 $\pm$ 2.4
Dwarf Planets |  |
Eris | (1.660 $\pm$ 0.020) $\cdot 10^{22}$ | 1200 $\pm$ 50
Pluto | (1.3029 $\pm$ 0.0027) $\cdot 10^{22}$ | 1188.3 $\pm$ 1.6
Planets |  |
Mercury | (3.30103 $\pm$ 0.00021) $\cdot 10^{23}$ | 2439.4 $\pm$ 0.1
Venus | (4.86731 $\pm$ 0.00023) $\cdot 10^{24}$ | 6051.8 $\pm$ 1.0
Earth | (5.97217 $\pm$ 0.00028) $\cdot 10^{24}$ | 6371.0080 $\pm$ 0.0001
Mars | (6.41691 $\pm$ 0.00030) $\cdot 10^{23}$ | 3389.5 $\pm$ 0.2
Jupiter | (1.898125 $\pm$ 0.000088) $\cdot 10^{27}$ | 69911 $\pm$ 6
Saturn | (5.68317 $\pm$ 0.00026) $\cdot 10^{26}$ | 58232 $\pm$ 6
Uranus | (8.68099 $\pm$ 0.0004) $\cdot 10^{25}$ | 25362 $\pm$ 7
neptune | (1.024092 $\pm$ 0.000048) $\cdot 10^{26}$ | 24622 $\pm$ 19