# On the source counts of VLBI-detected radio sources and the prospects of
all-sky surveys with current and next generation instruments

S. Rezaei,1,2 J. P. McKean3,4, A. T. Deller5 and J. F. Radcliffe6,7
1Leiden Observatory, Leiden University, Box 9513, 2300 RA Leiden, The
Netherlands
2Leiden Institute of Advanced Computer Science (LIACS), Leiden University,
P.O. Box 9512, 2300 RA Leiden, The Netherlands
3Kapteyn Astronomical Institute, University of Groningen, Postbus 800, NL-9700
AV Groningen, The Netherlands
4ASTRON, Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD
Dwingeloo, The Netherlands
5Centre for Astrophysics and Supercomputing, Swinburne University of
Technology, Mail H30, PO Box 218, Hawthorn, VIC 3122, Australia
6Department of Physics, University of Pretoria, Lynnwood Road, Hatfield,
Pretoria, 0083, South Africa
7Jodrell Bank Centre for Astrophysics, University of Manchester, Oxford Road,
Manchester M13 9PL, UK E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present an analysis of the detection fraction and the number counts of
radio sources imaged with Very Long Baseline Interferometry (VLBI) at 1.4 GHz
as part of the mJIVE–20 survey. From a sample of 24 903 radio sources
identified by FIRST, 4 965 are detected on VLBI-scales, giving an overall
detection fraction of $19.9\pm 2.9$ per cent. However, we find that the
detection fraction falls from around 50 per cent at a peak surface brightness
of 80 mJy beam-1 in FIRST to around 8 per cent at the detection limit, which
is likely dominated by the surface brightness sensitivity of the VLBI
observations, with some contribution from a change in the radio source
population. We also find that compactness at arcsec-scales is the dominant
factor in determining whether a radio source is detected with VLBI, and that
the median size of the VLBI-detected radio sources is 7.7 mas. After
correcting for the survey completeness and effective sky area, we determine
the slope of the differential number counts of VLBI-detected radio sources
with flux densities $S_{\rm 1.4~{}GHz}>1$ mJy to be $\eta_{\rm VLBI}=-1.74\pm
0.02$, which is shallower than in the cases of the FIRST parent population
($\eta_{\rm FIRST}=-1.77\pm 0.02$) and for compact radio sources selected at
higher frequencies ($\eta_{\rm JBF}=-2.06\pm 0.02$). From this, we find that
all-sky ($3\pi$ sr) surveys with the EVN and the VLBA have the potential to
detect $(7.2\pm 0.9)\times 10^{5}$ radio sources at mas-resolution, and that
the density of compact radio sources is sufficient (5.3 deg-2) for in-beam
phase referencing with multiple sources (3.9 per primary beam) in the case of
a hypothetical SKA-VLBI array.

###### keywords:

galaxies: active – radio continuum: galaxies – surveys – techniques: high
angular resolution

††pubyear: 2023††pagerange: On the source counts of VLBI-detected radio
sources and the prospects of all-sky surveys with current and next generation
instruments–A

##  1 Introduction

Very Long Baseline Interferometry (VLBI; Broten et al. 1967) is an
observational technique where the signals from individual radio telescopes are
combined coherently to produce a synthesised unfilled aperture. As the radio
telescopes can be widely separated across the Earth (and even located in
space), this technique currently provides the highest angular resolution
imaging in astronomy (with a synthesized beam size of typically 1 to 10 mas at
cm-wavelengths, and reaching 0.025 mas at sub-mm wavelengths; Event Horizon
Telescope Collaboration et al. 2019b). VLBI allows a wide range of science
goals to be realised, for example, detecting super-luminal motion in radio
jets (e.g., Cohen et al., 1977), imaging gravitational lenses (e.g., Porcas et
al., 1979), mapping the accretion disks of supermassive black holes (e.g.,
Miyoshi et al., 1995), tracing the expansion of stellar explosions (e.g.,
O’Brien et al., 2006), imaging outflows of atomic hydrogen from supermassive
black holes (e.g., Morganti et al., 2013), localising Fast Radio Bursts (e.g.,
Marcote et al., 2020), and making the first images of the shadow of a black
hole (e.g., Event Horizon Telescope Collaboration et al., 2019a).

It is over 50 years since the first fringes were produced between two
unconnected antennas separated by over 3000 km (Broten et al., 1967), yet the
unique applications of VLBI at cm-wavelengths have been restricted to studying
only $\sim 25\,000$ radio sources with very high brightness temperatures
($>10^{5}$ K) and over a very small fraction of the observable sky ($\sim 200$
deg2; e.g., Deller & Middelberg 2014; Herrera Ruiz et al. 2017; Petrov 2021).
This is significantly lower than the up to 5 million radio sources that have
been catalogued from recent all-sky surveys at arcsec-resolution (e.g., Becker
et al. 1995; Condon et al. 1998; Intema et al. 2017; Shimwell et al. 2019;
Lacy et al. 2020). This is because the effective field-of-view of a VLBI
experiment is typically only a few tens of arcsec in diameter, which is due to
historical computational limitations and data-rates that required significant
averaging of the visibility data. This means that large numbers of sources can
only be observed through many short observations. Coupled with the sparseness
of the telescopes forming the synthesised unfilled aperture, this results in
only the brightest sources being detectable, that is, those that are
sufficiently compact to have a measurable correlated flux on the available
baselines of the array.

This limitation of VLBI means that radio sources at relatively high flux
densities ($>1$ mJy) tend to make up the target population at cm-wavelengths,
which is mainly dominated by radio-loud Active Galactic Nuclei (AGN), or so-
called jetted AGN (Padovani, 2017). At lower flux densities ($<1$ mJy) a
higher proportion of radio sources have emission associated with star-
formation processes or weak jets (e.g.,Condon et al. 2012; Prandoni et al.
2018). Studying such objects at high angular resolution requires much deeper
observations, which further limits the number of objects that can be studied
(Garrett et al., 2001).

However, due to advances in computing and the development of new techniques
for correlating the signals from the different radio telescopes, it is now
possible to form multiple phase centres for a given observation (Deller et
al., 2011b). Using multiple phase centres allows the observer to image small
areas around known sources within the same observation. It works by shifting
the phase centre to the location of known sources and averaging the obtained
visibilities to obtain manageable-sized data sets. For example, the GOODS-N
(Radcliffe et al., 2018) and COSMOS fields had 699 and 750 phase centres per
observation. However, these are somewhat special cases and routine wide-field
VLBI experiments are currently limited to around 100 phase centres from a
typical observation. Nevertheless, this technique has revolutionised our view
of the mas-scale radio Universe, both through shallow wide-field surveys, like
the mJy Imaging VLBI Exploration at 20 cm (mJIVE–20; Deller & Middelberg 2014)
survey and through deep narrow-field observations of well-studied deep fields
(e.g., Middelberg et al., 2013; Herrera Ruiz et al., 2017; Radcliffe et al.,
2018).

The results from these studies have demonstrated a proof-of-concept for all-
sky VLBI surveys with current instruments, like the European VLBI Network
(EVN) or the Very Long Baseline Array (VLBA). Also, in the future, the
sensitivity of VLBI arrays are expected to improve dramatically with the
construction of the Square Kilometre Array (SKA-VLBI) and the next generation
Very Large Array (ngVLA). These instruments will lower the current sensitivity
limits, and given their wide fields-of-view and wide frequency bandwidths,
they can provide a significant increase in the numbers of sources detected on
arcsec- to mas-scales. However, this will require determining the number of
detectable sources on VLBI-scales as a function of flux density, so that
survey strategies can be properly developed.

In this paper, we determine the source counts of VLBI-detected radio sources
from the wide-area mJIVE–20 survey (Deller & Middelberg, 2014). With these
data, we estimate the number of radio sources that can be detected by the next
generation of wide-field VLBI surveys with the VLBA and the EVN in the short
term ($<10$ yr) and with a hypothetical SKA-VLBI array in the longer term.
This paper is arranged as follows. In Section 2, we present an overview of the
final mJIVE–20 survey catalogue that we use for our study. We determine the
source counts of VLBI-detected radio sources in Section 3, which includes
calculating the mJIVE–20 survey completeness from simulations and determining
the effective area of the survey. Using these data and the (expected)
capabilities of current and next generation radio interferometers, we also
calculate the number of likely detectable radio sources on VLBI-scales for
given survey strategies in Section 4. Finally, in Section 5, we present our
conclusions.

##  2 The mJIVE–20 survey final catalogue

In this section, we provide an overview of the mJIVE–20 survey, and discuss
the properties and detection rate of the final catalogue that we use to
determine the source counts of VLBI-detected radio sources.

###  2.1 Overview of the mJIVE–20 survey

The mJIVE–20 survey was carried out with the VLBA at 1.4 GHz with dual
polarization and a bandwidth of 64 MHz (recording rate of 512 Mbit s-1). The
main aim of the project was to better understand the radio source population
at mas resolution by exploiting recent developments in wide-field multi-phase
centre correlation techniques (Deller et al., 2011b). The survey strategy and
initial results were reported by Deller & Middelberg (2014); here we summarise
the main properties of the survey and discuss the final catalogue.

Each mJIVE–20 survey observation consisted of four individual pointings of 2
mins duration around a known VLBA calibrator (so that in-beam calibration
could be used; see Fig. 1) at 7 different hour angles to improve the uv-
coverage. This resulted in a typical synthesised beam-size of $16\times 6$ mas
for target fields at Declinations between $-7$ and 60 deg. The average rms
noise is around 150 $\mu$Jy beam-1 at each pointing centre. In total, there
were 410 separate observations that targeted 24 903 radio sources from the
Faint Images of the Radio Sky at Twenty-cm (FIRST) survey (Becker et al.,
1995), making the mJIVE–20 survey the largest targeted sample of radio sources
with VLBI to date at L band. From these observations, 4 965 sources were
detected above a threshold of 6.75$\sigma$, where $\sigma$ is the local rms
noise of the imaging data. For our analysis, we use the final source catalogue
from 2014 March 31, which was obtained using the blobcat source detection
software (Hales et al., 2012).

Figure 1: An example rms map from simulating four separate pointings of the
VLBA 25-m antennas at 1.4 GHz. The VLBA calibrator source is denoted by the
(red) star in the centre of the image, such that it is always in the primary
beam of each antenna for each pointing position. A sky area of about $1\times
1$ deg2 is surveyed with this pointing strategy, for each observation of a
VLBA calibrator.

###  2.2 Detection rate as a function of radio source surface brightness,
compactness and size

We now investigate the detectability of radio sources on VLBI-scales, based on
their observed properties at lower angular resolution. Note that this
comparison is not intended to be a thorough study of the astrophysical
characteristics of the radio source population, but instead aims to compare
the properties between the parent sample catalogue obtained with the VLA and
whether a detection would be made with a snapshot observation with the VLBA.

From a straight comparison of the number of sources observed and detected, we
see that the mJIVE–20 survey has an overall detection rate of $19.9\pm 2.9$
per cent, where the uncertainty is calculated from Poisson statistics.
However, to investigate the detection rate of radio sources on VLBI-scales
further, we show in Fig. 2 the number of radio sources detected and not
detected, and the detection fraction for the mJIVE–20 survey, as a function of
the FIRST peak surface brightness. We see that above a peak surface brightness
of 80 mJy beam-1 (for a 5.5 arcsec beam size), the majority of radio sources
detected in FIRST also have a VLBI counterpart in the mJIVE–20 survey. Below
this surface brightness limit, the detection fraction steadily drops from 50
to around 31 per cent at 5 mJy beam-1, before steeply falling to 8 per cent at
the FIRST detection limit (5$\sigma$) of 0.8 mJy beam-1.

This change in the detection fraction as a function of peak surface brightness
in the parent sample could result from a number of factors. First, it could be
due to a change in the intrinsic compactness of radio sources, where the
brightest objects are relativistically beamed towards the observer, that is, a
selection effect. This would explain why the majority of radio sources at the
highest surface brightness are also detected on VLBI-scales. Also, even for
those cases that are not strongly beamed, the depth of the mJIVE–20 survey
data is sufficient to detect a weak radio core in the brightest radio sources
in FIRST. This is consistent with the results of Herrera Ruiz et al. (2017),
who detected with the VLBA a large number of low flux density radio sources
($\leq$1 mJy) that were previously identified using the VLA. Moreover, for
those observations, the VLBA recovered between 60 and 80 per cent of the VLA
flux density. This suggests that a large portion of sources at and below the
mJIVE–20 survey detection limit may still be dominated by AGN activity.

