# Search for galactic axions with a high-Q dielectric cavity

D. Alesini  INFN, Laboratori Nazionali di Frascati, Frascati, Roma, Italy  D.
Babusci  INFN, Laboratori Nazionali di Frascati, Frascati, Roma, Italy  C.
Braggio  INFN, Sezione di Padova, Padova, Italy  Dipartimento di Fisica e
Astronomia, Padova, Italy  G. Carugno  INFN, Sezione di Padova, Padova, Italy
Dipartimento di Fisica e Astronomia, Padova, Italy  N. Crescini  Present
address: University Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 38000
Grenoble, France  INFN, Laboratori Nazionali di Legnaro, Legnaro, Padova,
Italy  Dipartimento di Fisica e Astronomia, Padova, Italy  D. D’Agostino
Dipartimento di Fisica E.R. Caianiello, Fisciano, Salerno, Italy  INFN,
Sezione di Napoli, Napoli, Italy  A. D’Elia<EMAIL_ADDRESS>INFN, Laboratori Nazionali di Frascati, Frascati, Roma, Italy  D. Di
Gioacchino  INFN, Laboratori Nazionali di Frascati, Frascati, Roma, Italy  R.
Di Vora<EMAIL_ADDRESS>INFN, Sezione di Padova, Padova, Italy
Dipartimento di Scienze Fisiche, della Terra e dell’Ambiente, Università di
Siena, 53100 Siena, Italy  P. Falferi  Istituto di Fotonica e Nanotecnologie,
CNR Fondazione Bruno Kessler, I-38123 Povo, Trento, Italy  INFN, TIFPA, Povo,
Trento, Italy  U. Gambardella  Dipartimento di Fisica E.R. Caianiello,
Fisciano, Salerno, Italy  INFN, Sezione di Napoli, Napoli, Italy  C. Gatti
INFN, Laboratori Nazionali di Frascati, Frascati, Roma, Italy  G. Iannone
Dipartimento di Fisica E.R. Caianiello, Fisciano, Salerno, Italy  INFN,
Sezione di Napoli, Napoli, Italy  C. Ligi  INFN, Laboratori Nazionali di
Frascati, Frascati, Roma, Italy  A. Lombardi  INFN, Laboratori Nazionali di
Legnaro, Legnaro, Padova, Italy  G. Maccarrone  INFN, Laboratori Nazionali di
Frascati, Frascati, Roma, Italy  A. Ortolan  INFN, Laboratori Nazionali di
Legnaro, Legnaro, Padova, Italy  R. Pengo  INFN, Laboratori Nazionali di
Legnaro, Legnaro, Padova, Italy  A. Rettaroli  INFN, Laboratori Nazionali di
Frascati, Frascati, Roma, Italy  G. Ruoso  INFN, Laboratori Nazionali di
Legnaro, Legnaro, Padova, Italy  L. Taffarello  INFN, Sezione di Padova,
Padova, Italy  S. Tocci  INFN, Laboratori Nazionali di Frascati, Frascati,
Roma, Italy

###### Abstract

A haloscope of the QUAX–$a\gamma$ experiment, composed of an high-Q resonant
cavity immersed in a 8 T magnet and cooled to $\sim 4.5$ K is operated to
search for galactic axion with mass $m_{a}\simeq 42.8~{}\mu\text{eV}$. The
design of the cavity with hollow dielectric cylinders concentrically inserted
in a OFHC Cu cavity, allowed us to maintain a loaded quality-factor Q $\sim
300000$ during the measurements in presence of magnetic field. Through the
cavity tuning mechanism it was possible to modulate the resonance frequency of
the haloscope in the region $10.35337-10.35345$ GHz and thus acquire different
dataset at different resonance frequencies. Acquiring each dataset for about
50 minutes, combining them and correcting for the axion’s signal estimation-
efficiency we set a limit on the axion-photon coupling
$g_{a\gamma\gamma}<0.731\times 10^{-13}$ GeV-1 with the confidence level set
at $90\%$.

##  I Introduction

The axion is an hypothetical particle that was theorized to solve the strong
CP problem. It arises from the spontaneous breaking of the Peccei-Quinn
symmetry of QCD [1, 2, 3]. In addition, the properties predicted for the
axion, charge neutrality, spin 0 and negligible interaction with the ordinary
matter, make this particle a strong candidate for the dark matter [4].
Cosmology and astrophysical considerations, suggest an axion mass range
$1~{}\mu\text{eV}<m_{a}<10~{}\text{meV}$ [5]. The hunt for axion is now world
spread and most of the experiments involved in this search use detectors based
on the haloscope design proposed by Sikivie [6, 7]. Among them are ADMX [8, 9,
10, 11], HAYSTAC [12, 13], ORGAN [14], CAPP-8T [15, 16], CAPP-9T [17], CAPP-
PACE [18], CAPP-18T [19], GrAHal [20], RADES [21, 22, 23], TASEH [24], QUAX
[25, 26, 27, 28, 29], and KLASH [30, 31]. Dielectric and plasma haloscopes
have also been proposed, like MADMAX [32] and ALPHA [33], respectively. The
haloscope concept is based on the immersion of a resonant cavity in a strong
magnetic field in order to stimulate the inverse Primakoff effect, converting
an axion into an observable photon [34]. To maximise the power of the
converted axion, it’s necessary to maximise the cavity quality-factor ($Q$)
and to tune the resonance frequency to match the axion mass. Different
solutions have being adopted to maximize the signal-to-noise ratio, facing the
problem from different angles. Resonant cavities of superconductive and
dielectric materials are becoming increasingly popular because of their high Q
[35, 36, 37, 38]. In this work we describe the results obtained operating the
haloscope of the QUAX–$a\gamma$ experiment using an high-Q dielectric cavity
immersed in a static magnetic field of 8 T and cooled down to $\sim 4.5$ K.
The results obtained allow us to exclude values of
$g_{a\gamma\gamma}>0.729\times 10^{-13}$ GeV-1 at 90% confidence level (C.L.)
in a region of mass 1.32 neV wide, centered at $42.8216\,\mu$eV.

##  II Experimental Setup

###  II.0.1 General description

The core of the haloscope is an extremely high-Q resonant cavity. The cavity
is extensively described in [39]: it is based on a right circular copper
cavity with hollow sapphire cylinders that confine higher order modes around
the cylinder axis. The useful mode is a TM030, which has an effective volume
$V\cdot C_{030}=3.4\times 10^{-2}$ liters at the resonant frequency of 10.3
GHz, where $C_{030}$ is a geometrical factor entering in the signal-power
estimation in Eq. (4). Under an 8 T-field we measured an internal quality
factor of more than $9\times 10^{6}$. The cavity and the magnet are hosted
inside a liquid-He cryostat at a temperature of about 4 K.

The principle scheme of the measurement set up is shown in Figure 1. The
microwave cavity is immersed in a 8 T maximum magnetic field, not shown in the
figure, generated by a 150 mm diameter bore of 500 mm length superconducting
magnet. When the magnet is driven by a 92 A current, the effective squared
field over the cavity length amount to 50.8 T2. The microwave cavity is read
by a tunable monopole antenna with coupling $\beta$. This is obtained by
acting on a manually-controlled mechanical feed-through, that allows for
$\beta$ values in the range 0.01 to 20. A weakly coupled port (coupling about
0.01) is used for calibration purposes and is connected to the room
temperature electronics by means of line L1. To avoid thermal power inputs
from room temperature, a 20 dB attenuation is inserted on L1. Cavity tuning
was obtained by displacing triplets of 2 mm-diameter sapphire rods relative to
the top and bottom cavity endcaps [39]. Again, independent motion of the two
triplets is obtained by manually controlled mechanical feed-throughs.

Figure 1: Schematics of the experimental apparatus. The microwave cavity
(orange) is immersed in the uniform magnetic field (not shown) generated by
the magnet. C1, C2 and C3 are circulators, HP is a 8 GHz high pass filter,
cryo HEMT and HEMT are High Electron Mobility Transistors. Attenuators are
shown with their reduction factor in decibels. The horizontal blue line
identifies the boundaries of the cryogenic stages of the apparatus. KF 40 are
rf feed-through with ISO-KF 40 flanges.

The power collected by the tunable antenna is amplified by a cryogenic high
electron mobility transistor amplifier (cryo HEMT in figure), isolated from
the cavity by means of the circulators C1 and C2. The output of the cryo HEMT
is filtered and then transmitted along line L3 to the room temperature
electronics, where it is first amplified by a room temperature HEMT and then
processed for data storage. The room temperature chain is the same used in our
previous measurements [26]: the HEMT output is frequency down-converted using
a mixer with the local oscillator frequency set to a value about 500 kHz below
the cavity resonance. The low frequency in phase and quadrature outputs of the
mixer are amplified and then sampled with a 2 Ms/s analog to digital converter
(ADC) and stored on a computer for off line data analysis. Data storage is
done with blocks of about 4 s of sampled data for both output channels of the
mixer. An auxiliary line L2 is used for calibration purposes: it is connected
to the line L3 by means of the circulator C1, and 32 dB of attenuation
prevents thermal leakage from room temperature components.

The room temperature electronic features also a Vector Network Analyser (VNA)
for measurement of the scattering parameters S12 (input from line L2 - output
from line L1), S31 and S32. From these scattering parameters it is possible to
derive the loaded quality factor $Q_{L}$, resonance frequency $f_{c}$ and
coupling $\beta$ of the tunable antenna. A diode noise source, having an
equivalent noise temperature of about $10^{4}$ K, can be fed to line L1 for
testing after being amplified in such a way to have an equivalent noise
temperature inside the microwave cavity slightly in excess of the
thermodynamic temperature. A microwave signal generator and a microwave
spectrum analyser are used for the measurement of the system noise temperature
as described below. All rf generators, the VNA and the spectrum analyser are
frequency locked to a GPS disciplined reference oscillator.

Following the figure, all components below the horizontal blue line sectioning
the 4K region are enclosed in a vacuum chamber immersed in a liquid helium
cryostat. A Ruthenium Oxide thermometer measures the temperature of the
cavity.

###  II.0.2 Data taking

