2023

[1]Xinyang Ge

[1]Sergio Maldonado

1]Faculty of Engineering and Physical Sciences, University of Southampton,
Southampton, SO16 7QF, UK

2]Photonics Laboratory, ETH Zurich, Zurich, CH-8093, Switzerland

# High-precision interferometric measurement of slow and fast temperature
changes in static fluid and convective flow

<EMAIL_ADDRESS>Joanna A. Zielińska<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [

###### Abstract

We explore the strengths and limitations of using a standard Michelson
interferometer to sample line-of-sight-averaged temperature in water via two
experimental setups: slow-varying temperature in static fluid and fast
temperature variations in convective flow. The high precision of our
measurements (a few mK) is enabled by the fast response time and high
sensitivity of the interferometer to minute changes in the refractive index of
water caused by temperature variations. These features allow us to detect the
signature of fine fluid dynamical patterns in convective flow in a fully non-
intrusive manner. For example, we are able to observe an asymmetry in the
rising thermal plume (i.e. an asynchronous arrival of two counter-rotating
vortices at the measurement location), which is not possible to resolve with
more traditional (and invasive) techniques, such as RTD (Resistance
Temperature Detector) sensors. These findings, and the overall reliability of
our method, are further corroborated by means of Particle Image Velocimetry
and Large Eddy Simulations. While this method presents inherent limitations
(mainly stemming from the line-of-sight-averaged nature of its results), its
non-intrusiveness and robustness, along with the ability to readily yield
real-time, highly accurate measurements, render this technique very attractive
for a wide range of applications in experimental fluid dynamics.

###### keywords:

Interferometry, Convective flow, Optical metrology, Temperature measurement

##  1 Introduction

Reliable measurement of temperature in fluids is important for a wide range of
applications, from industrial (e.g. power plants, heat exchangers, chemical
industries) to basic scientific research (e.g. Ross et al, 2001; Abram et al,
2018; Laffont et al, 2018). Fluid temperature is conventionally measured via
intrusive instrumentation. Thermocouples, thermistors and resistance
temperature detectors (RTDs) are popular examples of instruments which measure
fluid temperature by immersion of a physical probe. While this may not be
problematic for many industrial applications, when studying fluid dynamics
from a scientific perspective, this intrusiveness may be undesirable since the
probes may disturb local flow patterns. Also, these sensors typically have
moderate to poor accuracy. While RTDs are generally considered accurate for
most applications, their slow response time prevents them from capturing fast
changes in fluid temperature (Goumopoulos, 2018), of the type expected in
turbulent flows. Another class of temperature sensors are based on optical
fibers (Childs et al, 2000; Fernandez-Valdivielso et al, 2003; Rizzolo et al,
2017; Kim et al, 2018), which offer advantages such as small dimensions,
capability of multiplexing, chemical inertness, and immunity to
electromagnetic fields (Roriz et al, 2020). It must be noted, however, that
while optical fibre based sensors can be as small as a fraction of a
millimetre in diameter (Drusová et al, 2021), they must be immersed in the
test fluid, making them invasive instruments. Another intrusive technique is
acoustic thermometry, which relies on the measurement of the speed of sound in
a given fluid (Moldover et al, 2014; Wang et al, 2018). The response time of
this technique is determined by the time traveled by a probe sound wave, in
turn given by the thickness of the sample volume. This method can attain a
precision of up to $\pm$0.015 K with a sampling frequency of 5 Hz (Wang et al,
2018).

When it is not possible or desirable for sensors to be in direct contact with
the test fluid, a non-intrusive method may be used, such as laser-induced
fluorescence (LIF) or interferometry. LIF (originally proposed by Tango et al
1968) works by using laser light to excite a fluorescent dye, which
subsequently fluoresces at a different wavelength, revealing temperature
changes (Banks et al, 2019). LIF has been employed in numerous studies to
understand heat transfer in convective flow (e.g. Grafsrønningen and Jensen,
2012; Park et al, 2019; Kashanj and Nobes, 2023), and can be used to obtain
both two- and three-dimensional images (Taylor and Lai, 2021), with accuracy
in thermofluid applications of up to 0.17 K (Sakakibara and Adrian, 2004).

Interferometric techniques for visualisation and analysis of transparent fluid
flows are non-invasive, fast and extremely sensitive. These include
holographic interferometry (HI) and standard interferometry (SI).
Inteferometric methods extract fluid properties (e.g. temperature) from
spatial and/or temporal interference between two laser beams, one of which
passes through the test fluid (the other one being the reference). In HI, the
interference pattern between two laser beams is recorded on a camera and used
to reconstruct the optical phase distribution in the observed area, which
directly maps onto the 2D temperature distribution within the field of view of
the instrument. Several authors have employed HI to understand convection-
based heat transfer phenomena around horizontally placed heating rods (Herraez
and Belda, 2002; Ashjaee et al, 2007; Ashjaee and Yousefi, 2007; Narayan et
al, 2017) as well as temperature profiles during liquid cooling and heating
processes (Wu et al, 2013; Guerrero-Mendez et al, 2016; Wang et al, 2022).
These investigations have successfully extracted temperature variations in the
fluid, with reported accuracies ranging from 0.20 K (Guerrero-Mendez et al,
2016) to 0.42 K (Wu et al, 2013). Most studies have been limited to laminar
flows, and the reported experimental data are recorded under steady-state
conditions after the initial transients have decayed. This may be due to the
fact that, in dynamic measurements, HI performance is limited, as large phase
gradients across the beam (which may occur in turbulent flow settings) may
lead to non-distinguishable interference structure and errors in the
interpretation of fringe patterns (Wilkie and Fisher, 1963; Lira and Vest,
1987).

Unlike HI, SI does not provide spatial information about the distribution of
fluid properties such as temperature. Instead, in SI the temporal interference
between probe and reference beams is observed using a single-pixel detector,
which provides necessary information to measure temporal variations of the
optical path length in one location. Continuous modulation of the reference
arm length enables monitoring of the interference visibility, which provides
information about the performance of the instrument (which may be affected by
phase gradients perpendicular to the probe beam’s direction of propagation).
This is advantageous for highly non-uniform flow settings, permitting
measurements in the presence of large phase gradients (albeit with reduced
sensitivity). What is more, SI is relatively simple and economic, with data
processing requirements being significantly less computationally demanding
than those of HI. These characteristics are desirable for real-time
measurements of rapid temperature variations in fluids.

While SI is a well-established technique with a wide range of applications in
metrology and other branches of science (e.g. Smith and Dobson, 1989; Freise
et al, 2009; Ito et al, 2020), to the best of the authors’ knowledge, very few
researchers have employed this technique to sample temperature in fluid
dynamical setups. For example, Tomita et al (2006) and Mahdieh and Nazari
(2013) implemented SI to measure very slow water temperature changes in
uniformly heated water. However, their detection method (involving the imaging
of interference fringes using a camera) enables only a modest measurement
accuracy; e.g. of $\sim$6% in the case of Mahdieh and Nazari (2013). Zhang et
al (2019) used SI to measure salinity and temperature in slowly cooling water
(over a 10 hour period). Their method, necessitating empirical expressions
(relating temperature, salinity and refractive index), yielded measurements
with an accuracy of 0.12 K relative to the reference instrument (a
thermistor).

In this paper, we explore the potential and limitations of SI in experimental
fluid dynamics by applying it to the study of two different problems: slow-
varying temperature in static water and convective flow in water induced by a
horizontal heating rod. The latter, in particular, is a novel application of
SI. Results are compared against and complemented by well-established
experimental and numerical methods, such as RTD sensors, Particle Image
Velocimetry and Large Eddy Simulations. The structure of the paper is as
follows. In Section 2, we explain the experimental setup and the procedure to
extract temperature changes from shifts in interference patterns. Then,
Section 3 presents the results of a slowly varying temperature experiment in a
static fluid, which establishes the method’s reliability and sensitivity.
Next, in Section 4, we apply our technique to a convective flow, and compare
the obtained results against those from RTD sensors, further enhancing our
analysis by the use of Particle Image Velocimetry and Large Eddy Simulations.
Section 5 presents a brief discussion on the key merits and limitations of the
proposed technique, as well as an outlook for potential future applications
and improvements. Concluding remarks are presented in Section 6.

##  2 Methodology

###  2.1 Interferometry: basic theory

Figure 1: a) Experimental setup of the interferometer. The reference arm
length is controlled using a piezoeletric actuator driven by a voltage signal
$V(t)$. The probe arm passes twice through the water tank. Abbreviations: PZT
- piezoeletric actuator, M - mirror, BS - non-polarising 50:50 beamsplitter,
PD - photodetector, DAQ - data acquisition device. b) Oscilloscope view of the
PZT driving voltage, $V(t)$, and signal recorded on the detector PD, $S(t)$.
c) 3D render of the setup, showing approximate position of the heating rod and
RTD sensors inside the water tank, as well as a close-up schematic view of the
overlap between probe and reference beams at the PD. d) Illustration of the
method for extracting temperature changes from interference pattern shifts;
$S(t)$ and $S(t+\tau)$ are PD-detected signals during subsequent piezo
strokes.

In most interferometers, light from a single coherent source is split into two
beams that travel via different optical paths, which are then recombined to
produce an interference pattern. The resulting interference fringes provide
information about the difference in optical path lengths with a precision of a
fraction of the wavelength of light.

Here, we measure water temperature using a Michelson interferometer. The probe
arm of the interferometer passes through the water tank, while the second
(reference) arm path bypasses the tank. The interferometer measures the phase
difference between the two arms, which depends on the refractive index of
water, in turn given by its temperature.

The experimental setup is shown in Fig. 1(a). A laser beam is split into two
equal-power components with the help of a beam-splitter. One of the beams
passes through the water tank and is reflected off a stationary mirror, while
the other one is reflected off a mirror mounted on a piezoelectric actuator,
fed by a voltage signal shown in Fig. 1(b). The interference signal $S(t)$ is
detected by a photodetector. The signal $S(t)$ consists of fringes, resulting
from linear phase ramps (i.e. linear piezo displacement in time; see $V(t)$
curve in Fig. 1b) introduced by the piezoelectric transducer (PZT). Each time
the PZT reverses direction, we observe a ‘turning point’ in the fringe pattern
(see times 100, 200 and 300 in Fig. 1b). It must be noted that the probe beam
passes twice through the test fluid, which effectively increases the
instrument’s sensitivity by a factor of $2$ (compared to a single-pass
configuration).

Interferometry does not provide absolute temperature measurements, but rather
variations with respect to time. We detect time variations of the temperature
(refractive index) of water by comparing the fringe patterns recorded in
different moments in time (consecutive PZT sweeps). A shift of the fringe
pattern with respect to the PZT drive signal allows us to deduce the phase
difference between the arms, $\phi(t)$, which is expressed as:

$\phi(t)=k_{0}z_{\rm pzt}(t)-2k_{0}\int_{0}^{l}n_{\rm w}(t)\,dz+\phi_{0}(t),$
(1)

where $k_{0}$ is the wave-vector in vacuum, $k_{0}z_{\rm pzt}(t)$ is the phase
shift introduced by the displacement of the PZT (proportional to $V(t)$), the
$z$ coordinate is aligned with the direction of beam propagation through the
water tank along a path of length $l$, and $n_{\rm w}(t)$ is the refractive
index of water (which may vary in time). The term $\phi_{0}$ contains
mechanical drifts which affect the lengths of the interferometer arms, as well
as phase difference introduced by variations in air temperature.

