# Precision and accuracy of acoustic gunshot location in an urban environment

Robert B. Calhoun<EMAIL_ADDRESS>Clark Dunson  Murphey L. Johnson
Scott R. Lamkin  William R. Lewis  Robert L. Showen  Mark A. Sompel  Lester P.
Wollman  ShotSpotter, Inc., Newark, CA 94560, USA

###### Abstract

The muzzle blast caused by the discharge of a firearm generates a loud,
impulsive sound that propagates away from the shooter in all directions. The
location of the source can be computed from time-of-arrival measurements of
the muzzle blast on multiple acoustic sensors at known locations, a technique
known as multilateration. The multilateration problem is considerably
simplified by assuming straight-line propagation in a homogeneous medium, a
model for which there are multiple published solutions. Live-fire tests of the
ShotSpotter gunshot location system in Pittsburgh, PA were analyzed off-line
under several algorithms and geometric constraints to evaluate the accuracy of
acoustic multilateration in a forensic context. Best results were obtained
using the algorithm due to Mathias, Leonari and Galati under a two-dimensional
geometric constraint. Multilateration on random subsets of the participating
sensor array show that 96% of shots can be located to an accuracy of 15 m or
better when six or more sensors participate in the solution.

##  I  Introduction

The muzzle blast resulting from discharge of a firearm generates a loud,
impulsive noise that propagates away from the shooter in all directions. The
location of the shooter can be estimated by minimizing the difference between
the measured and predicted arrival time of the muzzle blast on multiple time-
synchronized sensors, a technique known as acoustic multilateration. Public
safety agencies use networks of acoustic sensors to locate illegal gunfire in
near real-time, 10–15 seconds after discharge (Aguilar, 2015; Lawrence _et
al._ , 2018; Watkins _et al._ , 2002; National Law Enforcement and Corrections
Technology Center, 2007).

Acoustic data are also used in forensic context (Brustad and Freytag, 2005;
Maher and Routh, 2015; Maher and Hoerr, 2018; Beck, 2019; Begault _et al._ ,
2019). Acoustic multilateration cannot identify a shooter, but it can
establish the discharge location and discharge time for each individual shot
fired in a repeatable and objective manner. The efficacy and accuracy of an
acoustic multilateration system depends on the number and position of
reporting sensors, the accuracy to which the initial arrival time can be
measured, and the validity of the propagation model, with the simplest
acoustic propagation model assuming straight-line propagation. The present
work outlines a procedure for computing shot locations from acoustic
measurements of the muzzle blast using acoustic multilateration, and
quantifies the precision and accuracy of the results obtained using data
collected under controlled but realistic conditions.

The data derive from a series of live fire tests of the ShotSpotter gunshot
location system conducted in Pittsburgh, PA. These tests were analyzed offline
to identify the combination of algorithm and geometric constraint that offers
the greatest precision when locating gunshots in a forensic context.
Pittsburgh was selected as a study area for a number of reasons: first, a
large number of tests were conducted (three firearms discharged from nine
firing positions) covering neighborhoods with a density of structures similar
to that of many ShotSpotter coverage areas; second, the hilly terrain offered
the opportunity to investigate the importance of locating in three dimensions;
and finally because the Pittsburgh array features an unusually high sensor
density, allowing for simulation of the effects of deployed sensor density on
detection rates and location accuracy.

Source data from the Pittsburgh tests comparable to what is submitted as
evidence in a criminal trial are provided as supplemental materials. These
materials comprise time-stamped acoustic recordings, time- and location-
stamped output of the pulse processing algorithms, and tagged pulse sets
suitable for acoustic multilateration. These data permit other investigators
to review or extend the current work.

##  II Theory of Gunshot Location

###  II.1 Source signal

The muzzle blast of a gunshot is an impulsive acoustic source that emits power
over a broad range of frequencies. (Weber, 1939) developed an expression for
the pressure of spark explosions that has applicability to muzzle blasts; the
resulting expression is similar to that introduced by (Friedlander, 1946). The
single parameter of the Weber model is the Weber radius $r_{w}$, the radius at
which the supersonic combustion products decelerate to sonic velocity. Unlike
explosions, muzzle blasts are highly directional (Maher and Shaw, 2010; Beck
_et al._ , 2011) but they can be modeled as a Weber source in which the size
of $r_{w}$ depends on the angle between the gun barrel and the receiver
(Hirsch and Bertels, 2013).

Muzzle blasts are simple impulses (Maher, 2009) and recordings made in open
environments (such as Kruger National Park in South Africa, Figure 1, top)
show reasonable agreement with simple models (Don and Cramond, 1987; Embleton,
1996). In built-up areas, propagation distances are shorter and the received
waveforms are more complicated (Figure 1, bottom.) (Hornikx, 2016) notes that
“sound pressure levels in urban areas depend at one hand on the actual sound
sources, and at the other hand on the propagation of sound from these sources
in the environment”, a statement that is particularly applicable to urban
gunshots, where a simple source signal becomes a complicated received signal
as a result of the environment. In urban environments the acoustic energy of
the muzzle blast is transferred from source to a receiver via many different
paths, including specular and diffuse reflection (echoes and reverberation),
diffraction (transmission around structures), refraction (bending due to
speed-of-sound gradient), as attenuated by absorption in the atmosphere and at
surfaces (Kapralos _et al._ , 2008). It is possible to model acoustic
propagation in an urban environment using ray-tracing and other computational
modeling techniques (Le Bot and Bocquillet, 2000; Parker _et al._ , 2007;
Albert and Liu, 2010; Mehra _et al._ , 2014; Pasareanu _et al._ , 2018;
Remillieux _et al._ , 2012; Stevens and Murphy, 2014) but these models require
a structural model of the area and are in some cases computationally
demanding.

A simpler approach (Showen, 1997; Showen _et al._ , 2008) is to deploy more
than the minimum number of sensors required for multilateration and assume a
sufficient number of sensors will receive detectable levels of muzzle blast
energy via paths well-approximated by a straight-line propagation model. These
“near line-of-sight” pulses can be identified by constructing sets of pulses
from appropriate sensors, computing a prospective location from the pulse set,
and comparing arrival time predictions of the model with the measured arrival
time on each sensor. The goal of this combinatorial optimization problem is to
find the largest set of pulses that meet a mutual consistency test. Specular
reflections (echoes) are rejected because the path between source and each
sensor is distinct, resulting in different echo patterns in each received
signal; only pulses from “near line-of-sight” paths pass the mutual
consistency test. This principle is illustrated in Figure 2, which shows the
envelopes of three successive shots detected on five different sensors. Strong
echoes are detected on three of the sensors, but only the initial impulse from
each shot is consistent with a straight-line propagation model on all the
fives sensors. The mutual consistency method can only be used when the number
of reporting sensors exceeds the minimum number of time differences of arrival
(TDoAs) required for multilateration. The method works best when the sensor
density is high enough to ensure gunshots are detected on at least two more
sensors than required for multilateration, but not so high that the system
generates false triggers on unwanted sources of impulsive noise such as
construction noise.

Figure 1: Firearms discharged at similar ranges in an open environment (top:
Kruger National Park, South Africa) and in an urban environment (bottom: city
in the midwestern United States); the vertical scale is with respect to ADC
full scale of 93 dB SPL = 1.0. Figure 2: Waveform envelopes for a sequence of
3 rounds from a 9mm pistol as recorded by sensors at $79\text{\,}\mathrm{m}$,
$223\text{\,}\mathrm{m}$, $350\text{\,}\mathrm{m}$, $474\text{\,}\mathrm{m}$
and $775\text{\,}\mathrm{m}$ showing that range alone is not a good predictor
of signal strength. Echo patterns vary from location to location, allowing the
near line-of-sight paths to be identified through their greater self-
consistency. The signals have been aligned in the time domain; the vertical
scale is the same for all plots, with amplitude of 1.0 corresponding to ADC
full scale of 93 dB SPL.

###  II.2 Pulse localization

In built environments, the acoustic energy from a muzzle blast arrives via
multiple paths. Energy transmitted via the shortest acoustic path arrives
first; additional energy from the blast is delivered via reflections and other
indirect paths over a period of time that may extend to several seconds. In
the present work, signal refers to the acoustic pressure change induced by the
muzzle blast as measured by a remote microphone as a function of time; the
goal of the pulse detection algorithm is to identify the arrival time of
pulses—fast rise time portions of the measured signal—that are plausibly due
to muzzle blast energy arriving via the shortest acoustic path. Because
shortest acoustic paths are often not line-of-sight paths, significant
attenuation may occur and the acceptance criteria for muzzle blast impulses
must necessarily be broad.

One effective strategy for detecting plausible muzzle blasts in built
environments is to convolve the envelope of an acoustic signal with a step-
detection kernel. This acknowledges that muzzle blast power will be broadly
distributed in the time domain, while giving the most weight to the earliest-
arriving signal.

Let $f(t)$ be the amplitude of an acoustic signal as a function of time $t$.
The instantaneous magnitude (or envelope) of $f(t)$ is given by the magnitude
of the analytic associate of $f(t)$

$\|f\|=\|f(t)+i\mathcal{H}(f(t))\|$

where $\mathcal{H}(f(t))$ is the Hilbert transform of $f(t)$ (Boashash, 2015).
Define a weighted step-detection kernel $g(t)$ as

$g(t)=\begin{cases}t/\tau+1&\textrm{for\quad}-\tau<t\leq 0\\\ \
-1/2&\textrm{for\quad}0<t\leq\tau\\\ 0&\textrm{otherwise}\\\ \end{cases}$ (1)

where $\tau$ is a length scale over which the signal should be averaged.

Define $x(t)$ as the convolution the signal envelope $\|f(t)\|$ and kernel
$g(t)$:

$x(t)=(\|f\|*g)(t).$ (2)

The pulse arrival time is defined as the time $t$ at which $x(t)$ is locally
maximized. Pulses with values of $x(t)$ that fall below a threshold are
discarded.

The kernel $g$ gives as much weight to the absence of signal prior to pulse
arrival as it does to the amplitude of the signal that follows; this provides
consistent estimates of the arrival time of the most direct (earliest-
arriving) path, even when the energy received via that path is a small
fraction of the total. Increasing the length scale $\tau$ reduces spurious
triggering on echoes and reverberation, but an excessively long window reduces
the ability to detect closely-spaced shots. Empirically, a value of
$50\text{\,}\mathrm{ms}$ has been found to work well on both isolated shots
and the range of semi-automatic (300 rpm) and fully automatic (450–950 rpm)
firearms encountered in a public-safety context. Note this is 10–20 times the
duration of a muzzle blast measured on an open range (Beck _et al._ , 2011;
Maher, 2011; Routh and Maher, 2016), which is on the order of
$2\text{\,}\mathrm{ms}$.

Figure 3: Pulse arrival time (red dot) plotted on the original waveform for
sensor at $223\text{\,}\mathrm{m}$ distance from Figure 2. Dotted line shows
output of pulse localization function $x(t)$ from (2). The full scale
amplitude approximately 93 dB SPL.

###  II.3 Pulse Selection

In urban environments, true line-of-sight paths are rare. “Near line-of-sight”
paths, where the received signal is line-of-sight after a single diffraction
over a barrier, can be used for location because the typical height of
buildings is small compared with the distances between acoustic sensors. Use
of non-line of sight paths with significant time delays, such as those
generated by specular reflections from buildings and other hard surfaces, will
cause significant location error if used in a multilateration calculation.
Associating pulses from one shot with pulses from a different shot will also
result in large multilateration error.

The construction of sets of pulses that are appropriate for multilateration
becomes more difficult as the number of shots and number of shooters
increases. ShotSpotter’s pulse selection algorithms use a combination of
combinatorial optimization, pattern-matching and numerical optimization
methods to search for the largest possible set of mutually-consistent pulses
for each shot. The selection of appropriate pulses for multilateration is an
involved topic (Calhoun, 2011; Showen _et al._ , 2008; Beck, 2019) and a full
discussion of the methods used is outside the scope of the present work. To
facilitate review of the present multilateration results, the output of the
ShotSpotter pulse selection routines (comprising sets of acceptably “near line
of sight” pulses for each shot) are provided as part of the supplemental
materials, along with sufficient audio data to verify that the selected pulses
are associated with the shortest acoustic path. (Supplementary materials
containing audio files, pulse properties, and tuples of near line-of-sight
pulses for each shot that are suitable for multilateration are currently
available at https://github.com/ShotSpotter/research.accuracy-of-gunshot-
location.) Because the search for pulse sets is not exhaustive, the tuples
provided are not guaranteed to be an optimal subset for each test. However,
all pulses in a tuple will have an error of less than $40\text{\,}\mathrm{ms}$
between the measured arrival time and the predictions of a straight-line
propagation model.

