# Critical scaling law for the deposition efficiency of inertia-driven
particle collisions with a cylinder in high Reynolds number air flow

Matthew R. Turner<EMAIL_ADDRESS>http://personal.maths.surrey.ac.uk/st/M.Turner/ School of Mathematics and
Physics, University of Surrey, Guildford, GU2 7XH, United Kingdom  Richard P.
Sear<EMAIL_ADDRESS>https://richardsear.me/ School of Mathematics and
Physics, University of Surrey, Guildford, GU2 7XH, United Kingdom

###### Abstract

The Earth’s atmosphere is an aerosol, it contains suspended particles. When
air flows over an obstacle such as an aircraft wing or tree branch, these
particles may not follow the same paths as the air flowing around the
obstacle. Instead the particles in the air may deviate from the path of the
air and so collide with the surface of the obstacle. It is known that particle
inertia can drive this deposition, and that there is a critical value of this
inertia, below which no point particles deposit. Particle inertia is measured
by the Stokes number, St. We show that near the critical value of the Stokes
number, Stc, the amount of deposition has the unusual scaling law of
exp(-1/(St-Stc)1/2). The scaling is controlled by the stagnation point of the
flow. This scaling is determined by the time for the particle to reach the
surface of the cylinder varying as 1/(St-Stc)1/2, together with the distance
away from the stagnation point (perpendicular to the flow direction)
increasing exponentially with time. The scaling law applies to inviscid flow,
a model for flow at high Reynolds numbers. The unusual scaling means that the
amount of particles deposited increases only very slowly above the critical
Stokes number. This has consequences for applications ranging from rime
formation and fog harvesting to pollination.

##  I Introduction

The Earth’s atmosphere is an aerosol, in that it contains suspended particles,
with sizes up to tens of micrometres [1]. For example, clouds and fogs are
aerosols of water droplets suspended in air. When air flows over an obstacle
such as an aircraft wing or tree branch, the suspended particles may follow
the air flow around the obstacle, or they may deposit on the surface of the
obstacle. For particles tens of micrometres in diameter, deposition on the
surface of the obstacle is largely due to the particle’s inertia. While the
air flow curves to move around the obstacle, the particle’s inertia means the
particle tends to move in straight lines, and so collides with the obstacle,
see Fig. 1. Here we study these inertia-driven collisions.

Deposition of particles from flowing air onto obstacles occurs in many
contexts. These include filtration [2, 3, 4], harvesting water from fog [5, 6,
7], pollination [8, 9], and rime (ice) formation [10, 11, 12]. It has been
studied for approximately a hundred years. A lot of the earliest work
considered applications where the air was flowing rapidly, at large Reynolds
numbers, which is relevant to applications such as rime forming on the leading
edges of aircraft wings. The large Reynolds number limit is also the limit we
consider here. Rime forms by water droplets below the freezing temperature
($T=0$$~{}^{\circ}$C) depositing on a surface and then freezing. This can coat
a wing with ice, risking a crash.

In 1931 Albrecht [13] found a threshold in the inertia, below which no (point)
particles deposited on the obstacle. A minimum amount of inertia is needed
before any particles are deposited. Then in the 1940s, first G.I. Taylor [14]
(also in the scientific papers of G.I. Taylor[15]), and then Langmuir and
Blodgett [16] calculated this critical value of the inertia. Here we confirm
this result, and build on it. We determine the critical scaling of the
deposition efficiency, and show how this relates to the flow field around the
obstacle.

There has been considerable work on this problem since the 1940s [17, 18],
motivated by its many applications, but the critical scaling has not been
studied before, for the high Reynolds number flow field considered by Taylor
[14, 15], and by Langmuir and Blodgett [16]. How the threshold varies with
Reynolds number has been studied by Phillips and Kaye[19], and Araújo et
al.[20] determined the critical scaling for zero-Reynolds-number flow.

Most of the work has been computational or theoretical, as experiments on
collisions in high- Reynolds-number flows are challenging. However, Wong and
coworkers did obtain some experimental data on deposition efficiency at
Reynolds numbers of hundreds [21]. This was for an aerosol with particles with
a narrow distribution of sizes. They found no measurable deposition below a
threshold near that predicted by Taylor [14, 15], and by Langmuir and Blodgett
[16]. So these experiments agree with theoretical predictions (for a
simplified model) that a threshold exists. Makkonen and coworkers [22, 23, 24]
have measured ice deposition on cylinders. This is for the typical case in the
environment, where the droplets have a broad range of sizes, which complicates
comparison with theoretical predictions. Here we mostly consider aerosols of
identical particles but we do look at the affect of a distribution of particle
sizes, in order to compare with this work.

Figure 1: Plot of the cross-section of the cylinder (yellow), the flow field
(streamlines in blue), and three trajectories. The green and black
trajectories collide with the cylinder, while the red trajectory misses. Here
$\mathrm{St}=0.7$ and the initial conditions are at Cartesian coordinate
$x_{C}(t=0)=-10$, with initial $y_{C}$ values given in the key. The black
trajectory is the one that defines the edge of the region where particles
collide, and so its value of $y_{C}=0.2793=\eta_{d}$. N.B. Numerical errors
mean we cannot determine $\eta_{d}$ to four significant figures, so the
$y_{C}$ value for black curve should not be taken as accurate to four figures.

###  I.1 Obstacles in high Reynolds number flow

The air flow over an aircraft wing is fast, speeds $U$ of order
$100\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$
($\simeq$360\text{\,}\mathrm{km}\text{\,}{\mathrm{h}}^{-1}$$). This speed,
combined with a wing leading edge radius $R$ of order
$10\text{\,}\mathrm{cm}$, means that the Reynolds number of the flow over the
wing is $\mathrm{Re}=UR/\nu\sim 10^{6}\gg 1$ with $\nu\simeq
10^{-5}$\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-2}$$, the kinematic
viscosity of air. So, following Langmuir and Blodgett and many others, we
approximate the airflow over a wing by inviscid, incompressible flow over an
infinite cylinder, where we have a simple analytic expression for the flow
field [25]. This is a simple model of high-Reynolds-number flow.

The cylinder is of radius $R$, with its axis along the $z_{C}$-axis, directed
into the page in Fig. 1. The cylinder is taken to be at rest in the frame of
reference. Far from the cylinder, flow is in the $x_{C}$ direction: $U{\bf
i}$. Note that we denote the Cartesian coordinates as $(x_{C},y_{C},z_{C})$,
with corresponding unit vectors ${\bf i},~{}{\bf j},~{}{\bf k}$, because below
we will use $x$ and $y$ to indicate distances from the stagnation point. In
cylindrical polar coordinates, the flow field
$\bm{\mathbf{u}}(\bm{\mathbf{r}})$ is given by

$\displaystyle\frac{\bm{\mathbf{u}}}{U}=\left(1-\frac{R^{2}}{r^{2}}\right)\cos(\theta)\widehat{\bm{\mathbf{r}}}-\left(1+\frac{R^{2}}{r^{2}}\right)\sin(\theta)\widehat{\bm{\mathbf{\theta}}},$
(1)

where $(r,\theta)$ are plane polar coordinates in the $(x_{C},y_{C})$-plane.
Streamlines of the flow field are illustrated in Fig. 1.

In the remainder of this paper, we will set the radius of the cylinder $R=1$
and the flow field speed $U=1$. So for example, both the Cartesian coordinate
$x_{C}$ and the distance along the Cartesian coordinate from the stagnation
point, $x$, are in units such that the cylinder radius $R=1$. This choice also
means that time is measured in units such that it takes unit time to move a
distance of the cylinder’s radius, when moving at speed $U$.

###  I.2 Particles in air flowing around an obstacle: The effect of inertia on
deposition

We approximate the particles by point particles, they are about ten thousand
times smaller than the wing leading edge. Point particles that follow the
streamlines of fluid flow perfectly never collide with the obstacle. However,
if the particles have inertia then when the air flow changes direction to flow
around the obstacle, the particle’s inertia may cause it to go straight on
potentially leading it to crash into the obstacle. So particle inertia can
cause collisions.

The inertia of a particle is quantified by its Stokes number $\mathrm{St}$,
defined for a particle of mass $m_{p}$ by

$\mathrm{St}=m_{p}UB\big{/}R$ (2)

with $B$ the particle mobility and $R$ the lengthscale of the obstacle. The
Stokes expression for the mobility is $B=1/(6\pi\eta a_{p})$ with $\eta$ the
dynamic viscosity of air, and $a_{p}$ the radius of the particle. The Stokes
number is the dimensionless ratio between the inertia of a particle – which
tends to cause the particle to move in straight lines – and the friction
between the particle and the surrounding air – which tends to make the
particle follow streamlines. For a wing of width $0.1\text{\,}\mathrm{m}$, the
Stokes number is of order one for droplets micrometres in diameter, so cloud
droplets of size micrometres, and tens of micrometres will deposit on the
wing.

In the $\mathrm{St}\to\infty$ limit, particles move in straight lines and so
in that limit a cylinder of radius $R$ sweeps out a strip of air of thickness
$2R$, collecting all the particles in this strip. This allows us to define the
deposition or collection efficiency

$\eta_{d}=\left.\begin{array}[]{l}\mbox{maximum displacement from cylinder
axis}\\\ \mbox{perpendicular to the flow direction, for}\\\ \mbox{which a
particle collides with cylinder surface}\end{array}\middle/R\right.$

which varies from zero when no particles collide, to one when
$\mathrm{St}\to\infty$. The displacement is taken far upstream of the
cylinder. Calculation of $\eta_{d}$ is done by starting particle trajectories
far upstream of the cylinder at varying values of the displacement $y_{C}$
normal to the flow direction. Then $\eta_{d}$ is defined by the largest
initial displacement along $y_{C}$, for which the particle collides. This is
illustrated by the black trajectory in Fig. 1.

Numerical results for $\eta_{d}$ as a function of $\mathrm{St}$, together with
the fit of Langmuir and Blodgett [16], are shown in Fig. 2. Note that at small
values of St, the collection efficiency increases very slowly with increasing
inertia. This is what we will explain here. It follows directly from the high
Reynolds number flow field. In the low-Reynolds-number limit, where the flow
field is very different, $\eta_{d}$ increases much more rapidly above the
critical Stokes number [20, 3, 26]. Langmuir and Blodgett’s fit is in Appendix
A, and the details of our numerical calculations are in Appendix B. See
Appendix B.1 for the details of the fit in Fig. 2.

Figure 2: Plot of the deposition or collection efficiency $\eta_{d}$ as a
function of the Stokes number. Shown are numerical results (blue circles), the
function of Langmuir and Blodgett [16] (green curve) and our fit to the region
of small $\delta=\mathrm{St}-\mathrm{St}_{c}$ (dashed red curve).

##  II Newton’s equation for a particle in flowing air

For a particle suspended in flowing air we assume that the only force on the
particle is the friction with the surrounding air, which is taken to be
proportional to the difference between particle’s velocity $\bm{\mathbf{v}}$
and the local flow velocity $\bm{\mathbf{u}}$. This force is taken to act on
the particle’s centre of mass. Then Newton’s equation for the particle motion
is

$\mathrm{St}\frac{d\bm{\mathbf{v}}}{dt}=-(\bm{\mathbf{v}}-\bm{\mathbf{u}}).$
(3)

In cylindrical polar coordinates, this is

$\displaystyle\mathrm{St}\left(\ddot{r}-r\left(\dot{\theta}\right)^{2}\right)=-\left(\dot{r}-u_{r}\right),$
(4a)
$\displaystyle\mathrm{St}\left(\ddot{\theta}+\frac{2\dot{r}\dot{\theta}}{r}\right)=-\left(\dot{\theta}-u_{\theta}\right),$
(4b)

where
$\bm{\mathbf{u}}=u_{r}\hat{\bm{\mathbf{r}}}+u_{\theta}\hat{\bm{\mathbf{\theta}}}$.

##  III The cylinder’s forward stagnation point

We are interested in the behaviour near the critical Stokes number of Langmuir
and Blodgett [16]. Here, particle trajectories pass close to the stagnation
point at the front of the cylinder. This stagnation point is at $r=1$ and
$\theta=\pi$, see Fig. 1. So we will study behaviour near this stagnation
point. We start by changing variables to the distance to contact with the
cylinder $x=r-1\ll 1$, and the angle from the angle of the stagnation point
$y=\pi-\theta\ll 1$. Note that $x$ and $y$ are not the conventional Cartesian
coordinates, which we denote by $x_{C}$ and $y_{C}$.

In these new coordinates, Newton’s equations for the particle, Eq. (4), become

$\displaystyle\ddot{x}-(1+x)\dot{y}^{2}$
$\displaystyle=-\frac{\dot{x}-u_{r}}{\mathrm{St}},$ (5a)
$\displaystyle\ddot{y}+\frac{2}{1+x}\dot{x}\dot{y}$
$\displaystyle=-\frac{\dot{y}+u_{\theta}}{\mathrm{St}}.$ (5b)

Note that $\dot{y}=-\dot{\theta}$ and $\ddot{y}=-\ddot{\theta}$ etc. Near the
stagnation point, we can expand the flow field in Eq. (1) as a series in $x$
and $y$