Second, it could be that those sources with a lower surface brightness do have
radio emission on VLBI-scales, but the mJIVE–20 survey observations were not
deep enough. This would explain the sharp drop in the detection rate towards a
lower surface brightness, that is, this is an observational effect. This
conclusion is also consistent with extremely deep EVN observations of the
GOODS-N field (rms 9 $\mu$Jy beam-1; Radcliffe et al. 2018) when an
e-MERLIN/VLA parent sample is used (Muxlow et al., 2020), which finds a
similar detection rate of $25\pm 6$ per cent (24 out of 94 objects) to the
mJIVE–20 survey. Even though the radio source population at this depth is
expected to be dominated by star-forming galaxies, as opposed to those
harbouring an AGN, it is interesting that providing deep enough observations
are carried out, a high fraction of VLBI detections is still found. Therefore,
the steep change in the detection fraction towards lower surface brightness
for the mJIVE–20 sample is likely a combination of the observing depth and a
change in the properties of the underlying source population.

Figure 2: Upper panel: The number of objects detected (purple) and not-
detected (blue) at VLBI-scales during the mJIVE–20 survey, as a function of
peak surface brightness measured from FIRST (5.5 arcsec beam size). Above a
peak brightness of around 90 mJy beam-1, the number of detections is greater
than the number of non-detections. Lower panel: The detection fraction at
VLBI-scales during the mJIVE–20 survey, as a function of peak surface
brightness measured from FIRST. The detection rate changes from around 50 per
cent for the highest surface brightness sources within FIRST, to below 10 per
cent at the detection threshold. This is due to the compactness of the radio
sources changing as a function of flux density.

To further investigate the latter, we have looked at several properties of the
source population to determine their effect on the mJIVE–20 survey detection
fraction. In particular, we have analysed the ratio of VLBI peak surface
brightness to VLBI integrated flux density (VLBA compactness), the ratio of
FIRST peak surface brightness to FIRST integrated flux density (VLA
compactness), and the VLBI and FIRST deconvolved major axis (taken from the
mJIVE–20 survey catalog obtained using blobcat), which are defined as $V_{\rm
comp}$, $F_{\rm comp}$, $V_{\rm size}$ and $F_{\rm size}$, respectively.
Distributions for the VLA compactness and size are presented in Fig. 3. The
same distributions for the VLBA properties, and also the ratio of peak VLBI
surface brightness to the peak FIRST surface brightness, which we refer to as
$P_{F/V}$ (radio source compactness), are presented in Fig. 4.

Figure 3: The compactness ($F_{\rm comp}$; upper panel) and major-axis size
($F_{\rm size}$; lower panel) in the FIRST survey for the detection (purple)
and non-detection (blue) samples in the mJIVE–20 survey. Figure 4: The
compactness ($V_{\rm comp}$; upper panel) and major-axis size ($V_{\rm size}$;
middle panel) for objects detected in the mJIVE–20 survey. Also shown is the
ratio of the FIRST and mJIVE–20 survey peak surface brightness ($P_{F/V}$;
lower panel).

As expected, those objects with a higher compactness and a smaller size on
VLA-scales are more likely to also have a detection on VLBI-scales. For the
VLA compactness, the mean (with standard deviation) and median of the
distributions for the VLBI detections and non-detections are 1.05
($\sigma=0.24$) and 1.00 beam-1, and 0.92 ($\sigma=0.38$) and 0.92 beam-1,
respectively. In the case of the size, the mean (with standard deviation) and
median of the distributions for the VLBI detections and non-detections are
1.60 ($\sigma=2.01$) and 1.24 arcsec, and 4.10 ($\sigma=4.10$) and 3.27
arcsec, respectively. Also, we find that those objects detected on VLBI-scales
by the mJIVE–20 survey tend to have a higher compactness and smaller intrinsic
size; the mean (with standard deviation) and median for the compactness and
size on VLBI-scales is 0.85 ($\sigma=0.27$) and 0.9 beam-1, and 11.5
($\sigma=13.0$) and 7.7 mas. Finally, we find that $P_{F/V}$ has a mean of
3.7, with a standard deviation of 13.4, and a median of 1.4.

Although the mean values of the above parameters tend to be consistent for the
samples of detections and non-detections on VLBI-scales, given the large
scatter for each parameter, the median values point toward compact sources
being more likely to be detected. Therefore, we now investigate the
distribution of sources detected, as a function of their VLBI size, $V_{\rm
size}$. In Fig. 5, we show the mJIVE–20 survey detection fraction as a
function of FIRST peak surface brightness, where each bin is also divided into
bins of $V_{\rm size}$. We consider all objects with a major axis larger than
16 mas as extended on VLBI-scales, since this is equivalent to the average
beam size of the mJIVE–20 survey. Due to the selection bias in the low surface
brightness regime ($<2$ mJy beam-1 in FIRST), only objects with a $V_{\rm
size}$ of $<8$ mas have been detected. Overall, extremely compact objects
($V_{\rm size}<1$ mas) have the least contribution to the detection fraction
compared to the more extended objects, implying that almost all radio sources
are partially-resolved with the VLBA, that is, they have some level of VLBI
structure. Another point from Fig. 5 is that for a surface brightness $>8$ mJy
beam-1, the contribution of objects with $1<V_{\rm size}<16$ mas is almost
constant. However, the contribution of those objects with $V_{\rm size}>16$
mas increases. This is likely due to the signal-to-noise ratio of the
brightest radio sources also being high, and therefore, we are more sensitive
to any possible extended emission. Similarly, extended objects make up a
smaller fraction of the lowest surface brightness part of the distribution
because here we are mainly sensitive to the most compact emission from an
object.

Overall, we find that VLBI detections are more likely from the most compact
radio sources observed at arcsec-scales, and that the emission observed on
VLBI-scales for around 65 to 80 per cent of the objects has a size in the
range of 1 to 16 mas.

Figure 5: The detection fraction for the mJIVE–20 survey, binned by peak
surface brightness in FIRST. The different colours show the detection fraction
in bins of $V_{\rm size}$ (deconvolved major axis of a single 2-dimensional
elliptical Gaussian fit to the imaging data).

##  3 VLBI-detected radio source counts

In this section, we determine the source counts of radio sources detected
during the wide-field mJIVE–20 survey. Our aim is to estimate the number of
VLBI-detected radio sources as a function of integrated flux density, so that
we can make robust estimates for the likely source counts from future surveys
with current (EVN and VLBA) and next generation (SKA-VLBI and ngVLA) VLBI
arrays.

We first calculate the sky-area observed during the mJIVE–20 survey, which due
to the primary beam attenuation of the VLBA antennas is also a function of
source surface brightness. Next we determine the completeness of the catalogue
by making simulations of mock mJIVE–20 survey data and running the source
detection procedure. With the sky area and an estimate of the completeness in
hand, we then determine the normalized Euclidean and differential source
counts for VLBI-detected radio sources. Here, we refer to the entire sample of
24 903 radio sources observed during the mJIVE–20 survey as the (FIRST) parent
population, and the 4 965 radio sources detected with the VLBA as the mJIVE–20
population.

###  3.1 The mJIVE–20 survey sky area

Calculating the radio source density on the sky requires some knowledge of the
sky area being investigated. As can be seen in Fig. 1, each mJIVE–20 survey
observation consisted of four different pointing positions around a central
calibrator source (Deller & Middelberg, 2014). This resulted in the effective
rms noise across the field being non-uniform, due to the primary beam
attenuation of the individual antennas and the combination of data from
different individual pointings and repeated observations. Also, two different
pointing configurations were employed during the mJIVE–20 survey. The first
used an offset from the central calibrator of $\pm 12$ arcmin in RA and Dec,
to include all parent population radio sources within 20 arcmin of each
pointing centre, and the second used offsets of $\pm 9.6$ arcmin in RA and
Dec, and included all parent population radio sources within 17 arcmin of each
pointing centre. A total of 306 unique calibrators were observed, with 252
observations using the set-up with $\pm 12$ arcmin offsets, and 54
observations using the set-up with $\pm 9.6$ arcmin offsets from the
calibrator.

To measure the sky area, we use a Monte Carlo integration by stone throwing
technique, which is based on the mean value theorem designed for numerical
integration using random numbers. In this method, the sky area is calculated
by the average number of hits in different iterations. The error in measuring
the sky area with this approach goes as ${1}/{\sqrt{N}}$, where $N$ is the
number of thrown stones. We first define a box around the intersected primary
beam footprints (circles), shown for example in Fig. 1, and measure the area
of that box. Then, we throw 1 million points into the box to estimate the
size. The ratio of the number of points that are inside any of the circles to
the total number of points in the box gives an estimation of the area inside
the circles. Considering the area within the individual observation
footprints, and the different pointing configurations, the total sky area of
the mJIVE–20 survey is found to be 237.95 deg2.

In Fig. 6, we show how the effective area changes as a function of the minimum
peak surface brightness of an mJIVE–20 population radio source (assuming a
local signal-to-noise ratio of 6.75). We see that at a surface brightness of
around 4 mJy beam-1, there is a knee in the effective sky area (corresponding
to 227 deg2; 95 per cent of the total sky observed), which is roughly
equivalent to 25 times the typical rms map noise of 150 $\mu$Jy beam-1 at each
pointing centre.

Figure 6: The effective area of the mJIVE–20 survey, as a function of minimum
detectable radio source surface brightness (for a detection threshold of
$6.75\sigma$). The total observed area in the mJIVE–20 survey is 237.95 deg2,
but changes as a function of source brightness due to the primary beam
attenuation and the combination of the visibility data from different
pointings and repeated observations.

###  3.2 The mJIVE–20 survey completeness

Figure 7: The cumulative completeness for a set of simulated mock VLBA
observations that have the same characteristics as the mJIVE–20 survey data.
These have been analyzed using the blobcat object detection algorithm, with
the same parameters for the detection thresholds ($T_{d}$ and $T_{f}$; see
text for details) used for the mJIVE–20 survey.

In order to calculate the source counts of the mJIVE–20 population to the
lowest flux densities, we also have to consider the completeness of the
catalogue, that is, estimate how many sources may have been missed due to
observational effects. To do this, we have made mock visibility data sets,
which then go through the same imaging and source detection process as the
mJIVE–20 survey data. By comparing the input and output mock catalogues, we
determine the completeness as a function of flux density.

The simulated data set needed to have similar properties to the radio sources
in the mJIVE–20 population catalogue. Therefore, we implemented a Monte Carlo
approach to define the simulated radio source size (either delta or
2-dimensional elliptical Gaussian function), peak surface brightness, position
angle, and a random $x,y$ position with respect to the phase centre, given the
actual distribution of these parameters in the mJIVE–20 catalogue (see Rezaei
et al. (2022) for a detailed overview of generating realistic mock sources for
mJIVE–20 survey). The generated sources were then converted to mock
visibilities using the Common Astronomy Software Applications (casa; McMullin
et al. 2007) package. The simulation tool in casa allows as an input an
existing interferometric visibility data set that can be used to obtain the
position of the antennas and other observational settings (such as frequency
and time sampling, total integration time). Therefore, we used the actual
mJIVE–20 survey interferometric visibility data sets when converting the
simulated model radio sources into mock visibility data. The final step in the
process used the tclean task within casa to produce de-convolved clean images.
As with the mJIVE–20 survey data, the images were centred on the position of
the surface brightness peak, had a size of $1024\times 1024$ pixels and a
pixel-scale of 0.75 mas pixel-1. The dirty images were cleaned for a maximum
of 1000 iterations or until a threshold of 0.2 mJy beam-1 was reached, while
masking was applied to the centre of the image, with a radius of 25 pixels.

To identify and measure the flux density of the sources in our mock imaging
data, we have used blobcat (Hales et al., 2012), the same post-processing
object detection package used as part of the mJIVE–20 survey. blobcat searches
for islands of pixels that could represent a possible object considering the
signal-to-noise ratio of a given pixel. This is accomplished using the surface
brightness distribution of the 2-dimensional imaging data and the background
rms noise as the input parameters. It requires the user to set a threshold for
the minimum detection signal-to-noise ratio (known as $T_{d}$) as well as a
cutoff ratio (known as $T_{f}$) for flooding the islands. Based on simulations
carried out by Deller & Middelberg (2014), $T_{d}$ and $T_{f}$ were set to
$6.5\sigma$ and $5\sigma$, respectively, during the mJIVE–20 survey. However,
only those sources with a peak surface brightness above $6.75\sigma$ were
added to the final mJIVE–20 survey catalogue, so we adopt that value for
$T_{d}$ in our simulations. In Fig. 7, we show the cumulative completeness of
the mJIVE–20 survey as a function of signal-to-noise ratio from 13 500
simulated radio sources that have a signal-to-noise ratio between 3 and 15.
The magenta dashed line shows the $6.75\sigma$ detection threshold of the
mJIVE–20 survey. We find from our simulations that blobcat is 97 per cent
complete at a threshold of $6.75\sigma$ and reaches full completeness at a
threshold of $7.8\sigma$ for this data set.

