# Comparing Different Approaches for Stellar Intensity Interferometry

Sebastian Karl,1 Andreas Zmija,2 Stefan Richter,1,3 Naomi Vogel,2 Dmitry
Malyshev,2 Adrian Zink,2 Thilo Michel,2 Gisela Anton,2 Joachim von
Zanthier,1,3 and Stefan Funk2
1Institute for Optics, Information and Photonics, Friedrich-Alexander-
Universität Erlangen-Nürnberg, Staudtstr. 1, Erlangen D 91058, Germany
2Erlangen Centre for Astroparticle Physics, Friedrich-Alexander-Universität
Erlangen-Nürnberg, Erwin-Rommel-Str. 1, Erlangen D 91058, Germany
3Erlangen Graduate School in Advanced Optical Technologies, Paul-Gordan-Str.
6, Erlangen D 91052, Germany
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Stellar intensity interferometers correlate photons within their coherence
time and could overcome the baseline limitations of existing amplitude
interferometers. Intensity interferometers do not rely on phase coherence of
the optical elements and thus function without high grade optics and light
combining delay lines. However, the coherence time of starlight observed with
realistic optical filter bandwidths ($>$0.1\text{\,}\mathrm{nm}$$) is usually
much smaller than the time resolution of the detection system
($>$10\text{\,}\mathrm{ps}$$), resulting in a greatly reduced correlation
signal. Reaching high signal to noise in a reasonably short measurement time
can be achieved in different ways: either by increasing the time resolution,
which increases the correlation signal height, or by increasing the photon
rate, which decreases statistical uncertainties of the measurement. We present
laboratory measurements employing both approaches and directly compare them in
terms of signal to noise ratio. A high time-resolution interferometry setup
designed for small to intermediate size optical telescopes and thus lower
photon rates (diameters $<\,$some meters) is compared to a setup capable of
measuring high photon rates, which is planned to be installed at Cherenkov
telescopes with dish diameters of $>$10\text{\,}\mathrm{m}$$. We use a Xenon
lamp as a common light source simulating starlight. Both setups measure the
expected correlation signal and work at the expected shot-noise limit of
statistical uncertainties for measurement times between $10\text{\,}\min$ and
$23\text{\,}\mathrm{h}$. We discuss the quantitative differences in the
measurement results and give an overview of suitable operation regimes for
each of the interferometer concepts.

###### keywords:

instrumentation: high angular resolution – instrumentation: interferometers –
telescopes – stars: imaging

††pubyear: 2021††pagerange: Comparing Different Approaches for Stellar
Intensity Interferometry–Comparing Different Approaches for Stellar Intensity
Interferometry

##  1 Introduction

In recent years, interest in intensity interferometry for astronomical
observations has grown. Prospects of precise measurements of angular diameters
of stars and resolving structures on stellar surfaces (Nuñez et al., 2012;
Nuñez & Domiciano de Souza, 2015; Dravins et al., 2013) resulted in a revival
of Hanbury-Brown Twiss interferometry in astrophysics.

High optical resolutions of less than $1\text{\,}\mathrm{\text{mas}}$ require
interferometers with baselines of several hundred metres when measuring at
optical wavelengths. Despite the problems due to atmospheric turbulences and
the challenge of combining the light paths from the telescopes with sub-
wavelength precision, stellar interferometers based on amplitude
interferometry provide successful measurements at such baselines, see for
example Roettenbacher et al. (2016); Che et al. (2011); Baines et al. (2007);
Mourard, D. et al. (2015). The Center for High Angular Resolution Astronomy
(CHARA) - one of the largest amplitude interferometers - has baselines of
around $300\text{\,}\mathrm{m}$ providing angular resolution capabilities up
to $0.4\text{\,}\mathrm{m}\mathrm{a}\mathrm{s}$ (ten Brummelaar et al., 2005).

The concept of intensity interferometry is promising for even longer
baselines. Correlating the light intensities rather than amplitudes provides a
robust way to overcome atmospheric disturbances at the cost of the achievable
limiting magnitudes. Pioneering research at the Narrabri Stellar Intensity
Interferometer was carried out by Hanbury Brown and Twiss (HBT) (Brown &
Twiss, 1956; Brown et al., 1967; Hanbury Brown et al., 1974) in the second
half of the 20th century. The challenge in this technique is the low amount of
thermal photons that can actually interfere compared to the total number of
arriving photons, increasing the difficulty in retrieving the signals in HBT
interferometry. For this kind of photon correlation measurement, high
electronic stability, photo-detectors with high time resolution and large
light collection areas (telescopes) are desirable. These are preconditions
that are well achieved by today’s telescopes and photo-detection equipment.

The performance of intensity interferometers is determined by a wide range of
parameters, such as the aperture of the telescopes, the optical bandwidth of
the correlated light, the timing resolution of the detection system, the
maximum possible photon rate, and the projected baseline, all of which
influence the expected signal-to-noise-ratio (S/N) of the measurement.
However, when scanning the available parameter space for stellar intensity
interferometry measurements at existing or future observatories, it appears
that there are two promising methods for enabling successful measurements of
the correlation signal.

The first method attempts to reach high timing resolution
($<$1\text{\,}\mathrm{ns}$$) and narrow optical bandwidths
($<$1\text{\,}\mathrm{nm}$$), keeping the correlation signal as pronounced as
possible. This approach requires optical telescopes with photon timing
dispersion smaller than the timing resolution of the detection system. Such
telescopes typically enable precise control of the light path necessary for
sharp optical filtering. As a downside, they usually do not come with large
photon collection areas resulting in smaller photon count rates at the
detectors compared to possible intensity interferometers at the world’s
largest telescopes. Also, high timing resolution photon counting detectors
lose efficiency when operated at photon rates higher than some (tens of) MHz
due to the detector dead time, reducing the detected photon rate further. This
may lead to rather long measurement times of up to a few hours for measuring
significant intensity correlations even with relatively bright stars of
apparent magnitude 3. An option to overcome this downside is using arrays of
telescopes, where each subset measures correlations along the same baseline,
or spectral and polarization multiplexing (Trippe et al., 2014). In absence of
arrays, the main challenge is to obtain stable measurement conditions over
periods of up to several hours. In this context, some promising results have
been achieved in recent years. Temporal photon bunching peaks of $\alpha$ Boo,
$\alpha$ CMi and $\beta$ Gem, providing magnitudes between $-1.68$ and $-0.11$
in the near-infrared, were measured at the C2PU observatory in France with a
telescope diameter of not more than $1\,$m at integration times between $2$
and $7$ hours (Guerin et al., 2017). Using two of these telescopes enabled
spatial bunching measurements of $\alpha$ Lyr, $\beta$ Ori and $\alpha$ Aur
(Guerin et al., 2018). The same telescopes were recently used to estimate the
distance to P Cygni using intensity interferometry (Rivet et al., 2020).
Further, the AQUEYE+ and IQUEYE astronomical high-timing resolution
instruments can be used for stellar intensity interferometry at sub-nanosecond
time resolution (Zampieri et al., 2016). Finally intensity interferometry
observations and diameter estimations on 5 stars were performed by Horch et
al. (2021) at the Southern Connecticut Stellar Interferometer.

In the second method, very large telescopes or telescope arrays are used. In
this context, Imaging Air Cherenkov telescopes (IACTs) are promising
candidates for photon correlation measurements (Bohec & Holder, 2006), since
they provide light collection areas of about $100\text{\,}\mathrm{m}^{2}$
(e.g. VERITAS (Maier, 2007) or H.E.S.S. (Bernlöhr et al., 2003)). Depending on
the optical design Cherenkov telescopes have different levels of internal
unisochronicity that make it impossible to achieve timing resolutions below
these time differences. While Davies-Cotton telescopes, such as VERITAS
(Weekes et al., 2002) and H.E.S.S (Bernlöhr et al., 2003) incur time
differences in the light path of a few nanoseconds, parabolic IACTs such as
MAGIC (Shayduk et al., 2005) or the Cherenkov Telescope Array (CTA) Large-
Sized Telescope (LST) reduce it to $\approx 100-400\,$picoseconds, which is
still significant compared to imaging-grade telescopes. In addition, imperfect
optics make it hard to collimate the light for optical filters, requiring
bandwidths of $1\text{\,}\mathrm{nm}$ or more if interferometric filters are
used. Photo-detection for such telescopes needs to be designed to work at
several hundreds of MHz. Even though the correlation signal is very small in
this arrangement the high photon rates suppress the shot noise enabling high
S/N in relatively short measurements times of a few minutes. The main
challenge here is the electronic stability on a high statistical level. Using
the MAGIC telescopes Acciari et al. (2019) measured photon bunching signals of
a few stars within minutes (for single baselines). Abeysekara et al. (2020)
used the four VERITAS telescopes to measure the angular diameter of two stars
by taking measurements at many different projected telescope baselines with
averaged measurement times of $17\text{\,}\min$ per baseline, demonstrating
the potential of IACTs for stellar intensity interferometry.

In this paper we shortly recapitulate the theory of intensity correlations in
section 2. In section 3 we describe the setup used to compare our
complementary approaches in the laboratory. We then evaluate our measurement
results in section 4. Finally, we compare and discuss our measurement results
in section 5, illustrating the scope of both methods. Here we pay special
attention to the S/N of our measurements, since it is the central parameter to
quantify the quality of photon correlation measurements. High timing
resolutions as well as high photon rates favour a high S/N, however, both are
difficult to achieve simultaneously with current state-of-the-art technology.
We interpret our measurements in order to quantify telescope layouts for both
approaches which are promising for high S/N. Further we will discuss future
observation possibilities based on upcoming detector developments, which will
especially boost the HTR technology towards higher photon rate capabilities
and thus larger telescopes.

