# Acoustic Neutrino Detection In a Adriatic Multidisciplinary Observatory
(ANDIAMO)

Antonio Marinelli<EMAIL_ADDRESS>Pasquale Migliozzi
<EMAIL_ADDRESS>Andreino Simonelli
<EMAIL_ADDRESS>INFN - Sezione di Napoli, Complesso Univ. Monte
S. Angelo, I-80126 Napoli, Italy

###### Abstract

The existence of cosmic accelerators able to emit charged particles up to EeV
energies has been confirmed by the observations made in the last years by
experiments such as Auger and Telescope Array. The interaction of such
energetic cosmic-rays with gas or low energy photons, surrounding the
astrophysical sources or present in the intergalactic medium, guarantee an
ultra-high-energy neutrino related emission. When these energetic neutrinos
interact in a medium produce a thermo-acoustic process where the energy of
generated particle cascades can be conveyed in a pressure pulse propagating
into the same medium. The kilometric attenuation length as well as the well-
defined shape of the expected pulse suggest a large-area-undersea-array of
acoustic sensors as an ideal observatory. For this scope, we propose to
exploit the existing and no more operative offshore (oil rigs) powered
platforms in the Adriatic sea as the main infrastructure to build an acoustic
submarine array of dedicated hydrophones covering a surface area up to 10000
Km2 and a volume up to 500 Km3. In this work we describe the advantages of
this detector concept using a ray tracing technique as well as the scientific
goals linked to the challenging purpose of observing for the first time ultra-
high-energy cosmic neutrinos. This observatory will be complementary to the
dedicated radio array detectors with the advantages of avoiding any possible
thermo-acoustic noise from the atmospheric muons.

###### keywords:

neutrinos, radiation mechanism: non thermal, instrumentation: detectors ,
methods: analytical , techniques: underwater acoustic

††journal: Astroparticle Physics

##  1 Introduction

The ultra-high-energy cosmic-rays (UHECRs) observed by Pierre Auger apparatus
[1] and Telescope Array [2] suggest the presence of astrophysical accelerators
who can originate EeV neutrinos considering the possible cosmic-rays (CRs)
interaction [3] with gas or photons in the source environment [4] or along the
path to the Earth [5]. On the other hand, the astrophysical flux observed by
IceCube telescope during the last decade [6, 7] demonstrates that the former
interactions occur with CRs accelerated at least up to an energy of hundred of
PeVs [8]. Who are the main responsibles for these neutrino emissions is still
a matter of debate. However, while PeV-EeV neutrinos can hardly be produced in
our Galaxy we can expect the extragalactic accelerators to be the main
candidates for that emission. The muonic neutrino at 290 TeV observed by
IceCube on 22 of September 2017 from a direction compatible with the blazar
TXS 0506+056, in coincidence with gamma-ray flare observed by Fermi-LAT and
MAGIC telescopes, makes blazars to climb the charts of the very-high-energy
neutrino emitter candidates. Most of the hadronic models related to this class
of powerful sources suggest a favourite range of neutrino production between
several hundreds of TeV up to a few EeVs. This implies that the astrophysical
samples of IceCube events can be only partially explained by the blazar
emission. This hypothesis can be confirmed whenever we will have ultra-high-
energy (UHE) neutrino observations. Another important neutrino emitter
candidate who is expected to emit up to EeV energies is represented by the
gamma-ray burst (GRB) who are able to accelerate CRs up to these energies [9]
during the so-called prompt phase [10]. In addition to the UHE neutrinos that
could arrive directly from one of the mentioned astrophysical sources, at
energies above PeV, there is also the cosmogenic neutrino flux produced by the
interaction UHECRs with the astrophysical background radiation fields. This
photon target is represented by the extra-galactic background light (EBL) and
the cosmic microwave backgrounds (CMB) along the path between the source and
the Earth [11]. While the former signal is expected to produce a spotted
observation in correlation with UHE sources, the latter is expected to be a
uniform signal over the full sky. The possibility to identify single
astrophysical sources of UHE neutrinos can be challenging for such a detector
due to the number of events needed. The information about the astrophysical
spectral shape observed will be fundamental to identify a class of
astrophysical emitters in this range of energy.
The detection of neutrinos from PeV to EeV range could represent a major
breakthrough and go beyond the capabilities of Cherenkov neutrino telescopes
like IceCube [12], Antares [13], Baikal [14] and KM3NeT [15]. A larger
instrumented surface and a different approach to collect the neutrino debris
signature are needed. The two main techniques proposed for the neutrino
detection at UHEs are the radio detection in atmosphere and the acoustic
detection in water/ice thanks to the Askarian effect [16] happening when high-
energy ionizing particles pass through a dense medium. The radio detection
techniques improved during the last decades thanks to the installation of
several prototypes and telescopes as the pioneering RICE [17, 18], ARA [19] ,
ARIANNA [20, 21], RNO-G [22], GRAND [23] as well as the balloon experiments
like ANITA [24]. Even if up-to-date none of these radio arrays were able to
find an excess of signal statistically significant as recently reported by ARA
collaboration [25] the analyses done helped to better understand the time
dependent behaviour of the radio arrays as well as to improve the rejection of
noise.
The use of acoustic techniques for UHE neutrino telescope favours the
exploiting of natural water/ice reservoirs like seas and big lakes [26, 27].
Thanks to the accelerated expansion of the cylindrical volume heated by the
interaction of neutrino with a nucleon and the production of a hadronic
cascade [28, 29] an acoustic array can reconstruct the generating neutrino
event. The two experimental setups represents a complementary approach to
cover a larger solid angle, while the first favours the Earth-skimming
neutrino events, the second privileges the down-going events.
Despite the test setups done with small acoustic arrays like SPATS [30],
O$\nu$DE5[31], ACoRNE [32], AMADEUS [33] and SOUND [34] already installed for
monitoring activities or built as sub-detector of major Cherenkov telescopes,
at different latitudes, a large scale underwater/ice acoustic telescope has
not yet finalized. For a dedicated acoustic array, a large area coverage
$\mathcal{O}(1000\leavevmode\nobreak\ $Km${}^{2})$ and a correspondent
sizeable volume $\mathcal{O}(100\leavevmode\nobreak\ $Km${}^{3})$ [26] are
needed, even considering a sparse units distribution.
In this work we introduce the possibility of exploiting the ENI not operative
powered oil rigs in the offshore of the Adriatic sea. The oil upstream
activity is terminated but the infrastructures are still available for
scientific purposes. We show that the Mediterranean Sea certainly represents a
preferential environment for an acoustic telescope due to the effects of water
temperature, currents, and salinity [14]. Moreover, we prove that by
exploiting the shape and the propagation pattern of the acoustic wave-field
generated by UHE neutrino cascade we can cover an array of unprecedented area.
This description can be translated to a large-scale distribution of small
sized strings equipped with equidistant acoustic sensors. In Fig. 1 we report
the possible sensitivity of an extended acoustic array simulation [35], made
by the AMADEUS group, in comparison with the sensitivity of mentioned radio-
arrays and the major expected diffuse UHE neutrino signals. It is important to
say that the acoustic array proposed here, ANDIAMO (Acoustic Neutrino
Detection in a Multidisciplinary Observatory), can be used not only for UHE
neutrino detection but also for marine biology and geophysical studies. The
large amount of data that will be collected over a wide span of frequencies
can be used in full by many scientific communities. Monitoring the
oceanographic parameters of a large portion of the Adriatic sea is an
effective way to study climate changes, and the extremely high sensitivity of
the acoustic sensors is suitable for seismic monitoring of local and
teleseismic events. The shallow waters of the site permits the installation of
ocean bottom seismometer (OBS) platforms at low cost. All the equipment
installed for astroparticle physics purposes is needed for a better
understanding of the Sea in terms of medium where our expected signal
propagates. At the same time, the collected signals will serve a vast branch
of science with a large amount of data.

Figure 1: In this plot are reported the expected sensitivities of the main
radio array detectors for UHE neutrinos, comprising the future IceCube
Gen2-Radio [36], in comparison with the capabilities of an acoustic array and
the main UHE diffuse neutrino expected signals. In details, with the blu stars
it is reported the sensitivity of one of the few large scale undersea acoustic
detectors recently simulated, in this case by the AMADEUS group [35]. This
simulated acoustic array refers to a huge sparse array (1500 Km3) of sensors
with a density of $\sim 100$ per Km3 and a kinetic threshold of 5 mPa.

##  2 The Science case

The study of possible neutrino emission beyond the energies already observed
by IceCube telescope can answer important questions about the processes
happening on astrophysical accelerators, as well as the nature of local Dark
Matter (DM) scenarios. Even though the identifications of single UHE
accelerators, like the ones reported in Fig. 2 and Fig. 3 can be challenging,
the observation of that kind of neutrino signal and the estimation of an
isotropic flux represent in itself a major step for the astroparticle
community.

###  2.1 Cosmogenic neutrinos

The ankle observed in the cosmic-ray spectral energy distribution at around
$5\times 10^{18}$ eV is interpreted as the transition point between the
Galactic and the extragalactic origin of the observed comic rays [37, 38].
With the observations made by the Pierre Auger experiment [1] and Telescope
Array [2] we already collected hundreds of these events with energy above
$10^{18}$ eV. When these events interact with background photons, like CMB
[39] and EBL [40], they are producing gamma-rays and neutrinos through the
decay of neutral and charged pions ($\pi^{0}\rightarrow\gamma\gamma$,
$\pi^{+}\rightarrow e^{+}\nu_{e}\nu_{\mu}\overline{\nu}_{\mu}$) who are
produced trough the $\Delta^{+}$ channel as described in the following
equation:

$p+\gamma_{bg}\rightarrow\Delta^{+}\rightarrow\begin{cases}&p+\pi^{0},\\\
&n+\pi^{+}\end{cases}$ (1)

The energy threshold of this process is of 1.08 GeV in the reference system of
the interacting particles [41]. This is also the main reason of the distortion
of the proton spectrum above $3\times 10^{19}$ eV during propagation, known as
the Greisen–Zatsepin–Kuzmin cutoff ([42], [43]).
Considering that the CMB density increases with the distance as $(1+z)^{3}$ we
should expect also the cosmogenic neutrino production depending on the
redshift considered. An equivalent dependency should be expected for UHE
neutrino produced by the interaction of cosmic rays with EBL event though the
EBL parameters are less known in respect to the CMB, due to the spectral
evolution of optical, infrared and ultraviolet. In this case, being the EBL
photons more energetic that the CMB ones, the energy threshold for neutrino
production through the photopion process becomes lower.
Beta-decay can contribute to the cosmogenic neutrino production as well, due
to the decay of neutron obtained with the charged pion production:

$n\rightarrow p+e^{-}+\overline{\nu}_{e}$ (2)

Heavier nuclei with $Z>1$ can also produce beta-decay process, adding a
possible contribution to the UHE neutrino flux generated through photopion
interaction. Following the presented description, the cosmogenic neutrinos
originated by the interaction of UHECRs with thermal photons are expected to
follow a double peaks spectral energy distribution (SED) due to the targets
considered. The low energy peak, at around $10^{16}$ eV, can be associated to
the interaction of cosmic-rays with EBL and the neutron beta-decay.
Conversely, the high energy peak, is expected to be the result of cosmic-rays
interaction with the CMB photons [44]. On the other hand the multi-messenger
studies of UHECRs and diffuse gamma-rays cannot constrain the expected
cosmogenic neutrino spectrum [45].