In a data taking run performed over a two weeks period in June 2021 we
searched for an axion signal as a cavity excess-power in a small frequency
band about 10.353 GHz, i.e. about an axion mass of 42.8 $\mu$eV. As shown in
[39], the unloaded $Q$ factor of our microwave cavity is of several millions,
well in excess of the axion quality factor $Q_{a}=10^{6}$.

We decided to perform measurements for three different values of the cavity
loaded quality-factor $Q_{L}=Q_{0}/(1+\beta)$:

  1. 1)

$\beta\simeq 1$, i.e. $Q_{L}>Q_{a}$

  2. 2)

$\beta\simeq 6$, i.e. $Q_{L}\simeq Q_{a}$

  3. 3)

$\beta\geq 14$, i.e. $Q_{L}\ll Q_{a}$

The total data taking session comprised 8 sub-runs in regime 1, 33 sub-runs in
regime 2 and 11 sub-runs in regime 3. We performed the following steps for
each sub-run:

  1. a.

looking to the S32 spectra with a VNA, we moved the cavity frequency to the
desired value by acting (inserting) the sapphire triplet for tuning. Normally,
for each sub-run a shift of half the cavity linewidth with respect to the
previous sub-run was done.

  2. b.

we stored the S32 spectra for the resulting cavity-configuration. This
spectrum corresponds to a reflection-type measurement for the tunable port of
the cavity.

  3. c.

we injected along the line L1 a white noise produced with the amplified noise-
source and collected down-converted low-frequency I and Q spectra with the
ADC. For this run, thermal-input spectra were usually integrated for about 2-3
minutes.

  4. d.

we removed any input to the system and collected data with the ADC. This step
has been chosen to last 750 data blocks, for a total time of 3000 s. For each
sub-run both the cavity temperatures and the liquid-helium level in the
cryostat are recorded.

The total time needed for a single sub-run is about one hour. We performed
measurements only during the day, so that a typical day started with the
liquid helium refilling of the cryostat, followed by the magnet charging
lasting about 40 minutes. Considering that at the end of the day we needed to
discharge the magnet for safety reasons, in a typical day we recorded about
5-6 sub-runs.

In this paper we present the results obtained for the sub-runs with strong
cavity-coupling such that $Q_{L}\ll Q_{a}$ (regime 3), while the other data
will be subject of a different paper. Indeed, extracting the axion signal in a
regime where $Q_{L}\simeq Q_{a}$ or even higher poses a series of issues
regarding systematics that necessitate a dedicated study. For the measurements
in regime 3 we dealt with a loaded quality-factor $Q_{L}\sim 3\times 10^{5}$,
which is much larger than the typical values used by other running-haloscopes,
where such value was never in excess of $10^{5}$ (with the exemption of our
previous measurement [26]).

###  II.0.3 Noise temperature and gain

We measured the system noise-temperature at the beginning, at the end and in
the middle of the global data taking period. This procedure [40] consists in
measuring precisely the gains of the three lines L1, L2 and L3 from the point
A1 in Figure 1. The knowledge of the gains allows us to extract the system
noise-temperature from the measurement of the noise level at the output of
line L3. Gain measurement is obtained by feeding a calibrated power-level with
the signal generator either from L1 or L2, and reading the outputs at L3 for
inputs from both inputs lines or at L1, for an input from L2 only. The last
measurement is only possible when the tunable antenna has a significant
coupling to the cavity, therefore only at the cavity resonant frequency. On
the contrary, to exploit the cavity reflection, the measurement from L2 to L3
is done at a frequency just off the cavity resonance. Figure 2 shows the power
measured at the output of line L3 by feeding power into L1 (red points) or L2
(blue points). For each measure, we estimated the gains and the intercepts at
zero, $P_{0}$, with a fit with a linear function. The transmission coefficient
from input L2 to output L1 is done in a similar way. Combining these results,
we computed the gains from A1 to L1, L2, L3 as ($g_{1}$, $g_{2}$, $g_{3}$) =
(-48.3, -39.3, 52) dB, respectively. We computed the system noise-temperature
as

$T_{\rm sys}=\frac{P_{0}}{k_{B}Bg_{3}}$ (1)

obtaining $T_{\rm sys}=17.3\pm 1$ K. Here, $k_{B}$ is the Boltzmann constant
and $B$ the resolution bandwidth.

Figure 2: Power output at the line L3 with variable input at the lines L1
(red) and L2 (blue). For the L1 input a rf signal at the cavity resonance-
frequency is used, while for the L2 input the frequency is detuned by 1 MHz
from the cavity resonance. Measurements are performed with a spectrum analyser
taking 500 RMS averages of a 100 MHz window with a resolution bandwidth $B=1$
MHz.

This particular value of the noise temperature was obtained with the magnetic
field on at the end of the last day of run, at the end of the sessions with
regime 3. Previous measurements, performed with the magnet off, were in
agreement with this one. The cause of the large noise-temperature observed was
identified in a malfunction of the cryo HEMT exhibiting a quite high added
noise. After the run, we measured its added noise separately finding a value
about 10-12 K.

###  II.0.4 Raw data processing

As described in II.0.2, we measured, for each subrun, the values of the cavity
parameters by taking a reflection spectrum on the pick-up antenna, and a
transmission spectrum with a thermal-noise source feeding power into the
weakly coupled antenna. Parameters are extracted by fitting the spectra. A
standard Lorentzian line shape is used to fit the transmission spectrum:

$S31^{2}(\nu)=\frac{A}{2\pi}\frac{\Gamma}{(\nu-\nu_{c})^{2}+(\Gamma/2)^{2}}$
(2)

With this equation we extract from the fit the cavity resonance frequency
$\nu_{c}$ and the linewidth $\Gamma$, directly related to the loaded factor of
merit $Q_{L}=\nu_{c}/\Gamma$. $A$ is a normalization constant.

A modified reflection function is used for the reflection spectrum to take
into account some impedance mismatch between the cavity and the first-stage
amplifier:

$S32(\delta)=C\,\Bigg{|}\frac{\beta-1-iQ_{0}\delta}{\beta+1+iQ_{0}\delta}+ic\Bigg{|}$
(3)

where $C$ is a normalization constant, $\delta=\nu/\nu_{c}-\nu_{c}/\nu$, with
$\nu_{c}$ the cavity resonance frequency, $\nu$ the frequency, $c$ is a free
parameter related to the impedance mismatch, $Q_{0}$ is the cavity unloaded
quality factor. From this second fit we obtain the value of the coupling
$\beta$.

Figure 3: Reflection spectrum obtained with the VNA and fit with the function
(3). The fit results in the following values: $\beta=14.59\pm 0.01$,
$\nu_{c}=10353366689\pm 20$ Hz, $Q_{0}=5565000\pm 8000$, $c=0.0127\pm 0.0001$,
$C=3.6526\pm 0.0002$.

##  III Data Analysis and Results

By tuning the cavity resonance frequency ($\nu_{c}$) we acquired 11 different
dataset, one for each $\nu_{c}$ in the range $10.35337-10.35345$ GHz (Tab. 1).
For each dataset we calculated a power spectrum.

Table 1: cavity resonance frequency, quality factor and cavity-antenna coupling for each dataset. $\nu_{c}$ [GHz] | $Q_{L}$ | $\beta$
---|---|---
10.3533667 | 365730 | 14.59
10.3533711 | 337630 | 15.91
10.3533792 | 315100 | 17.00
10.3533874 | 288190 | 18.00
10.3533955 | 286620 | 17.87
10.3534036 | 284810 | 17.66
10.3534159 | 283410 | 17.61
10.3534150 | 354000 | 13.74
10.3534250 | 292510 | 16.20
10.3534354 | 290290 | 16.42
10.3534464 | 285760 | 17.25

In this section, we discuss the cumulative results obtained from the combined
spectra. For sake of simplicity, examples from a single dataset are reported
when necessary. The expected power generated by the axion conversion inside
the haloscope is given by [41, 34]:

$\begin{split}P_{a}=\left(\frac{g_{a\gamma\gamma}^{2}}{m_{a}^{2}}\,\hbar^{3}c^{3}\rho_{a}\right)\,\left(\frac{\beta}{1+\beta}\omega_{c}\frac{1}{\mu_{0}}B_{0}^{2}VC_{030}Q_{L}\right)&\\\
\times\left(\frac{1}{1+\left(2Q_{L}\,\Delta_{\omega}/\omega_{c}\right)^{2}}\right)\end{split}$
(4)

In the first set of parenthesis, $\rho_{a}\sim 0.45$ GeV/cm3 [42] is the local
dark matter density, $g_{a\gamma\gamma}$ is the coupling constant of the
axion-photon interaction, $m_{a}$ is the axion mass. The second set of
parenthesis contains the vacuum permeability $\mu_{0}$, the magnetic field
$B_{0}$ and the volume $V$ of the cavity. $\omega_{c}=2\pi\nu_{c}$ is the
resonance angular frequency of the cavity, $\beta$ and $Q_{L}$ are antenna
coupling and loaded quality factor as described above. $C_{030}$ is a
geometrical factor equal to about 0.028 for the TM030 mode of this cylindrical
dielectric cavity. In the third brackets, a Lorentzian function describes the
effect of the detuning $\Delta_{w}=\omega_{c}-\omega_{a}$ between the cavity
and an axion having angular frequency $\omega_{a}$. In presence of a signal
due to axion conversion a power excess would be observable in the residuals of
the power spectrum. The residuals are obtained subtracting a Savitzky-Golay
(SG) filter [43] of the fourth order to the cavity power spectrum. The dynamic
interval of the SG filter was optimized to 59.2 kHz (91 bins). For each
dataset, we applied the SG filter to a window
[$\nu_{c}-3\Gamma,\nu_{c}+3\Gamma$] of about 180 kHz, corresponding to six
linewidths $\Gamma$ and centered on the cavity resonance frequency $\nu_{c}$
as showed in Figure 4.

Figure 4: FFT cavity power spectrum (blue dots) and SG filter (black line).
$\nu_{c}=10.3534149$ GHz, $Q_{L}=354000$

In the laboratory frame, the axion signal is expected to have a width of about
10 KHz [6, 44]. With a power spectrum with bin width $\Delta\nu=651$ Hz we
expect the axion signal to be distributed over 16 consecutive bins. We
normalized the residuals of each dataset to the expected noise power
$\sigma_{\scriptscriptstyle\textup{Dicke}}$ calculated with the Dicke
radiometer equation [45]