##  2 Theory of intensity correlations

The observable of intensity interferometry is the normalized second order
correlation function

$g^{(2)}(\textbf{{r}},\tau)=\frac{\expectationvalue{I(\textbf{{R}},t)I(\textbf{{R}}+\textbf{{r}},t+\tau)}}{\expectationvalue{I(\textbf{{R}},t)}\expectationvalue{I(\textbf{{R}}+\textbf{{r}},t+\tau)}},$
(1)

where $I(\textbf{{R}},t)$ denotes the intensity measured at position R at a
time $t$ and the brackets denote a time average. For thermal light sources
(TLS) the correlations in space and time can be factorized into the respective
absolute value squared of the first order correlation functions (Mandel &
Wolf, 1995):

$g^{(2)}(\textbf{{r}},\tau)=1+\absolutevalue{g^{(1)}(\textbf{{r}})}^{2}\absolutevalue{g^{(1)}(\tau)}^{2}.$
(2)

Here $\absolutevalue{g^{(1)}(\textbf{{r}})}^{2}$ corresponds to the modulus
squared of the measured visibility of amplitude interferometry (Labeyrie et
al., 2006) and is related to the source geometry by the Fourier relationship
of the van Cittert Zernike theorem (Mandel & Wolf, 1995). In the following,
the correlation functions’ spatial part is assumed to be $1$, as we restrict
ourselves to correlations of two detectors having zero separation. In analogy
to the spatial part, the temporal first order correlation can be calculated
using a Fourier transform relationship - the Wiener Khintchine theorem (Mandel
& Wolf, 1995):

$g^{(1)}(\tau)=\int_{-\infty}^{\infty}s(\omega)\mathrm{e}^{-i\omega\tau}\differential{\omega}$
(3)

where $s(\omega)$ denotes the normalized source spectrum.

In general, TLS show bunching, meaning it is more likely to detect two photons
arriving coincidentally rather than individually for a sufficiently small
coincidence window. This is signified by the second order correlation function
having its’ maximum $g^{(2)}(\tau=0)=2$ at zero time delay. For a time delay
much larger than the source coherence time $\tau_{c}$, the second order
correlation function drops off to its baseline value of 1. The corresponding
time delay is proportional to the source coherence time $\tau_{c}$, which can
be defined as the integral over the correlation peak if the light is
considered to be fully polarized

$\tau_{c}=\int_{-\infty}^{\infty}\absolutevalue{g^{(1)}(\tau)}^{2}\differential{\tau}\approx\frac{\lambda^{2}}{c\Delta\lambda}.$
(4)

According to eq. 4 the coherence time is approximated in terms of the
spectrum’s central wavelength $\lambda$, the spectral bandwidth
$\Delta\lambda$ and the speed of light $c$ eq. 4 (Fox, 2006). For a very
narrow rectangular filter, the approximate equality in eq. 4 holds rigorously.

For a real-world correlation measurement, the theoretical correlation function
from eq. 2 has to be convolved with the detector response to yield the
measurement expectation. Thus, for an electronic time resolution
$\tau_{\textnormal{e}}$, the peak amplitude $A$ is reduced to a value
$A=\frac{\tau_{\textnormal{c}}}{\tau_{\textnormal{e}}}$ from its initial value
of 1 (Loudon, 2000). A further decrease by a factor of $\frac{1}{2}$ is
incurred if no polarization filtering is employed, as will be the case in our
measurements.

Assuming only statistical fluctuations the S/N of the HBT bunching peak of
randomly polarized thermal light is given by (Brown et al., 1957)

$S/N=A\,n(\omega)\,\eta\,\sqrt{\frac{T_{\mathrm{meas}}}{\tau_{\textnormal{e}}}},$
(5)

where $A$ denotes the telescope collection area, $n(\omega)$ denotes the
source spectral flux density, $\eta$ denotes the detector quantum efficiency,
and $T_{\mathrm{meas}}$ denotes the measurement time. Note that the spatial
correlation factor is 1 in our measurements and thus is omitted in further
calculations. Equation 5 shows that for a HBT experiment, the filter bandwidth
does not influence the S/N measurement. It should however be noted that eq. 5
is only valid if the measured photon rate is approximately constant over the
whole measurement time. This is not the case for both the single photon and
current correlation measurements that are presented in this paper. Instead we
estimate the theoretical S/N utilizing shot noise as a lower limit to the
noise that might be incurred in the noisy correlation measurement.

##  3 Methods

Our two systems will be denoted as

  * •

HTR : "High time resolution" - photon counting measurements aiming for very
high time resolutions, but restricted in photon rate capability.

  * •

HPR : "High photon rates" - Photomultiplier tube (PMT) current measurements
aiming for very high photon rates, but restricted in time resolution.

Figure 1 shows a sketch of the combined setup. To simulate light coming from a
distant star (black body spectrum, spatially coherent), a high pressure xenon
arc lamp (XBO) is used (Tan et al., 2014). Light from this black body spectrum
source is coupled into a single mode fiber placed at an appropriate distance
to the XBO, ensuring the collection of only one spatial mode. After the light
is transported through the single mode fiber, it is collimated using an
achromatic fiber port. The collimated beam can now be transferred to any of
the setups by installing the flippable
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$ mirror into the beam for the HTR
setup, and removing it from the beam for the HPR setup. Details of both setups
are described in the following subsections.

Figure 1: Measurement setup including both the HTR and the HPR parts, which
can be swapped in operation by either installing or removing the flippable
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$ mirror. The interference filters are
denoted as follows: central wavelength/FWHM bandwidth. The high pressure xenon
arc lamp is denoted by XBO, and the Photomultiplyer tubes and the hybrid photo
detectors are denoted as PMT and HPM respectively. ND 1.5 refers to a neutral
density filter of optical density 1.5 placed in the HTR setup.

###  3.1 The HTR setup

The HTR setup is designed for use on optical telescopes, imprinting an optical
time spread of not more than $100\text{\,}\mathrm{ps}$ onto the incident
photons, and yielding light at a low intensity of less than
$10\text{\,}\mathrm{MHz}$ when filtered down to bandwidths
$\leq$1\text{\,}\mathrm{nm}$$. In this case, the second order correlation
function can be measured using single photon counting detectors and calculated
from the measured arrival time differences.

We use two HPM 100-06 Hybrid Photo detectors (HPDs) manufactured by
Becker&Hickl to measure the temporal correlation function using the HBT setup.
The photon arrival signals are discretized and correlated on a SPC-130-EMN
time-to-amplitude-converter (TAC) based card manufactured by Becker& Hickl. As
these detectors have a maximum continuous count rate of
$15\text{\,}\mathrm{MHz}$, usage at higher count rates within this regime
results in undetected photons due to their $\approx$70\text{\,}\mathrm{ns}$$
dead time. We measured the full width at half maximum (FWHM) overall system
timing resolution to be $41.61(8)\text{\,}\mathrm{ps}$ and the detectors
quantum detection efficiency to be 21 percent in the relevant wavelength
range, in close agreement to the manufacturer specifications (Becker & Hickl
GmbH, 2017). Considering that the detector dead time is orders of magnitude
larger than the final histogram length, simple start stop histogram
correlation is sufficient to calculate the correlation function accurately and
thus TAC based correlation can be utilized. To shift the signal out of the TAC
switching systematics and eliminate the loss of correlation events due to the
asymmetry of the TAC input channels, a cable delay of
$35\text{\,}\mathrm{n}\mathrm{s}$ is introduced.

Since the HPM detectors have a maximal photon rate of
$15\text{\,}\mathrm{MHz}$, the light directed into the HTR setup is attenuated
by a neutral density (ND) filter of optical density 1.5, resulting in a 97.2
percent decrease in incident photon flux (Thorlabs, Inc, 2021). Afterwards the
attenuated light is spectrally filtered using a narrow interference filter of
central wavelength $416\text{\,}\mathrm{n}\mathrm{m}$ FWHM width
$0.5\text{\,}\mathrm{nm}$ manufactured by Alluxa and directed onto the HPM
cathodes using a non-polarizing 50:50 beam splitter.

In addition to pronounced switching systematics, the TAC imprints further
systematics on the correlation measurements, consisting of two superimposed
oscillations and an increasing slope. Since these systematics cannot be
neglected, the measurements have to be calibrated. The calibration measurement
is performed using light from a halogen light bulb, which is not spectrally
filtered and incoherently scattered into a multimode fibre connected to our
detectors. The measurement produces a nearly flat correlation histogram, with
the HBT peak suppressed by a factor
$\approx$5\text{\times}{10}^{-5}\text{\,}$$ compared to the visibility
measured with spectral and spatial filtering. Subtracting the calibration
measurement from the actual measurement yields a well resolved correlation
peak. This correlation peak is then shifted by
$\approx$35\text{\,}\mathrm{ns}$$ to account for the cable delay.

###  3.2 The HPR setup

The HPR setup is designed to be operated at Imaging Air Cherenkov Telescopes
(IACTs) and therefore deals with different parameters. Filter bandwidths of
not less then $2\text{\,}\mathrm{nm}$ and thus photon rates of more than
$100\text{\,}\mathrm{MHz}$ even for stars of intermediate magnitude are
realistic. Typical telescope timing precisions of a few ns due to the
unisochronicity of different light paths in the telescopes decrease the
measured correlation signals strongly, however, the high photon fluxes of
several hundreds of MHz at the detectors decrease shot noise maintaining a
reasonable S/N.