###  2.2 Extragalactic accelerators

As mentioned in the previous subsection, the acceleration of cosmic-rays above
$10^{18}$ eV can occur in astrophysical plasmas when large-scale motion, such
as shocks and turbulent flows, is transferred to individual charged particles.
The maximum energy of accelerated particles, $E_{max}$, can be estimated by
requiring that the gyroradius of the particle be contained in the acceleration
region. Therefore, for a given strength, B, and coherence length, L, of the
magnetic field embedded in an astrophysical region, $E_{max}=Z\cdot e\cdot
B\cdot L$, where $Z\cdot e$ is the charge of the particle. This condition
known as “Hillas limit” [46] can allow $E_{max}\gtrsim 10^{20}$ eV with $Z\sim
1$ for astrophysical environment like the ones present in neutron stars
($B\sim 10^{13}$ G, $L\sim 10$ Km) [47], active galactic nuclei (AGNs) ($B\sim
10^{4}$ G, $L\sim 10$ AU) [48], radio lobes of AGNs ($B\sim 0.1\mu$ G, $L\sim
10$ Kpc), and clusters of galaxies ($B\sim\mu$ G, $L\sim 100$ Kpc) [49]. The
recent detection of a high-energy neutrino event ($\sim 290$ TeV) in
coincidence with a gamma-ray flare, at hundred of GeVs, from the blazar
TXS0506+056 [50, 51] (at a distance of $z\sim 0.3365$ [52]) makes this class
of sources one of the main candidates[53] of the diffuse astrophysical flux
measured by IceCube [54]. Even though this class of sources is extremely
variable [55] the part of neutrinos, above 100 TeV, emitted during a flaring
period can potentially be observed [56]. For blazars (like the sample reported
in Fig. 2) we can expect most of the UHE neutrinos produced through
photohadronic interaction of accelerated protons with photons emitted by
electrons within the jet. For other classes of AGN, like radio galaxies (as
the sample reported in Fig. 3), we can expect an additional UHE neutrino
production inside the lobes, present at the end of the hundred Kpc jets,
through proton-proton inelastic collisions [57]. A recent analysis taking into
account the position of UHE cosmic-rays observed by Pierre Auger [58] and
Telescope Array [2] experiments, a known blazar gamma-ray catalog and the
position of IceCube astrophysical events shows hints of possible correlation
between the two messengers and highlights the fact that these class of
accelerators have the capability of reaching the energy of $\sim 10^{20}$ eV.
On this regards also clusters of galaxies can have the possibility to produce
UHE neutrinos since cosmic rays accelerated up to UHE by AGN activities can be
contained for a certain amount of time inside the cluster due to the presence
of high magnetic field. They can reach values ranging from a few microgauss
($\mu$G) on scales of order $\sim 10$ Kpc for normal clusters [59], up to
$10-40$ $\mu$G on scales of $3-5$ Kpc for cool core clusters [60]. Such strong
magnetic fields can confine cosmic-rays of UHE. While propagating inside the
cluster, cosmic-rays can thus have interactions with the present photonic and
baryonic backgrounds producing gamma-rays and neutrinos.
Another possible candidate to produce UHE neutrinos it is represented by the
newly born millisecond pulsars with a diffuse component reported in fig. 10.
The amount of this neutrino production can be traced though the measurements
of the extra-galactic star formation rate following the model reported by
[61].

Figure 2: Exposure Skymap (in Galactic coordinates) of the presented UHE
neutrino telescope. In the Skymap are reported also the sources of the 2FHL
catalog spatially correlated with a UHECR observed by Auger or Telescope Array
(within $5^{\circ}$). In Black are reported the blazars TXS0506+056 and
PKS1502+106 as the two main neutrino emitter candidates observed up to now. A
maximum angular acceptance of 45∘ from the local zenith is considered. Figure
3: Exposure Skymap (in Galactic coordinates) of the presented UHE neutrino
telescope. In the Skymap also the radio galaxies, with a radio luminosity $\nu
L_{\nu}>2\times 10^{40}$ erg s-1, expected to have the possibilities of
accelerates cosmic-rays up to UHE in the giant lobes [62]. A maximum angular
acceptance of 45∘ from the local zenith is considered.

###  2.3 Extreme Dark Matter candidates

The standard cosmological model is based on the existence of Dark Matter (DM)
which, until now, did not show direct evidences via non-gravitational
interaction. An alternative way to probe DM interactions, beyond the
gravitational one, is the search for indirect signatures of DM decay or
annihilation. These processes can indeed lead to the production of an
astrophysical signal of cosmic-rays, gamma-rays and neutrinos. In the last
decades, the rise of high-energy multi-messenger astronomy has led to huge
improvements in the indirect searches for dark matter. Dark matter candidates
concentrated in the galactic halo and distributed in the intergalactic space
can produce a flux of UHE neutrinos through their decay and annihilation into
ordinary matter. Due to the unitarity bound on DM cross-section [63, 64] we
can expect a higher flux for the decaying case with respect to the
annihilating one. The neutrino flux from decaying dark matter can be expressed
by:

$\phi_{\nu}^{dec}\sim\frac{\rho_{DM}}{m_{DM}}\frac{1}{\tau_{DM}},$ (3)

where $\rho_{DM}\simeq 0.4GeV/cm^{3}$ is the typical galactic density of DM
particles, $\tau_{DM}$ is their lifetime, $m_{DM}$ is the DM mass, and $L\sim
10kpc$ is the length scale of our galaxy. On the other hand, for annihilating
dark matter the corresponding flux is of the order of:

$\phi_{\nu}^{ann}\sim\left(\frac{\rho_{DM}}{m_{DM}}\right)^{2}\sigma v_{DM}L,$
(4)

where $\sigma$ is the annihilation cross-section, $v_{DM}\simeq 10^{-3}c$
corresponds to the typical velocity of dark matter particles and $L$ is the
length scale. The possible values of neutrino fluxes obtained for annihilating
DM result far beyond the sensitivity of next-generation neutrino radio
telescopes, producing a number of neutrino events in these observatories
negligibly small. An exception can be represented by peculiar scenarios
featuring very dense and/or very cold dark matter substructures [65].
Conversely, for decaying DM particles with $\tau_{DM}\simeq 10^{29}$ s and
$m_{DM}=10^{9}$ GeV we can expect a production of neutrino flux within the
reach of upcoming neutrino UHE telescopes [66].

##  3 Experimental setup

The design of acoustic neutrinos detectors has been up to the present
ancillary to the optical ones. Therefore, the larger attenuation length of the
sound in water in respect to the light was never exploited properly. Here we
explore the possibility of detecting a neutrino generated acoustic signal by
instrumenting the area covered by the existing oil rigs structures.
As a best case scenario we consider all the available platforms (see Fig. 4)
and based on the actual geometry we calculated the possible sensitivity of the
acoustic telescope. The study has been performed by accounting for sound
attenuation in the conditions of Adriatic Sea and wave field propagation
characteristics in the velocity structure scenario reported by [67]. A ray
tracing simulation is obtained using a realistic velocity profile for both the
calculations of amplitude attenuation and ray propagation. The dependence of
wavefield propagation respect to the angle of the incident neutrino respect to
the normal vector to the sea surface is studied here for the first time in
shallow waters. In exploration geophysics and in earthquake seismology, the
use of large arrays of sensors is routine since decades. In both fields, the
target is to study the structure of the soil to reveal potential geological
structures that can trap oil or gas. The detection of an acoustic signal with
an extended array of sensors (hydrophones) is a problem that is solved since
decades both in submarine sonar technology and in a similar fashion by array
techniques for seismic studies and earthquake location like f-k analysis and
beam-forming that are explained in e.g. [68]. The problem in our case is made
easier since the source is well modelled by previous studies [69, 70].
We plan thus to instrument 100 platforms, the ones reported with a red
triangle in Fig. 4. Each platform will be equipped with a multi-storey
structure of hydrophones whose distribution in depth will be optimized
according to the simulation of the wave field propagation studies. These will
be based on the experimental data collected at the first R&D stage. In all
platforms, a sub-array named single platform sub array (hereinafter SIPSA)
will be deployed to better solve the f-k analysis uncertainties connected to
spatial aliasing due to possible multiple reflection from the seafloor and
water surface. We therefore plan to install a single projector in any array
element (platform) to calibrate the SIPSA sub-arrays and the entire large
detector. A continuous calibration using known signals resembling the neutrino
expected one according to simulations see Fig. 5(b) will permit to develop
many detection templates that will help to discriminate the signal and avoid
false detections. In total, we will instrument 100 platforms (see Fig. 4) each
of them will be instrumented with the SIPSA unit that can be formed by four to
ten hydrophones. A large advantage of shallow sea depth is indeed the
simplicity of the hardware needed. In fact, a single hydrophone with an
integrated pre-amplifier like [71] has a self-noise level well below sea state
zero and the integrated preamplifier permits the use of a simple analog cable.
A cable no longer than 50 meters permits to use the analog output of the
transducer and digitize it on board of the platform, avoiding the installation
of an expensive underwater digitizer and the development of custom DAQ
hardware. The reduction of the costs allows to multiply the number of sensors,
increase the signal-to-noise ratio as well as the array sensitivity.

###  3.1 The existing infrastructure in Adriatic Sea

The platforms usable for the ANDIAMO experiment are located in the offshore of
Adriatic Sea at a distance from the shore ranging from 10 to 60 kilometres.
Most of them are placed in an elongated shape from north to south, following
roughly the natural shape of the continental margin where oil and gas were
extracted. The depth of the sea is typical of a continental shelf, with depth
ranging from 25 to 80 meters.

Figure 4: In the map the red triangles show the location of the available
platforms and colour coded the bathymetry of the Adriatic Sea.

The shallow depth of Adriatic Sea gives us several advantages, the most
important is in terms of geometrical spreading. The acoustic wave generated by
neutrino interaction never goes in the far field regime where spherical
spreading is characteristic, this is the case for acoustic detectors
implemented on optical arrays for exampele, where the array dimensions are
small compared to the extension of the channel where they are installed. In
our case, the cascade length is comparable to the depth of the channel, there
is no transition region from near to far field propagation. The amplitude of
the acoustic wave will decay as $1/\sqrt{r}$ respect to a spherical spreading
where the decay is $1/r$. In terms of attenuation, we have 10 dB/decade
against 20 dB/decade for spherical spreading. In other words, we deal with a
cylindrical wave instead of a spherical one. The relatively high temperature
of the shallow Adriatic Sea together with a lower salinity due to the fresh
water influx from the largest Italian river, the Po river, will contribute to
decrease the attenuation [72]. The Mediterranean Sea is for many physical and
geological reasons the best environment for an underwater neutrino telescope
as already highlighted by [69]. Another advantage is logistical since mostly
all services needed such power and connectivity are supposed to be already
available on the oil rigs and no complex high pressure waterproof equipment is
needed.