###  II.4 Source Localization

Acoustic multilateration, the location of a source from multiple time
difference of arrivals (TDoAs), has been in use since World War I when French
(Aubin and Goldstein, 2014), Canadian (Finan and Hurley, 1997) and British
forces (Van der Kloot, 2005) used sound ranging to locate the source of
artillery fire. Many multilateration algorithms have been documented since
that time. A multilateration solution in $d$ dimensions requires at least
$n=d+1$ reporting sensors. Some $n=d+1$ solutions have two mathematically
consistent solutions. Additional methods not discussed in the present work,
such as measurements of angle-of-arrival (AoA) or received signal power, are
required to resolve ambiguous solutions, so the present work focuses on
solutions for which $n\geq d+2$, for which the multilateration solution is
always unique given non-degenerate sensor placement.

The TDoA source location problem reduces to finding the point of intersection
of multiple hyperboloids, so the algorithms comprise variations on a common
theme. ShotSpotter uses multiple algorithms because they exhibit different
numerical stability and sensitivity to the choice of a reference sensor. Below
we analyze the accuracy of four algorithms and three constraints under a
variety of array densities to compare their relative performance under real-
world test conditions. A goal of the present work is to identify which
algorithm is best-suited for forensic work. All solutions discussed assume
line-of-sight propagation in a homogeneous, stationary medium.

####  II.4.1 Reddi

This implementation is based on (Reddi, 1993). With reference sensor at
$(0,0,0)$ and other sensors at $(x_{i},y_{i},z_{i})$, $i=1,2,3,\dots n$,
source location $(x_{*},y_{*},z_{*})$ can be found from relative arrival times
$t_{i}$ by equating the squared distance from the reference to the source with
the square distance of the reference to the other sensors:

$\displaystyle ct_{i}=S_{i}-\sqrt{x^{2}+y^{2}+z^{2}}$ (3)

for $i=(1,2,\dots n)$ where

$S_{i}^{2}=(x-x_{i})^{2}+(y-y_{i})^{2}+(z-z_{i})^{2}.$

This reduces to a quadratic equation that can be solved for source location
$(x_{*},y_{*},z_{*})$. Discharge time $t_{*}$ is obtained from the mean value
of $t_{i}-S_{i}/c$ at $(x_{*},y_{*},z_{*})$. This implementation requires a
minimum of $d+1$ TDoAs, where $d$ is the number of dimensions.

Figure 4: Graphical representation of the intersecting hyperbola approach
characteristic of Reddi algorithm. The reference sensor that defines the time
differences of arrival is the first reporting sensor. The sample data is from
test 1 at Firing Position 4.

####  II.4.2 MLG

(Mathias _et al._ , 2008) avoid use of a reference sensor by minimizing the
squared error on each sensor. ShotSpotter’s “MLG” implementation slightly
modifies this algorithm by solving for time-distance $(ct)$ instead of $t$,
which improves computational stability. For the simplified case of
overdetermined systems ($n>d+1$) with receiver matrix

$\mathbf{R}=\begin{bmatrix}x_{0}&x_{1}&\dots&x_{i}\\\
y_{0}&y_{1}&\dots&y_{i}\\\ z_{0}&z_{1}&\dots&z_{i}\\\
-ct_{0}&-ct_{1}&\dots&-ct_{i}\\\ -c/2&-c/2&\dots&-c/2\\\ \end{bmatrix},$ (4)

the solution vector is given by

$\begin{bmatrix}x\\\ y\\\ z\\\ ct\\\ v/c\\\
\end{bmatrix}=\frac{1}{2}\mathbf{R}^{\dagger}\begin{bmatrix}x_{0}^{2}+y_{0}^{2}+z_{0}^{2}-(ct_{0})^{2}\\\
x_{1}^{2}+y_{1}^{2}+z_{1}^{2}-(ct_{1})^{2}\\\ \dots\\\
x_{i}^{2}+y_{i}^{2}+z_{i}^{2}-(ct_{i})^{2}\\\ \end{bmatrix}$ (5)

where $v$ is an error term and $\mathbf{R}^{\dagger}$ is the pseudoinverse of
$\mathbf{R}$. See (Mathias _et al._ , 2008) for handling of the fully
determined case ($n=d+1$) , and for a quantification of error propagation in
multilateration.

####  II.4.3 LeastSquares

The “LeastSquares” implementation is another least-squares minimizer, based on
time difference of arrivals between pairs of sensors. Equating the squares of
the distance and the time-distance gives one equation for each sensor $i$:

$(x_{i}-x)^{2}+(y_{i}-y)^{2}+(z_{i}-z)^{2}=(ct_{i}-ct)^{2}.$ (6)

Subtracting the $i$th equation from the ${(i+1)}$-st eliminates quadratic
terms. Define $X_{i}$ to be the difference between the LHS (distance) and RHS
(time-distance) of the $i$th difference equation. When

$\chi^{2}=\sum_{i=1}^{n-1}X_{i}^{2},$

the minimum value of $\chi^{2}$ is the point at which

$\frac{\partial\chi^{2}}{\partial x}=\frac{\partial\chi^{2}}{\partial
y}=\frac{\partial\chi^{2}}{\partial z}=\frac{\partial\chi^{2}}{\partial
ct}=0.$

This gives a linear system in $d+1$ variables:

$\begin{bmatrix}a_{i}a_{i}&a_{i}b_{i}&a_{i}c_{i}&a_{i}d_{i}\\\
b_{i}a_{i}&b_{i}b_{i}&b_{i}c_{i}&b_{i}d_{i}\\\
c_{i}a_{i}&c_{i}b_{i}&c_{i}c_{i}&c_{i}d_{i}\\\
d_{i}a_{i}&d_{i}b_{i}&d_{i}c_{i}&d_{i}d_{i}\\\
\end{bmatrix}\begin{bmatrix}x\\\ y\\\ z\\\ ct\\\
\end{bmatrix}=\begin{bmatrix}a_{i}e_{i}\\\ b_{i}e_{i}\\\ c_{i}e_{i}\\\
d_{i}e_{i}\\\ \end{bmatrix}$ (7)

where

$\displaystyle a_{i}$ $\displaystyle=(x_{i+1}-x_{i}),$ $\displaystyle b_{i}$
$\displaystyle=(y_{i+1}-y_{i}),$ $\displaystyle c_{i}$
$\displaystyle=(z_{i+1}-z_{i}),$ $\displaystyle d_{i}$
$\displaystyle=-c(t_{i+1}-t_{i}),\>\textrm{and}$ $\displaystyle e_{i}$
$\displaystyle=\dfrac{1}{2}(x_{i+1}^{2}+y_{i+1}^{2}+z_{i+1}^{2}+c^{2}t_{i}^{2}-x_{i}^{2}-y_{i}^{2}-z_{i}^{2}-c^{2}t_{i+1}^{2}).$

This requires an additional TDoA compared with the Reddi and MLG
implementations.

####  II.4.4 Iterative Discharge Time

Combustion of the propellant used firearms generates an optical signal known
as the muzzle flash that is strongest in the medium-wave IR region (Pauli _et
al._ , 2004). If the discharge time $t_{*}$ is known from muzzle flash timing
(6) reduces to finding the point of intersection of circles (or spheres) of
known radius $c(t_{i}-t_{*})$. This is a simplified version of (7).

When no muzzle flash time is available the discharge time can be determined
iteratively. The “Iterative Discharge Time” (IDT) implementation uses simplex
gradient descent to find $t_{*}$ based on starting value of $t_{*}=t_{0}-1.0$,
where $t_{0}$ is the arrival time on the nearest (earliest) sensor and the
constant roughly equal to the mean distance between a shooter and the first
reporting sensor. The main advantage over pure TDoA algorithms is that the
range can be bounded, which is helpful when locating sources outside the
sensor array. An additional benefit of all intersecting circle solutions in a
forensic context is that the location tuple obtained $(x_{*},y_{*},t_{*})$ can
readily be checked using the standard tools of plane geometry, a ruler and a
compass. A graphical representation of an intersecting circles solution is
shown in Figure 5.

Figure 5: Graphical representation of the intersecting circle approach
characteristic of Iterative Discharge Time algorithm. The sample data is from
test 1 at Firing Position 4.

###  II.5 Geometric Constraints

Acoustic multilateration can be done in two or three dimensions. For two-
dimensional multilateration, the sensor positions are projected against a
plane and the multilateration solution returned is computed in this plane. At
least three TDoAs are required. For three-dimensional multilateration, the
sensor elevation data is used without modification and the multilateration
solution computed in three-dimensional space. Four TDoAs are required for
multilateration in three dimensions.

There are two potential advantages of three-dimensional multilateration.
First, a 3D straight-line distance propagation model can be used. This is
potentially more accurate. Second, the elevation of the shooter can be
estimated. The downside of 3D multilateration is that the additional degree of
freedom allowed by three-dimensional multilateration can lead to less accurate
2D results if the sensor position matrix is poorly conditioned, e.g. the range
of vertical elevations of reporting sensors is small compared with the
horizontal spacing between them. With poorly-conditioned arrays, small errors
in the measurements of arrival time, sensor position, or speed of sound may
result in large changes in the shooter’s putative elevation.

There are plausible use cases for detecting gunshots above ground level,
including locating shots fired from the upper floors of a structure and
assisting in the classification of aerial fireworks. Nevertheless the vast
majority of firearms are discharged at ground level, so a practical approach
is to do the bulk of pulse selection in two dimensions, and then refine the
solution using a 3D model when sufficient data are available.

An approach that enables use of a 3D propagation model while avoiding the
instability of a full 3D solution is to constrain the shooter location to a
two-dimensional plane. This “2.5D” method (Calhoun, 2013) refines an initial
2D solution by assuming the shooter is standing on the ground (or on the roof
of a building) and constructing $\mathbf{R}$ with a duplicate set of sensors
with identical arrival time but position reflected through the horizontal
plane containing the shooter, Figure 6. The method requires use of a digital
elevation and/or building model to provide an estimate for $z_{*}$, the
$z\textrm{-elevation}$ of the shooter. With ground plane normal to $\hat{z}$,
each TDoA tuple $(x_{i},y_{i},z_{i},t_{i})$ in (4) is supplemented with a
second TDoA tuple $(x_{i},y_{i},2z_{*}-z_{i},t_{i})$. This doubles the number
of columns in (4) and the number of rows in (5) but no changes in algorithm
are required. The solution will be constrained to elevation $z_{*}$ by the
symmetry of the input data.

Figure 6: Method for constructing a sensor array through a ground plane by
reflection. This permits use of elevation data from the remote sensors while
constraining the output to the specified plane.

###  II.6 Weather Effects

Atmospheric conditions affect sound propagation and therefore the accuracy of
an acoustic gunshot location system. The speed of sound varies only slowly
with temperature, and little with humidity or pressure, so the dry air
relationship for the speed of sound $c$ may be used with acceptable accuracy.
This is:

$c=20.03\sqrt{T+273.15}\quad$\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$ (8)

where $T$ is the temperature in $\mathrm{\SIUnitSymbolCelsius}$.

Sound energy is transmitted as a compression wave in a fluid; a moving fluid
carries the sound energy with it. The multilateration routines of Section II.4
assume a stationary medium, but the effects of a homogeneous moving fluid can
be incorporated by shifting the apparent position of each sensor to account
for air motion from when the shot is fired to when it is detected. With the
$z$-axis component of wind vector $\vec{v}$ assumed to be zero, the wind-
corrected position of the sensors is given by:

$\begin{bmatrix}x_{i}^{\prime}\\\
y_{i}^{\prime}\end{bmatrix}=\begin{bmatrix}x_{i}\\\
y_{i}\end{bmatrix}-(t_{i}-t_{*})\begin{bmatrix}v_{x}\\\ v_{y}\end{bmatrix}$
(9)