From Fig. 7 we are able to calculate the completeness of the mJIVE–20 survey,
given the observational set-up (signal-to-noise ratio, uv-coverage, image de-
convolution) and the ability of the object detection algorithm to identify a
representative sample of realistic radio sources. This allows us to define a
completeness correction factor for each flux density bin, where those bins
below 100 per cent completeness are divided by the completeness given in Fig.
7; for this, we assume that the radio sources are mostly unresolved, which is
consistent with the results presented in the previous section.

###  3.3 A note on the resolution bias of the parent and mJIVE–20 population
samples

As the mJIVE–20 survey is not a blind search on VLBI-scales, but is instead a
targeted search of known FIRST radio sources, the completeness of that parent
population can also affect the final source counts of VLBI-detected radio
sources. In fact, it is known that below a flux density of around 2 mJy, the
completeness of the FIRST survey becomes progressively worse (White et al.,
1997), and the observed source counts dramatically drop-off from what is
predicted from models of radio source evolution and observations of deep
fields at the same frequency (e.g., Prandoni et al. 2018). This results in the
source counts of FIRST being lower by a factor of between 0.9 and 0.4 towards
the faint end of the flux density distribution. However, as this effect is
thought to be due to faint and extended radio sources being resolved out by
the 5.5 arcsec synthesised beam of the VLA (B-configuration; see White et al.
1997 for further discussion), it is almost certainly the case that such
objects would also be resolved out at VLBI-scales (in the absence of
significant variability). Therefore, we do not correct for the resolution bias
of the FIRST survey.

However, there will be a resolution bias associated with the mJIVE–20 survey
data, as the detectability and measured flux density of a given object will be
dependent on the structure (compactness) and also the uv-coverage of the
observations. For this reason, the source counts from deep VLBI observations
have tended to use the total flux density from lower resolution imaging
(e.g.,the VLA; Middelberg et al. 2013; Herrera Ruiz et al. 2018) when
discussing VLBI-detected radio sources with respect to the radio source
population in general. Here, we do use the flux density recovered on VLBI-
scales, as we are not primarily interested in testing galaxy formation models
or to separate AGN and star-forming galaxies, but instead, our goal is to
understand the number of radio sources that could be detected with wide-field
VLBI surveys. In Section 4, we will discuss the effect of VLBI resolution bias
further and present ways to mitigate it in the future.

###  3.4 Euclidean-normalized and differential source counts

We now calculate the source counts of the radio sources detected during the
mJIVE–20 survey. We also compare these with the source counts obtained from
studies of compact radio sources selected at higher frequencies. In Table 1,
we present the number of sources per flux-density bin, and the correction
factor that is determined by taking the effective sky area and completeness of
the survey into account. Here, we also present the (corrected) Euclidean-
normalized source counts for the mJIVE–20 survey.

In Fig. 8, we present these source counts as a function of flux density, which
we see are well approximated by a power-law down to flux densities of around 2
mJy, below which there is evidence of a knee in the distribution. This
downturn is likely associated with the change in the parent population of
FIRST radio sources described in the previous section. Below 1 mJy, the
Euclidean-normalized source counts become flatter, but are also quite noisy.
At this stage, it is not clear if this is a systematic or a real effect,
similar to the flattening of the radio source counts seen in lower angular
resolution observations at this frequency; this is attributed to a change in
the radio source population from AGN to star-forming dominated systems (Condon
et al., 2012; Prandoni et al., 2018). From a power-law fit to the data
presented in Fig. 8, we find that the Euclidean-normalized source counts can
be described by,

$n(S)=({2.4\pm 0.1})\left(\frac{S_{1.4}}{\rm 1~{}mJy}\right)^{({0.695\pm
0.013})}\quad{\rm Jy^{1.5}~{}sr^{-1}},$ (1)

where $n(S)$ is the differential number of sources (Euclidean-normalized) and
$S_{1.4}$ is the flux density at 1.4 GHz.

Figure 8: The Euclidean-normalized source counts for VLBI-detected radio
sources from the mJIVE–20 survey.

In Fig. 9, we present the differential source counts as a function of flux
density for the FIRST parent population and the mJIVE–20 population. From a
fit to the data between 1 and 100 mJy, we find that the FIRST parent
population can be described by the power-law,

$n(S)=({80.6\pm 5.5})\left(\frac{S_{1.4}}{\rm 100~{}mJy}\right)^{-({1.77\pm
0.02})}\quad{\rm mJy^{-1}~{}sr^{-1}},$ (2)

and that the mJIVE–20 parent population can be described by the power-law,

$n(S)=({19.4\pm 1.4})\left(\frac{S_{1.4}}{\rm 100~{}mJy}\right)^{-({1.74\pm
0.02})}\quad{\rm mJy^{-1}~{}sr^{-1}}.$ (3)

As in the previous case, $n(S)$ is the differential number of sources and
$S_{1.4}$ is the flux density at 1.4 GHz. From these fits, and by inspecting
Fig. 9, we see that the power-laws for the parent and mJIVE–20 populations are
diverging toward lower flux densities, which would be consistent with the
change in the detection fraction seen in Fig. 2.

Taking the result for the mJIVE–20 survey, we compare these source counts with
those for a complete sample of 117 flat-spectrum radio sources selected and
observed at 4.85 GHz with the VLA (McKean et al., 2007). Such sources are
assumed to be compact (size $<170$ mas; Myers et al. 2003), given their flat-
radio spectra (due to the super-position of synchrotron self-absorption from
many homogeneous emitting regions) and the higher selection frequency tending
to favour core-dominated radio sources. Therefore, flat-spectrum radio sources
are prime targets for VLBI observing programmes as they are expected to
dominate the compact radio source population.

We find that the slope of the source counts is steeper [$\eta=-2.06\pm 0.01$;
where $n(s)=k\,S^{\eta}$], and the normalization is smaller ($k=6.91\pm 0.42$)
for flat-spectrum radio sources selected at 4.85 GHz (Jodrell Bank Flat-
spectrum radio source sample; JBF; McKean et al. 2007), when compared to the
results for the VLBI-detected radio sources at 1.4 GHz. However, by assuming a
mean spectral index of $\alpha_{1.4}^{4.85}=-0.09$ (where
$S_{\nu}\propto\nu^{\alpha}$, and $\nu$ is the observing frequency; McKean et
al. 2007) for the JBF sample, we have calculated the expected number of
compact radio sources at 1.4 GHz down to a flux density of 1 mJy.
Interestingly, we find that the sky-density of flat-spectrum radio sources
selected at 4.85 GHz is very similar to that of the VLBI-detected radio
sources from the mJIVE–20 survey, but is a factor $1.22\pm 0.12$ higher. Such
similar source counts suggests that both samples are drawn from a similar
population, even though one is selected based on the radio spectra, and the
other on compactness at VLBI-scales. This has implications when choosing an
appropriate observing frequency for an all-sky VLBI survey, which we discuss
in the next section.

Figure 9: The differential source counts for VLBI-detected radio sources using the mJIVE–20 survey (red), with the FIRST parent population (blue) for comparison. The dashed lines are the best-fit power-law models to the mJIVE–20 survey and FIRST data. $S_{l}$ | $S_{u}$ | $\overline{S}$ | Number | Corr. | Counts
---|---|---|---|---|---
(mJy) | (mJy) | (mJy) |  |  | (sr-1 Jy1.5)
0.51 | 0.56 | 0.534 | 11 | 11.858 | 1.107$\pm$0.334
0.57 | 0.62 | 0.595 | 19 | 11.858 | 1.787$\pm$0.410
0.63 | 0.68 | 0.657 | 20 | 7.836 | 1.140$\pm$0.255
0.69 | 0.75 | 0.725 | 40 | 5.553 | 1.007$\pm$0.159
0.76 | 0.83 | 0.799 | 54 | 4.416 | 1.082$\pm$0.147
0.84 | 0.92 | 0.884 | 86 | 3.853 | 1.700$\pm$0.183
0.93 | 1.02 | 0.976 | 102 | 3.138 | 1.741$\pm$0.172
1.03 | 1.13 | 1.080 | 106 | 2.460 | 1.373$\pm$0.133
1.14 | 1.24 | 1.190 | 131 | 2.178 | 1.715$\pm$0.150
1.25 | 1.38 | 1.313 | 163 | 1.966 | 2.244$\pm$0.176
1.39 | 1.52 | 1.454 | 185 | 1.798 | 2.603$\pm$0.191
1.53 | 1.68 | 1.603 | 224 | 1.616 | 3.198$\pm$0.214
1.69 | 1.86 | 1.774 | 194 | 1.541 | 3.087$\pm$0.222
1.87 | 2.06 | 1.961 | 230 | 1.427 | 3.927$\pm$0.259
2.07 | 2.27 | 2.164 | 203 | 1.342 | 3.785$\pm$0.266
2.28 | 2.51 | 2.394 | 213 | 1.283 | 4.420$\pm$0.303
2.52 | 2.78 | 2.652 | 207 | 1.231 | 4.875$\pm$0.339
2.79 | 3.07 | 2.930 | 190 | 1.168 | 4.852$\pm$0.352
3.08 | 3.40 | 3.227 | 188 | 1.140 | 5.394$\pm$0.393
3.41 | 3.75 | 3.582 | 153 | 1.118 | 5.083$\pm$0.411
3.76 | 4.15 | 3.959 | 186 | 1.086 | 6.940$\pm$0.509
4.16 | 4.59 | 4.360 | 179 | 1.070 | 7.578$\pm$0.566
4.61 | 5.07 | 4.854 | 134 | 1.057 | 6.632$\pm$0.573
5.08 | 5.61 | 5.349 | 143 | 1.045 | 8.073$\pm$0.675
5.62 | 6.20 | 5.902 | 154 | 1.036 | 9.964$\pm$0.803
6.21 | 6.85 | 6.527 | 115 | 1.033 | 8.633$\pm$0.805
6.86 | 7.57 | 7.220 | 119 | 1.031 | 10.441$\pm$0.957
7.58 | 8.37 | 7.982 | 88 | 1.032 | 8.939$\pm$0.953
8.38 | 9.25 | 8.791 | 85 | 1.030 | 9.920$\pm$1.076
9.26 | 10.23 | 9.797 | 94 | 1.029 | 13.135$\pm$1.355
10.24 | 11.29 | 10.736 | 81 | 1.028 | 12.787$\pm$1.421
11.33 | 12.47 | 11.896 | 76 | 1.028 | 14.038$\pm$1.610
12.51 | 13.80 | 13.073 | 66 | 1.026 | 13.986$\pm$1.722
13.86 | 15.25 | 14.573 | 64 | 1.026 | 16.320$\pm$2.040
15.28 | 16.85 | 15.982 | 60 | 1.026 | 17.161$\pm$2.215
16.90 | 18.65 | 17.780 | 55 | 1.026 | 18.571$\pm$2.504
18.71 | 20.59 | 19.624 | 51 | 1.022 | 19.879$\pm$2.784
20.68 | 22.78 | 21.701 | 49 | 1.023 | 22.077$\pm$3.154
22.85 | 25.15 | 23.926 | 38 | 1.022 | 19.873$\pm$3.224
25.31 | 27.84 | 26.594 | 44 | 1.020 | 27.246$\pm$4.108
27.98 | 30.57 | 29.148 | 32 | 1.020 | 22.553$\pm$3.987
30.80 | 34.00 | 32.246 | 30 | 1.018 | 24.583$\pm$4.488
34.13 | 37.52 | 35.756 | 36 | 1.017 | 34.863$\pm$5.811
37.66 | 41.52 | 39.615 | 29 | 1.016 | 32.614$\pm$6.056
41.62 | 45.72 | 43.285 | 24 | 1.015 | 30.452$\pm$6.216
45.97 | 50.07 | 48.109 | 21 | 1.015 | 31.741$\pm$6.926
50.93 | 55.93 | 53.434 | 20 | 1.012 | 35.774$\pm$7.999
56.82 | 61.93 | 59.833 | 18 | 1.013 | 38.466$\pm$9.066
62.29 | 68.49 | 65.218 | 21 | 1.011 | 50.615$\pm$11.045
69.47 | 74.11 | 71.762 | 13 | 1.008 | 36.042$\pm$9.996
75.97 | 83.62 | 79.061 | 19 | 1.006 | 60.303$\pm$13.834
84.83 | 90.96 | 87.906 | 9 | 1.005 | 34.991$\pm$11.664
92.70 | 100.87 | 96.787 | 9 | 1.007 | 39.090$\pm$13.030
104.92 | 112.86 | 109.495 | 6 | 1.005 | 32.393$\pm$13.224
125.63 | 137.61 | 132.954 | 9 | 1.003 | 63.024$\pm$21.008
140.59 | 151.53 | 147.273 | 9 | 1.003 | 75.359$\pm$25.120
153.07 | 166.88 | 160.351 | 8 | 1.003 | 74.221$\pm$26.241
169.90 | 184.55 | 179.800 | 7 | 1.003 | 79.416$\pm$30.017
210.57 | 225.90 | 217.562 | 6 | 1.004 | 87.566$\pm$35.749
310.54 | 339.24 | 321.093 | 6 | 1.001 | 149.152$\pm$60.891
Table 1: The Euclidean normalized source counts for the mJIVE–20 survey. The
columns from left to right are, the lower ($S_{l}$) and upper ($S_{u}$) bounds
of a flux density bin, the mean flux density of a bin ($\overline{S}$), the
number of sources in each bin ($N$), the correction factor (taking into
account the sky area and the completeness; Corr.) and the resulting Euclidean
normalized source counts (Counts). Note that 33 sources with flux densities
$>340$ mJy are not included here.