$\sigma_{\scriptscriptstyle\textup{Dicke}}=k_{B}T_{\rm
sys}\sqrt{\Delta\nu/\Delta t}\,,$ (5)

where $T_{\rm sys}$, is the system noise-temperature, $\Delta\nu$ is the bin
width (651 Hz) and $\Delta t$ is the integration time (3000 s). The
distribution of the cumulative normalized-residuals from all the datasets is
shown in Fig. 5 along with a Gaussian fit, showing a standard deviation
compatible with 1.

Figure 5: Distribution of the cumulative residuals from each dataset
normalized to the $\sigma_{\scriptscriptstyle\textup{Dicke}}$.

We use the Least-Squares method to estimate the best value
$\hat{g}_{a\gamma\gamma}$ for the axion-photon coupling, by minimizing

$\chi^{2}=\sum_{\alpha=1}^{N_{\text{scan}}}\sum_{i=1}^{N_{\text{bin}}}\,\left[\frac{R^{(\alpha)}_{i}-S^{(\alpha)}_{i}(m_{a},g_{a\gamma\gamma}^{2})}{\sigma^{(\alpha)}_{\text{Dicke}}}\right]^{2}\,,$
(6)

where the $\alpha$ index runs over the $N_{\text{scan}}$ datasets taken with
different cavity resonant-frequencies, the index $i$ runs over the frequency
bins of each power spectrum, $R_{i\alpha}$ and $S_{i\alpha}$ are the residuals
and the expected power signals for the scan $\alpha$ and frequency bin $i$,
respectively. $S_{i\alpha}$ is calculated as the integral in the frequency
domain of equation (4) multiplied by the spectrum of the full standard halo
model distribution [44].

We express the expected power as
$S_{\alpha,i}(m_{a},g_{a\gamma\gamma}^{2})=g_{a\gamma\gamma}^{2}\,T_{\alpha,i}(m_{a})$,
and analytically minimize Eq. (6) by solving $\partial\chi^{2}/\partial
g^{2}_{a\gamma\gamma}=0$, and calculating the uncertainty according to the
formula $(\xi=g^{2}_{a\gamma\gamma})$:

$\frac{1}{\sigma^{2}_{\hat{\xi}}}=\frac{1}{2}\,\frac{\partial^{2}\chi^{2}}{\partial\xi^{2}}\,.$
(7)

Solving this equation, we get:
($\sum\sum\equiv\sum_{\alpha=1}^{N_{\text{scan}}}\sum_{i=1}^{N_{\text{bin}}}$)

$\overline{g^{2}}=\sigma^{2}(\overline{g^{2}})\,\left[\sum\sum\,\frac{R^{(\alpha)}_{i}\,T^{(\alpha)}_{i}(m_{a})}{(\sigma^{(\alpha)}_{\text{Dicke}})^{2}}\right]\,$
(8)

where $\overline{g^{2}}$ is the average squared coupling-constant that
accounts for the contributions of all the frequency bins of all the datasets,
and

$\sigma^{2}(\overline{g^{2}})=\sum\sum\,\left[\frac{T^{(\alpha)}_{i}(m_{a})}{\sigma^{(\alpha)}_{\text{Dicke}}}\right]^{2}$
(9)

is its variance. We repeated this procedure for different values of $m_{a}$
and calculated $\overline{g^{2}}$ and $\sigma(\overline{g^{2}})$ for axions
masses in the range $42.8210-42.8223$ $\mu$eV.

A candidate discovery requires the detection of a power excess larger than
5$\sigma$ above the noise, hence in the distribution of
$\overline{g^{2}}$/$\sigma(\overline{g^{2}})$. We did not find any candidate
(see Fig. 6) and the result is interpreted as an exclusion test in this axion-
mass range.

Figure 6: Histogram of the $\overline{g^{2}}$/$\sigma(\overline{g^{2}})$
distribution calculated using Eqs. (8) and (9). No excess above $5\sigma$ was
observed.

.

So far, we did not consider the efficiency of the SG filter in estimating the
axion signal. In order to quantify it, we run a Monte Carlo simulation where a
fake axion signal, with a known $g_{\text{injected}}^{2}$, is numerically
inserted in simulated power-spectra with different $\nu_{c}$. We used Eq. (8)
to estimate the $g^{2}$ for each injected signal
(g$g_{\text{calculated}}^{2}$), and determined the efficiency from the
relation between $g_{\text{calculated}}^{2}$ and $g_{\text{injected}}^{2}$. We
simulated the cavity power-spectra by adding random Gaussian noise (mean=0,
sigma=$\sigma_{Dicke}$, random values extracted between 0 and
$\pm\,\sigma_{\text{Dicke}}$, according to a Gaussian PDF) to the SG filters.
For a given axion mass, $m_{a}$, the estimation of the efficiency works as
follows: 1) for each data set we calculate a simulated spectrum; 2) a fake
axion signal with a known $g_{\text{injected}}^{2}$ is injected in the
simulated spectra; 3) Eq. (8) is used to compute g${}_{calculated}^{2}$; 4)
points 2 and 3 are repeated for different values of $g_{\text{injected}}^{2}$;
5) points 2-4 are repeated for a new set of simulated spectra in order to
increase the statistics. The output of this procedure is a distribution of
$g_{\text{calculated}}^{2}$ for each value of $g_{\text{injected}}^{2}$. The
relation between $g_{\text{injected}}^{2}$ and the mean of
$g_{\text{calculated}}^{2}$, calculated accounting for the contribution of all
the datasets, is shown in Figure 7 for axion proper frequency
faxion=10.35341562 GHz,
$g_{\text{injected}}^{2}=[0.01,1.755,6.502,14.25,25]\times 10^{-26}$ GeV-2,
and for 100 simulated spectra.

Figure 7: Relation between g${}_{\text{injected}}^{2}$ and the mean of the
$g_{\text{calculated}}^{2}$ distribution (blue points), along with the best
linear fit parameters. The belt represents the standard deviation of the
distribution of $g_{\text{calculated}}^{2}$ obtained after 100 simulations.

The distribution of $g_{\text{injected}}^{2}$ vs $g_{\text{calculated}}^{2}$
shows a linear relation with slope very close to 3 and an intercept different
from zero. These features are valid for all the axion masses not in the
immediate proximity of the edge of the power spectrum, where the slope
deviates considerably from 3. The slope of this linear relation is interpreted
as the inverse of the estimation efficiency. If we neglect the intercept, then
$g_{\text{calculated}}^{2}$/$g_{\text{injected}}^{2}\sim 1/3$, i.e. an
estimation efficiency of about 0.33. The intercept represents a contribution
to the $g_{\text{calculated}}^{2}$ given by the average noise of the simulated
spectra for a given ma.

Once we corrected the estimated $\overline{g^{2}}$ by the filter efficiency,
we calculated the limit on the axion-photon coupling with a 90% confidence
level as in [26], using a power constrained procedure for the
$\overline{g^{2}}$ that under fluctuates below $-\sigma$ [46].

Figure 8: The 90% single-sided C.L. upper limit for the axion coupling
constant $g_{a\gamma\gamma}$ as a function of the axion mass. The red solid
curve represents the expected limit in the case of no signal. The yellow
region indicates the QCD axion model band. Image realized using:
https://github.com/cajohare/AxionLimits