where $t_{i}$ the arrival time at each sensor and $v_{x}$ and $v_{y}$ are the
$x$\- and $y-$components of the wind vector. The discharge time $t_{*}$ is not
known a priori but it can be estimated by subtracting an appropriate value
from the arrival time at the first reporting sensor. A suitable value is the
(historic) mean distance from the shooter to the first reporting sensor; a
value of $1.0\text{\,}\mathrm{s}$ is appropriate for arrays with sensor
densities in the range of
$2\text{\,}{sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}12\text{\,}{sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}$.
The wind-correction process can be repeated with a more accurate value after
initial computation of discharge time $t_{*}$.

Figure 7: Adjustment of sensor positions to take homogeneous fluid movement
into account. Energy is propagated from source to sensor via a combination of
bulk fluid flow and acoustic propagation. The effects of bulk fluid flow can
be incorporated by shifting the apparent position of the sensor to the
position that would have received the signal in the absence of fluid flow.

Other deleterious effects of wind are not so easily handled. Wind increases
the background noise detected by outdoor microphones (Raspet _et al._ , 2006),
an effect that can be ameliorated but not entirely eliminated with proper
windscreens (Hosier and Donavan, 1979; Lin _et al._ , 2014; Schomer _et al._ ,
1990; Wuttke, 1992). Wind noise reduces the signal-to-noise ratio of gunshot
impulses, making them harder to detect. Gradients in wind speed and
temperature with elevation refract the sound, bending it upwards in the upwind
direction and downwards in the downwind direction (Wiener and Keast, 1959;
Wilson, 2003; Franke _et al._ , 1988). Upwind refraction creates a “shadow
zone” where the gunshot signal is not detectable. Because of the potential
lack of participation from upwind sensors during windy conditions, additional
sensors must be to deployed to enable successful multilateration in all
weather conditions. The ShotSpotter system does not otherwise compensate for
the effects of wind or temperature gradients.

##  III System Design

###  III.1 Array design

A typical deployment density is
$2\text{\,}{sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}12\text{\,}{sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}$,
with the density dependent on the density of structures and foliage in the
covered area. The low end of this range is used in open areas where few
structures impede acoustic propagation, such as anti-poaching deployments;
higher sensor densities are used in areas with a high density of structures,
high background noise (traffic, rapid transit systems), difficult wind
conditions, or unfamiliar environments. Sensors are installed on elevated
locations such as rooftops and utility poles to increase the number of near
line-of-sight audio paths, and to minimize background noise. Wind noise is
mitigated with a $3.5\text{\,}\mathrm{cm}$ layer of open-cell foam around the
microphones.

###  III.2 Sensor processing

All signal processing and pulse detection operations are performed in the
field on remotely-deployed single-board computers based on the Texas
Instruments AM335x platform. Sensor position is provided by a u-blox MAX 7 GPS
chipset, and the Linux system clock is disciplined using the PPS signal from
the GPS. The u-blox MAX-7 series of GPS chipsets have a datasheet location
accuracy of 2.5m (50% circular error probability (CEP)) and timepulse accuracy
$30\text{\,}\mathrm{n}\sec$ (RMS). Actual performance depends on the satellite
constellation visible from the installed sensor position. Because sensors are
rarely moved, averaging over long periods of time is a feasible method for
improving location precision. Each hour, the mean position over the previous
hour is stored, and a long-term average position re-estimated based on values
from up to 180 days of hourly records, (Figure 8) with the multilateration
calculations making use of the long-term averages. It is expected that the
increasing availability of high-accuracy real-time kinematic (RTK) GPS
chipsets will reduce or eliminate the need for position averaging in future
multilateration systems.

Figure 8: Horizontal GPS performance of a typical stationary sensor, showing
the distance from the centroid of GPS measurements sampled once per minute, of
the hourly averages, and of the long-term average of hourly values. 50% CEP is
slightly worse than predicted by the datasheet but averaging reduces the RMS
distance from the centroid to less than $0.02\text{\,}\mathrm{m}$.

Audio sound pressure level is measured using Kingstate KECG2742WBL waterproof
electret microphones and digitized using a Cirrus Logic WM8737 analog-to-
digital converter, which samples 24-bit audio data at 48 kHz. Microphone
sensitivity is $-42\pm 3$ dB SPL, giving a full-scale digital output of
approximately 93 dB SPL. Network connectivity is provided by an on-board
cellular data module.

For performance reasons, the pulse localization algorithm described in II.2 is
implemented as a discrete-time algorithm. The digitized input stream is split
up into 2048-length chunks with 50% overlap so that a Hann window can be
applied without data loss. After windowing, data is converted into the
frequency domain for filtering and analysis. A typical pulse localization
response is shown in Figure 3.

system |  $\sigma$ ($\mathrm{\SIUnitSymbolMicro s}$)
---|---
GPS timepulse | $3\text{\times}{10}^{-2}$
Linux system clock | $1\text{\times}{10}^{-3}$
Pulse start measurement | $1\text{\times}{10}^{2}$
Table 1: Scale of contributions to pulse arrival time error

Pulse arrival time measurements are derived from audio synchronized to the GPS
PPS timepulse. As shown in Table 1, the timing errors associated with the GPS
system are negligible compared with the uncertainty in pulse arrival time
measurement. The net effect of all timing errors is estimated by running the
pulse detection algorithm on gunshot-like impulsive signals synchronized with
GPS PPS. This measurement takes into account multiple sources of error. The
measured timing error through the ShotSpotter sensors has a mean of
$\mu=$-316\text{\,}\mathrm{\SIUnitSymbolMicro s}$$ with a standard deviation
of $\sigma=$99.6\text{\,}\mathrm{\SIUnitSymbolMicro s}$$, with a typical
measurement shown in Figure 9. The pulse arrival time data provided and used
for multilateration in the current work are to millisecond precision.

Figure 9: Pulse start timing measurements using simulated gunshot impulses
triggered by GPS timepulses

The audio processing code analyzes the first $250\text{\,}\mathrm{ms}$
following pulse arrival time to determine the likelihood that the impulse is
due to a gunshot muzzle blast. Measured characteristics of the impulse include
the integrated power of the impulse, the shape of its envelope, and a
normalized power spectrum. These data are transmitted to a common endpoint for
the multilateration computation.

###  III.3 Classification

The original 48 kHz data is discarded after pulse processing; data downsampled
to 12 kHz is cached on the sensor for 30 hours to support forensic analysis.
Following detection of a potential gunshot, a 2-second subset of downsampled
(12 kHz) FLAC-encoded audio is retrieved from the two closest sensors and used
to build an image mosaic (Figure 10) comprising features characteristic of the
gunshot audio on each sensor and features characteristic of the incident as a
whole. A full description of image mosaic construction is outside of the scope
of the present work, but the tiles comprising the image can be summarized as
follows: A) waveform image; B) discrete wavelet spectrogram of the de-noised
waveform; C) pulse time-of-arrival markers aligned using the shooter-sensor
time delay computed via multilateration; D) discrete wavelet spectrogram of
the noise (complement of tile B above); E) wavelet transform of the first
impulse; F) waveform of the first impulse; G) participating area, as defined
by the Voronoi diagram of participating sensors; H) location of recent nearby
incidents I) plots of various per-sensor parameters as a function of distance;
J) version-tracking code; K) time-of-date counters; L) individual shot
locations; M) power spectrum; N, O, P) pulse feature counters; Q, R) recent
incident counters; S) logo. Discrete wavelet transforms are performed using
the biorthogonal 2.2 wavelet and scaling functions.

Classification of the mosaic images is performed using the ResNet (He _et al._
, 2016) image classifier trained using field-collected data labelled by human
reviewers. Each mosaic image (one from the nearest sensor, and one from the
second-nearest) is classified independently; the outputs of the two
classifiers are combined to minimize the likelihood of false negatives. In the
vast majority of cases, “ground truth” (e.g. shell casings, video recordings)
is not available, and it is to be expected that some training data are
misidentified.

Figure 10: The top image shows a typical mosaic image used for machine
classification with ResNet image classifier. This example is from a
$9\text{\,}\mathrm{mm}$ firearm discharged at Firing Position 8. Each image
used for classification is comprised of 19 individual tiles; some tiles are
derived from a 2-second subset of 12 kHz audio, while others are derived from
the characteristics of the incident as a whole. See text for tile
descriptions.

Incidents classified as gunfire are pushed to the public safety agency’s
dispatch system and/or mobile devices. All of the pulse data associated with
an incident is stored in a database so that the measurements made by the
remote sensor are available for later analysis. If a forensic report is
requested (typically in the course of a criminal investigation) within 24
hours, eight seconds of 12 kHz downsampled audio may be retrieved from sensors
participating in the incident. The downsampled audio is invaluable for
verifying that the pulse selection algorithms have grouped the input pulse
data into appropriate subsets for multilateration.

###  III.4 Weather Measurements

Reliable temperature data with a precision of $\pm
1\text{\,}\mathrm{\SIUnitSymbolCelsius}$ are available via NWS METAR readings
from major airports; measurements from Pittsburgh International Airport (KPIT)
and Allegheny County Airport (KAGC) were used in the present work. Wind is
highly dependent on local conditions and sensor mounting elevation so the
airport values can provide only a rough estimate of local conditions, but both
airports reported negligible wind of
$0\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$
during the live fire test, so the effects of wind were not investigated.

##  IV Description of Live Fire Tests

ShotSpotter sensor arrays are routinely tested prior to active use by law
enforcement. The live fire test of the ShotSpotter system examined here was
conducted in Pittsburgh, PA by the Pittsburgh Bureau of Police in December,
2018. Pittsburgh has mixture of hilly and flat terrain. The regions covered by
these tests primarily contain detached structures one to three stories in
height, comparable to the scale of structures in many other American cities.
Most sensors are installed on rooftops; the major component of the variance in
sensor elevation is due to the hilly terrain. Several of the firing positions
were adjacent to institutional buildings that generated loud, sharp echoes.

Officers secured the target area prior to each test for safety reasons. Live
rounds were discharged into a bullet trap at distance $\approx$
$3\text{\,}\mathrm{m}$ from the firearm, Figure (11). All shots at a giving
firing position were discharged with the muzzle of the firearm pointed in the
same direction, which was not otherwise specified. The survey location of each
firing position was determined by ShotSpotter personnel during the test using
a smartphone application and qualitatively verified using aerial photographs
available on Google Earth. Survey error $\epsilon$ was not well-characterized
during the test but typical performance of such devices (which make use of
signals from cellular and WiFi networks in addition to GPS) when used outdoors
with good sky view is on the order of $\epsilon=$ $2.5\text{\,}\mathrm{m}$.

Figure 11: Firing Position 1, showing the shooter and bullet trap, illustrating the source of the signals used for the test. FP | description | $N/M$ | $n$ | $\epsilon$ | $\rho$
---|---|---|---|---|---
|  |  |  | ($\mathrm{m}$) | (${sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}$)
FP1 | overpass ~$35\text{\,}\mathrm{m}$ above grade | 36/36 | 17.8 | 4.14 | 9.4
FP2 | parking lot near top of hill | 36/36 | 14.3 | 5.20 | 11.4
FP3 | intersection on high bluff | 36/36 | 9.9 | 3.49 | 7.5
FP4 | intersection near top of hill | 36/36 | 13.8 | 4.76 | 9.4
FP5 | playing field at head of valley | 33/36 | 11.6 | 2.21 | 8.5
FP6 | playing field near top of hill | 36/36 | 13.6 | 5.54 | 8.5
FP7 | playing field at end of valley | 36/36 | 14.9 | 7.88 | 12.4
FP8 | intersection at bottom of valley | 36/36 | 11.2 | 2.18 | 9.9
FP9 | parking lot on sloped region | 36/36 | 10.6 | 6.93 | 8.0
Table 2: Realtime performance at the nine firing positions (FP) analyzed in
present study. Buildings at most sites comprised detached residential
structures 1–3 stories in height, with some larger commercial, residential,
and institutional buildings. Ratio $N/M$ is number of shots $N$ detected out
of shots fired $M$; $n$ is mean participating sensor count per shot; survey
error $\epsilon$ is the distance from the acoustically-determined centroid to
the survey location determined by smartphone-assisted GPS. Sensor density
$\rho$ is in ${sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}$. This is
2D performance using the algorithm matching the most sensors for each shot.
Figure 12: Cumulative distribution function of survey error during live-fire
testings. For live fire incidents, 3-shot bursts are detected and scored as a
single incident. Figure 13: Sensor positions shown relative to Firing
Positions 1–9, which are located at the center of each plot. The filled
circles show those sensors that participated in at least one shot; the dashed
“detection limit” is the convex hull of those sensors reporting in at least
50% of shots.