##  4 Prospects for all-sky VLBI surveys

In this section, we use the differential source counts determined above to
estimate the number of radio sources that could be detected from an all-sky
survey at mas-scale resolution with VLBI. We focus primarily on what would be
achievable with the current VLBA and EVN (using only the smaller dish radio
telescopes of the array in a hypothetical survey mode), and briefly discuss
the expectations from next generation instruments, like SKA-VLBI. For our
calculations, we adopt a survey strategy that maximizes the number of radio
sources detected on VLBI-scales.

###  4.1 Characteristics of the arrays

Our ability to survey the sky is dependent on the effective field-of-view of
the individual dishes and their sensitivity. This requires some knowledge of
the primary beam models, which we assume are consistent with a circular
aperture that has an axisymmetric Gaussian illumination pattern, such that the
primary beam diameter at the 50 per cent attenuation point is given by,

$\theta_{\rm FWHM}=1.22\lambda/D,$ (4)

where $D$ is the physical diameter of the telescope and $\lambda$ is the
observing wavelength. The theoretical sensitivity of an interferometic array
is defined from the rms thermal noise fluctuations, given by

$\sigma_{\rm rms}=\frac{1}{\eta_{\rm cor}}\frac{2kT_{\rm sys}}{A_{\rm
eff}\sqrt{N_{p}\,N_{A}(N_{A}-1)\,\Delta\nu\,\Delta t_{\rm int}}},$ (5)

where $T_{\rm sys}$ is the system temperature (assuming identical antenna
systems), $k$ is the Boltzmann constant, $\eta_{\rm cor}$ is the correlator
efficiency (assumed to be 0.89 for 2-bit sampling), $A_{\rm eff}$ is the
effective collecting area of the filled aperture, $N_{p}$ is the number of
polarisations (assumed to be 2), $N_{A}$ is the number of antennas in the
array, $\Delta\nu$ is the bandwidth in frequency and $\Delta t_{\rm int}$ is
the integration time. Typically, the antenna forward gain is defined by the
System Equivalent Flux Density (SEFD), which is given by

${\rm SEFD}=\frac{2kT_{\rm sys}}{A_{\rm eff}}.$ (6)

In the cases where the antenna systems are not identical, the SEFD on a given
baseline between antennas 1 and 2 is expressed as $\sqrt{{\rm
SEFD_{1}}\times{\rm SEFD_{2}}}$. Note that the VLBA and EVN are assumed to
operate at central observing frequencies of 1.4 and 1.7 GHz, respectively,
given the available frequency coverage of the individual receiver systems.

We see from Equation 4 that the observable solid angle is proportional to
$\lambda^{2}$, and therefore, surveys of radio sources at longer wavelengths
are typically more efficient. This is also because radio sources tend to have
a higher flux density toward lower frequencies. However, we see from our
comparison of the source counts for the mJIVE–20 and JBF surveys above that
the expected number of detectable radio sources on VLBI-scales is likely quite
similar at 1.4/1.7 and 5 GHz. Therefore, either frequency is likely viable.
However, the level of radio frequency interference at 5 GHz is currently much
less than at 1.4/1.7 GHz, and the usable bandwidth is currently a factor of
two to four times better. This, coupled with a slightly better receiver
response at 5 GHz, results in around a factor of between 1.5 and 2 improvement
in overall sensitivity at 5 GHz. Unfortunately, as this comes at the expense
of having to carryout a factor of around 9 to 13 more observations to cover
the same area of sky, we focus on 1.4/1.7 GHz as our preferred observing
frequency for a potential all-sky survey with VLBI. We note that this choice
is based on maximising the number of radio sources detected on VLBI-scales,
but other science goals may be better served with the improved angular
resolution and frequency coverage afforded at 5 GHz.

In addition, we also consider a hypothetical VLBI array that includes 20
antennas with the same characteristics as the SKA-MID. We assume that they are
positioned in different locations in Botswana, Ghana, Kenya, Madagascar,
Mauritius, Mozambique, Namibia and Zambia (the partner African VLBI Network
countries) with two dishs per country, and a further four antennas spread
across South Africa. The exact locations of these antennas can be seen in the
right-hand panel of Figure 10 Also, we assume that 512 MHz of frequency
bandwidth will be usable within the SKA Band 2, operating between 0.95 and
1.76 GHz, due to radio frequency interference.

In Table 2, we summarise the properties of the three arrays that we consider
here.

Array | $N_{A}$ |  Antennas | Diameter | Freq. | $\Delta\nu$ | FoV | SEFD
---|---|---|---|---|---|---|---
|  |  | (m) | (GHz) | (MHz) | (deg2) | (Jy)
VLBA | 10 |  Sc, Hn, Nl, Fd, La, Kp, Pt, Ov, Br, Mk | 25 | 1.4 | 256 | 0.262 | 365
EVN/e-MERLIN | 18 |  Mc, On, Tr, W1, Nt, Sh, Ur, Hh, Sv, Zc, |  |  |  |  |
|  |  Bd, Ir, Jb2, Cm, Da, De, Kn, Pi | 25, 32 | 1.7 | 128 | 0.117 | 485
SKA-VLBI | 20 |  AVN1-BW, AVN2-BW, AVN3-GH, AVN4-GH, AVN5-KE, AVN6-KE, AVN7-MG, AVN8-MG, AVN9-MU, AVN10-MU, AVN11-MZ, AVN12-MZ, AVN13-NAM, AVN14-NAM, AVN15-ZM, AVN16-ZM, AVN17-ZA, AVN18-ZA, AVN19-ZA, AVN20-ZA | 15 | 1.4 | 512 | 0.728 | 235
Table 2: The properties of the two arrays considered here for an all-sky
survey with VLBI. The antenna names are given using their standard
abbreviation. The portion of EVN/e-MERLIN listed here has a combination of
$10\times 25$\- and $8\times 32$-m telescopes, but for our field-of-view (FoV)
calculations, we use a diameter of 32 m. The SEFD is the average for the given
array. We also include a hypothetical VLBI array that includes 20 antennas
with the characteristics of the SKA-MID dish design, which provides mas-scale
angular resolution (excluding the SKA-MID core).

###  4.2 Prospects for in-beam calibration

Our ability to calibrate VLBI observations is currently limited by the density
of objects that are suitable for correcting the complex antenna gains as a
function of time. Also, given the large changes in phase, due to the long
baselines involved with VLBI observations, calibrators are often required to
be very close to the target field ($<2$ deg separation). There are currently
17 432 compact radio sources that have been identified as calibrators for VLBI
experiments (1 in 2.37 deg2 over the whole sky; Radio Fundamental Catalog;
Petrov 2021). However, astrometric accuracy scales with the calibrator
distance (Pradel et al., 2006), so additional calibrators will improve the
overall astrometric and imaging quality. It is this rationale that made in-
beam calibration ($<0.4$ deg calibrator-to-target separation) a key component
of the mJIVE–20 survey.

With the development of wide-bandwidth observations, radio sources at lower
flux densities can now be used as calibrators, when the amplitude and phase of
the visibilities are well calibrated as a function of frequency. For example,
the recommended phase referencing time at 1 to 8 GHz in good weather is 300 s.
For the VLBA and EVN arrays presented in Table 2, this is equivalent to a
baseline sensitivity of $\sigma_{\rm rms}=0.74$ and 0.98 mJy beam-1 at 1.4 and
1.7 GHz, respectively. Therefore, compact radio sources with flux-densities
$>7.4$ and $>9.8$ mJy (for $10\sigma_{\rm baseline}$) can be regarded as good
phase reference sources for the VLBA and EVN, respectively. For our
hypothetical SKA-VLBI array, the baseline sensitivity is $\sigma_{\rm
rms}=0.39$ mJy beam-1, which corresponds to a minimum calibrator flux density
of $>3.9$ mJy (for $10\sigma_{\rm baseline}$).

Using the differential source counts of VLBI-detected radio sources given by
equation 3, we calculate the expected sky density of radio sources needed for
phase referencing from pointed observations, and determine the fraction that
can be used for in-beam calibration. Note that the latter requires increasing
the lower limit on the flux density by a factor of 2, due to the attenuation
of the primary beam of the antennas. For the VLBA, we find that the sky
density of radio sources with flux densities $>7.4$ (pointed) and $>14.8$ mJy
(in beam) is 5.5 and 3.3 deg-2, respectively. This corresponds to an in-beam
calibrator density of 0.87 sources beam-1. In the case of the EVN, we expect
that the sky density of radio sources with flux densities $>9.8$ (pointed) and
$>19.6$ mJy (in beam) is 4.5 and 2.6 deg-2, respectively, with an in-beam
density of 0.30 sources beam-1 (for the 32-m antennas). Finally, for our
hypothetical SKA-VLBI array, the sky density of radio sources with flux
densities $>3.9$ (pointed) and $>7.8$ mJy (in beam) is 8.8 and 5.3 deg-2,
respectively. Therefore, the number of radio sources that can be used for in-
beam calibration with the SKA-VLBI array would be 3.9 sources beam-1.

Overall, we find that 1 in 3 pointings with the EVN will likely have an in-
beam calibrator, whereas this is likely the case for almost all pointings with
the VLBA (Ding et al., 2023; Deller et al., 2011a). In the case of the SKA-
VLBI array, we predict that there will be multiple in-beam calibrators
available, which can be important for extremely precise astrometric
applications of VLBI (e.g., Dodson & Rioja 2018).

###  4.3 Expected number of detected sources

Figure 10: Left-hand panel: The rms sensitivity against declination for a
theoretical AVN plus EVN array (top) and an AVN-only array (bottom) for a 10
min observation and 512 MHz bandwidth. Right-hand panel: The locations of the
theoretical 20 element AVN array plotted using basemap (Elson et al., 2022).
The positions of 16 of these telescopes were selected based upon the two
largest cities of SKA African partner countries while the remaining are the
four largest South African cities.

We now calculate the number of radio sources that can be found from an all-sky
survey with VLBI. For this, we integrate equation 3 to an appropriate flux-
density limit. We see from Fig. 8 that our differential source counts of VLBI-
detected radio sources is robust to around 1 mJy, below which the source
counts become much more noisy. In principle, we could extrapolate our
calculations below 1 mJy, but given the expected change in the radio source
population at these flux densities, this would likely result in an over-
estimate of the total number of radio sources that we would detect. Therefore,
we use a 1.4/1.7 GHz flux-density limit of 1 mJy for our calculations.

We also assume that, given the short amount of time used per observation, any
survey would have similar noise properties to the mJIVE–20 survey, and that
the same detection threshold would likely be needed ($6.75\sigma$); this would
result in a similar completeness of 97 per cent, as we determined from Fig. 7,
which we fold into our calculations. We also assume that the individual
pointings are separated by a beam width, in a similar configuration to what
was presented in Fig. 1. In a future work, we will simulate the effect of
different pointing strategies and noise realizations. However, for our current
pointing strategy, we see that the absolute value of the noise changes across
the field. Therefore, to make a $6.75\sigma$ detection at a flux density
(point-source) limit of 1 mJy, means having an average thermal noise of 150
$\mu$Jy beam-1. This corresponds to a thermal noise at the pointing centre
(where the primary beam response is maximum) of 130 $\mu$Jy beam-1
($1\sigma_{\rm rms}$). From the properties of the three arrays given in Table
2, and using equations 5 and 6, we find that the on-source integration time
needed to reach this sensitivity would be 215, 225 and 15 s for the VLBA, EVN
and a hypothetical SKA-VLBI array, respectively. This highlights the
improvement in sensitivity that is potentially available with SKA-VLBI, given
the advances in antenna and receiver design for this next generation
instrument. We note that the short integration times needed to reach the
required noise levels would allow for flexible scheduling of such
observations. Therefore, it would also be straightforward to define optimum
scheduling blocks when the target fields are visible to all of the available
telescopes. This is more challenging to achieve for long-track observations
that are typically used for deep-field science cases.