At photon rates of hundreds of MHz, time stamping of single photons becomes
difficult. Instead, the photon current of the photo-detector is measured in
both observation channels. The interferometer consists of two Photomultiplier
Tubes of type Hamamatsu R11265U, whose output current is amplified by a factor
of ten and then digitized with $1.6\text{\,}\mathrm{ns}$ sampling time using
the M4i.2212-x8 digitizer card from Spectrum Instrumentation.

The incoming beam from the single mode fiber is expanded to a diameter of
$\approx$1\text{\,}\mathrm{cm}$$ using a Keplerian beam expander, in order to
not damage the photo-cathodes due to the high light density of the former
beam.

If $\textbf{{A}}(t)$ and $\textbf{{B}}(t)$ are the vectors of the time-sampled
photo-currents with subtracted mean from the two PMT channels, then the cross
correlation with time shift $\tau$ is determined as:

$G^{(2)}(\tau)=\textbf{{A}}(t)\cdot\textbf{{B}}(t+\tau)$ (6)

Normalization is done by dividing every value by the mean value of the
correlation baseline, where the baseline is defined as
$G^{(2)}(\tau)|_{\tau\gg\tau_{C}}$.

$g^{(2)}(\tau)=\frac{G^{(2)}(\tau)}{\overline{G^{(2)}}|_{\tau\gg\tau_{c}}}$
(7)

This correlation procedure and the fact, that single photon pulse shapes are
extended over $>30$ sampling time bins induces correlation effects in the
$g^{(2)}$ of neighbouring bins. The baseline in part a) of fig. 3 shows that
neighbouring $g^{(2)}$ values are not independent of each other. This
correlation also results in a decreased value of the root mean square error
(RMSE) expectation level compared to Poissonian shot noise calculations. A
correction factor is obtained by averaging over experimental RMSE observations
at short measurement times and confirmed by waveform simulations. In
particular, the entire $10\text{\,}\min$ of measurement time are divided in
$175$ chunks each of $\approx$3.4\text{\,}\mathrm{s}$\,$ length. At these
short time scales the RMSE is assumed to be purely statistical for
$\tau\gg\tau_{\textnormal{c}}$. We further cross check the calculation of the
correction factor by waveform simulations which then are correlated. These
simulations consist of an "empty" waveform baseline to which photons are added
via a photon pulse shape template created by averaging many single photon
pulse shapes measured at low rates. The simulated photon heights also follow
measured photon pulse height distributions at low rates where the heights of
single photons can be determined. No photon bunching is simulated in the
waveforms, but the RMSE of the resulting correlation can be compared to the
expectation of the single photon time tagging shot noise expectation of the
simulated rate and directly yields the RMSE correction factor of $0.31$ in
this measurement. The RMSE correction factor is then multiplied to the
Poissonian shot noise expectation to yield the current correlation RMSE
expectation.

A detailed description of the general HPR setup and its measurement analysis
chain can be found in Zmija et al. (2021)

##  4 Experiments and Results

The correlation results are analysed in two aspects. First, the measurement of
the correlation time is analysed to investigate the coherence time. Second,
the fluctuation of the correlation function baseline is analysed in the region
around the peak to investigate possible background systematics in the
measurement. A good index is the RMSE value of the $g^{(2)}$ baseline.
Assuming no systematics, the RMSE value should be given by the Poissonian
noise statistics.

Figure 2: Measurement results for the HTR setup. Part a) shows the measured
second order correlation function after an accumulation time of
$22.7\text{\,}\mathrm{h}$, as well as a measurement expectation rescaled to
the measured coherence time and a Gaussian fit. Part b) shows the HTR baseline
RMSE time evolution over the accumulation time. Figure 3: Measurement results
for the HPR setup. Part a) shows the measured second order correlation
function after an accumulation time of $10\text{\,}\mathrm{min}$, as well as
the measurement expectation calculated numerically and a fit of that numerical
model to the measured correlation function. Part b) shows the HPR baseline
RMSE time evolution over the accumulation time.

Part a) of fig. 2 shows the correlation function measured by the HTR setup
($0.5\text{\,}\mathrm{nm}$ optical bandwidth, centred at
$416\text{\,}\mathrm{nm}$) after a measurement time of
$22.7\text{\,}\mathrm{h}$, where the error bars correspond to the baseline
RMSE. Evaluating the coherence time by numerical integration we recover a
value of about 85 percent of the expected coherence time ($\tau_{\text{c,
exp}}=$0.425\text{\,}\mathrm{ps}$$, $\tau_{\text{c,
meas}}=$0.35(3)\text{\,}\mathrm{ps}$$). The expected coherence time is
calculated from a manufacturer-supplied filter transmission spectrum. A
Gaussian fit to the measurement yields a coherence time of $\tau_{\text{c,
fit}}=$0.36(3)\text{\,}\mathrm{ps}$$, and is also shown in fig. 2. The
coherence times recovered by fit and numerical integration agree well with
each other, leading to the conclusion that a coherence time about 15 percent
lower than the expectation has been measured. Part a) of Figure 2 also shows
an expectation curve calculated from measurements of the HTR timing jitter and
a Gaussian filter, rescaled to yield a coherence time matching the measured
value. This curve is in good agreement with the measured correlation function,
and only nearly exceeds the $1\sigma$ boundary of the measurement values for
values centred around zero time delay. The reduction in measured coherence
time could be explained by slight errors in the manufacturer supplied filter
transmission spectrum, which could fully explain the observed deviations even
for a spectrometer FWHM resolution $<$0.1\text{\,}\mathrm{nm}$$.

Part b) of fig. 2 shows the RMSE behaviour of the HTR $g^{(2)}$ function
evaluated cumulatively over the whole measurement time. The RMSE is evaluated
directly on the measured correlation histogram baseline for each time period.
In order to estimate the error of our calculated RMSE, it is evaluated in 20
independent samples of size more than 20 times larger than the measured
correlation peak FWHM. The theoretical predictions for the expected RMSE is
directly extracted from the mean number of counts per bin of both the raw and
the normalization measurement, and are approximated very well by our
measurements.

Part a) of fig. 3 shows the correlation function obtained with the HPR setup
using the XBO lamp filtered at $465\text{\,}\mathrm{nm}$ with a FWHM of
$2\text{\,}\mathrm{nm}$. The error bars correspond to the baseline RMSE. Since
timing resolution and sampling time are far larger than the correlation time,
the expected $g^{(2)}$ signal shape is dominated by the correlation pulse
shape of PMT pulses. The expected pulse shape is derived by correlating the
average single photon pulse shapes of both channels obtained from a
calibration measurement, as previously described. The resulting height of the
correlation signal can then be scaled to fit the requirement of its integral
being the coherence time as defined (labelled "expectation" in fig. 3).
Furthermore, the shape can also be fitted to the data with the height and the
peak position being free fit parameters, and allows to derive the experimental
value of the coherence time. Both are plotted in fig. 3 and show good
agreement with each other. The theoretical coherence time stems from numerical
calculations using the optical filter function measured with a high-resolution
spectrometer.

Part b) of fig. 3 shows the RMSE of the analysed $g^{(2)}$ function evaluated
cumulatively over the whole measurement time. The uncertainty is computed by
taking $8$ independent parts of $g^{(2)}|_{\tau\gg\tau_{\textnormal{c}}}$ each
of $25$ time bins length. The dashed line indicates the Poissonian shot noise
statistics following nearly a $1/\sqrt{T}$ curve only slightly affected by the
non-constant photon rates. It is apparent that the measured RMSE values follow
the theory curve quite well for the whole measurement time of
$10\text{\,}\min$. Other than shot noise, contributions are well below the
$10^{-6}$ level.

##  5 Discussion

Figure 4: Illustration of parameter definitions for the S/N calculation. The
bunching peak FWHM is shown in green and the RMSE interval around the second
order correlation function baseline is shaded in orange.

In this section, we compare the two methods with respect to S/N in a given
measuring time.

It can be seen from the RMSE plots in fig. 2 and fig. 3 that both methods work
at the minimum noise level, which means that they are dominated by the photon
number Poisson noise (shot noise). As presented in fig. 2 and fig. 3 for both
methods we can determine the height $h$ of the $g^{(2)}$-correlation peak, the
FWHM time width $\Delta t$ of the peak and the noise RMSE of the
$g^{(2)}$-curve per time bin of bin width $t_{\text{bin}}$ The parameter
definitions are also displayed in fig. 4.

The signal-to-noise ratio of a correlation measurement can be approximated as

$\frac{S}{N}=\frac{h}{\text{RMSE}}\,\sqrt{\frac{\Delta t}{t_{\text{bin}}}},$
(8)

multiplying the peak height and RMSE error by the ratio of the correlation
peak FWHM and correlation bin size to yield an effective integral of the
correlation peak and its RMSE error.

As an example we compare the S/N values for the two setups at a measurement
time of 10 minutes. The corresponding RMSE in the HTR (HPR) setup is the first
(last) data point in fig. 2 (fig. 3). The S/N ratios are

$\left(\frac{S}{N}\right)_{\text{HTR}}=\frac{0.0056}{0.0069}\,\sqrt{\frac{60\,\text{ps}}{12\,\text{ps}}}=1.8$
(9) $\left(\frac{S}{N}\right)_{\text{HPR}}=\frac{1.27\cdot 10^{-5}}{3.51\cdot
10^{-7}}\,\sqrt{\frac{8\,\text{ns}}{1.6\,\text{ns}}}=80.9$ (10)