The steady structure of the platform will permit to install the SIPSA mini
array in every platform in a solution similar to the one adopted in [14]. From
another side, we need to carefully study the underwater noise that could be a
limiting factor. The stationary noise floor for our system is mostly dependent
on sea state in the frequency range of detection for the predicted acoustic
signal (1-20 kHz). Transient noises as dolphins and propellers can of course
spoil the sensitivity locally and in an unpredictable manner but in a long
term perspective can be filtered and recognized since they have a different
time /frequency signature compared to the neutrino generate ones. Another
possible issue with the shallow sea can be a complex ray path propagation and
possible reflections and absorption, especially due to the interaction with
sandy or clay seafloor that we will face in a specific paragraph of this
paper. Sound propagation in very shallow waters has been studied largely in
literature by [73, 74, 75, 76] converging to stating that shallow waters
propagation can be seen as a small scale problem respect to deep waters, for
this reason we will apply ray tracing to simulate our environment. Refraction
and reflections are the physical phenomena that guide the direction of the
wave field across the path from the source to the various receivers. The
refraction is connected to the change of the vertical velocity profile that
drives an acoustic impedance change. This profile changes across the year, and
it depends on water temperature, salinity pH, and depth (pressure). For these
reasons, a preliminary study using a realistic and appropriate velocity
profile is necessary. The structure of the velocity profile is indeed a game
changer in the possibility to design an effective acoustic telescope, in
particular we want to exploit the presence of the SOFAR channel to stretch the
acoustic detection range to the maximum. In fact, changing the physical
properties of the medium changes the sensitivity of the detector itself. We
will show that the available geometry of platforms distribution is suitable
for the realization of the largest acoustic array for UHE neutrinos detection
with a sensitive area of more than 10000 Km2

###  3.2 A novel reconstruction technique for acoustic signal

When a UHE neutrino interacts with matter at such high energies, it can decay
according to different channels and generate an hadronic cascade that in turn
gives rise to the thermo-acoustic effect.

(a) Peak pressure versus cascade energy calculated at 100 meters from the
shower.

(b) The waveform of the acoustic pulse at 100 m for a 1e11 GeV shower.

Figure 5: Acoustic parameters of the source.

The extension in depth of the cascade gives the length of the cylinder where
the thermo-acoustic effect occurs. The intensity of this cylindrical wave,
also called in literature pancake, depends on the angle of the cascade with
respect to the vertical and of course on the generating phenomena, i.e. the
possible decay channels. Once the acoustic wave is generated with a given
intensity, the problem is shifted to a proper reconstruction and validation of
the signal using the state-of-art acoustic techniques. The waveform generated
by a UHE neutrino interaction resembles a bipolar signal, whose amplitude at
100 meters from the cascade can range from 10 Pascals to some milliPascal. The
central frequency is expected to be around 50 kHz at 100 meters with a
duration of the order of several microseconds according to simulations
proposed by [70]. This energy deposition is similar to a delta function,
exciting thus a white kind of spectrum in the frequency domain. When this
pulse propagates, it undergoes attenuation due to both geometrical spreading
and absorption. This causes a broadening of the wavelet in time and an obvious
decrease in amplitude. Nowadays, state-of-the-art hydrophones have a self-
noise level that is lower than the sea state level zero of the Wenz sea state
plots [77]. This means that, supposed to have the most sensitive hydrophones,
we are left with the problem of a realistic calculation of the attenuation and
frequency content change in the signal versus distance. This study is
fundamental to assess how far can we stretch the sensitive area of the
detector. The sea noise is a limiting factor to the sensitivity of the array.
This will be discussed in detail in the section about the amplitude spectra.
We will account for geometrical spreading and attenuation in the frequency
domain. The generated signal creates in the far field (> 100 m) a cylindrical
wave that propagates and attenuates with an $1/\sqrt{r}$ attenuation law. The
frequency content of the bipolar signal ranges from 1kHz to 100 kHz and the
peak pressure is expected at a frequency of approximately 50 kHz at 100 meters
distance. The signal undergoes all the acoustic phenomena involved in wave
propagation in a vertically inhomogeneous medium, including refraction and
dispersion/attenuation. The propagating pressure wave field changes the energy
and energy distribution in the frequency domain with respect to distance from
the event that generated it. It is no mystery that attenuation is a frequency
dependent phenomenon and roughly speaking can be summarized with this
statement: the longer wavelengths travels longer from the source. The previous
considerations were made in order to answer to the basic question of this
paper, i.e.: Can an array created by using the existing oil rigs (or at least
a group of them) be suitable for the detection of UHE neutrinos ? The results
showing the feasibility of the ANDIAMO project can be grouped in two parts,
the first regarding the attenuation/dispersion and the second one about the
ray paths simulations addressing the following points:

  * 1.

Establish how far can we ear an acoustic signal generated by a neutrino
induced shower in a relaistic velocity profile scenario.

  * 2.

Determine which are the angular limits for detection in our setup.

####  3.2.1 Amplitude spectra

We calculate the amplitude spectrum of the acoustic pressure reference signal
of Fig. 5(b) i.e. the one intended as the near source signal.

Figure 6: Sound pressure level versus frequency and distance from the source
for different shower energies in a 50 meters deep sea assuming isovelocity
conditions.

In a second step we calculate, for a fixed channel depth and for several
distances and frequencies, the attenuated spectra for several source levels.
The result is reported in Fig. 6. As a comparison, for a fixed energy, we
calculate the dispersion of the spectrum for a 3000 meters deep channel to
highlight the advantage of shallow waters (see Fig. 7). From the spectra, it
is evident that above a cascade energy of $10^{19}$ GeV signal can be measured
by any single array element (oil rig) that stays in a 10-kilometre radius from
the source event.

Figure 7: Calculated spectra for different sea depths with fixed shower
energy, on the left panel a 30 meters deep sea and on the right panel a 3000
meters deep sea, we assume isovelocity conditions.

We assume a linear path propagation from source to receiver induced by an
isovelocity condition for simplicity. We calculate the transmission losses
(TL) for every spectral component of the source signal reported in Fig. 5(b)
and using the following equation:

$TL=10\log\frac{I_{s}}{I(r)}$ (5)

we obtain in the end the sound pressure level spectra at different distances
for different source energies and incidence angle of 0 deg, the angular
dependence of attenuation will be studied in the second part where we will
account for seafloor reflection. We use as source level the sound pressure
level at 100 meters because it can be considered far field compared to the
extension of the neutrino induced cascade (that is not a point-like source)
but still much smaller than the inter distance occurring between the array
elements. An extrapolation of the results shown by [70] permitted to produce a
plot (see Fig. 5(a)) of peak-pressure versus cascade energy at a fixed
distance of 100 m which is more suitable for our study. The results of the
simulated spectra for a 50 meters deep channel for several energies and
distances from the source are reported in Fig. 6.

####  3.2.2 Ray tracing and signal intensity versus cascade angle

The problem of propagation in shallow waters, i.e. inside the continental
shelf, can be complex, and it depends strongly on the velocity structure of
the water column.

Figure 8: Ray tracing using a typical velocity profile for a 50 m deep sea at
different neutrino incidence angle and a shower energy of $10^{19}$ eV.

In our setup, the seafloor is pretty homogeneous, and we have a general mild
increase of the seafloor depth from 20 to 50 meters according to the natural
slope of the continental shelf. Future studies based on the geophysical and
oceanographic acquired data of the detector site will permit to produce more
accurate simulations of the propagation of the signals.

Figure 9: Map illustrating the covered sensitive area for vertical neutrino
induced showers at $10^{19}$ eV , color coded is the SIPSA sub array receiver
density

The typical shallow water conditions are studied largely in underwater
acoustics [73, 74] we can summarize the part of our interest with the
following scenario that is confirmed by experimental data e.g. [67] for
Mediterranean sea. In shallow waters we have a first mixed surface layer where
the speed of sound is constant, followed by a thermocline and an increase due
to pressure as shown in the left top panel of Fig. 8. During winter the entire
water column is completely mixed, and the sun exposure is lower. We have a
less steep thermal gradient and a larger mixing due to the heavy sea
conditions. This creates a condition that can be more similar to isovelocity
conditions. The latter implies the straight ray propagation given the
homogeneity of the medium. By simulating a cylindrical source of 20 meters of
extension in depth (cascade length), we could ray-trace the direction of the
wave-fronts generated at different shower angle for a 5 kHz tone which is,
according to spectral simulations, the peak level at the distances of interest
given the actual platforms geometry. Additional signal attenuation induced by
transmission losses at the sea floor are accounted for rays that are not
confined inside the SOFAR channel. The transmission losses at the sea floor
constitute a severe limiting factor for angles larger than 45 degrees since
after this value we have a very steep drop in the reflection coefficient. We
consider the sea surface as a perfect reflector, hence we account only for sea
bottom losses. We report the results in Fig. 8, we can easily note that at
after 3000 meters for a shower induced at 45∘ the signal drops below the sea
state level zero value for a 5kHz tone. On the contrary for angles between
$\pm 10^{\circ}$ the channeling effect permits the focusing of the energy
inside the the SOFAR channel and extending the range of the detector. This
aspect of wave propagation is limiting the effective sky solid angle that we
can observe with the current setup. Of course a denser array can push further
the angular span that we can cover. In summary the attenuation, geometrical
spreading and reflection/transmission losses are carefully calculated in the
simulation of the acoustic environment of the ANDIAMO experiment. From the
acoustical point of view the sensitive area changes according to angle and
energy of the incident neutrino an example of the covered area is reported in
Fig. 9 where the color map represent the SIPSA sub arrays density. In this
case for a $10^{19}$ eV cascade vertically originated the sensitive area is of
the order of 10000 km2.

##  4 Possible sensitivity to UHE neutrinos

Since the calculation of a detailed exposure for the proposed ANDIAMO array go
beyond the scope of this work, mainly because of a non-defined geometry and
the absence of a effective area obtained though a dedicated Monte Carlo
simulation, a possible approximation is obtained through a semi-analytic
analysis. The limits for the angle and area are set by the results shown in
the previous section where acoustic propagation is studied in detail. The
maximal solid angle around the zenith of the detector is set by the abrupt
fall of the reflection coefficient which in turn limits the total distance
travelled by the acoustic signal generated by an inclined cascade. The
sensitive area is then angle and energy dependent, for a $10^{19}$ eV cascade
it ranges from 10000 Km2 for a vertical cascade to about 2800 Km2 for a 45∘
inclined one. This area variability is accounted in the calculation of the
sensitivity shown in Fig. 10.

For downward-going neutrinos, the calculation of the exposure involves the
ANDIAMO array aperture, the neutrino interaction probability, an ideal
identification efficiency, and the integration in time. A sum over time and
integration in solid angle would yield the exposure ($\mathcal{E}$) to UHE
neutrinos. Assuming a $1:1:1$ flavour ratio (as expected due to the effects of
neutrino oscillations during propagation from the sources to the earth), the
total exposure can be expressed as:

$\mathcal{E}(E_{\nu})=\frac{2\pi}{m}\sum_{i}\left[\sigma_{i}(E_{\nu})\int
dtd\theta dD\sin\theta\cos\theta A^{i}_{eff}(\theta,D,E_{nu},t)\right],$ (6)