Testing followed the standard firing sequence used by ShotSpotter. This
comprises twelve shots (three single shots and three groups of three shots)
for each firing position and tested weapon. This test used three different
handguns at each firing position (.45 cal, .40 cal, $9\text{\,}\mathrm{mm}$).
The test sequence was repeated at each of nine firing positions, for a total
of 162 tests comprising 324 shots. A general description of each firing
position is given in Table 2. Temperature during the test was $-1\pm
1^{\circ}\textrm{C}$ with light (0–2 m/s) wind. During the real-time tests,
the system detected 100% of incidents and located 96.9% of incidents within a
$25\text{\,}\mathrm{m}$ benchmark criterion. The cumulative distribution
function of survey error below $25\text{\,}\mathrm{m}$ for all tests is shown
Figure 12. The real-time ShotSpotter system bundles multiple closely-spaced
shots into a single incident, so scoring is per-test rather than per-shot.

Following each activation, the system requests two-second samples of 12 kHz
audio to facilitate machine classification. The results of the on-sensor pulse
analysis are stored in a database for later use. Following the test,
ShotSpotter personnel retrieved up to eight seconds of additional audio from
selected sensors located within $800\text{\,}\mathrm{m}$ of the firing
position. The range limit is in alignment with the audio retrieval policy
adopted by the company when servicing forensic requests. The present study
limits performance analysis to sensors located within an
$800\text{\,}\mathrm{m}$ radius of each firing position.

Selected data from these tests are available as supplementary materials. These
materials include pulse data (impulse arrival time, signal-to-noise ratio,
acoustic power and sensor position) and up to six 12 kHz audio files for each
test. Courts of law routinely redact sensor locations from trial records
because of the risk of retaliation against property owners; accordingly, the
sensor locations in the dataset are provided relative to an arbitrary (and
unpublished) origin that is different for each firing position.

The pulse data include a record of impulsive sounds triggered during the time
range around each test. These include muzzle-blast sounds with near line-of-
sight propagation paths that are useful for multilateration location as well
as multipath signals and other impulsive noise sources that are unrelated to
the gunfire test. ShotSpotter uses various pattern-matching and combinatorial
consistency algorithms to pick sets of pulses suitable for multilateration.
While these routines are not revealed in the present work, the output of the
pulse selection algorithms (as tuples of pulses for each shot of each test)
are provided as supplementary materials.

##  V Analysis

###  V.1 Multilateration Implementations

The nine firing positions (Table 2) were re-analyzed offline to evaluate the
performance of the four 2D implementations available (Reddi, LeastSquares,
MLG, and IDT). The input data were tuples of shortest acoustic path pulses
within $800\text{\,}\mathrm{m}$. The mean number of participating sensors per
shot ranged from 9.6 to 17.5, well over the four needed for unambiguous 2D
multilateration. Results by implementation are shown in Table 3. The survey
error $\epsilon$ is the two-dimensional distance from the acoustic
multilateration centroid to the survey location determined by smartphone GPS;
$\sigma_{1}$ and $\sigma_{2}$ are the square roots of the eigenvalues of the
two-dimensional covariance matrix. Two-dimensional distance is considered the
appropriate error measure in a public safety context.

All four routines solve the same mathematical problem—minimizing the
difference between the predicted and measured arrival times at the several
sensors assuming a straight-line propagation model—so unsurprisingly similar
results are obtained from all four algorithms. The Reddi implementation uses
the first reporting sensor as a reference, which gives more weight to the
nearest sensor than to all of the rest. The LeastSquares and related IDT
routines yield good accuracy and stability, but the LeastSquares requires an
extra TDoA and the IDT routine has a running time roughly ten times longer
than Reddi or MLG for the same input data because it is iterative. The MLG
implementation requires only the minimum number of TDoAs for a solution, does
not arbitrarily overweight a single sensor, works in two or three dimensions
and also has the fastest running time, since the computation required to set
up $\mathbf{R}$ is simple and $\mathbf{R}^{\dagger}$ can be computed using
SVD. The MLG implementation was used in the remainder of the study.

Algo | FP | $N/M$ | n | $\epsilon$ | $\sigma_{1}$ | $\sigma_{2}$
---|---|---|---|---|---|---
|  |  |  | $(\textrm{m})$ | $(\textrm{m})$ | $(\textrm{m})$
Iterative | FP1 | 36/36 | 17.5 | 3.14 | 0.62 | 0.41
Discharge | FP2 | 36/36 | 17.5 | 5.14 | 1.40 | 0.25
Time | FP3 | 36/36 | 17.5 | 0.61 | 1.40 | 0.83
(IDT) | FP4 | 36/36 | 17.5 | 4.09 | 2.33 | 0.56
| FP5 | 36/36 | 14.0 | 3.25 | 1.71 | 0.93
| FP6 | 36/36 | 14.0 | 5.21 | 1.29 | 0.40
| FP7 | 36/36 | 13.9 | 7.17 | 0.94 | 0.39
| FP8 | 36/36 | 14.0 | 2.67 | 1.39 | 0.47
| FP9 | 36/36 | 9.9 | 6.58 | 0.66 | 0.55
Reddi | FP1 | 36/36 | 9.9 | 3.55 | 1.03 | 0.62
| FP2 | 36/36 | 9.8 | 4.66 | 1.05 | 0.33
| FP3 | 36/36 | 9.9 | 3.56 | 2.40 | 0.74
| FP4 | 36/36 | 13.8 | 7.69 | 5.52 | 1.86
| FP5 | 36/36 | 13.6 | 2.13 | 6.60 | 0.77
| FP6 | 36/36 | 11.2 | 6.30 | 0.76 | 0.50
| FP7 | 26/36 | 13.6 | 18.37 | 5.49 | 3.21
| FP8 | 36/36 | 11.4 | 2.15 | 1.39 | 0.32
| FP9 | 36/36 | 11.4 | 6.78 | 5.06 | 0.84
LeastSquares | FP1 | 36/36 | 11.4 | 3.09 | 0.62 | 0.35
| FP2 | 36/36 | 11.4 | 4.01 | 1.68 | 0.20
| FP3 | 36/36 | 13.6 | 0.75 | 1.31 | 0.54
| FP4 | 31/36 | 13.6 | 3.93 | 1.35 | 0.71
| FP5 | 36/36 | 13.6 | 3.89 | 1.76 | 0.85
| FP6 | 36/36 | 13.6 | 5.62 | 1.47 | 0.50
| FP7 | 36/36 | 14.9 | 7.82 | 1.13 | 0.48
| FP8 | 35/36 | 9.6 | 2.48 | 1.95 | 1.60
| FP9 | 36/36 | 14.9 | 6.33 | 1.25 | 0.36
MLG | FP1 | 36/36 | 14.9 | 2.94 | 0.55 | 0.42
| FP2 | 36/36 | 11.2 | 4.03 | 1.87 | 0.31
| FP3 | 36/36 | 11.2 | 0.29 | 1.49 | 0.23
| FP4 | 36/36 | 10.8 | 4.55 | 1.40 | 0.60
| FP5 | 36/36 | 11.2 | 4.20 | 1.97 | 0.67
| FP6 | 36/36 | 10.6 | 5.28 | 1.25 | 0.38
| FP7 | 36/36 | 10.5 | 7.86 | 0.96 | 0.37
| FP8 | 36/36 | 10.1 | 2.77 | 1.89 | 1.51
| FP9 | 36/36 | 10.6 | 5.76 | 0.91 | 0.56
Table 3: Performance of several 2D multilateration implementations at nine
firing positions (FP). The ratio $N/M$ is number of shots $N$ detected out of
$M$ shots fired, while $n$ is the mean number of sensors used in the location.
The survey error $\epsilon$ is the 2D RMS distance from the survey location;
$\sigma_{1}$ and $\sigma_{2}$ are the square roots of the eigenvalues of the
covariance matrix, determined using PCA. These correspond to the standard
deviations $\sigma_{X}$, $\sigma_{Y}$ of a bivariate normal distribution.
Figure 14: Results of MLG implementation under 2D constraint for Firing
Positions 1–9.

###  V.2 Effects of geometric constraint

While a ShotSpotter array in California’s Central Valley has a vertical range
of just 10 m over its entire 5 km extent, the Pittsburgh array under study has
sensor elevation change (comprising terrain plus building height) of up to 160
m over 1 km. This makes the array suitable for a study of the effect of
geometric constraint on location accuracy.

All incidents were relocated under the three geometric constraints (2D, 2.5D,
3D) to quantify the effect of including elevation data in the multilateration
calculation. The system was permitted to use any algorithm in selecting pulse
sets for multilateration; the MLG algorithm was used for the final location
computation. The results (Table 4) show that 2D location results are very
similar, with no apparent improvement due to the use of 3D straight-line
distance as a propagation model.

A possible explanation is that a shot fired at ground level makes one
diffracted path at the nearest building and then travels above the roofline
along a path well-approximated by the 2D straight-line distance. 2D
multilateration is implicitly locating in the plane of the sensors, so the
system may be locating the best-fit location of a point slightly above the
shooter. Because the 2D distances between sensors are so much larger than
height of typical sensors above ground, the error introduced by this process
is small.

These results suggest that the 2D geometric constraint is appropriate for
acoustic multilateration of gunshots even in hilly terrain. The resulting
locations are of high accuracy with precision of a few meters even when the
sensors comprising the array contain a non-plane component. When a sufficient
number of sensors are available for 3D multilateration, the results obtained
show acceptable 2D location accuracy. However, the error in the $z-$axis is
high, and the results show that the largely planar arrays used in this study
do not provide sufficient elevation determination to identify the floor of a
multi-story structure from which a shot might be fired.

Figure 15: Schematic side view of how sound from a muzzle blast discharged
near ground level relies on sound diffraction to reach a roof-mounted sensor.
From the point of view of the sensor, the source appears to be an above-
ground-level source with attenuated amplitude and nearly 2D propagation path.
Figure 16: Schematic top view of the shortest acoustic path from a shot
discharged near ground level to several rooftop-mounted sensors. Diffraction
shifts the apparent position of the source from the point of view of each
reporting sensor, but errors from multiple sensors tend to cancel when sensors
surround the actual source.

Table 4 also shows the elevation survey error $\epsilon_{z}$, which is the RMS
difference between presumed shooter elevation of $1.5\text{\,}\mathrm{m}$
above grade and the output of the multilateration algorithm. The elevation
survey error is only meaningful under the 3D constraint, as the shooter
elevation is arbitrarily set to $1.5\text{\,}\mathrm{m}$ above grade under the
2D and 2.5D geometric constraints.