The difference between the two results of the different setups can be
explained by the time resolution (resp. peak width) $\Delta t$ and the
relative intensity $f_{\text{trans}}$ between the two setups. The relative
intensity is determined by the applied grey filter $f_{\text{trans,
grey}}=0.028$, the difference in quantum efficiency $f_{\text{trans, QE}}=0.5$
and the difference in source brightness at the measured wavelengths
$f_{\text{trans, spectrum}}=0.5$, all reducing the HTR photon flux relative to
$f_{\text{trans, HPR}}=1$. Further, one has to correct for the different
amount of temporal coherence intrinsically determined by the respective
filter’s central wavelength, $f_{\text{coh}}$. The achieved S/N-values should
be proportional to $f_{\text{trans}}f_{\text{coh}}/\sqrt{\Delta t}$. For the
HPR setup the effective time resolution differs from the peak width due to the
net width of photon pulse correlations. Assuming no change in S/N between the
concepts of time tagging and current correlation, $\Delta
t_{\text{HPR}}=c^{2}\cdot\Delta t_{\text{peak}}$ holds, where $c=0.31$ is the
previously mentioned RMSE correction factor.

With these values we expect

$\frac{\left(\frac{S}{N}\right)_{\text{HTR}}}{\left(\frac{S}{N}\right)_{\text{HPR}}}=\frac{(f_{\text{trans}}f_{\text{coh}}/\sqrt{\Delta
t})_{\text{HTR}}}{(f_{\text{trans}}f_{\text{coh}}/\sqrt{\Delta
t})_{\text{HPR}}}=\frac{0.0075\cdot 0.8/\sqrt{60\,\text{ps}}}{1\cdot
1/\sqrt{0.76\,\text{ns}}}=0.021$ (11)

which is in excellent agreement with the measured S/N values.

In our laboratory measurements, the HPR setup provides a much higher S/N in
the same measurement time than the HTR setup. However, this is mainly because
of the strong reduction of photon flux in the HTR setup, which is analogous to
measuring the same star with a smaller telescope. From the rates measured in
the laboratory, we extract the virtual telescope areas employed for our
autocorrelation measurement. As an exemplary star we choose a magnitude 2 star
and calculate that HPR is measuring with an equivalent of two $195$ m2
telescopes ($15.8\,$m diameter), while HTR uses two telescopes of size $1.17$
m2 ($1.2\,$m diameter). A total light collection and detection efficiency of
$0.1$ is assumed here for both systems.

Telescope | Area (m2) |  $\Delta t$ (ns) | Reference
---|---|---|---
H.E.S.S. I. | 100 | 5 | Bernlöhr et al. (2003)
CTA SST (req) | 8 | $1.5\times 2.35$ | Rulten et al. (2016)
CTA MST (req) | 88 | $0.8\sqrt{12}$ | Schlenstedt (2014)
CTA LST | 370 | 0.1 |
MAGIC | 227 | 0.4 | Shayduk et al. (2005)
VERITAS | 113 | 4 | Weekes et al. (2002)
PW CDK 20 | 0.19 | $\approx 0.01$ | Cavazzani et al. (2012)
PW CDK 1000 | 0.8 | $\approx 0.01$ | Cavazzani et al. (2012)
HTR (this work) | 1.17 | 0.041 |
HPR (this work) | 202.77 | 0.76 |
Table 1: Summary of values for telescope areas and time resolutions used for
fig. 5. The CTA requirements are used for the CTA SST and MST (we assume a
Gaussian distribution of photon arrival times for the SST and a rectangular
distribution for the MST), for the CTA LST we refer to a private communication
with D. Mazin. For the PW telescopes at least $\lambda/4$ precision is assumed
due to their imaging capability, and the timing resolution of photons
impinging on the detector is dominated by jitter introduced by atmospheric
fluctuations (Cavazzani et al., 2012).

Table 1 displays a few telescope types used in (upcoming) telescope arrays,
two smaller optical telescopes as well as the the virtual telescopes of this
work. The sizes and the telescope’s optical time spreads are given. The data
is visualized in fig. 5.

Telescopes suitable for HTR are located in the lower left corner of the plot
providing small collection areas but high isochronicity and thus small optical
time spread. Values for their isochronicity are assumed to be dominated by
atmospheric turbulence as described in Cavazzani et al. (2012), since the
telescope itself displays imaging capabilities with unisochronicity
$<\lambda/4$. The time resolution is thus dominated by the single photon
detection equipment rather than the optical time spread of the telescope
itself. The HPR setup can be operated at very high photon rates, which doesn’t
intrinsically limit the collection area, however, usually such telescopes are
parts of Cherenkov telescope arrays and have an unisochronicity of a few
nanoseconds (even though with MAGIC and the CTA LST also high time-resolution
Cherenkov telescopes exist/will be built). These telescopes are mainly located
in the upper right corner in fig. 5, their time resolutions are dominantly
given by the telescope’s unisochronicities. Due to the large point spread
function (PSF) of Cherenkov telescopes, intensity interferometry at such
telescopes would suffer from night sky background to a higher extent than for
telescope types with a smaller PSF. The large night sky background
contribution might set a limiting magnitude for observations much smaller than
the limiting magnitudes estimated from our background-free laboratory
measurements, especially for bright full-moon nights. Further cosmic ray
produced Cherenkov light background could influence S/N for faint limiting
magnitudes (Hanbury Brown, 1974).

A S/N map for a measurement time of 10 minutes of a mag 2 star is displayed In
fig. 5, see the right numbers of the colour bar. For a more intuitive view
when comparing the different telescopes the reader may also consider the left
numbers on the colour bar, denoting the maximum magnitude of a star which
satisfies $S/N\geq 5$ after $10$ minutes measurement time. However, please
note that the red data-points do not carry meaning in this interpretation.

It can be seen that with the HTR setup and two $1\text{\,}\mathrm{m}$ diameter
telescopes, S/N is comparable to the ones achievable with two small Cherenkov
telescopes such as the CTA SSTs. For mirrors on the order of
$3\text{\,}\mathrm{m}$ diameter or $>7\,$m2 of area, an HTR setup could
compete with measurements taken at CTA MSTs. Such larger mirror diameters are
feasible for use with single photon counting detectors for faint stars,by
utilizing multiple single photon detectors and beam splitters per telescope to
reduce the incident flux on each detector, or using detectors with dead times
on the order of $10\text{\,}\mathrm{ns}$. While such detectors with single
photon resolution are currently under development (Orlov et al., 2019a, b)
they are not yet widely available. Another promising detector type for
intensity interferometry are superconducting nanowire single photon detectors
(SNSPDs) due to their high quantum efficiency, low dead time, and high timing
resolution, for which efficient coupling to multimode fibers has been
demonstrated recently (see e.g. Chang et al. (2019)). While SNSPDs are
commercially avauliable, they would require effcient fiber coupling at the
telescope, requiring high-grade optics. This means, that currently the raw S/N
favours the current-correlating HPR approach for IACTs. However, an array of
$1\,$m telescopes used in conjunction with a HTR-like setup could
alternatively be used fruitfully in the future.

Assuming the time resolution of CTA-MSTs in combination with HPR intensity
interferometry electronics being not better than the MST’s requirement (an
admittedly conservative estimation), we compare use of 2 MSTs with HPR
electronics to an array of optical $1\text{\,}\mathrm{m}$ diameter telescopes.
For this comparison, the availability of different telescope baselines when
using an array of telescopes is not subject of investigation, instead the
individual telescopes add up to a combined light collection area. We calculate
an S/N value of $36.5$ for the HPR-MSTs, which corresponds to a measurement
with 26 pairs of $1\text{\,}\mathrm{m}$ telescopes (S/N = $36.0$) and thus
only 23 percent of the mirror area of two MSTs, demonstrating the potential of
high time resolution equipment for intensity correlation measurements.
Employing an HTR-like HBT instrumentation with higher quantum efficiency as
for example found in Si APDs, the number of $1\text{\,}\mathrm{m}$ telescopes
necessary to match the two MST S/N could be reduced by a factor of two. As
suggested in Trippe et al. (2014), the sensitivity of both intensity
interferometry setups could be enhanced by spectral and polarization
multiplexing, leading to accessible apparent magnitudes below 10. This
enhancement would probably scale especially favourably for a HTR-like approach
fielding narrow spectral filtering and a high quality beam that can be easily
controlled.

Note also that the small telescopes could be arranged in flexible geometries,
and might thus prove advantageous for the extraction of specific features of
astronomical objects. Further, they could realize larger baselines and thus
potentially obtain a higher angular resolution than CTA telescopes, again
since they are not constrained to the geometry of CTA. In particular, they
could be used to supplement the results gained from intensity interferometry
with CTA telescopes.

Given that CTA-S will be built, also the use of all 14 planned MSTs can be
considered. The HPR-S/N then corresponds to $255$ which would require $368$
HTR telescopes to achieve a similar S/N. Not even taking LSTs and SSTs with
the additive longer baselines into this consideration it demonstrates that
with good reason CTA-S is considered to be a most-promising concept for
successful intensity interferometry operations coming up in the future
decade(s).

Figure 5: Map of telescope properties. Two quantities are shown on the colour
bar. The right scale shows the $S/N$ as defined in this paper for a magnitude
2 star after $10\,$minutes of data accumulation using two telescopes of each
type at $100\,$percent spatial coherence. 5 contour lines correspond to a
factor of $10$ change. The red data points correspond to the measurements
discussed above, plotted as hypothetical star measurements with telescope
sizes matching the photon rates of the measurements. The left scale of the
colour bar shows the limiting magnitude of a star requiring $S/N\geq 5$ for a
$10\,$minutes measurement. Arrows pointing down denote an upper limit on the
time resolution, set by the current HTR time resolution in combination with
the Planewave telescopes. We assume a combined light collection and detection
efficiency of $0.1$ for the S/N and magnitude map. This is a rough assumption
not specific to any of the shown telescope arrays, thus there may be
differences to (possible) realizations for each of the arrays.