$\displaystyle\bm{\mathbf{u}}$ $\displaystyle=$
$\displaystyle\left(-2x+3x^{2}+\mathcal{O}({\rm
cubic~{}terms})\right)\hat{\bm{\mathbf{r}}}$ (6)
$\displaystyle+\left(-2y+2xy+\mathcal{O}({\rm
cubic~{}terms})\right)\hat{\bm{\mathbf{\theta}}}.$

It is worth noting that as we do not have stick boundary conditions here,
$u_{\theta}$ is not zero at the cylinder surface, except at the stagnation
point; $u_{r}$ is zero at the surface is because there is no flow into the
cylinder.

We substitute the flow field of Eq. (6) into the equation for the particle
trajectory (Eq. (5)). Then if we keep only linear and quadratic terms we
obtain

$\displaystyle\ddot{x}-\dot{y}^{2}$
$\displaystyle=-\frac{1}{\mathrm{St}}\left(\dot{x}+2x-3x^{2}\right),$ (7a)
$\displaystyle\ddot{y}+2\dot{x}\dot{y}$
$\displaystyle=-\frac{1}{\mathrm{St}}\left(\dot{y}-2y+2xy\right),$ (7b)

We now solve these equations to find out which are the trajectories of
particles that collide with the cylinder.

##  IV Particle trajectories on axis, and the critical Stokes number

We start by considering particle trajectories precisely on axis, where
$\theta=\pi$ and $y=0$. This will enable us to determine the critical value of
the Stokes number, below which no particles collide with the cylinder. We
follow Taylor [14, 15], and Langmuir and Blodgett [16] here. On axis the
system reduces to a one dimensional problem.

Near the stagnation point $x\ll 1$, and we can approximate the flow by keeping
only the leading order terms in $x$. Then Eq. (7a) is just

$\mathrm{St}\ddot{x}=-\dot{x}-2x,$ (8)

which as both Taylor, and Langmuir and Blodgett realised, is just the
differential equation for damped simple harmonic motion (SHM). It has
solutions

$x(t)=A_{0}\exp\left(\lambda_{1}t\right)+B_{0}\exp\left(\lambda_{2}t\right),$
(9) $\lambda_{1},\lambda_{2}=\frac{-1\pm\sqrt{1-8\mathrm{St}}}{2\mathrm{St}},$
(10)

where $A_{0}$ and $B_{0}$ are fixed by the initial conditions, for example
$x(t=0)=x_{0}$ and $\dot{x}(t=0)=u_{0}$. See Appendix C.

For $0<\mathrm{St}<1/8$ both $\lambda$ are real and negative (overdamped SHM
solutions), and so the particle approaches the cylinder surface at a speed
which decays exponentially. There is only a collision in the $t\to\infty$
limit, i.e., no collision at finite time. But if $\mathrm{St}>1/8$ then we
have complex $\lambda$ (underdamped SHM solutions), and the collision will
occur in finite time. The critical value of the Stokes number is therefore
$\mathrm{St}_{c}=1/8$. We can define the distance from this critical value as

$\delta=\mathrm{St}-\mathrm{St}_{c}.$ (11)

Note that there are initial conditions for Eq. (9), for which the particle
collides in finite time[27]. However, they do not appear to be physically
relevant, as Ingham et al.[27] discuss.

###  IV.1 Time to collision, on axis

The collision occurs when $x(t_{COLL})=0$ and is set by the angular frequency
— the imaginary part of Eq. (10) —
$\omega=\sqrt{1-8\mathrm{St}}/(2\mathrm{St})\simeq 8\sqrt{2}\delta^{1/2}$. As
expected (see Appendix C for details) the time to collide $t_{COLL}$ is half
the period

$t_{COLL}\simeq\frac{\pi}{8\sqrt{2}}\delta^{-1/2}~{}~{}~{}~{}~{}~{}\delta\ll
1.$ (12)

As the critical Stokes number is approached from above, the time to reach the
cylinder surface and collide diverges as $1/\delta^{1/2}$. Numerical
calculations using the full flow field agree with this observation, see
Appendix B.1.

Figure 3: Plot of the trajectories in Fig. 1 but in the $xy$ plane. The green
and black trajectories collide with the cylinder, while the red trajectory
misses. Here $\mathrm{St}=0.7$ and the initial conditions are at Cartesian
coordinate $x_{C}(t=0)=-10$, with initial $y_{C}$ values given in the key.
Note that the trajectory which just collides does so tangentially to the
surface.

##  V Particle trajectories off axis, and the critical scaling of the
deposition efficiency

We now consider the full two-dimensional case, near the stagnation point. Off
axis we also need an equation for $y$. Retaining only the leading order linear
terms in Eq. (7b) for $y$ leads to

$\mathrm{St}\ddot{y}=-\dot{y}+2y,$ (13)

in which there is no coupling with the $x$ direction. This equation is almost
the damped SHM equation again, but in SHM language, the ‘force’ term has the
opposite sign, so it is not a restoring force but drives exponential growth of
$y$. The solution is

$y(t)=C_{0}\exp\left(\mu_{1}t\right)+D_{0}\exp\left(\mu_{2}t\right)$ (14)

with

$\mu_{1},\mu_{2}=\frac{-1\pm\sqrt{1+8\mathrm{St}}}{2\mathrm{St}}$ (15)

Now, $\mu_{1}>0$ and $\mu_{2}<0$ so at long times the $\mu_{1}$ solution
dominates; it increases exponentially while the other solution decays to zero.
The constants $C_{0}$ and $D_{0}$ are related to the initial conditions. If
$y(t=0)=y_{0}$ and $\dot{y}(t=0)=v_{0}$, then

$C_{0}=\frac{v_{0}-\mu_{2}y_{0}}{\mu_{1}-\mu_{2}}~{}~{}~{}{\rm
and}~{}~{}~{}D_{0}=-\frac{v_{0}-\mu_{1}y_{0}}{\mu_{1}-\mu_{2}}.$ (16)

Having determined the leading order behaviour of $y$, we return to Eq. (7a)
for $x$. The leading order $y$ term in Eq. (7a) is the $\dot{y}^{2}$ term.
Retaining only this term, we have

$\ddot{x}+\frac{\dot{x}+2x}{\mathrm{St}}=\dot{y}^{2}.$ (17)

The new term acts like an effective force that pushes the particle away from a
collision.

The general solution to Eq. (17) is given in Eq. (37) in Appendix D. We take
this equation, expand for $\delta\ll 1$, keep only leading-order terms, and
then set $x=0$ at the collision time of Eq. (12). This yields

$\displaystyle x(t_{COLL})$ $\displaystyle=$ $\displaystyle-
A_{0}\exp\left[-\frac{\pi}{2\sqrt{2}}\delta^{-1/2}\right]$ (18)
$\displaystyle+E_{0}\exp\left[\frac{(\sqrt{2}-1)\pi}{\sqrt{2}}\delta^{-1/2}\right]=0,$

with $E_{0}=C_{0}^{2}(11-6\sqrt{2})/49$. This equation is a function of
$\delta$ and the initial boundary conditions on $x$ and $y$.

Equation (18) has two terms. The first term is negative and is the same for
the on-axis case. It is this term that results in a collision after a time
$\sim 1/\delta^{1/2}$. The second term is new for the off-axis case, it
increases exponentially with $\delta^{-1/2}$, and the prefactor scales with
$C_{0}^{2}\sim(v_{0}-\mu_{2}y_{0})^{2}$, i.e., with the initial conditions on
$y$.

We want to estimate $\eta_{d}$, which is set by the largest displacement
$y_{0}$ for which the particle collides with the cylinder, which is when Eq.
(18) is satisfied. Near the critical Stokes number $y_{0}$ will be small, and
the exponential variation of the terms in Eq. (18) suggests it needs to be
exponentially small, so we write $C_{0}=\zeta\exp[\alpha(\delta)]$ where
$\zeta$ is a constant to leading order in $\delta$, and $\alpha$ is a function
of $\delta$.

Then we insert $C_{0}=\zeta\exp(\alpha)$ into Eq. (18). For small $\delta$
solutions are only possible when the exponential exponents are equal. Equating
the exponents of the two terms gives us

$-\frac{\pi}{2\sqrt{2}}\delta^{-1/2}=2\alpha+\frac{(\sqrt{2}-1)\pi}{\sqrt{2}}\delta^{-1/2},$
(19)

and

$\alpha=-\frac{(4-\sqrt{2})}{8}\pi\delta^{-1/2}\simeq-1.015\delta^{-1/2}.$
(20)

Collisions only occur at small $\delta$ when $y_{0}$ scales as
$\exp(-1.015\delta^{-1/2})$. This sets the width of the strip of air over
which particles collide with the cylinder, and so the deposition efficiency is