In Fig. 8 we show the calculated upper-limit
$g_{a\gamma\gamma}^{\scriptscriptstyle\textup{CL}}$ in the axion mass range
$42.821-42.8223$ $\mu$eV, i.e. a mass window of about 1.32 neV centered in
42.8216 $\mu$eV. The reference upper-limit of this analysis is the value at
the maximum sensitivity (the minimum spectrum reported in Fig. 8,
$g_{a\gamma\gamma}^{\scriptscriptstyle\textup{CL}}<0.731\times 10^{-13}$ GeV-1
at 90% C.L.

##  IV Conclusions

We reported the results of the search of galactic axions using an high-Q
dielectric haloscope. The investigated mass range is $42.8210-42.8223$
$\mu$eV. We set a limit for the axion-photon coupling a factor about 4 from
the axion-QCD band. We demonstrated the robustness of this detection approach
and the importance of working with a high Q cavity in the search for axions.
In fact, we managed to reach a sensitivity that almost touches the QCD band,
even though the equivalent thermal noise of our system was very high (about 17
K) due to an experimental setback, and the detection efficiency was
considered. In future experiments of this kind the sensitivity could be
further improved, reducing the overall noise and improving the thermalization
of the cavity.

## acknowledgments

DA VERIFICARE

We are grateful to E. Berto, A. Benato, and M. Rebeschini for the mechanical
work; F. Calaon and M. Tessaro for help with the electronics and cryogenics.
We thank G. Galet and L. Castellani for the development of the magnet power
supply, and M. Zago who realized the technical drawings of the system. We
deeply acknowledge the Cryogenic Service of the Laboratori Nazionali di
Legnaro for providing us with large quantities of liquid helium on demand.
This work was supported by INFN and partially supported by EU through FET Open
SUPERGALAX project, Grant N.863313

## References

  * Weinberg [1978] S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek [1978] F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Peccei and Quinn [1977] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Irastorza and Redondo [2018] I. G. Irastorza and J. Redondo, Prog. Part. Nucl. Phys. 102, 89 (2018).
  * Sikivie [1983] P. Sikivie, Phys. Rev. Lett. 51, 1415 (1983).
  * Sikivie [1985] P. Sikivie, Phys. Rev. D 32, 2988 (1985).
  * Braine  _et al._ [2020] T. Braine _et al._ , Phys. Rev. Lett. 124, 101303 (2020).
  * Du  _et al._ [2018] N. Du _et al._ , Phys. Rev. Lett. 120, 151301 (2018).
  * Boutan  _et al._ [2018] C. Boutan _et al._ , Phys. Rev. Lett. 121, 261302 (2018).
  * Bartram  _et al._ [2021] C. Bartram _et al._ , Phys. Rev Lett. 127, 261803 (2021).
  * Backes  _et al._ [2021] K. Backes _et al._ , Nature 590, 238 (2021).
  * Zhong  _et al._ [2018] L. Zhong _et al._ , Phys. Rev. D 97, 092001 (2018).
  * McAllister  _et al._ [2017] B. T. McAllister, G. Flower, E. N. Ivanov, M. Goryachev, J. Bourhill, and M. E. Tobar, Phys. Dark Univ. 18, 67 (2017).
  * Choi  _et al._ [2021] J. Choi, S. Ahn, B. Ko, S. Lee, and Y. K. Semertzidis, Nucl. Inst. Meth. Phys. Res. A 1013, 165667 (2021).
  * Lee  _et al._ [2020] S. Lee, S. Ahn, J. Choi, B. Ko, and Y. K. Semertzidis, Phys. Rev. Lett. 124, 101802 (2020).
  * Jeong  _et al._ [2020] J. Jeong, S. Youn, S. Bae, J. Kim, T. Seong, J. E. Kim, and Y. K. Semertzidis, Phys. Rev. Lett. 125, 221302 (2020).
  * Kwon  _et al._ [2021] O. Kwon _et al._ , Phys. Rev. Lett. 126, 191802 (2021).
  * Lee  _et al._ [2022] Y. Lee, B. Yang, H. Yoon, M. Ahn, H. Park, B. Min, D. Kim, and J. Yoo, Phys. Rev. Lett. 128, 241805 (2022).
  * Grenet  _et al._ [2021] T. Grenet, R. Ballou, Q. Basto, K. Martineau, P. Perrier, P. Pugnat, J. Quevillon, N. Roch, and C. Smith, arXiv:2110.14406 (2021).
  * Álvarez Melcón  _et al._ [2020] A. Álvarez Melcón _et al._ , JHEP 07 (2020), 084.
  * Álvarez Melcón  _et al._ [2018] A. Álvarez Melcón _et al._ , JCAP 05 (2018), 040.
  * Álvarez Melcón  _et al._ [2021] A. Álvarez Melcón _et al._ , JHEP 10 (2021), 075.
  * Chang  _et al._ [2022] H. Chang _et al._ , arXiv:2205.05574 (2022).
  * Alesini  _et al._ [2019a] D. Alesini _et al._ , Phys. Rev. D 99, 101101 (2019a).
  * Alesini  _et al._ [2021a] D. Alesini _et al._ , Phys. Rev. D 103, 102004 (2021a).
  * Barbieri  _et al._ [2017] R. Barbieri, C. Braggio, G. Carugno, C. S. Gallo, A. Lombardi, A. Ortolan, R. Pengo, G. Ruoso, and C. C. Speake, Phys. Dark Univ. 15, 135 (2017).
  * Crescini  _et al._ [2018] N. Crescini _et al._ , Eur. Phys. J. C 78, 1 (2018).
  * Crescini  _et al._ [2020] N. Crescini _et al._ , Phys. Rev. Lett. 124, 171801 (2020).
  * Gatti  _et al._ [2018] C. Gatti _et al._ , arXiv:1811.06754 (2018).
  * Alesini  _et al._ [2019b] D. Alesini _et al._ , arXiv:1911.02427 (2019b).
  * Caldwell  _et al._ [2017] A. Caldwell _et al._ , Phys. Rev. Lett. 118, 091801 (2017).
  * Lawson  _et al._ [2019] M. Lawson, A. J. Millar, M. Pancaldi, E. Vitagliano, and F. Wilczek, Phys. Rev. Lett. 123, 141802 (2019).
  * Al Kenany  _et al._ [2017] S. Al Kenany _et al._ , Nucl. Instr. Meth. Phys. Res. A 854, 11 (2017).
  * Di Gioacchino  _et al._ [2019] D. Di Gioacchino _et al._ , IEEE Trans. App. Sup. 29, 1 (2019).
  * Ahn  _et al._ [2019] D. Ahn, O. Kwon, W. Chung, W. Jang, D. Lee, J. Lee, S. W. Youn, D. Youm, and Y. K. Semertzidis, arXiv:1904.05111 (2019).
  * Alesini  _et al._ [2020] D. Alesini _et al._ , Rev. Sci. Instr. 91, 094701 (2020).
  * Alesini  _et al._ [2021b] D. Alesini _et al._ , Nucl. Instr. Meth. Phys. Res. A 985, 164641 (2021b).
  * Di Vora  _et al._ [2022] R. Di Vora _et al._ , Phys. Rev. Applied 17, 054013 (2022).
  * Braggio  _et al._ [2022] C. Braggio _et al._ , arXiv:2205.02053 (2022).
  * Brubaker  _et al._ [2017] B. M. Brubaker _et al._ , Phys. Rev. Lett. 118, 061302 (2017).
  * Zyla  _et al._ [2020] P. A. Zyla _et al._ , (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 10.1093/ptep/ptaa104 (2020), 083C01, https://academic.oup.com/ptep/article-pdf/2020/8/083C01/34673722/ptaa104.pdf .
  * Savitzky and Golay [1964] A. Savitzky and M. J. Golay, Anal. Chem. 36, 1627 (1964).
  * Turner [1990] M. S. Turner, Phys. Rev. D 42, 3572 (1990).
  * Dicke [1946] R. H. Dicke, The measurement of thermal radiation at microwave frequencies, in  _Classics in Radio Astronomy_ (Springer, Dordrecht, 1946) pp. 106–113.
  * Cowan  _et al._ [2011] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, arXiv:1105.3166 (2011).


11institutetext: Department of Engineering Science, University of Oxford, UK

# Principled Ultrasound Data Augmentation for Classification of Standard
Planes

Lok Hin Lee  Yuan Gao and J. Alison Noble

###### Abstract

Deep learning models with large learning capacities often overfit to medical
imaging datasets. This is because training sets are often relatively small due
to the significant time and financial costs incurred in medical data
acquisition and labelling. Data augmentation is therefore often used to expand
the availability of training data and to increase generalization. However,
augmentation strategies are often chosen on an ad-hoc basis without
justification. In this paper, we present an augmentation policy search method
with the goal of improving model classification performance. We include in the
augmentation policy search additional transformations that are often used in
medical image analysis and evaluate their performance. In addition, we extend
the augmentation policy search to include non-linear mixed-example data
augmentation strategies. Using these learned policies, we show that principled
data augmentation for medical image model training can lead to significant
improvements in ultrasound standard plane detection, with an an average
F1-score improvement of 7.0% overall over naive data augmentation strategies
in ultrasound fetal standard plane classification. We find that the learned
representations of ultrasound images are better clustered and defined with
optimized data augmentation.

###### Keywords:

Data augmentation Fetal Ultrasound.

##  1 Introduction

The benefits of data augmentation for training deep learning models are well
documented in a variety of tasks, including image recognition [20, 5, 23] and
regression problems [12, 22]. Data augmentation acts to artificially increase
the size and variance of a given training dataset by adding transformed copies
of the training examples. This is particularly evident in the context of
medical imaging, where data augmentation is used to combat class imbalance
[10], increase model generalization [11, 18], and expand training data [28,
9]. This is usually done with transformations to the input image that are
determined based on expert knowledge and cannot be easily transferred to other
problems and domains. In ultrasound, this usually manifests as data
augmentation strategies consisting of small rotations, translations and
scalings [1]. However, whilst it is appealing to base augmentation strategies
on “expected” variations in input image presentation, recent work has found
that other augmentation strategies that generate “unrealistic looking”
training images [19, 24] have led to improvements in generalization
capability. There has therefore been great interest in developing data
augmentation strategies to automatically generate transformations to images
and labels that would lead to the greatest performance increase in a neural
network model. In this paper, inspired by the RandAugment [7] augmentation
search policy, we automatically look for augmentation policies that outperform
standard augmentation strategies in ultrasound imaging based on prior
knowledge and extend our algorithm to include mixed-example data augmentation
[24] in the base policy search. We evaluate the proposed method on second-
trimester fetal ultrasound plane detection, and find that a randomly
initialized network with our augmentation policy achieves performance
competitive with methods that require external labelled data for network pre-
training and self-supervised methods. We also evaluate our method on a fine-
tuning a pre-trained model, and find that using an optimized augmentation
policy during training improves final performance.

###  1.0.1 Contributions:

Our contributions are three fold: 1) We investigate the use of an augmentation
search policy with hyperparameters that does not need expensive reinforcement
learning policies and can be tuned with simple grid search; 2) We extend this
augmentation search policy to combinations that include mixed-example based
data augmentation and include common medical imaging transformations; 3) We
explain the performance of optimal augmentation strategies by looking at their
affinity, diversities and effect on final model performance.

###  1.0.2 Related Work

Medical image datasets are difficult and expensive to acquire. There has
therefore been previous work that seeks to artificially expand the breadth of
training data available in medical image classification [11, 16], segmentation
[21, 10] and regression [9].

Original Data Manipulation: Zeshan et al. [16] evaluate the performance of
eight different affine and pixel level transformations by training eight
different CNNs for predicting mammography masses and find that ensembling the
trained models improves the classification performance significantly. Nalepa
et al. [21] elastically deform brain MRI scans using diffeomorphic mappings
and find that tumour segmentation is improved. However, in the above works,
specific augmentations and parameters are selected arbitrarily and are task
and modality dependent. In contrast, we propose an automated augmentation
policy search method that can out perform conventional medical imaging
augmentation baselines.

Artificial Data Generation: Florian et al. [9] generates new training samples
in by linearly combining existing training examples in regression. Models
trained to estimate the volume of white matter hyperintensities had
performance comparable to networks trained with larger datasets. Zach et
al.[10] also linearly combine training examples and target labels linearly
inspired by mix-up [27] but focus on pairing classes with high and low
incidence together, which was found to be beneficial for tasks with high class
imbalance such as in brain tumor segmentation. Maayan et al. [11] train a
conditional generative adversarial network (cGAN) to generate different types
of liver lesions and use the synthesized samples to train a classification
network. Dakai et al. [18] use a cGAN to synthesize 3D lung nodules of
different sizes and appearances at multiple locations of a lung CT scan. These
generated samples were then used to finetune a pretrained lung nodule
segmentation network that improved segmentation of small peripheral nodules.
However, cGAN-based methods are difficult to train and have significant
computational costs during augmentation.