##  6 Outlook

With the presented laboratory measurements, we have proven that both portrayed
methods can work at the photon noise limit. The photon counting (HTR) ansatz
aims at highest time resolution with moderate photon rates. It achieved
$60\text{\,}\mathrm{ps}$ time resolution maintained over several hours. The
related noise level is $10^{-3}$. The photon current (HPR) ansatz aims at
highest PMT current and moderate time resolution. It achieved a rate of
$950\text{\,}\mathrm{MHz}$ and a stable noise level of $10^{-6}$ over several
minutes.
Given that both methods work at the ideal noise limits, they deliver
predictable signal-to-noise levels for astronomical observations with
reasonable telescope fields. In particular, it has been shown that small HTR
telescopes can compete with HPR at bigger IACTs. Based on the availability of
(future) IACT arrays HPR-implementations seem to be encouraging science cases.
However, further developments of HTR electronics to higher rate-capabilities
promise to stay competitive.

## Acknowledgements

This work was supported with a grant by the Deutsche Forschungsgemeinschaft
(‘Optical intensity interferometry with the H.E.S.S. gamma-ray telescopes’- FU
1093/3-1, Project number 426212122). The authors thank Heide Constantini and
Nolan Matthews for their helpful comments on the manuscript.

## Data Availability

The data that support the plots within this paper and other findings of this
study are available from the corresponding authors upon reasonable request.
For the HTR setup all intermediary correlation histograms are available in
time-intervals of $30\text{\,}\mathrm{s}$, however no raw photon event data
can be supplied. For the HPR all intermediary correlation histograms are
available in time-intervals of $3.436\text{\,}\mathrm{s}$. Due to the large
size of the digitized waveforms in excess of a TB, the raw data cannot be made
available online.

## References

  * Abeysekara et al. (2020) Abeysekara A. U., et al., 2020, Nature Astronomy
  * Acciari et al. (2019) Acciari V. A., et al., 2019, Monthly Notices of the Royal Astronomical Society, 491, 1540
  * Baines et al. (2007) Baines E. K., van Belle G. T., ten Brummelaar T. A., McAlister H. A., Swain M., Turner N. H., Sturmann L., Sturmann J., 2007, The Astrophysical Journal, 661, L195
  * Becker & Hickl GmbH (2017) Becker & Hickl GmbH 2017, HPM-100-06/07: Ultra-High Speed Hybrid Detectors for TCSPC, https://www.becker-hickl.com/wp-content/uploads/2018/11/db-hpm-06-07-v02.pdf
  * Bernlöhr et al. (2003) Bernlöhr K., et al., 2003, Astroparticle Physics, 20, 111
  * Bohec & Holder (2006) Bohec S. L., Holder J., 2006, The Astrophysical Journal, 649, 399
  * Brown & Twiss (1956) Brown R. H., Twiss R. Q., 1956, Nature, 177, 27
  * Brown et al. (1957) Brown R. H., Twiss R. Q., Lovell A. C. B., surName g., 1957, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 242, 300
  * Brown et al. (1967) Brown R. H., Davis J., Allen L., 1967, Monthly Notices of the Royal Astronomical Society, 137, 375
  * Cavazzani et al. (2012) Cavazzani S., Ortolani S., Barbieri C., 2012, Monthly Notices of the Royal Astronomical Society, 419, 2349
  * Chang et al. (2019) Chang J., et al., 2019, Appl. Opt., 58, 9803
  * Che et al. (2011) Che X., et al., 2011, The Astrophysical Journal, 732, 68
  * Dravins et al. (2013) Dravins D., LeBohec S., Jensen H., Nuñez P. D., 2013, Astroparticle Physics, 43, 331
  * Fox (2006) Fox M., 2006, Quantum optics: an introduction. Oxford Master Series in Physics, 6, Oxford University Press, USA, http://gen.lib.rus.ec/book/index.php?md5=AF074928C421227E6B7EBFCC30425DE2
  * Guerin et al. (2017) Guerin W., Dussaux A., Fouché M., Labeyrie G., Rivet J.-P., Vernet D., Vakili F., Kaiser R., 2017, Monthly Notices of the Royal Astronomical Society, 472, 4126
  * Guerin et al. (2018) Guerin W., Rivet J.-P., Fouché M., Labeyrie G., Vernet D., Vakili F., Kaiser R., 2018, Monthly Notices of the Royal Astronomical Society, 480, 245
  * Hanbury Brown (1974) Hanbury Brown R., 1974, The intensity interferometer; its application to astronomy. Taylor & Francis, Halsted Press, http://gen.lib.rus.ec/book/index.php?md5=d604cdb5cacf1ea40203b4b71256ea19
  * Hanbury Brown et al. (1974) Hanbury Brown R., Davis J., Allen L. R., 1974, Monthly Notices of the Royal Astronomical Society, 167, 121
  * Horch et al. (2021) Horch E. P., Weiss S. A., Klaucke P. M., Pellegrino R. A., Rupert J. D., 2021, Observations with the Southern Connecticut Stellar Interferometer. I. Instrument Description and First Results (arXiv:2112.07758)
  * Labeyrie et al. (2006) Labeyrie A., Lipson S. G., Nisenson P., 2006, An Introduction to Optical Stellar Interferometry. Cambridge University Press, Cambridge, doi:10.1017/CBO9780511617638, https://www.cambridge.org/core/books/an-introduction-to-optical-stellar-interferometry/2EA3ABDA8557CF3277063391C02E899D
  * Loudon (2000) Loudon R., 2000, The Quantum Theory of Light, 3rd ed. (Oxford Science Publications), 3 edn. Oxford University Press, USA, libgen.li/file.php?md5=ccc4545c9d7c0a3216fb446644bfcd47
  * Maier (2007) Maier G., 2007, VERITAS: Status and Latest Results (arXiv:0709.3654)
  * Mandel & Wolf (1995) Mandel L., Wolf E., 1995, Optical Coherence and Quantum Optics. Cambridge University Press, doi:10.1017/CBO9781139644105
  * Mourard, D. et al. (2015) Mourard, D. et al., 2015, A&A, 577, A51
  * Nuñez & Domiciano de Souza (2015) Nuñez P. D., Domiciano de Souza A., 2015, Monthly Notices of the Royal Astronomical Society, 453, 1999
  * Nuñez et al. (2012) Nuñez P. D., Holmes R., Kieda D., Rou J., LeBohec S., 2012, Monthly Notices of the Royal Astronomical Society, 424, 1006
  * Orlov et al. (2019a) Orlov D. A., Glazenborg R., Ortega R., Kernen E., 2019a, CEAS Space Journal, 11, 405
  * Orlov et al. (2019b) Orlov D. A., Glazenborg R., Ortega R., Kernen E., 2019b, in Itzler M. A., Bienfang J. C., McIntosh K. A., eds, Vol. 10978, Advanced Photon Counting Techniques XIII. SPIE, pp 93 – 100, https://doi.org/10.1117/12.2519061
  * Rivet et al. (2020) Rivet J.-P., et al., 2020, Monthly Notices of the Royal Astronomical Society, 494, 218
  * Roettenbacher et al. (2016) Roettenbacher R. M., et al., 2016, Nature, 533, 217
  * Rulten et al. (2016) Rulten C., Zech A., Okumura A., Laporte P., Schmoll J., 2016, Astroparticle Physics, 82, 36
  * Schlenstedt (2014) Schlenstedt S., 2014, in Ground-based and Airborne Telescopes V. p. 91450O
  * Shayduk et al. (2005) Shayduk M., Hengstebeck T., Kalekin O., Pavel N., Schweizer T., 2005, in 29th International Cosmic Ray Conference (ICRC29), Volume 5. p. 223
  * Tan et al. (2014) Tan P. K., Yeo G. H., Poh H. S., Chan A. H., Kurtsiefer C., 2014, The Astrophysical Journal, 789, L10
  * Thorlabs, Inc (2021) Thorlabs, Inc accessed: 06/2021, Unmounted Absorptive Neutral Density Filters - Transmission and Reflection Graphs, https://www.thorlabs.com/newgrouppage9.cfm?objectgroup_id=5011
  * Trippe et al. (2014) Trippe S., et al., 2014, Journal of The Korean Astronomical Society, 47, 235
  * Weekes et al. (2002) Weekes T., et al., 2002, Astroparticle Physics, 17, 221
  * Zampieri et al. (2016) Zampieri L., et al., 2016, in Malbet F., Creech-Eakman M. J., Tuthill P. G., eds, Vol. 9907, Optical and Infrared Interferometry and Imaging V. SPIE, pp 140 – 153, doi:10.1117/12.2233688, https://doi.org/10.1117/12.2233688
  * Zmija et al. (2021) Zmija A., Vogel N., Anton G., Malyshev D., Michel T., Zink A., Funk S., 2021, Monthly Notices of the Royal Astronomical Society, 509, 3113
  * ten Brummelaar et al. (2005) ten Brummelaar T. A., et al., 2005, The Astrophysical Journal, 628, 453–465


11institutetext: National University of Singapore. 21 Lower Kent Ridge Rd,
Singapore 119077 11email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
22institutetext: Agency for Science, Technology and Research (A*STAR),
Singapore Institute of Manufacturing Technology (SIMTech), Singapore 138634,
Singapore 22email<EMAIL_ADDRESS>
# Celestial Machine Learning