| FP | $N/M$ | n | $\epsilon$ | $\sigma_{1}$ | $\sigma_{2}$ | $\epsilon_{z}$
---|---|---|---|---|---|---|---
|  |  | $(\textrm{m})$ | $(\textrm{m})$ | $(\textrm{m})$ | $(\textrm{m})$ |
2D | FP1 | 36/36 | 17.53 | 2.96 | 0.67 | 0.28 | 34.95
| FP2 | 36/36 | 14.00 | 3.60 | 2.82 | 0.38 | 0.22
| FP3 | 36/36 | 9.86 | 0.36 | 1.43 | 0.58 | 0.54
| FP4 | 36/36 | 13.72 | 5.57 | 1.26 | 0.75 | 0.35
| FP5 | 36/36 | 11.44 | 3.30 | 2.18 | 0.76 | 0.00
| FP6 | 36/36 | 13.61 | 5.87 | 1.56 | 0.76 | 0.52
| FP7 | 36/36 | 14.94 | 4.57 | 0.48 | 0.23 | 0.37
| FP8 | 36/36 | 11.17 | 3.11 | 2.41 | 1.59 | 0.66
| FP9 | 36/36 | 10.58 | 6.44 | 1.14 | 0.46 | 0.71
2.5D | FP1 | 36/36 | 17.53 | 3.02 | 0.70 | 0.23 | 34.95
| FP2 | 35/36 | 13.58 | 3.62 | 2.72 | 0.35 | 0.22
| FP3 | 36/36 | 9.83 | 0.54 | 1.21 | 0.38 | 0.54
| FP4 | 36/36 | 13.78 | 5.66 | 1.43 | 1.08 | 0.35
| FP5 | 36/36 | 11.44 | 1.96 | 1.87 | 0.79 | 0.00
| FP6 | 36/36 | 13.61 | 6.77 | 1.07 | 0.77 | 0.52
| FP7 | 36/36 | 14.94 | 4.91 | 0.85 | 0.44 | 0.37
| FP8 | 36/36 | 11.17 | 4.96 | 2.42 | 1.64 | 0.66
| FP9 | 36/36 | 10.58 | 6.15 | 1.15 | 0.32 | 0.71
3D | FP1 | 36/36 | 17.53 | 3.76 | 0.51 | 0.32 | 42.40
| FP2 | 35/36 | 13.58 | 4.46 | 2.56 | 0.28 | 14.57
| FP3 | 36/36 | 9.83 | 2.51 | 1.95 | 0.72 | 45.55
| FP4 | 36/36 | 13.78 | 5.63 | 1.35 | 1.01 | 4.79
| FP5 | 36/36 | 11.44 | 2.29 | 3.17 | 0.73 | 3.82
| FP6 | 36/36 | 13.61 | 6.43 | 0.81 | 0.62 | 20.35
| FP7 | 36/36 | 14.94 | 4.84 | 1.73 | 0.81 | 3.73
| FP8 | 36/36 | 11.17 | 4.11 | 2.92 | 1.66 | 41.69
| FP9 | 35/36 | 10.28 | 5.68 | 0.86 | 0.38 | 18.75
Table 4: Performance of MLG multilateration implementation under 2D, 3D, and
“2.5D” (3D constrain-to-plane) geometric constraints. Table columns have the
same meaning as Table 3, with the addition of $\epsilon_{z}$, the distance
from multilateration elevation to survey elevation. Note that for 2D and 2.5D
solutions, the elevation coordinate is set arbitrarily at
$1.5\text{\,}\mathrm{m}$ above the ground elevation at the latitude, longitude
of the multilateration results.

###  V.3 Effects of sensor density

A reporting sensor is one on which at least one gunshot impulse is detected
programmatically by a sensor from shots fired at a particular location. A
participating sensor is one identified by the pulse selection algorithms as a
member of a subset suitable for multilateration.

This series of firing tests was conducted on a night with still air after
autumn leaf drop, resulting in good acoustic propagation of the gunshot
sounds. The reporting sensor count was as high as 21 for the most powerful
firearm (.45 cal) fired from the elevated position at FP1; the median number
of participating sensors across all weapons and firing positions was 10, well
above the minimum $d+2=4$ participating sensors required for unambiguous 2D or
2.5D multilateration.

To investigate the effects of reduced deployment density, the multilateration
performance of reduced-density arrays was evaluated using a Monte-Carlo
simulation. These simulations make use of the real-time pulse detection
performance; the simulation is in varying the set of sensors considered to
have been deployed during the test. 25 randomly-generated arrays were
constructed at each firing point for the desired number of participating
sensors. Simulations were run for 4, 5, 6, 8, 10, 12 and 15 sensors,
corresponding to a sensor density of approximately
$2\text{\,}{sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}7.5\text{\,}{sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}$.
In some cases the target array density exceeded the number of sensors
participating in a shot, so the maximum value of participating sensors was
used. Multilateration was performed using the MLG algorithm under the 2D
geometric constraint. Simulated incidents were filtered using the same pulse
selection and filtering criteria used in production systems, with the result
that some valid multilateration solutions were discarded because the signal
amplitude on the strongest participating sensor failed to meet the (arbitrary)
73 dB SPL amplitude threshold applied to all incidents.

The results show that detection rates and location accuracy increase as the
number of participating sensors increases. With six participating sensors,
96.3% of incidents can be located to an accuracy of $15\text{\,}\mathrm{m}$ or
better. Detection rates and location accuracy continue to improve as the
number of participating sensors increases above six, but additional reporting
sensors provide only marginal improvements. The six-sensor result may be
interpreted as “sufficient to provide good array geometry after removing non-
line-of-sight paths”.

As this technique varies the deployed sensor density rather than the
propagation distance, it is best used to estimate the precision and accuracy
obtainable as a function of the number of reporting sensors, a known quantity
in a forensic context. It is less suitable for predicting the sensor density
required to ensure good detection rates under non-ideal conditions. The
presence of wind, rain, traffic noise and similar background noise sources
reduce signal-to-noise ratio. Furthermore, while the simulation shows that any
set of six reporting sensors will yield accurate multilateration results,
implicit in the simulation is that the reporting sensors are well-distributed,
since the sensors were randomly drawn from a set of participating sensors that
fully surround the shooter, as shown in Figure 13. Wind may refract the sound
away from a given direction entirely, and some barriers (such as highway
embankments or sound walls) may attenuate the sound below the point of
detectability, potentially resulting in a situation where the shooter is far
outside the convex hull of reporting sensors. Because the effect is due to the
geometry of the receiver array, the error of a multilateration solution can be
estimated using the analytical approach of (Mathias _et al._ , 2008) or using
Monte-Carlo simulation, in which the shooter location is repeatedly re-
computed after perturbing the input measurements according to the measurement
error estimates in Section III. Because of the high number of participating
sensors in the Pittsburgh live-fire test, the data published in the present
work provide a good basis for future study of the effects of array geometry on
location accuracy.

Figure 17: Location accuracy of randomly-generated, reduced-density arrays,
determined using a Monte-Carlo simulation of 25 random arrays of participating
sensors at each firing position. These plots show performance with a maximum
of 4 (left column), 8 (middle column), or 12 (right columns) sensors,
corresponding to a sensor density of
$2\text{\,}{sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}6\text{\,}{sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}$.
Firing positions shown reflect those with highest sensor participation during
the live fire test (FP1, top row), medium participation (FP6, middle row) or
low participation (FP3, bottom). See 18 for detection rates and location
accuracy at all firing positions. These results show that the robustness of
multilateration results is improved by deploying a sufficient number of
sensors to ensure good participation. Figure 18: Cumulative distribution
function (CDF) of simulated detection accuracy with random arrays containing
varying numbers of participating sensors. The CDF for each simulations is
based on relocating all shots using 25 random arrays comprising 4, 5, 6, 8, or
10 participating sensors at each firing position, corresponding to a sensor
density of approximately
$2\text{\,}{sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}5\text{\,}{sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}$.
The CDF of the live-fire test from Figure 12 is included for comparison.

##  VI Conclusion

Source location via acoustic multilateration is a well-developed technique
that was reduced to practice during the First World War. Inexpensive
computation and cellular radio networks allow the use of these techniques by
automated systems to detect gunshots in near real time.

For the sensor arrays considered, the multilateration algorithm of (Mathias
_et al._ , 2008) had the greatest precision. Use of a 3D straight-line
propagation model (either unconstrained or constrained to the ground plane)
did not significantly improve location accuracy over the simpler 2D
propagation model. Under these conditions, acoustic multilateration provided a
2D accuracy between $0.36\text{\,}\mathrm{m}6.44\text{\,}\mathrm{m}$ with
respect to surveyed location and a 2D precision (or relative accuracy) of
$0.48\text{\,}\mathrm{m}2.82\text{\,}\mathrm{m}$. High relative accuracy is
essential when performing forensic analysis of situations that may involve
multiple shooters. The precision of elevation estimates was much lower, with
full 3D solutions yielding an RMS accuracy vs digital elevation model of
$17.8\text{\,}\mathrm{m}$. The array geometry tested was not adequate to
identify the elevation of shots fired from the upper stories of a building. In
coverage areas where firearm discharge from an elevated structure is of
concern, sensors should be positioned at both the roof and base of tall
structures in the area to provide better elevation precision.

Detection on the minimum number of sensors required (nominally three in two
dimensions) can produce an accurate result when acoustic paths are near line-
of-sight, but the fully-determined solution is not robust under the presence
of timing errors. Additional information (such as waveform analysis) must be
used to estimate the likelihood that the TDoAs are good estimates of the line-
of-sight travel time. With four sensors, the presence of a timing error can be
detected (since each 3-tuple produces a different solution) but not corrected.
With five or more sensors, error-inducing multipath inputs may be identified
and rejected. Simulations with random sensor arrays show that excellent
precision and accuracy are reliably obtained with arbitrary arrays of six
reporting sensors drawn at random from a spatially-distributed set of
participating sensors. Detection rate and accuracy continue to improve as the
number of participating sensors increases, with additional sensors above six
providing only marginal improvement. The six-sensor criterion corresponds to a
deployed sensor density of approximately
$3\text{\,}{sensors}\mathrm{/}\mathrm{k}\mathrm{m}\mathrm{{}^{2}}$ in the
arrays studied; a higher sensor density is needed to ensure sufficient sensor
participation in practical deployments. Higher density arrays help ensure an
adequate participating sensor count (preferably six, but at least three) in
highly built-up areas, at the edge of coverage areas, or in adverse weather
conditions. The ideal sensor deployment configuration is one that will detect
a high fraction of outdoor firearm discharges from any part of the coverage
area under all weather conditions, and do so using a minimal number of
sensors. A model that incorporates the combined effects of structures,
terrain, foliage, wind, and ground reflection would be of value in the design
of effective, economical sensor arrays for acoustic gunshot location, and is a
topic we seek to address in future work.

###### Acknowledgements.

The authors thank the Pittsburgh Bureau of Police for conducting a
comprehensive live-fire test of the ShotSpotter system in their community. We
also thank the ShotSpotter Operations Team, and Project Manager Athnaiel Bazi
in particular, for designing, deploying, and maintaining an effective array
for gunshot detection. The authors are employed by and hold an ownership
interest in ShotSpotter, Inc.