where the sum runs over the three neutrino flavours and the CC and NC
interactions, with $\sigma^{i}$ the corresponding $\nu-$nucleon interaction
cross-section [78] and $m$ the nucleon mass. The integral is performed over
the zenith angle $\theta$ the interaction depth $D$ of the neutrino (in units
of g cm-2), and the blind search period. $A^{i}_{eff}$ represents the
effective area. Without having a finalized geometry for the ANDIAMO array, we
did a analytic approach to obtain a possible effective area and a related
sensitivity of this future detector. In particular, we assume $A^{i}_{eff}$ to
be approximated by $Ap_{\nu\rightarrow\mu,e,\tau}(E_{\nu},\langle
E_{\mu,e,\tau}\rangle$) for an array of surface $A$ (from 2800 to 10000 Km2
for the ANDIAMO case). Where $p_{\nu\rightarrow\mu,e,\tau}(E_{\nu})$ can be
described by the following equation:

$p_{\nu\rightarrow\mu,e,\tau}(E_{\nu})=R_{\mu,e,\tau}(\langle
E_{\mu,e,\tau}\rangle)/\lambda_{\nu}(E_{\nu}),$ (7)

with $\lambda_{\nu}(E_{\nu})$ representing the interaction length of a
neutrino with energy $E_{\nu}$ and $R_{\mu,e,\tau}$ describing the leptonic
range in the considered medium. These quantities allow us to roughly obtain a
possible flux limit observable, or in other words, a possible sensitivity of
our detector. More in details, under the assumption that the UHE neutrino flux
behaves like $\Phi(E)=\frac{dN}{dE}=kE^{-2}$ the corresponding sensitivity can
be written as:

$k=\frac{N_{evt}}{\int^{E_{max}}_{E_{min}}E^{-2}\mathcal{E}(E)dE},$ (8)

where the $\mathcal{E}(E)$ correspond to the exposure obtained in equation 6
and $N_{evt}$ indicates the number of expected events. No background events
neither spectral noises were considered in the energy range where ANDIAMO
array can have the capability of reconstructing UHE neutrinos.

Figure 10: In this plot we report the expected sensitivity of the ANDIAMO
concept with a surface array equivalent to the one described in Fig. 4 and a
observational time of one decade. The three neutrino families are considered
with the range of sensitivity spanning the different sensitive areas covered
considering the possible beneficial incident angles. For comparison we show
also the limit on the possible UHE cosmogenic neutrino flux obtained from the
UHECR observations as well as possible UHE neutrino diffuse contributions from
AGNs, millisecond pulsars and heavy dark matter decays. IceCube full-sky
measurements are reported too.

The expected sensitivity, reported in Fig. 10, is compared with the main
expected UHE neutrino diffuse fluxes. Event though more accurate studies can
be obtained with the use of a dedicated Monte Carlo simulation chain whenever
the geometry will be finalized, the preliminary calculations obtained through
an analytical approach highlights the potential discoveries of such acoustic
array concept.

##  5 Conclusions

The advances of neutrino astronomy of the last decade surge the importance of
knowing the neutrino spectral features in the energy range from PeV to EeV.
Several radio arrays prototypes have been recently finalized with the idea of
covering larger surfaces in order to be sensible to the possible
extraterrestrial UHE neutrino flux. An alternative way to detect this flux it
is represented by the underwater/ice acoustic techniques carried on in the
last decades with small prototypes or through the possible use of hydrophones
installed for calibration purposes and not optimized for an acoustic detector
design. In this work we explore the possibility of building a dedicated
acoustic array in the Adriatic Sea using the platforms already installed, and
not operative, for oil and gas upstream.
Detailed calculations on the propagation and attenuation of the signal
generated by UHE neutrinos in shallow waters show the possibility to exploit
the favourable conditions of the Adriatic Sea. This advantage is given by the
physical and geographical conditions of the sea, i.e. mild temperature and low
salinity. Also, the existence of many infrastructures partially instrumented
is a plus. Indeed, it will permit the deployment of the largest acoustic array
for UHE neutrinos detection as well as for marine biology, oceanographic and
geophysical studies. The calculations performed show that it is feasible to
realize an acoustic telescope with an instrumented maximal area of $\sim$10000
Km2 exploiting all the available platforms all over the Adriatic Sea with an
expected sensitivity shown in Fig. 10. A first attempt to estimate the
possible sensitivity shows that, even though the possibility of resolving
single sources of UHE neutrinos is still a challenging purpose, constraining
the neutrino spectral features in the PeV-EeV range can be possible within a
decade of data taking.
Considering the moderate costs of this detection technique and the uniqueness
of this project i.e. complementary to the shower earth-skimming radio array,
we aim to make further studies with a dedicated Monte Carlo simulation
optimizing the geometry of hydrophone lines and analysing more in details the
possible backgrounds with in situ measurements. Moreover, we leave open the
future possibility of taking into account a bigger instrumental area with
additional lines equipped without pre-existing platforms.

## Acknowledgements

The authors would like to thank Giorgio Riccobene, Robert Lahmann, Kay Graf
and Giorgio Gratta for their helpful suggestions and discussions regarding the
detector feasibility study. This work was partially supported by the research
grant number 2017W4HA7S “NAT-NET: Neutrino and Astroparticle Theory Network”
under the program PRIN 2017 funded by the Italian Ministero dell’Universit‘a e
della Ricerca (MUR).

## References

  * [1] J. Abraham, et al., Measurement of the Energy Spectrum of Cosmic Rays above $10^{18}$ eV Using the Pierre Auger Observatory, Phys. Lett. B 685 (2010) 239–246.  arXiv:1002.1975, doi:10.1016/j.physletb.2010.02.013.
  * [2] T. Abu-Zayyad, et al., The Cosmic Ray Energy Spectrum Observed with the Surface Detector of the Telescope Array Experiment, Astrophys. J. Lett. 768 (2013) L1.  arXiv:1205.5067, doi:10.1088/2041-8205/768/1/L1.
  * [3] F. Halzen, D. Hooper, High-energy neutrino astronomy: The Cosmic ray connection, Rept. Prog. Phys. 65 (2002) 1025–1078.  arXiv:astro-ph/0204527, doi:10.1088/0034-4885/65/7/201.
  * [4] J. N. Bahcall, E. Waxman, High-energy astrophysical neutrinos: The Upper bound is robust, Phys. Rev. D 64 (2001) 023002.  arXiv:hep-ph/9902383, doi:10.1103/PhysRevD.64.023002.
  * [5] K. Kotera, D. Allard, A. Olinto, Cosmogenic Neutrinos: parameter space and detectabilty from PeV to ZeV, JCAP 10 (2010) 013.  arXiv:1009.1382, doi:10.1088/1475-7516/2010/10/013.
  * [6] M. Aartsen, et al., Observation of High-Energy Astrophysical Neutrinos in Three Years of IceCube Data, Phys. Rev. Lett. 113 (2014) 101101.  arXiv:1405.5303, doi:10.1103/PhysRevLett.113.101101.
  * [7] M. Aartsen, et al., Evidence for Astrophysical Muon Neutrinos from the Northern Sky with IceCube, Phys. Rev. Lett. 115 (8) (2015) 081102.  arXiv:1507.04005, doi:10.1103/PhysRevLett.115.081102.
  * [8] M. Aartsen, et al., Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector, Science 342 (2013) 1242856.  arXiv:1311.5238, doi:10.1126/science.1242856.
  * [9] K. Murase, S. Nagataki, High energy neutrino emission and neutrino background from gamma-ray bursts in the internal shock model, Phys. Rev. D 73 (2006) 063002\.  arXiv:astro-ph/0512275, doi:10.1103/PhysRevD.73.063002.
  * [10] D. Biehl, D. Boncioli, A. Fedynitch, W. Winter, Cosmic-Ray and Neutrino Emission from Gamma-Ray Bursts with a Nuclear Cascade, Astron. Astrophys. 611 (2018) A101.  arXiv:1705.08909, doi:10.1051/0004-6361/201731337.
  * [11] V. Berezinsky, G. Zatsepin, Cosmic rays at ultrahigh-energies (neutrino?), Phys. Lett. B 28 (1969) 423–424.  doi:10.1016/0370-2693(69)90341-4.
  * [12] M. G. Aartsen, et al., The IceCube Neutrino Observatory: Instrumentation and Online Systems, JINST 12 (03) (2017) P03012.  arXiv:1612.05093, doi:10.1088/1748-0221/12/03/P03012.
  * [13] M. Ageron, et al., ANTARES: the first undersea neutrino telescope, Nucl. Instrum. Meth. A 656 (2011) 11–38.  arXiv:1104.1607, doi:10.1016/j.nima.2011.06.103.
  * [14] V. Aynutdinov, et al., Acoustic search for high-energy neutrinos in the Lake Baikal: Results and plans, Nucl. Instrum. Meth. A 662 (2012) S210–S215.  doi:10.1016/j.nima.2010.11.153.
  * [15] S. Adrian-Martinez, et al., Letter of intent for KM3NeT 2.0, J. Phys. G 43 (8) (2016) 084001.  arXiv:1601.07459, doi:10.1088/0954-3899/43/8/084001.
  * [16] G. Askar’yan, Excess negative charge of an electron-photon shower and its coherent radio emission, Sov. Phys. JETP 14 (2) (1962) 441–443.
  * [17] I. Kravchenko, et al., Performance and simulation of the RICE detector, Astropart. Phys. 19 (2003) 15–36.  arXiv:astro-ph/0112372, doi:10.1016/S0927-6505(02)00194-9.
  * [18] I. Kravchenko, S. Hussain, D. Seckel, D. Besson, E. Fensholt, J. Ralston, J. Taylor, K. Ratzlaff, R. Young, Updated Results from the RICE Experiment and Future Prospects for Ultra-High Energy Neutrino Detection at the South Pole, Phys. Rev. D 85 (2012) 062004.  arXiv:1106.1164, doi:10.1103/PhysRevD.85.062004.
  * [19] P. Allison, et al., Performance of two Askaryan Radio Array stations and first results in the search for ultrahigh energy neutrinos, Phys. Rev. D 93 (8) (2016) 082003.  arXiv:1507.08991, doi:10.1103/PhysRevD.93.082003.
  * [20] S. Barwick, et al., Design and Performance of the ARIANNA HRA-3 Neutrino Detector Systems, IEEE Trans. Nucl. Sci. 62 (5) (2015) 2202–2215.  arXiv:1410.7369, doi:10.1109/TNS.2015.2468182.
  * [21] S. Barwick, et al., A First Search for Cosmogenic Neutrinos with the ARIANNA Hexagonal Radio Array, Astropart. Phys. 70 (2015) 12–26.  arXiv:1410.7352, doi:10.1016/j.astropartphys.2015.04.002.
  * [22] J. Aguilar, et al., Design and Sensitivity of the Radio Neutrino Observatory in Greenland (RNO-G)arXiv:2010.12279.
  * [23] J. Álvarez-Muñiz, et al., The Giant Radio Array for Neutrino Detection (GRAND): Science and design, Science China Physics, Mechanics, and Astronomy 63 (1) (2020) 219501.  arXiv:1810.09994, doi:10.1007/s11433-018-9385-7.
  * [24] P. Gorham, et al., The Antarctic Impulsive Transient Antenna Ultra-high Energy Neutrino Detector Design, Performance, and Sensitivity for 2006-2007 Balloon Flight, Astropart. Phys. 32 (2009) 10–41.  arXiv:0812.1920, doi:10.1016/j.astropartphys.2009.05.003.
  * [25] P. Allison, et al., Constraints on the diffuse flux of ultrahigh energy neutrinos from four years of Askaryan Radio Array data in two stations, Phys. Rev. D 102 (4) (2020) 043021.  arXiv:1912.00987, doi:10.1103/PhysRevD.102.043021.
  * [26] N. Kurahashi, J. Vandenbroucke, G. Gratta, Search for Acoustic Signals from Ultra-High Energy Neutrinos in 1500 km3 of Sea Water, Phys. Rev. D 82 (2010) 073006.  arXiv:1007.5517, doi:10.1103/PhysRevD.82.073006.
  * [27] N. G. Lehtinen, S. Adam, G. Gratta, T. K. Berger, M. J. Buckingham, Sensitivity of an underwater acoustic array to ultra-high energy neutrinos, Astroparticle Physics 17 (3) (2002) 279–292.  doi:10.1016/s0927-6505(01)00158-x.