From Data to Mars and Beyond with AI Feynman

Zi-Yu Khoo 11 Abel Yang 11 Jonathan Sze Choong Low 22 Stéphane Bressan 11

###### Abstract

Can machine learning discover Kepler’s first law from data? We emulate
Johannes Kepler’s discovery of the equations of the orbit of Mars with the
Rudolphine tables using AI Feynman, a physics-inspired tool for symbolic
regression.

##  1 Introduction

In 2020, Silviu-Marian Udrescu and Max Tegmark introduced _AI Feynman_ [16], a
symbolic regression algorithm that could rediscover from data one hundred
equations from the _Feynman Lectures on Physics_ [3]. Although the authors
motivated their work with the example of Johannes Kepler’s successful
discovery of the orbital equations of Mars, to our knowledge, they did not
report an attempt to rediscover it with their algorithm. We show that AI
Feynman can emulate Kepler’s discovery of the orbital equation of Mars from
the Rudolphine tables.

The discovery of Kepler’s laws of planetary motion illustrates of the process
of science, encapsulating the principles of parsimony and physical
considerations. Prior to Kepler, early astrologers such as Nicolaus Copernicus
and Tycho Brahe hypothesized various models to explain the movement of
celestial bodies. Armed with the Prutenic tables, Copernicus modelled the
orbit of Mars as heliocentric, with a deferent having two epicycles [12].
However, Kepler, assistant to Tycho Brahe, had access to the best available
data collected in Europe. In 1627, Kepler compiled Brahe’s sightings of Mars
into a set of 180 heliocentric data of the position of Mars in the Rudolphine
tables. The translation of sightings to the Rudolphine tables only embeds
assumptions of heliocentrism and planarity of the orbit. Kepler could have
described the motions of Mars as an oval or added additional epicycles to the
Copernician model, but instead described it as elliptical in Astronomia nova
in 1609 [8]. We use AI Feynman to emulate Kepler’s discovery of the elliptical
orbital equation of Mars, from the Rudolphine tables.

The Rudolphine tables already embed assumptions regarding planarity and
heliocentricity of Mars’ orbit. To make further inferences using AI Feynman,
we can add biases to the data regarding its physical units. We have four
experiments. In the first experiment, AI Feynman is oblivious to biases. In
the second and third experiments, AI Feynman is biased through transforming
data which are angles and limiting the search space for orbital equations
respectively. The fourth experiment combines the biases of the second and
third. AI Feynman benefits from these biases, and produces the best result in
the fourth experiment. Information regarding the physical units of the data
likely also guided Kepler in his discovery of the elliptical orbit of Mars,
and in this way, AI Feynman emulates Kepler’s discovery from data in the
Rudolphine tables. In this paper, we present the design, results and
discussion of our experiments with AI Feynman.

##  2 Related Work

Finding an equation describing the orbit of Mars is a combinatorial challenge.
The task is NP-hard in principle due to its exponentially large search space
of equations [16]. To circumvent this, one may use universal function
approximators such as multilayer perceptron neural networks [5].
Alternatively, symbolic regressions search for a parsimonious and elegant form
of the unknown equation.

There are three main classes of symbolic regression methods [10]: regression-
based, expression tree-based and physics- or mathematics-inspired. We use AI
Feynman, a machine learning and physics-inspired algorithm [16].

Regression-based symbolic regression methods [10], given solutions to the
unknown equation, find the coefficients of a fixed basis that minimise the
prediction error. As the basis grows, the fit improves, but the functional
form of the unknown equation becomes less sparse or parsimonious. Sparse
regressions promote sparsity through regularisation, as proposed by Robert
Tibshirani [15] who used the $\mathit{l}_{1}$ norm, thus inventing the Lasso
regression. A state-of-the-art sparse symbolic regression approach is the
Sparse Identification of Nonlinear Dynamics by Steven Brunton et al. in [1].
It leverages regularisation and identifies equations of motion of a system
using a sparse regression over a chosen basis.

Committing to a basis limits the applicability of regression-based methods.
Expression tree-based symbolic regression methods based on genetic programming
[10] can instead discover the form and coefficients of the unknown equation.

Seminal work by John Koza et al. [9] represented each approximation of an
unknown equation as a genetic programme with a tree-like data structure, with
traits (or nodes in the tree) representing functions or operations, and
variables representing real numbers. The fitness of each genetic programme is
its prediction error. Fitter genetic programmes undergo a set of transition
rules comprising selection, crossover and mutation to iteratively find the
optimal equation form.

Genetic programmes may greedily mimic nuances of the unknown equation [14],
limiting generalisability. David Goldberg [4] used Pareto optimisation to
balance the objectives of fit and parsimony in symbolic regression. In each
iteration, the fittest genetic programmes lie on the non-dominated Pareto-
frontier. State-of-the-art symbolic regression using genetic programming
include `Eureqa` by Michael Schmidt and Hod Lipson [13] and `PySR` by Miles
Cranmer [2].

Expression tree-based methods do not guarantee that more accurate
approximations of an equation are symbolically closer to the truth. If an
expression tree-based method finds a reasonably accurate equation with wrong
functional form, it risks getting stuck near a local optimum [16]. Functions
of practical interest in physics exhibit simplifying properties such as
symmetry or separability [16]. Physics-inspired symbolic regression methods
leverage these simplifying properties to guarantee taking a step in the right
direction. Udrescu and Tegmark [16] use a neural network to test data
describing the unknown equation for simplifying properties and recursively
break the unknown equation into simpler unknown equations with fewer variables
to solve [16]. Each simpler unknown equation can be solved by regression with
a basis set of non-linear functions. AI Feynman then outputs a sequence of
increasingly complex equations that provide progressively better accuracy,
along a Pareto front, based on work by Goldberg [4] and Guido Smits [14]. We
use AI Feynman to rediscover the orbital equation of Mars.

##  3 Methodology

In publishing the Rudolphine tables, Kepler had already assumed planarity and
heliocentricity of the orbit of Mars. We experiment if AI Feynman performs
better with biases, based on our knowledge of the physical units of the data.

Informing a learning algorithm of physics amounts to introducing appropriate
biases that can steer the learning process towards identifying physically
consistent solutions according to George Karnadiakis [6]. Karniadakis
identifies three types of bias: observational, inductive and learning biases.
Observational biases are introduced directly through data that embody the
underlying physics or carefully crafted data augmentation procedures. Training
a machine learning model on such data allows it to learn an output that
reflects the physical structure of the data. Inductive biases correspond to
prior assumptions incorporated by tailored interventions to a machine learning
model architecture, so predictions are guaranteed to satisfy a set of given
physical laws. Learning biases can be introduced by appropriate choice of loss
functions, constraints and inference algorithms that can modulate the training
phase of an machine learning model to explicitly favour convergence towards
solutions that adhere to the underlying physics. We consider the introduction
of observational and inductive biases.

With knowledge that data is known to be an angle, only trigonometric functions
can transform it. We introduce an observational bias by applying the sine and
cosine functions to inputs of the unknown equation that are known to be
angles. The resulting numerical values, hence embodying the underlying
periodicity of the data, are input to AI Feynman. The observational bias
guides AI Feynman in finding an equation that reflects the periodic nature of
Mars’ orbit.

With knowledge that data are physical quantities, they cannot be transformed
by exponential and logarithmic functions, as these only transform
dimensionless quantities. We introduce an inductive bias by eliminating
candidate functions. For each simpler recursive unknown equation AI Feynman
has to solve, it transforms equations in the current Pareto front by one of
several non-linear functions. These non-linear functions include exponential,
logarithmic, trigonometric, polynomial and radical functions. The inductive
bias limits the search space to trigonometric, polynomial and radical
functions.

We conduct four experiments corresponding to four possible combinations of
observational and inductive biases with the AI-Feynman algorithm. Experiment 1
does not use any bias. Experiment 2 and 3 only use observational and inductive
bias respectively. Experiment 4 uses both observational and inductive bias.

While AI Feynman explores the Pareto front, Kepler may have instead made use
of thought experiments to hypothesize an elliptical orbit. Fitting data from
the Rudolphine tables to the equation of an ellipse using non-linear least
squares returns the coefficients representing the eccentricity and semi-major
axis. These are $0.0926$ and $1.5235$ respectively. For reference, the
National Aeronautics and Space Administration suggest an $0.0934$ and $1.5237$
respectively [11].

##  4 Performance Evaluation

In the Rudolphine tables [7], the table titled Tabula Aequationum MARTIS, or
Table of Corrections for Mars, contains four columns of data: Anomalia
eccentri, Intercolumnium, Anomalia coaequata, and Intervallu. These columns
represent the eccentric anomaly, an interpolating factor, the coequated or
true anomaly, and the distance between the Sun and Mars respectively. The full
Rudolphine tables (snippet found in Figure 1) was digitised for this
experiment.

Figure 1: The four columns of data provided in the Rudolphine Tables.

We apply AI Feynman to the data of Intervallu and Anomalia coaequata to
recover the equation of orbit for Mars. Anomalia coaequata is an angle in
degrees minutes seconds, which we convert to decimal degrees. Intervallu is
the distance between the Sun and Mars, which we scale from a magnitude of
$E+05$ to $E+00$. The code for AI Feynman, with minor modifications to embed
observational and inductive biases, is at https://github.com/zykhoo/AI-
Feynman.

We compare equations along the AI Feynman Pareto frontier with the orbital
equation of Mars in Equation  ‣ 4. $r$ is the Intervallu and $\theta$ is the
Anomalia coaequata. $\epsilon$ is the eccentricity of the ellipse, and $a$ is
the semi-major axis.