Automated Augmentation Policy Search: There are augmentation policy search
methods in the natural image analysis [6, 14, 19] that learn a series of
transformations which are parameterized by their magnitudes and probability.
However, these searches are expensive, and cannot be run on the full training
dataset as the hyperparameter search for each transformation require
significant computational resources. RandAugment (RA) [7] finds that
transformations can have a shared magnitude and probability of application and
achieve similar performance, without expensive reinforcement learning.
However, RA is limited to single-image transformations. We therefore explore
the use of an extended RA policy search with additional transformations that
are more specific to medical imaging, and expand its capabilities to include
mixed-example image examples to include artificial data in model training.

##  2 Methods

Figure 1: Overview of our proposed learning framework.

In this section we describe our proposed framework for augmentation policy
search, depicted in Figure 1 which consists of three key procedures i) data
generation, ii) data augmentation, iii) policy searching and interpretation.

###  2.0.1 Mixed-Example Data Augmentation

Figure 2: The procedure for non-linear mixed-example data augmentation using
an image pair and the final artificial mixed-example image.

The original dataset $D=\\{(X_{i},Y_{i})\\}$ consists of a series of $i$
ultrasound frames $X$ and their associated classes $Y$. We first generate a
paired dataset
$D_{paired}=\\{(x_{1},x_{2})_{\frac{i}{2}},(y_{1},y_{2})_{\frac{i}{2}}\\}$ by
pairing examples from different classes. Examples of artificial data are then
generated using non-linear methods [25], which are found to be more effective
than linear intensity averaging (mix-up)[27]. As illustrated in Figure 2,
instead of pixel-wise averaging, the bottom $\lambda_{1}$ fraction of image
$x_{1}$ is vertically concatenated with the top $1-\lambda_{1}$ fraction of
image $x_{2}$. Similarly, the right $\lambda_{2}$ fraction of image $x_{1}$ is
horizontally concatenated with the left $1-\lambda_{2}$ fraction of image
$x_{2}$. After the concatenations, the resulted images are combined to produce
an image $\tilde{x}$ in which the top-left is from $x_{1}$, the bottom right
is from $x_{2}$, and the top-right and bottom-left are mixed between the two.
Moreover, instead of linear pixel averaging, we treat each image as a zero-
mean waveform and normalize mixing coefficients with image intensity energies
[26]. Formally, given initial images $x_{1,2}$ with image intensity means and
standard deviations of $\mu_{1,2}$ and $\sigma_{1,2}$, the generated
artificial mixed-example image $\tilde{x}$ is:

$\tilde{x}=\begin{cases}x_{1}(i,j)-\mu_{1}&\quad\text{if}\;i\leq\lambda_{1}H\;\text{and}\;j\leq\lambda_{2}W\\\
\frac{c}{\phi}[x_{1}(i,j)-\mu_{1}]+\frac{1-c}{\phi}[x_{2}(i,j)-\mu_{2}]&\quad\text{if}\;i\leq\lambda_{1}H\;\text{and}\;j>\lambda_{2}W\\\
\frac{1-c}{\phi}[x_{1}(i,j)-\mu_{1}]+\frac{c}{\phi}[x_{2}(i,j)-\mu_{2}]&\quad\text{if}\;i>\lambda_{1}H\;\text{and}\;j\leq\lambda_{2}W\\\
x_{2}(i,j)-\mu_{2}&\quad\text{if}\;i>\lambda_{1}H\;\text{and}\;j>\lambda_{2}W\\\
\end{cases}$

where $c$ is the mixing coefficient
$(1+\frac{\sigma_{1}}{\sigma_{2}}\cdot\frac{1-\lambda_{3}}{\lambda_{3}})^{-1}$
and $\phi$ is the normalization term defined as $\sqrt{c^{2}+(1-c)^{2}}$. The
row index and column index is represented by $i,j$ and the height and width of
the images are represented by $H,W$.

We sample $\lambda_{1,2,3}\sim Beta(m/10,m/10)$ where $m$ is a learnt
hyperparameter varied from 0-10. As m approaches 10, $\lambda$ values are more
uniformly distributed across 0-1 and artificial images are more interpolated.
The ground truth label after interpolation is determined by the mixing
coefficients and can be calculated with:

$\tilde{y}=(\lambda_{3}\lambda_{1}+(1-\lambda_{3})\lambda_{2})y_{1}+(\lambda_{3}(1-\lambda_{1})+(1-\lambda_{3})(1-\lambda_{2}))y_{2}$

###  2.0.2 Original Data Augmentation

Augmentation transformations are then applied to the mixed images. Inspired by
[7], we do not learn specific magnitudes and probabilities of applying each
transformation in a given transformation list. Each augmentation policy is
instead defined only by $n$, which is the number of transformations from the
list an image undergoes, and $m$, which is the magnitude distortion of each
transformation. Note that $m$ is a shared hyperparameter with the mixed-
example augmentation process.

Figure 3: Examples of how each transformation and augmentation policy affect
input images. Each color represents a transformation, and augmented image is
transformed using a number of transformations (n) at a magnitude of (m).
$\dagger$: our additional transformations.

We investigate the inclusion in the transformation list transformations
commonly used in ultrasound image analysis augmentation: i) grid distortions
and elastic transformation [4] and ii) speckle noise [2]. The transformation
list then totals 18 transformations, examples of which can be seen in Figure
3.

###  2.0.3 Optimization

We define $f$ and $\theta$ as a convolutional neural network (CNN) and its
parameters. As depicted in Figure 1, we train a CNN with the augmented mini-
batch data $\tilde{x}^{i}$ and obtain the predicted output class scores
$f_{\theta}(\tilde{x}^{i})$. These are converted into class probabilities
$p(\tilde{x}^{i})$ with the softmax function. The KL-divergence between
$f_{\theta}(\tilde{x}^{i})$ and $\tilde{y}^{i}$ is then minimized with back-
propagation and stochastic gradient descent

$L=\frac{1}{B}D_{KL}(\tilde{y}^{i}\parallel
p(\tilde{x}^{i}))=\frac{1}{B}\sum_{i=1}^{B}\sum_{j=1}^{C}\tilde{y}_{j}^{i}log\frac{\tilde{y}_{j}^{i}}{\left\\{p(\tilde{x}^{i})\right\\}_{j}}$

where $B$ is the batch size, $C$ is the number of classes and $L$ is the loss.

Due to the limited search space, the hyperparameters $n$ and $m$ that produce
the optimum classification performance can then be found using grid search as
seen in Figure 1. The best performing $m,n$ tuple is then used during final
model evaluation.

###  2.0.4 Quantifying Augmentation Effects

Next, we investigate how augmentation improves model generalization and
quantify how different augmentation policies affect augmented data
distributions and model performance. We adopt a two dimensional metric -
affinity and diversity [13] to do this. Affinity quantifies the distribution
shift of augmented data with respect to the unaugmented distribution captured
by a baseline model; the diversity quantifies complexity of the augmented
data. Given training and validation datasets, $D_{t}$ and $D_{v}$, drawn from
the original dataset $D$, we can generate an augmented validation dataset
$D(m,n)_{v}^{{}^{\prime}}$ derived from $D_{v}$ using $m,n$ as hyperparameters
for the augmentation policy. The affinity $A$ for this augmentation policy is
then:

$A=\mathbb{E}[L({D_{v}^{{}^{\prime}}})]-\mathbb{E}[L(D_{v})]$

where $\mathbb{E}[L(D)]$ represents the expected value of the loss computed on
the dataset $D$ loss of a model trained on $D_{t}$.

The diversity, $D$, of the augmentation policy $a$ is computed on the
augmented training dataset $D_{t}^{{}^{\prime}}$ with respect to the expected
final training loss, $L_{t}$, as:

$D=\mathbb{E}[L(D_{t}^{{}^{\prime}})]$

Intuitively, the greater the difference in loss between an augmented
validation dataset and the original dataset on a model trained with
unaugmented data, the greater the distribution shift of the augmented
validation dataset. Similarly, the greater the final training loss of a model
on augmented data, the more complexity and variation there is in the final
augmented dataset.

##  3 Experiments and Results

We use a clinically acquired dataset consisting of ultrasound second-trimester
fetal examinations. A GE Voluson E8 scanner was used for ultrasound image
acquisition. For comparison with previous work [8, 17], fetal ultrasound
images were labelled into 14 categories. Four cardiac view classes (4CH, 3VV,
LVOT, RVOT) corresponding to the four chamber view, three vessel view, left
and right ventricle tracts respectively; the brain transcerebellar and
transventricular views (TC, TV); two fetal spine sagittal and coronal views
(SpineSag, SpineCor); the kidney, femur, abdominal circumference standard
planes, profile view planes and background images. The standard planes from
135 routine ultrasound clinical scans were labelled, and 1129 standard plane
frames were extracted. A further 1127 background images were also extracted
and three-fold cross validation was used to verify the performance of our
network.

###  3.0.1 Network Implementation

The performance of the SE-ResNeXt-50 [15] backbone is well validated on
natural images and therefore used. Networks were trained with an SGD optimizer
with learning rate of $10^{-3}$, a momentum of 0.9 and a weight decay of
$10^{-4}$. Networks were trained for a minimum of 150 epochs, and training was
halted if there was 20 continuous epochs without improvement in validation
accuracy. Models were implemented with PyTorch and trained on a NVIDIA GTX
1080 Ti. Random horizontal and vertical flipping were used in all RA policies
as a baseline augmentation. Models were trained with a batch size of 50. We
evaluated the performance of networks trained with augmentation policies with
values of $m,n$ where $m,n=\\{1,3,5,7,9\\}$ and used a simple grid search for
augmentation strategies to find optimal $m,n$ values.

###  3.0.2 Results on CNNs with Random Initialization

Table 1: Results for standard plane detection (mean $\pm$ std $\%$). The best
performing augmentation strategies are marked in bold for each metric.

| Random Initialization | Initialized with external data
---|---|---
| No Aug. | SN Pol. |  RA [7] |  Mix. RA (ours) |  Siam. Init. [17] |  Saliency [8] |  SonoNet [8, 3] |  SonoNet + Mix.RA (ours)
Precision | 56.5$\pm$1.2 | 70.4$\pm$2.3 | 74.7$\pm$1.8 | 75.1$\pm$1.8 | 75.8$\pm$1.9 | 79.5$\pm$1.7 | 82.3$\pm$1.3 | 86.3$\pm$1.3
Recall | 55.1$\pm$1.2 | 64.9$\pm$1.6 | 72.2$\pm$2.3 | 73.4$\pm$1.9 | 76.4$\pm$2.7 | 75.1$\pm$3.4 | 87.3$\pm$1.1 | 85.1$\pm$1.3
F1-Score | 55.4$\pm$1.2 | 67.0$\pm$1.3 | 72.8$\pm$1.8 | 74.0$\pm$1.8 | 75.7$\pm$2.0 | 76.6$\pm$2.6 | 84.5$\pm$0.9 | 85.4$\pm$1.5