URL https://doi.org/10.1016/s0927-6505(01)00158-x

  * [28] G. Askarian, B. Dolgoshein, A. Kalinovsky, N. Mokhov, ACOUSTIC DETECTION OF HIGH-ENERGY PARTICLE SHOWERS IN WATER, Nucl. Instrum. Meth. 164 (1979) 267–278.  doi:10.1016/0029-554X(79)90244-1.
  * [29] J. G. Learned, Acoustic Radiation by Charged Atomic Particles in Liquids: An Analysis, Phys. Rev. D 19 (1979) 3293.  doi:10.1103/PhysRevD.19.3293.
  * [30] Y. Abdou, et al., Design and performance of the South Pole Acoustic Test Setup, Nucl. Instrum. Meth. A 683 (2012) 78–90.  arXiv:1105.4339, doi:10.1016/j.nima.2012.03.027.
  * [31] G. Riccobene, et al., Long-term measurements of acoustic background noise in very deep sea, Nucl. Instrum. Meth. A 604 (2009) S149–S157.  doi:10.1016/j.nima.2009.03.195.
  * [32] S. Danaher, First data from ACoRNE and signal processing techniques, J. Phys. Conf. Ser. 81 (2007) 012011.  doi:10.1088/1742-6596/81/1/012011.
  * [33] J. Aguilar, et al., AMADEUS - The Acoustic Neutrino Detection Test System of the ANTARES Deep-Sea Neutrino Telescope, Nucl. Instrum. Meth. A 626-627 (2011) 128–143.  arXiv:1009.4179, doi:10.1016/j.nima.2010.09.053.
  * [34] J. Vandenbroucke, G. Gratta, N. Lehtinen, Experimental study of acoustic ultrahigh - energy neutrino detection, Astrophys. J. 621 (2005) 301–312.  arXiv:astro-ph/0406105, doi:10.1086/425336.
  * [35] F. Simeone, A. Capone, Acoustic detection of UHE neutrinos in the Mediterranean sea: status and perspective, in: European Physical Journal Web of Conferences, Vol. 135 of European Physical Journal Web of Conferences, 2017, p. 06002.  doi:10.1051/epjconf/201713506002.
  * [36] M. G. Aartsen, et al., IceCube-Gen2: the window to the extreme Universe, J. Phys. G 48 (6) (2021) 060501.  arXiv:2008.04323, doi:10.1088/1361-6471/abbd48.
  * [37] V. Berezinsky, A. Z. Gazizov, S. I. Grigorieva, On astrophysical solution to ultrahigh-energy cosmic rays, Phys. Rev. D 74 (2006) 043005.  arXiv:hep-ph/0204357, doi:10.1103/PhysRevD.74.043005.
  * [38] R. Aloisio, V. Berezinsky, A. Gazizov, Transition from galactic to extragalactic cosmic rays, Astropart. Phys. 39-40 (2012) 129–143.  arXiv:1211.0494, doi:10.1016/j.astropartphys.2012.09.007.
  * [39] P. de Bernardis, et al., A Flat universe from high resolution maps of the cosmic microwave background radiation, Nature 404 (2000) 955–959.  arXiv:astro-ph/0004404, doi:10.1038/35010035.
  * [40] A. Franceschini, G. Rodighiero, M. Vaccari, The extragalactic optical-infrared background radiations, their time evolution and the cosmic photon-photon opacity, Astron. Astrophys. 487 (2008) 837.  arXiv:0805.1841, doi:10.1051/0004-6361:200809691.
  * [41] A. Mucke, J. P. Rachen, R. Engel, R. J. Protheroe, T. Stanev, On photohadronic processes in astrophysical environments, Publ. Astron. Soc. Austral. 16 (1999) 160.  arXiv:astro-ph/9808279, doi:10.1071/AS99160.
  * [42] K. Greisen, End to the cosmic ray spectrum?, Phys. Rev. Lett. 16 (1966) 748–750.  doi:10.1103/PhysRevLett.16.748.
  * [43] G. T. Zatsepin, V. A. Kuzmin, Upper limit of the spectrum of cosmic rays, JETP Lett. 4 (1966) 78–80.
  * [44] S. Das, S. Razzaque, N. Gupta, Modeling the spectrum and composition of ultrahigh-energy cosmic rays with two populations of extragalactic sources, Eur. Phys. J. C 81 (1) (2021) 59.  arXiv:2004.07621, doi:10.1140/epjc/s10052-021-08885-4.
  * [45] M. Ackermann, et al., The spectrum of isotropic diffuse gamma-ray emission between 100 MeV and 820 GeV, Astrophys. J. 799 (2015) 86.  arXiv:1410.3696, doi:10.1088/0004-637X/799/1/86.
  * [46] A. M. Hillas, The Origin of Ultrahigh-Energy Cosmic Rays, Ann. Rev. Astron. Astrophys. 22 (1984) 425–444.  doi:10.1146/annurev.aa.22.090184.002233.
  * [47] E. Waxman, J. N. Bahcall, High-energy neutrinos from cosmological gamma-ray burst fireballs, Phys. Rev. Lett. 78 (1997) 2292–2295.  arXiv:astro-ph/9701231, doi:10.1103/PhysRevLett.78.2292.
  * [48] C. M. Urry, P. Padovani, Unified schemes for radio-loud active galactic nuclei, Publ. Astron. Soc. Pac. 107 (1995) 803.  arXiv:astro-ph/9506063, doi:10.1086/133630.
  * [49] V. S. Berezinsky, P. Blasi, V. S. Ptuskin, Clusters of Galaxies as a Storage Room for Cosmic Rays, Astrophys. J. 487 (1997) 529–535.  arXiv:astro-ph/9609048, doi:10.1086/304622.
  * [50] M. G. Aartsen, et al., Multimessenger observations of a flaring blazar coincident with high-energy neutrino IceCube-170922A, Science 361 (6398) (2018) eaat1378.  arXiv:1807.08816, doi:10.1126/science.aat1378.
  * [51] M. G. Aartsen, et al., Neutrino emission from the direction of the blazar TXS 0506+056 prior to the IceCube-170922A alert, Science 361 (6398) (2018) 147–151.  arXiv:1807.08794, doi:10.1126/science.aat2890.
  * [52] S. Paiano, R. Falomo, A. Treves, R. Scarpa, The redshift of the BL Lac object TXS 0506+056, Astrophys. J. Lett. 854 (2) (2018) L32.  arXiv:1802.01939, doi:10.3847/2041-8213/aaad5e.
  * [53] K. Murase, Y. Inoue, C. D. Dermer, Diffuse Neutrino Intensity from the Inner Jets of Active Galactic Nuclei: Impacts of External Photon Fields and the Blazar Sequence, Phys. Rev. D 90 (2) (2014) 023007.  arXiv:1403.4089, doi:10.1103/PhysRevD.90.023007.
  * [54] M. G. Aartsen, et al., All-sky Search for Time-integrated Neutrino Emission from Astrophysical Sources with 7 yr of IceCube Data, Astrophys. J. 835 (2) (2017) 151.  arXiv:1609.04981, doi:10.3847/1538-4357/835/2/151.
  * [55] A. Marinelli, J. R. Sacahui, A. Sharma, M. Osorio-Archila, Analyzing the gamma-ray activity of neutrino emitter candidates: comparing TXS 0506+056 with other blazars, Montly Notices of the Royal Astronomical Societydoi:10.1093/mnras/stab1312.
  * [56] F. Oikonomou, K. Murase, P. Padovani, E. Resconi, P. Mészáros, High energy neutrino flux from individual blazar flares, Mon. Not. Roy. Astron. Soc. 489 (3) (2019) 4347–4366.  arXiv:1906.05302, doi:10.1093/mnras/stz2246.
  * [57] N. Fraija, E. Aguilar-Ruiz, A. Galván-Gámez, A. Marinelli, J. A. de Diego, Study of the PeV neutrino, $\gamma$-rays, and UHECRs around the lobes of Centaurus A, Mon. Not. Roy. Astron. Soc. 481 (4) (2018) 4461–4471.  arXiv:1709.05766, doi:10.1093/mnras/sty2561.
  * [58] A. Aab, et al., The Pierre Auger Cosmic Ray Observatory, Nucl. Instrum. Meth. A 798 (2015) 172–213.  arXiv:1502.01323, doi:10.1016/j.nima.2015.06.058.
  * [59] T. E. Clarke, Faraday rotation observations of magnetic fields in galaxy clusters, J. Korean Astron. Soc. 37 (5) (2004) 337–342.  arXiv:astro-ph/0412268, doi:10.5303/JKAS.2004.37.5.337.
  * [60] T. A. Ensslin, C. Vogt, Magnetic turbulence in cool cores of galaxy clusters, Astron. Astrophys. 453 (2006) 447.  arXiv:astro-ph/0505517, doi:10.1051/0004-6361:20053518.
  * [61] K. Fang, K. Kotera, K. Murase, A. V. Olinto, Testing the Newborn Pulsar Origin of Ultrahigh Energy Cosmic Rays with EeV Neutrinos, Phys. Rev. D 90 (10) (2014) 103005, [Erratum: Phys.Rev.D 92, 129901 (2015)].  arXiv:1311.2044, doi:10.1103/PhysRevD.90.103005.
  * [62] J. H. Matthews, A. R. Bell, K. M. Blundell, A. T. Araudo, Ultrahigh energy cosmic rays from shocks in the lobes of powerful radio galaxies, Mon. Not. Roy. Astron. Soc. 482 (4) (2019) 4303–4321.  arXiv:1810.12350, doi:10.1093/mnras/sty2936.
  * [63] K. Griest, M. Kamionkowski, Unitarity Limits on the Mass and Radius of Dark Matter Particles, Phys. Rev. Lett. 64 (1990) 615.  doi:10.1103/PhysRevLett.64.615.
  * [64] J. Smirnov, J. F. Beacom, TeV-Scale Thermal WIMPs: Unitarity and its Consequences, Phys. Rev. D 100 (4) (2019) 043029.  arXiv:1904.11503, doi:10.1103/PhysRevD.100.043029.
  * [65] J. Zavala, Galactic PeV neutrinos from dark matter annihilation, Phys. Rev. D 89 (12) (2014) 123516.  arXiv:1404.2932, doi:10.1103/PhysRevD.89.123516.
  * [66] M. Chianese, D. F. G. Fiorillo, R. Hajjar, G. Miele, S. Morisi, N. Saviano, Heavy decaying dark matter at future neutrino radio telescopes, JCAP 05 (2021) 074.  arXiv:2103.03254, doi:10.1088/1475-7516/2021/05/074.
  * [67] S. Salon, A. Crise, P. Picco, E. de Marinis, O. Gasparini, Sound speed in the mediterranean sea: an analysis from a climatological data set, Annales Geophysicae 21 (3) (2003) 833–846.  doi:10.5194/angeo-21-833-2003.