We also see from equation 3 that the slope of the differential number counts
of VLBI-detected radio sources is quite shallow ($\eta=-1.74\pm 0.02$).
Therefore, maximizing the number of detected radio sources is best achieved
through observing the largest possible sky area, which for our calculations we
assume is $3\pi$ sr (equivalent to about 31 000 deg2). Given a flux-limited
survey to 1 mJy, and a completeness of 97 per cent at this limit, such a
survey could potentially detect $(7.2\pm 0.9)\times 10^{5}$ radio sources on
VLBI-scales, a factor of about 30 times more than is currently known.

Finally, we estimate the total time such a survey would take with the current
VLBA and EVN (using only the 32-m antennas), and a next generation instrument,
like SKA-VLBI. From Table 2, we see the effective field-of-view of each array.
From this we estimate that to survey $3\pi$ sr would take $1.18\times 10^{5}$,
$2.64\times 10^{5}$ and $1.18\times 10^{4}$ pointings for the VLBA, EVN and a
hypothetical SKA-VLBI array, respectively. Given the integration time to reach
the required noise level, this equates to about 7000, 16500 and 180 h of on-
source observing time for the VLBA, EVN and SKA-VLBI, respectively. Note that
about 2 in 3 of the EVN pointings will also require additional phase
referencing calibration (see above), which will add up to around another 30
per cent to the total time needed to complete such a survey with that array.
For this reason, such surveys with the EVN are likely prohibitive, whereas
with the VLBA and SKA-VLBI they would be very much feasible.

##  5 Conclusions

We have analyzed the final catalogue of the mJIVE–20 survey, which observed 24
903 radio sources in the FIRST survey using the VLBA at 1.4 GHz to an rms
noise level of about 150 $\mu$Jy beam-1, detecting 4 965 radio sources on
VLBI-scales. Through comparing the number of detections and non-detections in
the mJIVE–20 catalogue, we found that the detection fraction is a strong
function of the peak surface brightness (compactness and radio source size) of
the objects in the FIRST survey; the detection fraction is over 50 per cent at
a peak surface brightness at 80 mJy beam-1 and falls steadily to about 31 per
cent at 5 mJy beam-1. Below this, the detection fraction falls sharply to 8
per cent at the detection threshold of both surveys. This is likely due in
part to a change in the composition of the radio source population, from AGN
to star-formation dominated objects, but is also due to the VLBI observations
not being sensitive enough to detect low-level compact emission toward lower
flux densities. We found an overall detection fraction of $19.9\pm 2.9$ per
cent. Given the limited uv-coverage of the mJIVE–20 survey observations, we
found that those radio sources that are detected tend to be unresolved (median
VLBI compactness 0.9 beam-1), with a median size of 7.7 mas (about half a
beam-size). Finally, we found that 20 to 35 per cent of the radio sources
detected were resolved, with sizes $>16$ mas.

From an analysis of the VLBI-detected source counts, we see hints of a similar
behaviour in the Euclidean-normalized distribution that has been reported on
arcsec-scales, that is, a downturn around a flux density of 2 mJy and then a
potential flattening below 1 mJy. This might arise from the considerable
challenge in accurately determining the size and, by extension, the integrated
flux density of sources with low signal-to-noise ratios. During the
deconvolution process, if the calculated size is smaller than the actual beam
size, the resulting integrated flux could be significantly lower than the peak
flux density. As a result, these measurements are accompanied by substantial
error bars. We also determine the differential number counts for VLBI-detected
radio sources, finding that the total number of objects is similar to those of
flat-spectrum radio sources selected at higher frequencies, suggesting that
they come from a similar population.

From our analysis of the differential source counts, we found that the sky
density of suitable phase reference sources is of order 2.6 to 3.3 deg-2. This
should be sufficient for in-beam phase referencing in the case of around 30
and 90 per cent of the observations carried out with the EVN and VLBA,
respectively. However, we found that for a VLBI facility that includes
antennas with a similar specification to the SKA-MID design, the expected sky
density of phase reference sources is about 5.3 deg-2, which equates to
multiple in-beam calibration sources for each observation.

Finally, we investigated the number of sources that could be found from all-
sky surveys carried out with the VLBA and EVN. From this analysis, we found
that a factor 30 more VLBI-detected radio sources could be identified with
around 7000 h of observations with the VLBA. However, in the case of the EVN,
such surveys would be rather expensive, given the smaller field of view of the
antennas in that array. For a hypothetical SKA-VLBI array, such a survey would
take a fraction of the time, and could be completed with just 180 h of
observations. Our analysis is currently limited by our knowledge of the radio
source counts below 1 mJy. More focused surveys of the radio sky with VLBI,
down to a limiting sensitivity of around 100 $\mu$Jy beam-1 ($6\sigma$
detection threshold) would be extremely informative in testing models for
radio source populations and making robust predictions for the expectations
with SKA-VLBI.

## Acknowledgements

This paper is based on research developed in the DSSC Doctoral Training
Programme, co-funded through a Marie Skłodowska-Curie COFUND (DSSC 754315).
JPM acknowledges support from the Netherlands Organization for Scientific
Research (NWO) (Project No. 629.001.023) and the Chinese Academy of Sciences
(CAS) (Project No. 114A11KYSB20170054). The National Radio Astronomy
Observatory is a facility of the National Science Foundation operated under
cooperative agreement by Associated Universities, Inc.

## Data Availability

The data used in this research is publicly available in the mJIVE–20 survey
and FIRST databases.

## References

  * Becker et al. (1995) Becker R. H., White R. L., Helfand D. J., 1995, ApJ, 450, 559
  * Broten et al. (1967) Broten N. W., et al., 1967, Science, 156, 1592
  * Cohen et al. (1977) Cohen M. H., et al., 1977, Nature, 268, 405
  * Condon et al. (1998) Condon J. J., Cotton W. D., Greisen E. W., Yin Q. F., Perley R. A., Taylor G. B., Broderick J. J., 1998, AJ, 115, 1693
  * Condon et al. (2012) Condon J. J., et al., 2012, ApJ, 758, 23
  * Deller & Middelberg (2014) Deller A. T., Middelberg E., 2014, AJ, 147, 14
  * Deller et al. (2011a) Deller A. T., et al., 2011a, in 20th Meeting of the European VLBI Group for Geodesy and Astronomy. pp 178–182 (arXiv:1110.1979), doi:10.48550/arXiv.1110.1979
  * Deller et al. (2011b) Deller A. T., et al., 2011b, PASP, 123, 275
  * Ding et al. (2023) Ding H., et al., 2023, MNRAS, 519, 4982
  * Dodson & Rioja (2018) Dodson R., Rioja M., 2018, in 14th European VLBI Network Symposium & Users Meeting (EVN 2018). p. 86
  * Elson et al. (2022) Elson P., et al., 2022, SciTools/cartopy: v0.20.3, doi:10.5281/zenodo.6775197, https://doi.org/10.5281/zenodo.6775197
  * Event Horizon Telescope Collaboration et al. (2019a) Event Horizon Telescope Collaboration et al., 2019a, ApJ, 875, L1
  * Event Horizon Telescope Collaboration et al. (2019b) Event Horizon Telescope Collaboration et al., 2019b, ApJ, 875, L2
  * Garrett et al. (2001) Garrett M. A., et al., 2001, A&A, 366, L5
  * Hales et al. (2012) Hales C. A., Murphy T., Curran J. R., Middelberg E., Gaensler B. M., Norris R. P., 2012, MNRAS, 425, 979
  * Herrera Ruiz et al. (2017) Herrera Ruiz N., et al., 2017, A&A, 607, A132
  * Herrera Ruiz et al. (2018) Herrera Ruiz N., et al., 2018, A&A, 616, A128
  * Intema et al. (2017) Intema H. T., Jagannathan P., Mooley K. P., Frail D. A., 2017, A&A, 598, A78
  * Lacy et al. (2020) Lacy M., et al., 2020, PASP, 132, 035001
  * Marcote et al. (2020) Marcote B., et al., 2020, Nature, 577, 190
  * McKean et al. (2007) McKean J. P., Browne I. W. A., Jackson N. J., Fassnacht C. D., Helbig P., 2007, MNRAS, 377, 430
  * McMullin et al. (2007) McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, in Shaw R. A., Hill F., Bell D. J., eds, Astronomical Society of the Pacific Conference Series Vol. 376, Astronomical Data Analysis Software and Systems XVI. p. 127
  * Middelberg et al. (2013) Middelberg E., et al., 2013, A&A, 551, A97
  * Miyoshi et al. (1995) Miyoshi M., Moran J., Herrnstein J., Greenhill L., Nakai N., Diamond P., Inoue M., 1995, Nature, 373, 127
  * Morganti et al. (2013) Morganti R., Fogasy J., Paragi Z., Oosterloo T., Orienti M., 2013, Science, 341, 1082
  * Muxlow et al. (2020) Muxlow T. W. B., et al., 2020, MNRAS, 495, 1188
  * Myers et al. (2003) Myers S. T., et al., 2003, MNRAS, 341, 1
  * O’Brien et al. (2006) O’Brien T. J., et al., 2006, Nature, 442, 279
  * Padovani (2017) Padovani P., 2017, Nature Astronomy, 1, 0194
  * Petrov (2021) Petrov L., 2021, AJ, 161, 14
  * Porcas et al. (1979) Porcas R. W., Booth R. S., Browne I. W. A., Walsh D., Wilkinson P. N., 1979, Nature, 282, 385
  * Pradel et al. (2006) Pradel N., Charlot P., Lestrade J. F., 2006, A&A, 452, 1099
  * Prandoni et al. (2018) Prandoni I., Guglielmino G., Morganti R., Vaccari M., Maini A., Röttgering H. J. A., Jarvis M. J., Garrett M. A., 2018, MNRAS, 481, 4548
  * Radcliffe et al. (2018) Radcliffe J. F., et al., 2018, A&A, 619, A48
  * Rezaei et al. (2022) Rezaei S., McKean J. P., Biehl M., Javadpour A., 2022, MNRAS, 510, 5891
  * Shimwell et al. (2019) Shimwell T. W., et al., 2019, A&A, 622, A1
  * White et al. (1997) White R. L., Becker R. H., Helfand D. J., Gregg M. D., 1997, ApJ, 475, 479

##  Appendix A Hypothetical AVN locations

Table 3: The geocentric coordinates of a hypothetical African VLBI Network using 20 SKA-MID dishes. Geocentric coordinates |  |
---|---|---
$x$ | $y$ | $z$ | Station code | Country
(m) | (m) | (m) |  |
$5217030.16853176$ | $2534221.5108732083$ | $-2644582.0758971698$ | AVN1-BW | Botswana
$5277247.922210887$ | $2748624.8644503574$ | $-2289328.600036085$ | AVN2-BW | Botswana
$6347825.285561881$ | $-20717.88800793007$ | $618657.4567346348$ | AVN3-GH | Ghana
$6332560.529802635$ | $-178731.63840587594$ | $737359.7953479206$ | AVN4-GH | Ghana
$5104748.927307364$ | $3821230.688393553$ | $-142231.03817239657$ | AVN5-KE | Kenya
$5049932.002818037$ | $3895996.353310751$ | $5528.713103446548$ | AVN6-KE | Kenya
$4075522.7134878747$ | $4450989.035238711$ | $-2056922.0894522057$ | AVN7-MG | Madagascar
$4090030.2866659025$ | $4391134.402974519$ | $-2153824.1004977995$ | AVN8-MG | Madagascar
$3218155.6429566364$ | $5051492.497668295$ | $-2185029.212682493$ | AVN9-MU | Mauritius
$3213616.693183972$ | $5047594.365627962$ | $-2200557.8274866333$ | AVN10-MU | Mauritius
$4837312.3501758$ | $3090783.6593368473$ | $-2770624.9059733003$ | AVN11-MZ | Mozambique
$4841119.498520394$ | $3080179.155624241$ | $-2775745.617128389$ | AVN12-MZ | Mozambique
$5632553.87753738$ | $1731035.1721495497$ | $-2432816.736449252$ | AVN13-NAM | Namibia
$5688944.726890755$ | $1471992.094229581$ | $-2471560.196273584$ | AVN14-NAM | Namibia
$5415867.714182049$ | $2914108.406276784$ | $-1684593.4154433303$ | AVN15-ZM | Zambia
$5456235.737166995$ | $2979000.8406833992$ | $-1421987.8923449807$ | AVN16-ZM | Zambia
$5066809.931828817$ | $2715444.8969223998$ | $-2753782.8449464357$ | AVN17-ZA | South Africa
$5026325.59405273$ | $1674472.17663092$ | $-3539546.651210291$ | AVN18-ZA | South Africa
$5003701.377073657$ | $2463938.124620917$ | $-3083590.4930300973$ | AVN19-ZA | South Africa
$5054121.780624786$ | $2692033.063658318$ | $-2799397.734311496$ | AVN20-ZA | South Africa