The effectiveness of our mixed-example augmentation search policy algorithm
(Mix. RA) on SE-ResNeXt-50 models that are randomly initialized is compared
with models trained with the baseline RandAugment (RA) augmentation search
policy; a commonly used augmentation strategy (SN Pol.) in line with that
found in [8], where images are augmented with random horizontal flipping,
rotation $\pm 10^{\circ}$, aspect ratio changes $\pm 10\%$, cropping and
changing brightness $\pm 25\%$ and image cropping $95-100\%$; and no
augmentation (No. Aug.).

Figure 4: Confusion matrix for Mix. RA (left) and the difference in precision
between Mix. RA and SN Pol.

From Table 1 we can see that the proposed method Mix. RA outperforms all other
augmentation methods on every metric with random network initialization,
including the baseline RA augmentation policy search.

To better understand how Mix. RA outperforms naive augmentation, we show the
confusion matrix for the best performing model Mix. Aug and the difference in
confusion matrix between it and naive augmentation SN Pol. We find that in
general, heart plane classification is improved with a mean increase in macro
F1-Score of $4.0\%$. Other anatomical planes with the exception of the femur
plane also show consistent increases in performance with movement of
probability mass away from erroneously classified background images to the
correct classes suggesting the model is able to recognize greater variation in
each anatomical class.

Figure 5: t-SNE embeddings for different augmentation policies. The boxes
represent the regions from which example images are taken from. The blue and
purple boxes contain examples of the (SpineCor, SpineSag, and BG) and (3VT,
RVOT, BG) classes respectively taken from the highlighted positions in the
t-SNE embedding. Best viewed in color.

The t-SNE embeddings of the penultimate layer seen in Figure 5 can also be
used to visualize the differences in feature spaces in trained networks with
different augmentation policies. Compared to the model trained with no
augmentation, our best performing policy leads to a better separation of the
abdominal and profile standard planes from the background class as well as
clearer decision boundaries between anatomical classes. The two brain views
(TC, TV) and the demarcation of the boundary between the kidney view and
abdominal plane view is also better defined.

Between the best performing policy $m,n=(5,3)$ and an underperforming policy
$m,n=(7,7)$, we find that profile planes are better separated from the
background class and the abdominal planes better separated from the kidney
views, which suggests that the optimum $m,n$ value increases network
recognition of salient anatomical structures. However, in all three cases, the
cardiac views remain entangled. This can be attributed to the difficulty of
the problem, as even human expert sonographers cannot consistently
differentiate between different cardiac standard plane images. We also find
that the background class also contains examples of the anatomies in each
class, but in sub-optimal plane views, which leads to confusion during
classification. This difficulty is illustrated in example background images in
Figure 5 where the heart and spine are visible in the BG class.

###  3.0.3 Pre-trained Networks

We also compare our work to methods where networks were initialized with
external data as seen in the right of Table 1. Baseline methods of self-
supervised pre-training using video data [17] (Siam. Init.), multi-modal
saliency prediction [8] (Saliency) and Sononet (Sononet) [3] were used to
initialize the models and the models fine-tuned on our dataset. Using our
augmentation policy during fine-tuning of a pre-trained SonoNet network
further increased the performance of standard plane classification with an
increase in final F1-score of 0.9% when $m,n=(5,1)$. This reduction in optimum
transformation magnitude may be due to the change in network architecture from
SE-ResNeXt-50 to a Sononet, as the smaller Sononet network may not be able to
capture representations the former is able to. Furthermore, we find that
augmentation policy with a randomly initialized network Mix. RA approaches the
performance of the Siam. Init. and Saliency pre-trained networks. This is
despite the fact that the Siam. Init. requires additional 135 US videos for
network self-supervised initialization, and Saliency required external multi-
modal data in the form of sonographer gaze.

###  3.0.4 Ablation Study

To better understand the individual contributions to the Mix. RA augmentation
search policy, we show the results of an ablation study on the components of
Mix. RA in Table 2.

Table 2: Ablation study on the individual components of our Mix. RA policy
search algorithm on training of a randomly initialized CNN for ultrasound
standard plane detection. All metrics are macro-averages due to the class
imbalance. The Linear Mix. RA is included as a baseline mixed-example
augmentation strategy.

| No Aug. | SN Pol. | RA |  | RA + Speckle
---
\raisebox{-.9pt} {1}⃝
| RA + Deform.
---
\raisebox{-.9pt} {2}⃝
| Ext. RA
---
\raisebox{-.9pt} {1}⃝ $+$ \raisebox{-.9pt} {2}⃝
| Linear Mix. RA
---
\raisebox{-.9pt} {1}⃝ $+$ \raisebox{-.9pt} {2}⃝
| Non-Linear Mix. RA
---
\raisebox{-.9pt} {1}⃝ $+$ \raisebox{-.9pt} {2}⃝
Precision | 56.5$\pm$1.2 | 70.4$\pm$2.3 | 73.9$\pm$1.7 | 74.6$\pm$1.7 | 74.6$\pm$1.7 | 74.0$\pm$2.4 | 74.6$\pm$1.6 | 75.1$\pm$1.8
Recall | 55.1$\pm$1.2 | 64.9$\pm$1.6 | 72.9$\pm$2.0 | 72.1$\pm$1.8 | 72.9$\pm$1.8 | 74.5$\pm$1.3 | 73.2$\pm$1.5 | 73.4$\pm$1.9
F1-Score | 55.4$\pm$1.2 | 67.0$\pm$1.3 | 72.8$\pm$1.8 | 72.9$\pm$1.7 | 73.2$\pm$1.2 | 73.6$\pm$1.3 | 73.7$\pm$1.6 | 74.0$\pm$1.8

It can be seen that both including Speckle noise transformations and
deformation (Grid, Elastic) transformations lead to increased classifier
performance for standard plane classification of +0.1% and 0.3% respectively
with further improvement when both are combined together with Ext. RA. We find
that both RA and Ext. RA had an optimal $m,n=(5,3)$, suggesting that the
magnitude ranges for our additional transformations are well matched to the
original transformation list. This performance increase is further boosted
when mixed-example augmentation is introduced on top of Ext. RA, with non-
linear mixed-example augmentations outperforming a linear mix-up based method.

###  3.0.5 Affinity and Diversity

Figure 6: Affinity and diversity metrics for RA, Ext. RA and Mix. RA
augmentation policy search algorithms. Each point represents a (m, n) value in
hyperparameter space. Best performing $m,n$ values are highlighted in red for
each policy search method. Color represents final F1-Score on a randomly
initialized CNN.

The affinity and diversity of the augmentation policies is shown in Fig. 6. We
find that there exists a “sweet spot” of affinity and diversity using non-
mixed class augmentation strategies at an affinity distance of $\sim$3.8 and
diversity of $\sim$0.25 which maximized model performance, corresponding to
$m,n=(5,3)$. At high $m,n$ values, affinity distance is too high and the
distribution of the augmented data is too far away from the validation data,
decreasing model performance. However, at low $m,n$ values, the diversity of
the augmented data decreases and the model sees too little variation in input
data.

It can also be seen that the Mix. RA augmented dataset showed a reduced
affinity distance to the original dataset than the Ext. RA dataset at the same
$m,n=(5,3)$ value, implying that our proposed transforms shifts augmented
images to be more similar to the original images. Moreover, using a mixed-
example data augmentation strategy drastically increased diversity for any
given value of data distribution affinity, which improved final model
performance. The best performing mixed-example augmentation policy $m,n=(3,3)$
reduced the magnitude of each transformation compared to the optimal non-
linear augmentation policy. This suggests that mixed-example augmentation acts
to increase the diversity of training images which reduces the magnitude
required during further processing.

##  4 Conclusion

The results have shown that we can use a simple hyper-parameter grid search
method to find an augmentation strategy that significantly outperforms
conventional augmentation methods. For standard plane classification, the best
performing augmentation policy had an average increase in F1-Score of 7.0%
over that of a standard ultrasound model augmentation strategy. Our
augmentation policy method is competitive with the Siam. Init. [17] despite
the latter needing additional external data for pre-training. Our method also
improves the performance of a Sononet pre-trained model when fine-tuned using
our augmentation policy search method. From t-SNE plots and confusion matrix
differences, we can see that the performance increase is from better
classification of background-labelled planes, despite the difficulty in heart
plane classification. It should be noted that a large degree of
misclassification was due to standard planes being mis-classified into
background images or vice-versa, and qualitative evaluation of t-SNE clusters
show that this was due to background images being examples of class anatomies.
The ablation study also shows that our additional transformations improve
model performance, and non-linear mixed-example augmentation further improves
classification. The evaluation using affinity and diversity indicate that the
hyperparameter search involves a trade-off between diversity and affinity. We
find that using non-linear mixed-class data augmentation before
transformations drastically increases diversity without further increasing
affinity distance between the training data and augmented data, which helps
explain the increase in model performance. In conclusion, we have shown that
our augmentation policy search method outperforms standard manual choice of
augmentation policy. The augmentation policy search method presented does not
have any inference-time computational cost, and has the potential to be
applied in other medical image settings where training data is insufficient
and costly to acquire.