## References

  * Aguilar (2015) Aguilar, J. R. (2015). “Gunshot detection systems in civilian law enforcement,” Journal of the Audio Engineering Society 63(4), 280–291.
  * Albert and Liu (2010) Albert, D. G., and Liu, L. (2010). “The effect of buildings on acoustic pulse propagation in an urban environment,” The Journal of the Acoustical Society of America 127(3), 1335–1346.
  * Aubin and Goldstein (2014) Aubin, D., and Goldstein, C. (2014). _The War of Guns and Mathematics: Mathematical Practices and Communities in France and Its Western Allies Around World War I_ , 42, 144–149 (American Mathematical Society).
  * Beck (2019) Beck, S. D. (2019). “Who fired when: Associating multiple audio events from uncalibrated receivers,” in _Audio Engineering Society Conference: 2019 AES International Conference on Audio Forensics_.
  * Beck _et al._ (2011) Beck, S. D., Nakasone, H., and Marr, K. W. (2011). “Variations in recorded acoustic gunshot waveforms generated by small firearms,” The Journal of the Acoustical Society of America 129(4), 1748–1759.
  * Begault _et al._ (2019) Begault, D. R., Beck, S. D., and Maher, R. C. (2019). “Overview of forensic audio gunshot analysis techniques,” in _Audio Engineering Society Conference: 2019 AES International Conference on Audio Forensics_.
  * Boashash (2015) Boashash, B. (2015). _Time-frequency signal analysis and processing: a comprehensive reference_ (Academic Press).
  * Brustad and Freytag (2005) Brustad, B. M., and Freytag, J. C. (2005). “A survey of audio forensic gunshot investigations,” in _Audio Engineering Society Conference: 26th International Conference: Audio Forensics in the Digital Age_.
  * Calhoun (2011) Calhoun, R. B. (2011). “Systems and methods of processing impulses including bullet pulses and/or muzzle pulses in association with time domain representations” US Patent 7,961,550.
  * Calhoun (2013) Calhoun, R. B. (2013). “Systems and methods with improved three-dimensional source location processing including constraint of location solutions to a two-dimensional plane” U.S. patent 8,369,184 (Feb 5,2013).
  * Don and Cramond (1987) Don, C. G., and Cramond, A. J. (1987). “Impulse propagation in a neutral atmosphere,” The Journal of the Acoustical Society of America 81(5), 1341–1349.
  * Embleton (1996) Embleton, T. F. W. (1996). “Tutorial on sound propagation outdoors,” The Journal of the Acoustical Society of America 100(1), 31–48.
  * Finan and Hurley (1997) Finan, J. S., and Hurley, W. J. (1997). “McNaughton and Canadian operational research at Vimy,” Journal of the Operational Research Society 48(1), 10–14.
  * Franke _et al._ (1988) Franke, S. J., Raspet, R., and Liu, C. H. (1988). “Numerical predictions of atmospheric sound-pressure levels in shadow zones,” The Journal of the Acoustical Society of America 83(2), 816–820.
  * Friedlander (1946) Friedlander, F. G. (1946). “The diffraction of sound pulses i. diffraction by a semi-infinite plane,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 186(1006), 322–344.
  * He _et al._ (2016) He, K., Zhang, X., Ren, S., and Sun, J. (2016). “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 770–778.
  * Hirsch and Bertels (2013) Hirsch, K.-W., and Bertels, W. (2013). “Estimation of the directivity pattern of muzzle blasts,” Proc. AIA-DAGA 961–963.
  * Hornikx (2016) Hornikx, M. (2016). “Ten questions concerning computational urban acoustics,” Building and Environment 106, 409–421.
  * Hosier and Donavan (1979) Hosier, R. N., and Donavan, P. R. (1979). “Microphone windscreen performance,” Final Report National Bureau of Standards, Washington, DC. Acoustical Engineering Div. .
  * Kapralos _et al._ (2008) Kapralos, B., Jenkin, M., and Milios, E. (2008). “Sonel mapping: a probabilistic acoustical modeling method,” Building Acoustics 15(4), 289–313.
  * Lawrence _et al._ (2018) Lawrence, D. S., La Vigne, N. G., Goff, M., and Thompson, P. S. (2018). “Lessons learned implementing gunshot detection technology: Results of a process evaluation in three major cities,” Justice Evaluation Journal 1(2), 109–129.
  * Le Bot and Bocquillet (2000) Le Bot, A., and Bocquillet, A. (2000). “Comparison of an integral equation on energy and the ray-tracing technique in room acoustics,” The Journal of the Acoustical Society of America 108(4), 1732–1740.
  * Lin _et al._ (2014) Lin, I.-C., Hsieh, Y.-R., Shieh, P.-F., Chuang, H.-C., and Chou, L.-C. (2014). “The effect of wind on low frequency noise,” in _43rd International Congress and Exposition on Noise Control Engineering (Internoise 2014)_ , Institute of Noise Control Engineering, Vol. 249, pp. 1137–1148.
  * Maher (2009) Maher, R. C. (2009). “Audio forensic examination,” IEEE Signal Processing Magazine 26(2), 84–94.
  * Maher (2011) Maher, R. C. (2011). “Acoustical modeling of gunshots including directional information and reflections,” in _Audio Engineering Society Convention 131_ , Audio Engineering Society.
  * Maher and Hoerr (2018) Maher, R. C., and Hoerr, E. (2018). “Audio forensic gunshot analysis and multilateration,” in _Audio Engineering Society Convention 145_ , Audio Engineering Society.
  * Maher and Routh (2015) Maher, R. C., and Routh, T. (2015). “Advancing forensic analysis of gunshot acoustics,” in _Audio Engineering Society Convention 139_ , Audio Engineering Society.
  * Maher and Shaw (2010) Maher, R. C., and Shaw, S. R. (2010). “Directional aspects of forensic gunshot recordings,” in _Audio Engineering Society Conference: 39th International Conference: Audio Forensics: Practices and Challenges_ , Audio Engineering Society.
  * Mathias _et al._ (2008) Mathias, A., Leonardi, M., and Galati, G. (2008). “An efficient multilateration algorithm,” in _Digital Communications-Enhanced Surveillance of Aircraft and Vehicles, 2008. TIWDC/ESAV 2008. Tyrrhenian International Workshop on_ , IEEE, pp. 1–6.
  * Mehra _et al._ (2014) Mehra, R., Raghuvanshi, N., Chandak, A., Albert, D. G., Keith Wilson, D., and Manocha, D. (2014). “Acoustic pulse propagation in an urban environment using a three-dimensional numerical simulation,” The Journal of the Acoustical Society of America 135(6), 3231–3242.
  * National Law Enforcement and Corrections Technology Center (2007) National Law Enforcement and Corrections Technology Center (2007). “Spot the Shot,” TechBeat 1–2.
  * Parker _et al._ (2007) Parker, M., Ketcham, S., and Cudney, H. (2007). “Acoustic wave propagation in urban environments,” in _2007 DoD High Performance Computing Modernization Program Users Group Conference_ , IEEE, pp. 233–237.
  * Pasareanu _et al._ (2018) Pasareanu, S. M., Burdisso, R. A., and Remillieux, M. C. (2018). “A numerical hybrid model for outdoor sound propagation in complex urban environments,” The Journal of the Acoustical Society of America 143(3), EL218–EL224.
  * Pauli _et al._ (2004) Pauli, M., Seisler, W., Price, J., Williams, A., Maraviglia, C., Evans, R., Moroz, S., Ertem, M. C., Heidhausen, E., and Burchick, D. A. (2004). “Infrared detection and geolocation of gunfire and ordnance events from ground and air platforms,” Technical Report ADA460225 .
  * Raspet _et al._ (2006) Raspet, R., Webster, J., and Dillion, K. (2006). “Framework for wind noise studies,” The Journal of the Acoustical Society of America 119(2), 834–843.
  * Reddi (1993) Reddi, S. S. (1993). “An exact solution to range computation with time delay information for arbitrary array geometries,” IEEE Transactions on Signal Processing 41(1), 485–486.
  * Remillieux _et al._ (2012) Remillieux, M. C., Corcoran, J. M., Haac, T. R., Burdisso, R. A., and Svensson, U. P. (2012). “Experimental and numerical study on the propagation of impulsive sound around buildings,” Applied Acoustics 73(10), 1029–1044.
  * Routh and Maher (2016) Routh, T. K., and Maher, R. C. (2016). “Recording anechoic gunshot waveforms of several firearms at 500 kilohertz sampling rate,” in _Proceedings of Meetings on Acoustics 171 ASA_ , ASA, Vol. 26, p. 030001.
  * Schomer _et al._ (1990) Schomer, P. D., Raspet, R., Wagner, M., Walker, D., and Marshall, D. (1990). “Methods for detecting low-frequency signals in the presence of strong winds,” Technical Report ADA223980 .
  * Showen (1997) Showen, R. (1997). “Operational gunshot location system,” in _Surveillance and Assessment Technologies for Law Enforcement_ , International Society for Optics and Photonics, Vol. 2935, pp. 130–140.
  * Showen _et al._ (2008) Showen, R. L., Calhoun, R. B., Chu, W. C., and Dunham, J. W. (2008). “Acoustic gunshot location in complex environments—concepts and results,” in _Sensors, and Command, Control, Communications, and Intelligence (C3I) Technologies for Homeland Security and Homeland Defense VII_ , S.P.I.E., Vol. 6943, pp. 694318–1–694318–11.
  * Stevens and Murphy (2014) Stevens, F., and Murphy, D. (2014). “Spatial impulse response measurement in an urban environment,” in _Audio Engineering Society Conference: 55th International Conference: Spatial Audio_ , Audio Engineering Society.
  * Van der Kloot (2005) Van der Kloot, W. (2005). “Lawrence Bragg’s role in the development of sound-ranging in World War I,” Notes and Records of the Royal Society 59(3), 273–284.
  * Watkins _et al._ (2002) Watkins, C., Green Mazerolle, L., Rogan, D., and Frank, J. (2002). “Technological approaches to controlling random gunfire: Results of a gunshot detection system field test,” Policing: An International Journal of Police Strategies & Management 25(2), 345–370.
  * Weber (1939) Weber, W. (1939). “Das schallspektrum von knallfunken und knallpistolen mit einem beitrag über anwendungsmöglichkeiten in der elektroakustischen meßtechnik,” Ph.D. thesis, University of Berlin.
  * Wiener and Keast (1959) Wiener, F. M., and Keast, D. N. (1959). “Experimental study of the propagation of sound over ground,” The Journal of the Acoustical Society of America 31(6), 724–733.
  * Wilson (2003) Wilson, D. K. (2003). “The sound-speed gradient and refraction in the near-ground atmosphere,” The Journal of the Acoustical Society of America 113(2), 750–757.
  * Wuttke (1992) Wuttke, J. (1992). “Microphones and wind,” Journal of the Audio Engineering Society 40(10), 809–817.


# Demonstration of nanosized high-impedance superconducting NbRe Dayem bridges

S. Battisti<EMAIL_ADDRESS>NEST, Istituto Nanoscienze-CNR and
Scuola Normale Superiore, Pisa, Italy  J. Koch  Department of Physics,
University of Konstanz, 78457 Konstanz, Germany  A. Paghi  NEST, Istituto
Nanoscienze-CNR and Scuola Normale Superiore, Pisa, Italy  L. Ruf  Department
of Physics, University of Konstanz, 78457 Konstanz, Germany  A. Gulian
Advanced Physics Laboratory, Institute for Quantum Studies, Chapman
University, Burtonsville, MD 20866, USA  S. Teknowijoyo  Advanced Physics
Laboratory, Institute for Quantum Studies, Chapman University, Burtonsville,
MD 20866, USA  C. Cirillo  CNR-Spin, c/o Università degli Studi di Salerno,
I-84084 Fisciano (Sa), Italy  Z. Makhdoumi Kakhaki  Dipartimento di Fisica “E.
R. Caianiello”, Università degli Studi di Salerno, via Giovanni Paolo II 132,
I-84084, Fisciano, Salerno, Italy  C. Attanasio  Dipartimento di Fisica “E. R.
Caianiello”, Università degli Studi di Salerno, via Giovanni Paolo II 132,
I-84084, Fisciano, Salerno, Italy  E. Scheer  Department of Physics,
University of Konstanz, 78457 Konstanz, Germany  A. Di Bernardo  Department of
Physics, University of Konstanz, 78457 Konstanz, Germany  2\. Dipartimento di
Fisica “E. R. Caianiello”, Università degli Studi di Salerno, via Giovanni
Paolo II 132, I-84084, Fisciano, Salerno, Italy  G. De Simoni  NEST, Istituto
Nanoscienze-CNR and Scuola Normale Superiore, Pisa, Italy  F. Giazotto
<EMAIL_ADDRESS>NEST, Istituto Nanoscienze-CNR and Scuola Normale
Superiore, Pisa, Italy

###### Abstract

Here we demonstrate nano-sized superconducting Dayem-bridge weak-links made of
different stoichiometric compositions of NbRe. Our devices possess a
relatively high critical temperature, normal-state resistance, and kinetic
inductance. In particular, the high kinetic inductance makes this material a
good alternative to more conventional niobium-based superconductors (e.g., NbN
or NbTiN) for the realization of superinductors and high-quality factor
resonators, whereas the high normal-state resistance yields a large output
voltage in superconducting switches and logic elements realized upon this
compound. Moreover, out-of-plane critical magnetic fields exceeding 2 T ensure
that possible applications requiring high magnetic fields can also be
envisaged. Altogether, these features make this material appealing for a
number of applications in the framework of quantum technologies.

##  1  Introduction

NbRe is a relatively-high critical temperature ($T_{c}^{bulk}\sim 9$K 1) non-
centrosymmetric superconductor 2, which was proposed as a suitable material
for superconducting single-photon detectors (SSPD) 3, 4, 5, 6. Thanks to a
coherence length as short as $\xi\sim 5$nm, this material might provide a
platform for the realization of ultra-thin superconducting films (thickness
$t_{NbRe}\simeq 3.5$nm 7) with large critical temperatures ($T_{c}\simeq
5.3$K) easily accessible in 4He closed-cycle refrigerators. Such a short
coherence length is the consequence of a sizable normal-state resistivity
($\rho_{N}\sim 3.1x10^{-5}$ $\Omega$m deduced by the normal state resistance
and the dimensions of our structures), which turns out in a high $I_{S}R_{N}$
product, where $I_{S}$ is the switching current and $R_{N}$ is the normal-
state resistance of the film. This peculiar characteristic, although somewhat
rare for metals, can be achieved in dirty superconductors such as NbN.
Differently, in polycrystalline NbRe this feature stems from the highly
disordered nature of the film 7. Another crucial figure of merit of NbRe is
its high kinetic inductance, which makes the devices based on this material
appealing for all those applications requiring a high inductance such as
superinductors 8, 9, 10, fluxonium qubits 11, high Q-factor resonators 12, and
radiation sensors 13, 14. Here, we report for the first time on devices based
on a Dayem-bridge (DB) geometry made of NbRe films with different
stoichiometry. The DB structure can serve as a tool to realize several
monolithic devices such as superconducting quantum interference devices
(SQUIDs), magnetometers 15, and logic gates 16. Furthermore, the opportunity
to suppress the supercurrent of the DB via an external gate voltage 17 makes
this kind of structure appealing for the realization of variable kinetic
inductors with high dynamic range or nanosized superconducting switches with
high output voltage due to the high $I_{S}R_{N}$ product.