For small temperature variations $\delta T=T_{1}-T_{0}$, we can approximate
the dependence of water refractive index on temperature $T_{1}$ as

$n_{\rm w}(T_{1})\approx n_{\rm w}(T_{0})+\left.\left(\frac{\partial
n}{\partial T}\right)\right\rvert_{T_{0}}\delta T\,,$ (2)

where $\partial n_{\rm w}/\partial T=-9.6\times 10^{-5}~{}\rm{K}^{-1}$ within
the range of 285 to 305 K at 101,325 Pa (Fernández-Prini and Dooley, 1997).
Note that $\partial n/\partial T$ for air is approximately $-9\times
10^{-7}~{}{\rm K}^{-1}$ (Ciddor, 1996) (dry air, 450 ppm CO2), such that the
experiments described in this paper are much more sensitive to water than to
air temperature changes (by two orders of magnitude) that may occur in the
lab. Indeed, this is why we select water as the test fluid.

Note that the temperature variation deduced from $\phi(t)$ corresponds to the
average change along the beam path (inside the water).

The signal $S(t)$ detected by the photodiode (see PD in Fig. 1) and $\phi(t)$
are related via:

$S(t)=\eta\left[I_{1}+I_{2}+2\sqrt{I_{1}I_{2}}\gamma\cos(\phi(t))\right],$ (3)

where $I_{1}$ and $I_{2}$ are the optical intensities corresponding to the
probe and reference beams, respectively, and $\gamma$ is the absolute value of
the first order normalized correlation function between probe and reference
fields (Mandel and Wolf, 1995). The calibration factor $\eta$ converts the
signal to Volts. The two beams have approximately equal power, such that:

$I_{1}\approx I_{2}=\frac{S_{0}}{2\eta}\,,$ (4)

where $S_{0}$ is the DC part of the signal $S(t)$ (see Fig. 1b). The amplitude
of the interference fringes relative to $S_{0}$ (interference visibility) is
dictated by $\gamma$. Unit visibility ($\gamma=1$) occurs when the spatial
coherence (overlap between probe and reference beams) is perfect and the path
length difference between probe and reference beams is much smaller than the
coherence length. In this work the latter is always true, i.e. our anticipated
optical path length differences are always well within the coherence length,
which for our laser source is equal to the path length difference resulting
from a water temperature change of 59 K. The interference visibility in our
experiment is therefore limited by spatial misalignment between the two beams
and wavefront distortion introduced by water, which results in imperfect
overlap of probe and reference beams at the detector (see inset in Fig. 1c).

The method for extracting temperature changes from the fringe pattern is
depicted in Fig. 1d. We calibrate the PZT using the fact that a change of
(double-pass) path length difference of one wavelength is equivalent to the
shift of the sinusoidal signal pattern by one period. Therefore, a voltage
change of $|V_{1}-V_{0}|$ (see Fig. 1d) corresponds to a single-pass path
length difference of $\lambda/2$, yielding a calibration factor
$\alpha=\frac{\lambda}{2|V_{1}-V_{0}|}$. The temporal change in the refractive
index (integrated over $l$) corresponding to a voltage change $V_{2}-V_{1}$
can then be obtained from:

$\int_{0}^{l}\left(n_{\rm w}(t+\tau)-n_{\rm
w}(t)\right)\,dz=\alpha(V_{2}-V_{1}),$ (5)

where $V_{1}$ is the ‘$S_{0}$-crossing’ of the signal at time $t$ and $V_{2}$
corresponds to the subsequent piezo stroke at time $t+\tau$ (see Fig. 1d).
Finally, we can calculate the temperature change as:

$\widehat{\delta
T}=\widehat{T}(t+\tau)-\widehat{T}(t)=\frac{\alpha}{l}(V_{2}-V_{1})\left(\frac{\partial
n_{w}}{\partial T}\right)^{-1},$ (6)

where $\widehat{T}(t)$ denotes average water temperature along the beam
(corresponding to the average refractive index along the beam) and
$\widehat{\delta T}$ is the corresponding line-of-sight-average change in
temperature (but, for clarity, we drop the hat notation hereinafter). In our
experiments, we sample temperature changes sufficiently fast so that
$|V_{2}-V_{1}|<|V_{1}-V_{0}|$ is always satisfied.

##  3 Static fluid experiment

As a first assessment of the proposed technique, we compare its performance
against that of state-of-the-art RTD (Resistance Temperature Detector) sensors
for a long experiment where water temperature varies slowly. In order to
quantify the stability and sensitivity of the technique, we also analyse its
Allan variance.

###  3.1 Experimental setup

A schematic of the experimental setup is illustrated in Fig. 1. The outer
dimensions of the glass tank are $0.20\times 0.10\times 0.30$ m in length,
width and height, respectively. We use a collimated laser with a wavelength of
532 nm (Thorlabs, CPS532) and orient the tank such that the laser travels
along its length. The beam diameter is approximately 3 mm. The beam splitter
is a non-polarising cube beam splitter with a 50:50 (Reflectance :
Transmission) beam splitting ratio and the moveable mirror is powered by a
discrete piezo ring stack (Thorlabs, PK44LA2P2) with a maximum displacement of
9 $\mu$m when driven with no load. The piezo actuator is powered by a waveform
generator (ISO-TECH, GFG-8216A), which can produce different types of
waveforms (e.g. triangular, sinusoidal) with varying frequencies. Four Class
1/10 DIN RTD sensors (OMEGA, P-M-1/10-1/4-6-0-P-3) are evenly spaced inside
the water with their probes located near the laser path for comparison against
the line-of-sight-averaged results obtained from the interferometer; for the
same reason, all RTD measurements reported in this paper refer to the average
of all four sensors. The RTDs are connected to a data logger (DataTaker DT85M)
for recording temperature with a 4-wire configuration. Their accuracy is
$\pm(0.03+0.0005T_{s})\,{}^{\circ}C$, where $T_{s}$ is the surrounding fluid
temperature.

The tank is filled with 20 cm of water and the laser beam propagates
horizontally at 5 cm from the bottom of the tank. For this experiment, the
piezo is driven with a triangular signal at 40 Hz. The experiment runs for 13
hrs in order to compare the results from the proposed technique against
measurements from RTD sensors in a simple setting where water temperature
varies slowly due to ambient temperature changes. The experiments have been
performed remotely at night hours to minimise mechanical vibrations induced by
human movement. A vibration absorption mat is also placed underneath the
optical table to further reduce possible noise.

###  3.2 Results

The temperature variations measured by the RTD sensors and inferred from
interferometry are compared in Fig. 2. The observed drop in water temperature
is entirely due to the drop in room temperature over night (no additional
cooling system is employed). The sampling frequency of the data logger for the
RTD sensors is 1 Hz, and the reported value is the average of the four
sensors. For the interferometry results, we captured a triggered result based
on the voltage from the waveform generator every 6 seconds. The data
transmission rate of the oscilloscope (Tektronix, TDS-2004C) limits this
interval. While a different data logger could have been employed for
continuous data recording (see Section 4.1), the volume of data corresponding
to this 13 hrs long experiment would have been unnecessarily large for post-
processing purposes. Recalling that the proposed technique yields temperature
changes rather than absolute values, we have fixed the reference temperature
at $t=0$ as equal to that from the RTD sensors (i.e. their average). This is
the only calibration we carry out. Fig. 2 shows that the two curves agree
extremely well; the difference between the methods has an average of only
$2.47\times 10^{-4}$ K. Most noticeable differences occur right before the
auto-recalibration process of the data logger (see inset in Fig. 2). During
operation, the RTD data logger measures the amplifier’s internal ‘offset
voltage’. If it is found to have drifted by a specified amount ($1\mu V$), a
calibration cycle is performed (dataTaker, 2013). Regardless of the error
caused by the data logger drift, the difference between the two methods is
within the expected accuracy of the RTD sensors ($\sim 0.038$ K). This
comparison against state-of-the-art temperature sensors in a simple setting
serves as a confirmation of the high accuracy achieved by the proposed method.

Figure 2: Comparison of measurements from RTD sensors and interferometry for
slowly varying water temperature in static fluid. Differences are barely
noticeable (see inset).

###  3.3 Interferometer stability

We quantify the stability of the interferometer under hydrostatic conditions,
which determines the precision of our measurements for a given integration
time. A standard method to characterise a sensor’s stability is the two-sample
variance, also known as Allan variance (Allan, 1966), which quantifies noise
in the system and resulting measurement limitations. For the derivation of the
two-sample variance $\sigma^{2}_{T}(\tau)$ of a time trace $T(t)$ with length
$\mathcal{T}$, the time trace is split into $N$ equal blocks of length
$\tau=\mathcal{T}/N$. For each of these blocks an average value is denoted as
$\bar{T}_{k}^{(\tau)}$, where $k\in(1,...,N)$. The Allan variance is the
expectation value of the two consecutive values of the block average
$\bar{T}_{k}^{(\tau)}$ as a function of integration time $\tau$:

$\displaystyle\sigma^{2}_{T}(\tau)$
$\displaystyle=\left<\frac{1}{2}\left(\bar{T}_{k+1}^{(\tau)}-\bar{T}_{k}^{(\tau)}\right)^{2}\right>$
$\displaystyle\approx\frac{1}{2(N-1)}\sum_{k=1}^{N-1}\left(\bar{T}_{k+1}^{(\tau)}-\bar{T}_{k}^{(\tau)}\right)^{2}.$
(7)

Figure 3: Allan deviation plot for the static fluid setting; slope $\tau^{-1}$
shown as dashed line.

The result of the Allan variance analysis is shown in Fig. 3, which has been
calculated by acquiring data for 1 min with a $40$ Hz sampling frequency. We
observe that for integration times shorter than approximately 2 s, the signal
behaviour is dominated by high frequency noise. Increasing the averaging
window beyond $\sim 2$ s reveals the influence of other noise sources,
including long-term mechanical drifts in the optical path lengths. We estimate
that the precision of our instrument is approximately $1$ mK per piezo cycle
(at 40 Hz), and up to $5~{}{\rm\mu K}$ for 10 s integration time. Under the
assumption that the water temperature remains constant over a $1/40$ s
interval (a reasonable assumption given the experimental conditions), $1$ mK
can also be taken as a good estimate of the accuracy of the interferometer (at
a 40 Hz sampling rate).

##  4 Convective flow experiment

In this section we describe the measurements of rapid temperature variations
due to natural convection induced by a heating rod in the fluid.
Interferometry results are compared against those from RTD sensors. Our
interpretation of the measurements obtained is enhanced via numerical
simulations and Particle Image Velocimetry.

###  4.1 Experimental setup