URL https://angeo.copernicus.org/articles/21/833/2003/

  * [68] J. Schweitzer, J. Fyen, S. Mykkeltveit, S. J. Gibbons, M. Pirli, D. Kühn, T. Kværna, Seismic arrays, New Manual of Seismological Observatory Practice 2 (NMSOP2)doi:10.2312/GFZ.NMSOP-2_CH9.
URL https://gfzpublic.gfz-potsdam.de/pubman/item/item_43213

  * [69] V. Niess, V. Bertin, Underwater acoustic detection of ultra high energy neutrinos, Astroparticle Physics 26 (4-5) (2006) 243–256.  doi:10.1016/j.astropartphys.2006.06.005.
URL https://doi.org/10.1016/j.astropartphys.2006.06.005

  * [70] S. Bevan, A. Brown, S. Danaher, J. Perkin, C. Rhodes, T. Sloan, L. Thompson, O. Veledar, D. Waters, Study of the acoustic signature of uhe neutrino interactions in water and ice, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 607 (2) (2009) 398–411.  doi:https://doi.org/10.1016/j.nima.2009.05.009.
URL https://www.sciencedirect.com/science/article/pii/S0168900209009401

  * [71] B. . K. S. . V. Measurement, Type 8106 low-noise hydrophone with preamplifier, in: Modeling and Measurement Methods for Acoustic Waves and for Acoustic Microdevices, InTech, 2021.
URL https://www.bksv.com/en/transducers/acoustic/microphones/hydrophones/8106

  * [72] L. E. Kinsler, A. R. Frey, A. B. Coppens, J. V. Sanders, Fundamentals of Acoustics, 4th Edition, 1999.
  * [73] W. A. Kuperman, J. F. Lynch, Shallow-water acoustics, Physics Today 57 (10) (2004) 55–61.  doi:10.1063/1.1825269.
URL https://doi.org/10.1063/1.1825269

  * [74] E. de Sousa Costa, E. Bauzer, J. B. C. Filardi, Underwater acoustics modeling in finite depth shallow waters, in: Modeling and Measurement Methods for Acoustic Waves and for Acoustic Microdevices, InTech, 2013.  doi:10.5772/56132.
URL https://doi.org/10.5772/56132

  * [75] F. Chen, Y. Li, L. Lv, Analysis of shallow water sound velocity profile impact on detection performance of active sonar signal, in: 2016 IEEE/OES China Ocean Acoustics (COA), IEEE, 2016.  doi:10.1109/coa.2016.7535676.
URL https://doi.org/10.1109/coa.2016.7535676

  * [76] M. B. Porter, M. Siderius, Acoustic propagation in very shallow water, in: 2010 International WaterSide Security Conference, IEEE, 2010.  doi:10.1109/wssc.2010.5730274.
URL https://doi.org/10.1109/wssc.2010.5730274

  * [77] G. M. Wenz, Acoustic ambient noise in the ocean: Spectra and sources, The Journal of the Acoustical Society of America 34 (12) (1962) 1936–1956.  arXiv:https://doi.org/10.1121/1.1909155, doi:10.1121/1.1909155.
URL https://doi.org/10.1121/1.1909155

  * [78] A. Cooper-Sarkar, S. Sarkar, Predictions for high energy neutrino cross-sections from the ZEUS global PDF fits, JHEP 01 (2008) 075.  arXiv:0710.5303, doi:10.1088/1126-6708/2008/01/075.


# Square-well model for superconducting pair-potential

Erkki Thuneberg  QTF Centre of Excellence, Aalto University, Espoo, FI-00076,
Finland

###### Abstract

We study Andreev reflection in a one-dimensional square-well pair-potential.
We discuss the history of the model. The current-phase relation is presented
as a sum over Matsubara frequencies. How the current arises from bound and
continuum levels is found by analytic continuation. We discuss two limiting
cases of the square-well potential, the zero-length well and the infinite
well. The model is quantitatively valid in some cases but forms the basis for
understanding a wide range of problems in inhomogeneous superconductivity and
superfluidity.

###### keywords:

Andreev bound state, SNS junction, square-well pair-potential, Josephson
current

Andreev reflection is a fundamental process in superconductivity. The process
consists of a particle type excitation being reflected as hole type
excitation. In contrast to usual reflection, the momentum change in Andreev
reflection is small compared to the particle’s momentum, as the hole type
excitation moves in the opposite direction than its momentum. Andreev
reflection treats the transport of electricity and heat differently. It is the
essential process that allows electric current to flow in any inhomogeneous
superconductor but it resists the heat current.

The purpose of this article is to study Andreev reflection in a simple model.
We model the pair potential by a square-well form, see Fig. 1. This could be
compared with the model of square-well potential  and its special case the
infinite potential well, that are studied in practically all textbooks of
quantum mechanics. The study of such models has a twofold justification. 1)
Although the model is not precisely realized in Nature, it still forms a
reasonable approximation for several applications. 2) A model that can be
solved precisely allows one to learn physics, which is useful also in cases
where the model ceases to be correct quantitatively.

Figure 1: The square-well potential for the pair potential $\Delta(x)$: it
equals $\Delta_{0}e^{-{\rm i}\phi/2}$ for $x<-L/2$, $0$ for $|x|<L/2$, and
$\Delta_{0}e^{{\rm i}\phi/2}$ for $x>L/2$. The figure also shows bound state
energies (dashed lines). These are evaluated with parameter values
$L\Delta/\hbar v_{F}=4.5$ and $\phi=1$ rad. At energies $E>\Delta_{0}$ the
states form a continuum.

We specify the parameters of the square-well model (Fig. 1). We consider the
pair potential $\Delta(x)$ that depends on a single coordinate $x$. The
potential is assumed to vanish in the well of width $L$. Outside of this
range, the amplitude is constant, $|\Delta|=\Delta_{0}$. We allow a phase
difference $\phi$ between the two sides. Because the absolute phase is not
important, it is convenient to choose the phase symmetrically,
$\Delta(x)=\Delta_{0}e^{-{\rm i}\phi/2}$ for $x<-L/2$, and
$\Delta=\Delta_{0}e^{{\rm i}\phi/2}$ for $x>L/2$. It is convenient to call the
middle region as normal (N), and the surrounding regions as superconducting
banks (S). We can think that $x$ is a coordinate along a quasiparticle
trajectory in a bulk metal. Alternatively, we can consider the system to
represent a single quantum channel. (For an introduction to quantum channels
see Ref. [1].) We neglect all scattering (by interfaces, walls, impurities or
lattice vibrations) except the one caused by the pair potential. We count the
energies starting from the Fermi level. For brevity, the bank gap amplitude
$\Delta_{0}$ is denoted by $\Delta$ in the following.

The motivation for Andreev was to understand experimental observations that
the heat conductivity in the intermediate state was found much smaller than in
the normal or zero-field superconducting state. The intermediate state of type
I superconductors, which occurs in the magnetic field, was known to have a
lamellar structure consisting of alternating normal and superconducting
regions. This naturally leads to thinking that the normal-superconducting
interfaces could be responsible for the reduced heat conductivity. The
problem, however, is that the force on the electron at such an interface is
too weak to cause ordinary scattering of electrons that have momentum near the
Fermi momentum. Andreev solved this puzzle by inventing a new type of
reflection, where an incident electron is reflected as a hole-type excitation.
The latter travels in the opposite direction at a momentum that is only
slightly changed from that of the incident electron. In the first publication
on the topic [2], he used Gor’kov’s equations to solve the reflection problem
for a pair-potential step. More generally, he showed that the result for a
realistic N-S interface is not essentially different for electron energies
slightly above the gap, the relevant ones for heat conduction at a low
temperature. This allowed him to estimate the heat conductivity of the
electrons at low temperature, in good agreement with experiment.

Andreev studied the problem of the square-well potential (Fig. 1) in Ref. [3].
At energies below $\Delta$, an electron is reflected from one step as a hole,
which is then reflected back from the other step as an electron. This leads to
the formation of discrete energy levels, which we now call Andreev bound
states. For low-energy levels, he found equally spaced energy values

$\displaystyle E_{n}=\frac{\pi\hbar v_{F}}{L}(n+\gamma).$ (1)

Here $v_{F}$ is the Fermi velocity, $n$ is an integer and $\gamma$ is a
constant, and it is assumed $L\gg\hbar v_{F}/\Delta$. Note that the level
spacing is the same as one would obtain in a deep ordinary potential well at
energies close to the Fermi energy. Andreev applied this to the intermediate
state, where the discreteness of the levels should be visible in thermodynamic
properties at temperatures $T\lesssim\hbar v_{F}/L$.

The next important development was by Kulik [4]. He had the idea that the
square-well pair-potential could be used as a model of the SNS Josephson
junction. For that, he took into account that there could be a phase
difference between the superconducting banks. For the bound state energies, he
found the equation

$\displaystyle E_{n,\pm}=\frac{\hbar v_{F}}{2L}\left(\pm\phi+2\pi
n+2\arccos\frac{E_{n,\pm}}{\Delta}\right).$ (2)

This formula is valid for all bound state energies and arbitrary $L$. For low-
energy levles, the last term in the parenthesis approaches $\pi$. Thus the
energy (2) reduces to the Andreev’s result (1) but now gives a phase-dependent
$\gamma(\phi)=\frac{1}{2}\pm\phi/2\pi$. The levels that increase in energy
with increasing $\phi$ correspond to momentum to the right. This means that
there is an electron moving to the right and a hole to the left. For levels
that decrease in energy with increasing $\phi$, it is the opposite. Thus,
depending on the phase difference and the occupations of these levels, there
is a current through the system. Thus the SNS system forms a Josephson
junction.

Kulik went on to calculate the current. His method was complicated. In
particular, the contributions from the bound states and continuum were
calculated separately. The result reported in Ref. [4] turned out to be
incorrect. In particular, the continuum contribution was found to have a phase
dependence of $\sin\phi$, contrary to later calculations to be discussed
below.

The Josephson current was next calculated by Ishii [5]. His method was also
rather complicated. An essential advantage in his calculation was that he used
the imaginary energies $\epsilon_{m}=\pi T(2m+1)$ (also known as Matsubara
frequencies or energies), where $m$ is an integer. This allowed him to write a
single compact expression that included both bound state and continuum
contributions. His general expression [Eq. (3.5) in [5]] can be written in the
form