## References

  * [1] Asaduz, Z., Sang, H.P., et al.: Generative approach for data augmentation for deep learning-based bone surface segmentation from ultrasound images. ICARS (2020)
  * [2] Bargsten, L., Schlaefer, A.: Specklegan: a generative adversarial network with an adaptive speckle layer to augment limited training data for ultrasound image processing. ICARS (2020)
  * [3] Baumgartner, C.F., et al.: SonoNet: Real-Time Detection and Localisation of Fetal Standard Scan Planes in Freehand Ultrasound. IEEE TMI (2017)
  * [4] Buslaev, A., Iglovikov, V.I., et al.: Albumentations: Fast and flexible image augmentations. MDPI (2020)
  * [5] Cecilia, S., Michael, J.D.: Improved mixed-example data augmentation. In: IEEE Winter Conference on Applications of Computer Vision (WACV) (2019)
  * [6] Cubuk, E.D., Zoph, B., Mane, D., Vasudevan, V., Le, Q.V.: Autoaugment: Learning augmentation strategies from data. In: CVPR (2019)
  * [7] Cubuk, E.D., Zoph, B., Shlens, J., Le, Q.V.: Randaugment: Practical automated data augmentation with a reduced search space. In: CVPR (2020)
  * [8] Droste, R., et al.: Ultrasound Image Representation Learning by Modeling Sonographer Visual Attention. IPMI (2019)
  * [9] Dubost, F., Bortsova, G., et al.: Hydranet: Data Augmentation for Regression Neural Networks. MICCAI (2019)
  * [10] Eaton-Rosen, Z., Bragman, F., et al.: Improving data augmentation for medical image segmentation. In: MIDL (2018)
  * [11] Frid-Adar, M., Diamant, I., et al.: Gan-based synthetic medical image augmentation for increased cnn performance in liver lesion classification. Neurocomputing (2018)
  * [12] Gan, Z., Henao, R., Carlson, D., Carin, L.: Learning deep sigmoid belief networks with data augmentation. Proceedings of Machine Learning Research (PMLR (2015)
  * [13] Gontijo-Lopes, R., Smullin, S.J., Cubuk, E.D., Dyer, E.: Affinity and diversity: Quantifying mechanisms of data augmentation (2020)
  * [14] Ho, D., Liang, E., et al.: Population based augmentation: Efficient learning of augmentation policy schedules. In: ICML (2019)
  * [15] Hu, J., et al.: Squeeze-and-Excitation Networks. IPAMI (2020)
  * [16] Hussain, Z., Gimenez, F., Yi, D., Rubin, D.: Differential data augmentation techniques for medical imaging classification tasks. AMIA (2017)
  * [17] Jiao, J., et al.: Self-supervised representation learning for ultrasound video. In: IEEE 17th International Symposium on Biomedical Imaging (2020)
  * [18] Jin, D., Xu, Z., et al.: CT-realistic lung nodule simulation from 3D conditional generative adversarial networks for robust lung segmentation. in MICCAI (2018)
  * [19] Lim, S., Kim, I., Kim, T., Kim, C., Kim, S.: Fast autoaugment. in Advances in Neural Information Processing Systems (NeurIPs) (2019)
  * [20] Luke, T., Geoff, N.: Improving deep learning using generic data augmentation. IEEE Symposium Series on Computational Intelligence (2018)
  * [21] Nalepa, J., et al.: Data augmentation via image registration. In: ICIP (2019)
  * [22] Ohno, H.: Auto-encoder-based generative models for data augmentation on regression problems (2019)
  * [23] Ryo, T., Takashi, M.: Data augmentation using random image cropping and patches for deep CNNs. IEEE TCSVT (2020)
  * [24] Summers, C., Dinneen, M.J.: Improved mixed-example data augmentation. In: IEEE Winter Conference on Applications of Computer Vision (WACV) (2019)
  * [25] Summers, C., Dinneen, M.J.: Improved mixed example data augmentation. In: IEEE Winter Conference on Applications of Computer Vision (WACV) (2019)
  * [26] Tokozume, Y., Ushiku, Y., Harada, T.: Between-class learning for image classification. In: CVPR (2018)
  * [27] Zhang, H., Cisse, M., Dauphin, Y.N., Lopez-Paz, D.: mixup: Beyond empirical risk minimization. In: ICLR (2018)
  * [28] Zhao, A., Balakrishnan, G., et al.: Data augmentation using learned transformations for one-shot medical image segmentation. CVPR (2019)


Furthermore, the (proof) complexity $f$ of A is defined to be

$f(n)\coloneqq\max_{\begin{subarray}{c}\enskip\mathsf{A}\in\textbf{{A}}\colon\\\
\lvert\mathsf{A}\rvert\leq
n\end{subarray}}\;\min_{\begin{subarray}{c}p\in\mathbb{N}\colon\\\ p\text{
$\operatorname{\mathit{pv}}$-decides }\mathsf{A}\end{subarray}}\lvert p\rvert$

In plain English, the complexity of a proof problem the minimum number of
symbols required to decide every sentence up to length $n$.

Onto the result.

[Proof Hierarchy Theorem] Let T be a consistent and computably enumerable
theory which represents every computable function. Let $g$ be computable,
monotone, and $g(n)=\Omega(n)$. And let $f(n+1)=o(g(n))$. Then for any
canonical proof verifier $\operatorname{\mathit{pv}}$ there exists a proof
problem with complexity $O(g(n))$ but not $O(f(n))$.

###### Proof.

Consider the following algorithm.
$Y(\text{\T1\guilsinglleft$f$\T1\guilsinglright},x)$
$m\coloneqq\bigl{\lvert}\neg\,f(\raise
1.46208pt\hbox{\footnotesize$\ulcorner$}\hbox{$\text{\T1\guilsinglleft$f$\T1\guilsinglright}$}\raise
1.46208pt\hbox{\footnotesize$\urcorner$},\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$x$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\raise
0.69225pt\hbox{\footnotesize$\ulcorner$}\hbox{$\text{\small\char 51}$}\raise
0.69225pt\hbox{\footnotesize$\urcorner$}\bigr{\rvert}$ for
$p\in\mathbb{N}:\lvert p\rvert\leq g(m)$ if $p$
$\operatorname{\mathit{pv}}$-proves $\neg\,f(\raise
1.46208pt\hbox{\footnotesize$\ulcorner$}\hbox{$\text{\T1\guilsinglleft$f$\T1\guilsinglright}$}\raise
1.46208pt\hbox{\footnotesize$\urcorner$},\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$x$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\raise
0.69225pt\hbox{\footnotesize$\ulcorner$}\hbox{$\text{\small\char 51}$}\raise
0.69225pt\hbox{\footnotesize$\urcorner$}$ return ✓ if $p$
$\operatorname{\mathit{pv}}$-proves $\phantom{\neg\,}f(\raise
1.46208pt\hbox{\footnotesize$\ulcorner$}\hbox{$\text{\T1\guilsinglleft$f$\T1\guilsinglright}$}\raise
1.46208pt\hbox{\footnotesize$\urcorner$},\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$x$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\raise
0.69225pt\hbox{\footnotesize$\ulcorner$}\hbox{$\text{\small\char 51}$}\raise
0.69225pt\hbox{\footnotesize$\urcorner$}$ return ✗ return ✗ This gives rise to
the proof problem $\textbf{{B}}\enskip\coloneqq\enskip\Bigl{\\{}\neg\,Y(\raise
1.46208pt\hbox{\footnotesize$\ulcorner$}\hbox{$\text{\T1\guilsinglleft$Y$\T1\guilsinglright}$}\raise
1.46208pt\hbox{\footnotesize$\urcorner$},\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$x$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\raise
0.69225pt\hbox{\footnotesize$\ulcorner$}\hbox{$\text{\small\char 51}$}\raise
0.69225pt\hbox{\footnotesize$\urcorner$}:x\in\mathbb{N}\Bigr{\\}}$ where the
$x$ just serves to turn our single undecidable sentence from the First
Incompleteness Theorem into an infinitely large family. To start with, B can’t
have complexity $g(m)$ by reasoning exactly as in the First Incompleteness
Theorem. Any proofs this short would be found by $Y$ and swiftly lead to an
inconsistency. That means however that _eventually_ $Y$ will arrive at its
last line and output ✗. By the Runtime Lemma that makes every sentence in B
provable and gives this problem a well defined proof complexity $h$ such that
$g(m)<h(m)$. Now let $g^{-1}(j)\coloneqq\operatorname*{arg\,min}_{i}j\leq
g(i)$ and define
$\textbf{{A}}\enskip\coloneqq\enskip\Bigl{\\{}\mathsf{B}\land\neg\,\mathsf{S}\dots\mathsf{S}\mathsf{0}\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\mathsf{0}:\mathsf{B}\in\textbf{{B}},\;n=g^{-1}(h(m))\Bigr{\\}}$
where $m$ is the length of the original sentence $\mathsf{B}\in\textbf{{B}}$
and $n$ is the length of the new sentence $\mathsf{A}\in\textbf{{A}}$. In
other words, we add $\mathsf{S}$’s until we hit $n=g^{-1}(h(m))$. The aim here
is to dilute the sentences in B with padding until we reach complexity
$O(g(n))$ and no lower. Also, recall that by virtue of T representing at least
_some_ function we have the axioms of successorship available. First, let’s
verify that A has proof complexity $O(g(n))$. To that end note that $i\leq
g(g^{-1}(i))$ for any $i\in\mathbb{N}$. So we can say that it takes $h(m)\leq
g(g^{-1}(h(m)))=g(n)$ symbols to prove any $\mathsf{B}\in\textbf{{B}}$. And it
takes $O(n)$ symbols to prove the padding. Since $g$ is monotone and
$g(n)=\Omega(n)$ this gives us proofs of length $O(g(n))$ overall. And second,
we need to rule out $O(f(n))$ as a proof complexity for A. This time, opposite
to the previous case, note that $g(g^{-1}(i)-1)<i$ for any $i\in\mathbb{N}$.
Hence, $O(f(n))=o(g(n-1))=o(g(g^{-1}(h(m))-1))=o(h(m)).$ This leads us to
conclude that if A and thereby B had proofs of length $O(f(n))$, then it would
have proofs of length $o(h(m))$. But that contradicts the assumption that
$h(m)$ denoted the length of the shortest proofs for B. ∎

As already mentioned, this theorem is quite similar to the Nondeterministic
Time Hierarchy. Here’s a reminder.

[Nondeterministic Time Hierarchy Theorem] Let $g$ be time constructible (one
can compute $g(n)$ in $O(g(n))$ steps). And let $f(n+1)=o(g(n))$. Then
$\mathsf{NTIME}(f(n))\subsetneq\mathsf{NTIME}(g(n))$.

###### Proof.

Omitted. ∎

There are a couple of lines along which I’d like to compare these theorems.

  * •

Fascinatingly, even though the proof of the Proof Hierarchy Theorem runs along
very different lines than the proof of the Nondeterministic Time Hierarchy
Theorem (at least the one I’m familiar with by Žák Stanislav [11]), the exact
same $+1$ shows up in the formulation.

  * •

Whereas the Nondeterministic Time Hierarchy Theorem requires $g$ to be time
constructible, the Proof Hierarchy Theorem is already satisfied with mere
computability.

  * •

The Proof Hierarchy Theorem doesn’t make any demands regarding the
computational complexity of checking proofs. In particular checking whether a
given sentence is an axiom might be arbitrarily difficult.

  * •

Fundamentally, the nondeterminstic time hierarchy is a hierarchy _among_
verifiers, while the proof hierarchy is a hierarchy among the inputs for a
_fixed_ verifier.

  * •

The main sticky point that prevents this proof from being turned into a
simpler demonstration of the Nondeterministic Time Hierarchy Theorem is that,
as Manuel Blum demonstrated [2], computational problems might not have a
_minimum_ complexity $h$.

##  8 Oracles

Last but not least, I want to take a brief look at a topic near and dear to
computer scientists of the present: oracles. This section was born from the
following line of thought.

Take Peano Arithmetic and add a set of uncomputable axioms to it. The simplest
example would be an axiomatization of the Halting Problem by means of a
predicate $\mathsf{H}$ so that either $\mathsf{H}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$x$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})$ or $\neg\mathsf{H}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$x$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})$ is an axiom depending on whether or not
$x$ describes a halting computation. This theory, call it
$\textbf{{PA}}+\mathsf{H}$, is no longer computably enumerable so the First
Incompleteness Theorem won’t apply anymore. The step that breaks in the proof
is that the algorithm $X$ we constructed now needs to run an uncomputable
proof verifier as a subroutine. Sure, we could make $X$ “computable” if we
gave it access to an oracle for the Halting Problem. But then we still face
the challenge of representing $X$. We only know how to represent computable
functions, not functions that call on a Halting Problem oracle. To represent
such functions it seems that we’d also need a way to represent the Halting
Problem. Except, that is precisely what $\textbf{{PA}}+\mathsf{H}$
conveniently provides for us!

To turn this idea into a proof we’re going to need a lemma. This one:

[Oracle Representation Lemma] Let T be an extension of Q which also represents
$g$. Then T represents any function $f$ that can be computed using oracle
access to $g$.

###### Proof.

First we set up the following pair of algorithms.
$\operatorname{\mathit{check}}(\bm{x},p)$ simulate $f$ running on $\bm{x}$
using $p$ in place of the oracle for $g$ if $\ p$ isn’t of the form
$[[q_{1},r_{1}],\dots,[q_{k},r_{k}]]$ or some $q_{i}$ isn’t a query $f$
would’ve made  return ✗ return ✓ $\operatorname{\mathit{sim}}(\bm{x},p)$
simulate $f$ running on $\bm{x}$ using $p$ in place of the oracle for $g$
return $f(\bm{x})$ It’ll also be helpful to use a shorthand:
$\mathsf{\operatorname{\textsf{{Ck}}}(\bm{\mathsf{a}},p)\enskip\coloneqq\enskip\operatorname{\mathit{check}}(\bm{\mathsf{a}},p)\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\raise
0.69225pt\hbox{\footnotesize$\ulcorner$}\hbox{$\text{\small\char 51}$}\raise
0.69225pt\hbox{\footnotesize$\urcorner$}\land\forall i<\lVert
p\rVert\colon\text{$g$}(p[i][0])\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}p[i][1]}.$
Then we can give this representation of $f$:
$\text{\T1\guilsinglleft$f$\T1\guilsinglright}(\bm{\mathsf{a}},\mathsf{b})\enskip\coloneqq\enskip\mathsf{\exists
p\colon\operatorname{\mathit{sim}}(\bm{\mathsf{a}},p)\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}b\land\operatorname{\textsf{{Ck}}}(\bm{\mathsf{a}},p)\land\forall
q<p\colon\neg\operatorname{\textsf{{Ck}}}(\bm{\mathsf{a}},q)}.$ To verify that
‹$f$› works as advertised suppose that $f(\bm{x})=y$. Then there exists a
shortest corresponding list $p_{y}$ of oracle inputs and outputs for $g$. This
allows us to confirm the following two required statements.
_$\textbf{{T}}\vdash\text{$f$}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$y$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})$_ : This easily follows from
$\mathsf{\operatorname{\mathit{sim}}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$p_{y}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$y$}\raise
0.0pt\hbox{\footnotesize$\urcorner$}\land\operatorname{\textsf{{Ck}}}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$p_{y}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})\land\forall q<\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$p_{y}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$}\colon\neg\operatorname{\textsf{{Ck}}}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},q)}.$ The bounded quantifiers can be
converted into a finite number of cases using the Bounded Quantification Lemma
(or the non-existence of negative objects if necessary). And then it’s just a
matter of proving via the Runtime Lemma that a bunch of representable
functions when run on the given inputs produce the given outputs.
_$\textbf{{T}}\vdash\mathsf{\exists!y\colon\text{$f$}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},y)}$_ : Existence has already been
handled. For uniqueness let $\mathtt{y}$ be arbitrary and assume
$\text{\T1\guilsinglleft$f$\T1\guilsinglright}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},\mathtt{y})$. Then we can obtain a
$\mathtt{p}$ such that $\mathsf{\text{$\operatorname{\mathit{sim}}$}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},\mathtt{p})\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\mathtt{y}\land\operatorname{\textsf{{Ck}}}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},\mathtt{p})\land\forall
q<\mathtt{p}\colon\neg\operatorname{\textsf{{Ck}}}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},q)}.$ Of course, we also still have
$\mathsf{\text{$\operatorname{\mathit{sim}}$}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$p_{y}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$y$}\raise
0.0pt\hbox{\footnotesize$\urcorner$}\land\operatorname{\textsf{{Ck}}}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$p_{y}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})\land\forall q<\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$p_{y}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$}\colon\neg\operatorname{\textsf{{Ck}}}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},q)}$ which allows us to easily rule out
$\mathtt{p}<\raise 0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$p_{y}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$}$ and $\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$p_{y}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$}<\mathtt{p}$. By the totality of ordering
w.r.t. numerals this leaves
$\mathtt{p}\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$p_{y}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$}$ and hence
$\mathsf{\mathtt{y}\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\operatorname{\mathit{sim}}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},\mathtt{p})\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\operatorname{\mathit{sim}}(\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$\bm{x}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$},\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$p_{y}$}\raise
0.0pt\hbox{\footnotesize$\urcorner$})\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\raise
0.0pt\hbox{\footnotesize$\ulcorner$}\hbox{$y$}\raise
0.0pt\hbox{\footnotesize$\urcorner$}}.\qed$