The experimental setup is depicted in Figure 1. The data acquisition device
used for this experiment is a Red Pitaya 125-10 (a 10-bit data logger), which
facilitates continuous high-frequency data streaming for rapidly fluctuating
water temperature measurements. The tank is positioned on a jack lift, with a
rectangular opening in the optical table allowing for horizontal and vertical
adjustments to sample data from various locations without necessitating laser
realignment. The tank geometry is based on that reported by Park et al (2019)
and is shown in Fig. 4. The material is common glass with a nominal thickness
of 10 mm. The tank is filled with water up to 20 cm and the power output from
the heater is set to 12 W. For analysis and comparison against RTD
measurements, six different locations were chosen as illustrated in Fig. 4.
Some of these points are right above the heater, while others are slightly
off. This is done to capture qualitatively different parts of the experiment;
e.g. unlike locations 4-6, the area right above the heater (locations 1-3) is
expected to be strongly dominated by the rising plume (see Fig. 7).

Figure 4: Sketch of the water tank used in the convective flow experiments
showing geometrical details and location of all six locations sampled (in
red). All dimensions in mm.

As we show later, a short-duration experiment is not appropriate to compare
interferometric measurements against those from RTD sensors due to the
latter’s slow response time. For this reason, experiments were conducted over
a 480 s time span. Data collection began right after switching on the heater,
which was then switched off at $t=80$ s, and data recording continued for
another 400 s. Therefore, while instantaneous temperature is not expected to
be the same for both methods (given their different response times), the final
recorded temperature should in theory match (by design, we take the initial
temperature to be the same in both methods).

###  4.2 Results

Figure 5: Line-of-sight-average temperature change relative to initial
temperature, $T-T_{0}$, for the convective flow experiment. Vertical dashed
line represents the time at which the heater was turned off. RTD data (black
line) denotes the average of all 4 RTD sensors. Insets in a) and b) are
further analysed in Section 4.2.2 and Fig. 8.

####  4.2.1 Comparison against RTD sensors

Comparison between the RTDs and interferometry in six locations within the
tank for the convective flow experiment is illustrated in Fig. 5. Overall,
both methods agree well (especially if the RTDs accuracy is taken into
account, which is of $\pm$0.04 K). However, it is worth analysing several
aspects in detail. First of all, all six sampled locations exhibit a high
degree of agreement between both techniques at the end of the experiment. This
is important because the proposed interferometric method works (i.e. yields
instantaneous temperature) by accumulating temperature changes relative to the
previous instant in time (piezo stroke), such that if major errors occurred,
these would be evident from the final temperature measurement. We are thus
confident that such error accumulation does not take place. Secondly, in all
locations a lag between both methods is evident, with temperature peaks being
smoothed out and recorded later by the RTDs (this is particularly clear in
location 4). This is an expected behaviour given the well-known slow response
time of RTDs. We can further verify this by applying a simple moving average
function to the original interferometric measurements (in this case, we report
at time $t$ the average temperature over the previous $30$ s), which closely
replicates the RTD data, as illustrated in Fig. 6 (evidently, this agreement
could improve if a different moving average filter, better mimicking the RTDs
behaviour, were employed). Finally, notice that rapid temperature variations
are not captured by the RTDs, which is especially evident in location 1, where
interferometry reports a sharp rise of about 0.23 K within 1.0 s, followed by
a decrease in temperature. This phenomenon, also observed in location 2, is
presumably caused by the rising plume, and we further analyse it in Section
4.2.2. Temperature fluctuations are also present in locations 3 and 6, which
we attribute to the plume’s transition to turbulent flow as it approaches the
free surface. Locations 4 and 5 do not exhibit the aforementioned discrepancy
between techniques due to the the primary mode of heat transfer in these
locations, which is conduction rather than convection, leading to the absence
of the sharp temperature increases observed in other locations.

Figure 6: Measurements from RTDs and interferometry in location 4 after
applying a simple moving average function to the interferometry results. Here,
the moving-mean temperature reported (red dashed line) at time $t$ represents
the average of the previous $30$ s.

####  4.2.2 Comparison against PIV and LES

The temperature variations measured in locations 1 and 2 at around $t=30$ s
are particularly intriguing, as interferometry not only reveals a rapid
increase in temperature due to the rising plume, but also uncovers three peaks
within it (see insets in Figs. 5a and 5b). To investigate further the origin
of these peaks, we analyse the setup using Particle Image Velocimetry (PIV)
and Large Eddy Simulations (LES), focusing our attention on a cross-section
normal to the heating rod passing through locations (beam paths) 1, 2 and 3,
as illustrated in Fig. 7 (see Appendices A and B for technical details on the
PIV and LES settings, respectively).

Figure 7: a) Illustration of a Large Eddy Simulation (temperature field),
indicating the location of the cross-section of interest (purple vertical
line), which has also been analysed by means of Particle Image Velocimetry
(see Fig. 8). b) Non-uniform temperature distribution observed in the heating
rod using an infrared (IR) camera. Images provided for qualitative analysis
only (hence, no colour bar).

Fig. 8 shows velocity magnitude and temperature contour plots obtained from
PIV and LES in the plane of interest. There are several reasons to expect
quantitative discrepancies between PIV and LES. For example, LES assumes a
uniform heat flux distribution from the heating rod; however, we observed a
marked non-uniformity in temperature employing an infrared camera in the lab
(see Fig. 7b). Nevertheless, a qualitative interpretation of these results is
sufficient to substantiate our forthcoming arguments.

Comparison of velocity magnitude obtained by LES and PIV shows a good
qualitative agreement, particularly regarding the time of arrival of the
plume’s front to location 1. In the LES temperature plot, this hot water front
is very clear (see Fig. 8a3), such that we can reliably associate it with the
first peak detected via interferometry (see Fig. 5a). At this point, it is
worth recalling that our method captures the average temperature change along
the laser path in the water. Therefore, as the plume progresses upwards, the
laser at location 1 travels through the thin central ‘neck’ and lateral ‘arms’
of the plume (Fig. 8b3), detecting a decrease of temperature, until the
counter-rotating vortices at the edge of the arms reach this location (Fig.
8c3), leading to the detection of a second peak. To explain then the third
peak in Fig. 5a (inset), we notice the asymmetry in the rising plume observed
via PIV (Figs. 8b1 and 8c1), which is to be expected given the system’s
sensitivity to initial conditions (Narayan et al, 2017), as well as
imperfections, such as the heat flux non-uniformity discussed above. Because
of this asymmetry, a third peak is detected, corresponding (in this case) to
the left-arm eddy, which reaches location 1 with a time delay relative to the
right-arm eddy (compare positions of dotted white circles in Figs. 8c1 and
8c2). Additional support for this interpretation comes from the fact that
multiple repetitions of this experiment consistently revealed the presence of
the aforementioned three peaks at locations 1 and 2 almost every time, but,
crucially, not always. Sometimes only two peaks would be detected, which can
be attributed to conditions when the plume’s asymmetry was negligible. This
method’s sensitivity thus allows us to capture the signature of fine fluid
dynamical patterns in a fully non-intrusive manner.

Figure 8: Velocity magnitude obtained via Particle Image Velocimetry (left
panels) and Large Eddy Simulations (centre panels), as well as temperature
contour plots from LES (right panels), at three different times after start of
the experiment: a) 25 s, b) 32 s, and c) 40 s. The green horizontal lines
represent the laser path at locations 1, 2 and 3. Note, the grey semicircle in
the PIV images represents the ‘mask’ used for image processing, not the cross-
sectional area of the heating rod (which is shown in white in the LES images).
From PIV results, an asymmetry in the rising plume is evident.

###  4.3 Thermal gradient effect on interference visibility

During the convective flow experiment, the interference visibility $\gamma$
(see Fig. 1b) decreases substantially when the thermal plume reaches the
measurement location, as indicated by the blue circle in Fig. 9a. We claim
that this effect is caused by the large vertical temperature gradient across
the laser beam induced by the thermal plume’s front, causing significant
refraction of the probe beam and subsequent misalignment with the reference
beam at the photodetector. This hypothesis is confirmed by recording the
beam’s position using a camera, the results of which are shown in Fig. 9b.
Note that despite the reduced visibility periods observed in Fig. 9a, we are
still able to track the same interference peak and measure temperature change
reliably in the convective flow setting. In other words, the visibility drop
does not lead to accumulation of measurement errors, as discussed in Section
4.2.1. However, in other experimental settings (e.g. a stronger plume
presenting a much larger temperature gradient), visibility reduction can lead
to temporary data loss. This in turn makes visibility reduction due to thermal
gradients a useful indicator of measurement performance. If the thermal
gradient is excessive, our method (unlike HI) will provide no data, instead of
yielding a distorted measurement. This feature distinguishes the present
approach (SI) from holographic interferometry, where refraction (due to the
thermal gradients in the fluid) can lead to phase reconstruction errors (Lira
and Vest, 1987) (see also Section 1), and renders SI especially suitable for
probing rapid temperature changes in three-dimensional flows like the one
herein considered.

Figure 9: a) Interference signal detected by the PD at location 1, indicating
the drop in visibility $\gamma$ (blue dashed circle) caused by the arrival of
the plume’s front (see inset in Fig. 5a). b) We use a camera to record the
temporary change in the probe beam’s position due to refraction, shown by the
red (original/permanent position) and yellow (temporary/displaced position)
circles; compare with inset in Fig. 1c.

##  5 Discussion

We now provide a summary of the main strengths and limitations of the proposed
technique, as well as an outlook of potential future applications and
improvements.

While using relatively simple and economic instrumentation, standard
interferometry has the potential to yield highly precise temperature
measurements in a fully non-intrusive fashion. Moreover, the robustness of
this technique, along with its small data processing requirements, render it
attractive for real-time applications. However, the temperature measurements
obtained represent the average along the beam path within the test fluid (the
line of sight), which may limit application of this method to a set of
specific problems (see below). What is more, the interferometer is very
sensitive to ambient disturbances (such as mechanical vibrations) and, for
optimal performance, may need to be realigned before the start of an
experiment (which can be time-consuming and requires specific skills from the
experimentalist).

Despite current limitations, the present technique could be employed in small-
scale experiments where high precision and non-intrusiveness are desirable.
For example, it may be used to characterise quasi-two-dimensional flows (where
the dimension traversed by the probe beam is much smaller than the other two),
yielding data that might potentially serve as benchmark for validation of
Computational Fluid Dynamics models. Moreover, in 3D flows (like that in
Section 4), several interferometers may be arranged so as to span a plane,
potentially enabling the recovery of more localised temperature data. Due to
the increase of precision with integration time (up to a certain point), this
technique can also be helpful in long experiments with slow-evolving
temperature where very high precision measurements (in the order of ${\rm\mu
K}$) might be beneficial.

In this paper, we have purposely applied this method to two problems with a
significant disparity in nature and complexity (i.e. slow varying temperature
in static water and a convective 3D flow), with the goal of better
illustrating the strengths and weaknesses of our approach. However, our
overarching aim is to promote within the community the adoption and adaptation
of standard interferometry for the study of a wider range of problems in
experimental fluid mechanics.

##  6 Conclusions

We explore some of the merits and shortcomings of standard interferometry by
applying it to two fluid problems: slow-varying temperature in static water
and natural convective flow. We detect slow and rapid line-of-sight-average
temperature changes in water at mK sensitivity in a fully non-intrusive
manner. This is enabled by the interferometer’s high sensitivity to minute
temperature-induced changes in the fluid’s refractive index. We achieve a
precision of about 1 mK per measurement with a sampling rate of 40 Hz. The
high sensitivity and fast response time of our method allow us to capture the
signature of fine fluid dynamical patterns that cannot be resolved with
traditional methods such as RTD sensors. For example, in convective flow
induced by a horizontal heating rod, we observe an asynchronous arrival of two
counter-rotating vortices at the measurement location, revealing an asymmetric
thermal plume, which we further corroborate by means of Particle Image
Velocimetry and Large Eddy Simulations.