Figure 1:  (a): False color scanning electron micrograph (SEM) of one
representative Dayem-bridge device with composition A. The 4-wire measurement
setup and the characteristic nanobridge lateral dimensions are also indicated.
$I$ indicates the current flowing through the device, while $V$ the voltage
drop occurring across it. (b): False color SEM image of a representative
device with composition B (the same 4-wire setup was used for this device).
(c): Resistance versus temperature curve ($R$ vs $T$) of the device shown in
panel (a). (d): $R$ vs $T$ of the device shown in panel (b). The dashed lines
indicate the critical temperatures $T_{c}^{exp}$ of the shown devices.

##  2  Film deposition, device fabrication, and measurement setup

In this work, we have investigated two different types of composition of NbRe
thin films, i.e., composition A and B, which have a stoichiometric ratio of
Nb0.1Re0.9 and Nb0.18Re0.82, respectively.
We start by describing the fabrication process for the samples with
composition A. The Re-Nb films were prepared via magnetron sputtering using an
ATC series UHV Hybrid deposition system (AJA International, Inc.) with a base
pressure of $1\times 10^{-8}$ Torr. The Re target (ACI Alloys, Inc., $99.99\%$
purity) and the Nb target (Testbourne, Ltd., $99.95\%$ purity) were both
accommodated inside $1.5"$ diameter DC guns. The sapphire substrate (AdValue
Technology, thickness $650\mu$m, C-cut) was cleaned thoroughly with isopropyl
alcohol before it was mounted on the holder. In the chamber’s configuration,
the substrate holder is at the center of the chamber facing upwards, while the
sputtering guns are located at the top. The substrate is rotated in plane
throughout the whole deposition process to ensure a homogeneous deposition
layer over the whole surface. A pre-deposition in-situ cleaning of the
substrate involved heating it up to $700$ ∘C for 10 minutes followed by a
gentle bombardment of Ar+ ions at $600$ ∘C for 5 minutes using a Kauffman
source at $0.75$ mTorr. Plasma striking for guns was performed at $30$ mTorr
with shutters closed and with the initial power $30$ W for each target. Then
the gun powers were ramped up to $150$ W for Nb and $500$ W for Re, with an
increase of the chamber pressure to $5$ mTorr in Ar flow. Afterward, the
deposition of Nb and Re took place. For better adhesion, the Re shutter was
opened for $10$ seconds later than the Nb shutter. Simultaneous Nb and Re
sputtering lasted $150$ s, after that the Re shutter was closed, and in $10$
seconds the Nb shutter was also closed. After this, the temperature was again
raised to $900$ ∘C for in-situ annealing during $30$ minutes and then cooled
down to ambient temperature. All the heating/cooling protocols consistently
used a $30$ ∘C/min ramp rate. To pattern the DB structures via reactive ion
etching (RIE), a positive electron beam lithography (EBL) was performed with a
PMMA mask on the sputtered films. Then $60$ nm of aluminum (Al) was deposited
with an e-beam evaporator in order to have an Al hard mask, which is resistant
to fluorine etching, on top of our NbRe film. The films were etched using a
CF4, Ar and O2 gas mixture of 40/6/4 sccm at a power of $150$ W . Finally, the
removal of the Al hard mask in an alkaline solution provides the finished
device shown in Fig.1(a).
The composition B films were sputtered via a stochiometric NbRe (Nb0.18Re0.82)
target ((Testbourne Ltd., 99.95% purity)) with a power of 250 W and an Ar
pressure of 3 mTorr at a base pressure of $1\times 10^{-8}$ mTorr, which
resulted in a deposition rate of $0.33$ nm/s. The 20-nm-thick NbRe films were
sputtered on Al2O3 substrates. To pattern the DBs via RIE, a negative
lithography was performed with an etching-resistant negative resist. The
etching was performed with an inductively-coupled plasma – reactive ion etcher
(ICP-RIE), where the samples were etched for 30 s (at first 15 s and than
three steps with 5 s after each checking the resistance in the wafer prober)
in $26$ sccm Ar + $26$ sccm Cl2 at $10$ mTorr with an RIE power of $20$ W and
an ICP power of $750$ W (the sample plate was cooled to 10°C). A
representative device with composition B, produced with this process is shown
in Fig.1(b).

##  3  Experimental results

Figure 2:  (a)[(b)]: Current vs voltage ($I-V$) characteristics for selected values of the bath temperature for a representative DB device with composition A[(B)]. The gray areas correspond to the dissipationless region. The curves are horizontally shifted for clarity. Figure 3:  (a)[(b)]: Switching ($I_{S}$) and retrapping ($I_{R}$) currents as a function of bath temperature $T$ for a representative device with composition A[B]. The error bars are smaller than the size of the circles. The red line shows a fit with the Bardeen model in Eq.1, whereas the green line shows the KO1 model in Eq. 2a calculated with the measured parameters. Dashed lines are guides to the eyes. Table 1: The table summarises all the parameters of the measured devices. $T_{c}^{exp}$ indicates the measured critical temperature, while $T_{c_{theo}}$ and $I_{S_{theo}}^{0}$ indicate respectively the estimated critical temperature and switching current obtained with a Bardeen fit as in Eq. 1. The product $I_{S}R_{N}$ is calculated for the maximum experimental switching current available, i.e., for the lower temperatures in Fig.2(a-b) |  composition A (Re0.9Nb0.1) |  composition B (Re0.82Nb0.18)
---|---|---
|  Device 1A |  Device 2A |  Device 1B |  Device 2B
L [nm] |  180 |  170 |  200 |  200
W [nm] |  80 |  75 |  64 |  56
t [nm] |  30 |  30 |  20 |  20
$T_{c}^{exp}$ [K] |  $7.4\pm 0.1$ |  $7.4\pm 0.1$ |  $6.1\pm 0.1$ |  $6.1\pm 0.1$
$T_{c_{theo}}$ [K] |  $7.1\pm 0.1$ |  $7.0\pm 0.1$ |  $5.82\pm 0.03$ |  $5.91\pm 0.06$
$I_{S_{theo}}^{0}$ [$\mu$A] |  $158\pm 7$ |  $136\pm 2$ |  $33\pm 3$ |  $18\pm 3$
$R_{N}$ [$\Omega$] |  $2276\pm 49$ |  $2301\pm 38$ |  $1136\pm 5$ |  $1171\pm 6$
$I_{S}R_{N}$ [mV] |  $250\pm 5$ |  $232\pm 4$ |  $24.4\pm 0.4$ |  $14.4\pm 0.3$

The devices were then wired, as displayed in Fig1(a), with a 4-wire standard
setup in order to take the full characterization of their properties. The
samples were cooled down in 4He cryogen-free cryostats with a base temperature
of $\ 3$K. In Fig.1(c) the characteristic $R$ vs $T$ shows that the critical
temperature of the devices with composition A is $T_{c_{1A}}^{exp}\sim 7.4$K,
while for composition B is $T_{c_{1B}}^{exp}\sim 6.1$K [see Fig.1(d)]. The
$T_{c}^{exp}$ value for both samples is taken when $dR/dT$ is maximum [dashed
lines in Figs.1(c,d)].
Figure 2(a) shows the current-voltage (I-V) curves of device 1A for different
temperatures horizontally shifted for clarity. It can be seen that at
$T_{bath}=3$K the switching current is $I_{S}=107\mu$A and then rising the
temperature it lowers monotonically, becoming vanishing small at $T\sim 7.4$K.
Figure 2(a) marks also the hysteretical behavior of these DBs, stemming from
the Joule heating occurring in the nanobridge when in the normal state. Figure
2(b) reports the same characterization for the device 1B, which displays lower
switching (and retrapping) currents, and critical temperature. In fact, we
have that for $T_{bath}=3$K the switching current is $I_{S}=11\mu$A. This much
lower value with respect to device 1A one is due to the differences in the
geometry of the two structures (see Tab.LABEL:tab): device 1B is narrower and
thinner, resulting in both a weaker superconductivity and in a smaller section
of the superconductor. Indeed, the critical temperature of device 1B is
$T_{c_{1B}}\simeq 6.1$K. This leads to a reduction of the superconducting
energy gap $\Delta$ and of the switching current 18. An analogous discussion
can be done for the other devices investigated in this study (see
Tab.LABEL:tab). From the curves shown in Fig.1(c,d) the device’s normal-state
resistance $R_{N}$ can be extracted. This parameter is particularly relevant
for applications of these structures since it determines the output voltage
($I_{S}R_{N}$) of the device when it switches from the superconducting to the
normal state. The measured values are collected in Tab.LABEL:tab, and one can
observe that $R_{N}^{A}\simeq 2R_{N}^{B}$. Taking into account the
nanobridge’s geometrical parameters, we argue that this behavior might
originate from the higher percentage of Re in the composition of A devices.
For this reason, we believe that concerning applications that require a high
$I_{S}R_{N}$ product such as the recently investigated superconducting gate-
tunable transistors 17, 19, 20 or other superconducting switching devices such
as cryotrons 21 or nanocrytrons 22, the composition A promises better
performances. Figure 3 shows the trend of the switching current $I_{S}$ and
the retrapping current $I_{R}$ as a function of bath temperature $T$ for both
samples 1A and 1B. The curves were fitted with the Bardeen equation 23 that
holds for a thin superconducting wire:

$I_{S}(T)=I_{S_{theo}}^{0}\left[1-\left(\frac{T}{T_{c_{theo}}}\right)^{2}\right]^{3/2}.$
(1)

Furthermore, we verified the agreement with the Kulik-Omelyanchouk model in
the diffusive regime (dirty limit) KO1 24,

$I_{S}(T)=\text{max}_{\phi}(I(\phi)),$ (2a) $I_{S}(\phi)=\frac{4\pi
k_{B}T}{eR_{N}}\sum_{\omega>0}\frac{\Delta(T)\text{cos}(\phi/2)}{\delta}\text{arctan}[\frac{\Delta(T)\text{sin}(\phi/2)}{\delta}],$
(2b)

calculated with the measured parameters for $R_{N}$ and $T_{c}$. In Eq.2b,
$\delta=\sqrt{\Delta(T)^{2}\text{cos}^{2}(\phi/2)+\omega^{2}}$, $\omega=\pi
k_{B}T(2n+1)$ are the Matsubara frequencies, $k_{B}$ is the Boltzmann
constant, $\phi$ is the phase difference across the weak-link, and
$\Delta(T)=1.764k_{B}T_{c}^{exp}\textrm{tanh}[1.74\sqrt{T_{c}^{exp}/T-1}]$ is
the formula used for the temperature dependence of the superconductive gap. It
is worth noting that it was not possible to obtain a reasonable agreement
between the KO1 model and the experimental data, thereby confirming that the
length of the wire is much greater than the superconducting coherence length,
i.e., $L\gg\xi$. We note, indeed, that the KO1 model holds for $L\ll\xi$ 24.
On the other hand, the Bardeen formula fits well the experimental data over
the full temperature range, yielding a coefficient of determination ${\cal
R}^{2}>0.99$. This highlights the BCS nature of the NbRe thin films.

Figure 4(a) shows several I-V curves measured at 2K for increasing values of
the external perpendicular-to-plane magnetic field for a device with
composition B. The full trend of the switching and retrapping currents as a
function of the magnetic field is highlighed in Fig.4(b). Up to $2$T (i.e.,
the maximum magnetic field that our magnet could support), the NbRe DB still
shows a sizable switching current of the order of $\sim 6\mu$A. This suggests
that NbRe has a perpendicular-to-plane critical magnetic field at least
comparable to other high-kinetic inductance superconductors 25, 26. We
emphasize that the robustness of this material to magnetic fields makes it an
optimal candidate for experimental situations and applications where high
magnetic fluxes cannot be avoided or are even required such as e. g., quantum
Hall experiment in hybrid platforms 27.