# Enhancement of stability of metastable states in the presence of Lévy noise

A. A<EMAIL_ADDRESS>C<EMAIL_ADDRESS>B.
<EMAIL_ADDRESS>1Radiophysical Department,
Lobachevsky State University, Gagarin ave. 23, 603950 Nizhni Novgorod, Russia
2aDipartimento di Fisica “E.R. Caianiello”, Università di Salerno, Via
Giovanni Paolo II 132, Fisciano (SA), I-84084, Campania, Italy
2bINFN, Sezione di Napoli Gruppo Collegato di Salerno, Complesso Universitario
di Monte S. Angelo, Napoli, I-80126, Campania, Italy
3Dipartimento di Fisica e Chimica “E. Segrè”, Group of Interdisciplinary
Theoretical Physics, Università degli Studi di Palermo, Viale delle Scienze,
Edificio 18, I-90128 Palermo, Italy

###### Abstract

The barrier crossing event for superdiffusion in the form of symmetric Lévy
flights is investigated. We derive from the fractional Fokker-Planck equation
a general differential equation with the corresponding conditions useful to
calculate the mean residence time of a particle in a fixed interval for an
arbitrary smooth potential profile, in particular metastable, with a sink and
a Lévy noise with an arbitrary index $\alpha$. A closed expression in
quadrature of the nonlinear relaxation time for Lévy flights with the index
$\alpha=1$ in cubic metastable potential is obtained. Enhancement of the mean
residence time in the metastable state, analytically derived, due to Lévy
noise is found.

###### pacs:

05.40.Fb, 05.10.Gg, 02.50.Ey

The presence of anomalous diffusion can be explained as a deviation of the
statistics of fluctuations from the Gaussian law, giving rise to the
generalization of the Central Limit Theorem by Lévy and Gnedenko Lev37-Gne54 .
One of the manifestations of anomalous diffusion are Lévy flights. These are
stochastic processes characterized by the occurrence of extremely long jumps,
so that their trajectories are not continuous anymore. The length of these
jumps is distributed according to a Lévy stable statistics with a power law.
The relevance of the Lévy motion appears in many physical, natural, and social
complex systems. The Lévy-type statistics, in fact, is observed in various
scientific areas, where scale invariance phenomena take place, see Che06 ;
Dub08 ; Pad19  and references therein. The Lévy flights are a special class of
Markovian processes, therefore the Markovian analysis can be used to derive
the generalized Kolmogorov equation directly from Langevin equation with Lévy
noise Dub05 .
Metastability, as well as the transition process between metastable states, is
a ubiquitous phenomenon in nature affecting different fields of natural
sciences and advancing in its understanding it is a key challenge in complex
systems Fia05-Val18 ; Men13 ; Ser16 ; Gho20 ; Luc22 ; Hru18 ; Mor21 .
Experimental Nos09 ; Hru18  and theoretical Bia09 ; Mor21  results show that
long-lived metastable states, even if observed in different areas of physics,
were not explained.
Recently, the first passage time properties of Lévy flights for random search
strategies have been investigated Pal19 . Furthermore, noise-induced escape
from a metastable state in the presence of Lévy noise governs a plethora of
transition phenomena in complex systems, of physical, chemical and biological
nature, ranging from the motion of molecules to climate signals (see Che06 ;
Ser16 ; Luc22 ; Che07 ; Dub09 ; Dub20 ; Bau23 ; Jia23  and references
therein). The main focus in these papers is to understand how the barrier
crossing behavior in different potential profiles $U(x)$, is modified by the
presence of the Lévy-stable noise $L_{\alpha}(t)$ with index $\alpha$. This is
achieved by particle displacement analysis which obeys the following Langevin
equation

$\frac{dx}{dt}=-U^{\prime}\left(x\right)+L_{\alpha}\left(t\right),$ (1)

where $\alpha$ is the stability index of the Lévy distribution, with
$0<\alpha<2$. The main tools to investigate the barrier crossing problem for
Lévy flights in these above-mentioned papers are the first passage times and
residence times. Concerning these time characteristics, there are a lot of
numerical results and some analytical approximations Che06 ; Che07 ; Bau23 .
However, exact analytical results for barrier crossing problems in metastable
systems in the presence of Lévy noise are lacking. Obtaining exact analytical
results remains an open problem. This Letter aims to answer this question by
investigating the barrier crossing process in a system with a metastable state
driven by Lévy noise. Starting from the fractional Fokker-Planck equation
corresponding to Eq. (1), we investigate the barrier crossing event, by
focusing on the nonlinear relaxation time (NLRT). Specifically, we look at the
average time spent by the particle into the potential well (see Fig. 1). This
time is, in other words, the mean residence time in the metastable state,
which characterizes its stability.
Here, we address the following open questions: (_i_) the exact results of the
NLRT for a particle moving in an arbitrary smooth potential profile with a
sink and in the presence of Lévy noise with arbitrary index $\alpha$; (_ii_) a
closed expression in quadrature of the NLRT for Lévy flights with the index
$\alpha=1$ (Cauchy noise) in cubic metastable potential; (_iii_) the
enhancement of stability of metastable states, analytically derived, due to
Lévy noise.

The anomalous diffusion in the form of Lévy flights, for a particle moving in
a potential profile $U(x)$, is described by the following Fokker-Planck
equation with fractional space derivative

$\frac{\partial P}{\partial t}=\frac{\partial}{\partial
x}\left[U^{\prime}\left(x\right)P\right]+D_{\alpha}\frac{\partial^{\alpha}P}{\partial\left|x\right|^{\alpha}}\,,$
(2)

where $P\left(x,t|x_{0},0\right)$ is the transition probability density, and
$D_{\alpha}$ is the noise intensity parameter, in the sense that the size of a
cloud of particles undergoing Lévy motion increases with time as
$\left(D_{\alpha}t\right)^{1/\alpha}$. Equation (2) can be easily obtained
directly from Eq. (1) Dub05 .

The NLRT for a diffusion in a potential profile with a sink, see Fig. 1, is
defined as

$T\left(x_{0}\right)=\int_{0}^{\infty}\mathrm{Pr}\left(t,x_{0}\right)dt,$ (3)

where

$\mathrm{Pr}\left(t,x_{0}\right)=\int_{L_{1}}^{L_{2}}P\left(x,t|x_{0},0\right)dx$
(4)

represents the probability to find a particle in the interval
$\left(L_{1},L_{2}\right)$ at the time $t$, if it starts from some internal
point $x_{0}$. Substituting Eq. (4) into Eq. (3) and changing the order of
integration, we arrive at

$T\left(x_{0}\right)=\int_{L_{1}}^{L_{2}}Y\left(x,x_{0}\right)dx,$ (5)

where

$Y\left(x,x_{0}\right)=\int_{0}^{\infty}P\left(x,t|x_{0},0\right)dt.$ (6)

Integrating Eq. (2) with respect to $t$ from $0$ to $\infty$ and taking into
account the initial condition
$P\left(x,0|x_{0},0\right)=\delta\left(x-x_{0}\right)$ and the asymptotic
condition $P\left(x,\infty|x_{0},0\right)=0$ (for a potential with a sink), we
obtain the following equation for the function $Y\left(x,x_{0}\right)$

$\frac{d}{dx}\left[U^{\prime}\left(x\right)Y\right]+D_{\alpha}\frac{d^{\alpha}Y}{d\left|x\right|^{\alpha}}=-\delta\left(x-x_{0}\right).$
(7)

To solve Eq. (7) it is better to consider the Fourier transform of the
function $Y\left(x,x_{0}\right)$, i.e.,

$\widetilde{Y}\left(k,x_{0}\right)=\int_{-\infty}^{\infty}Y\left(x,x_{0}\right)\,e^{ikx}\,dx.$
(8)

For a smooth potential profiles $U\left(x\right)$, after Fourier transform,
Eq. (7) can be written in the differential form

$ik\,U^{\prime}\left(-i\frac{d}{dk}\right)\widetilde{Y}+D_{\alpha}\left|k\right|^{\alpha}\widetilde{Y}=e^{ikx_{0}}.$
(9)

It is convenient to introduce the new function, namely, the derivative of the
function $\widetilde{Y}\left(k,x_{0}\right)$ with respect to $x_{0}$

$G\left(k,x_{0}\right)=\frac{\partial}{\partial
x_{0}}\,\,\widetilde{Y}\left(k,x_{0}\right).$ (10)

Differentiating both parts of Eq. (9) with respect to $x_{0}$ we find

$U^{\prime}\left(-i\frac{d}{dk}\right)G-iD_{\alpha}\left|k\right|^{\alpha-1}\mathrm{sgn}\,k\,G=e^{ikx_{0}},$
(11)

where $\mathrm{sgn}\,x$ is the sign function.

Substituting $Y\left(x,x_{0}\right)$ from the backward Fourier transformation
into Eq. (5) and changing the order of integration, we have