Resolving potentially large thermal gradients induced by a rising plume is
particularly challenging for any temperature measurement technique. For
example, RTD sensors, in addition to being intrusive, are unable to capture
this phenomenon due to their relatively slow response time, whereas
Holographic Interferometry based methods cannot reliably operate in the
presence of large thermal gradients. Conversely, we are able to operate in the
presence of said steep gradients thanks to the fast response and robustness of
our approach.

Despite its current limitations (mainly stemming from the line-of-sight-
average nature of its measurements), the technique herein described may be
utilised to investigate slow and rapid fluid temperature variations in certain
types of experiments (e.g. quasi-2D flows), generating high-quality data that
might be used, for example, for validation of Computational Fluid Dynamics
models. Moreover, as illustrated here, combination of this with other
techniques (such as PIV), can enhance the analysis of a wide variety of
problems in experimental fluid dynamics.

##  Appendix A PIV settings

A green laser (532 nm) with an output power of 4.5 mW is used as the source,
and a cylindrical lens is employed to generate the light sheet necessary for
Particle Image Velocimetry. A high-speed CMOS camera (Baumer VCXU-15M) was
utilised to capture particle images at 100 frames per second. A bandpass
filter was placed in front of the camera to eliminate ambient light and ensure
that only the scattered green light from the tracer particles was captured.
Polyamide particles are used as tracer particles, with a diameter of around 55
$\mu$m and a density of 1.016 g/cm3. Our field of view (FOV) is a square
region that extends from the centerline of the heating rod to the uppermost
water level, with a width equal to the inner width of the tank. Preliminary
tests indicated that the maximum velocity magnitude was approximately 10 mm/s.
As suggested by Lu and Sick (2013), particle displacements between subsequent
frames should not exceed 1/4 of the interrogation window size (8 pixels in our
configuration) to reduce the number of pairing loss (loss of particle images
within the interrogation window between subsequent frames). The time required
for a particle to move 8 pixels is 0.07 s for the maximum expected velocity,
which is much larger than the time between adjacent frames for our
configuration (0.01 s), indicating that our chosen frame rate is suitable.

Image processing was performed using PIVlab (Thielicke and Stamhuis, 2014). A
multi-pass interrogation algorithm with window deformation was implemented to
improve the calculation accuracy of the velocity field. Three passes were used
in total, which are 128$\times$128, 64$\times$64 and 32$\times$32 pixels with
a 50% overlap. To investigate the influence of the number of passes, we
compare the velocity magnitude in a selected point (along location 1 and just
above the heating rod) obtained by also using 2 passes only (128$\times$128
and 64$\times$64 pixels) and 4 passes (128$\times$128, 64$\times$64,
32$\times$32 and 16$\times$16 pixels). Results are shown in Fig. 10. We
observe that 2-passes produces lower velocities compared to the other two
settings, while 4-passes introduces the most noise. Thus, in this paper, we
report results using the 3-pass interrogation algorithm. The interrogation
area size for each pass are 11.5 mm, 5.8 mm and 2.9 mm.

Figure 10: Velocity magnitude at a point along location 1 and just above the
heater, obtained by using three different pass settings when processing PIV
images.

##  Appendix B LES settings

The computational domain in the simulation is set to match our experimental
settings, as detailed in Section 4.1. The Navier-Stokes equations are solved
using a Large Eddy Simulation (LES) approach with Boussinesq approximation,
and the subgrid-scale (SGS) turbulence was modelled using the dynamic
Smagorinsky model. The mesh grid has a total number of $40\times 10^{6}$
hexagonal cells. The minimum grid size is 0.03 mm, and the non-dimensional
wall distance $y^{+}$ on the heating rod is 0.023. For the boundary
conditions, we assign no-slip conditions to the side and bottom walls, while
the top wall is designated as a free-slip boundary, considering the fact that
any free surface level fluctuations that may occur will be negligible in
magnitude compared with the water depth (McSherry et al, 2017). We set the
heater’s boundary condition as a no-slip wall and define a constant heat flux
of 1504 Wm-2, in line with the experiment’s power output of 12 W. The initial
temperature is set to 22∘C everywhere. Simulations are carried out in ANSYS
Fluent 21.2. The numerical algorithm and solver settings were adapted from Ma
and He (2021), who also analysed a similar case study using LES. A pressure-
based solver was employed, and SIMPLEC scheme was used for the pressure-
velocity coupling scheme. Discretisation of the gradients and pressure
equation was conducted by least squares cell based and body force weighted
scheme, respectively. For the momentum and energy equations, boundary central
differencing and second order upwind schemes are adopted, respectively. Time
marching is achieved by a bounded second order implicit scheme. The
simulations were run on the University of Southampton’s Iridis 5
supercomputer.

To determine the integration time step, Fig. 11 illustrates the mean velocity
magnitude along location 1 (that is, averaged over the beam path at location
1) for three different values of the time step. The decrease in time step from
0.01 to 0.005 s yields a marginal mean difference of only 0.24%. In contrast,
a reduction in the time step from 0.05 to 0.01 s results in a larger mean
difference of 1.79%. Therefore, we conclude that the benefits of the
additional computational time and resources required for smaller time steps
beyond 0.01 s are not justified. Hence, we determine that a time step of 0.01
s provides sufficient precision for our needs and was consequently selected
for our simulations.

Figure 11: LES-obtained mean velocity magnitude along location 1 (i.e.
averaged over the beam path at location 1) for three different time steps.

The quality of the LES computational grid was judged based on a Reynolds-
Averaged Navier-Stokes (RANS) simulation of the same case. The turbulent
integral length scale, $l_{1}=k^{3/2}/\varepsilon$ (where $k$ and
$\varepsilon$ are the turbulent kinetic energy and energy dissipation rate,
respectively), was derived from the RANS simulation and served as a benchmark
for assessing the quality of the mesh. The LES results are deemed satisfactory
if more than 80% of the turbulent kinetic energy is resolved. This goal can be
achieved by ensuring the filter width, $\Delta$, is smaller than $l_{1}/12$,
as suggested by Pope (2000), where $\Delta$ is computed according to the
volume of the computational cell, $V$, using $\Delta=V^{1/3}$. Figure 12 shows
the ratio $l_{1}/\Delta$ computed for the grid used in our LES model. The
contour range is confined between 0 and 12, such that the grid size in the
white area complies with the $l_{1}/\Delta>12$ rule, and therefore, the core
region of the domain possesses an acceptable resolution. The non-dimensional
wall distance, or $y^{+}$, on the heater and the walls is 0.023, which is
sufficiently low to ensure a well-resolved LES near the wall.

Figure 12: Contour of the ratio between the RANS integral length scale $l_{1}$
and the filter width $\Delta$ at 40 s after start of the experiment; the white
area complies with the $l_{1}/\Delta>12$ rule. Same cross-section as in Fig.
8, but spanning the whole height of the tank.

## Declarations

Ethical Approval not applicable

Competing interests The authors have no competing interests to declare that
are relevant to the content of this article.

Authors’ contributions SM and JAZ conceived the research idea. XG performed
the experimental work and LES, and processed all data. All authors discussed
the results and contributed to the final manuscript.

Funding Partial funding was received from the University of Southampton. XG
was partly funded by the University of Southampton through a Presidential
Scholarship.

Availability of data and materials Data available from the authors upon
request

Acknowledgements The authors acknowledge the use of the IRIDIS High
Performance Computing Facility, and associated support services at the
University of Southampton, in the completion of this work. The authors also
thank two anonymous reviewers, whose thorough feedback has helped improve the
overall quality of this paper.