Figure 4(c) shows the kinetic inductance per square ($L_{K}^{\square}$)
calculated within the BCS theory for a superconducting stripe (indeed in our
devices we have $W>10\xi$) 28

$L_{K}^{\square}(T)=\frac{R_{N}\hbar}{\pi\Delta(T)}\frac{1}{\textrm{tanh}(\Delta(T)/k_{B}T)},$
(3)

using for $R_{N}$ and $T_{c}$ the experimentally measured values reported in
Tab.LABEL:tab. We stress that, in principle, $L_{K}^{\square}$ obtains values
as large as several hundreds of pH/$\square$ at a few Kelvins, which is
comparable or larger than that obtained in NbN $\simeq 82$pH/$\square$ 29 and
NbTiN $\simeq 115$pH/$\square$ 30 thin films for $T>2$K. For $T\gtrsim 4$K,
the device with composition B has a higher $L_{K}^{\square}$ with respect to
that with composition A. This might be due to the lower $T_{c}$ and
$\Delta(T)$ of B-devices. Indeed, although $L_{K}^{\square}$ is proportional
to $R_{N}$, Eq. 3 reminds us that the kinetic inductance is also exponentially
dependent on $\Delta(T)$. However, it is worth highlighting that by reducing
the DB cross-section, in order to intentionally lower the critical temperature
$T_{c}$ and the gap $\Delta(T)$ of devices with composition A, it would be
possible to take advantage of the larger resistivity of the films with higher
Re concentration to achieve improved performances in terms of kinetic
inductance.
In general, we can conclude that by varying the Re percentage of these films
one can explore the parameters space and find the better performances in terms
of $L_{K}^{\square}$ and $I_{S}R_{N}$ for a given geometry.

Figure 4:  (a) $I-V$ characteristics as a function of perpendicular-to-plane
magnetic field $B$ applied to a representative device with composition B. (b):
Switching and retrapping currents of the curves in (a) as a function of $B$.
Dashed lines are guides to the eye. (c): Kinetic inductance per square
$L_{K}^{\square}$ as a function of temperature $T$ for the two compositions.
Dashed lines indicate the asymptotes of $L_{K}^{\square}$ due to $\Delta(T)$
approaching $0$ at $T=T_{c}$. All the measurements shown in panels (a) and (b)
were taken at $T_{bath}2K$.

##  4  Conclusions

In summary, we have reported the fabrication and characterization of NbRe
superconducting Dayem nanobridges. The better performance, in terms of
$T_{c}$, $L_{K}$, and $I_{S}R_{N}$, of the devices with composition A suggests
that a higher percentage of Re is beneficial for applications requiring these
characteristics. The achievable values of the kinetic inductance, which are
higher than the ones reported in literature for NbN 29, suggest that the NbRe
can be used in several applications requiring high $L_{K}$ such as
superinductors 8, 9, 10, radiation sensors 3, 4, 5, and high-quality factor
resonators 12. Moreover, the sizable critical magnetic field highlights its
possible exploitation in qubit systems 11. Given all these features, we
believe that this superconducting metal might be considered as a promising
alternative to more conventional Nb compounds (e.g., NbN and NbTiN) for
quantum technology applications 31.

We acknowledge the EU’s Horizon 2020 Research and Innovation Framework
Programme under Grant No. 964398 (SUPERGATE), No. 101057977 (SPECTRUM), and
the PNRR MUR project PE0000023- NQSTI for partial financial support. We also
acknowledge partial support of this research by US ONR grants No.
N00014-21-1-2879, No. N00014-20-1-2442 and No. N00014-23-1-2866.

## References

  * Cirillo et al. 2015 Cirillo, C.; Fittipaldi, R.; Smidman, M.; Carapella, G.; Attanasio, C.; Vecchione, A.; Singh, R. P.; Lees, M. R.; Balakrishnan, G.; Cuoco, M. Evidence of Double-Gap Superconductivity in Noncentrosymmetric Nb 0.18 Re 0.82 Single Crystals. _Phys. Rev. B_ 2015, _91_ , 134508
  * Shang et al. 2018 Shang, T. et al. Time-Reversal Symmetry Breaking in Re-Based Superconductors. _Phys. Rev. Lett._ 2018, _121_ , 257002
  * Caputo et al. 2017 Caputo, M.; Cirillo, C.; Attanasio, C. NbRe as Candidate Material for Fast Single Photon Detection. _App. Phys. Lett._ 2017, _111_ , 192601
  * Cirillo et al. 2020 Cirillo, C.; Chang, J.; Caputo, M.; Los, J. W. N.; Dorenbos, S.; Esmaeil Zadeh, I.; Attanasio, C. Superconducting Nanowire Single Photon Detectors Based on Disordered NbRe Films. _App. Phys. Lett._ 2020, _117_ , 172602
  * Cirillo et al. 2021 Cirillo, C.; Granata, V.; Spuri, A.; Di Bernardo, A.; Attanasio, C. NbReN: A Disordered Superconductor in Thin Film Form for Potential Application as Superconducting Nanowire Single Photon Detector. _Phys. Rev. Mat._ 2021, _5_ , 085004
  * Ejrnaes et al. 2022 Ejrnaes, M.; Cirillo, C.; Salvoni, D.; Chianese, F.; Bruscino, C.; Ercolano, P.; Cassinese, A.; Attanasio, C.; Pepe, G. P.; Parlato, L. Single Photon Detection in NbRe Superconducting Microstrips. _App. Phys. Lett._ 2022, _121_ , 262601
  * Cirillo et al. 2016 Cirillo, C.; Carapella, G.; Salvato, M.; Arpaia, R.; Caputo, M.; Attanasio, C. Superconducting Properties of Noncentrosymmetric Nb 0.18 Re 0.82 Thin Films Probed by Transport and Tunneling Experiments. _Phys. Rev. B_ 2016, _94_ , 104512
  * Bell et al. 2012 Bell, M. T.; Sadovskyy, I. A.; Ioffe, L. B.; Kitaev, A. Yu.; Gershenson, M. E. Quantum Superinductor with Tunable Nonlinearity. _Phys. Rev. Lett._ 2012, _109_ , 137003
  * Pop et al. 2014 Pop, I. M.; Geerlings, K.; Catelani, G.; Schoelkopf, R. J.; Glazman, L. I.; Devoret, M. H. Coherent Suppression of Electromagnetic Dissipation Due to Superconducting Quasiparticles. _Nature_ 2014, _508_ , 369–372
  * Manucharyan et al. 2009 Manucharyan, V. E.; Koch, J.; Glazman, L. I.; Devoret, M. H. Fluxonium: Single Cooper-Pair Circuit Free of Charge Offsets. _Science_ 2009, _326_ , 113–116
  * Hazard et al. 2019 Hazard, T. M.; Gyenis, A.; Di Paolo, A.; Asfaw, A. T.; Lyon, S. A.; Blais, A.; Houck, A. A. Nanowire Superinductance Fluxonium Qubit. _Phys. Rev. Lett._ 2019, _122_ , 010504
  * Doyle et al. 2008 Doyle, S.; Mauskopf, P.; Naylon, J.; Porch, A.; Duncombe, C. Lumped Element Kinetic Inductance Detectors. _J Low. Temp. Phys._ 2008, _151_ , 530–536
  * Kerman et al. 2007 Kerman, A. J.; Dauler, E. A.; Yang, J. K. W.; Rosfjord, K. M.; Anant, V.; Berggren, K. K.; Gol’tsman, G. N.; Voronov, B. M. Constriction-Limited Detection Efficiency of Superconducting Nanowire Single-Photon Detectors. _App. Phys. Lett._ 2007, _90_ , 101110
  * Santavicca et al. 2016 Santavicca, D. F.; Adams, J. K.; Grant, L. E.; McCaughan, A. N.; Berggren, K. K. Microwave Dynamics of High Aspect Ratio Superconducting Nanowires Studied Using Self-Resonance. _J. of App. Phys._ 2016, _119_ , 234302
  * Paolucci et al. 2019 Paolucci, F.; Vischi, F.; De Simoni, G.; Guarcello, C.; Solinas, P.; Giazotto, F. Field-Effect Controllable Metallic Josephson Interferometer. _Nano Lett._ 2019, _19_ , 6263–6269
  * De Simoni et al. 2020 De Simoni, G.; Puglia, C.; Giazotto, F. Niobium Dayem Nano-Bridge Josephson Gate-Controlled Transistors. _App. Phys. Lett._ 2020, _116_ , 242601
  * De Simoni et al. 2018 De Simoni, G.; Paolucci, F.; Solinas, P.; Strambini, E.; Giazotto, F. Metallic Supercurrent Field-Effect Transistor. _Nat. Nano._ 2018, _13_ , 802–805
  * Tinkham 2004 Tinkham, M. _Introduction to Superconductivity_ ; Courier Corporation, 2004
  * Ruf et al. 2023 Ruf, L.; Elalaily, T.; Puglia, C.; Ivanov, Yu. P.; Joint, F.; Berke, M.; Iorio, A.; Makk, P.; De Simoni, G.; Gasparinetti, S.; Divitini, G.; Csonka, S.; Giazotto, F.; Scheer, E.; Di Bernardo, A. Effects of Fabrication Routes and Material Parameters on the Control of Superconducting Currents by Gate Voltage. _APL Mat._ 2023, _11_ , 091113
  * Paolucci et al. 2019 Paolucci, F.; De Simoni, G.; Solinas, P.; Strambini, E.; Puglia, C.; Ligato, N.; Giazotto, F. Field-Effect Control of Metallic Superconducting Systems. _AVS Quant. Sci._ 2019, _1_ , 016501
  * Buck 1956 Buck, D. The Cryotron-A Superconductive Computer Component. _Proc. IRE_ 1956, _44_ , 482–493
  * McCaughan and Berggren 2014 McCaughan, A. N.; Berggren, K. K. A Superconducting-Nanowire Three-Terminal Electrothermal Device. _Nano Lett._ 2014, _14_ , 5748–5753
  * Bardeen 1962 Bardeen, J. Critical Fields and Currents in Superconductors. _Rev. Mod. Phys._ 1962, _34_ , 667–681
  * Golubov et al. 2004 Golubov, A. A.; Kupriyanov, M. Yu.; Il’ichev, E. The Current-Phase Relation in Josephson Junctions. _Rev. Mod. Phys._ 2004, _76_ , 411–469
  * Suzuki et al. 1987 Suzuki, M.; Baba, M.; Anayama, T. Critical Magnetic Fields of Superconducting NbN Films Prepared by Reactive Sputtering. _Jpn. J. Appl. Phys._ 1987, _26_ , 947
  * Makhdoumi Kakhaki et al. 2023 Makhdoumi Kakhaki, Z.; Leo, A.; Chianese, F.; Parlato, L.; Pepe, G. P.; Nigro, A.; Cirillo, C.; Attanasio, C. Upper Critical Magnetic Field in NbRe and NbReN Micrometric Strips. _Beilstein J. Nanotechnol._ 2023, _14_ , 45–51
  * Rickhaus et al. 2012 Rickhaus, P.; Weiss, M.; Marot, L.; Schönenberger, C. Quantum Hall Effect in Graphene with Superconducting Electrodes. _Nano Lett._ 2012, _12_ , 1942–1945
  * Annunziata et al. 2010 Annunziata, A. J.; Santavicca, D. F.; Frunzio, L.; Catelani, G.; Rooks, M. J.; Frydman, A.; Prober, D. E. Tunable Superconducting Nanoinductors. _Nanotechnology_ 2010, _21_ , 445202
  * Niepce et al. 2019 Niepce, D.; Burnett, J.; Bylander, J. High Kinetic Inductance Nb N Nanowire Superinductors. _Phys. Rev. Applied_ 2019, _11_ , 044014
  * Harvey-Collard et al. 2020 Harvey-Collard, P.; Zheng, G.; Dijkema, J.; Samkharadze, N.; Sammak, A.; Scappucci, G.; Vandersypen, L. M. K. On-Chip Microwave Filters for High-Impedance Resonators with Gate-Defined Quantum Dots. _Phys. Rev. Appl._ 2020, _14_ , 034025
  * Polini et al. 2022 Polini, M. et al. Materials and Devices for Fundamental Quantum Science and Quantum Technologies. 2022