$\displaystyle
J=\frac{2e}{\hbar}T\sum_{m=-\infty}^{\infty}\frac{\Delta^{2}\sin\phi}{(2\epsilon_{m}^{2}+\Delta^{2})\cosh\frac{2\epsilon_{m}L}{\hbar
v_{F}}+2\epsilon_{m}\sqrt{\epsilon_{m}^{2}+\Delta^{2}}\sinh\frac{2\epsilon_{m}L}{\hbar
v_{F}}+\Delta^{2}\cos\phi}.$ (3)

The same result has been obtained in several publications using different
approaches [6, 7, 8, 9, 10, 11, 12, 13]. We comment on some of the
calculations. Svidzinskii, Antsygina, and Bratus’ solved the problem both by
using the Gorkov equations and by using the Bololiubov-de Gennes equations,
and discussed the problem in the earlier work by Kulik [6, 7]. Bezuglyi, Kulik
and Mitsai already went on to study more general problems, but commented on
the earlier calculations and provided alternative expressions for the current
in limiting cases [8]. Kulik and Omel’yanchuk used the quasiclassical Green’s
function satisfying the Eilenberger equation to study the special case $L=0$
[14, 15]. Kupriyanov used the same method in the case of general $L$ and
derived the current in the same form as given in (3) [9]. He also discussed
several limiting cases. A calculation aimed for conducting channels and based
on Bololiubov-de Gennes equations is given by Furusaki, Takayanagi and Tsukada
[10]. A same type of derivation with plenty of explanations was given by
Bagwell [11]. The derivation using the Eilenberger equation is presented in
the review [12] together with numerous studies of the Josephson current in
various models. Recently the derivation of (3) has been presented once more in
response to false claims about its invalidity [13]. In summary, formula (3)
was discovered over 50 years ago, and it has formed the starting point for
numerous other studies since.

Equation (3) is a highly convenient expression to calculate the current. The
problem in (3) is its rather abstract form as a sum over imaginary energies.
In order to gain physical understanding, one should do an analytic
continuation of the imaginary energies $E=i\epsilon_{m}$ to the real $E$ axis.
We notice that the denominator of (3) vanishes on the real axis at points
$E=E_{n}$ satisfying

$\displaystyle\phi(E)=\arccos\left[-(1-2\frac{E^{2}}{\Delta^{2}})\cos\frac{2EL}{\hbar
v_{F}}+2\frac{E}{\Delta}\sqrt{1-\frac{E^{2}}{\Delta^{2}}}\sin\frac{2EL}{\hbar
v_{F}}\right].$ (4)

This is just another expression for the bound state energies given in (2). We
express the sum over $\epsilon_{m}$ as an integral and deform the integration
path as explained in many textbooks, for example in Ref. [16], and also in
Ref. [5]. The result is

$\displaystyle J$ $\displaystyle=$
$\displaystyle\frac{2e}{\hbar}\sum_{n}\frac{dE_{n}}{d\phi}f(E_{n})+\frac{2e}{\hbar}\left(\int_{-\infty}^{-\Delta}+\int_{\Delta}^{\infty}\right)$
(5)
$\displaystyle\times\frac{\frac{2}{\pi}|E|\sqrt{E^{2}-\Delta^{2}}\sin\frac{2EL}{\hbar
v_{F}}\Delta^{2}\sin\phi}{[(\Delta^{2}-2E^{2})\cos\frac{2EL}{\hbar
v_{F}}+\Delta^{2}\cos\phi]^{2}+4E^{2}(E^{2}-\Delta^{2})\sin^{2}\frac{2EL}{\hbar
v_{F}}}f(E)dE.$

The first term is the contribution from the bound states at $|E|<\Delta$. It
is the sum over the energy values satisfying (4). The second term is the
contribution from the branch cut of the square root in the continuum
$|E|>\Delta$. In both terms a multiplier is the Fermi function
$f(E)=1/(e^{E/T}+1)$. The derivation given here proves formula (5) for
equilibrium occupations of the levels, but the equation is valid also for
distributions $f(E)$ differing from $1/(e^{E/T}+1)$.

A noteworthy feature of the expression (5) is that it is a sum and integral
over both positive and negative energies. This is a consequence of the Green’s
function technique, where a hole-type excitation (a particle is removed from
the ground state) appears as a negative-energy pole and a particle-type
excitation (a particle added to the ground state) as a positive-energy pole.
We call the presentation of both negative and positive energies the
semiconductor picture, since it makes the superconducting gap look like the
gap between the valence and conduction bands in a semiconductor. An
alternative to the semiconductor picture is the excitation picture, which
shows only the excitation energies, which are positive for both types of
excitations. This is used in Fig. 1. In the following, we prefer the
excitation picture for two reasons. Firstly, consider a level of a given wave
vector that decreases in energy from positive to negative. In the
semiconductor picture this goes smoothly, in contrast to the excitation
picture, where an excitation of the given momentum disappears and an
excitation of the opposite momentum is created. Secondly, the equilibrium
current (the supercurrent) can be interpreted to arise mainly from filled
negative energy levels since $f(E)$ in (5) is small for the positive energy
levels (and vanishes at $T=0$). The semiconductor picture is a direct
generalization of the simple model of the normal state, where the free
electron levels are filled up to the Fermi energy.

We study expression (5) for the current by evaluating the different terms
numerically at zero temperature, $T=0$. The results are presented in Figures
2-4. The difference between the three figures is the length of the normal
region. The left panels in the figures show the energies of the bound states
as functions of the phase difference $\phi$. These are conveniently obtained
by plotting $\phi$ as a function of $E$ (4) and then exchanging the axes.
Because of the symmetry of the energy levels with respect to $E=0$, we discuss
only the levels at negative energies, which are filled in equilibrium at
$T=0$. With increasing $\phi$ the levels degenerate at $\phi=0$ shift up or
down in energy depending on their direction of momentum, as illustrated in
Ref. [13]. In the case of Fig. 2, there are two bound states for small $\phi$.
At a critical $\phi_{c}=2$ the negative-momentum level merges with the
continuum, so that only the positive-momentum bound state survives in the
range $\phi_{c}<\phi<\pi$.

Figure 2: Numerical results in the square-well pair-potential model with a
short N region. (left) The bound state energies as functions of the phase
difference $\phi$ from (4). (middle) The contributions of the filled bound
states ($E<0$) to the current are given by the derivatives of the bound state
energies (5). (right) The total current (blue) (3) and the contribution of the
continuum (orange) [integral term in (5)]. The parameters are $L\Delta/\hbar
v_{F}=1$ and $T=0$.

The contribution of a bound state to the current (5) is given by the
derivative of its energy with respect to the phase. These are evaluated in the
middle panels of Figs. 2-4. In the case of Fig. 2, we see that the currents of
the two filled bound states are in opposite directions. Their currents cancel
at $\phi=0$. With increasing $\phi$ the negative contribution decreases
leaving a net positive current (in units of $2e\Delta/\hbar$, where $e$ is the
electron charge). At $\phi>\phi_{c}$ the negative-momentum level has moved to
the continuum so that the bound-state current comes from the single positive-
momentum level.

The right panels of Figs. 2-4 show the continuum current (5) by orange lines.
We see that it vanishes only at special points $\phi=n\pi$ with integer $n$.
In the case of Fig. 2, it has a maximal negative value at the cusp at
$\phi=\phi_{c}$, where the negative-momentum bound state disappears. The total
current is given by the blue line. This equals the sum of the bound-state
currents (middle panel) and the continuum contribution (orange line). We see
that the cusp in the continuum current just cancels the vanishing of bound-
state current of the negative-momentum level at $\phi_{c}$. The total current
does not exhibit any special behavior at $\phi_{c}$. Similar conclusions are
obtained in all Figs. 2-4. The insensitivity of the current on the bound
states near the gap edge was related to scattering theory in Refs. [6, 7]. It
can also be argued that the current is determined by levels near the Fermi
level and thus should be independent of the details at the relatively high
energy of the gap edge [13]. Further plots of the different contributions to
the current can be found in Refs. [11, 17].

Figure 3: Same as Fig. 2 but for a longer N region, $L\Delta/\hbar v_{F}=10$.
There are more bound states. Their currents alternate between values $\pm
0.045$ approximately, except the energy level that merges with the continuum
at $\phi_{c}=1.150$.

Figure 4: Same as Figs. 2 and 3 but $L\Delta/\hbar v_{F}=5.7$. This is an
example of a bound state emerging from the continuum with increasing $\phi$.

We see from the plots above that the zero-temperature current is discontinuous
at $\phi=\pi$. The current in the neighborhood of this point is

$\displaystyle J(\phi=\pi\mp 0^{+},T=0)=\pm\frac{ev_{F}}{L^{*}},\quad
L^{*}=L+\frac{\hbar v_{F}}{\Delta}.$ (6)

This current comes solely from the two levels that cross the Fermi level at
$\phi=\pi$, as the contributions from all other bound states cancel each other
and the continuum contribution vanishes. The quantity $L^{*}$ can be
interpreted as the effective length of a bound state, where the latter term in
$L^{*}$ describes the penetration of the bound state in the superconducting
banks [18].

The square-well pair-potential model has a few limiting cases. The case of the
zero-length potential well, $L=0$, was first studied by Kulik and Omel’yanchuk
[14, 15]. It is easy to calculate from equations (4) and (5) the bound state
energies $E=\pm\Delta\cos(\phi/2)$ and the current

$J(L=0)=\frac{e}{\hbar}\Delta\sin\frac{\phi}{2}\tanh\frac{\Delta\cos(\phi/2)}{2T}.$
(7)

The current comes from the two bound states only as the continuum contribution
vanishes.

The second limiting case we consider is the infinite potential well. For an
ordinary potential, this model is well presented in books of quantum mechanics
because it has a fully analytic solution. In the case of the pair potential we
take the limit $\Delta\rightarrow\infty$. It is easy to calculate from
equations (4) and (3) that the bound states have energies $E_{n}=\frac{\hbar
v_{F}}{2L}[\pm\phi+(2n+1)\pi]$ and the current [8]

$\displaystyle
J(\Delta\rightarrow\infty)=\frac{2e}{\hbar}T\sum_{m=-\infty}^{\infty}\frac{\sin\phi}{\cosh(2\epsilon_{m}L/\hbar
v_{F})+\cos\phi}.$ (8)

The infinitely deep potential well is a good approximation for low-energy
levels, $E_{n}\ll\Delta$ in a long junction, $L\gg\hbar v_{F}/\Delta$ because
in this limit the bound state energies (2) become independent of $\Delta$.
This model is equivalent to the one used by Bardeen and Johnson [18] (except
that they use $L^{*}$ instead of $L$).

We study the current (8). At $T=0$ the sum goes to an integral. Evaluating
this gives

$\displaystyle J(\Delta\rightarrow\infty,T=0)=\frac{ev_{F}}{\pi
L}\phi\quad{\rm for}\ |\phi|<\pi$ (9)