We’re now ready to answer the question regarding the completeness of
$\textbf{{PA}}+\mathsf{H}$ we raised earlier. After all, this theory can
represent the Halting Problem through the formula
$(\mathsf{H}(\mathsf{x})\land\mathsf{y}\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\raise
0.69225pt\hbox{\footnotesize$\ulcorner$}\hbox{$\text{\small\char 51}$}\raise
0.69225pt\hbox{\footnotesize$\urcorner$})\lor(\neg\mathsf{H}(\mathsf{x})\land\mathsf{y}\mathrel{\leavevmode\resizebox{1.15pt}{1.15pt}{=}}\raise
0.76207pt\hbox{\footnotesize$\ulcorner$}\hbox{$\text{\small\char 55}$}\raise
0.76207pt\hbox{\footnotesize$\urcorner$})$. So using the Oracle Representation
Lemma it can represent the algorithm $X$ used in the proof of the First
Incompleteness Theorem and the argument goes through.

Unquestionably, there’s a general theorem to be proved here. Actually, many of
them are. For the assumption of computable enumerability can now be relaxed
throughout our entire discussion. But relaxed to what? What exactly
characterizes theories such as $\textbf{{PA}}+\mathsf{H}$? One option might be
“a theory which represents its own proof verifier”. But that assumption is
already doing half the work towards establishing our generalization. There has
to be something more natural.

In any case, I need to rest the pen eventually. And this seems like a good
place to take a break.

> The final answer I badly wanted to blurt out was: “No, no, no—Gödel _did_
> invent Lisp!”
>
> Douglas Hofstadter

## References

  * [1] Paul Bernays, David Hilbert Grundlagen der Mathematik, Volume II, Springer-Verlag, 1939
  * [2] Manuel Blum, A Machine-Independent Theory of the Complexity of Recursive Functions, Journal of the ACM, Apr. 1967
  * [3] Stephen Arthur Cook, A hierarchy for nondeterministic time complexity, Proceedings of the fourth annual ACM symposium on Theory of computing, STOC ’72, pp. 187–192
  * [4] Kurt Gödel, Über formal unentscheidbare Sätze der Principia Mathematica und verwandter Systeme I, Monatshefte für Mathematik und Physik, Issue 38, Dec. 1931, pp. 173–198
  * [5] Kurt Gödel, Über die Länge von Beweisen, Ergebnisse eines mathematischen Kolloquiums, Issue 7, 1936, pp. 23–24
  * [6] Juris Hartmanis, Richard Edwin Stearns, On the computational complexity of algorithms, Transactions of the American Mathematical Society, Issue 117, May 1. 1965, pp. 285–306
  * [7] Robert G. Jeroslow, Redundancies in the Hilbert-Bernays Derivability Conditions for Gödel’s Second Incompleteness Theorem, Journal of Symbolic Logic, Volume 38, No. 3, Sep. 1973, pp. 359–367
  * [8] Sebastian Oberhoff, Incompleteness Ex Machina, Bulletin of the European Association for Theoretical Computer Science, No 128, June 2019, pp. 91–105
  * [9] John Barkley Rosser, Extensions of Some Theorems of Gödel and Church, Journal of Symbolic Logic, Volume 1, No. 3, Sep. 1936, pp. 87–91
  * [10] Robert Irving Soare, Computability and Recursion, Bulletin of Symbolic Logic 2, 1996, pp. 284–321
  * [11] Žák Stanislav, A Turing machine time hierarchy, Theoretical Computer Science, Elsevier Science B.V. 26, Oct. 1983, pp. 327–333
  * [12] Alan Turing, On Computable Numbers, with an Application to the Entscheidungsproblem, Proceedings of the London Mathematical Society, Volume s2–42, Issue 1, 1937, pp. 230–265