$T\left(x_{0}\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\widetilde{Y}\left(k,x_{0}\right)\frac{e^{-ikL_{1}}-e^{-ikL_{2}}}{ik}\,\,dk.$
(12)

After differentiation of both sides of Eq. (12) with respect to $x_{0}$, in
accordance with Eq. (10), we find

$T^{\prime}\left(x_{0}\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}G\left(k,x_{0}\right)\frac{e^{-ikL_{1}}-e^{-ikL_{2}}}{ik}\,\,dk.$
(13)

One can easily check that, after replacing $k$ with $-k$, Eq. (11) coincides
with the equation for the complex conjugate function
$G^{*}\left(k,x_{0}\right)$, i.e.
$G\left(-k,x_{0}\right)=G^{*}\left(k,x_{0}\right)$. As a result, Eq. (13) can
be rearranged into a simpler form

$T^{\prime}\left(x_{0}\right)=\int_{0}^{\infty}\mathrm{Re}\left\\{G\left(k,x_{0}\right)\frac{e^{-ikL_{1}}-e^{-ikL_{2}}}{\pi
ik}\right\\}dk,$ (14)

where $\mathrm{Re}\left\\{...\right\\}$ denotes the real part of the
expression.

If a sink of the potential profile $U(x)$ is located at the point $x=\infty$
we have

$\lim_{x_{0}\rightarrow\infty}T\left(x_{0}\right)=0.$ (15)

After integrating Eq. (14) with respect to $x_{0}$ and taking into account the
condition (15), we find

$T\left(x_{0}\right)=\int_{x_{0}}^{\infty}\mathrm{Re}\left\\{\int_{0}^{\infty}G\left(k,z\right)\frac{e^{-ikL_{2}}-e^{-ikL_{1}}}{\pi
ik}\,\,dk\right\\}dz.$ (16)

Thus, it is sufficient to solve Eq. (11) only in the region $k>0$, obtaining

$U^{\prime}\left(-i\frac{d}{dk}\right)G-iD_{\alpha}k^{\alpha-1}G=e^{ikx_{0}}.$
(17)

Equations (16) and (17), which are among the main results of the paper, give
the exact relations useful to calculate the NLRT of the symmetric Lévy flights
with arbitrary index $\alpha$ in a smooth potential profile with a sink at
$x=\infty$. Now, we focus on a particle moving in a metastable cubic potential
profile (see Fig. 1)

Figure 1:  Cubic potential with metastable state, archetype model for any
metastable state. $L_{1}$ and $L_{2}$ are the interval boundaries, $x_{0}$ is
the initial condition, $a$ is a potential parameter.

$U\left(x\right)=-\frac{x^{3}}{3}+a^{2}x,$ (18)

and driven by a Cauchy-stable noise with Lévy index $\alpha=1$. Using Eq. (18)
and placing $\alpha=1$ in Eq. (11), we get

$\frac{d^{2}G}{dk^{2}}+\left(a^{2}-iD_{1}\,\mathrm{sgn}\,k\right)G=e^{ikx_{0}},$
(19)

which for $k>0$ becomes

$\frac{d^{2}G}{dk^{2}}+\left(a^{2}-iD_{1}\right)G=e^{ikx_{0}}.$ (20)

The general solution of Eq. (20) is the sum of the general solution of the
homogeneous equation and its particular solution. Under the condition of its
limitation, not divergent for arbitrary $k$, it takes the form

$G\left(k,x_{0}\right)=C\,e^{-\lambda
k}+\frac{e^{ikx_{0}}}{a^{2}-x_{0}^{2}-iD_{1}}\,,$ (21)

where $C$ is an unknown complex constant and $\lambda$ is one of the complex
roots

$z=\pm\sqrt{iD_{1}-a^{2}},$

having a positive real part, $\lambda=\lambda_{1}+i\lambda_{2}$, where

$\displaystyle\lambda_{1}=\left(a^{4}+D_{1}^{2}\right)^{1/4}\sin\left[\frac{1}{2}\,\arctan\left(\frac{D_{1}}{a^{2}}\right)\right],$
$\displaystyle\lambda_{2}=\left(a^{4}+D_{1}^{2}\right)^{1/4}\cos\left[\frac{1}{2}\,\arctan\left(\frac{D_{1}}{a^{2}}\right)\right].$
(22)

To find the unknown constant $C$ we use the continuity conditions for the
function $G\left(k,x_{0}\right)$ and its first derivative at the point $k=0$

$\displaystyle\lim_{k\rightarrow
0^{+}}G\left(k,x_{0}\right)=\lim_{k\rightarrow 0^{-}}G\left(k,x_{0}\right),$
$\displaystyle\lim_{k\rightarrow
0^{+}}\frac{dG\left(k,x_{0}\right)}{dk}=\lim_{k\rightarrow
0^{-}}\frac{dG\left(k,x_{0}\right)}{dk}\,.$ (23)

For $k<0$, Eq. (19) transforms into

$\frac{d^{2}G}{dk^{2}}+\left(a^{2}+iD_{1}\right)G=e^{ikx_{0}},$ (24)

and its solution under the condition of its limitation reads

$G\left(k,x_{0}\right)=C^{*}e^{\lambda^{*}k}+\frac{e^{ikx_{0}}}{a^{2}-x_{0}^{2}+iD_{1}}\,.$
(25)

Using Eqs. (21), (25) and conditions (23) we get

$\displaystyle C+\frac{1}{a^{2}-x_{0}^{2}-iD_{1}}$ $\displaystyle=$
$\displaystyle C^{*}+\frac{1}{a^{2}-x_{0}^{2}+iD_{1}}\,$
$\displaystyle-\lambda C+\frac{ix_{0}}{a^{2}-x_{0}^{2}-iD_{1}}$
$\displaystyle=$
$\displaystyle\lambda^{*}C^{*}+\frac{ix_{0}}{a^{2}-x_{0}^{2}+iD_{1}}.\,\,$
(26)

The final expression for the constant $C$, obtained from the system (26),
reads

$C=-\frac{D_{1}\left(x_{0}+\lambda_{2}+i\lambda_{1}\right)}{\lambda_{1}[\left(a^{2}-x_{0}^{2}\right)^{2}+D_{1}^{2}]}\,.$
(27)

Substituting Eqs. (21) and (27) into Eq. (16) and calculating the internal
integral we arrive at

$\displaystyle T\left(x_{0}\right)$ $\displaystyle=$
$\displaystyle\frac{D_{1}}{\pi}\int_{x_{0}}^{\infty}\left(A\frac{z+\lambda_{2}}{\lambda_{1}}+B\right)\frac{dz}{\left(z^{2}-a^{2}\right)^{2}+D_{1}^{2}}$
(28) $\displaystyle+$
$\displaystyle\frac{D_{1}}{\pi}\int_{x_{0}}^{\infty}\ln\left|\frac{z-L_{1}}{z-L_{2}}\right|\frac{dz}{\left(z^{2}-a^{2}\right)^{2}+D_{1}^{2}}$
$\displaystyle+$
$\displaystyle\int_{x_{0}}^{L_{2}}\frac{\left(z^{2}-a^{2}\right)dz}{\left(z^{2}-a^{2}\right)^{2}+D_{1}^{2}}\,,$

where

$\displaystyle A$ $\displaystyle=$
$\displaystyle\arctan\frac{\lambda_{2}+L_{2}}{\lambda_{1}}-\arctan\frac{\lambda_{2}+L_{1}}{\lambda_{1}}\,,$
$\displaystyle B$ $\displaystyle=$
$\displaystyle\frac{1}{2}\ln\frac{\lambda_{1}^{2}+\left(L_{2}+\lambda_{2}\right)^{2}}{\lambda_{1}^{2}+\left(L_{1}+\lambda_{2}\right)^{2}}\,.$
(29)

The exact quadrature formula of Eq. (28) is the other main result of the
paper. For unstable initial conditions of the particle just beyond the
potential barrier and within the interval $a<x_{0}<L_{2}$, the normalized mean
residence time in the metastable state $\left<T(x_{0})\right>/T_{d}(x_{0})$ as
a function of the noise intensity parameter $D_{1}$ has a nonmonotonic
behaviour with a maximum (see curves in Figs. 2 and 3). This is the noise
enhanced stability (NES) phenomenon, already investigated with Gaussian noise
sources Fia05-Val18 .

$T_{d}(x_{0})$ is the dynamical time obtained by setting $D_{\alpha}=0$ in Eq.
(28)

Figure 2:  The results of analytical calculations according to Eq. (28) are
shown. (a) Normalized NLRT $\left<T(x_{0})\right>/T_{d}(x_{0})$ for a
metastable cubic potential as a function of the noise intensity parameter
$D_{1}$ for different positions $L_{1}$ of the left boundary ranging from $0$
to $-6$ with steps $0.2$. (b) Density plot of
$\left<T(x_{0})\right>/T_{d}(x_{0})$ versus $L_{1}$ and $D_{1}$. The white
dotted line marks the position of the NES maxima. The parameter values are:
$x_{0}=3.0$, $L_{2}=4$, $a=1$.

$T_{d}=\int_{x_{0}}^{L_{2}}\frac{dz}{z^{2}-a^{2}}\,.$ (30)

For $x_{0}<a<L_{2}$ the integral in Eq. (30) diverges, which means the
impossibility for a particle to cross the potential barrier, located at the
point $x=a$, in the absence of driving noise. For $a<x_{0}<L_{2}$ we obtain
the finite dynamical time

$T_{d}(x_{0})=\frac{1}{2a}\,\ln\frac{\left(L_{2}-a\right)\left(x_{0}+a\right)}{\left(L_{2}+a\right)\left(x_{0}-a\right)}\,.$
(31)

The normalized NLRT $\left<T(x_{0})\right>/T_{d}(x_{0})$ for a metastable
cubic potential as a function of the noise intensity parameter $D_{1}$ for
different positions $L_{1}$ of the left boundary and a fixed value of the
right boundary $L_{2}=4$ in a semilog plot is shown in Fig. 2. The different
values of $L_{1}$ range from $0$ to $-6$ with steps $0.2$. A nonmonotonic
behavior of the normalized NRLT with a maximum as a function of the noise
intensity is observed for all values of $L_{1}$ analyzed, that is the particle
is temporarily trapped in the metastable state.

Figure 3:  The results of analytical calculations according to Eq. (28) are
shown. (a) NLRT $\left<T(x_{0})\right>$ for a metastable cubic potential as a
function of the noise intensity parameter $D_{1}$ for different positions
$L_{2}$ of the right boundary ranging from $3$ to $7$ with steps $0.2$. (b)
Density plot of $\left<T(x_{0})\right>$ versus $L_{2}$ and $D_{1}$. The white
dotted line marks the position of the NES maxima. The parameter values are:
$x_{0}=3.0$, $L_{1}=-4$, $a=1$.

In Fig. 2b a density plot of $\left<T(x_{0})\right>/T_{d}(x_{0})$ versus
$L_{1}$ and $D_{1}$ is shown. The white dotted line marks the position of the
NES maxima.

The maxima and all curves increase as the value of the left limit $L_{1}$
decreases. This gives rise to an increasing size of the basin of attraction of
the metastable state Ser16 , responsible for the increase in the normalized
NRLT. Furthermore, in the limit $D_{1}\rightarrow 0$ there are three different
asymptotic values of the normalized NRLT, the value of which increases with
the size of the basin of attraction when $L_{1}$ varies from $0$ to $-6$. We
note that in the limit $D_{1}\rightarrow 0$ and for unstable initial
conditions of the particle there is, with a Gaussian noise source, a divergent
behavior of $T(x_{0})$ Fia05-Val18 . For Lévy flights, instead, we obtain for
$T(x_{0})$ a finite nonmonotonic behavior as a function of the noise intensity
parameter $D_{1}$. Because of heavy tails, a particle spends a finite time in
the metastability area even in the limit $D_{1}\rightarrow 0$. For very large
noise intensity, in the limit $D_{1}\rightarrow\infty$, the normalized NRLT
has a power law behavior as a function of the noise intensity Che06 ; Dub08 .
In Fig. 3, the NRLT $\left<T(x_{0})\right>$ versus $D_{1}$ in a semilog plot
for different positions $L_{2}$ of the right boundary at a fixed value of the
left boundary $L_{1}=-4$ is shown. The different values of $L_{2}$ range from
$3$ to $7$ with steps $0.2$. A nonmonotonic behavior of
$\left<T(x_{0})\right>$ versus $D_{1}$ for all values of $L_{2}$ investigated
is observed. In Fig. 3b a density plot of $\left<T(x_{0})\right>$ versus
$L_{2}$ and $D_{1}$ is shown. The white dotted line marks the position of the
NES maxima.
The asymptotic behaviors shown in Figs. 2 and 3 reproduce the asymptotic
expressions of Eq. (28) in the limits $D_{1}\rightarrow 0$ and
$D_{1}\rightarrow\infty$. In particular, for $D_{1}\rightarrow 0$ we have Sup

$T\left(x_{0}\right)\simeq
C_{1}\left[\frac{a}{x_{0}-a}+\frac{1}{2}\ln\frac{x_{0}-a}{x_{0}+a}\right]+T_{d}(x_{0})\,,$
(32)

where

$\displaystyle C_{1}\simeq\left\\{\begin{array}[]{ll}1,&L_{1}<-a,\\\
1/2,&L_{1}=-a,\\\ 0,&L_{1}>-a,\end{array}\right.$ (36)

giving rise to three different asymptotic values of the NLRT for
$D_{1}\rightarrow 0$ in the case of fixed $L_{2}$ (see Fig. 2) and different
asymptotic values depending on the different values of $L_{2}$ in the case of
fixed $L_{1}$ (see Fig. 3). And for $D_{1}\rightarrow\infty$ Sup  we have

$T\left(x_{0}\right)\sim\frac{1}{D_{1}}\,,$ (37)

that is a power law behavior of the NLRT as a function of the noise intensity
(see Figs. 2 and 3).

In conclusion, we obtain the general equations useful to calculate the NLRT
for superdiffusion in the form of symmetric Lévy flights, for an arbitrary
Lévy index $\alpha$ and an arbitrary smooth potential profile with a sink. For
a Cauchy driven noise ($\alpha=1$) we find the closed expression in
quadratures of the NLRT as a function of the noise intensity parameter, the
initial position, and the parameters of the potential. The interplay between
trapping in the metastable state, at small noise intensities, and long jumps
of Lévy flights produces a finite nonmonotonic enhancement of the mean
residence time in the metastable state. Our general equations provide a useful
tool to describe different dynamical behaviors in complex systems,
characterized by anomalous diffusion and non-exponential relaxation phenomena,
such as spatial extended systems Jia23 .

###### Acknowledgements.

This work was partially supported by Italian Ministry of University and
Research (MUR).

## References

  * (1) P. Lévy, _Theory de l’addition des variables Aléatoires_ (Gauthier–Villars, Paris, 1937); B. V. Gnedenko and A. N. Kolmogorov, _Limit Distributions for Sums of Independent Random Variables_ (Addison–Wesley, Cambridge, 1954).
  * (2) A. V. Chechkin, V. Yu. Gonchar, J. Klafter, and R. Metzler, Adv. Chem. Phys. 133, 439 (2006).
  * (3) A. A. Dubkov, B. Spagnolo, and V. V. Uchaikin, Int. J. Bifur. Chaos 18, 2649 (2008).
  * (4) A. Padash, A. V. Chechkin, B. Dybiec, I Pavlyukevich, B. Shokri, R. Metzler, J. Phys. A: Math. Theor. 52, 454004 (2019).
  * (5) A. Dubkov and B. Spagnolo, Fluct. Noise Lett. 5, L267 (2005).
  * (6) A. Fiasconaro, B. Spagnolo, and S. Boccaletti , Phys. Rev. E 72, 061110 (2005); D. Valenti, A. Carollo, B. Spagnolo, Phys. Rev. A 97, 042109 (2018); I.A. Surazhevsky, V.A. Demin _et al._ , Chaos, Soliton. Fract. 146, 110890 (2021); C. Guarcello Chaos, Soliton. Fract. 153, 111531 (2021).
  * (7) P. J. Menck, J. Heitzig, N. Marwan, J. Kurths, Nat. Phys. 9, 89–92 (2013).
  * (8) L. Serdukova, Y. Zheng, J. Duan, J. Kurths, Chaos 26, 073117 (2016); _ibidem_ Sci. Rep. 7, 9336 (2017).
  * (9) D. K. Ghosh, A. Ranjan, Protein. Science 29, 1559–1568 (2020).
  * (10) V. Lucarini, L. Serdukova, G. Margazoglou, Nonlinear Proc. Geoph. 29, 183-205 (2022).
  * (11) L. Hruby, N. Dogra, M. Landini, T. Donner, T. Esslinger, Proc. Natl. Acad. Sci. USA 115, 3279-3284 (2018); G. Zeng, J. Gao _et al._ , Proc. Natl. Acad. Sci. USA 117, 17528-17534 (2020); J. Ravnik, M. Diego _et al._ , Nat. Commun. 12, 2323 (2021).
  * (12) T. Mori, Phys. Rev. Res. 3, 043137 (2021); X.-D. Zheng, C. Li, S. Lessard, Y. Tao, Phys. Rev. Lett. 120, 218101 (2018).
  * (13) V. Nosenko, S. K. Zhdanov, A. V. Ivlev, C. A. Knapek, and G. E. Morfill, Phys. Rev. Lett. 103, 015001 (2009); J. J. L. Morton, A. Tiwari, G. Dantelle, K. Porfyrakis, A. Ardavan, G. A. Briggs, Phys. Rev. Lett. 101, 013002 (2008); S. Hayami, D. Urakami, Y. Kojima _et al._ , Inorg. Chem. 49, 1428–1432 (2010); B. Tu, R. Si, Y. Shen, J. Wang, B. Wei, C. Chen, K. Yao, Y. Zou _et al._ , Phys. Rev. Res. 5, 043014 (2023).
  * (14) I. Biazzo, F. Caltagirone, G. Parisi, and F. Zamponi, Phys. Rev. Lett. 102, 195701 (2009); R. Giachetti, E. Sorace, Phys. Rev. Lett. 101, 190401 (2008); L. M. Krauss and J. Dent, Phys. Rev. Lett. 100, 171301 (2008); Astron. Astrophys 432, 117–138 (2005).
  * (15) V. V. Palyulin, G. Blackburn _et al._ , New J. Phys. 21, 103028 (2019); A. Padash, T. Sandev _et al._ , Fractal Fract. 6, 260 (2022); M. R. Evans, S. N. Majumdar, G. Schehr, J. Phys. A: Math. Theor. 53, 193001 (2020).
  * (16) A. V. Chechkin, O. Yu. Sliusarenko, R. Metzler, and J. Klafter, Phys. Rev. E 75, 041101 (2007).
  * (17) A. A. Dubkov, A. La Cognata, and B. Spagnolo, J. Stat. Mech-Theory E., P01002 (2009).
  * (18) A. A. Dubkov, B. Dybiec, B. Spagnolo, A. Kharcheva, C. Guarcello, D. Valenti, Phys. Rev. E 102, 042142 (2020).
  * (19) A. Baule, P. Sollich, Sci. Rep. 13, 3853 (2023).
  * (20) C. Jiao and G. A. Gottwald, Proc. Roy. Soc. A-Math. Phys. 479, 20230349 (2023).
  * (21) See Supplemental Material.

Supplemental Material
Enhancement of stability of metastable states in the presence of Lévy noise

##  S-I Investigation of the NLRT in the limits $D_{1}\rightarrow 0$ and
$D_{1}\rightarrow\infty$

To investigate the NLRT in the asymptotic limits $D_{1}\rightarrow 0$ and
$D_{1}\rightarrow\infty$, we start from the expressions of the parameters
$\lambda_{1}$ and $\lambda_{2}$ of Eq. (22) of the main text

$\displaystyle\lambda_{1}=\left(a^{4}+D_{1}^{2}\right)^{1/4}\sin\left[\frac{1}{2}\,\arctan\left(\frac{D_{1}}{a^{2}}\right)\right],$
(S1)
$\displaystyle\lambda_{2}=\left(a^{4}+D_{1}^{2}\right)^{1/4}\cos\left[\frac{1}{2}\,\arctan\left(\frac{D_{1}}{a^{2}}\right)\right],$

and, using trigonometry formulas, we rewrite Eq. (S1) in a simpler form

$\displaystyle\lambda_{1}=\frac{a}{\sqrt{2}}\,\sqrt{\sqrt{1+\frac{D_{1}^{2}}{a^{4}}}-1}\,,$
$\displaystyle\lambda_{2}=\frac{a}{\sqrt{2}}\,\sqrt{\sqrt{1+\frac{D_{1}^{2}}{a^{4}}}+1}\,.$
(S2)

1\. Asymptotics for $D_{1}\rightarrow 0$

For small values of $D_{1}$, we can use the approximate expansion:
$\sqrt{1+x}\simeq 1+x/2-x^{2}/8$, where $x=D_{1}^{2}/a^{4}\ll 1$. As a result,
we obtain

$\lambda_{1}\simeq\frac{D_{1}}{2a}\left(1-\frac{D_{1}^{2}}{8a^{4}}\right),\,\,\lambda_{2}\simeq
a\left(1+\frac{D_{1}^{2}}{8a^{4}}\right).$ (S3)

First of all, it is better to write the expression for $A$ in another form.
Using the well-known relation

$\displaystyle\arctan{\frac{1}{x}}=\frac{\pi}{2}-\mathrm{arccot}{\frac{1}{x}}=\frac{\pi}{2}-\arctan{x},$

we can rewrite expressions (29) of the main text for $A$ and $B$ in the
following form

$\displaystyle
A=\arctan\frac{\lambda_{1}}{\lambda_{2}+L_{1}}-\arctan\frac{\lambda_{1}}{\lambda_{2}+L_{2}}+\pi\cdot\,1\left(-\lambda_{2}-L_{1}\right)\,,$
$\displaystyle
B=\frac{1}{2}\ln\frac{\lambda_{1}^{2}+\left(L_{2}+\lambda_{2}\right)^{2}}{\lambda_{1}^{2}+\left(L_{1}+\lambda_{2}\right)^{2}}\,,$
(S4)

where $1(x)$ is the step function.

Substituting Eq. (S3) into Eq. (S4) for $A$ and $B$ we have

$\displaystyle
A\simeq\left\\{\begin{array}[]{ll}\pi+D_{1}(L_{2}-L_{1})/\left[2a(a+L_{1})(a+L_{2})\right],&L_{1}<-a,\\\
\pi/2-D_{1}(3a+L_{2})/\left[4a^{2}(a+L_{2})\right],&L_{1}=-a,\\\
D_{1}(L_{2}-L_{1})/\left[2a(a+L_{1})(a+L_{2})\right],&L_{1}>-a.\end{array}\right.$
(S8) $\displaystyle
B\simeq\left\\{\begin{array}[]{ll}\ln\left[\left(a+L_{2}\right)/\left|a+L_{1}\right|\right],&L_{1}\neq-a,\\\
\ln\left[2a\left(a+L_{2}\right)/D_{1}\right],&L_{1}=-a.\end{array}\right.$
(S11)

As seen from Eq. (S8), the value $A$ does not go to zero in the limit
$D_{1}\rightarrow 0$ in the case $L_{1}\leq-a$.

Substitution of Eqs. (S8) and (S11) into Eq. (28) of the main text gives in
the limit $D_{1}\rightarrow 0$

$T\left(x_{0}\right)\simeq\frac{2aC_{1}}{\pi}\int_{x_{0}}^{\infty}\frac{dz}{\left(z-a\right)^{2}(z+a)}+T_{d}(x_{0})\,,$
(S12)

where

$\displaystyle C_{1}\simeq\left\\{\begin{array}[]{ll}\pi,&L_{1}<-a,\\\
\pi/2,&L_{1}=-a,\\\ 0,&L_{1}>-a\end{array}\right.$ (S16)

and $T_{d}(x_{0})$ is the dynamical time.

The integral in Eq. (S12) can be calculated in the analytical form. As a
result, we obtain in the limit $D_{1}\rightarrow 0$

$T\left(x_{0}\right)\simeq
C_{1}\left[\frac{a}{x_{0}-a}+\frac{1}{2}\ln\frac{x_{0}-a}{x_{0}+a}\right]+\frac{1}{2a}\,\ln\frac{\left(L_{2}-a\right)\left(x_{0}+a\right)}{\left(L_{2}+a\right)\left(x_{0}-a\right)}\,.$
(S17)

This gives rise to three different asymptotic values of the NLRT for
$D_{1}\rightarrow 0$ in the case of fixed $L_{2}$ (see Fig. 2 of the main
text) and different asymptotic values depending on the different values of
$L_{2}$ in the case of fixed $L_{1}$ (see Fig. 3 of the main text).

2\. Asymptotics for $D_{1}\rightarrow\infty$

Now we consider the case of very large $D_{1}$. From Eq. (S2) we easily find

$\lambda_{1}\simeq\sqrt{\frac{D_{1}}{2}}\,,\quad\lambda_{2}\simeq\sqrt{\frac{D_{1}}{2}}\,.$
(S18)

Substituting Eq. (S18) into Eq. (29) of the main text we obtain the following
approximate expressions for the constants $A$ and $B$

$A=B\simeq\frac{L_{2}-L_{1}}{\sqrt{2D_{1}}}.$ (S19)

Substitution of Eqs. (S18) and (S19) into Eq. (28) of the main text gives

$\displaystyle
T\left(x_{0}\right)\simeq\frac{L_{2}-L_{1}}{\pi}\int_{x_{0}}^{\infty}\frac{(z+\sqrt{2D_{1}})\,dz}{\left(z^{2}-a^{2}\right)^{2}+D_{1}^{2}}$
(S20)
$\displaystyle+\frac{D_{1}}{\pi}\int_{x_{0}}^{\infty}\ln\left|\frac{z-L_{1}}{z-L_{2}}\right|\frac{dz}{\left(z^{2}-a^{2}\right)^{2}+D_{1}^{2}}+\int_{x_{0}}^{L_{2}}\frac{\left(z^{2}-a^{2}\right)dz}{\left(z^{2}-a^{2}\right)^{2}+D_{1}^{2}}\,.$

The first integral in Eq. (S20) can be calculated in analytical form and for
the large $D_{1}$ gives

$\frac{L_{2}-L_{1}}{\pi}\int_{x_{0}}^{\infty}\frac{(z+\sqrt{2D_{1}})\,dz}{\left(z^{2}-a^{2}\right)^{2}+D_{1}^{2}}\simeq\frac{3\left(L_{2}-L_{1}\right)}{4D_{1}}\,.$
(S21)

The second integral in Eq. (S20) can be estimated for large $D_{1}$ using the
mean value theorem for a definite integral, namely

$\displaystyle\frac{D_{1}}{\pi}\int_{x_{0}}^{\infty}\ln\left|\frac{z-L_{1}}{z-L_{2}}\right|\frac{dz}{\left(z^{2}-a^{2}\right)^{2}+D_{1}^{2}}\simeq$
$\displaystyle\frac{D_{1}}{\pi}\ln\left(\frac{\sqrt{D_{1}}-L_{1}}{\sqrt{D_{1}}-L_{2}}\right)\int_{x_{0}}^{\infty}\frac{dz}{\left(z^{2}-a^{2}\right)^{2}+D_{1}^{2}}\simeq$
(S22)
$\displaystyle\frac{(L_{2}-L_{1})\sqrt{D_{1}}}{\pi}\int_{x_{0}}^{\infty}\frac{dz}{\left(z^{2}-a^{2}\right)^{2}+D_{1}^{2}}\simeq\frac{(L_{2}-L_{1})\sqrt{2}}{4D_{1}}\sim\frac{1}{D_{1}}\,.$

The last integral in Eq. (S20), due to the finite limits, can be easily
estimated

$\int_{x_{0}}^{L_{2}}\frac{\left(z^{2}-a^{2}\right)dz}{\left(z^{2}-a^{2}\right)^{2}+D_{1}^{2}}\simeq\frac{1}{D_{1}^{2}}\left[\frac{L_{2}^{3}-x_{0}^{3}}{3}-a^{2}(L_{2}-x_{0})\right]\sim\frac{1}{D_{1}^{2}}\,.$
(S23)

Taking into account Eqs. (S21), (S22), and (S23), we find finally for large
$D_{1}$

$T\left(x_{0}\right)\sim\frac{1}{D_{1}}\,,$ (S24)

that is a power law behavior in agreement with previous investigations (see
Refs. [2,3] of the main text).