$\eta_{d}\sim\left\\{\begin{array}[]{cc}0&\delta\leq 0\\\
\exp\left[-\frac{(4-\sqrt{2})}{8}\pi\delta^{-1/2}\right]&0<\delta\ll
1\end{array}\right.$ (21)

Above the critical Stokes number, the collection efficiency has the unusual
scaling $\exp(-1/\delta^{1/2})$. This means that the deposition efficiency
increases only slowly above the critical Stokes number, see Fig. 2. Numerical
calculations with the full flow field confirm the result. They are in Appendix
B.1. Limited experimental data [21] near $\mathrm{St}_{c}$ means that a
quantitative comparison with our scaling result is not possible.

The scaling is a direct consequence of the form of the flow field near the
stagnation point, Eq. (6). From the fact that the flow field varies linearly
with $x$ and $y$. The flow field leads to the time to collision scaling as
$1/\delta^{1/2}$, and to the $y$ coordinate increasing exponentially with
time, which in turn gives us Eq. (21).

Numerical results for particle trajectories as functions of $x$ and $y$ are
plotted in Fig. 3. The trajectory that just collides, and so defines
$\eta_{d}$, is in black. Note the rapidly increasing $y$ as the collision is
approached, and that the collision occurs tangential to the cylinder surface.

##  VI Comparison with particle deposition in Stokes flow

Our results are for a flow field that neglects viscosity, in effect the
infinite-Reynolds-number limit. The opposite limit is where inertia is
negligible and viscosity dominates. This is the zero-Reynolds-number or
Stokes-flow limit, and it is relevant to the filtration of aerosol particles
from air by fibrous filters, where the Reynolds number is small [3, 2, 28].
For particles in Stokes flow, there is also a critical value of the Stokes
number, but the scaling of the deposition efficiency is
$\eta_{d}\sim\delta^{1/2}$, which is completely different scaling. Langmuir
and Blodgett suspected that there were critical Stokes numbers for a number of
different flow fields including spheres in Stokes flow, but this scaling was
first found for spheres in Stokes flow by Araújo et al.[20].

The particle trajectories are very different for the inviscid and Stokes flow
fields. For example, with the slip boundary conditions in the inviscid limit,
the trajectories where the particle just collides with the cylinder surface do
so tangentially, see Fig. 3. This is because the particle’s velocity normal to
the cylinder surface tends to zero at the collision, while the tangential
component is non-zero. Particles in Stokes flow do not collide tangentially
because here the stick boundary conditions mean that the tangential fluid-flow
velocity is zero at contact [26].

At intermediate Reynolds numbers, there will be a boundary layer of thickness
$\sim 1/$Re1/2. Within this boundary layer the flow field is approximately
viscous flow, outside it is closer to inviscid flow. So as the Reynolds nunber
is increased the boundary layer thins and the behaviour changes continuously
from the Stokes to the inviscid limit. Work by Robinson and coworkers suggests
that the deposition behaviour changes smoothly between the limits [26].

Figure 4: Plot of the mean deposition or collection efficiency,
$\overline{\eta_{d}}$, as a function of the median Stokes number. Shown are
numerical results for monodisperse droplets (blue circles), and for
polydisperse droplets whose Stokes numbers obey a log-normal distribution with
width parameter $\sigma=1/2$ (orange curve).

##  VII Deposition from an aerosol of particles with a range of diameters

In natural aerosols, the particles are almost never all of the same diameter.
A broad distribution of particle sizes is typical [22, 24, 1]. Within our
simplified model, each particle is characterised by a single parameter: the
Stokes number. A distribution of particle sizes results in a distribution of
Stokes numbers. The Stokes number of a particle is given by Eq. (2). When mass
$m_{p}\propto a_{p}^{3}$ and the Stokes mobility $B\propto 1/a_{p}$, the
scaling of the Stokes number is St $\propto a_{p}^{2}$, i.e., the Stokes
number scales with the square of the particle radius. Thus a distribution of
particle radii gives a distribution of Stokes numbers. An aerosol of particles
with a range of radii will have particles with a distribution of Stokes
numbers: $p(\mathrm{St})$.

We need a model probability distribution function, $p$. We use the standard
log-normal distribution

$p(\mathrm{St};\mathrm{St}_{\mathrm{MED}},\sigma)=\frac{1}{\mathrm{St}\sigma(2\pi)^{1/2}}\exp\left[-\frac{(\ln\mathrm{St}-\mu)^{2}}{2\sigma^{2}}\right]$
(22)

Here, the median Stokes number $\mathrm{St}_{\mathrm{MED}}=\exp(\mu)$ provides
a measure of the typical Stokes number. The width of the distribution is set
by $\sigma$: the ratio of the standard deviation to the median value is
$\exp(\sigma^{2}/2)(\exp(\sigma^{2})-1)^{1/2}$.

Within our model, the particles are independent. So the mean deposition
efficiency, $\overline{\eta_{d}}$, for the aerosol is simply

$\overline{\eta_{d}}(\mathrm{St}_{\mathrm{MED}},\sigma)=\int\eta(\mathrm{St})p(\mathrm{St};\mathrm{St}_{\mathrm{MED}},\sigma){\rm
dSt}$ (23)

In Fig. 4, we compare the mean deposition efficiency as a function of median
Stokes number, of monodisperse particles, and particles with a log-normal
distribution with width parameter $\sigma=1/2$. This distribution has a ratio
of standard deviation to the median of $0.60$. Near the critical value of the
Stokes number, the convolution of deposition efficiency with the distribution
of Stokes numbers of the particles smooths out the transition at
$\mathrm{St}_{c}$. For a distribution of particle sizes and hence and Stokes
numbers, the deposition efficiency is never zero, unless all particles have
Stokes numbers less than $\mathrm{St}_{c}$. Thus when the median Stokes number
is below $\mathrm{St}_{c}$ polydisperse aerosols have larger deposition
efficiencies than monodisperse ones. But regardless of polydispersity the
deposition efficiency tends to one at larges values of the median Stokes
number. So at large median Stokes numbers a broad distribution of particle
sizes has little effect on the mean deposition efficiency.

Figure 5: A quantity $M$, which is proportional to the mass deposited, plotted
as a function of the median Stokes number. Shown are numerical results for
monodisperse droplets (dark blue circles) and for polydisperse droplets whose
Stokes numbers obey a log-normal distribution with width parameter
$\sigma=1/2$ (magenta curve).

In experiment, the total mass deposited on an obstacle can be measured [22,
24]. As the mass of a particle scales as $a_{p}^{3}$, the mass of a particle
scales as St3/2. So the three-halves moment of the efficiency times the
distribution of particle Stokes numbers gives us a quantity proportional to
the total mass deposited. We call this quantity $M$:

$M(\mathrm{St}_{\mathrm{MED}},\sigma)=\int\eta(\mathrm{St})p(\mathrm{St};\mathrm{St}_{\mathrm{MED}},\sigma)\mathrm{St}^{3/2}{\rm
dSt}$ (24)

We have plotted this deposited mass as a function of median Stokes number in
Fig. 5. As with the deposition efficiency in Fig. 4, when a range of particle
sizes are present, the transition is smoothed over. Mass is deposited at all
values of the median Stokes number. This is consistent with the work of
Makkonen and coworkers [22, 24], who find that some mass is deposited under
all their experimental conditions. However, the mass deposited is always
larger than for monodisperse droplets with the same median Stokes number. The
total mass of aerosol particles is larger in the polydisperse case. (at the
same median Stokes number). In addition, the large Stokes number tail of the
distribution contributes large amounts to the mass deposited as here the
deposition efficiency is high and these large particles contribute large
amounts to the total mass deposited.

##  VIII Conclusion

We have built on Taylor’s [14, 15], and Langmuir and Blodgett’s work [16] to
quantitatively understand the behaviour of the deposition efficiency
$\eta_{d}$ near the critical value of the Stokes number. Just above the
critical value, $\eta_{d}\sim\exp(-1/\delta^{1/2})$. This unusual scaling
comes from the fact that for an off-axis particle trajectory, the particle’s
displacement from the cylinder’s axis increases exponentially with time, while
the time to reach the cylinder scales as $1/\delta^{1/2}$. Thus unless the
initial displacement from the axis is exponentially small, the particle misses
the cylinder. The scaling follows directly from the (inviscid) flow field we
used; the behaviour of particles in Stokes flow with stick boundary conditions
is very different [20, 3, 26]. The lesson here is that particle deposition
from flowing air is very sensitive to details of the flow field, especially
near the stagnation point. The slow increase in deposition efficiency follows
from the slip boundary conditions. Particle deposition has varied
applications, from rime formation on aircraft wings v to pollination [8, 9].
The careful analysis of particle trajectories near the forward stagnation
point will greatly contribute to the understanding of this issue.

Here, we have studied a cylinder, but Langmuir and Blodgett [16] also found a
critical value of the Stokes number for spheres in inviscid, axisymmetric
flow. This was at the slightly lower value of $\mathrm{St}_{c}=1/12$. It is
straightforward to show that the scaling laws found here for a cylinder are
the same (up to multiplicative constants) for a sphere. The details are in
Appendix E. Our results for the scaling were obtained by expanding around the
stagnation point, it is possible that they are general to at least most
inviscid flows impinging on a locally parabolic surface in two or three
dimensions.

##  IX Supplementary Material

The supplementary material is a Python Jupyter notebook that performs all
numerical calculations, and produces all figures in this work.

###### Acknowledgements.

It is a pleasure to thank Josh Robinson and Patrick Warren for inspiring and
useful discussions. The data that support the findings of this study are
available within the article and in a Python Jupyter notebook in supplementary
material.

##  Appendix A Langmuir and Blodgett expression for deposition efficiency

The expression used by Langmuir and Blodgett [16] is

$\eta_{d}=\left\\{\begin{array}[]{cc}0&\mathrm{St}\leq 1/8\\\
0.466\left[\log_{10}(8\mathrm{St})\right]^{2}&1/8<\mathrm{St}<1.1\\\
\mathrm{St}/(\mathrm{St}+\pi/2)&1.1<\mathrm{St}\end{array}\right.$ (25)

Note that the functional form near $\mathrm{St}=1/8$ is incorrect but that it
is quite close to numerical calculations. This expression, and closely related
ones, have been and are used extensively to model ice accretion on the leading
edges of surfaces moving through air [10, 17, 18, 7]. Finstad et al.[18, 17]
summarise work up to 1998 on improving the approximation of Eq. (25). None of
the improvements has the correct functional form for small
$\delta=\mathrm{St}-1/8$ but it should be said the work was focused on
producing empirical expressions that are accurate over a wide range of values
of the Stokes number, not on obtaining the correct value near the dynamical
transition.

##  Appendix B Details of numerical calculations

The numerical calculations were done with a Python Jupyter notebook available
in the supplementary material. Particle trajectories were calculated as
follows. Particles were started at Cartesian coordinates upstream of the
cylinder $(x_{C0},y_{C0})$, i.e., a distance $x_{C0}$ along the $x_{C}$-axis
(the flow direction), and displaced a perpendicular distance $y_{C0}$ off
axis. All calculations here are for $x_{C0}=-10$; we checked that increasing
the distance upstream further to $-20$ had very little effect. Trajectories
were obtained by integrating Newton’s equation of motion ((3)), and checking
for collisions, which occur whenever $r=1$.

Figure 6: Plot of time to collision, as a function of
$\delta=\mathrm{St}-\mathrm{St}_{c}$. The points are numerical results, and
the lines are fits of function $A+B\delta^{-1/2}$. Blue points and cyan line
are for on axis collisions, with best fit values $A=9.78$ and $B=0.309$.
Orange points and red line are for collisions at the edge of the deposition
zone, with $A=9.92$ and $B=0.309$. Note that both values for $B$ are close to
our prediction of $t_{COLL}=\pi/(8\sqrt{2}\delta^{1/2})=0.278/\delta^{1/2}$.

###  B.1 Verification of scaling found by analytical mathematics, by
comparison to numerical calculations

On axis, our analytical results predict that as the critical Stokes number is
approached from above, the time to reach the surface of the cylinder and
collide, scales as $1/\delta^{1/2}$. There is another contribution, which is
the time for the particle to go from its initial position upstream to a point
where $x\ll 1$. That contribution is approximately equal to initial
displacement$/U$.

So, we fit a function of the form $A+B\delta^{-1/2}$ to numerical results for
the time to collision. The numerical results and fit are shown in Fig. 6. We
see excellent agreement.

We perform the same analysis for the trajectories that define $\eta_{d}$, at
the edge of the strip of air over which collisions occur. The results are in
Fig. 6, and again the agreement is excellent.

We also predict that the deposition efficiency has the scaling form

$\ln\eta_{d}=A+B\frac{1}{\delta^{1/2}}~{}~{}~{}~{}\delta\ll 1.$ (26)

In Fig. 7 we compare the results of numerical calculations with a fit of this
form. Again the agreement is excellent.

Figure 7: Plot of $\ln(\eta_{d})$ as a function of $1/\delta^{1/2}$, to show
$\exp(-1/\delta^{1/2})$ scaling of the deposition efficiency. Points are
numerical results, line is fit of function $A+B(1/\delta^{1/2})$ with best fit
values $A=-0.193$ and $B=-1.118$. Note that value for $B$ is close to our
prediction of $\alpha=-1.015/\delta^{1/2}$ (Eq. (20)).

A final check is on the match between the numerics and the analytics; the
analytical mathematics relies on an expansion in $x$ and $y$ and so is only
valid near the stagnation point where $x,y\ll 1$. We are interested in the
deposition efficiency $\eta_{d}(\delta)$. The deposition efficiency is set by
the trajectory with the largest initial Cartesian $y_{C}$ coordinate, call it
$y_{C0}$, at large and negative Cartesian $x_{C}$ (far upstream of the
cylinder) for which a collision occurs. There the particle velocity is taken
to be that in flow field, which is close to $U$ along the Cartesian $x_{C}$
axis.

What the analytic mathematics uses is the initial boundary condition on $y$,
as it appears as $C_{0}$, which has the scaling

$C_{0}\sim v_{0}-\mu_{2}y_{0}\sim\exp(\alpha)$ (27)

and in effect shows that $C_{0}\sim\exp(-1/\delta^{1/2})$. Here $y_{0}$ and
$v_{0}$ are the initial conditions on $y$ and $\dot{y}$, respectively. To
complete the link between the scaling we obtained for $C_{0}$ and $\eta_{d}$
we need to show that $C_{0}$ varies smoothly with the value of $y_{C0}$. We do
this in Fig. 8. There we have plotted $y$, $\dot{y}$ and $\dot{y}-\mu_{2}y$ at
$x=0.1$, as functions of $y_{C0}$. These are all obtained from numerical
calculations. We see that they all vary smoothly with $y_{C0}$, so the scaling
found for $C_{0}$ will translate $\eta_{d}$. In Fig. 8 we selected $x=0.1$ as
a value less than one but not $\ll 1$, plots with $x=0.05$ and $x=0.2$, are
similar.

Figure 8: Plot of $y$, $\dot{y}$ and their combination $\dot{y}-\mu_{2}y$
which sets $C_{0}$, as functions of the initial off-axis deviation in
Cartesian coordinates. All trajectories start at $x_{C}=-10$, and end at
$x=0.1$. The Stokes number is $\mathrm{St}=1/4$.

##  Appendix C Solution for on-axis case

The solution of Eq. (8) is

$x(t)=A_{0}\exp\left(\lambda_{1}t\right)+B_{0}\exp\left(\lambda_{2}t\right)$
(28)
$\lambda_{1},\lambda_{2}=\frac{-1\pm\sqrt{1-8\mathrm{St}}}{2\mathrm{St}},$
(29)

where $A_{0}$ and $B_{0}$ are fixed by the initial conditions. With
$x(t=0)=x_{0}$ and $\dot{x}(t=0)=u_{0}$, $A_{0}$ and $B_{0}$ are given by

$A_{0}=\frac{\lambda_{2}x_{0}-u_{0}}{\lambda_{2}-\lambda_{1}}~{}~{}~{}{\rm
and}~{}~{}~{}B_{0}=\frac{u_{0}-\lambda_{1}x_{0}}{\lambda_{2}-\lambda_{1}}.$
(30)

For $\mathrm{St}<1/8$ both $\lambda$ are real (overdamped SHM solutions), and
so the particle approaches the cylinder surface $(x=0)$ at a speed which
decays as $\exp(-t)$, and hence there is only a collision in the $t\to\infty$
limit, i.e., no collision at finite time. A similar phenomena occurs in the
$\mathrm{St}=1/8$ case. But if $\mathrm{St}>1/8$ then we have complex
$\lambda$ (underdamped SHM solution), and the collision will occur in finite
time. In this case the solution can be written as

$\displaystyle x(t)$ $\displaystyle=$
$\displaystyle\exp\left[-t/(2\mathrm{St})\right]\left(A_{0}\cos\left[\frac{\sqrt{8\mathrm{St}-1}}{2\mathrm{St}}t\right]\right.$
(31)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.+B_{0}\sin\left[\frac{\sqrt{8\mathrm{St}-1}}{2\mathrm{St}}t\right]\right),$

or by using trigonometric identities

$x(t)=\psi\exp\left[-t/(2\mathrm{St})\right]\cos\left[\frac{\sqrt{8\mathrm{St}-1}}{2\mathrm{St}}t+\phi\right],$
(32)

where $\psi$ is an amplitude and $\phi$ is a phase which are related to the
constants $A_{0}$ and $B_{0}$, which are given by

$A_{0}=x_{0}~{}~{}~{}{\rm
and}~{}~{}~{}B_{0}=\frac{2\mathrm{St}\,u_{0}+x_{0}}{\sqrt{8\mathrm{St}-1}}.$
(33)

Thus

$\phi=\tan^{-1}\left(\frac{2\mathrm{St}u_{0}/x_{0}+1}{\sqrt{8\mathrm{St}-1}}\right).$
(34)

And rewriting using $\delta$

$x(t)=\psi\exp\left[-t/(2\mathrm{St})\right]\cos\left[\frac{8\sqrt{2}\delta^{1/2}}{1+8\delta}t+\phi\right].$
(35)

The collision occurs when $x(t_{COLL})=0$ and is set by the cosine term, and
so by the angular frequency, which scales as $\delta^{1/2}$. In the $\delta\to
0$ limit, $\phi$ in Eq. (34) tends to $-\pi/2$ (as $u_{0}<0$) and the
collision time is set by $8\sqrt{2}\delta^{1/2}t_{COLL}=\pi$, which leads to
Eq. (12).

##  Appendix D Solution of Eq. (17)

In order to derive the scaling argument on the coupling term in Eq. (17), we
consider only the exponentially growing term on the RHS from the solution Eq.
(14). Hence the differential equation becomes

$\ddot{x}+\frac{\dot{x}+2x}{\mathrm{St}}=C_{0}^{2}\mu_{1}^{2}\exp(2\mu_{1}t).$
(36)

The general solution to this equation consists of a complementary function
which satisfies the homogeneous equation Eq. (8) and a particular integral
which satisfies the in-homogeneous RHS.

The general solution is then

$\displaystyle x(t)$ $\displaystyle=$
$\displaystyle\exp\left[-t/(2\mathrm{St})\right]\left(A_{0}\cos\left[\frac{\sqrt{8\mathrm{St}-1}}{2\mathrm{St}}t\right]\right.$
(37)
$\displaystyle+\left.B_{0}\sin\left[\frac{\sqrt{8\mathrm{St}-1}}{2\mathrm{St}}t\right]\right)+E_{0}\exp(2\mu_{1}t),$

where

$E_{0}=\frac{\mu_{1}^{2}C_{0}^{2}}{4\mu_{1}^{2}+\frac{2\mu_{1}}{\mathrm{St}}+\frac{2}{\mathrm{St}}},$

and

$\displaystyle A_{0}$ $\displaystyle=$ $\displaystyle x_{0}-E_{0},$ (38)
$\displaystyle B_{0}$ $\displaystyle=$
$\displaystyle\frac{2\mathrm{St}\,u_{0}+x_{0}-E_{0}(1+4\mathrm{St}\,\mu_{1})}{\sqrt{8\mathrm{St}-1}}.$
(39)

##  Appendix E Inertia-driven collisions with a sphere

For a sphere in an inviscid axisymmetric flow, the flow field is[25]

$\frac{\bm{\mathbf{u}}_{\rm
sphere}}{U}=\left(1-\frac{R^{3}}{r^{3}}\right)\cos(\theta)\widehat{\bm{\mathbf{r}}}-\left(1+\frac{R^{3}}{2r^{3}}\right)\sin(\theta)\widehat{\bm{\mathbf{\theta}}},$
(40)

which when expanded about the forward stagnation point with $U=R=1$ gives

$\displaystyle\bm{\mathbf{u}}_{\rm sphere}$ $\displaystyle=$
$\displaystyle\left(-3x+6x^{2}+\mathcal{O}({\rm
cubic~{}terms})\right)\hat{\bm{\mathbf{r}}}$ (41)
$\displaystyle+\left(-\frac{3}{2}y+\frac{3}{2}xy+\mathcal{O}({\rm
cubic~{}terms})\right)\hat{\bm{\mathbf{\theta}}}.$

As Eq. (41) is the same as the equivalent for a cylinder, Eq. (6), apart from
numerical prefactors, there is also a critical value of the Stokes number for
a sphere. The different numerical factors shift it to the lower value of
$\mathrm{St}_{c}=1/12$, which agrees with the result of Langmuir and Blodgett
[16]. The change from a $-2x$ term in the flow field for a cylinder to a $-3x$
term for a sphere just shifts the critical value down by a factor of $2/3$.
The equivalent of Eq. (8) for a sphere just has a $-3x$ term in place of the
$-2x$ term. This means that we have the same $1/\delta^{1/2}$ scaling of the
time to collide.

In spherical polar coordinates, the acceleration terms along a constant
azimuthal angle trajectory, in Eq. (3) are the same as those for the two-
dimensional cylindrical polar coordinates. Hence the left-hand-sides of Eq.
(4), and thus of Eq. (7a) and Eq. (7b), are the same as presented here, i.e.
the $\dot{y}^{2}$ term is the significant term, and the analysis presented for
the cylinder can be repeated for the sphere. The result is that the scaling
law for the deposition efficiency has the same $\delta$ dependence as for the
cylinder except with modified multiplicative constants.

## References

  * Pruppacher [1978] H. R. Pruppacher,  _Microphysics of clouds and precipitation_ (D. Reidel Pub. Co, Dordrecht, Holland, 1978).
  * Wang and Otani [2013] C.-s. Wang and Y. Otani, Removal of nanoparticles from gas streams by fibrous filters: A review, Industrial & Engineering Chemistry Research 52, 5 (2013).
  * Robinson  _et al._ [2021] J. F. Robinson, I. Rios de Anda, F. J. Moore, J. P. Reid, R. P. Sear, and C. P. Royall, Efficacy of face coverings in reducing transmission of COVID-19: Calculations based on models of droplet capture, Physics of Fluids 33, 043112 (2021).
  * Robinson  _et al._ [2022a] J. F. Robinson, I. Rios de Anda, F. J. Moore, F. K. A. Gregson, J. P. Reid, L. Husain, R. P. Sear, and C. P. Royall, How effective are face coverings in reducing transmission of COVID-19?, Aerosol Science and Technology 56, 473 (2022a).
  * Parker and Lawrence [2001] A. R. Parker and C. R. Lawrence, Water capture by a desert beetle, Nature 414, 33 (2001).
  * Shahrokhian  _et al._ [2020] A. Shahrokhian, J. Feng, and H. King, Surface morphology enhances deposition efficiency in biomimetic, wind-driven fog collection, J. Roy. Soc. Interface 17, 20200038 (2020).
  * Azeem  _et al._ [2020] M. Azeem, A. Guérin, T. Dumais, L. Caminos, R. E. Goldstein, A. I. Pesci, J. de Dios Rivera, M. J. Torres, J. Wiener, J. L. Campos, and J. Dumais, Optimal design of multilayer fog collectors, ACS App. Mat. Int. 12, 7736 (2020).
  * Niklas [1985] K. J. Niklas, Wind pollination—a study in controlled chaos: Aerodynamic studies of wind-pollinated plants reveal a high degree of control in the apparently random process of pollen capture, American Scientist 73, 462 (1985).
  * Paw U and Hotton [1989] K. T. Paw U and C. Hotton, Optimum pollen and female receptor size for anemophily, Am J Botany , 445 (1989).
  * Makkonen [1984] L. Makkonen, Modeling of ice accretion on wires, Journal of Applied Meteorology and Climatology 23, 929 (1984).
  * Makkonen [2000] L. Makkonen, Models for the growth of rime, glaze, icicles and wet snow on structures, Philosophical Transactions: Mathematical, Physical and Engineering Sciences 358, 2913 (2000).
  * Gao  _et al._ [2021] L. Gao, T. Tao, Y. Liu, and H. Hu, A field study of ice accretion and its effects on the power production of utility-scale wind turbines, Renew. Energ. 167, 917 (2021).
  * Albrecht [1931] F. Albrecht, Theoretische Untersuchungen über die Ablagerung von Staub aus strömender Luftund ihre Anwendung auf die Theorie Staubfilter, Physik. Zeithschift 32, 48 (1931).
  * Taylor [1940] G. I. Taylor, Notes on possible equipment and technique for experiments on icing on aircraft, Reports and Memoranda of the Aeronautical Research Committee 2024 (1940).
  * Taylor [1963] G. I. Taylor,  _The scientific papers of Sir Geoffrey Ingram Taylor. Volume Three. Aerodynamics and Mechanics of Projectiles and Explosions_ , Vol. 3 (Cambridge University Press, 1963).
  * Langmuir and Blodgett [1946] I. Langmuir and K. B. Blodgett,  _A Mathematical Investigation of Water Droplet Trajectories_, AAF technical report (Army Air Forces Headquarters, Air Technical Service Command, 1946).
  * Finstad  _et al._ [1988a] K. J. Finstad, E. P. Lozowski, and E. M. Gates, A computational investigation of water droplet trajectories, Journal of Atmospheric and Oceanic Technology 5, 160 (1988a).
  * Finstad  _et al._ [1988b] K. J. Finstad, E. P. Lozowski, and L. Makkonen, On the median volume diameter approximation for droplet collision efficiency, Journal of Atmospheric Sciences 45, 4008 (1988b).
  * Phillips and Kaye [1999] C. Phillips and S. Kaye, The influence of the viscous boundary layer on the critical stokes number for particle impaction near a stagnation point, Journal of aerosol science 30, 709 (1999).
  * Araújo  _et al._ [2006] A. D. Araújo, J. S. Andrade, and H. J. Herrmann, Critical Role of Gravity in Filters, Phys. Rev. Lett. 97, 138001 (2006).
  * Wong  _et al._ [1955] J. B. Wong, W. E. Ranz, and H. F. Johnstone, Inertial impaction of aerosol particles on cylinders, Journal of Applied Physics 26, 244 (1955).
  * Makkonen and Stallabrass [1987] L. Makkonen and J. Stallabrass, Experiments on the cloud droplet collision efficiency of cylinders, Journal of Applied Meteorology and Climatology 26, 1406 (1987).
  * Makkonen [1992] L. Makkonen, Analysis of rotating multicylinder data in measuring cloud-droplet size and liquid water content, Journal of Atmospheric and Oceanic Technology 9, 258 (1992).
  * Makkonen  _et al._ [2018] L. Makkonen, J. Zhang, T. Karlsson, and M. Tiihonen, Modelling the growth of large rime ice accretions, Cold Regions Science and Technology 151, 133 (2018).
  * Acheson [1990] D. J. Acheson,  _Elementary fluid dynamics_ (Clarendon Press, Oxford, 1990).
  * Robinson  _et al._ [2022b] J. F. Robinson, P. B. Warren, M. R. Turner, and R. P. Sear, private communication (2022b).
  * Ingham  _et al._ [1990] D. Ingham, L. Hildyard, and M. Hildyard, On the critical stokes’ number for particle transport in potential and viscous flows near bluff bodies, Journal of Aerosol Science 21, 935 (1990).
  * Rios de Anda  _et al._ [2022] I. Rios de Anda, J. W. Wilkins, J. F. Robinson, C. P. Royall, and R. P. Sear, Modeling the filtration efficiency of a woven fabric: The role of multiple lengthscales, Physics of Fluids 34, 033301 (2022).


# On the Optimal Configuration of a Square Array Group Testing Algorithm

Ugnė Čižikovienė  Institute of Applied Mathematics, Faculty of Mathematics and
Informatics, Vilnius University, Naugarduko 24, Vilnius LT-03225, Lithuania
Viktor Skorniakov corresponding author; e-mail<EMAIL_ADDRESS>Institute of Applied Mathematics, Faculty of Mathematics and Informatics,
Vilnius University, Naugarduko 24, Vilnius LT-03225, Lithuania

###### Abstract

Up to date, only lower and upper bounds for the optimal configuration of a
Square Array (A2) Group Testing (GT) algorithm are known. We establish exact
analytical formulae and provide a couple of applications of our result. First,
we compare the A2 GT scheme to several other classical GT schemes in terms of
the gain per specimen attained at optimal configuration. Second, operating
under objective Bayesian framework with the loss designed to attain minimum at
optimal GT configuration, we suggest the preferred choice of the group size
under natural minimal assumptions: the prior information regarding the
prevalence suggests that grouping and application of A2 is better than
individual testing. The same suggestion is provided for the Minimax strategy.

##  1 Introduction

The task of identification of infected patients in a given cohort is the
frequent one. Though the plain consecutive testing of all individuals is an
obvious solution, there are many other ways to approach that problem. The term
Group Testing (GT) refers to the testing strategy when the testing of distinct
specimens is replaced by the testing of groups of pooled specimens. It appears
that the idea was first described in the famous paper of Dorfman [12]. He
looked for a cost saving way to screen the U.S. soldiers for syphilis during
the period of the World War II and suggested the following scheme. Instead of
testing each individual blood sample, pool $N$ samples and test the group; if
the group tests positive, retest each single individual; if the group tests
negative, then all individuals in the group are healthy and no retesting is
needed. It is clear that, when the prevalence of the disease is low, a small
fraction of pools needs retesting thereby leading to significant cost savings.

Since the appearance of [12], the idea came to stay to many other fields
(quality control, informational sciences, environmental sciences, etc.) as
well. In biomedical context, GT is widely applied to screen for infectious
diseases like HIV, hepatitis and, most recently, COVID-19 ([42], [40], [5],
[25], [35], [38], [1], [10], [11], [16], [22], [30], [27]). It also appears to
be a very useful technique in genetics ([13], [9], [28], [7]).

In this paper, we focus on the Square Array (A2) GT algorithm introduced by
Phatarfod and Sudbury [37] and later generalized by Berger, Mandell and
Subrahmanya [36]. The A2 operates as follows. Given $n^{2}$ specimens, one
places them on $n\times n$ matrix and tests the pools defined by subsets
corresponding to rows and columns. The cases lying on the intersections of
rows and columns exhibiting positive responses are further retested and,
assuming that the test is perfect, all infected are identified.

One of the most important characteristics of each GT algorithm is the _optimal
configuration_. To introduce the concept, consider an arbitrary GT scheme and
let $T_{N}$ denote the total (random) number of tests performed over the
cohort spanning $N$ individuals. The optimal configuration is the size of the
cohort which minimizes function $N\mapsto\frac{\operatorname{E}T_{N}}{N}$,
i.e., the expected number of tests per individual. Though, for a given
prevalence, numerical solution of the optimal configuration is always
possible, analytical formulae provide much more insights allowing, in
particular, analytical comparisons of different GT schemes. For the case of
the Dorfman scheme described above, the optimal configuration was established
by Samuels [31] quite long ago. However, for a couple of its modifications,
namely, the modified Dorfman scheme [33] and Sterrett scheme [34], the
analytical formulae, though conjectured, were unknown [24] and established
only recently [44]. Pretty much the same situation is with A2. To our best
knowledge, only Hudgens and Kim [17] addressed this problem and did not
succeed in providing the final solution. Namely, in [17] the authors obtained
lower and upper bounds for the optimal A2 configuration leaving the exact
formulae undiscovered. In this paper, we fill in the gap by providing exact
solution similar to that obtained by Samuels [31] and Skorniakov and
Čižikovienė [44] for the case of classical Dorfman scheme and its
modifications.

The paper is organized as follows. In Section 2, we state our results and
provide a couple of applications. Section 3 is devoted to discussion. Finally,
there are two appendices containing proofs, figures, and link to the web
repository with tabular data.

##  2 Results and Applications

###  2.1 Statement of the Results

Before proceeding to the statement of results, we first formulate Binomial
Testing Assumptions (BTA) which are assumed to hold in the remaining part of
the paper by default.

  * (BTA1)

The tested cohort consists of independent individuals. Each individual is
infected with the same constant probability $p\in(0,1)$ (termed prevalence in
the sequel).

  * (BTA2)

The test under consideration is perfect and there is no dilution effect. That
is, pooling does not affect the performance of the test.

The way A2 operates was already described in the introductory Section 1.
However, we need an explicit expression for the average number of tests
$\operatorname{E}T_{N}$ applied under A2 to the cohort spanning $N=n^{2},n\geq
2$, individuals. The latter was derived by Phatarfod and Sudbury [37] and is
equal to

$2n+n^{2}\left(1-2q^{n}+q^{2n-1}\right),$ (2.1)

where $q:=1-p$. Therefore, in this parametrization, an average number of tests
per person

$t(q,n):=\frac{1}{N}\operatorname{E}T_{N}\Big{|}_{N=n^{2}}=\frac{2}{n}+1-2q^{n}+q^{2n-1},$
(2.2)

and the corresponding optimal configuration

$n_{opt}=n_{opt}(q):=\mathrm{argmin}_{n\in\\{2,3,\ldots\\}}t(q,n).$

The existence, uniqueness and bounds on $n_{opt}$ for various values of
$q\in(0,1)$ were established by Hudgens and Kim [17]. Our reconsideration
(given in Theorem 2.1 and Corollary 2.1 below) aimed to sharpen their results
(see discussion in Section 3) and provide the complete theoretical
characterization of A2.

######  Theorem. 2.1.

Let $g(q,n)=\frac{2}{n}-2q^{n}+q^{2n-1}=t(q,n)-1$.

  * (i)

For $(q,n)$ ranging in $(1/2,1)\times(2,\infty)$, system of equations

$\begin{cases}&\ln q=nq^{n}\left(1-\frac{q^{n-1}}{2}\right)\ln q\\\ &n\ln
q=-\frac{\left(1-\frac{q^{n-1}}{2}\right)}{\left(1-{q^{n-1}}\right)}\end{cases}$
(2.3)

has a unique solution $(q_{*},n_{*})\approx(0.748416,4.453524)$.

  * (ii)

For any fixed $q\in(q_{*},1)$ and with respect to $n$, equation $g(q,n)=0$
admits two solutions $n_{L},n_{U}:2<n_{L}<n_{*}<n_{U}<\infty$. On
$(n_{L},n_{U})$, $n\mapsto g(q,n)$ attains values in $(-\infty,0)$ whereas on
$(2,\infty)\setminus[n_{L},n_{U}]$ it attains values in $(0,\infty)$.

  * (iii)

For any fixed $q\in(q_{*},1)$, the region $(n_{L},n_{U})$ is the one where A2
is efficient, i.e., $t(q,n)<1$ for $n\in(n_{L},n_{U})$. In that region, there
exists a unique (and, therefore, global) minimizer $n_{min}$ of $(2,\infty)\ni
n\mapsto t(q,n)$. For $q\in[0.755,1)$, it is given by

$n_{min}=\frac{1}{p^{\frac{2}{3}}}+\frac{1}{2p^{\frac{1}{3}}}+0.2+3p^{2}+t_{*}$
(2.4)

for some $t_{*}\in[0,1]$.

$(2,\infty)\ni n\mapsto t(q,n)$ also has a unique (and, therefore, global)
maximizer located in the region $(n_{U},\infty)$. For any fixed
$q\in(0,q_{*})$, A2 is never optimal, i.e., $(2,\infty)\ni n\mapsto t(q,n)$
attains values in $(1,\infty)$.

######  Corollary. 2.1.

Let $g(q,n)$ be as in Theorem (2.1). Then $g(q,5)=0$ has a unique solution
$q_{5}\approx 0.750209961$. For all $q\in(q_{5},1)$, $n_{opt}(q)$ belongs to
the set

$\left\\{\left\lfloor\frac{1}{p^{\frac{2}{3}}}+\frac{1}{2p^{\frac{1}{3}}}+3p^{2}+0.2\right\rfloor+i:i=0,1,2\right\\}.$
(2.5)

######  Remark 2.1.

In the Corollary 2.1, the region $(q_{5},1)\subset(q_{*},1)$. An explanation
for this truncation stems from the fact that $(q_{5},1)$ is the region where
practical application of A2 makes sense. To be more precise, applying Theorem
2.1 for a fixed $q\in(q_{*},q_{5})$, we have that $t(q,n_{min})<1$ with
$n_{min}=n_{min}(q)$ given by (2.4). However, this $n_{min}(q)\in(4,5)$ and
$\min(t(q,4),t(q,5))>1$ when $q\in(q_{5},q_{*})$. We touch this question
briefly in the discussion Section 3 when talking about relation of our results
to those of Hudgens and Kim [17]. Though we do not provide a separate proof of
this fact, technical details can be filled in after inspection of the proofs
presented in the Appendix A.

###  2.2 Examples of Applications

####  2.2.1 Comparison to other GT schemes

In order to shed the light on to performance of A2, we compare it to several
other GT schemes in terms of magnitude of optimal configuration $N_{opt}$ and
gain across the range $p\in(0,0.249790)$ where application of A2 seems
reasonable. In this comparative analysis, for a fixed $q=1-p$, we define the
gain as

$G(q)=(1-t(q,N_{*})).$

Here, $t(q,N)=\frac{\operatorname{E}T_{N}}{N}$ is an average number of tests
per person when the tested group size is equal to $N$ whereas $N_{*}$ stands
for the unrounded optimal configuration, i.e., the minimizer of the continuous
argument function $[1,\infty)\ni N\mapsto t(q,N)$. Defined this way, the gain
multiplied by 100 has a meaning of an average number of tests saved per 100
persons in comparison to usual one by one testing. Preceded by several
remarks, below comes a description of the a fore mentioned GT schemes. To
distinguish between the schemes, we assign one letter abbreviations to each
and, by making use of these letters, superscript all related quantities. E.g.,
$t^{(A2)}(q,N_{*}^{(A2)})$ denotes the value of $t(q,N_{*})$ when A2 is the
scheme under consideration.

######  Remark 2.2.

Recall that, in section 2.1, we have introduced reparametrization of $N$ by
equality $N=n^{2}$, where $n$ is the number of rows (columns) in the squared
array used in the definition of A2. In this example and in all what follows
after, this reparametrization remains in force: writing $t^{(A2)}(q,N)$, we
factually mean $t$ given in equation (2.4) with $n$ ranging continuously
unless stated otherwise. When we want to emphasize reference to (2.4), $n$ is
used instead of $N$. For all other schemes considered, reparametrizations of
similar kind do not apply.

######  Remark 2.3.

The behaviour of quantities compared is more naturally interpreted in terms of
the prevalence $p=1-q$. Therefore, $p$ but not $q$ is used in the accompanying
graphs. Also, because of the same reason and for the sake of convenience, we
quite often denote an argument of the function considered by $q$ and write an
explicit formula in terms of $p$.

######  Remark 2.4.

In this example, we have chosen to operate on the continuous scale because it
is much easier to perceive visually in comparison to the discrete one.
However, keeping in a view the practical aspect, comprehensive numerical
results (see Appendix B) are given on the discrete scale. Comparing both one
can find out that the discrepancies are small.

_Dorfman scheme D_. $\rhd$ The Dorfman scheme was described in the
introductory Section 1. For this scheme,

$t^{(D)}(q,N)=\frac{1}{N}+1-q^{N}$ (2.6)

and $N_{*}^{(D)}$ solves equation $1/N^{2}=-q^{N}\ln q$. Samuels [31] has
shown that rounded optimal configuration $N_{opt}^{(D)}$ is either $1+\lfloor
p^{-1/2}\rfloor$ or $2+\lfloor p^{-1/2}\rfloor$. Hence,
$N_{*}^{(D)}=\sqrt{p^{-1}}(1+o(1)),\ p\to 0+$. Approximation
$N_{*}^{(D)}\approx\sqrt{p^{-1}}$ is accurate enough (see [32] for tabulated
numerical results), and one can further show that

$t^{(D)}(q,N_{*}^{(D)})=2\sqrt{p}(1+o(1)),\ p\to 0+.$

Consequently, $G^{(D)}(q)=(1-2\sqrt{p})(1+o(1)),\ p\to 0+$. $\lhd$

_Sterrett scheme S_. $\rhd$ Sterrett [34] suggested the following modification
of the Dorfman111in fact, his suggestion was built on the already modified
Dorfman scheme scheme: one should retest initial positive pool sequentially
one by one until appearance of the first positive case and then again apply
pool testing to the remaining tail. If the remaining untested set tests
positive, one should proceed recursively as previously until the remaining set
tests negative or the whole set of individuals gets tested. Malinovsky and
Albert [24] conjectured that, for $p\in(0,(3-\sqrt{5})/2)$, rounded optimal
configuration $N_{opt}^{(S)}$ lies in the set

$\left\\{\lfloor\sqrt{2p^{-1}}\rfloor,\lfloor\sqrt{2p^{-1}}\rfloor+1,\lfloor\sqrt{2p^{-1}}\rfloor+2\right\\}.$

Skorniakov and Čižikovienė [44] affirmed their conjecture showing along the
way that $N_{*}^{(S)}\in\left[\sqrt{2p^{-1}}-1,\sqrt{2p^{-1}}+1\right]$ for
$p\in(0,(3-\sqrt{5})/2)$. From the latter result it then follows that

$\displaystyle N_{*}^{(S)}=\sqrt{2p^{-1}}(1+o(1)),\quad
t^{(S)}(q,N_{*}^{(S)})=\sqrt{2p}(1+o(1)),\ p\to 0+.\ $

Hence, $G^{(S)}(q)=(1-\sqrt{2p})(1+o(1)),\ p\to 0+$. $\lhd$

_Halving scheme H_. $\rhd$ This scheme resembles divide and conquer sorting
and can be described by the following steps.

  * Step 1.

Test initial pooled cohort. If it tests negative, finish; if it tests
positive, proceed to Step 2.

  * Step 2.

Divide the cohort into two approximately equal parts consisting of the first
and second halves and apply the whole algorithm (starting from Step 1) to the
two obtained parts recursively.

It is difficult to trace back the first reference discussing this scheme in
detail. To our best knowledge, it was treated already by Johnson et. al. [19].
However, it appears that its asymptotic analysis was first accomplished not so
long ago by Zamman and Pippenger [43] whereas in [32] it was discussed a fresh
without a strict focus on asymptotic regime when $p\to 0+$. There it was shown
that rounded optimal configuration

$N_{opt}^{(H)}\in\\{\lfloor 1/(2\log_{2}(1/q))\rfloor,\lfloor
1/(2\log_{2}(1/q))\rfloor+1\\}$

leading thereby to the following asymptotic relationships:

$\displaystyle N_{*}^{(H)}=-\frac{1}{2\log_{2}q}(1+o(1)),\quad
t^{(H)}(q,N_{*}^{(H)})=-({2p}\log_{2}p)(1+o(1)),$ $\displaystyle
G^{(H)}(q)=(1+{2p}\log_{2}p)(1+o(1)),\ p\to 0+.\ \lhd$

Figure 1 shows the behaviour of $N_{*}$ and gain $G$ for the case of A2 and
the three schemes described above. Note that, in case of A2, total optimal
pool size $N_{*}^{(A2)}=n_{min}^{2}$ with $n_{min}$ given in Theorem 2.1. Due
to that raise to the square, $N_{*}^{(A2)}$ grows to infinity much faster than
the counterparts of the remaining schemes, and, because of this, in the top
left sub-figure, the range of $p$ starts quite far from the origin and the
accompanying bottom left sub-figure on the log scale is given. It clearly
depicts the relationships

$\ln\left(N_{*}^{(D)}\right),\ln\left(N_{*}^{(S)}\right)\sim-\frac{1}{2}\ln(p),\ln\left(N_{*}^{(H)}\right)\sim-\ln(p),\ln\left(N_{*}^{(A2)}\right)\sim-\frac{4}{3}\ln(p),$
(2.7)

following from the formulae given above and clearly showing that the
asymptotic slope of $\ln\left(N_{*}^{(A2)}\right)$ on the $-\ln p$ scale is
the largest one.

Turning to sub-figures on the right, one sees that there are ranges where A2
outperforms other schemes. Numerical solutions are as follows:

  * •

$R_{A2,D}=\\{p\in(0,0.249790):G^{(A2)}(q)>G^{(D)}(q)\\}=(0,0.115589)$;

  * •

$R_{A2,S}=\\{p\in(0,0.249790):G^{(A2)}(q)>G^{(S)}(q)\\}=(0,0.028071)$;

  * •

$R_{A2,H}=\\{p\in(0,0.249790):G^{(A2)}(q)>G^{(H)}(q)\\}=(0.012936,0.220788)$.

Within these ranges, numerically estimated maximal differences are

  * •

$\max_{p\in R_{A2,D}}100(G^{(A2)}(q)-G^{(D)}(q))=6.2179$ at $p=0.017128$;

  * •

$\max_{p\in R_{A2,S}}100(G^{(A2)}(q)-G^{(S)}(q))=1.9342$ at $p=0.003984$;

  * •

$\max_{p\in R_{A2,H}}100(G^{(A2)}(q)-G^{(H)}(q))=6.5951$ at $p=0.104908$.

Finally, comparing behaviour of optimal pool sizes one has that
$\sqrt{N_{*}^{(A2)}}$ (and, therefore, $N_{*}^{(A2)}$ as well) exceeds
$\max(N_{*}^{(D)},N_{*}^{(S)})$ for all $p\in(0,0.249790)$. For scheme H,
however, the following holds true:

$\displaystyle\sqrt{N_{*}^{(A2)}}<N_{*}^{(H)}\quad\text{for}\quad
p\in(0,0.023178];$
$\displaystyle\sqrt{N_{*}^{(A2)}}>N_{*}^{(H)}\quad\text{for}\quad
p\in[0.023179,0.249790);$
$\displaystyle{N_{*}^{(A2)}}>N_{*}^{(H)}\quad\text{for all}\quad
p\in(0,0.249790).$

Figure 2 provides a zoomed in visual illustration. Regarding the role and
distinction between $\sqrt{N_{*}^{(A2)}}$ and ${N_{*}^{(A2)}}$, consult the
discussion Section 3.

####  2.2.2 Optimal configuration when the prevalence is unknown

In reference [23], the authors sought for the pool size leading to optimal
testing by making use of scheme D when the prevalence is unknown. In their
work, two approaches were used. Both (approaches) were based on the following
loss function. Given scheme $X$, define

$L^{(X)}(q,N)=t^{(X)}(q,N)-t^{(X)}(q,N_{opt}^{(X)}),$ (2.8)

where
$N^{(X)}_{opt}=N^{(X)}_{opt}(q)=\operatorname*{argmin}_{N\in\\{1,2,\ldots\\}}t^{(X)}(q,N)$
is the optimal configuration when the prevalence $p$ (and hence $q=1-p$) is
known. It is clear that $L^{(X)}(q,N)\geq 0\
\forall(q,N)\in(0,1)\times\mathbb{N}$ and, for a given $q$, $L^{(X)}(q,N)=0$
precisely when $N=N_{opt}^{(X)}(q)$.

In what follows, to distinguish between $N_{opt}^{(X)}(q)$ and optimal
configuration suitable for unknown $q$’s, the latter configuration is denoted
by $N_{\star}^{(X)}$.

The first approach in [23] was to make use of mini–max strategy and take
$N_{\star}^{(X)}$ as a minimizer of

$\\{1,2,\ldots\\}\ni N\mapsto\sup_{q}L^{(X)}(q,N).$ (2.9)

The second approach was to make use of Bayesian paradigm and, after putting
the prior $\pi$ on $q$, to take $N_{\star}^{(X)}$ as a minimizer of

$\\{1,2,\ldots\\}\ni
N\mapsto\operatorname{E}_{\pi}L^{(X)}(q,N)=\operatorname{E}_{\pi}t^{(X)}(q,N)-c(\pi),\quad\\\
c(\pi)=\operatorname{E}_{\pi}t^{(X)}(q,N_{opt}^{(X)}).$ (2.10)

In this example, we have adopted both approaches to the case of A2. When using
the Bayesian one, $\pi$ was taken uniform over $(0.750210,1)$. Thereby, we
have modeled situation when the only prior information is that application of
A2 makes sense (see Corollary 2.1). Also, we have modified (2.10) and used

$\\{1,2,\ldots\\}\ni
N\mapsto\operatorname{E}_{\pi}(L^{(X)}(q,N))^{2}=\operatorname{E}_{\pi}\left(t^{(X)}(q,N)-t^{(X)}(q,N_{opt}^{(X)})\right)^{2}$
(2.11)

instead. To justify our choice, note that, in (2.10), $c(\pi)$ does not depend
on $N$. Therefore, minimization of the target function amounts to minimization
of $\\{1,2,\ldots\\}\ni N\mapsto\operatorname{E}_{\pi}t^{(X)}(q,N)$. This way
important information carrying function $q\mapsto t^{(X)}(q,N_{opt}^{(X)})$
remains unutilized.

Figures 3–4 show graphs of (2.9) and (2.11) for the case of $X=\mathrm{A2}$
and the previously mentioned prior $\pi$. Note that, adopting the above to our
case, in (2.8), (2.9), (2.10), we have used function $t$ defined by (2.2).
Numerical estimation yielded the following values:

  * •

$N_{\star}^{(A2)}=12^{2}$ for the case of mini–max approach;

  * •

$N_{\star}^{(A2)}=7^{2}$ for the case of Bayesian approach.

Finishing, it is important to note that, though the strategy discussed above
leads to suboptimal testing in a stable environment where reliable estimation
of prevalence is possible, it appears to be a reasonable strategy when the
prevalence is varying rapidly and is difficult to capture by data at hand.
Therefore, at least in the initial stage, it can be considered as a good
alternative for optimal testing during pandemics like COVID–19.

##  3 Discussion

It was already mentioned in the introductory Section 1 that, to our best
knowledge, [17] is the only reference where A2 was treated in the same way
like we did here. Therefore, we first discuss our input in comparison with
[17] and then turn to a more general setting.

###  3.1 Comparison with the previous work

Figure 5 illustrates the behaviour of the function $n\mapsto q_{n}$ appearing
in the proof of Theorem 2.1. The minimum of this function, denoted by $q_{*}$
is an exact lower bound of the region where A2 is efficient _on a continuous
scale_ in a sense that, given prevalence $p\in(0,1-q_{*})$ (or, alternatively,
$q\in(q_{*},1)$), the achievable minimum of $(2,\infty)\ni n\mapsto
t^{(A2)}(q,n)$ is strictly smaller than $1$. This always holds true for the
unrounded optimal configuration, i.e., the minimizer of $(2,\infty)\ni
n\mapsto t^{(A2)}(q,n)$ on the continuous scale. In [17], the authors also
provide a region of this kind. We utilize that region in Corollary 2.1 and
denote it $(q_{5},1)$. Remark 2.1 explains why it can be thought of as a
region where A2 makes sense from the practical point of view and how this can
be derived from our results. Linking our work to theirs, it is important to
mention that, deriving their proofs, Hudgens and Kim [17] operated on the
discrete scale and obtained this region as the one where $\\{5,6,\ldots\\}\ni
n\mapsto t^{(A2)}(q,n)$ is efficient. Having proved that $n=2,3,4$ are never
optimal, they have also verified that $(q_{5},1)$ is the region where
practical application of A2 makes sense.

Summing up, in this direction, our input adds the missing part to the
theoretical characterization of A2 yet does not bring novelty from the applied
point of view. Turning to formulae (2.4), (2.5), the situation is essentially
different. Kim and Hudgens [17] have given lower and upper bounds on the
configuration $\sqrt{N_{opt}^{(A2)}}=n_{opt}^{(A2)}\in\\{5,6,\ldots\\}$ and
ascertained that the optimal configuration is unique. These bounds facilitate
numerical computations and allow some insights on the analytical scale yet
they do not accommodate all the benefits brought by formulae (2.4)–(2.5):

  1. 1.

(2.5) provides _a ready to use_ expressions for $N_{opt}^{(A2)}$ which are
practically important when it comes to $p$’s ranging near the origin where
numerical solutions might cause computational problems;

  2. 2.

(2.4) and (2.5) enable to obtain analytical comparisons of A2 with other
schemes. An example of such analysis is given in Subsection 2.2.1 when
deriving relationships (2.7). Without knowing analytical expressions of
$N_{opt}^{(A2)}$ it would be impossible to contrast the asymptotic slope of
$-\ln p\mapsto\ln N_{*}^{(A2)}(p)$ to the analogous asymptotic slopes of other
schemes considered there. One can go even further. Namely, if $X$ is another
scheme of interest, one can compare analytically the behaviour of
$t^{(A2)}(q,N_{opt}^{(A2)})$ with $t^{(X)}(p,N_{opt}^{(X)})$ as $p\to 0+$ (see
a comment in the forthcoming Subsection 3.2 regarding an importance of
analysis of this kind). It is again worth to stress up that, in this
asymptotic setting, validity of pure numerical analysis is always questionable
whereas analytical result leads to accurate and definite analysis.

###  3.2 Discussion of other aspects

Several authors have already noted that array based algorithms can be more
efficient in certain settings [36], [21], [20]. Our findings confirm these
observations: comparisons given in Subsection 2.2.1 demonstrated that there
were regions where A2 performed better than other considered schemes. The same
behaviour is expected for other unconsidered schemes satisfying binomial
testing assumptions. Therefore, though initially applied in a frame of genetic
a screening [4], [3], [6], A2 is an appropriate candidate for other GT
applications as well. Since GT is most effective when the prevalence $p$ is
low, for any scheme $X$, asymptotic behaviour of $t^{(X)}(p,N_{opt}^{(X)})$ as
$p\to 0+$, is an important characteristic. Among the schemes contrasted in
Subsection 2.2.1, A2 takes an intermediate position (with respect to this
characteristic) because222the relationships are not difficult to justify by
making use of the formulae (2.5)

$\displaystyle\lim\limits_{p\to
0+}\frac{t^{(A2)}(q,N_{opt}^{(A2)})}{t^{(X)}(q,N_{opt}^{(X)})}=0\quad\text{for}\quad
X\in\\{D,S\\},$ $\displaystyle\text{and}\quad\lim\limits_{p\to
0+}\frac{t^{(H)}(q,N_{opt}^{(H)})}{t^{(A2)}(q,N_{opt}^{(A2)})}=0.$

This, however, comes at a cost of a quickly increasing _maximal tested pool
size_ at optimal configuration with relationships inverse to those above:

$\displaystyle\lim\limits_{p\to
0+}\frac{N_{opt}^{(X)}(q)}{\sqrt{N^{(A2)}_{opt}(q)}}=0\quad\text{for}\quad
X\in\\{D,S\\},$ (3.1) $\displaystyle\text{and}\quad\lim\limits_{p\to
0+}\frac{\sqrt{N_{opt}^{(A2)}(q)}}{N_{opt}^{(H)}(q)}=0.$ (3.2)

We would like to emphasize that, in case of A2,
$\sqrt{N_{opt}^{(A2)}}=n_{opt}^{(A2)}$ (not $N_{opt}^{(A2)}$) is the maximal
tested pool size (corresponding to the test of a distinct row/column) whereas
in case of all the rest schemes it is $N_{opt}^{(X)},X\in\\{D,H,S\\}$. Thus
the appearance of $\sqrt{N_{opt}^{(A2)}}$ in formulae (3.1)–(3.2) above.
$N_{opt}^{(A2)}$ is the batch size needed to organize optimal testing since
$N_{opt}^{(A2)}=n_{opt}^{(A2)}$ samples never get tested as a single specimen.
This distinction is important since in the applied setting one has to take
into account the fact that the operating characteristics of the test might
depend on the size of the tested pool. That is, sensitivity and/or specificity
of the test might become unacceptably low when the pooled sample is formed out
of $N\geq N_{critical}$ individual samples. Also, subjects might be dependent
and it might be reasonable to take into account imperfectness of the test even
on the single individual basis. These facts limit usage of many schemes
satisfying binomial testing assumptions. E.g., scheme H, also known under the
name of optimal testing algorithm and giving best asymptotic performance as
$p\to 0+$, requires large pool sizes. Inspecting figure 2, one sees that A2
also exhibits such behaviour. Nonetheless, investigations giving theoretical
characterization of the GT scheme considered are important because of the
following reasons.

  1. Convenience.

When choosing between two realizable schemes, the preference does not always
fall on the optimal one. If it turns out that the optimal scheme yields a
small surplus in gain, it can be exchanged to a more simply realizable
competitor. E.g., scheme D was recently employed in Lithuania [2] for massive
COVID-19 testing and a large number of other countries do the same [8].
Turning to A2, it appears that one of the reasons of its emergence in genetic
applications was an operational convenience.

  2. Tolerable errors.

Though there exist a lot of generalizations allowing imperfectness of the test
(e.g., [21], [20], [14]), as noted by several authors [36], [17], [24], the
schemes assuming perfect tests can be quite accurate since modern tests
exhibit very small errors.

  3. Benchmarking.

BTA based schemes, being more simple to treat analitically, provide
theoretically justified benchmark thresholds for more elaborated schemes which
assume imperfectness of the test and/or other specific conditions (e.g.,
testing outcome dependence on subject specific characteristics).

  4. Basement.

More elaborated schemes may emerge on the basis of the simpler ones. An
example is the scheme S introduced in Subsection 2.2.1. It was built on the
top of scheme D. In [21] and [14], A2 serves as a basis for extensions
incorporating imperfectness of the test and dilution effect.

## References

  * [1] Baha Abdalhamid, Christopher R Bilder, Emily L McCutchen, Steven H Hinrichs, Scott A Koepsell, and Peter C Iwen.  Assessment of Specimen Pooling to Conserve SARS CoV-2 Testing Resources.  American Journal of Clinical Pathology, 153(6):715–718, May 2020\.
  * [2] Erika Alonderytė.  Lithuanian firms to get 30 million euros to test employees, February 2021\.
  * [3] C.T. Amemiya, M.J. Alegria-Hartman, C. Aslanidis, C. Chen, J. Nikolic, J.C. Gingrich, and P.J. de Jong.  A two-dimensional YAC pooling strategy for library screening via STS and Alu -PCR methods.  Nucleic Acids Research, 20(10):2559–2563, 1992.
  * [4] Emmanuel Barillot, Bruno Lacroix, and Daniel Cohen.  Theoretical analysis of library screening using a N-dimensional pooling strategy.  Nucleic Acids Research, 19(22):6241–6247, 1991.
  * [5] Christopher R. Bilder, Joshua M. Tebbs, and Peng Chen.  Informative Retesting.  Journal of the American Statistical Association, 105(491):942–955, September 2010.
  * [6] W.J. Bruno, E. Knill, D.J. Balding, D.C. Bruce, N.A. Doggett, W.W. Sawhill, R.L. Stallings, C.C. Whittaker, and D.C. Torney.  Efficient pooling designs for library screening.  Genomics, 26(1):21–30, March 1995.
  * [7] Chang Chang Cao and Xiao Sun.  Combinatorial pooled sequencing: experiment design and decoding.  Quantitative Biology, 4(1):36, 2016.
  * [8] Wikipedia contributors.  COVID-19 testing.  https://en.wikipedia.org/wiki/COVID-19_testing.
  * [9] David J. Cutler and Jeffrey D. Jensen.  To pool, or not to pool?  Genetics, 186(1):41–43, 2010.
  * [10] Timo de Wolff, Dirk Pflüger, Michael Rehme, Janin Heuer, and Martin-Immanuel Bittner.  Evaluation of pool-based testing approaches to enable population-wide screening for COVID-19.  PLOS ONE, 15(12):e0243692, December 2020.
  * [11] Andreas Deckert, Till Bärnighausen, and Nicholas NA Kyei.  Simulation of pooled-sample analysis strategies for COVID-19 mass testing.  Bulletin of the World Health Organization, 98(9):590–598, September 2020.
  * [12] R. Dorfman.  The detection of defective members of large populations.  The Annals of Mathematical Statistics, 14(4):436–440, 1943.
  * [13] Dingzhu Du and Frank Hwang.  Pooling designs and nonadaptive group testing: important tools for DNA sequencing.  Number v. 18 in Series on applied mathematics. World Scientific, New Jersey, 2006.  OCLC: ocm70407992.
  * [14] Y.G. Habtesllassie, Linda M. Haines, H.G. Mwambi, and J.W. Odhiambo.  Array-based schemes for group screening with test errors which incorporate a concentration effect.  Journal of Statistical Planning and Inference, 167:41–57, December 2015.
  * [15] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [16] Catherine A. Hogan, Malaya K. Sahoo, and Benjamin A. Pinsky.  Sample Pooling as a Strategy to Detect Community Transmission of SARS-CoV-2.  JAMA, 323(19):1967, May 2020.
  * [17] Michael G. Hudgens and Hae-Young Kim.  Optimal Configuration of a Square Array Group Testing Algorithm.  Communications in Statistics - Theory and Methods, 40(3):436–448, January 2011.
  * [18] J. D. Hunter.  Matplotlib: A 2d graphics environment.  Computing in Science & Engineering, 9(3):90–95, 2007.
  * [19] Norman L. Johnson, Samuel Kotz, and Xizhi Wu.  Inspection Errors for Attributes in Quality Control.  Springer US, Boston, MA, 1991.
  * [20] Hae-Young Kim and Michael G. Hudgens.  Three-Dimensional Array-Based Group Testing Algorithms.  Biometrics, 65(3):903–910, September 2009.
  * [21] Hae-Young Kim, Michael G. Hudgens, Jonathan M. Dreyfuss, Daniel J. Westreich, and Christopher D. Pilcher.  Comparison of Group Testing Algorithms for Case Identification in the Presence of Test Error.  Biometrics, 63(4):1152–1163, December 2007.
  * [22] Stefan Lohse, Thorsten Pfuhl, Barbara Berkó-Göttel, Jürgen Rissland, Tobias Geißler, Barbara Gärtner, Sören L Becker, Sophie Schneitler, and Sigrun Smola.  Pooling of samples for testing for SARS-CoV-2 in asymptomatic people.  The Lancet Infectious Diseases, 20(11):1231–1232, November 2020\.
  * [23] Yaakov Malinovsky and Paul S. Albert.  A note on the minimax solution for the two-stage group testing problem.  The American Statistician, 69(1):45–52, 2015.
  * [24] Yaakov Malinovsky and Paul S. Albert.  Revisiting Nested Group Testing Procedures: New Results, Comparisons, and Robustness.  The American Statistician, 73(2):117–125, April 2019.
  * [25] S. May, A. Gamst, R. Haubrich, C. Benson, and D. M. Smith.  Pooled nucleic acid testing to identify antiretroviral treatment failure during HIV infection.  JAIDS Journal of Acquired Immune Deficiency Syndromes, 53(2):194–201, 2010.
  * [26] Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondřej Čertík, Sergey B. Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K. Moore, Sartaj Singh, Thilina Rathnayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta, Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, Štěpán Roučka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and Anthony Scopatz.  Sympy: symbolic computing in python.  PeerJ Computer Science, 3:e103, January 2017.
  * [27] Leon Mutesa, Pacifique Ndishimye, Yvan Butera, Jacob Souopgui, Annette Uwineza, Robert Rutayisire, Ella Larissa Ndoricimpaye, Emile Musoni, Nadine Rujeni, Thierry Nyatanyi, Edouard Ntagwabira, Muhammed Semakula, Clarisse Musanabaganwa, Daniel Nyamwasa, Maurice Ndashimye, Eva Ujeneza, Ivan Emile Mwikarago, Claude Mambo Muvunyi, Jean Baptiste Mazarati, Sabin Nsanzimana, Neil Turok, and Wilfred Ndifon.  A pooled testing strategy for identifying SARS-CoV-2 at low prevalence.  Nature, 589(7841):276–280, January 2021.
  * [28] Amir Najafi, Damoun Nashta-ali, Seyed Abolfazl Motahari, Mehrdad Khani, Babak H. Khalaj, and Hamid R. Rabiee.  Fundamental limits of pooled-DNA sequencing.  https://arxiv.org/abs/1604.04735, 2016.
  * [29] The pandas development team.  pandas-dev/pandas: Pandas, February 2020.
  * [30] Christopher D Pilcher, Daniel Westreich, and Michael G Hudgens.  Group Testing for Severe Acute Respiratory Syndrome– Coronavirus 2 to Enable Rapid Scale-up of Testing and Real-Time Surveillance of Incidence.  The Journal of Infectious Diseases, 222(6):903–909, August 2020\.
  * [31] S. M. Samuels.  The Exact Solution to the Two-Stage Group-Testing Problem.  Technometrics, 20(4):497–500, November 1978.
  * [32] Viktor Skorniakov, Remigijus Leipus, Gediminas Juzeliūnas, and Kęstutis Staliūnas.  Group testing: Revisiting the ideas.  Nonlinear Analysis: Modelling and Control, 26(3):534–549, May 2021\.
  * [33] M. Sobel and P. A. Groll.  Group testing to eliminate efficiently all defectives in a binomial sample.  Bell System Technical Journal, 38:1179–1252, 1959.
  * [34] Andrew Sterrett.  On the Detection of Defective Members of Large Populations.  The Annals of Mathematical Statistics, 28(4):1033–1036, December 1957.
  * [35] Susan L. Stramer, Ulrike Wend, Daniel Candotti, Gregory A. Foster, F. Blaine Hollinger, Roger Y. Dodd, Jean-Pierre Allain, and Wolfram Gerlich.  Nucleic Acid Testing to Detect HBV Infection in Blood Donors.  New England Journal of Medicine, 364(3):236–247, January 2011.
  * [36] Toby Berger; James W. Mandell; P. Subrahmanya.  Maximally efficient two-stage screening.  Biometrics, 56, 2000.
  * [37] R. M. Phatarfod; Aidan Sudbury.  The use of a square array scheme in blood testing.  Statistics in Medicine, 13, 1994.
  * [38] Joshua M. Tebbs, Christopher S. McMahan, and Christopher R. Bilder.  Two-Stage Hierarchical Group Testing for Multiple Infections with Application to the Infertility Prevention Project.  Biometrics, 69(4):1064–1073, December 2013.
  * [39] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [40] Lawrence M. Wein and Stefanos A. Zenios.  Pooled Testing for HIV Screening: Capturing the Dilution Effect.  Operations Research, 44(4):543–569, August 1996.
  * [41] Wes McKinney.  Data Structures for Statistical Computing in Python.  In Stéfan van der Walt and Jarrod Millman, editors, Proceedings of the 9th Python in Science Conference, pages 56 – 61, 2010\.
  * [42] Eugene Litvak Xin Ming Tu and Marcello Pagano.  On the informativeness and accuracy of pooled testing in estimating prevalence of a rare disease: Application to hiv screening.  Biometrika, 82:287–297, 06 1995.
  * [43] N. Zaman and N. Pippenger.  Asymptotic analysis of optimal nested group-testing procedures.  Probability in the Engineering and Informational Sciences, 30:547–552, 2016.
  * [44] Ugnė Čižikovienė and Viktor Skorniakov.  On a couple of unresolved group testing conjectures.  Communications in Statistics - Theory and Methods, 2021.  (in press).

##  Appendix A Proofs

######  Remark A.1.

In the proof of Theorem 2.1, Step 1 is a repetition of Lemma 1 in [17]. We
have decided to rewrite it here because it is very short and, along the way,
some notions used in the sequel appear.

_Proof of Theorem 2.1._ _Step 1._ We first show that, for any fixed
$n\in(2,\infty)$, there exists a unique $q_{n}\in(0,1)$ such that

$g(q_{n},n)=0,\quad g(q,n)<0\ \forall q\in(q_{n},1)\quad\text{and}\quad
g(q,n)>0\ \forall q\in(0,q_{n}).$

To this end, note that

$\frac{\partial}{\partial q}g(q,n)=-2nq^{n-1}+(2n-1)q^{2n-2}=\\\
-2nq^{n-1}\left(1-\frac{2n-1}{2n}q^{n-1}\right)<0\ \forall n\in(2,\infty).$

Thus, given $n\in(2,\infty)$, $q\mapsto g(q,n)$ is decreasing on $(0,1)$.
Since $g(0+,n)=\frac{2}{n}>0$ and $g(1-,n)=\frac{2}{n}-1<0$, the claim holds
true.

_Step 2._ From _Step 1_ it follows that, for any fixed $q\in(0,1)$, we have a
well defined function $n\mapsto q_{n}$ which is given implicitly by equation
$g(q_{n},n)=0$. Since this function is continuous, its range $I\subset(0,1)$
is an interval. Further, for $\varepsilon\in(0,1)$ and $n=2+\varepsilon$,

$g(1-\varepsilon,n)=\frac{2}{n}-2(1-n\varepsilon+O(\varepsilon^{2}))+(1-(2n-1)\varepsilon+O(\varepsilon^{2}))=\\\
\frac{-\varepsilon}{2+\varepsilon}+\varepsilon+O(\varepsilon^{2})>0,$

provided $\varepsilon$ is small enough. Hence, analysis accomplished in _Step
1_ implies that $q_{2+\varepsilon}\in(1-\varepsilon,1)$. Therefore,
$I=(q_{*},1)$ for some $q_{*}\in(1/2,1)$. To justify the lower bound $1/2$,
note that

$g(1/2,n)>0\Longleftrightarrow 2^{n}>n\left(1-\frac{1}{2^{n}}\right).$

Since the right hand side holds true for all $n>2$, it follows that
$g(1/2,n)>0\ \forall n>2$. Also note that $g(q,n)>0\
\forall\,(q,n)\in(0,q_{*})\times(2,\infty)$ since an opposite contradicts the
definition of $q_{*}$.

_Step 3._ Fix $q\in(q_{*},1)$. By _Step 1–Step 2_ , $g(q,n)=0$ has at least
one solution $n=n(q)$ (suffices it to take $n$ such that $q_{n}=q$). To show
that there are two solutions, put $c=\frac{1}{2q}$, make a change of variable
$q^{n}=x$, and rewrite $g(q,n)=0$ in a form

$-\ln q=-\ln x(x-cx^{2}).$ (A.1)

Consider function $h(x)=-(x-cx^{2})\ln x$ for $x\in(0,1)$. Note that

$\frac{\mathrm{d}}{\mathrm{d}x}h(x)=-\left((1-cx)+(1-2cx)\ln
x\right)=0\Longleftrightarrow-\ln x=\frac{1-cx}{1-2cx}.$ (A.2)

Since $1-cx>1-c>0$ and $-\ln x>0$ for all $x\in(0,1)$, it follows that
$1-2cx>0$ as well, provided $x$ solves (A.2). Therefore, the range of possible
solutions of (A.2) shrinks to $(0,q)$. Moreover, relationships
$\lim\limits_{x\to 0+}\frac{\mathrm{d}}{\mathrm{d}x}h(x)=\infty$,
$\lim\limits_{x\to 1-}\frac{\mathrm{d}}{\mathrm{d}x}h(x)=c-1<0$ imply that
(A.2) has at least one solution. Since $\frac{1-cx}{1-2cx}+\ln
x=1+\frac{cx}{1-2cx}+\ln x$ increases on $(0,q)$, the solution is unique.
Denote it $x_{0}$. Based on the sign of the derivative, we have that
$h\uparrow$ on $(0,x_{0})$ and $h\downarrow$ on $(x_{0},1)$. Hence, at
$x_{0}$, $h$ attains its maximum and (A.1) admits exactly two solutions if
$h(x_{0})>-\ln q$, one solution if $h(x_{0})=-\ln q$, and has no solutions if
$h(x_{0})<-\ln q$. By the choice of $q$ (recall that $q>q_{*}$), the last case
can not hold. To exclude the second one, note that the function
$\left(\frac{1}{2},\frac{1}{2q_{*}}\right)\ni c\mapsto x_{0}(c)$ is well
defined and decreasing since

$-\ln x_{0}=\frac{1-cx_{0}}{1-2cx_{0}}\Rightarrow\\\
-\frac{\mathrm{d}}{\mathrm{d}c}\ln
x_{0}(c)=-\frac{\frac{\mathrm{d}}{\mathrm{d}c}x_{0}(c)}{x_{0}(c)}=\frac{\mathrm{d}}{\mathrm{d}c}\frac{1-cx_{0}}{1-2cx_{0}}=\frac{x_{0}+c\frac{\mathrm{d}}{\mathrm{d}c}x_{0}(c)}{(1-2cx_{0})^{2}}\Rightarrow\\\
\frac{\mathrm{d}}{\mathrm{d}c}x_{0}(c)=\frac{-x_{0}^{2}}{(1-2cx_{0})^{2}+cx_{0}}<0.$

Therefore, $(q_{*},1)\ni q\mapsto x_{0}(q)$ is increasing. Taking into account
that $q\mapsto-\ln q$ is decreasing, we finally deduce that $h(x_{0})>-\ln q$
for all $q\in(q_{*},1)$. The monotonicity of $q\mapsto x_{0}(q)$ and
$q\mapsto-\ln q$ also leads to conclusion that $q_{*}$ can be solved from
equation

$-\ln q_{*}=h(x_{0}(q_{*}))$

along with a unique $n_{*}\in(2,\infty)$. Hence (i).

_Step 4._ Assume the setting of _Step 3_. Let
$0<x_{L}=q^{n_{U}}<x_{U}=q^{n_{L}}<q$ denote two solutions of (A.1). By above,
$h(x)>-\ln q\Longleftrightarrow x\in(x_{L},x_{U})$. Reverting to
$(0,\infty)\ni n\mapsto g(q,n)$, this reads as $g(q,n)<0\Longleftrightarrow
n\in(n_{L},n_{U})$. Note that $n\mapsto g(q,n)>0$ in the neighborhood of
$\infty$. Also, from _Step 1_ – _Step 2_ , we have that $n\mapsto g(q,n)>0$ in
the right neighborhood of 2 and that $n_{U}\in(2,\infty)$. Therefore,
continuity of $n\mapsto g(q,n)>0$ yields that $n_{L}\in(2,\infty)$ as well.
Finally, it is clear that $n_{L}<n_{*}<n_{U}$ (see figure 5 for a graphical
illustration). Hence (ii).

_Step 5._ In this step, we identify the number and location of zeroes of the
derivative of $(2,\infty)\ni n\mapsto g(q,n)$ having fixed $q\in(q_{*},1)$.
From analysis given in _Step 4_ , it follows that $(2,\infty)\ni n\mapsto
g(q,n)$ has at least two extremes: there must be a minimum in $(n_{L},n_{U})$
(since function is negative here), and maximum in $(n_{U},\infty)$ (since the
function is positive here and $\lim_{n\to\infty}g(q,n)=0+$. To see that there
are no other extremes, consider equation

$\frac{\partial}{\partial n}g(q,n)=-\frac{2}{n^{2}}-2q^{n}\ln q+2q^{2n-1}\ln
q=0,$ (A.3)

and rewrite it by making use of notions introduced in _Step 3_ as follows:

$-x\ln x=\sqrt{-\ln
q}\sqrt{\frac{x}{1-2cx}},x\in\left(0,\frac{1}{2c}\right)=(0,q).$ (A.4)

Next, note that:

  * •

$h_{1}(x)=-x\ln x$ is strictly convex–up and positive on $(0,1)$ with
$\lim\limits_{x\to 0+}h_{1}(x)=\lim\limits_{x\to 1-}h_{1}(x)=0+$;

  * •

$h_{2}(x)=\sqrt{\frac{x}{1-2cx}}$ is strictly positive and increasing on
$\left(0,\frac{1}{2c}\right)=(0,q)$, it has one inflection point and
$\lim\limits_{x\to 0+}h_{2}(x)=0$, $\lim\limits_{x\to q-}h_{2}(x)=\infty$.

Taking this information into account, we conclude that (A.4) can have at most
two solutions and confirm thereby the assertion stated above.

_Step 6._ It remains to justify expression (2.4). Let

$n(q,t)=\frac{1}{p^{\frac{2}{3}}}+\frac{1}{2p^{\frac{1}{3}}}+0.2+3p^{2}+t,\quad
t\in[0,1].$ (A.5)

It suffices to prove that, for all $q\in[0.755,1)$, the following statements
hold true:

  * (a)

$\max(g(q,n(q,0)),g(q,n(q,1))<0$; and

  * (b)

$\exists t\in[0,1]:\frac{\partial}{\partial n}g(q,n)\Big{|}_{n=n(q,t)}=0$.

Analytical calculations behind (a) and (b) are standard yet very lengthy and
tedious. Therefore, we omit the details and end up with a graphical proof and
a sketch of the analytical one.

Figures 6–7 show graph of $q\mapsto\max(g(q,n(q,0)),g(q,n(q,1)))$ from which
it is evident that (a) holds. Analytical proof consists of the following
steps.

  * (s1)

Calculate $\frac{\partial}{\partial q}g(q,n(q,i)),i=0,1$.

  * (s2)

Check that $\frac{\partial}{\partial q}g(q,n(q,i))<0,i=0,1$ on $[0.755,1)$ and
deduce that $g(q,n(q,i))$ decrease on [0.755,1).

  * (s3)

Conclude that (a) indeed holds since

$g(0.755,n(0.755,0))\approx-0.002258\quad\text{and}\quad
g(0.755,n(0.755,1))\approx-0.013690.$

Turning to (b), first rewrite (A.3) as follows:

$-{n^{2}}q^{n}\ln q(1-q^{n-1})=1.$

Next, consider function $h(t,q)=-{n^{2}(q,t)}q^{n(q,t)}\ln
q(1-q^{n(q,t)-1})-1$ with $n(q,t)$ given by (A.5) and $q\in[0.755,1)$. Since
$t\mapsto h(t,q)$ is continuous, it suffices to show that, for any
$q\in[0.755,1)$, $h(0,q)<0$ and $h(1,q)>0$. Figure 8 shows graphs of $q\mapsto
h(0,q),q\mapsto h(1,q)$. These confirm (b). Considering analytical part, the
following is the suggested route.

  * (s1)

Calculate $\frac{\partial^{2}}{\partial q^{2}}h(i,q),i=0,1$.

  * (s2)

Check that $\frac{\partial^{2}}{\partial q^{2}}h(0,q)>0$ whereas
$\frac{\partial^{2}}{\partial q^{2}}h(1,q)<0$ on $[0.755,1)$ and deduce that
$h(0,q)$ is convex downwards whereas $h(1,q)$ is convex upwards on [0.755,1).

  * (s3)

By making use of Taylor’s expansion, check that $\lim\limits_{q\to
1-}h(0,q)=0-$ and $\lim\limits_{q\to 1-}h(1,q)=0+$.

  * (s4)

Conclude that (b) indeed holds since

$h(0,0.755)\approx-0.2645889\quad\text{and}\quad h(1,0.755)\approx 0.081749.$

∎

_Proof of Corollary 2.1._ Uniqueness of $q_{5}$ was established in Step 1 of
the proof of Theorem 2.1. Hudgens and Kim [17] (Lemmas 2, 7, and 14) have
demonstrated that $n_{opt}(q)\not\in\\{2,3,4\\}\ \forall q\in(0,1)$. From
their results we also have that $n_{opt}(q)=5\ \forall q\in(q_{*},0.755]$. It
is straightforward to verify that

$\forall\ q\in(q_{5},0.755]\
\left\lceil\frac{1}{p^{\frac{2}{3}}}+\frac{1}{2p^{\frac{1}{3}}}+3p^{2}+1.2\right\rceil=5.$

Hence the claim in the region $(q_{5},0.755]$. For $q\in[0.755,1)$, it follows
from Theorem 2.1 by noting that at least one of numbers in the set (2.5)
belongs to $\\{n\in(2,\infty):t(q,n)<1\\}$ because of (a) in Step 6 of proof
of Theorem 2.1. ∎

##  Appendix B Tables and Figures

The material contained in this appendix was produced by making use of an open
source computer algebra system SymPy [26] and the python packages constituting
the core kit for scientific numerical programming with python: NumPy [15],
SciPy [39], Matplotlib [18] and Pandas [41], [29]. The table accompanying
Subsection 2.2.1 is placed in our Open Science Framework repository (file
A2_table.csv). For the prevalence $p$ ranging over $(0,0.249790)$ at a step
equal to $10^{-6}$ and for each scheme considered there, the following
information is given: optimal (rounded) pool size $N_{opt}(q)$ and gain
$G(q)=1-t(q,N_{opt}(q))$. In case of A2, $n^{(A2)}_{opt}(q)$ and
$t^{(A2)}(q,n^{(A2)}_{opt}(q))$ are reported (see Subsection 2.2.1, Remark 2.2
for the details).

Figure 1: Graph showing behaviour of optimal pool sizes and gains on original
and log–log scales. Figure 2: Graph showing behaviour of maximal tested pool
size for different schemes. Figure 3: Graph of $\\{2,3,\ldots\\}\ni
n\mapsto\sup_{q}(t^{(A2)}(q,n)-t^{(A2)}(q,n_{opt}^{(A2)}))$ Figure 4: Graph of
$\\{2,3,\ldots\\}\ni
n\mapsto\operatorname{E}_{\pi}(t^{(A2)}(q,n)-t^{(A2)}(q,n_{opt}^{(A2)}))^{2}$
Figure 5: Graph illustrating relationships of $n\mapsto q_{n}$ and related
quantities. The lower curve corresponds to $q=0.86$. Figure 6: Graph of
$q\mapsto\max(g(q,n(q,0)),g(q,n(q,1)))$ for ${q\in[0.765,1)}$. For reference,
function identically equal to 0 is plotted in red. Figure 7: Graph of
$q\mapsto\max(g(q,n(q,0)),g(q,n(q,1))$ for ${q\in[0.755,0.765]}$. For
reference, function identically equal to 0 is plotted in red. Figure 8: Graphs
of $q\mapsto h(0,q),q\in[0.755,1)$ (plotted in blue) and $q\mapsto
h(1,q),q\in[0.755,1)$ (plotted in green).