## References

  *   * Abram et al (2018) Abram C, Fond B, Beyrau F (2018) Temperature measurement techniques for gas and liquid flows using thermographic phosphor tracer particles. Progress in energy and combustion science 64:93–156. https://doi.org/10.1016/j.pecs.2017.09.001
  * Allan (1966) Allan D (1966) Statistics of atomic frequency standards. Proceedings of the IEEE 54(2):221–230. 10.1109/PROC.1966.4634
  * Ashjaee and Yousefi (2007) Ashjaee M, Yousefi T (2007) Experimental study of free convection heat transfer from horizontal isothermal cylinders arranged in vertical and inclined arrays. Heat transfer engineering 28(5):460–471. https://doi.org/10.1080/01457630601165822
  * Ashjaee et al (2007) Ashjaee M, Eshtiaghi A, Yaghoubi M, et al (2007) Experimental investigation on free convection from a horizontal cylinder beneath an adiabatic ceiling. Experimental Thermal and Fluid Science 32(2):614–623. https://doi.org/10.1016/j.expthermflusci.2007.07.004
  * Banks et al (2019) Banks D, Robles V, Zhang B, et al (2019) Planar laser induced fluorescence for temperature measurement of optical thermocavitation. Experimental Thermal and Fluid Science 103:385–393. https://doi.org/10.1016/j.expthermflusci.2019.01.030
  * Childs et al (2000) Childs PR, Greenwood J, Long C (2000) Review of temperature measurement. Review of scientific instruments 71(8):2959–2978. https://doi.org/10.1063/1.1305516
  * Ciddor (1996) Ciddor PE (1996) Refractive index of air: new equations for the visible and near infrared. Appl Opt 35(9):1566–1573. 10.1364/AO.35.001566, URL https://opg.optica.org/ao/abstract.cfm?URI=ao-35-9-1566
  * dataTaker (2013) dataTaker (2013) DT80 Range User’s Manual. Thermo Fisher Scientific, Australia
  * Drusová et al (2021) Drusová S, Bakx W, Doornenbal PJ, et al (2021) Comparison of three types of fiber optic sensors for temperature monitoring in a groundwater flow simulator. Sensors and Actuators A: Physical 331:112,682. https://doi.org/10.1016/j.sna.2021.112682
  * Fernández-Prini and Dooley (1997) Fernández-Prini R, Dooley R (1997) Release on the refractive index of ordinary water substance as a function of wavelength, temperature and pressure. International Association for the Properties of Water and Steam pp 1–7
  * Fernandez-Valdivielso et al (2003) Fernandez-Valdivielso C, Matias IR, Arregui FJ, et al (2003) Thermochromic-effect-based temperature optical fiber sensor for underwater applications. Optical Engineering 42(3):656–661. https://doi.org/10.1117/1.1541620
  * Freise et al (2009) Freise A, Chelkowski S, Hild S, et al (2009) Triple michelson interferometer for a third-generation gravitational wave detector. Classical and Quantum Gravity 26(8):085,012
  * Goumopoulos (2018) Goumopoulos C (2018) A high precision, wireless temperature measurement system for pervasive computing applications. Sensors 18(10):5. https://doi.org/10.3390/s18103445
  * Grafsrønningen and Jensen (2012) Grafsrønningen S, Jensen A (2012) Simultaneous piv/lif measurements of a transitional buoyant plume above a horizontal cylinder. International Journal of Heat and Mass Transfer 55(15-16):4195–4206
  * Guerrero-Mendez et al (2016) Guerrero-Mendez C, Anaya TS, Araiza-Esquivel M, et al (2016) Real-time measurement of the average temperature profiles in liquid cooling using digital holographic interferometry. Optical Engineering 55(12):121,730. 10.1117/1.OE.55.12.121730
  * Herraez and Belda (2002) Herraez J, Belda R (2002) A study of free convection in air around horizontal cylinders of different diameters based on holographic interferometry. temperature field equations and heat transfer coefficients. International journal of thermal sciences 41(3):261–267. https://doi.org/10.1016/S1290-0729(01)01314-X
  * Ito et al (2020) Ito Y, Veysset D, Kooi SE, et al (2020) Interferometric and fluorescence analysis of shock wave effects on cell membrane. Communications physics 3(1):124
  * Kashanj and Nobes (2023) Kashanj S, Nobes DS (2023) Application of 4d two-colour lif to explore the temperature field of laterally confined turbulent rayleigh–bénard convection. Experiments in Fluids 64(3):46
  * Kim et al (2018) Kim H, Byun H, Song Y, et al (2018) Multi-channel fiber-optic temperature sensor system using an optical time- domain reflectometer. Results in Physics 11:743–748. https://doi.org/10.1016/j.rinp.2018.10.024
  * Laffont et al (2018) Laffont G, Cotillard R, Roussel N, et al (2018) Temperature resistant fiber bragg gratings for on-line and structural health monitoring of the next-generation of nuclear reactors. Sensors 18(6):1791. https://doi.org/10.3390/s18061791
  * Lira and Vest (1987) Lira IH, Vest CM (1987) Refraction correction in holographic interferometry and tomography of transparent objects. Appl Opt 26(18):3919–3928. 10.1364/AO.26.003919, URL https://opg.optica.org/ao/abstract.cfm?URI=ao-26-18-3919
  * Lu and Sick (2013) Lu L, Sick V (2013) High-speed particle image velocimetry near surfaces. JoVE (Journal of Visualized Experiments) 76(50559). 10.3791/50559
  * Ma and He (2021) Ma H, He L (2021) Large eddy simulation of natural convection heat transfer and fluid flow around a horizontal cylinder. International Journal of Thermal Sciences 162:106,789. https://doi.org/10.1016/j.ijthermalsci.2020.106789
  * Mahdieh and Nazari (2013) Mahdieh M, Nazari T (2013) Measurement of impurity and temperature variations in water by interferometry technique. Optik 124(20):4393–4396
  * Mandel and Wolf (1995) Mandel L, Wolf E (1995) Optical Coherence and Quantum Optics. Cambridge University Press, United States of America, 10.1017/CBO9781139644105
  * McSherry et al (2017) McSherry RJ, Chua KV, Stoesser T (2017) Large eddy simulation of free-surface flows. Journal of Hydrodynamics, Ser B 29(1):1–12. https://doi.org/10.1016/S1001-6058(16)60712-6
  * Moldover et al (2014) Moldover MR, Gavioso RM, Mehl JB, et al (2014) Acoustic gas thermometry. Metrologia 51(1):R1. 10.1088/0026-1394/51/1/R1
  * Narayan et al (2017) Narayan S, Singh AK, Srivastava A (2017) Interferometric study of natural convection heat transfer phenomena around array of heated cylinders. International Journal of Heat and Mass Transfer 109:278–292. https://doi.org/10.1016/j.ijheatmasstransfer.2017.01.106
  * Park et al (2019) Park H, Park J, Jung SY (2019) Measurements of velocity and temperature fields in natural convective flow. International Journal of Heat and Mass Transfer 139:292–302. https://doi.org/10.1016/j.ijheatmasstransfer.2019.05.022
  * Pope (2000) Pope SB (2000) Turbulent flows. Cambridge university press, New York
  * Rizzolo et al (2017) Rizzolo S, Périsse J, Boukenter A, et al (2017) Real time monitoring of water level and temperature in storage fuel pools through optical fibre sensors. Scientific Reports 7(8766):1–10. https://doi.org/10.1038/s41598-017-08853-7
  * Roriz et al (2020) Roriz P, Silva S, Frazao O, et al (2020) Optical fiber temperature sensors and their biomedical applications. Sensors 20(7):2113. 10.3390/s20072113
  * Ross et al (2001) Ross D, Gaitan M, Locascio LE (2001) Temperature measurement in microfluidic systems using a temperature-dependent fluorescent dye. Analytical chemistry 73(17):4117–4123. https://doi.org/10.1021/ac010370l
  * Sakakibara and Adrian (2004) Sakakibara J, Adrian R (2004) Measurement of temperature field of a rayleigh-bénard convection using two-color laser-induced fluorescence. Experiments in fluids 37(3):331–340. 10.1007/s00348-004-0821-3
  * Smith and Dobson (1989) Smith LM, Dobson CC (1989) Absolute displacement measurements using modulation of the spectrum of white light in a michelson interferometer. Applied Optics 28(16):3339–3342
  * Tango et al (1968) Tango WJ, Link JK, Zare RN (1968) Spectroscopy of k2 using laser-induced fluorescence. The Journal of Chemical Physics 49(10):4264–4268. https://doi.org/10.1063/1.1669869
  * Taylor and Lai (2021) Taylor AT, Lai EP (2021) Current state of laser-induced fluorescence spectroscopy for designing biochemical sensors. Chemosensors 9(10):275. https://doi.org/10.3390/chemosensors9100275
  * Thielicke and Stamhuis (2014) Thielicke W, Stamhuis E (2014) Pivlab–towards user-friendly, affordable and accurate digital particle image velocimetry in matlab. Journal of open research software 2(1). 10.5334/jors.bl
  * Tomita et al (2006) Tomita E, Kawahara N, Toda Y (2006) A new sensor for temperature measurement of water by laser interferometry technique. Proceedings of 13th International Symposium on Applications of Laser Techniques to Fluid Mechanics pp 26–29
  * Wang et al (2022) Wang C, Xie X, Zhang H, et al (2022) Single-element real-time interferometric system for measuring dynamic temperature field of liquid medium. AIP Advances 12(4):045,010. https://doi.org/10.1063/5.0087196
  * Wang et al (2018) Wang Y, Zou F, Cegla FB (2018) Acoustic waveguides: An attractive alternative for accurate and robust contact thermometry. Sensors and Actuators A: Physical 270:84–88. https://doi.org/10.1016/j.sna.2017.12.049
  * Wilkie and Fisher (1963) Wilkie D, Fisher SA (1963) Measurement of temperature by mach-zehnder interferometry. Proceedings of the Institution of Mechanical Engineers 178(1):461–470. 10.1177/002034836317800166
  * Wu et al (2013) Wu B, Zhao J, Wang J, et al (2013) Visual investigation on the heat dissipation process of a heat sink by using digital holographic interferometry. Journal of Applied Physics 114(19):193,103. https://doi.org/10.1063/1.4832479
  * Zhang et al (2019) Zhang H, Xu X, Zhao H, et al (2019) Water temperature and salinity measurement using frequency comb. Applied Sciences 9(23):5043


[1]Konstantin Batygin

[1]Division of Geological and Planetary Sciences, California Institute of
Technology, 1200 E. California Blvd., Pasadena, 91125, CA, USA

2]Department of Space Studies, Southwest Research Institute, 1050 Walnut St.,
Suite 300, Boulder, 80302, CO, USA

# Self-Gravitational Dynamics Within the Inner Oort Cloud