and repeated periodically with period $2\pi$. This result was found by Ishii
(the linear dependence on $\phi$ in Ref. [5] and the coefficient in Ref.
[19]). We see that the linear dependence is rather well satisfied in the
selected examples of the general model, Figs. 2-4.

An interesting feature of the current (9) is that it corresponds to motion at
a velocity $v_{s}=(\hbar/2m)\phi/L$, which can be interpreted as the
superfluid velocity. Thus the normal region seems not to differ from a
superconductor, as long as the critical velocity $v_{c}=\pi\hbar/2mL$ is not
exceeded. Expressed in other words, the system behaves like it would be
Galilean invariant, although it is not because the locations of the N-S
interfaces are fixed [18].

The temperature dependence of the current can be calculated from (3) in the
general case. Examples in two limiting cases are shown in Figure 5.

Figure 5: The current-phase relations at different temperatures. (left) Zero-
length pair-potential well (7) evaluated using the BCS gap function
$\Delta(T)$ at reduced temperatures $T/T_{c}=0,0.1,\ldots,0.9$. (right)
Infinitely-deep pair-potential well (8) evaluated at temperatures 0, 1,…, 5 in
units of $\hbar v_{F}/2\pi L$. In both panels, an increasing temperature
corresponds to lower current, but the scales of the current and the
temperature differ essentially.

Let us consider the application of the square-well pair-potential model to
experimentally realizable systems. The main deficiency of the model in several
applications is that the pair-potential is not self-consistent. For example,
consider a normal section in a uniform-thickness superconducting wire. Because
ordinary scattering was neglected, the normal region induces a proximity
effect on the superconducting banks. As a result, the pair potential is
reduced in the banks near the normal region. This affects especially the
levels near the gap edge, which can penetrate deeper into the banks. There
are, however, some cases where the square-well model can be justified
quantitatively.

  * •

The limit of a long junction, $L\gg\hbar v_{F}/\Delta$. Here the current is
determined by low-energy levels, $E_{n}\ll\Delta$. The energies (2) of these
levels are not dependent on $\Delta$. Thus we expect this to be the case also
for $\Delta(x)$ modified by the proximity effect.

  * •

The case of a weak link separating two wider superconducting leads. In this
case the order parameter in the banks is determined by quasiparticle
trajectories staying in the leads with a minimal contribution from the
transmitting trajectories [14].

  * •

The case that the Fermi velocity in the banks is much larger than in the
normal region. This means that only a few of the trajectories incident on the
junction will be transmitted, and thus the proximity effect on the banks is
small [9].

Another possible point of worry in the square-well pair-potential model could
be that current conservation is strictly valid only in the normal region.
There is no current deep in the superconducting regions of constant phase.
This problem can be corrected by adding a spatially varying phase in the
superconducting regions. In the three cases where the proximity effect in the
banks can be neglected, the current in the normal region is always small
compared to the critical current in the superconducting regions. Thus the
induced phase gradient in the superconducting regions is small and does not
affect the calculation of the bound states and the current. This means that
current conservation is not a problem that is separate from the self-
consistency problem. Thus attempts to restore current conservation without
solving the self-consistency of the pair potential are misguided.

The square-well pair-potential model can be used as a qualitative model in
cases where its assumptions are not strictly valid. For example, consider a
flux line in a type II superconductor and a quasiparticle trajectory passing
the vortex line at distance $b$ (impact parameter). There is phase change
along the trajectory. We can model this using the square-well model. The phase
along the trajectory changes by $2\pi$ as a function of $b$ changing from
$-\infty$ to $+\infty$. There has to be a bound state whose energy changes
from positive to negative energy, and an opposite momentum state that shifts
from negative to positive energy. The zero crossing takes place as the
trajectory goes through the vortex line, $b=0$. This gives a qualitative
semiclassical description of the Caroli-de Gennes-Matricon vortex-core states.

The model can also be used to demonstrate the difference between an ideal
normal metal and a superconductor. Consider an occupation of the levels that
the positive momentum levels are filled up to an energy $E>0$ and
correspondingly, there is a deficit in filling the negative momentum levels.
Such a state would not decay if we strictly assumed the pure limit. The
current in such a state can be considerably higher than that obtained with
equilibrium occupation. This is not supercurrent since, for a filled positive
momentum level, there are empty negative momentum levels at lower energy.
Inelastic scattering can reduce the distribution and the current to their
equilibrium values.

There are a great number of generalizations of the square-well pair-potential
model. The most obvious is to allow for a self-consistent pair potential, see
Refs. [20, 21] for some early calculations. A second important generalization
is to consider the normal scattering of electrons. This couples the levels
with opposite momentum and results in that the crossings of the opposite
momentum levels at $\phi=n\pi$ (integer $n$) are replaced by avoided crossings
[8, 22, 23, 11, 24]. This is illustrated in Fig. 6. There are other
generalizations to non-equilibrium phenomena, different geometries,
unconventional superconductors/superfluids etc. Some of the topics are
considered in Refs. [12, 25, 26]. In fact, all studies of inhomogeneous
superconductors can be seen as generalizations of the square-well pair-
potential model.

Figure 6: Same as Fig. 4 but including barriers at both NS interfaces. The
noteworthy difference to Fig. 4 is the opening of gaps in the energy spectrum
at $\phi=0$ and $\pi$. The figure is based on the calculation by Galaktionov
and Zaikin [24]. The parameters are the transmission probability 0.9 at both
barriers, $L\Delta/\hbar v_{F}=5.7$, zero temperature, and the phase acquired
in one round trip between the two barriers equal to $\pi/2$. [Note the
misprint of an extra factor $\sin\chi$ in formula (16) of Ref. [24].]

I will finish this article with a personal view of how I learned and used the
square-well pair-potential model. First, an impurity in the vortex core acts
similarly as in a SNS junction at phase difference $\phi=\pi$: it increases
the superconductor’s condensation energy. This leads to a pinning force on the
vortex to the impurity [27, 28]. This is a single-impurity version of the
effect known earlier that impurity scattering makes the superconducting
coherence length shorter. Second, an impurity or a wall in a p-wave superfluid
acts like a SNS junction as the phase on the quasiparticle trajectory jumps in
scattering [29, 30, 31]. Third, we have studied both equilibrium and non-
equilibrium currents in pinholes, the analog of point contacts in
superconductors, of superfluid 3He [32, 33, 34]. This allows one to explain
$\pi$ states observed experimentally in an array of pinholes.

In summary, similarly to the role of the square-well potential model to learn
quantum mechanics, the square-well pair-potential model is an analytically
solvable model that allows us to learn inhomogeneous superconductivity and
fermion superfluidity.

## References

  * [1] S. Datta, Electronic Transport in Mesoscopic Systems (Cambridge, 1995).
  * [2] A.F. Andreev, J. Eksp. Teor. Fiz. 46, 1823 (1964) [Sov. Phys. JETP 19, 1228 (1964)].
  * [3] A.F. Andreev, J. Eksp. Teor. Fiz. 49, 655 (1965) [Sov. Phys. JETP 22, 455 (1966)].
  * [4] I.O. Kulik, Zh. Eksp. Teor. Fiz. 57, 1745 (1969) [Sov. Phys. JETP 30, 944 (1970)].
  * [5] C. Ishii, Prog. Theor. Phys. 44, 1525 (1970).
  * [6] A.V. Svidzinskii, T.N. Antsygina, and E.N. Bratus’, Zh. Eksp. Teor. Fiz. 61, 1612 (1971) [Sov. Phys. JETP 34, 860 (1972)].
  * [7] A.V. Svidzinsky, T.N. Antsygina, and E.N. Bratus’, J. Low Temp. Phys. 10, 131 (1973).
  * [8] A.I. Bezuglyi, I.O. Kulik and Yu.N. Mitsai, Fiz. Nizk. Temp. 1, 57 (1975) [Sov. J. Low Temp. Phys. 1, 27 (1975)].
  * [9] M.Yu. Kupriyanov, Fiz. Nizk. Temp. 7, 700 (1981) [Sov. J. Low Temp. Phys. 7, 342 (1981)].
  * [10] A. Furusaki, H. Takayanagi, and M. Tsukada, Phys. Rev. Lett. 67, 132 (1991).
  * [11] P.F. Bagwell, Phys. Rev. B 46, 12573 (1992).
  * [12] A.A. Golubov, M.Yu. Kupriyanov, and E. Il’ichev, Rev. Mod. Phys. 76, 411 (2004).
  * [13] E.V. Thuneberg, Phys. Rev B 108, 176501 (2023).
  * [14] I.O. Kulik and A.N. Omel’yanchuk, Fiz. Nizk. Temp. 3, 945 (1977) [Sov. J. Low Temp. Phys. 3, 459 (1977)].
  * [15] I.O. Kulik and A.N. Omel’yanchuk, Fiz. Nizk. Temp. 4, 296 (1978) [Sov. J. Low Temp. Phys. 4, 142 (1978)].
  * [16] A.L. Fetter and J.D. Walecka, Quantum theory of many-particle systems (McGraw-Hill, New York 1971).
  * [17] U. Gunsenheimer, U. Schüssler, and R. Kümmel, Phys. Rev. B 49, 6111 (1994).
  * [18] J. Bardeen and J.L. Johnson, Phys. Rev. B 5, 72 (1972).
  * [19] C. Ishii, Prog. Theor. Phys. 47, 1464 (1972).
  * [20] A. Levy Yeyati, A. Martin-Rodero, and F. J. Garcia-Vidal, Phys. Rev. B 51, 3743 (1995).
  * [21] R.A. Riedel, L. Chang, and P.F. Bagwell, Phys. Rev. B 54, 16082 (1996).
  * [22] W. Haberkorn, H. Knauer and J. Richter, Phys. Stat. Sol (a) 47, K161 (1978).
  * [23] C.W.J. Beenakker, Phys. Rev. Lett. 67, 3836 (1991).
  * [24] A.V. Galaktionov and A.D. Zaikin, Phys. Rev. B 65, 184507 (2002).
  * [25] W. Belzig, F.K. Wilhelm, C. Bruder, G. Schön and A.D. Zaikin, Superlattices and microstructures 25, 1251 (1999).
  * [26] Andreev bound states, Phil. Trans. R. Soc. A 376, Issue 2125, editors L.H. Greene and J.A. Sauls (2018).
  * [27] E.V. Thuneberg, J. Kurkijärvi, and D. Rainer, Phys. Rev. B 29, 3913 (1984).
  * [28] E.V. Thuneberg, J. Low Temp. Phys. 54, 415 (1984).
  * [29] E.V. Thuneberg, J. Kurkijärvi, and D. Rainer, J. Phys. C 14, 5615 (1981).
  * [30] W. Zhang, J. Kurkijärvi, and E.V. Thuneberg, Phys. Rev. B36, 1987 (1987).
  * [31] R. Hänninen and E.V. Thuneberg, Phys. Rev. B 67, 214507 (2003).
  * [32] J.K. Viljas and E.V. Thuneberg, Phys. Rev. B 65, 64530 (2002).
  * [33] J.K. Viljas and E.V. Thuneberg, Phys. Rev. Lett. 93, 205301 (2004).
  * [34] E.V. Thuneberg, AIP conference proceeding 850, 103 (2006).