$\displaystyle r=\frac{a}{1+\epsilon\times cos\theta}$ (0)

We also present the mean description length loss[16, 17] (DL) computed between
each predicted and true Intervallu. It minimises the geometric mean instead of
the arithmetic mean, which encourages improving already well-fit points [17].

For Experiments 1 and 3, the inputs to AI Feynman are $\theta$ and $r$. For
Experiments 2 and 4, the inputs to AI Feynman are $\cos(\theta)$,
$\sin(\theta)$ and $r$. We traverse the equations along the Pareto frontier
returned by AI Feynman in increasing goodness of fit and increasing complexity
(or equivalently decreasing parsimony) and present them. The results of
Experiments 1, 2, 3 and 4 are presented in Tables 1, 2, 3 and 4 respectively.
We omit results independent of $\theta$.

Eqn No. | Equation | DL
---|---|---
(1a) |  $\displaystyle r=\frac{4}{3}-0.09\times\theta^{2}$ | 24.976
(1b) |  $\displaystyle r=(2.78-0.26\times\theta^{2})^{0}.5$ | 24.926
(1c) |  $\displaystyle r=\arccos{(-0.02\times\theta^{3}+0.09\times\theta^{2}-0.1)}$ | 23.577
(1d) |  $\displaystyle r=\frac{1}{-0.01\times\theta^{3}+0.04\times\theta^{2}+0.6}$ | 22.515
(1e) |  $\displaystyle r=(0.01\times\theta^{3}-0.04\times\theta^{2}+1.29)^{2}$ | 22.273
(1f)0 |  $\displaystyle r=\arccos{(-0.02\times\theta^{3}+0.09\times\theta^{2}+0.01\times\theta-0.1)}$ | 21.356
(1g)1 |  $\displaystyle r=\log(0.09\times\theta^{3}-0.38\times\theta^{2}-\frac{1}{9}\times\theta+5.3)$ | 20.841
(1h)2 |  $\displaystyle r=0.02\times\theta^{3}-0.09\times\theta^{2}-0.01\times\theta+1.67$ | 20.238
Table 1: Results of Experiment 1, which took 748 seconds.

Eqn No. | Equation | MSE
---|---|---
(2a)4 |  $\displaystyle r=\log{\cos\theta+5}$ | 26.006
(2b)5 |  $\displaystyle r=\frac{1}{7}\times\cos\theta+1.5$ | 24.053
(2c)qn:exp2exp1 |  $\displaystyle r=1.5\times\exp{0.1\times\cos\theta}$ | 23.512
(2d)qn:exp2cor1 |  $\displaystyle\mathbf{r=\frac{1}{\frac{2}{3}-0.0556244812357114\times\cos\theta}}$ | 22.857
(2e)qn:exp2exp2 |  $\displaystyle r=1.5119670200057298\times\exp{(0.1\times\cos\theta)}$ | 22.457
(2f)0 |  $\displaystyle r=1.510965630582+(\cos\theta/(\sin\theta+6))$ | 21.070
(2g)qn:exp2exp3 |  $\displaystyle r=1.51366746425629\times exp(0.0931480601429939\times\cos\theta)$ | 20.762
(2h)qn:exp2cor2 |  $\displaystyle\mathbf{r=\frac{1}{0.662428796291351-0.0612906403839588\times\cos\theta}}$ | 19.781
(2i)qn:exp2cor3 |  $\displaystyle\mathbf{r=(0.662428796291351-0.0612906403839588\times\cos\theta)^{-1.00133872032166}}$ | 12.211
Table 2: Results of Experiment 2, which took 1451 seconds.

Eqn No. | Equation | MSE
---|---|---
(3a)0 |  $\displaystyle r=\frac{4}{3}-0.09\times\theta^{2}$ | 24.976
(3b)1 |  $\displaystyle r=(2.78-0.25\times\theta^{2})^{0}.5$ | 24.842
(3c)2 |  $\displaystyle r=\arccos{(-0.02\times\theta^{3}+0.09\times\theta^{2}-0.1)}$ | 23.577
(3d)3 |  $\displaystyle r=\frac{1}{-0.01\times\theta^{3}+0.04\times\theta^{2}+0.6}$ | 22.515
(3e)4 |  $\displaystyle r=(0.01\times\theta^{3}-0.04\times\theta^{2}+1.29)^{2}$ | 22.273
(3f)5 |  $\displaystyle r=\arccos{(-0.02\times\theta^{3}+0.09\times\theta^{2}+0.01\times\theta-0.1)}$ | 21.356
(3g)7 |  $\displaystyle r=0.02\times\theta^{3}-0.09\times\theta^{2}-0.01\times\theta+1.67$ | 20.238
Table 3: Results of Experiment 3, which took 621 seconds.

Eqn No. | Equation | MSE
---|---|---
(4a)0 |  $\displaystyle r=\frac{1}{7}\times\cos\theta+1.5$ | 24.053
(4b)qn:exp4x0x11 |  $\displaystyle r=\cos\theta/(\sin\theta+6)+1.5$ | 23.617
(4c)1 |  $\displaystyle r=\arccos(0.0420224035468255-0.142857142857143\times\cos\theta)$ | 23.392
(4d)qn:exp4cor1 |  $\displaystyle\mathbf{r=\frac{1}{\frac{2}{3}-0.0566732120453772\times\cos\theta}}$ | 22.575
(4e)qn:exp4x0x12 |  $\displaystyle r=1.511006320056+(\cos\theta/(\sin\theta+6))$ | 21.089
(4f)2 |  $\displaystyle r=\tan(0.0425049090340329\times\cos\theta+0.986141372332807)$ | 20.057
(4g)3 |  $\displaystyle r=\tan(0.0427569970488548\times\cos\theta+0.98658412694931)$ | 20.021
(4h)qn:exp4cor2 |  $\displaystyle\mathbf{r=\frac{1}{0.662420213222504-0.0612917765974998\times\cos\theta}}$ | 19.747
(4i)qn:exp4cor3 |  $\displaystyle\mathbf{r=(0.662420213222504-0.0612917765974998\times\cos\theta)^{-1.00130701065063}}$ | 12.208
Table 4: Results of Experiment 4, which took 1184 seconds.

Experiment 1 does not use any bias. The search space is big and AI Feynman
does not find an equation form that matches the orbital equation of Mars along
its Pareto frontier. In Experiment 2, AI Feynman makes use of an observational
bias. As the input data to AI Feynman embodies the underlying periodicity of
the data, AI Feynman can use this information to guide its search for an
equation that reflects the periodic structure of Mars’ orbit. Therefore, three
out of nine equations along the Pareto front have an equation form that
matches the true orbit of Mars. These are Equations LABEL:eqn:exp2cor1,
LABEL:eqn:exp2cor2 and LABEL:eqn:exp2cor3. In Experiment 3, AI Feynman makes
use of an inductive bias. While the search space for AI Feynman is smaller, an
inductive bias does not guide its search for an equation that reflects the
periodic structure of Mars’ orbit. Therefore AI Feynman does not find an
equation form that matches the orbital equation of Mars along its Pareto
frontier. In Experiment 4, AI Feynman makes use of both an observational and
an inductive bias. AI Feynman makes use of the observational bias to guide its
search for an equation that reflects the periodic structure of Mars’ orbit. It
also makes use of the inductive bias to limit the search space, resulting in
fewer equations along the Pareto front. Therefore, three out of ten equations
along the Pareto front have an equation form that matches the true orbit of
Mars. These are Equations LABEL:eqn:exp4cor1, LABEL:eqn:exp4cor2 and
LABEL:eqn:exp4cor3.

Experiments 1 and 2 highlight the importance of an observational bias in
guiding AI Feynman. In Experiment 1, none of the equations along the Pareto
front match the orbital equation for Mars, compared to three out of eleven
equations along the Pareto front for Experiment 2. However, as the
observational bias doubles the number of inputs to AI Feynman, it takes
approximately twice as long to run. This is because AI Feynman recurses
through each input. The depth of the recursion is doubled when the number of
inputs is doubled.

Experiments 2 and 4 highlight the importance of an inductive bias in limiting
the search space for AI Feynman. In Experiment 2, we can observe many
equations have one of two common forms. Three of them utilise an exponential
function applied to $\cos\theta$ (Equations LABEL:eqn:exp2exp1,
LABEL:eqn:exp2exp2 and LABEL:eqn:exp2exp3). Another three utilise an inverse
function applied to $\cos\theta$ (Equations LABEL:eqn:exp2cor1,
LABEL:eqn:exp2cor2 and LABEL:eqn:exp2cor3), which matches the true orbit of
Mars. However, in Experiment 4, the equation form with an inverse function
applied to $\cos\theta$ in (Equations LABEL:eqn:exp4cor1, LABEL:eqn:exp4cor2
and LABEL:eqn:exp4cor3) matches the true obit of Mars, and is also the most
prevalent. Therefore, AI Feynman, augmented with both an observational and
inductive bias, is best able to rediscover Kepler’s first law for the orbit of
Mars. The inductive bias also reduces the time taken for AI Feynman to run.
This is because the search space for AI Feynman is limited, and time is saved
from having to search a smaller search space.

We also observe that Equations LABEL:eqn:exp2cor2, LABEL:eqn:exp2cor3,
LABEL:eqn:exp4cor2, and LABEL:eqn:exp4cor3, with forms that match the true
orbit of Mars, have the lowest mean description length loss of less than $20$.

Lastly, we observe that Equations LABEL:eqn:exp2cor2 and LABEL:eqn:exp4cor2
suggest $a=1.52$ and $\epsilon=0.0925$ similar to the values suggested in the
Rudolphine tables.