<EMAIL_ADDRESS>David Nesvorný<EMAIL_ADDRESS>*  [

###### Abstract

The formation of the Inner Oort Cloud (IOC) — a vast halo of icy bodies
residing far beyond Neptune’s orbit — is an expected outcome of the solar
system’s primordial evolution within a stellar cluster. Recent models have
shown that the process of early planetesimal capture within the trans-
Neptunian region may have been sufficiently high for the cumulative mass of
the Cloud to approach several Earth masses. In light of this, here we examine
the dynamical evolution of the IOC, driven by its own self-gravity. We show
that the collective gravitational potential of the IOC is adequately
approximated by the Miyamoto-Nagai model, and use a semi-analytic framework to
demonstrate that the resulting secular oscillations are akin to the von
Zeipel-Lidov-Kozai resonance. We verify our results with direct $N$-body
calculations, and examine the effects of IOC self-gravity on the long-term
behavior of the solar system’s minor bodies using a detailed simulation.
Cumulatively, we find that while the modulation of perihelion distances and
inclinations can occur within an observationally relevant range, the
associated timescales vastly surpass the age of the sun, indicating that the
influence of IOC self-gravity on the architecture of the solar system is
negligible.

###### keywords:

Trans-Neptunian objects, Orbital dynamics, Perturbation theory, N-body
simulations

##  1 Introduction

The prevailing consensus holds that our solar system originated within a
stellar cluster – a conclusion bolstered by multiple lines of evidence that
extend beyond the statistical observation that most stars do not form in
isolation [2]. A key piece of supporting evidence for this notion is the
prevalence of decay products from short-lived radionuclides (notably 26Al),
within primitive bodies of the solar system. The widespread presence of such
elements is commonly viewed as a strong indicator that the proto-solar cloud
emerged in a relatively densely populated stellar association [1, 3].

The presence of a stellar cluster around the nascent solar system had a number
of important implications, influencing not just the isotopic makeup of solar
system bodies but also their dynamics. In fact, the formation of the Inner
Oort Cloud (IOC) — a population of long-period minor bodies that do not
strongly interact with the giant planets or the Galaxy — necessitates
extrinsic perturbations [15, 12]. Particularly, the assembly of the IOC is
envisioned to have taken place through the following sequence of events.
First, as giant planets accreted and migrated within the proto-solar nebula,
they scattered primordial planetesimals onto highly eccentric orbits. Over
time, these orbits chaotically diffused outward, where the combined effects of
stellar flybys and the cluster’s tidal forces began to alter their
trajectories. For a subset of these objects, external perturbation led to
orbital circularization, dynamically detaching them from the giant planets.
Eventually, as the sun exited the cluster, these orbits became fossilized in
the trans-Neptunian region of the solar system, generating a population of
objects with heliocentric distances that span hundreds to thousands of AU [8,
9, 17].

In light of the expected existence of a vast population of IOC objects,
several authors have explored the potential role of self-gravitational effects
in shaping the large-scale architecture of the trans-Neptunian region. For
instance, Madigan and McCourt [18], Zderic and Madigan [29] have examined the
emergence of the so-called inclination instability, which occurs due to the
self-interaction of a multitude of highly eccentric orbits initially confined
to a common plane. Their findings suggest that non-trivial dynamical behavior
could arise on a timescale far exceeding the orbital period, provided the
combined mass of the orbiting bodies exceeds $\sim 10$ Earth masses. However,
the recent study of Das and Batygin [11] showed that the inclusion of self-
consistent modeling of the giant planets – particularly Neptune’s scattering
effects – fully suppresses the inclination instability.

In a separate theoretical approach, Sefilian and Touma [25] investigated the
collective effects of a massive, lopsided trans-Neptunian disk of
planetesimals extending to hundreds of AU. The authors proposed that such a
disk could account for the observed clustering of longitudes of perihelion
among long-period, dynamically detached TNOs [4]. Yet, the physical mechanisms
responsible for the formation and maintenance of such an asymmetric structure
in the outer solar system remain elusive.

Recently, Nesvorný et al. [22] presented a detailed study of the IOC’s
formation, which incorporates sophisticated modeling of the early migration of
the giant planets, along with a self-consistent treatment of the cluster’s
effects. Through subsequent $N-$body simulations, they showed that forming the
solar system in a dense stellar environment produces perihelion and semimajor
axis distributions among long-period trans-Neptunian objects that closely
mirrors the observational data. Intriguingly, these simulations also indicate
that the total mass captured in the IOC can exceed the mass of the Earth by as
much as a factor of three, a revelation that underscores the efficiency with
which cluster-driven dynamics could have shaped the outer reaches of the solar
system.

In a distinct effort, Saillenfest et al. [24] considered the secular evolution
of the IOC under the combined perturbations from the giant planets and the
Galactic tide. Remarkably, they demonstrated that much of the parameter space
associated with the IOC is perforated by chaos, insinuating that the long-term
evolution of the IOC is considerably less inert than previously thought
(although the rate of chaotic diffusion is exceptionally slow).

Inspired by the combined works of Saillenfest et al. [24] and Nesvorný et al.
[22], here we explore previously overlooked dynamics within the IOC that are
driven by its own self-gravity. To this end, in Section 2, we introduce the
Miyamoto-Nagai potential-density pair as a simple model for the collective
gravity of the IOC. Then, in Section 3, we demonstrate from semi-analytic
grounds that the emergent secular interactions are akin to von Zeipel-Lidov-
Kozai cycles. Subsequently, we show that while self-gravity can, in principle,
drive dynamically detached orbits to join the Scattered Disk, the timescale to
achieve this greatly exceeds the lifetime of the solar system, unless the mass
of the IOC is assumed to be unjustifiably large.

Figure 1: Visualization of the particle distribution within the IOC as derived
from the $t=300\,$Myr snapshot of the Nesvorný et al. [22] cluster2
simulation. Only objects with semi-major axes in excess of 250 AU and
perihelia greater than the scattering threshold given by equation (1) are
shown. The left panel displays a top-down view (x-y projection) of the
particle distribution, highlighting the approximate axisymmetry of the
generated cloud. The middle panel presents the side view (x-z projection) of
the same distribution, effectively illustrating the cloud’s dual composition:
a central flattened disk is embedded within a more extended, quasi-spherical
halo. The right panel depicts the distribution of orbital elements of the IOC:
the inclination (top) and eccentricity (bottom) are shown as functions of the
semi-major axis, with the bounding (critical) perihelion distance (equation 1)
marked with a red curve. Additionally, the semi-major axis – eccentricity
panel depicts contours of constant perihelion distance, with the values of $q$
labeled accordingly.

##  2 Miyamoto-Nagai Model of the Inner Oort Cloud

As a first step in our analysis, let us consider the structure of the IOC. To
construct the physical model, we derive the initial conditions from the
cluster2 simulation presented in Nesvorný et al. [22], adopting the
distribution of IOC objects at the 300 Myr mark. This simulation provides a
nuanced model of the solar system’s primordial evolution, incorporating
detailed and well-tested prescriptions of giant planet migration and the
dynamical influence of the Sun’s birth cluster. Specifically, the cluster2
simulation models stellar encounters self-consistently, assuming a Plummer
model for the Sun’s birth cluster characterized by a mass of 1200$M_{\odot}$
and a Plummer radius of 0.35 pc. This environment approximates conditions akin
to the Orion Nebular Cluster (see e.g, Batygin et al. 6), and the product of
the local stellar density and the sun’s residence time within the cluster –
taken to be $\eta\,\tau\sim 10^{4}$Myr pc-3 – is tuned to reproduce the semi-
major axis distribution of high-perihelion scattered disk objects.

The choice of the 300 Myr snapshot is informed by the timeline of early solar
system evolution and cluster dynamics. Despite being early in the the solar
system’s lifetime, by this epoch, the giant planets’ migration — modeled in
Nesvorný et al. [22] with a characteristic migration timescale of
$\tau_{\rm{mig}}=10$ Myr has largely concluded. Additionally, this point marks
a phase where the solar system’s birth cluster had already dissipated, and the
formation of the IOC is complete.

Objects on unstable111See e.g., Batygin et al. [5], Huang and Gladman [16] for
long-term integrations of unstable large-semi-major axis TNOs. (i.e., Neptune-
crossing and scattering) orbits will diffuse out of the solar system on a sub-
Gyr timescale, and will therefore not contribute to the cumulative
gravitational potential of the (long-term stable) “inert” inner Oort Cloud.
Indeed, in the limit where $q\rightarrow a_{\rm{N}}$, the semi-major axis
diffusion approaches $\mathcal{D}_{a}\sim
8/(5\pi)\,(m_{\rm{N}}/M_{\odot})\,\sqrt{\mathcal{G}\,M_{\odot}\,a_{\rm{N}}}$
[7, 14]. In an effort to filter these particles out of our analysis and
isolate the relevant population within the broader trans-Neptunian region, we
exclude particles within the active Neptune-scattering zone, as well as those
with (comparatively) short periods. This is achieved by removing objects with
semi-major axis smaller than $250\,$AU and perihelion distance below the
scattering threshold [7]:

$\displaystyle
q_{\rm{crit}}=\mathrm{sup}\Bigg{[}q_{0},a_{\rm{N}}\sqrt{\log\left(\frac{24^{2}}{5}\frac{m_{\rm{N}}}{M_{\odot}}\left(\frac{a}{a_{\rm{N}}}\right)^{5/2}\right)}\Bigg{]},$
(1)

where quantities with the subscript “N" refer to Neptune, and $q_{0}=40\,$AU.
It is important to note that this expression for the critical perihelion
distance is derived under the assumption of planar motion ($i\rightarrow 0$)
and small semi-major axis ratio ($a_{\rm{N}}/a\ll 1$) (Batygin et al. 7; see
also Hadden and Tremaine 14), and hence, our approach may be somewhat
restrictive. That is, chaotic diffusion of orbits is reduced at high
inclinations, suggesting that our cut-off might be more stringent than
necessary, given that it is derived in the planar approximation. However, for
the purposes of our initial exploration into the IOC’s dynamics, this
simplification provides a suitable starting point.

Figure 2: Density profile of the IOC. The left panel showcases the
azimuthally-averaged particle density, $\rho_{\rm{num}}$ derived from detailed
$N-$body simulation of Nesvorný et al. 22 (see also Fig. 1). The middle panel
shows the Miyamoto-Nagai density profile, $\rho$ (given by equation 2), with
scaling parameters $\tilde{a}=200\,$AU and $\tilde{b}/\tilde{a}=5$. In both
panels, the density is normalized to unity at $r=250$ AU and $z=0$, and a
sequence of contours is labeled, depicting $\log_{10}(\rho)$. The right panel
depicts the error associated with the analytic fit,
$\epsilon=(\rho-\rho_{\rm{num}})/\rho_{\rm{num}}$. As can be seen, the error
within the relevant regions is limited to tens of percent, which is
satisfactory for the purposes of our study.

Examination of the particle distribution in physical space (as shown in Fig 1)
reveals that the particle distribution within the IOC is, to a good
approximation, axisymmetric222To ascertain the degree of axial symmetry, we
computed the distribution of azimuthal angles of particles across cylinders
with a radial thickness of 100 AU at radial intervals of 100 AU, and applied
the Kolmogorov-Smirnov test for uniformity. This calculation yielded an
average $p$-value of $\langle p\rangle=0.45$, indicating that axial symmetry
is statistically consistent with the simulation data.. However, the
distribution is markedly not spherical: when viewed from the side, there is a
discernible disk-like structure embedded within a more diffuse halo. To model
this distribution, we adopt the Miyamoto-Nagai potential-density pair [19].

The Miyamoto-Nagai potential is a versatile model that allows for the
simulation of a broad range of astrophysical systems, with adjustable
thickness and radial distribution of mass. It is defined by two crucial
parameters: the radial scale length ($\tilde{a}$) and the scale height
($\tilde{b}$). These parameters govern the shape and density distribution of
the particle population. As $\tilde{b}$ approaches zero, the model converges
to the Kuzmin potential, representing a razor-thin disk. Conversely, when
$\tilde{a}$ approaches zero, the model reduces to a Plummer sphere. The ratio
$\tilde{b}/\tilde{a}$ is thus a key factor in determining the disk’s flatness:
a high $\tilde{b}/\tilde{a}$ ratio results in a mass distribution that is
nearly spherical, while a low $\tilde{b}/\tilde{a}$ ratio leads to a
significantly flattened distribution. This flexibility allows the Miyamoto-
Nagai potential to adeptly model various astrophysical structures by adjusting
$\tilde{a}$ and $\tilde{b}$ accordingly.

Written in terms of cylindrical coordinates, the potential $\Psi$ and the
corresponding density have the functional form:

$\displaystyle\Psi$
$\displaystyle=-\frac{\mathcal{G}\,M_{\rm{IOC}}}{\sqrt{\big{(}\tilde{a}+\sqrt{\tilde{b}^{2}+z^{2}}\big{)}^{2}+r^{2}}},$
(2) $\displaystyle\rho$
$\displaystyle=\frac{\nabla^{2}\,\Psi}{4\,\pi\,\mathcal{G}}=\frac{\tilde{b}^{2}M_{\rm{IOC}}\,\big{(}\big{(}\tilde{a}+\sqrt{\tilde{b}^{2}+z^{2}}\big{)}^{2}\big{(}\tilde{a}+3\sqrt{b^{2}+z^{2}}\big{)}+\tilde{a}\,r^{2}\big{)}}{4\,\pi\,\big{(}\tilde{b}^{2}+z^{2}\big{)}^{3/2}\big{(}\big{(}\tilde{a}+\sqrt{\tilde{b}^{2}+z^{2}}\big{)}^{2}+r^{2}\big{)}^{5/2}},$

where $M_{\rm{IOC}}$ is the total mass of the cloud.

To determine the appropriate parameter values, we began by fitting the
analytic model to particle distribution shown in Fig 1, which yielded
$\tilde{a}=155$ AU and $\tilde{b}/\tilde{a}=7$. We subsequently fine-tuned
these parameters to $\tilde{a}=200$ AU and $\tilde{b}/\tilde{a}=5$. This
adjustment slightly reduced the discrepancy between the analytic model and the
numerical results within the disk region, albeit at the cost of accuracy in
the much more sparsely populated regions at very large heliocentric distances.
Figure 2 depicts a comparison between the azimuthally-averaged particle
distribution obtained from the detailed $N-$body simulation and that given by
equation 2.

##  3 Secular Dynamics Within the IOC

With the potential of the IOC now defined, we are in a position to write down
the model Hamiltonian that will govern the dynamics within this system. The
chosen potential, characterized by its smoothness and axisymmetry, implies
that the IOC will predominantly drive secular evolution. This class of
dynamics corresponds to long-term orbital changes, smoothing over short-term
perturbations. It is important to acknowledge, however, that this phase-
averaged description, only provides an adequate approximation for bodies
outside the scattered disk, and more generally, for objects not trapped in
high-order MMRs with Neptune [20].

Accounting for secular perturbations from Jupiter, Saturn, Uranus, and Neptune
to quadrupolar order, we have:

$\displaystyle\mathcal{H}=\frac{\mathcal{G}\,\mathcal{C}\,\big{(}3\cos^{2}(i)-1\big{)}}{8\,a^{3}\,\big{(}1-e^{2}\big{)}^{3/2}}+\frac{1}{2\,\pi}\oint\Psi\,d\mathcal{M},$
(3)

where $\mathcal{M}$ is the mean anomaly of the IOC particle, and
$\mathcal{C}=\sum_{i}m_{i}\,a_{i}^{2}$ is the sum of the moments of inertia of
the planetary orbits. Although the integral in the above expression cannot be
evaluated analytically, general properties of $\mathcal{H}$ can nevertheless
be deduced from symmetry. First, $\mathcal{H}$ is, by construction,
independent of the orbital phase, meaning that the semi-major axis, $a$, is
rendered a constant of motion. Similarly, axial symmetry implies that
$\mathcal{H}$ is independent of the longitude of ascending node, $\Omega$,
meaning that the conjugate action, $\mathcal{J}$, – which corresponds to the
$\hat{z}$-component of the specific angular momentum vector, is conserved. In
other words, $\mathcal{H}=f(a,\mathcal{J},e,\omega)$ depends only on a single
angle – the argument of perihelion, $\omega$. Moreover, because the difference
between ascending and descending nodes disappears upon phase-averaging, the
dependence of $\mathcal{H}$ on $\omega$ is bound to be $\pi-$periodic.

These properties place $\mathcal{H}$ into the broad category of von Zeipel-
Lidov-Kozai type Hamiltonians, i.e., systems with only a single degree of
freedom – related to the $(e,\omega)$ variable pair – parameterized by $a$ and
$\mathcal{J}$ [27, 26]. Thus, upon an appropriate choice of the constants of
motion, the phase-space portrait of an IOC particle can be fully quantified by
projecting level curves of $\mathcal{H}$ onto the $(\omega,q)$ plane.

Figure 3: Level curves of Hamiltonian (3) projected onto the $(\omega,q)$
plane. These contours depict the secular dynamics within the IOC for particles
with semi-major axes of $a=1000\,$AU (left panel) and $a=5000\,$AU (right
panel), for a total IOC mass of $M_{\rm{IOC}}=3M_{\oplus}$. Both panels
illustrate the characteristic structure of von Zeipel-Lidov-Kozai (vZLK)
resonance. A series of contours are labeled according to
$\log_{10}(-\mathcal{H})$, evaluated in units of $M_{\odot}$, AU, and yr. The
dotted yellow line corresponds to the Neptune-scattering threshold in the
coplanar limit, and the white curves depict output from an $N-$body
simulation, for equivalent values of $\mathcal{J}$.

Though the IOC as a whole covers a broad range of parameters, here, our
selection of constants of motion is motivated by the observational census of
known long-period TNOs, which primarily consists of high-eccentricity objects
situated near the plane of the solar system. Considering the perturbative
nature of the model, the maximum allowable eccentricity corresponds to a
perihelion distance equal to Neptune’s semi-major axis, beyond which orbits
would cross and fundamentally disrupt the model’s accuracy. Thus, for
definitiveness, we define the action variable to correspond to an orbit with
zero inclination and $q=a_{\rm{N}}$:

$\displaystyle\mathcal{J}=\sqrt{1-e^{2}}\,\cos(i)=\sqrt{\frac{a_{\rm{N}}\,(2\,a-a_{\rm{N}})}{a^{2}}}=\mathrm{const.}$
(4)

We reiterate, however, that even for orbits with larger perihelion distances,
the secular Hamiltonian can yield an inaccurate representation of the
dynamics, if the IOC particle becomes embedded within a dense network of mean
motion resonances that facilitate Neptune scattering [7]. Despite these
limitations, this semi-analytic approach is still effective in exploring
whether orbital evolution driven by IOC’s self-gravity, can elevate the
eccentricities of initially high-perihelion objects to the scattering
threshold, thereby cycling them between the IOC and scattered disk.

Figure 4: Time-series of perihelion distance and inclination from idealized
N-body simulations of IOC particles at $a=1000\,$AU (left panel) and
$a=5000\,$AU (right panel) for $M_{\rm{IOC}}=1M_{\oplus}$. Driven by IOC’s
self-gravity, the particles execute vZLK dynamics, where perihelion distance
and inclination oscillate in concert, preserving the vertical component of the
specific angular momentum vector. These numerically computed orbits are
projected onto the $(\omega,q)$ in Figure 3, providing a direct comparison
between the simulation results and the analytical model. While the agreement
between the two is satisfactory, these time-series also underscore the
exceedingly long timescales over which IOC-driven secular dynamics operate. As
depicted, the orbital evolution unfolds over durations that far exceed the age
of the solar system.

Figure 3 shows the contours of $\mathcal{H}$ on the $(q,\omega)$ plane for IOC
particles at $a=1000$ and $5000\,$AU, where we have set
$M_{\rm{IOC}}=3M_{\oplus}$, in agreement with the simulation results of
Nesvorný et al. [22]. In both panels, the familiar structure of the vZLK
resonance can be seen, and intriguingly, secular IOC-driven dynamics appears
to modulate the perihelia over an observationally relevant (i.e., $q\sim
100\,$AU) range. Moreover, the reduction of the phase-space portrait towards
trivial circulation for $q$ approaching $30$ AU can be understood as emanating
from the growing contribution of the first term in Hamiltonian (3). More
specifically, our analysis demonstrates that the contribution of the first
term, representing precession driven by the giant planets, to the overall
secular dynamics of the IOC, significantly varies with perihelion distance. At
$a=1000$ AU, this term accounts for less than 2% of the Hamiltonian’s value at
$q=75$ AU, increasing to 7% and 25% at $q=50$ and $30$ AU, respectively.
Meanwhile, at $a=5000$ AU, its influence is relatively subdued, constituting
approximately 10% at $q=30$ AU and diminishing to 2% at $q=50$ AU, thereby
allowing the vZLK cycle to penetrate to smaller perihelion distances.

To further quantify these self-gravitational dynamics, we supplement our
analytic calculations with direct numerical orbit propagation of test
particles. To carry out the integrations, we used the mercury6 gravitational
software package [10], utilizing the MVS algorithm [28]. To ensure consistency
with our analytical model, we include the accelerations due to the IOC via
$\mathbf{a}=-\nabla\,\Psi$ and model the quadrupolar fields of the giant
planets as an oblateness of the Sun, represented by
$J_{2}=\mathcal{C}/(2\,M_{\odot}\,\mathcal{R}^{2})$ moment, setting the
central body’s radius (which also acts as an absorbing boundary condition to
$\mathcal{R}=a_{\rm{N}}$). The initial conditions for the IOC objects, treated
as test particles, are set with a perihelion of 200 AU, an argument of
perihelion of $\omega=90\deg$, and an initial inclination determined by the
constant of motion, $\mathcal{J}$.

The time-series of the perihelion distance and inclination are shown on Figure
4. The numerically computed evolution of the orbits is also projected onto the
($\omega,e$) diagrams shown in Figure 3. Overall, these simulations show
excellent agreement with the semi-analytic model. However, the time-series
also highlight the staggeringly long timescales associated with the evolution
driven by the self-gravity of the IOC. In fact for the adopted parameters,
evolution unfolds on a timescale comparable to the age of the universe or even
longer.

Given the extraordinarily long integration period, we clarify that the intent
behind these simulations is not to present a realistic portrayal of the Solar
System’s future evolution. Instead, this calculation merely serves as an
illustrative approach to fully capture the dynamical cycle of particles
evolving under the influence of the inner Oort Cloud’s self-gravitational
potential. Indeed, owing to external perturbations such as stellar flybys, the
inner Oort cloud is expected to be depleted on a timescale much shorter than
that depicted in Figure 4. As importantly, the dynamical lifetime of the outer
solar system itself is limited to mere tens of Gyr [30]. Accordingly, these
simulations should be viewed as nothing more than a theoretical exposition of
the isolated dynamics facilitated by Hamiltonian (3).

##  4 N-body simulations

Figure 5: Time-series of perihelion distance (top panel) and inclination
(middle panel) of IOC objects, derived from the $N=10^{4}$ GENGA simulation. A
subset of trajectories that are entrained in vZLK cycles are highlighted,
emphasizing long-term evolution of such particles. Notably, particles
entrained in cycles of $\omega-$libration are rare (comprising $\sim 2\%$ of
all objects), but their existence draws a parallel between direct $N$-body
simulations and the semi-analytic model based upon the Miyamoto-Nakai
potential. The time-series of these objects are further illustrated in the
$q-\omega$ diagram (bottom panel), confirming the libration of the argument of
perihelion and showcasing the vZLK behavior expected from the semi-analytic
model.

In an effort to explore the self-gravitational dynamics of the IOC further, we
conducted a pair of additional N-body simulations. First, as a test of the
secular framework outlined in the previous section, we used the GENGA GPU-
accelerated code [13] to directly model the gravitational interactions,
without assuming the Miyamoto-Nagai potential. The principal goal of this
calculation was to assess if the vZLK-type behavior insinuated by the semi-
analytic model can be recovered within the context of a direct (“brute force")
calculation.

The initial conditions for our simulation were directly derived from the
distributions shown in the right panel of Figure 1. To expedite the dynamical
evolution, however, we implemented two significant adjustments: firstly, we
augmented the cloud’s mass and secondly, we consolidated its distribution into
a denser orbital configuration. In particular, we excluded all objects with
semi-major axes greater than 1000 AU and perihelia below $q<35\,$AU, while
elevating the total mass of the remaining distribution to 10 Earth masses.
This reduction in the range of semi-major axes (and orbital periods), coupled
with the substantial enhancement in total mass, significantly intensified the
gravitational interactions within the system. The perturbations of Jupiter,
Saturn, and Uranus were encapsulated within an effective $J_{2}$ oblateness
coefficient as before, while Neptune was explicitly included on its current
orbit. The time-step for the simulation was set at 5 years.

The simulation spanned a total integration time of 4.5 billion years,
mirroring the age of the solar system. The results, presented in Figure 5,
elucidate the long-term evolution of the particles. In particular, the top two
panels of Figure 5 display the time-series of perihelion distances and
inclinations, respectively. The bottom panel of Figure 5 shows the perihelion
distance as a function of the argument of perihelion for a subset of the
simulated particles. This diagram serves as a direct comparison to the semi-
analytic model presented earlier in Figure 3, highlighting the libration zones
characteristic of vZLK resonance. Though the depicted evolution is notably
stochastic (owing in part to the relatively coarse representation of the IOC),
the behavior of these particles aligns with the anticipated vZLK dynamics
driven by the IOC’s self-gravity.

Figure 6: Comparison of observed and simulated cumulative distributions for
orbital elements and the Tisserand parameter of Jupiter Family Comets (JFCs).
Panels (a) to (d) show the cumulative fraction $f$ of JFCs for semi-major
axis, inclination, perihelion distance, and Tisserand parameter with respect
to Jupiter ($T_{\rm{J}}$), respectively. The red dots indicate observations
(from the minor planet center), while the solid black line represents results
from our simulation incorporating self-gravity with a 3 Earth mass IOC. This
figure closely resembles the self-gravity-free simulation presented in
Nesvorný et al. [23] underscoring the negligible influence of self-gravity on
the orbital evolution of JFCs and the solar system as a whole.

Having not observed any obvious contradictions between the GENGA simulations
and the dynamical behavior expected within the context of the orbit-averaged
Miyamoto-Nagai framework, we revisited the calculations of Nesvorný et al.
[22], incorporating the potential (2) as an auxiliary effect, and retaining
the IOC mass at $3M_{\odot}$. Broadly speaking, this simulation revealed no
significant difference in any of the readily-observable trans-Neptunian
populations of minor bodies within the solar system, we additionally
interrogated the cometary population for any hints of difference between runs
with and without self-gravity. While the perturbations associated with IOC
self-gravity are confined to large heliocentric distances, short period comets
within the solar system constitute an observationally well-quantified
aggregate of minor bodies that are derived from the trans-Neptunian region
(namely, they originate at large orbital periods but scatter inwards). As a
result, perturbations to the outer solar system naturally propagate to the
cometary populations, and their distribution has previously been shown to be
affected by perturbations from dwarf planets [21] as well as the hypothetical
Planet 9 [23].

We examined the generation of Jupiter Family Comets (JFCs) in the last billion
years as a proxy for the effects of self-gravity. The results, presented in
Figure 6, allow for direct comparison with the findings of Nesvorný et al.
[23], which did not account for the cluster potential. As expected, inclusion
of IOC self-gravity did not meaningfully alter the production of JFCs compared
to the published results of Nesvorný et al. [23]. The fact that this
simulation does not violate the JFC constraint is consistent with the
aforementioned notion that the timescales for self-gravitational dynamics
within the IOC are too long to substantially impact the solar system’s
evolution. Indeed, the muted effect of self-gravity aligns with theoretical
expectations outlined in the preceding section, and underscores its ultimately
negligible role in shaping the dynamical architecture of the trans-Neptunian
region.

##  5 Discusssion

In our study, we have considered a self-consistent model for the formation of
the IOC to analyze the secular dynamics of objects within the IOC, under the
influence of their collective gravitational potential. By aligning a Miyamoto-
Nagai potential-density pair with the structure observed in simulation
outputs, we quantified the von Zeipel-Lidov-Kozai-type dynamics among
dynamically detached orbits. Our findings even suggest that the range of
perihelion modulation facilitated by the vZLK cycles can reach perihelion
distances of $q\lesssim 100$ AU, where the “detached" population of scattered
disk objects is observed. However, our simulations also established that the
timescale for these self-gravitational dynamics to manifest significantly
exceeds not only the solar system’s age but also the expected dynamical
lifetime of the giant planets (see Zink et al. 30).

Although our analysis is predicated on a specific scenario for the IOC’s
formation, and alternative model parameters could theoretically yield
quantitatively different outcomes, this is unlikely because the particular
simulation results of Nesvorný et al. [22] adopted in this study already
invoke cluster parameters and solar residence time that are close to the upper
bound of the constraints set by the dynamical structure of the cold classical
belt [6]. This means that more massive models of the IOC are likely to be
incompatible with other solar system constraints. Our results thus suggest
that any realistic parameter range would likely lead to similar conclusions.
Consequently, our research joins the ranks of other unsuccessful attempts to
attribute the dynamical architecture of the outer solar system to self-
gravitating effects.

Acknowledgements

We are thankful to Fred Adams, Alessandro Morbidelli, Cristian Beauge, and
Mike Brown for insightful discussions. We thank the anonymous referees for
providing careful and insightful reviews of the manuscript. K.B. is grateful
to Caltech, the David and Lucile Packard Foundation, and the National Science
Foundation (grant number: AST 2109276) for their generous support.

Conflict of Interest The authors declare that they have no competing interests
or other interests that might be perceived to influence the results and/or
discussion reported in this paper.

## References

  * Adams and Laughlin [2001] Adams, F. C., Laughlin, G. 2001. Constraints on the Birth Aggregate of the Solar System. Icarus 150, 151–162. doi:10.1006/icar.2000.6567
  * Adams [2010] Adams, F. C. 2010. The Birth Environment of the Solar System. Annual Review of Astronomy and Astrophysics 48, 47–85. doi:10.1146/annurev-astro-081309-130830
  * Arakawa and Kokubo [2023] Arakawa, S., Kokubo, E. 2023. Number of stars in the Sun’s birth cluster revisited. Astronomy and Astrophysics 670. doi:10.1051/0004-6361/202244743
  * Batygin and Brown [2016] Batygin, K., Brown, M. E. 2016. Evidence for a Distant Giant Planet in the Solar System. The Astronomical Journal 151. doi:10.3847/0004-6256/151/2/22
  * Batygin et al. [2019] Batygin, K., Adams, F. C., Brown, M. E., Becker, J. C. 2019. The planet nine hypothesis. Physics Reports 805, 1–53. doi:10.1016/j.physrep.2019.01.009
  * Batygin et al. [2020] Batygin, K., Adams, F. C., Batygin, Y. K., Petigura, E. A. 2020. Dynamics of Planetary Systems within Star Clusters: Aspects of the Solar System’s Early Evolution. The Astronomical Journal 159. doi:10.3847/1538-3881/ab665d
  * Batygin et al. [2021] Batygin, K., Mardling, R. A., Nesvorný, D. 2021. The Stability Boundary of the Distant Scattered Disk. The Astrophysical Journal 920. doi:10.3847/1538-4357/ac19a4
  * Brasser et al. [2006] Brasser, R., Duncan, M. J., Levison, H. F. 2006. Embedded star clusters and the formation of the Oort Cloud. Icarus 184, 59–82. doi:10.1016/j.icarus.2006.04.010
  * Brasser et al. [2012] Brasser, R., Duncan, M. J., Levison, H. F., Schwamb, M. E., Brown, M. E. 2012. Reassessing the formation of the inner Oort cloud in an embedded star cluster. Icarus 217, 1–19. doi:10.1016/j.icarus.2011.10.012
  * Chambers [1999] Chambers, J. E. 1999. A hybrid symplectic integrator that permits close encounters between massive bodies. Monthly Notices of the Royal Astronomical Society 304, 793–799. doi:10.1046/j.1365-8711.1999.02379.x
  * Das and Batygin [2023] Das, A., Batygin, K. 2023. Suppression of the inclination instability in the trans-Neptunian Solar system. Monthly Notices of the Royal Astronomical Society 523, 6103–6113. doi:10.1093/mnras/stad1687
  * Fernández [1997] Fernández, J. A. 1997. The Formation of the Oort Cloud and the Primitive Galactic Environment. Icarus 129, 106–119. doi:10.1006/icar.1997.5754
  * Grimm et al. [2022] Grimm, S. L., Stadel, J. G., Brasser, R., Meier, M. M. M., Mordasini, C. 2022. GENGA. II. GPU Planetary N-body Simulations with Non-Newtonian Forces and High Number of Particles. The Astrophysical Journal 932. doi:10.3847/1538-4357/ac6dd2
  * Hadden and Tremaine [2024] Hadden, S., Tremaine, S. 2024. Scattered disc dynamics: the mapping approach. Monthly Notices of the Royal Astronomical Society 527, 3054–3075. doi:10.1093/mnras/stad3478
  * Hills [1981] Hills, J. G. 1981. Comet showers and the steady-state infall of comets from the Oort cloud.. The Astronomical Journal 86, 1730–1740. doi:10.1086/113058
  * Huang and Gladman [2024] Huang, Y., Gladman, B. 2024. Primordial Orbital Alignment of Sednoids. The Astrophysical Journal 962. doi:10.3847/2041-8213/ad2686
  * Kaib and Quinn [2008] Kaib, N. A., Quinn, T. 2008. The formation of the Oort cloud in open cluster environments. Icarus 197, 221–238. doi:10.1016/j.icarus.2008.03.020
  * Madigan and McCourt [2016] Madigan, A.-M., McCourt, M. 2016. A new inclination instability reshapes Keplerian discs into cones: application to the outer Solar system. Monthly Notices of the Royal Astronomical Society 457, L89–L93. doi:10.1093/mnrasl/slv203
  * Miyamoto and Nagai [1975] Miyamoto, M., Nagai, R. 1975. Three-dimensional models for the distribution of mass in galaxies.. Publications of the Astronomical Society of Japan 27, 533–543. https://ui.adsabs.harvard.edu/abs/1975PASJ…27..533M
  * Morbidelli [2002] Morbidelli, A. 2002. Modern celestial mechanics : aspects of solar system dynamics. Modern celestial mechanics : aspects of solar system dynamics, by Alessandro Morbidelli. London: Taylor & Francis, 2002, ISBN 0415279399.
  * Muñoz-Gutiérrez et al. [2019] Muñoz-Gutiérrez, M. A., Peimbert, A., Pichardo, B., Lehner, M. J., Wang, S.-Y. 2019. The Contribution of Dwarf Planets to the Origin of Jupiter Family Comets. The Astronomical Journal 158. doi:10.3847/1538-3881/ab4399
  * Nesvorný et al. [2023] Nesvorný, D., Bernardinelli, P., Vokrouhlický, D., Batygin, K. 2023. Radial distribution of distant trans-Neptunian objects points to Sun’s formation in a stellar cluster. Icarus 406. doi:10.1016/j.icarus.2023.115738
  * Nesvorný et al. [2017] Nesvorný, D., Vokrouhlický, D., Dones, L., Levison, H. F., Kaib, N., Morbidelli, A. 2017. Origin and Evolution of Short-period Comets. The Astrophysical Journal 845. doi:10.3847/1538-4357/aa7cf6
  * Saillenfest et al. [2019] Saillenfest, M., Fouchard, M., Ito, T., Higuchi, A. 2019. Chaos in the inert Oort cloud. Astronomy and Astrophysics 629. doi:10.1051/0004-6361/201936298
  * Sefilian and Touma [2019] Sefilian, A. A., Touma, J. R. 2019. Shepherding in a Self-gravitating Disk of Trans-Neptunian Objects. The Astronomical Journal 157. doi:10.3847/1538-3881/aaf0fc
  * Tremaine [2023] Tremaine, S. 2023. The Hamiltonian for von Zeipel-Lidov-Kozai oscillations. Monthly Notices of the Royal Astronomical Society 522, 937–947. doi:10.1093/mnras/stad1029
  * von Zeipel [1910] von Zeipel, H. 1910. Sur l’application des séries de M. Lindstedt à l’étude du mouvement des comètes périodiques. Astronomische Nachrichten 183, 345. doi:10.1002/asna.19091832202
  * Wisdom and Holman [1991] Wisdom, J., Holman, M. 1991. Symplectic maps for the N-body problem.. The Astronomical Journal 102, 1528–1538. doi:10.1086/115978
  * Zderic and Madigan [2023] Zderic, A., Madigan, A.-M. 2023. Steeper Scattered Disks Buckle Faster. The Astrophysical Journal 948. doi:10.3847/2041-8213/accde2
  * Zink et al. [2020] Zink, J. K., Batygin, K., Adams, F. C. 2020. The Great Inequality and the Dynamical Disintegration of the Outer Solar System. The Astronomical Journal 160. doi:10.3847/1538-3881/abb8de