##  5 Conclusion

We have successfully shown that AI Feynman can rediscover from the Rudolphine
Tables the equation of Kepler’s first law for the planet Mars, given
information regarding physical quantities of the data in the form of
observational and inductive biases. The discovery of physical laws is a bi-
optimisation problem of parsimony and accuracy, that can be guided by physical
units of the data available. AI Feynman is able to emulate Kepler’s discovery
of the orbital equation of Mars because it implements an optimisation of both
parsimony and accuracy, and can be guided by information regarding the
physical units of the data.

As future work, we are looking into how AI Feynman can repeat this discovery
process directly from sightings of Mars and the Sun from Earth. We use a
modern reproduction of these sightings from the National Aeronautics and Space
Administration’s Horizons system. This challenges AI Feynman to perform a
change from the geocentric to heliocentric reference frame, and we are
investigating how this change can be incorporated within its algorithm.
Lastly, we are experimenting with the planet Mercury, which has a precessing
orbit.

## References

  * [1] Brunton, S.L., Proctor, J.L., Kutz, J.N.: Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences 113(15), 3932–3937 (2016)
  * [2] Cranmer, M.: Pysr: Fast & parallelized symbolic regression in python/julia (Sep 2020), http://doi.org/10.5281/zenodo.4041459
  * [3] Feynman, R.P.: The Feynman lectures on physics. Reading, Mass. : Addison-Wesley Pub. Co., c1963-1965. (1965 c1963)
  * [4] Goldberg, D.E.: Genetic Algorithms in Search, Optimization and Machine Learning. Addison-Wesley Longman Publishing Co., Inc., USA, 1st edn. (1989)
  * [5] Hornik, K., Stinchcombe, M., White, H.: Multilayer feedforward networks are universal approximators. Neural Networks 2(5), 359–366 (Jan 1989)
  * [6] Karniadakis, G.E., Kevrekidis, I.G., Lu, L., Perdikaris, P., Wang, S., Yang, L.: Physics-informed machine learning. Nature Reviews Physics 3(6), 422–440 (Jun 2021)
  * [7] Kepler, J., Brahe, T., Eckebrecht, P.: Tabulæ Rudolphinæ, quibus astronomicæ scientiæ, temporum longinquitate collapsæ restauratio continetur (1627)
  * [8] Kepler, J., Donahue, W.H.: New Astronomy. Cambridge University Press (1992)
  * [9] Koza, J.R.: Genetic Programming: On the Programming of Computers by Means of Natural Selection. MIT Press, Cambridge, MA, USA (1992)
  * [10] Makke, N., Chawla, S.: Interpretable Scientific Discovery with Symbolic Regression: A Review (11 2022)
  * [11] National Aeronautics and Space Administration: Mars fact sheet. https://nssdc.gsfc.nasa.gov/planetary/factsheet/marsfact.html (2022)
  * [12] Rosen, E.: The commentariolus of copernicus. Osiris 3, 123–141 (1937). https://doi.org/10.1086/368473, https://doi.org/10.1086/368473
  * [13] Schmidt, M., Lipson, H.: Distilling free-form natural laws from experimental data. Science 324(5923), 81–85 (2009)
  * [14] Smits, G.F., Kotanchek, M.: Pareto-Front Exploitation in Symbolic Regression, pp. 283–299. Springer US, Boston, MA (2005)
  * [15] Tibshirani, R.: Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, Series B 58, 267–288 (1994)
  * [16] Udrescu, S.M., Tegmark, M.: AI Feynman: A physics-inspired method for symbolic regression. Science Advances 6(16) (2020)
  * [17] Wu, T.: Intelligence, physics and information – the tradeoff between accuracy and simplicity in machine learning (2020)

## Appendix

Eqn No. | Eqn | MSE
---|---|---
1 | $asin(-666.000000000000\times((sin(pi)+1)-1))$ | 2.33
2 | $1.50000000000000$ | 0.0106
3 | $pi/2$ | 0.0123
4 | $1.65306122448980$ | 0.0268
5 | $1.66666666666667-0.09\times x_{0}^{2}$ | 0.048
6 | $(2.78-0.26\times x_{0}^{2})^{0}.5$ | 0.0767
7 | $acos(-0.02\times x_{0}^{3}+0.09\times x_{0}^{2}-0.1)$ | 0.000169
8 | $1/(-0.01\times x_{0}^{3}+0.04\times x_{0}^{2}+0.6)$ | 0.000706
9 | $(0.01\times x_{0}^{3}-0.04\times x_{0}^{2}+1.29)^{2}$ | 0.000458
10 | $acos(-0.02\times x_{0}^{3}+0.09\times x_{0}^{2}+0.01\times x_{0}-0.1)$ | 4.7e-05
11 | $log(0.09\times x_{0}^{3}-0.38\times x_{0}^{2}-0.111111111111111\times x_{0}+5.3)$ | 1.07e-05
12 | $0.02\times x_{0}^{3}-0.09\times x_{0}^{2}-0.01\times x_{0}+1.67$ | 4.41e-05
Table 5: Experiment 1 Results Figure 2: A plot of the results of Experiment 1. The equations correspond to those presented in Table 1. The y-axis represents the Intervallu and x-axis represents the Anomalia coaequata. The true values for the Intervallu from the Rudolphine tables are also plotted and labeled ”original”. Eqn No. | Eqn | MSE
---|---|---
1 | $1.50000000000000$ | 0.0106
2 | $log(x_{0}+5)$ | 0.0092
3 | $0.142857142857143\times x_{0}+1.5$ | 0.000309
4 | $1.5\times exp(0.1\times x_{0})$ | 0.00021
5 | $acos(0.0418258837357514-0.142857142857143\times x_{0})$ | 0.000166
6 | $1/(0.666666666666667-0.0571196591636008\times x_{0})$ | 0.000212
7 | $1.5082674607662332\times exp(0.1\times x_{0})$ | 7.44e-05
8 | $1.510965630582+(x_{0}/((((((x_{1}+1)+1)+1)+1)+1)+1))$ | 0.000179
9 | $1.51366746425629\times exp(0.0931480601429939\times x_{0})$ | 5.39e-06
10 | $tan(0.0427570976316929\times x_{0}+0.986583888530731)$ | 1.98e-06
11 | $tan(0.0428397443727006\times x_{0}+0.986126406475687)$ | 4.21e-06
12 | $1/(0.662416338920593-0.0612923018634319\times x_{0})$ | 7.26e-07
13 | $(0.662416338920593-0.0612923018634319\times x_{0})^{(}-1.0012925863266)$ | 6.01e-10
Table 6: Experiment 2 Results Figure 3: A plot of the results of Experiment 2. The equations correspond to those presented in Table 2. The y-axis represents the Intervallu and x-axis represents the Anomalia coaequata. The true values for the Intervallu from the Rudolphine tables are also plotted and labeled ”original”. Eqn No. | Eqn | MSE
---|---|---
1 | $0$ | 2.33
2 | $1.50000000000000$ | 0.0106
3 | $pi/2$ | 0.0123
4 | $1.65306122448980$ | 0.0268
5 | $1.66666666666667-0.09\times x_{0}^{2}$ | 0.048
6 | $(2.78-0.26\times x_{0}^{2})^{0}.5$ | 0.0767
7 | $acos(-0.02\times x_{0}^{3}+0.09\times x_{0}^{2}-0.1)$ | 0.000169
8 | $1/(-0.01\times x_{0}^{3}+0.04\times x_{0}^{2}+0.6)$ | 0.000706
9 | $(0.01\times x_{0}^{3}-0.04\times x_{0}^{2}+1.29)^{2}$ | 0.000458
10 | $acos(-0.02\times x_{0}^{3}+0.09\times x_{0}^{2}+0.01\times x_{0}-0.1)$ | 4.7e-05
11 | $(0.06\times x_{0}^{3}-0.26\times x_{0}^{2}-0.05\times x_{0}+2.78)^{0}.5$ | 5.61e-06
12 | $0.02\times x_{0}^{3}-0.09\times x_{0}^{2}-0.01\times x_{0}+1.67$ | 4.41e-05
Table 7: Experiment 3 Results Figure 4: A plot of the results of Experiment 3. The equations correspond to those presented in Table 3. The y-axis represents the Intervallu and x-axis represents the Anomalia coaequata. The true values for the Intervallu from the Rudolphine tables are also plotted and labeled ”original”. Eqn No. | Eqn | MSE
---|---|---
1 | $0$ | 2.33
2 | $1.50000000000000$ | 0.0106
3 | $1.56250000000000$ | 0.0115
4 | $0.142857142857143\times x_{0}+1.5$ | 0.000309
5 | $x_{0}/(x_{1}+6)+1.5$ | 0.000416
6 | $1/(0.666666666666667-0.0557172402393568\times x_{0})$ | 0.000265
7 | $1.510957104465+(x_{0}/((((((x_{1}+1)+1)+1)+1)+1)+1))$ | 0.000179
8 | $1/(0.662428081035614-0.0612907484173775\times x_{0})$ | 7.75e-07
9 | $0.140863761305809\times x_{0}-0.0146051803603768\times x_{1}+1.52623379230499$ | 1.1e-06
10 | $(0.662428081035614-0.0612907484173775\times x_{0})^{(}-1.001335978508)$ | 6.01e-10
Table 8: Experiment 4 Results Figure 5: A plot of the results of Experiment 4.
The equations correspond to those presented in Table 4. The y-axis represents
the Intervallu and x-axis represents the Anomalia coaequata. The true values
for the Intervallu from the Rudolphine tables are also plotted and labeled
”original”.