# Eddy-Viscous Modeling and the Topology of Extreme
Circulation Events in Three-Dimensional Turbulence

G.B. Apolinário1, L. Moriconi2, R.M. Pereira3, and V.J.
Valadão2111Corresponding author<EMAIL_ADDRESS>1ENS de Lyon, CNRS,
Laboratoire de physique, F-69342 Lyon, France  2Instituto de Física,
Universidade Federal do Rio de Janeiro, C.P. 68528, CEP: 21941-972, Rio de
Janeiro, RJ, Brazil,  3Instituto de Física, Universidade Federal Fluminense,
24210-346 Niterói, RJ, Brazil.

###### Abstract

We discuss the role of particular velocity field configurations – instantons,
for short – which are supposed to dominate the flow during the occurrence of
extreme turbulent circulation events. Instanton equations, devised for the
stochastic hydrodynamic setup of homogeneous and isotropic turbulence, are
applied to the interpretation of direct numerical simulation results. We are
able in this way to model the time evolution of extreme circulation events for
a broad range of scales, through the combined use of eddy viscosity
phenomenology and exact creeping instantons. While this approach works well
for the core of circulation instantons, it fails to describe their tails. In
order to overcome this difficulty, we put forward a numerical treatment of the
axisymmetric instanton equations. Circulation instantons are then found to
have a surprising topological structure, which consists of a system of paired
counter-rotating vortex rings centered around the symmetry axis of a
background axisymmetric vortical flow.

##  I Introduction

Predicting extreme events is a fundamental problem in completely distinct
domains of science and technology, as in astrophysics aschwanden2019self ;
larrodera2021estimation , climatology bouchet2014stochastic ;
laurie2015computation , wind power generation couto2021identification , and a
further broad list of other important contemporary fields such as finance
rocco2014extreme , meteorology yiou2004extreme ; grotjahn2016north  and
seismology soloviev2008transformation . Extreme events often take place in
relatively very short time scales, but their consequences can be perennially
disastrous as in the case of earthquakes, hurricanes and financial market
crashes.

Navier-Stokes dynamics, our specific focus in this work, is a particularly
interesting instance for the study of extreme events. The formation of intense
energy-carrying localized structures in turbulent flows is a fingerprint of
the non-local and nonlinear aspects of the underlying equations of motion.
Fast spatial variations of the dissipation field or of the entangled system of
vortex tubes which confine strong vorticity yoshida2000statistical ;
kaneda2003energy  are, however, difficult to observe, and decades went by
until it was realized that the self-similar cascade picture of turbulence,
originally proposed by Kolmogorov (K41) frisch , should be revised to include
the existence of such extreme events. As it is well known, the predicted K41
scaling exponents for the statistical moments of velocity differences fail to
reproduce the respective numerically and experimentally measured values
batchelor1949nature ; anselmet1984high ; vincent1991spatial  at high moment
orders. Commonly understood as an important feature of the intermittency
phenomenon frisch ; kaneda2003energy , these deviations of the K41 scaling
exponents have been accurately obtained through multifractal formulations
kolmogorov1962refinement ; obukhov1962 ; novikov1994intermittency ;
she1994universal ; apolinario2020vortex ; pereira2022hard , but more
fundamental derivations are still in order.

A systematic analytical strategy to address the problem of intermittent
velocity fluctuations was initially introduced in the framework of Burgers
turbulence burgers1948mathematical . Asymptotic results were found for the
right tail of the velocity gradient probability distribution functions
gurarie1996instantons , based on the steepest-descent approximation as
implemented in the Martin-Siggia-Rose-Janssen-de-Dominicis
martin1973statistical ; janssen1976lagrangean ; dominicis1976technics  (MSRJD)
functional formalism. Subsequent theoretical and numerical efforts have been
made to describe the left tail of Burgers velocity gradient and velocity
difference PDFs balkovsky1997intermittency ; chernykh2001large ;
moriconi2009instanton ; grafke2015relevance ; apolinario2019onset . Despite
the relevant progress achieved for the Burgers and effective Lagragian
turbulent models johnson2016closure ; grigorio2017instantons ;
apolinario2019instantons ; alqahtani2021extreme , the description of
intermittency in three-dimensional turbulence is still very open to
investigation along the lines of the MSRJD approach.

In this work, we use the MSRJD functional formalism to study the scaling
properties of velocity circulation in three-dimensional homogeneous and
isotropic turbulence, motivated by the fact that velocity circulation
fluctuations are known to provide clear signatures of turbulence intermittency
umeki1993probability ; migdal1994loop ; cao1996properties ; benzi1997self .
Velocity circulation is a natural tool to probe localized vortex structures,
which contain most of the turbulent kinetic energy in high Reynolds number
flows kaneda2003energy ; frisch . The first formally structured attempt to
explore velocity circulation as the main actor of turbulent fluctuations dates
back to 1994 migdal1994loop . Promising results about the circulation
probability distribution functions (cPDFs) were then advanced, but further
developments and validation studies faced serious obstacles in the following
years, due to computational limitations. More recently, once hardware
shortcomings were overcome, the subject of circulation statistics was
vigorously revisited under the light of much improved high-resolution
simulations, not only in classical iyer2019circulation ; iyer2021area  but
also in quantum turbulence muller2021intermittency ; polanco2021vortex .
Interesting perspectives on the multifractal and structural views of classical
turbulence cascades have been put forward since then apolinario2020vortex ;
pereira2022hard ; moriconi2021multifractality ; moriconi2022statistics .

Throughout our considerations, the circulation variable is defined as

$\Gamma_{R}=\oint_{\mathcal{C}_{R}}v_{i}dx_{i}=\iint_{\mathcal{D_{R}}}\omega_{i}dA_{i}\
,\ $ (1)

where $\mathcal{D_{R}}$ is a circular domain of radius $R$ and
${\mathcal{C}_{R}}$ is the circular contour that encloses $\mathcal{D_{R}}$.
The fields $v_{i}$ and $\omega_{i}$ are, respectively, the velocity and
vorticity fields related through $\omega_{i}=\epsilon_{ijk}\partial_{j}v_{k}$
(we use Einstein notation for the summation over repeated indices, unless
explicitly stated). We are interested in discussing velocity field
configurations – instantons, for short – which are associated to extreme
circulation events, that is, large deviations of $\Gamma_{R}$. On one front of
analysis, we model them as exact creeping flow structures creeping  related to
turbulent eddy-viscosity boussinesq1877essai  time scales. The approach proves
to be phenomenologically meaningful, once we validate it by means of extensive
turbulence databases. On another front, we numerically solve the axisymmetric
nonlinear instanton equations, to find that the circulation instantons have a
remarkable topological structure described by the superposition of two
counter-rotating vortex rings that share the same symmetry axis and carry
opposite helicity in the presence of a background axial vorticity field.

This paper is organized as follows. Sec. II briefly reviews the main technical
aspects of the MSRJD functional formalism. We, then, obtain closed analytical
solutions for the circulation instantons in the viscous limit of the Navier-
Stokes equations and show their somewhat unexpected usefulness to model the
time evolution of extreme circulation events, from the analysis of direct
numerical simulation (DNS) data. The nonlinear instanton equations are
introduced in Sec. III, and their numerical solutions are carefully discussed,
with emphasis placed on the qualitative differences between them and the
creeping solutions derived in Sec. II. Finally, in Sec. IV, we summarize our
findings and indicate directions of further research.

##  II The functional approach and creeping instantons

###  II.1 MSRJD Formalism

Many physical systems which exhibit extreme events among their dynamical
states are suitable for the application of large deviation techniques
weiss1995introduction ; evertsz1995large ; ragone2018computation . It is often
possible to address in these cases the MSRJD methodology martin1973statistical
; janssen1976lagrangean ; dominicis1976technics , where one typically deals
with differential equations that lead to solutions – the so-called instantons
– which dominate the probability distribution tails that describe the
occurrence of extreme events.

We are interested in applying these ideas to three-dimensional turbulence in
its stochastic hydrodynamics formulation brenig1976stochastic . The stochastic
incompressible Navier-Stokes equations are written as

$\partial_{t}v_{i}+v_{j}\partial_{j}v_{i}-\nu\partial^{2}v_{i}=-\partial_{i}P+f_{i}\
,\ $ (2) $\partial_{i}v_{i}=0\ ,\ $ (3)

where $v_{i}=v_{i}(\vec{x},t)$ is the velocity field, $P=P(\vec{x},t)$ is the
pressure field ensuring the incompressibility constraint (Eq. 3), $\nu$ is the
kinematic viscosity, and $f_{i}=f_{i}(\vec{x},t)$ is a zero-average Gaussian
random field used to model large-scale forcing, correlated as

$\langle
f_{i}(\vec{x},t)f_{j}(\vec{x}^{\prime},t^{\prime})\rangle=\delta_{ij}\chi(\vec{x}-\vec{x}^{\prime})\delta(t-t^{\prime})\
.\ $ (4)

The exact form of $\chi(\vec{x}-\vec{x}^{\prime})$ is not phenomenologically
relevant insofar as its correlation length $L$ is much larger than the
investigated length scales. Assuming homogeneity and isotropy, the correlation
kernel $\chi(\vec{x}-\vec{x}^{\prime})$ depends only on
$|\vec{x}-\vec{x}^{\prime}|$ and is parametrized by its amplitude
$D_{0}=\chi(\vec{0})$ and its correlation length $L$, which can be
conveniently defined as $L^{2}=\chi(0)/\chi^{\prime\prime}(0)$.

In our particular application of the MSRJD formalism, we deal with a
functional measure which is path-integrated over $v_{i}(\vec{x},t)$ and an
auxiliary field $p_{i}(\vec{x},t)$ ivashkevich1997symmetries ;
grafke2015instanton . Lagrange multipliers $Q(\vec{x},t)$ and $\lambda$ are
also introduced, to ensure incompressibility and to fix the circulation
$\Gamma_{R}$ at a given time instant, respectively. More concretely, the cPDF
is evaluated in the MSRJD formalism as

$\displaystyle\rho(\Gamma_{R})=\left\langle\delta\left(\Gamma_{R}-\oint_{\mathcal{C}}v_{i}(\vec{x},0)dx_{i}\right)\right\rangle=$
$\displaystyle=\mathcal{N}^{-1}\int
D[\vec{v}]D[P]D[\vec{p}]D[Q]\int_{-\infty}^{\infty}d\lambda\
\exp\left\\{-\frac{S[\vec{v},P,\vec{p},Q,\lambda]}{g^{2}}\right\\}\ ,\ $ (5)

where the pseudo-statistical weight (it is a complex number) contributions for
$\rho(\Gamma_{R})$ are given as the exponential of the MSRJD action,
$S[\vec{v},P,\vec{p},Q,\lambda]$, specified below, and $\mathcal{N}$ is an
unimportant normalization constant (to be suppressed hereafter). A single
dimensionless parameter $g^{2}=\chi(\vec{0})L^{4}\nu^{-3}$ can be factored out
from the action and its role is to control the noise strength and to set what
should be considered an extreme event: these are the ones associated to MSRJD
actions which are much larger than $g^{2}$.

The above functional integration represents the probability density functional
to find velocity circulation $\Gamma_{R}$ as defined in Eq. (1), at a
particular time instant, namely $t=0$, due to the cumulative effect of forcing
since the remote past, say $t\rightarrow-\infty$. All fields are assumed to
vanish in the remote past. Moreover, instantons are restricted to $t<0$, such
that auxiliary fields are imposed to vanish at $t\rightarrow 0^{+}$ as a
boundary condition chernykh2001large . Measuring length in units of $L$ and
time in units of the viscous time scale

$\tau_{\nu}=\frac{L^{2}}{\nu}\ ,\ $ (6)

a dimensionless form of the MSRJD action can be written as

$\displaystyle\tilde{S}[\vec{v},P,\vec{p},Q,\lambda]=-\frac{1}{2}\int
d^{3}xd^{3}x^{\prime}dt\
p_{i}(\vec{x},t)\tilde{\chi}(\vec{x}-\vec{x}^{\prime})p_{i}(\vec{x}^{\prime},t)+\int
d^{3}xdt\ Q\partial_{i}v_{i}+$ $\displaystyle+\int d^{3}xdt\
p_{i}(\partial_{t}v_{i}+v_{j}\partial_{j}v_{i}+\partial_{i}P-\partial^{2}v_{i})-\lambda\left(\Gamma_{R/L}-\oint_{\mathcal{C}_{R/L}}\\!\\!dx_{i}v_{i}(\vec{x},t=0)\right)\
,\ $ (7)

where, now, $\tilde{\chi}(\vec{x}-\vec{x}^{\prime})$ is the force-force
correlator with unit correlation length. The connection between $g$ and the
Reynolds number Re is made by the assumption that a well defined inertial
range develops when the dissipation length
$\eta_{K}=(\nu^{3}/\chi(\vec{0}))^{1/4}$ is much smaller than the correlation
length $L$. Since the characteristic velocity scale is
$U=(\chi(\vec{0})L)^{1/3}$, we obtain, as a consequence, that Re
$=UL/\nu=g^{2/3}$.

Instantons correspond to the extrema of the MSRJD action and are assumed to
provide, from Eq. (II.1), the dominant saddle-point contributions in the path-
integral evaluation of extreme event probabilities. They are found as
solutions of the variational principle $\delta\tilde{S}=0$, which amounts here
to the instanton equations,

$(\partial_{t}-\partial^{2})v_{i}+v_{j}\partial_{j}v_{i}+\partial_{i}P=\int
d^{3}x^{\prime}\
p_{i}(\vec{x}^{\prime},t)\tilde{\chi}(\vec{x}-\vec{x}^{\prime})\ ,\ $ (8)
$(\partial_{t}+\partial^{2})p_{i}+v_{j}\partial_{j}p_{i}+v_{j}\partial_{i}p_{j}+\partial_{i}Q=\lambda\delta(t)\oint_{\mathcal{C}_{R/L}}\\!\\!dx^{\prime}_{i}\delta^{3}(\vec{x}-\vec{x}^{\prime})\
,\ $ (9)

subject to the constraints,

$\partial_{i}v_{i}=0\ ,\ \partial_{i}p_{i}=0\ ,\ $ (10)
$\Gamma_{R/L}=\oint_{\mathcal{C}_{R/L}}\\!\\!dx^{\prime}_{i}v_{i}(x^{\prime},t=0)\
.\ $ (11)

By solving the above equations, one gets asymptotic expressions for the cPDF
tails, from the mapping of the values of the Lagrange multiplier $\lambda$ to
the final velocity circulation $\Gamma_{R/L}[\vec{v}(\lambda)]$. The instanton
calculus yields, in this way,

$\rho(\Gamma_{R/L})\simeq
C\exp\left[-\frac{1}{g^{2}}\tilde{S}_{c}(\lambda)\right]$ (12)

as an asymptotic approximation that holds for $g\ll\tilde{S}_{c}(\lambda)$,
where $\lambda=\lambda(\Gamma_{R/L})$ and $C$ is a normalization constant. A
more mathematically rigorous discussion of the arguments leading to results
which are analogous to (12) is the essential subject of large deviation theory
dembo1998large . As a direct consequence of the Gärtner-Ellis theorem, a
probability density function necessarily follows the Large Deviation Principe
(Eq. (12), in our case) if the associated action is strictly convex and the
map between $\lambda$ and the observable of interest is differentiable for an
arbitrary $\lambda$. Non-convex actions or ill-defined $\lambda$-maps can be
nevertheless adapted to large deviation evaluations through the use of
alternative integration measures alqahtani2021instantons .

###  II.2 Creeping Instantons and Eddy Viscous Modeling

Instanton equations are rarely amenable to analytical treatment without the
help of phenomenological or empirical inputs about the underlying dynamics.
Approximate solutions of (8-11) and of analogous magnetohydrodynamic equations
were formerly carried out under restrictive assumptions about the roles of
symmetries and the strain rate field moriconi1998circulation ;
moriconi2002circulation . Somewhat surprisingly, as we will discuss in the
following, interesting information on extreme events can be obtained from the
exact form of viscous (creeping) instantons.

To work out the creeping instantons, we introduce the diffusion Green’s
function with vanishing boundary condition in $\mathbb{R}^{3}$,

$G(\vec{x},t)=\frac{1}{(4\pi|t|)^{-\frac{3}{2}}}\exp\left(-\frac{|\vec{x}|^{2}}{4|t|}\right)\
.\ $ (13)

Neglecting the nonlinear terms in Eqs. (8-11), the velocity and the conjugate
solutions read

$v^{(1)}_{i}(\vec{x},t)=\int_{-\infty}^{t}dt^{\prime}\int
d^{3}x^{\prime}d^{3}x^{\prime\prime}G(\vec{x}-\vec{x}^{\prime},t-t^{\prime})\tilde{\chi}(\vec{x}^{\prime}-\vec{x}^{\prime\prime})p_{i}^{(1)}(\vec{x}^{\prime\prime},t^{\prime})\
,\ $ (14)
$p^{(1)}_{i}(\vec{x},t)=\lambda\oint_{\mathcal{C}_{R/L}}dx_{i}^{\prime\prime}\int_{t}^{0}dt^{\prime}\int
d^{3}x^{\prime}G(\vec{x}-\vec{x}^{\prime},t-t^{\prime})\delta(t^{\prime})\delta^{3}(\vec{x}^{\prime}-\vec{x}^{\prime\prime})\
.\ $ (15)

The MSRJD action, Eq. (II.1), additionally reduces to

$S^{(1)}_{0}[v_{i},p_{i}]=\frac{\lambda}{2g^{2}}\int_{\mathcal{C}_{R/L}}dx_{i}v_{i}^{(1)}(\vec{x},0)=\frac{\lambda\Gamma^{(1)}_{R/L}(\lambda)}{2g^{2}}\
.\ $ (16)

Defining, at this point, the explicit form of the force-force correlator as

$\tilde{\chi}(\vec{x}-\vec{x}^{\prime})=\frac{1}{(2\pi)^{3/2}}\exp\left(-\frac{|\vec{x}-\vec{x}^{\prime}|^{2}}{2}\right)\
,\ $ (17)

we obtain from (14) and (15) the creeping instanton solution by setting
$x_{3}$ to be perpendicular to the contour $\mathcal{C}_{R/L}$

$v^{(1)}_{i}(\vec{x},t)=\lambda\pi\sqrt{\frac{R/L}{x_{\perp}}}\epsilon_{3ji}\frac{x_{j}}{x_{\perp}}\int_{0}^{\sqrt{\frac{x_{\perp}R/L}{1-2t}}}du\
e^{-\left(\frac{(R/L)^{2}+x_{\perp}^{2}+x_{3}^{2}}{2x_{\perp}(R/L)}\right)u^{2}}I_{1}(u^{2})\
,\ $ (18)

where $I_{1}(\cdot)$ is the modified Bessel function of the first kind and
$x_{\perp}^{2}=x_{1}^{2}+x_{2}^{2}$. We remark that part of the spatio-
temporal dependence of the velocity field is encoded in the upper limit of the
above integration. It is interesting to note that this integral can be exactly
computed only at points that belong to the contour ${\mathcal{C}_{R/L}}$,
defined by $x_{3}=0$ and $|x_{\perp}|=R/L$. We find, in cylindrical
coordinates,

$\vec{v}^{(1)}(\mathcal{C}_{R/L},t)=\lambda\frac{\pi}{6}\frac{R^{3}/L^{3}}{(1-2t)^{3/2}}\
{}_{2}F_{2}\left(\frac{3}{2},\frac{3}{2},\frac{5}{2},3;-2\frac{R^{2}/L^{2}}{1-2t}\right)\hat{\theta}\
,\ $ (19)

where ${}_{2}F_{2}$ is the hypergeometric function of $2+2$ entries. The time-
dependent instanton circulation is readily obtained as

$\Gamma^{visc}_{R}(t)\equiv\oint_{\mathcal{C}_{R/L}}\vec{v}^{(1)}(\vec{x},t)\cdot
d\vec{x}=\lambda\frac{\pi^{2}}{3}\frac{R^{4}/L^{4}}{(1-2t)^{3/2}}\
{}_{2}F_{2}\left(\frac{3}{2},\frac{3}{2},\frac{5}{2},3;-2\frac{R^{2}/L^{2}}{1-2t}\right)\
.\ $ (20)

Since the circulation is a linear function of $\lambda$, the resulting cPDF is
Gaussian. Recovering the original circulation units, its variance
$\sigma_{R}^{2}$ is found to satisfy

$\frac{\sigma_{R}^{2}}{\nu^{2}{\hbox{Re}^{3}}}=\frac{\pi^{2}}{3}\frac{R^{4}}{L^{4}}\
{}_{2}F_{2}\left(\frac{3}{2},\frac{3}{2},\frac{5}{2},3;-2\frac{R^{2}}{L^{2}}\right)\
.\ $ (21)

We have, for $R\ll L$,

${}_{2}F_{2}\left(\frac{3}{2},\frac{3}{2},\frac{5}{2},3;-2\frac{R^{2}}{L^{2}}\right)=1-\frac{3}{5}\left(\frac{R}{L}\right)^{2}+{\cal{O}}((R/L)^{4})\
,\ $ (22)

such that, in this asymptotic limit,

$\tilde{\sigma}_{R}^{2}\equiv\frac{3\sigma_{R}^{2}L^{4}}{(\pi
R^{2})^{2}\nu^{2}{\hbox{Re}}^{3}}\simeq 1\ .\ $ (23)

Although (23) is expected to hold only for creeping flows, it may be possibily
relevant even for higher Reynolds number solutions, once scales are probed and
the MSRJD path-integration (II.1) is dominated by smooth velocity field
configurations. The educated guess that generalizes (23) to turbulent flows
sustained by alternative forcing mechanisms, where $L$ and Re are not
necessarily defined in terms of Gaussian force-force correlators, reads

$\frac{3\sigma_{R}^{2}L^{4}}{(\pi
R^{2})^{2}\nu^{2}{\hbox{Re}}^{3}}\simeq{\cal{O}}(1)\ .\ $ (24)

To investigate the correctness of (24), we have worked with four DNS datasets
of homogeneous and isotropic turbulent flows, publicly available from the
Johns Hopkins Turbulence Database (JHTDB) platform222For more information
about the public datasets see http://turbulence.pha.jhu.edu/. perlman2007data
; JHTD2 ; JHTD3 ; JHTD4 . The simulations were developed in periodic cubic
lattices with dimensions $1024^{3}$ (dataset I), $4096^{3}$ (datasets II), and
$8192^{3}$ (datasets III and IV), corresponding to Taylor-based Reynolds
numbers 418, 610, 613 and 1280 for datasets from I to IV. We have produced
from these datasets ensembles which contain respectively $1\times 10^{7}$,
$8\times 10^{6}$, $9\times 10^{6}$, and $4.5\times 10^{7}$ circulation
samples, computed by employing the first equality of Eq. (1). Denoting by
$\eta_{K}$ the Kolmogorov dissipation length frisch , the results reported in
Fig. 1 fully corroborate our expectations: relation (24) is in fact verified
for $R<\eta_{K}$ (it should be noted that small scale evaluations of
$\sigma_{R}$ are subject to relatively larger error bars, due to the
inaccuracy associated to the representation of a circular contour in a square
lattice, as closed polygonal line).

Fig. 1 also shows that the LHS of (24) leads, incidentally, to an excellent
collapse of data across the inertial range scales $30<R/\eta_{K}<300$.
Furthermore, the circulation variance has, in this range, a scaling dependence
with $R$ which is very well approximated by the K41 scaling exponent $8/3$
iyer2019circulation , as indicated by a dashed line Note that at fixed
integral scale and fixed energy injection rate per unit mass,
$\epsilon=U^{3}/L$, relation (24) implies that
$\sigma^{2}_{R}\propto\nu^{-1}R^{4}$. Insisting that an analogous result
should hold in the inertial range, we replace $\nu$ in this last relation by
the scale-dependent eddy viscosity smith1998renormalization ,

$\nu_{R}=\nu_{0}\left(\frac{R}{R_{0}}\right)^{\frac{4}{3}}\ ,\ $ (25)

where $\nu_{0}$ and $R_{0}$ are some reference viscosity and length scale
parameters, to obtain the observed scaling
$\sigma^{2}_{R}\propto\nu_{R}^{-1}R^{4}\propto R^{8/3}$.

As a brief historical digression, we point out that the concept of eddy
viscosity, long ago conceived by Boussinesq boussinesq1877essai  and much
later revived through the distant interaction approximation (DIA) formalism
kraichnan1964decay  and dynamical renormalization group techniques
forster1977large ; yakhot1986renormalization , is essentially an effective
transport parameter that allows for a coarse-grained description of the
turbulent cascade. This is the main physical motivation for a number of
numerically tractable models of turbulence, as the celebrated Smagorinsky sub-
grid formulation of the Navier-Stokes equations smagorinsky1963general ;
scotti1993generalized , which underlies the whole field of large eddy
simulations germano1991dynamic ; mason1994large ; piomelli1999large ;
piomelli2014large .

Figure 1: The dimensionless variance $\tilde{\sigma}_{R}^{2}$, defined in
(23), as a function of the contour radius $R$ for various Reynolds numbers
(symbols). Dashed lines represent $\tilde{\sigma}_{R}^{2}=1$ and
$\tilde{\sigma}_{R}^{2}\propto R^{-4/3}$.

The replacement of the molecular viscosity $\nu$ by the eddy viscosity
$\nu_{R}$ in (24) sounds at this point like a purely rhetoric remark. It is
possible, however, that creeping instantons, when combined with eddy viscosity
ideas, can in fact be used to model the time evolution of high Reynolds number
circulation instantons. The argument goes as follows. Consider the set of all
flow realizations which have evolved from the remote past and ended with
extreme circulation

$\Gamma_{R}\in[\bar{\Gamma}_{R}(0)-\delta\Gamma,\bar{\Gamma}_{R}(0)+\delta\Gamma]\,\
$ (26)

at time $t=0$, where $\bar{\Gamma}_{R}(0)\gg\delta\Gamma>0$ ($\delta\Gamma$ is
just a measurement bin). Let now $\bar{\Gamma}_{R}(t)$ represent the time
averaged circulation taken over the ensemble of all of these flow
realizations, which we refer to as the filtered instanton circulation.
Resorting to dynamic similarity and to the fact that $\eta_{K}/L\ll 1$ (i.e.,
Reynolds number is high), we write down the dimensionless circulation ratio
(time units are restored to the original ones, from now on),

$\hat{\Gamma}_{R}(t)\equiv\frac{\bar{\Gamma}_{R}(t)}{\bar{\Gamma}_{R}(0)}=f(t/\tau_{\nu},R/L,\eta_{K}/L)\simeq
f(t/\tau_{\nu},R/L,0)\equiv\tilde{f}(t/\tau_{\nu},R/L)\ .\ $ (27)

It follows from Eq. (27) that

$\hat{\Gamma}_{R}(\nu_{R}t/\nu)\simeq\tilde{f}(t/\tau_{\nu_{R}},R/L)\ .\ $
(28)

In consonance with eddy viscosity phenomenology, we assume that Eq. (28) is
scale invariant, viz.,

$\frac{d}{dR}\tilde{f}(t/\tau_{\nu_{R}},R/L)=0\ .\ $ (29)

We note that Eq. (29) is just a renormalization group equation
smith1998renormalization , which has, as general solution, the functional
relationship

$\tilde{f}(t/\tau_{\nu_{R}},R/L)\equiv h(t/\tau_{\nu})\ .\ $ (30)

This leads us, from (27), to

$\hat{\Gamma}_{R}(t)\simeq\tilde{f}(t/\tau_{\nu},R/L)=h(t^{\prime}/\tau_{\nu})\
,\ $ (31)

with

$t^{\prime}=t\frac{\nu}{\nu_{R}}=t\frac{\nu}{\nu_{0}}R^{-\frac{4}{3}}\ ,\ $
(32)

where we have used (25) to make the $R$-dependence explicit in the definition
of $t^{\prime}$.

In an analogous fashion, one has for the exact viscous solution (20)

$\hat{\Gamma}^{visc}_{R}(t)\equiv\frac{\Gamma^{visc}_{R}(t)}{\Gamma^{visc}_{R}(0)}=g(t/\tau_{\nu},R/L)\equiv\frac{1}{(1-2t/\tau_{\nu})^{3/2}}\frac{{}_{2}F_{2}\left(\frac{3}{2},\frac{3}{2},\frac{5}{2},3;-2\frac{R^{2}/L^{2}}{1-2t/\tau_{\nu}}\right)}{{}_{2}F_{2}\left(\frac{3}{2},\frac{3}{2},\frac{5}{2},3;-2\frac{R^{2}}{L^{2}}\right)}\
.\ $ (33)

It is natural to conjecture, furthermore, that in the limit
$R\rightarrow\eta_{K}$, the filtered and the viscous instanton circulations
match, that is,

$\displaystyle\hat{\Gamma}_{\eta_{K}}(t)=\tilde{f}(t/\tau_{\nu},\eta_{k}/L)=\hat{\Gamma}^{visc}_{\eta_{K}}(t)=g(t/\tau_{\nu},\eta_{K}/L)$
$\displaystyle\simeq
g(t/\tau_{\nu},0)\equiv\tilde{g}(t/\tau_{\nu})=\frac{1}{(1-2t/\tau_{\nu})^{3/2}}\frac{{}_{2}F_{2}\left(\frac{3}{2},\frac{3}{2},\frac{5}{2},3;0\right)}{{}_{2}F_{2}\left(\frac{3}{2},\frac{3}{2},\frac{5}{2},3;0\right)}\
.\ $ (34)

Pushing the validity of Eq. (30) down to dissipative scales, we get

$\tilde{f}(t/\tau_{\nu_{R}},R/L)=\tilde{f}(t/\tau_{\nu},\eta_{k}/L)\ ,\ $ (35)

which, in view of (34), leads to

$\tilde{f}(t/\tau_{\nu_{R}},R/L)\simeq\tilde{g}(t/\tau_{\nu})\ .\ $ (36)

Comparing (36) with (30), it follows that

$h(t/\tau_{\nu})\simeq\tilde{g}(t/\tau_{\nu})\ ,\ $ (37)

so that $\hat{\Gamma}_{R}(t)$ can be modeled, from (31), as the viscous
instanton solution $\tilde{g}(t^{\prime}/\tau_{\nu})$, that is,

$\hat{\Gamma}_{R}(t)=\hat{\Gamma}^{visc}_{0}(t^{\prime})=\frac{1}{(1-2t^{\prime}/\tau_{\nu})^{3/2}}\frac{{}_{2}F_{2}\left(\frac{3}{2},\frac{3}{2},\frac{5}{2},3;0\right)}{{}_{2}F_{2}\left(\frac{3}{2},\frac{3}{2},\frac{5}{2},3;0\right)}\
.\ $ (38)

In order to perform validation tests of Eqs. (31) and (38), we have
implemented a numerical filtering procedure to extract, from the raw data of
the JHTDB, the past time evolutions of extreme circulation events, which we
detail below.

###  II.3 Numerically Filtered Circulation Instantons

We have taken Dataset I for statistical analyses, since this is the dataset
which allows us to work with the largest ensemble of circulation time series,
defined for circular contours centered at equally spaced grid points. Letting
$\ell$ be the lattice parameter ($\ell\approx 2.1\eta_{K}$), we have studied
circulation fluctuations for four different radii, namely, $R=8\ell$,
$16\ell$, $32\ell$, and $64\ell$. For each given radius, we considered
$32^{3}$ circular contours oriented normally to each Cartesian direction, for
a total of $3\times 32^{3}$ contours per radius.

We conventionally define the set $\Lambda_{n}$ of extreme events as the
ensemble of circulations events whose absolute values $|\Gamma_{R}|$ reach
some multiple $n$ of the circulation standard deviation,
$\sigma_{R}=\sqrt{\langle\Gamma_{R}^{2}\rangle}$, within a small tolerance
window. In the language of Eq. (26), we consider the interval defined by
$\bar{\Gamma}_{R}(0)=n\sigma_{R}$ with a 0,5% tolerance, that is,
$\delta\Gamma=5\times 10^{-3}\bar{\Gamma}_{R}(0)$.

Figure 2: Filtered circulation instanton
$\hat{\Gamma}_{R}(t)\equiv\bar{\Gamma}_{R}(t)/\bar{\Gamma}_{R}(0)$ as a
function of the rescaled time $t^{\prime}=t(R/16\ell)^{-4/3}$ for a number of
radii and circulation sets $\Lambda_{n}$ at $R_{\lambda}=433$ (solid lines).
The eddy viscous modeling of the collapsed data is attained with the help of
Eq. (38) (dash-dotted line). Inset: the corresponding non-collapsed plots, as
functions of the original time variable $t$.

Once an extreme circulation event belonging to a given set $\Lambda_{n}$ is
identified, we assign it the observation time instant $t=0$ and save its
earlier time evolution. Then, a time dependent average $\bar{\Gamma}_{R}(t)$
over all saved series for each set $\Lambda_{n}$ is computed. Similar
filtering procedures have been applied in instanton studies of Burgers
turbulence grafke2013instanton , Lagrangian turbulence models
grigorio2017instantons ; apolinario2019instantons  and rogue wave formation
dematteis2018rogue ; dematteis2019experimental .

Our results for $\hat{\Gamma}_{R}(t)$ are shown in Fig. 2. Taking
$R_{0}=16\ell$ as an arbitrary reference length scale, we find that plots of
$\hat{\Gamma}_{R}(t)$ as a function $t(R/16\ell)^{-4/3}$ collapse reasonably
well for all the investigated radii and sets $\Lambda_{n}$, as predicted by
(31). Following now Eq. (38), we carry out an $L^{2}$-norm minimization of

$\sum_{R}||\hat{\Gamma}_{R}(\nu_{R}t/\nu)-\hat{\Gamma}^{visc}_{0}(t)||\ ,\ $
(39)

to adjust $\nu_{0}=4.5\times 10^{-5}\nu$ as the reference viscosity in (25).
Meaningful comparisons are then also reported in the same Fig. 2 for the eddy
viscous modeling of $\hat{\Gamma}_{R}(t)$, as suggested by (32) and (38).

Figure 3: Normalized filtered azimuthal velocity $\bar{v}_{\theta}(r,0,t)$
obtained from DNS data for different radii $R$ (vertically shifted for
clarity). Time evolution is color-coded (in units of the simulation time step)
and vertical lines indicate $r=R$ in each simulation. For comparison purposes,
part of a viscous instanton solution (dashed curve) is displayed for
$R=64\ell$.

In addition to the eddy viscous modeling of the filtered circulations
$\hat{\Gamma}_{R}(t)$ discussed above, one may wonder whether the functional
forms of filtered instanton velocity fields could be effectively modeled by
viscous instantons as well. With this aim, we take the symmetry axis of each
prescribed circular contour as the $z$-axis of a local Cartesian reference
system, and orient it in such a way that velocity circulations are rendered
positive at the event occurrence time. Then, we perform a similar time-
dependent averaging procedure over all events in $\Lambda_{n}$, as we did for
$\hat{\Gamma}_{R}(t)$, but now for velocity field configurations in
cylindrical coordinates,
$\vec{v}=v_{r}\hat{r}+v_{\theta}\hat{\theta}+v_{z}\hat{z}$. Due to isotropy,
the mean velocity components are functions of $r$, $z$, and $t$ only, i.e.,
$\langle v_{r}\rangle\equiv\bar{v}_{r}(r,z,t)$, $\langle
v_{\theta}\rangle\equiv\bar{v}_{\theta}(r,z,t)$, and $\langle
v_{z}\rangle\equiv\bar{v}_{z}(r,z,t)$.

Unfortunately, the radial and axial components $\bar{v}_{r}(r,z,t)$ and
$\bar{v}_{z}(r,z,t)$ turn out to develop relatively small intensities, which
prevents us from extracting a clear behavior out of the noise. The mean
azimuthal component $\bar{v}_{\theta}(r,0,t)$, on the other hand, was found to
be well-resolved in all studied cases. Fig. 3 shows comparisons between the
numerically filtered $\bar{v}_{\theta}(r,t)$ (normalized by
$\bar{v}_{\theta}(R,0)$) and the analogous normalized azimuthal velocity
obtained from the viscous solution (18), where the dimensionless time
$t/\tau_{\nu}$ is substituted by $t^{\prime}/\tau_{\nu}$, as defined in Eq.
(32). It is seen that the eddy viscous modeling only provides an adequate
description of $\bar{v}_{\theta}(r,t)$ for $r\leq R$, which is just the core
region of the filtered instantons. Notwithstanding such a limitation, extreme
circulations events are in fact well accounted by eddy viscous modeling,
since, as the filtering results indicate, $\bar{v}_{\theta}(r,0)$ gets its
maximum value at $r=R$.

To further investigate the detailed structure of the circulation instantons
moving beyond the viscous solutions, one must deal with the full nonlinear
Euler-Lagrange equations (9-11), an issue we address in the next section.

##  III nonlinear Instantons

A numerical scheme to solve the hydrodynamic instanton equations, similar to
Eqs. (8-9), was introduced by Chernykh and Stepanov chernykh2001large  in the
context of Burgers turbulence, being later extended to a variety of models
dematteis2018rogue ; dematteis2019experimental . In our case, the solution
algorithm inserts, at “iteration step $n$”, a given velocity field in Eq. (9)
to solve it backward in time. Its solution is defined as the auxiliary field
$\vec{p}(\vec{r},t)$ at iteration step $n$. This, in turn, is used in Eq. (8)
to provide the velocity field at iteration step $n+1$. An iteration cycle
proceed recursively, until a prescribed convergence for the instanton fields
is attained. The numerical procedure may start off with a vanishing velocity
field or one may alternatively use some improved guess, as we do by using the
creeping instanton described by (18).

We have worked out solutions for the circulation instantons in cylindrical
coordinates with Chebyshev collocation points (see the appendix for details)
for a broad range of $\lambda$ values and contour radii. The numerical
solutions were obtained through the Chernykh-Stepanov method outlined above,
enhanced, at larger values of $\lambda$, by the convexification procedure
discussed in alqahtani2021instantons , in order to avoid possible spoiling
inflection effects in the derivation of the cPDF tails. This amounts, in
practical terms, to perform the replacement of $\Gamma$ by the nonlinear
tilted measure ${\hbox{sign}}(\Gamma)\ln|\Gamma|$ in the MSRJD action (II.1).

Figure 4: Vector streamlines of a circulation instanton at its final evolution
time on the $(r,z)$ plane. Values of $\sqrt{v_{r}^{2}(r,z)+v^{2}_{z}(r,z)}$
are qualitatively represented by the gray level of the streamlines, while the
intensity of the azimuthal velocity $v_{\theta}(r,z)$ by a (red-colored)
density plot. Values vary from 0 (white) to approximately $11.9$ (black) and
$79.3$ (red), in units of $\nu/L$. This flow configuration is obtained as a
solution of Eqs. (8-11) for $\lambda=3\times 10^{6}$ and $R/L\approx 0.32$.
The vertical dashed lines give the radial positions $r=R$ and $r=L$.

Instability issues are known to affect the convergence of the Chernykh-
Stepanov solutions at large values of the Lagrange multiplier $\lambda$
alqahtani2021instantons . Although a number of technical improvements have
been implemented to solve instanton equations in more efficient ways
grafke2014arclength ; grafke2019numerical ; grigorio2020parametric , three-
dimensional instanton equations are far more computationally expensive than
their dimensionally reduced counterparts.

Fig. 4 yields the typical velocity profile of a circulation instanton. A
threefold vortex structure emerges, with interesting topological properties.
It consists of two paired counter-rotating vortex rings, as they can be
clearly identified from the streamline portrait there depicted in the $(r,z)$
plane. These rings define regions of opposite helicities, a fact associated to
the background axisymmetric velocity field that circulates around the symmetry
axis of the flow.

Figure 5: The instanton profiles of $v_{\theta}(r,z=0)$ at the final evolution
time for several values of $\lambda$, and $R/L\approx 0.32$. The vertical
dashed lines give the radial positions $r=R$ and $r=L$, while the inverted
triangles locate the peak positions of $v_{\theta}(r,z=0)$.

It is not difficult to see, from Fig. 4, that the azimuthal component of the
velocity field is not peaked at $r=R$ as we have found in our analysis of the
filtered DNS instantons (extreme circulation events), previously summarized in
Fig. 3. This is not a completely casual remark. In fact, it turns out that the
coupled instanton equations (8-11) become very stiff at large values of
$\lambda$, precisely where we expect to model asymptotically large circulation
fluctuations.

In Fig. 5 we show that, as $\lambda$ grows, the peak position of
$v_{\theta}(r,z=0)$ slowly drifts to the left (smaller values of $r$). A rough
estimate indicates that a peak at $r=R$ would be reached for $\lambda\sim
10^{10}$, which is several orders of magnitude beyond the domain of
convergence achieved in our applications of the Chernykh-Stepanov method.

We find, thus, that direct numerical schemes to solve the instanton equations
(8-11), like the ones we have applied, are able to address the parabolic cores
of cPDFs, while the description of their tails remains challenging. In order
to circumvent the stiffness of the instanton equations, and to model the far
cPDF tails, space-time reparametrizations chernykh2001large ;
grafke2014arclength ; grigorio2020parametric  and further large deviation
techniques are likely to be necessary, as illustrated by the use of hybrid
Monte Carlo algorithms, so far only applied to Burgers turbulence
margazoglou2019hybrid .

##  IV Conclusions

We have investigated the occurrence of extreme circulation events in three-
dimensional homogeneous and isotropic turbulence. We show that the time
evolution of circulation, conditionally averaged to the observation of
prescribed large deviation events can be effectively modeled with the help of
eddy-viscosity phenomenology combined with instanton functional techniques.
Our modeling analysis is corroborated from a careful treatment of large DNS
databases, at various Reynolds numbers.

Proceeding with the same set of extreme circulation events, we have
furthermore inspected time-dependent conditionally averaged velocity field
configurations. Resulting axisymmetric vortex structures clearly arise from
these filtering procedures. The comparison with the eddy-viscous (creeping)
instantons is only reasonable at the core of the flow configurations
(circulation instantons), which, nevertheless, is all one needs to get
accurate averages of the time-dependent instanton circulations ending in
extreme events.

An extensive numerical effort to solve the instanton equations suggests that
extreme circulation events have an interesting underlying topological
structure. Axisymmetric circulation instantons are composed of a main vortex,
which is at the center of a surrounding pair of counter-rotating vortex rings.
The existence of this triple vortex structure drives our attention to possible
formal connections with superfluidity. It is known that the motion of a
superfluid vortex ring is followed by a companion pair of normal vortex rings
kivotides2000triple . We note, then, that the instanton Eqs. (8) and (9) are
actually analogous to the HVBK hydrodynamic equations describing the self-
induced propagation of a quantum vortex ring.

Although the employed numerical databases allowed us to validate a viscous-
eddy modeling of extreme circulation events, they are not large enough to
disclose the topological features of the circulation instantons. Larger
databases are required, in order to accomplish this challenging task. Another
related important issue, deserved for further studies, has to do with the
implementation of algorithmic improvements in the numerical solution of the
instanton equations, in order to model the far cPDF tails iyer2019circulation
; iyer2021area  and the off-core azimuthal velocity component of the
circulation instantons. It is worth emphasizing that some alternatives have
been already applied with success to the paradigmatic example of Burgers
turbulence grafke2014arclength ; grafke2015relevance .

Intense vortex structure is known to have an important role in turbulence at
high Reynolds numbers. Such structures can be modeled by an interacting
molecular gas when statistical properties of circulation in 2D slices of
turbulent flows are analyzed apolinario2020vortex ;
moriconi2021multifractality ; moriconi2022statistics ; pereira2022hard . A
possible source of future investigation is to explore a single or a bunch of
modeled structures in order to compare its spatial distribution to the
instanton solutions.

## Acknowledgement

This work has been partially supported by CAPES via Grant No.
88887.336246/2019-00 (V.J.V.) and by the Simons Foundation Award ID 651475
(G.B.A.).

## References

  * (1) M.J. Aschwanden, ApJ 880, 105 (2019).
  * (2) C. Larrodera, L. Nikitina, and C. Cid, Space Weather 19, e2021SW002902 (2021).
  * (3) F. Bouchet, C. Nardini, and T. Tangarife , Fluid Dyn. Res. 46, 061416 (2014).
  * (4) J. Laurie and F. Bouchet, New. J. Phys. 17, 015009 (2015).
  * (5) A. Couto, P. Costa, and T. Simões, Energies 14, 3944 (2021).
  * (6) M. Rocco, J. Econ. Surv. 28, 82 (2014).
  * (7) P. Yiou and M. Nogaj, Geophys. Res. Lett. 31, (2004).
  * (8) R. Grotjahn, R. Black, R. Leung, M.F. Wehner, M. Barlow, M. Bosilovich, A. Gershunov, W.J. Gutowski, J.R. Gyakum, R.W. Katz, et al, Clim. Dyn. 46, 1151 (2016).
  * (9) A. Soloviev, Nonlinear Process Geophys 15, 209 (2008).
  * (10) K. Yoshida and N. Hatekeyama, J. Phys. Soc. Japan 69, 1661 (2000).
  * (11) Y. Kaneda, T. Ishihara, M. Yokokawa, K. Itakura, and A. Uno, Phys. Fluids 15, L21 (2003).
  * (12) U. Frisch, Turbulence, (New York, Cambridge University Press, 1995).
  * (13) G.K. Batchelor and A.A. Townsend, Proc. R. Soc. A: Math. Phys. Sci. 199, 238 (1949).
  * (14) F. Anselmet, Y. Gagne, E.J. Hopfinger, and R.A. Antonia, J. Fluid Mech. 140, 63 (1984).
  * (15) A. Vincent and M. Meneguzzi, J. Fluid Mech. 225, 1 (1991).
  * (16) A.M. Obukhov, J. Fluid Mech. 13, 77 (1962).
  * (17) A.N. Kolmogorov, J. Fluid Mech. 13, 82 (1962).
  * (18) E. Novikov and R. Stewart, Izv. Geophys. Ser 3, 408 (1964).
  * (19) Z.S. She and E. Leveque, Phys. Rev. Lett. 72, 336 (1994).
  * (20) G.B. Apolinário, L. Moriconi, R.M. Pereira and V.J. Valadão, Phys. Rev. E 102, 041102 (2020).
  * (21) L. Moriconi, R.M. Pereira, and V.J. Valadão, arXiv:2205.10470.
  * (22) J.M. Burgers, Adv. Appl. Mech. 1, 171 (1948).
  * (23) V. Gurarie and A. Migdal, Phys. Rev. E 54, 4908 (1996).
  * (24) P.C. Martin, E.D. Siggia, and H.A. Rose, Phys. Rev. A 8, 423 (1973).
  * (25) H.K. Janssen, Z. Phys. B: Cond. Mat. 23, 377 (1976).
  * (26) C. de Dominicis, J. Phys. Colloque 37, C1–247 (1976).
  * (27) E. Balkovsky, G. Falkovich, I. Kolokolov, and V. Lebedev, Phys. Rev. Lett. 78, 1452 (1997).
  * (28) A.I. Chernykh and M.G. Stepanov, Phys. Rev. E 64, 263061 (2001).
  * (29) L. Moriconi, Phys. Rev. E 79, 046324 (2009).
  * (30) T. Grafke, R. Grauer, T. Schäfer, and E. Vanden-Eijnden, Eur. Phys. Lett. 109, 34003 (2015).
  * (31) G.B. Apolinário, L. Moriconi, and R.M. Pereira, Phys. Rev. E 99, 033104 (2019).
  * (32) P.L. Johnson and C. Meneveau, J. Fluid Mech. 804, 387 (2016).
  * (33) L.S. Grigorio, F. Bouchet, R.M. Pereira, and L. Chevillard, J. Phys. A: Math. Theor. 50, 055501 (2017).
  * (34) G.B. Apolinário, L. Moriconi, and R.M. Pereira, Phys. A: Stat. Mech. Appl. 514, 741 (2019).
  * (35) M. Alqahtani, L. Grigorio, and T. Grafke, arXiv:2108.02103.
  * (36) M. Umeki, J. Phys. Soc. Japan 62, 3788 (1993).
  * (37) A. Migdal, Int J Mod Phys A 9, 1197 (1994).
  * (38) N. Cao, S. Chen, and K.R. Sreenivasan, Phys. Rev. Lett. 76, 616 (1996).
  * (39) R. Benzi, L. Biferale, M.V. Struglia, and R. Tripiccione, Phys. Rev. E 55, 3739 (1997).
  * (40) K.P. Iyer, K.R. Sreenivasan, and P.K. Yeung, Phys. Rev. X 9, 041006 (2019).
  * (41) K.P. Iyer, S.S. Bharadwaj, and K.R. Sreenivasan, Proc. Natl. Acad. Sci. U.S.A. 118, e2114679118 (2021).
  * (42) N.P. Müller, J.I. Polanco, and G. Krstulovic, Phys. Rev. X 11, 011053 (2021).
  * (43) J.I. Polanco, N.P. Müller, and G. Krstulovic, Nat. Commun. 12, 7090 (2021).
  * (44) L. Moriconi, Phys. Rev. E 103, 062137 (2021).
  * (45) L. Moriconi and R.M. Pereira, arXiv:2204.12862.
  * (46) P.A. Durbin and G. Medic, Fluid Dynamics with a Computational Perspective, (Cambridge, Cambridge University Press, 2010).
  * (47) J. Boussinesq, Mem. Pres. Par Diu. Sav. Acad. Sci. Paris 23, 46 (1877).
  * (48) A. Weiss, IEEE JSAC 13, 938 (1995).
  * (49) C. Evertsz and K. Berkner, Chaos Solit. Fractals 6, 121 (1995).
  * (50) F. Ragone, J. Wouters, and F. Bouchet, Proc. Natl. Acad. Sci. U.S.A. 115, 24 (2018).
  * (51) L. Brenig, M. Malek-Monsour, and W. Horsthemke, Phys. Lett. A 59, 341 (1976).
  * (52) E. Ivashkevich, J. Phys. A: Math. Gen. 30, L525 (1997).
  * (53) T. Grafke, R. Grauer, and T. Schäfer, J. Phys. A: Math. Theor. 48, 333001 (2015).
  * (54) A. Dembo and A.D.O Zeitouni, Large deviations techniques and applications, (Berlin, Springer, 2010).
  * (55) M. Alqahtani and T. Grafke, J. Phys. A: Math. Theor. 54, 175001 (2021).
  * (56) L. Moriconi and F. Takakura, Phys. Rev. E 58, 3187 (1998).
  * (57) L. Moriconi and F. Nobre, Phys. Rev. E 65, 036302 (2002).
  * (58) E. Perlman, R. Burns, Y. Li, and C. Meneveau, Proc. ACM/IEEE, 1 (2007).
  * (59) P.K. Yeung, D.A. Donzis, and K.R. Sreenivasan, J. Fluid Mech. 700, 5 (2012).
  * (60) P.K. Yeung, X.M. Zhai, and K.R. Sreenivasan, Proc. Natl. Acad. Sci. U.S.A. 112, 12633-12638 (2015).
  * (61) P.K. Yeung, K.R. Sreenivasan, and S.B. Pope, Phys. Rev. Fluids 3, 064603 (2018).
  * (62) L.M. Smith and S.L. Woodruff, Annu. Rev. Fluid Mech. 30, 275 (1998).
  * (63) R.H. Kraichnan, Phys. Fluids 7, 1030 (1964).
  * (64) D. Forster, D.R. Nelson, and M.J. Stephen, Phys. Rev. A 16, 732 (1977).
  * (65) V. Yakhot and S.A. Orszag, J. Sci. Comput. 1, 3 (1986).
  * (66) J. Smagorinsky, Mon. Weather Rev. 91, 99 (1963).
  * (67) A. Scotti, C. Meneveau, and D.K. Lilly, Phys. Fluids A: Fluid Dyn. 5, 2306 (1993).
  * (68) M. Germano, U. Piomelli, P. Moin, and W.H. Cabot, Phys. Fluids A: Fluid Dyn. 3, 1760 (1991).
  * (69) P.J. Mason, Q. J. R. Meteorol. Soc. 120, 1 (1994).
  * (70) U. Piomelli, Prog. Aerosp. Sci. 35, 335 (1999).
  * (71) U. Piomelli, Philos. Trans., Math. Phys. Eng. Sci. 372, 20130320 (2014).
  * (72) T. Grafke, R. Grauer, and T. Schäfer, J. Phys. A: Math. Theor. 46, 062002 (2013).
  * (73) G. Dematteis, T. Grafke,and E. Vanden-Eijnden, Proc. Natl. Acad. Sci. U.S.A. 115, 855 (2018).
  * (74) G. Dematteis, T. Grafke, M. Onorato, and E. Vanden-Eijnden, Phys. Rev. X 9, 041057 (2019).
  * (75) T. Grafke, R. Grauer, T. Schäfer, and E. Vanden-Eijnden, Multiscale Model. Simul. 12, 566 (2014).
  * (76) T. Grafke and E. Vanden-Eijnden, Chaos 29, 063118 (2019).
  * (77) L. Grigorio, J. Phys. A: Math. Theor. 53, 445001 (2020).
  * (78) G. Margazoglou, L. Biferale, R. Grauer, K. Jansen, D. Mesterházy, T. Rosenow, and R. Tripiccione, Phys. Rev. E 99, 053303 (2019).
  * (79) D. Kivotides, C.F. Barenghi, and D.C. Samuels, Science 290, 777 (2000).
  * (80) C. Canuto, M. Hussaini, A. Quarteroni, and T. Zang, Spectral methods in fluid dynamics, (Berlin, Springer, 1988).
  * (81) R. Peyret, Spectral methods for incompressible viscous flow, (Berlin, Springer, 2013).
  * (82) K.J. Arrow, L. Hurwicz and, H. Uzawa, Studies in linear and nonlinear programming, (California, Stanford University Press, 1958).
  * (83) I. Raspo, S. Hugues, E. Serre, A. Randriamampianina, and P. Bontoux, Comput. Fluids 31, 745 (2002).
  * (84) J. Lopez and J. Shen, J. Comput. Phys. 139, 308 (1998).
  * (85) N. Peres, S. Poncet, and E. Serre, J. Comput. Phys. 231, 6290 (2012).
  * (86) T. Schorlepp, T. Grafke, S. May, and R. Grauer, arXiv:2107.06153.

*

##  Appendix A Spectral Approach to the Instanton Equations in Cylindrical
Coordinates

In order to numerically solve the set of Eqs. (8-11) we rewrite them in
cylindrical coordinates. In light of the nonlinear nature of hydrodynamic
interactions, there would be no a priori reason to consider axisymmetric
solutions in detriment of any other. Nevertheless, we expect axisymmetric
solutions to dominate extreme events since the forcing term is invariant under
rotations around the $z$-axis. Taking this symmetry into account, Eqs. (8-11)
read

$\begin{cases}\hat{\mathcal{L}}^{-}v_{r}+v_{r}/r^{2}-v_{\theta}^{2}/r+\partial_{r}P=(\tilde{\chi}\star
p_{r})\ ,\\\
\hat{\mathcal{L}}^{-}v_{\theta}+v_{\theta}/r^{2}+v_{\theta}v_{r}/r=(\tilde{\chi}\star
p_{\theta})\ ,\\\ \hat{\mathcal{L}}^{-}v_{z}+\partial_{z}P=(\tilde{\chi}\star
p_{z})\ ,\\\ \partial_{r}(rv_{r})/r+\partial_{z}v_{z}=0\ ,\\\
\hat{\mathcal{L}}^{+}p_{r}-p_{r}/r^{2}-v_{\theta}p_{\theta}/r+v_{r}\partial_{r}p_{r}+v_{\theta}\partial_{r}p_{\theta}+v_{z}\partial_{r}p_{z}+\partial_{r}Q=0\
,\\\
\hat{\mathcal{L}}^{+}p_{\theta}-p_{\theta}/r^{2}+(2v_{\theta}p_{r}-v_{r}p_{\theta})/r=\lambda\delta(t)\delta(r-R)\delta(z)\
,\\\
\hat{\mathcal{L}}^{+}p_{z}+v_{r}\partial_{z}p_{r}+v_{\theta}\partial_{z}p_{\theta}+v_{z}\partial_{z}p_{z}+\partial_{z}Q=0\
,\\\ \partial_{r}(rp_{r})/r+\partial_{z}p_{z}=0\ ,\ \end{cases}$ (40)

where
$\hat{\mathcal{L}}^{\pm}\psi=\partial_{t}\psi+v_{z}\partial_{z}\psi+v_{r}\partial_{r}\psi\pm(\partial_{z}^{2}\psi+\partial_{r}(r\partial_{r}\psi)/r)$
and $(\phi\star\psi)$ is the spatial convolution between $\phi$ and $\psi$. We
implement a pseudo-spectral approach where the fields are expanded in a
truncated Fourier-Chebyshev series of orders $N$ and $M$, respectively. The
collocation points are defined as

$z_{j}=\frac{2\pi j}{N}\ ,\ $ (41)

and

$r_{m}=\cos{\left(\frac{\pi m}{M-1}\right)}\ ,\ $ (42)

where, $j$ and $m$ are integers in the ranges $j\in[-N/2,N/2-1]$ and
$m\in[0,M-1]$. Along the longitudinal direction, usual collocation points for
Fourier series are used, while for the radial direction we employ Gauss-
Lobatto collocation points, corresponding to the extrema of the $M^{th}$ order
Chebyshev polynomials canuto1988spectral . This choice of a non-homogeneous
grid hinders a direct comparison among simulations with different resolutions.
However, it allows one to adopt a Fast Fourier Transform (FFT) algorithm to
compute nonlinear terms, which are otherwise computationally expensive
peyret2013spectral . Moreover, for even $M$ the coordinate singularity at
$r=0$ is explicitly avoided on the collocation, with the drawback of requiring
the radial inverval to be duplicated, as $r\in[0,1]\rightarrow r\in[-1,1]$, as
a requirement of the Chebyshev polynomial expansion.

All fieds in Eq. (40) are given by a truncated polynomial series of order $N$
in the Fourier basis and $M$ in the Chebyshev basis,

$\phi(r_{m},z_{j})\equiv\phi_{mj}=\sum_{l=0}^{M-1}\sum_{k=-N/2}^{N/2-1}\hat{\phi}_{lk}T_{l}(r_{m})e^{ikz_{j}}\
,\ $ (43)

where $T_{l}(x)$ is the $l^{th}$ Chebyshev polynomial. Exceptions are the
pressure fields $Q(\vec{x},t)$ and $P(\vec{x},t)$, which are of order $M-2$ in
the Chebyshev expansion. This guarantees the stability of the numerical
method, an approach known as the $\mathbb{P}_{N}-\mathbb{P}_{N-2}$
approximation peyret2013spectral , further discussed ahead. Radial derivatives
of these fields are readily computed using recursion relations of the
Chebyshev polynomials, namely,

$(\partial_{r}\phi)_{mj}=\sum_{l=0}^{M-1}D_{ml}\phi_{lj}\ ,\ $ (44)

where the $M\times M$ matrix $D$ has, for all fields but the pressure related
ones, the entries

$D_{mj}=\begin{cases}(-1)^{m+j}c_{m}/(c_{j}(r_{m}-r_{j}))&\ ,\ m\neq j\\\
-r_{m}/(2(1-r_{m}^{2}))&\ ,\ m=j,\ m,j\neq[0,(M-1)]\\\ \pm(2(N-1)^{2}+1)/6&\
,\ m=j=0\text{ or }m=j=(M-1),\end{cases}$ (45)

with $c_{0}=c_{M-1}=2$ and $c_{i}=1$ otherwise. As for the pressure fields,
which are 2 degrees lower in the Chebyshev expansion, they must be
interpolated from the $M-2$ Gauss-Lobatto grid to the $M$ grid where the other
fields are defined. This can be done through a Lagrange interpolation (see
Ref. peyret2013spectral  for details), giving

$D^{PQ}_{mj}=\begin{cases}(-1)^{j+m}(1-r_{j}^{2})/((r_{m}-r_{j})(1-r_{m}^{2}))&\
,\ m\neq j,\ m,j\in[1,M-2]\\\ 3r_{m}/(2(1-r_{m}^{2}))&\ ,\ m=j,\
m\in[1,M-2].\end{cases}$ (46)

By virtue of the explicit form of these matrices, $r$-derivatives are
evaluated in physical space while $z$-derivatives are better performed in
Fourier space. Ultimately, in order to have consistently represent all fields
in the mirrored radial domain, they must satisfy the following reflection
properties:

$\begin{cases}&\phi_{mj}=\phi_{(-m)j}\text{ for scalar fields,}\\\
&\phi_{z,mj}=\phi_{z,(-m)j}\text{ for the }z\text{-component of vector
fields}\\\ &\phi_{\beta,mj}=-\phi_{\beta,(-m)j}\text{ for the radial and
azimuthal components (}\beta\text{ is either }r\text{ or }\theta).\end{cases}$
(47)

The time discretization follows a combined Adams-Bashforth/Implicit Backward
Differentiation method of second order (AB/BDI2) peyret2013spectral , which
consists on implicit evaluations of the linear terms and explicit evaluations
of the nonlinear terms. As an illustration, consider the equation

$\partial_{t}\phi=L(\phi)+N(\phi)+J\ ,\ $ (48)

where, $L(\phi)$ and $N(\phi)$ stand for the linear and nonlinear terms of the
differential equation, respectively, and $J$ accounts for a forcing term or
pressure gradient. The above equation is discretized at regularly spaced time
instants $t_{n}=ndt$ as

$\frac{3\phi^{(n+1)}-4\phi^{(n)}+\phi^{(n-1)}}{2dt}=L\phi^{(n+1)}+2N(\phi^{(n)})-N(\phi^{(n-1)})+J^{(n+1)}\
,\ $ (49)

where $\phi^{(n)}=\phi(t_{n})$. In the first time step, we set
$\phi^{(-1)}=\phi^{(0)}$ and change $dt\rightarrow 3dt/2$, reducing this step
to the usual Euler scheme. The choice of discretizing the pressure gradient
terms as $\partial_{i}P^{(n+1)}$ leads to a Stokes problem which is solvable
by the Uzawa method uzawa1958iterative . We note that similar discretization
setups were successfully applied to the Navier-Stokes equations in cylindrical
coordinates with a few different boundary conditions raspo2002spectral ;
lopez1998efficient ; peres20123d . We also remark that the nonlinear terms are
properly de-aliased following a standard 3/2-rule.

Applying the above discretization procedure to the system of Eqs. (40)
supplemented by Dirichlet boundary conditions, one finds a system of coupled
equations for each independent Fourier mode $k\in[N/2,N/2-1]$ for the forward
and backward time integration of the velocity and conjugate field,
respectively,

$\begin{cases}\hat{L}_{r,k}\ket{v_{r,k}^{n}}=-\hat{D}^{PQ}\ket{P_{k}^{n}}+\ket{f_{r,k}^{n}}&\
,\ j\in[1,M-2]\ ,\\\
\hat{L}_{r,k}\ket{v_{\theta,k}^{n}}=\ket{f_{\theta,k}^{n}}&\ ,\ j\in[1,M-2]\
,\\\ \hat{L}_{z,k}\ket{v_{z,k}^{n}}=-ik\ket{P_{k}^{n}}+\ket{f_{z,k}^{n}}&\ ,\
j\in[1,M-2]\ ,\\\
\hat{R}^{-1}\hat{D}\hat{R}\ket{v_{r,k}^{n}}+ik\ket{v_{z,k}^{n}}=0&\ ,\
j\in[0,M-1]\ ,\\\
(\ket{v_{r,k}^{n}},\ket{v_{\theta,k}^{n}},\ket{v_{z,k}^{n}})=0&\ ,\ j=0\text{
and }j=M-1\ ,\\\ \\\
\hat{L}_{r,k}\ket{p_{r,k}^{n}}=\hat{D}^{PQ}\ket{Q_{k}^{n}}+\ket{g_{r,k}^{n}}&\
,\ j\in[1,M-2]\ ,\\\
\hat{L}_{r,k}\ket{p_{\theta,k}^{n}}=\ket{g_{\theta,k}^{n}}&\ ,\ j\in[1,M-2]\
,\\\ \hat{L}_{z,k}\ket{p_{z,k}^{n}}=ik\ket{Q_{k}^{n}}+\ket{g_{z,k}^{n}}&\ ,\
j\in[1,M-2]\ ,\\\
\hat{R}^{-1}\hat{D}\hat{R}\ket{p_{r,k}^{n}}+ik\ket{p_{z,k}^{n}}=0&\ ,\
j\in[0,M-1]\ ,\\\
(\ket{p_{r,k}^{n}},\ket{p_{\theta,k}^{n}},\ket{p_{z,k}^{n}})=0&\ ,\ j=0\text{
and }j=M-1\ ,\\\ \end{cases}$ (50)

where the ket notation represents the $N$-dimensional vector
$\ket{\phi_{i,k}^{n}}=[\phi_{i,jk}^{n}]$ with $j\in[0,M-1]$. The $M\times M$
matrices $\hat{R}=\mathrm{diag}(r_{0},r_{1},\cdots,r_{M-1})$,
$\hat{L}_{r,k}=(3/(2|dt|)+k^{2})\hat{\mathbb{I}}-\hat{D}^{2}-\hat{R}^{-1}\hat{D}+\hat{R}^{-2}$
and $\hat{L}_{z,k}=\hat{L}_{r,k}-\hat{R}^{-2}$ can be efficiently inverted and
stored in a pre-processing stage. The $M$ vectors
$(\ket{f_{r,k}^{n}},\ket{f_{\theta,k}^{n}},\ket{f_{z,k}^{n}})$ and
$(\ket{g_{r,k}^{n}},\ket{g_{\theta,k}^{n}},\ket{g_{z,k}^{n}})$ are,
respectively, the explicit part of the discretized Eqs. 8 and 9. For instance,
in Eq. 49 one has
$\ket{f_{k}^{n}}=(4\phi^{(n-1)}-\phi^{(n-2)}+2N(\phi^{(n-1)})-N(\phi^{(n-2)}))/3$.

The algebraic system defined by Eq. 50 has a unique solution for every $k\neq
0$. The $\mathbb{P}_{N}-\mathbb{P}_{N-2}$ approximation prevents zero
eigenvalues of the Uzawa operator333The Uzawa operator is obtained by solving
formally the pressure field by setting the momentum equations into the
incompressibility constraint, in this case
$\hat{Z}_{k}=\hat{R}^{-1}\hat{D}\hat{R}\hat{L}_{r}^{-1}\hat{D}^{PQ}-k^{2}\hat{L}_{z}^{-1}$.
for $k=0$ and avoids the requirement of prescribing boundary conditions to the
pressure field. Indeed, non-uniqueness of the solution is related to the fact
that pressure fields are defined up to a constant and a simple calculation
shows that the unique solution consistent with the boundary conditions and the
incompressibility constraint for $k=0$ is $\ket{v_{r,0}^{n}}=0$.

The convolutions in Eq. (40) and in the action integral (Eq. II.1) can be
efficiently computed using the explicit form of the Fourier transformed
correlation $\tilde{\chi}$ and the inverted Chebyshev derivative matrix,

$\mathrm{FFT}_{z}\big{[}(\tilde{\chi}\star
p_{\beta})\big{]}_{jk}^{(n)}=(2\pi)^{3/2}Le^{-\frac{r_{j}^{2}}{2L^{2}}}e^{-\frac{k^{2}L^{2}}{2}}\sum_{l}\left(\hat{D}^{-1}_{M/2,l}-\hat{D}^{-1}_{0,l}\right)U_{jl,\beta}f^{(n)}_{l,\beta}\
,\ $ (51)
$S=(2\pi)^{7/2}L\sum_{n,k,l,j,\beta}e^{-\frac{k^{2}L^{2}}{2}}\left(\hat{D}^{-1}_{M/2,l}-\hat{D}^{-1}_{0,l}\right)\left(\hat{D}^{-1}_{M/2,j}-\hat{D}^{-1}_{0,j}\right)f^{(n)}_{lk,\beta}U_{jl,\beta}\left(f^{(n)}_{jk,\beta}\right)^{\star}\
,\ $ (52)

where $\beta=(r,\theta,z)$,
$f^{(n)}_{lk,\beta}=r_{l}\exp{(-r_{l}^{2}/2L^{2})}FFT_{z}[p_{\beta}]^{(n)}_{lk}$,
$U_{jk,\beta}=I_{0}(r_{j}r_{k}/L^{2})$ for $\beta=z$, and
$U_{jk,\beta}=I_{1}(r_{j}r_{k}/L^{2})$ for $\beta=(r,\theta)$, with
$I_{0},I_{1}$ being modified Bessel functions of the first kind.

In order to validate the numerical method, we performed several numerical
experiments of the linear instanton. Fig. 6 compares the $v_{\theta}$
component of the numerical solution with the analytical result from Eq. 18. We
first note that the Dirichlet boundary condition has a relevant influence on
the solution for $r\gtrsim 0.5$. This is not a surprising though, since Eq. 18
holds for unbounded domains, thus a slower decay is expected. Variations of
about $\pm 5\%$ in the peak position are also observed depending on how
distant to the boundaries one sets $L$ and $R$. We fixed $L=1/\pi$ to minimize
such boundary effects.

Figure 6: Difference between normalized numerical and analytical solutions of
the azimuthal velocity component at $t=0$ as functions of $r$ and $z$, for
$R\approx 0.09$ and $L=1/\pi$. Inset: profiles for $z=0$.

The viscous action was calculated by means of Eq. 52, and it was found to be
compatible with $S=\Gamma^{\alpha_{\Gamma}}/(2\sigma^{2})$ with
$\alpha_{\Gamma}=1.99999993(5)$ and $\sigma=(\beta_{\sigma}/N)R^{1.95(2)}$,
with errors estimated by averaging results obtained with different grid
resolutions, both in time and space. Small variations of $\beta_{\sigma}$ are
seen when $dt$ and/or $M$ are changed, but one must keep in mind that the
radial collocation is not regular, and hence direct comparisons among
different resolutions in $M$ are not perfect. As for the $dt$ dependence,
finite time effects are in play, since decreasing $dt$ for a fixed number of
timesteps $N_{t}$ also causes the total simulation time to decrease, so the
boundary conditions $p_{i}(\vec{x},T)=v_{i}(\vec{x},T)=0$ are effectively
imposed on different time instants $T=-N_{t}dt$. In out tests, we worked with
all combinations of $N_{t}=100$, $200$, and $400$, with $dt=0.0025$, $0.005$,
and $0.01$.

The conclusion drawn from this set of numerical experiments is that both
spatial and statistical properties are accurately captured by the numerical
algorithm used to solve the instanton equations, at least in the linear
approximation.


Gleb I. PosobinJohn Q. Open and Joan R. Acces 2 42nd Conference on Very
Important Topics (CVIT 2016) CVIT 2016 CVIT 2016 December 24–27, 2016 Little
Whinging, United Kingdom 42 23

# Computing majority with low-fan-in majority queries111This work was
partially supported by the Russian Academic Excellence Project ’5-100’.

Gleb I. Posobin  National Research University, Higher School of Economics,
Moscow, Russia
<EMAIL_ADDRESS>
###### Abstract.

In this paper we examine the problem of computing majority function $MAJ_{n}$
on $n$ bits by depth-two formula, where each gate is a majority function on at
most $k$ inputs. We present such formula that gives the first nontrivial upper
bound for this problem, with $k=\frac{2}{3}n+4$. This answers an open question
in [7].

We also look at this problem in adaptive setting — when we are allowed to
query for value of $MAJ_{k}$ on any subset, and wish to minimize the number of
such queries. We give a simple lower bound for this setting with $\lceil
n/k\rceil$ queries, and we present two algorithms for this model: the first
one makes $\approx 2\frac{n}{k}\log k$ queries in the case when we are limited
to the standard majority functions, and the second one makes $\frac{n}{k}\log
k$ queries when we are allowed to change the threshold of majority function.

###### Key words and phrases:

majority, upper bound, threshold, circuit complexity

###### 1991 Mathematics Subject Classification:

F.1.3 Complexity Measures and Classes

##  1\. Introduction

We study the problem of computing majority function $MAJ_{n}$ if we are only
allowed to query values of $MAJ_{k}$ (that is, majority on $k$ variables)
functions.

Majority functions and circuits consisting of them arise in various areas of
computational complexity (see e.g. [9, 3, 4]). Particularly, iterated (or
recursive) majority that consists of iterated application of majority of small
number of inputs to itself, provides an example of a function with interesting
complexity properties in various models ([5, 8, 10, 6]), and helps in various
constructions.

Another motivation for this problem comes from the studies of boolean circuits
of constant depth. $\widehat{\textsf{TC}}{}^{0}$ is defined as the class of
functions computable by constant-depth polynomial-sized circuits consisting of
majority gates, and it plays one of the major roles in this area ([4]).

The first model we look at is a standard boolean circuit of depth two, but
where all the gates are threshold functions of at most $k$ variables. It was
shown in [1] that to be able to compute such circuits $k$ must be
$\Omega(n^{4/5})$. There is a trivial upper bound $k=n$: just take all the
gates on the first level to match inputs (that is, gate
$f_{i}=x_{i}=MAJ_{\\{i\\}}(x)$), and the output gate is a standard majority
$MAJ_{n}(f_{1},\ldots,f_{n})$. [7] did not provide any nontrivial upper bound
on $k$ and posed an open question whether such bound exists. We answer this
question positively and present the depth-two circuit with
$MAJ_{\frac{2}{3}n+4}$ gates that computes $MAJ_{n}$ function. This result was
independently obtained by Bruno Bauwens by computer search, as told to us in
personal communication. For the sake of completeness, we also give a simple
proof of lower bound $k\geq\frac{1}{\sqrt[3]{2}}n^{2/3}$.

The second model is an adaptive one — algorithm can query the $MAJ_{k}$
function on any subset of size at most $k$, and we are interested in the worst
case number of queries made to compute $MAJ_{n}$. This problem was motivated
by the proof in [7], which does not use any properties of the function in the
output gate except for its monotonicity, and also by [2], where they study
adaptive computation of $MAJ_{n}$ with discrepancy queries on $k$ elements. We
give a simple proof for the lower bound of $\lceil n/k\rceil$ queries, and we
also present two algorithms for adaptive setting: one that runs in fixed-
threshold setting and which requires $2(\frac{n}{k-4}+1)(\log k+4)\approx
2\frac{n}{k}\log k$ queries, and another that runs in adjustable-threshold
adaptive model (that is, when the algorithm can specify majority’s threshold)
and requires $\lceil n/k\rceil\lceil\log(k+1)\rceil$ queries.

##  2\. Definitions

By $\log$ we denote binary logarithm and $n$ always means the length of input
vector $x=(x_{1},\ldots,x_{n})\in\\{0,1\\}^{n}$.

We are going to study majority functions
$MAJ_{n}(x_{1},\ldots,x_{n})=[x_{1}+\ldots+x_{n}\geq n/2]$, where $[P]=1$ if
condition $P$ is true, and $[P]=0$ otherwise. We will also call $MAJ_{k}$
standard majority functions. By $sum_{S}(x)$ we denote $\sum_{i\in S}x_{i}$.
We define $MAJ_{S}(x)$ for $S\subset\\{1,\ldots,n\\}$ as
$[sum_{S}(x)\geq|S|/2]$. We also define $MAJ_{S}(x;t)=[sum_{S}(x)\geq t]$.

Threshold functions are a generalization of majority functions:
$f(x_{1},\ldots,x_{n})=[\alpha_{1}x_{1}+\ldots+\alpha_{n}x_{n}\geq t]$.

The first (“static“) model that we look at is $MAJ_{k}\circ MAJ_{k}$ — it is
defined as depth-two formula, consisting of threshold gates
$[\sum\alpha_{i}x_{i}\geq t]$ such that $\sum\alpha_{i}\leq k$ and
$\alpha_{i}$ are positive integers.

In the second (“adaptive“) model an algorithm is given access to an oracle
that knows vector $x$. There are two variations of this model that we look at.

In the adaptive model with adjustable threshold the oracle receives a pair
$(S,t)$, where $S$ is a set $S\subset\\{1,\ldots,n\\}$ such that $|S|\leq k$
and $t$ is an integer, and responds with the value of $MAJ_{S}(x;t)$.

In the adaptive model with fixed threshold the oracle receives a set
$S\subset\\{1,\ldots,n\\}$ such that $|S|\leq k$ and responds with the value
of $MAJ_{S}(x)$.

##  3\. Static model

###  3.1. Lower bound

######  Lemma 3.1.

For large enough $n$, $k$ should be at least $\frac{1}{\sqrt[3]{2}}n^{2/3}$ in
order to compute $MAJ_{n}$ in $MAJ_{k}\circ MAJ_{k}$ model.

This is not the best currently known lower bound, for a more involved proof of
$k\geq n^{13/19+o(1)}$ see [7]. This proof is simpler and is adapted from the
proof of a stronger claim in [7] that shows that for some function in
$MAJ_{k}\circ MAJ_{k}$ to differ from $MAJ_{n}$ on small fraction (less than
$1/10$) of inputs $x$ with $\sum_{i}x_{i}=\lfloor(n-1)/2\rfloor$, $k$ must be
$\Omega(n^{2/3})$. We provide this proof just for the sake of completeness.

######  Proof 3.2.

Let’s look at the boolean hypercube $H^{n}=\\{0,1\\}^{n}$. It is known ([11])
that among all monotone boolean functions majority is the one with the largest
influence — that is, its value changes on the largest possible number of edges
of the hypercube. The number of such edges for $MAJ_{k}$ is
$k\binom{k}{k/2}\leq k^{1/2}2^{k}$. So if some gate $f_{i}$ depends on $l\leq
k$ variables, then it changes its value on at most $l^{1/2}2^{l}$ edges of
$H^{l}$.

Obviously, on each edge $(u,v)$ in $H^{n}$ such that
$\sum_{i}u_{i}<n/2,\sum_{i}v_{i}\geq n/2$, there must be at least one function
$f_{i}$ in the first level of the formula, such that $f_{i}(u)\not=f_{i}(v)$.
There are $n\binom{n}{n/2}>\frac{n^{1/2}2^{n}}{\sqrt{2}}$ such edges.

Let $f_{i}$ be some gate in the first level of the formula, without loss of
generality assume that $f_{i}$ depends only on $x_{1},\ldots,x_{l}$ for $l\leq
k$. Then, for each $u\in H^{l}$ there are at most $2^{n-l}/(n-l)^{1/2}$
possible continuations $u_{l+1},\ldots,u_{n}$ of $u$ that have $\sum_{1\leq
i\leq n}u_{i}=\lfloor(n-1)/2\rfloor$. This is so, since $u_{l+1},\ldots,u_{n}$
must sum to $\lfloor(n-1)/2\rfloor-\sum_{1\leq i\leq l}u_{i}$, and the number
of such assignments of $u$ is a binomial coefficient. The largest binomial
coefficient $\binom{n-l}{z}$ for fixed $n-l$ is the central one:
$z=\lceil(n-l)/2\rceil$, for which there is a well-known bound
$\binom{n-l}{z}\leq 2^{n-l}/(n-l)^{1/2}$.

This means that for a fixed $f_{i}$ there are at most
$l\frac{2^{l}}{l^{1/2}}\frac{2^{n-l}}{(n-l)^{1/2}}\leq\sqrt{2}k^{1/2}\frac{2^{n}}{n^{1/2}}$
such edges $(u,v)$ that $\sum_{i}u_{i}<n/2,\sum_{i}v_{i}\geq n/2$ and
$f_{i}(u)\not=f_{i}(v)$ (note that $(n-l)\geq\frac{1}{2}n$ for large enough
$n$, if $l\leq k=\mathcal{O}(n^{2/3})$). Summing this for all $f_{i}$ we get:

$\sqrt{2}k^{3/2}\frac{2^{n}}{n^{1/2}}\geq\frac{n^{1/2}2^{n}}{\sqrt{2}}\implies
k^{3}\geq\frac{n^{2}}{2}$

###  3.2. Upper bound

######  Theorem 3.3.

For any $n$ $MAJ_{n}$ is computable by some circuit in
$MAJ_{\frac{2}{3}n+4}\circ MAJ_{\frac{2}{3}n+4}$.

######  Proof 3.4.

Case 1. $n$ is divisible by $3$, $k=\frac{2}{3}n+1$.
We need to define $\frac{2}{3}n+1$ first-level threshold functions $f_{i}$ in
the formula. Let $f_{1}=x_{1},\ldots,f_{\frac{1}{3}n}=x_{\frac{1}{3}n}$, each
equivalently being a majority function on one corresponding element. Denote
$\frac{2}{3}n$ indices of yet unused inputs by
$M=\\{\frac{1}{3}n+1,\ldots,n\\}$. Then let
$f_{\frac{1}{3}n+i}=MAJ_{M}(x;\lceil\frac{1}{6}n\rceil+i-1)$. The output gate
is a standard majority on all $f_{i}$.

$\left[\sum
f_{i}\geq\frac{n}{3}+\frac{1}{2}\right]=\left[\Big{(}\sum_{\mathclap{i\leq
n/3}}x_{i}+\sum_{\mathclap{0\leq i\leq n/3}}[sum_{M}(x)\geq\lceil
n/6\rceil+i]\Big{)}\geq\frac{n}{3}+\frac{1}{2}\right]$

If
$sum_{M}(x)\geq\lceil\frac{1}{6}n\rceil+\frac{1}{3}n=\lceil\frac{1}{2}n\rceil$,
then $\sum_{0\leq
i\leq\frac{1}{3}n}[sum_{M}(x)\geq\lceil\frac{1}{6}n\rceil+i]=\frac{1}{3}n+1$,
that means that both output gate and $MAJ_{n}(x)$ will be equal to one.

If $\lceil\frac{1}{6}n\rceil\leq sum_{M}(x)\leq\lceil\frac{1}{2}n\rceil$, then
$\sum_{i\leq\frac{1}{3}n}x_{i}+\sum_{0\leq
i\leq\frac{1}{3}n}[sum_{M}\geq\lceil\frac{1}{6}n\rceil+i]=\sum_{i\leq\frac{1}{3}n}x_{i}+sum_{M}(x)-\lceil\frac{1}{6}n\rceil+1=\sum_{i}x_{i}-\lceil\frac{1}{6}n\rceil+1$.
Then the circuit attains value $1$ iff
$\sum_{i}x_{i}-\lceil\frac{1}{6}n\rceil+1\geq\frac{1}{3}n+\frac{1}{2}$.
Rearranging, we get that the value of the circuit is $1$ iff
$\sum_{i}x_{i}\geq\frac{1}{3}n+\lceil\frac{1}{6}n\rceil-\frac{1}{2}=\lceil\frac{1}{2}n\rceil-\frac{1}{2}$,
which is equivalent to $[\sum_{i}x_{i}\geq\frac{1}{2}n]=MAJ_{n}(x)$.

If $sum_{M}(x)<\lceil\frac{1}{6}n\rceil$, $MAJ_{n}(x)=0$ and the formula’s
value is also $0$, since $\sum_{i\leq\frac{1}{3}n}x_{i}\leq\frac{1}{3}n$.

Case 2. $n=3m+1$ for some integer $m$, $k=2m+3\leq\frac{2}{3}n+3$.
Let $f_{i}=x_{i}$ for $1\leq i\leq m+1$. Let $M=\\{m+2,\ldots,3m+1\\}$. Then
let $f_{m+1+i}=MAJ_{M}(x;\lceil\frac{1}{2}(m-1)\rceil+i-1)$ for $1\leq i\leq
m+2$, so there are $2m+3$ gates in the first level of the circuit. The output
gate is again a majority on all $f_{i}$.

If
$sum_{M}(x)\geq\lceil\frac{1}{2}(m-1)\rceil+m+1=\lceil\frac{3}{2}m+\frac{1}{2}\rceil$,
then both $MAJ_{n}(x)$ and the output gate are equal to one.

If $\lceil\frac{1}{2}(m-1)\rceil\leq
sum_{M}(x)\leq\lceil\frac{1}{2}(m-1)\rceil+m+1$, then
$\sum_{i}x_{i}=\sum_{i\leq
m+1}x_{i}+sum_{M}(x)-\lceil\frac{1}{2}(m-1)\rceil+1=\sum_{i}x_{i}-\lceil\frac{1}{2}(m-1)\rceil+1$.
Then the circuit attains value $1$ iff
$\sum_{i}x_{i}-\lceil\frac{1}{2}(m-1)\rceil+1\geq m+\frac{3}{2}$. Rearranging,
we get that the value of the circuit is $1$ iff $\sum_{i}x_{i}\geq
m+\lceil\frac{1}{2}(m-1)\rceil+\frac{1}{2}=\lceil\frac{3}{2}m+\frac{1}{2}\rceil-\frac{1}{2}$.
Rounding the last part to the smallest greater integer we get that the circuit
attains value $1$ iff $\sum_{i}x_{i}\geq\lceil\frac{3}{2}m+\frac{1}{2}\rceil$,
exactly when $MAJ_{n}(x)=1$.

If $sum_{M}(x)<\lceil\frac{1}{2}(m-1)\rceil$, then
$\sum_{i}x_{i}<m+1+\lceil\frac{1}{2}(m-1)\rceil=\lceil\frac{3}{2}m+1/2\rceil$,
so both $MAJ_{n}(x)$ and the circuit will have value $0$.

Case 3. $n=3m+2$ for some integer $m$, $k=2m+5\leq\frac{2}{3}n+4$.
Let $f_{i}=x_{i}$ for $1\leq i\leq m+2$. Let $M=\\{m+3,\ldots,3m+2\\}$. Then
let $f_{m+2+i}=MAJ_{M}(x;\lceil\frac{1}{2}m-1\rceil+i-1)$ for $1\leq i\leq
m+3$, so there are $2m+5$ gates in the first level of the circuit. The output
gate is again a majority on all $f_{i}$.

If
$sum_{M}(x)\geq\lceil\frac{1}{2}m-1\rceil+m+2=\lceil\frac{3}{2}m\rceil+1=\lceil\frac{3}{2}m+1\rceil$,
then both $MAJ_{n}(x)$ and the output gate are equal to one.

If $\lceil\frac{1}{2}m-1\rceil\leq
sum_{M}(x)\leq\lceil\frac{1}{2}m-1\rceil+m+2=\lceil\frac{3}{2}m+1\rceil$, then
$\sum_{i}x_{i}=\sum_{i\leq
m+2}x_{i}+sum_{M}(x)-\lceil\frac{1}{2}m-1\rceil+1=\sum_{i}x_{i}-\lceil\frac{1}{2}m\rceil+2$.
Then the circuit attains value $1$ iff
$\sum_{i}x_{i}-\lceil\frac{1}{2}m\rceil+2\geq m+\frac{5}{2}$. Rearranging, we
get that the value of the circuit is $1$ iff $\sum_{i}x_{i}\geq
m+\lceil\frac{1}{2}m\rceil+\frac{1}{2}=\lceil\frac{3}{2}m+1\rceil-\frac{1}{2}$.
Rounding the last part to the smallest greater integer we get that the circuit
attains value $1$ iff $\sum_{i}x_{i}\geq\lceil\frac{3}{2}m+1\rceil$, exactly
when $MAJ_{n}(x)=1$.

If $sum_{M}(x)<\lceil\frac{1}{2}m-1\rceil$, then
$\sum_{i}x_{i}<m+2+\lceil\frac{1}{2}m-1\rceil=\lceil\frac{3}{2}m+1\rceil$, so
both $MAJ_{n}(x)$ and the circuit will have value $0$.

##  4\. Adaptive model

###  4.1. Lower bound

######  Lemma 4.1.

Any algorithm for computing $MAJ_{n}$ in adaptive model requires $\lceil
n/k\rceil$ queries to oracle.

######  Proof 4.2.

We will set bits of vector $x$ the algorithm has not yet asked about right
before answering a query from the algorithm. Initialize a set $S=\varnothing$
— this will be the set of all the indices that appeared in algorithm’s
requests. When algorithm makes a request for $MAJ_{T}(x)$, we firstly set
$x_{i}$ for $i$ in $T\backslash S$: set any $z=\lfloor|T\cup S|/2\rfloor-
sum_{S}(x)$ variables $x_{i}$ to $1$ and the rest $|T\backslash S|-z$
variables to $0$. Then we update $S:=S\cup T$. After that we answer the query
with $MAJ_{T}(x)$.

Clearly, $sum_{S}(x)=\lfloor|S|/2\rfloor$ at each step. If algorithm has made
less than $\lceil n/k\rceil$ requests, then $|S|<n$, so there are some
variables that the algorithm knows nothing about.

There are $n-|S|$ such variables. If we set $x_{i}$ with $i\not\in S$ to $1$,
then $\sum_{i}x_{i}=\lfloor|S|/2\rfloor+n-|S|=n-\lceil|S|/2\rceil\geq
n-\lceil(n-1)/2\rceil\geq n-n/2=n/2$, and $MAJ_{n}(x)$ will be equal to $1$.
If we set $x_{i}$ with $i\not\in S$ to $0$, then
$\sum_{i}x_{i}=\lfloor|S|/2\rfloor\leq\lfloor(n-1)/2\rfloor\leq(n-1)/2<n/2$,
and $MAJ_{n}(x)$ will be equal to $0$.

So, whichever answer the algorithm chooses, there will be two vectors $x$ and
$y$ consistent with all our answers, such that $MAJ_{n}(x)\not=MAJ_{n}(y)$, so
on either $x$ or $y$ the algorithm will err.

###  4.2. Upper bound in adjustable-threshold setting

######  Lemma 4.3.

There exists an algorithm for adjustable-threshold adaptive model that
determines $MAJ_{n}(x)$ in $\lceil\frac{n}{k}\rceil\lceil\log(k+1)\rceil$
queries to the oracle.

######  Proof 4.4.

Split $\\{1,\ldots,n\\}$ into $\lceil\frac{n}{k}\rceil$ disjoint blocks
$B_{i}$ of size at most $k$ each. We know that
$MAJ_{B_{i}}(x;0)=1,MAJ_{B_{i}}(x;|B_{i}|+1)=0$, so using binary search we can
find such value $h_{i}$ that
$MAJ_{B_{i}}(x;h_{i})=1,MAJ_{B_{i}}(x;h_{i}+1)=0$. This means that
$sum_{B_{i}}(x)=h_{i}$. Since binary search works in
$\lceil\log(|B_{i}|+1)\rceil$ steps, in
$\lceil\frac{n}{k}\rceil\lceil\log(|B_{i}|+1)\rceil$ queries we will know
$\sum_{i}x_{i}=\sum_{i}sum_{B_{i}}(x)=\sum_{i}h_{i}$, and so we will easily
find $MAJ_{n}(x)=[\sum_{i}h_{i}\geq\frac{n}{2}]$.

The upper bound for the adaptive model with adjustable threshold turns out to
be at most two times better than the one we get in the case of the adaptive
model with fixed threshold.

###  4.3. Upper bound in fixed-threshold setting

######  Lemma 4.5.

There exists an algorithm such that it accepts two disjoint sets $A,B$,
$|A|=|B|\leq k$, $MAJ_{A}(x)\not=MAJ_{B}(x)$, and if given access to the
oracle from the definition of fixed-threshold adaptive model it finds in at
most $\lceil\log(|A|+1)\rceil$ queries to the oracle a set $S\subset A\cup B$
such that $sum_{S}(x)=|S|/2$ and $|S|\geq|A|-1$.

######  Proof 4.6.

We will define such procedure find_balanced_set(A, B, x).

Let $A_{1},\ldots,A_{|A|}$ be some enumeration of elements from $A$ and
$B_{1},\ldots,B_{|B|}$ be some enumeration of elements from $B$.

Define $c=|A|$, $i(m)=A_{m}$ if $m\leq c$, otherwise $i(m)=B_{m-c}$, and
$v(m)=x_{i(m)}$. Also define $f(m)=MAJ_{c}(v(m),\ldots,v(m+c-1))$. From the
specification of find_balanced_set(A, B, x) we know that $f(1)\not=f(c+1)$, so
we can use binary search to find position $h$ such that
$f(1)=f(h)\not=f(h+1)=f(c+1)$. This implies that $v(h)\not=v(h+c)$, and,
moreover, $v(h)=f(h)$.

If $c$ is even, this implies
$sum_{\\{i(h+1),\ldots,i(h+c-1)\\}}=\frac{c}{2}-1$, so either
$sum_{\\{i(h),\ldots,i(h+c-1)\\}}=\frac{c}{2}$, if $f(h)=1$, or
$sum_{\\{i(h+1),\ldots,i(h+c)\\}}=\frac{c}{2}$, if $f(h)=0$. The procedure
find_balanced_set(A, B, x) then returns the set
$\\{i(h+1-f(h)),\ldots,i(h+c-f(h))\\}$.

If $c$ is odd, $sum_{\\{i(h+1),\ldots,i(h+c-1)\\}}=\frac{c-1}{2}$, so the
procedure returns the set $\\{i(h+1),\ldots,i(h+c-1)\\}$.

Clearly, this procedure makes at most $\lceil\log(c+1)\rceil$ queries.

######  Claim 1.

Suppose that for some set $S\subset\\{1,\ldots,n\\}$ we know that
$sum_{S}(x)=|S|/2$. Then, $MAJ_{n}(x)=MAJ_{\\{1,\ldots,n\\}\backslash S}(x)$,
so we can forget about indices from $S$ and not take them into consideration,
as though vector $x$ has size $n-|S|$.

The basic idea of our adaptive algorithm is simple: on each iteration we find
some large set $S$ with $sum_{S}(x)=|S|/2$, and remove the variables with
indices in $S$ by claim 1. We find such set $S$ by splitting
$\\{1,\ldots,n\\}$ into $n/k$ subsets of similar sizes, querying each of them,
finding two sets of sizes $\approx|S|$ with different $MAJ$ values on them,
and running the find_balanced_set procedure from lemma 4.5. This is the main
idea, but we have to be careful with the sizes of subsets, since $k$ may not
divide $n$ evenly, and if we are too careless, after several iterations we may
have to split the set of current indices again and make queries for each of
new subsets, which could significatnly increase the total number of queries.

######  Theorem 4.7.

There exists an algorithm for computing $MAJ_{n}$ function in the fixed-
threshold adaptive setting using at most $2(\frac{n}{k-4}+1)(\log k+4)$
queries.

######  Proof 4.8.

Split the set $\\{1,\ldots,n\\}$ into $\lceil\frac{n}{k}\rceil$ disjoint sets
$S_{i}$ of sizes $\ell=\lceil\frac{n}{\lceil n/k\rceil}\rceil$ or $\ell-1$
each. Enumerate $S_{i}$ in such way that $S_{1},\ldots,S_{a}$ are of size
$\ell$ each, and the rest are of size $\ell-1$. $\ell\leq k$, since
$\frac{n}{\lceil n/k\rceil}\leq\frac{n}{n/k}=k$. For all $S_{i}$ query
$MAJ_{S_{i}}(x)$. If all the answers are the same and equal to
$y\in\\{0,1\\}$, then $MAJ_{n}(x)=y$ and we are done.

Case 1. Suppose that there are sets $S_{i}$ of both sizes $\ell$ and $\ell-1$.
Then, since not all the answers are the same, there are two indices $i$ and
$j$ such that $|S_{i}|=\ell$, $|S_{j}|=\ell-1$ and
$MAJ_{S_{i}}(x)\not=MAJ_{S_{j}}(x)$ (otherwise, all $S_{i}$ with size $\ell$
have $MAJ_{S_{i}}(x)=MAJ_{S_{j}}(x)$ for all $S_{j}$ of size $\ell-1$, so all
the subsets $S_{i},S_{j}$ have the same answer, a contradiction).

Now we will find a large (with size $\approx\ell$) subset
$S\subset\\{1,\ldots,n\\}$ of indices to exclude from our consideration (by
claim 1) in at most logarithmic in $\ell$ number of queries.

Choose any element $y$ from $S_{i}$ and query
$MAJ_{S_{i}\backslash\\{y\\}}(x)$. Suppose that
$MAJ_{S_{i}}(x)\not=MAJ_{S_{i}\backslash\\{y\\}}(x)$. Clearly, $x_{y}$ then is
equal to $MAJ_{S_{i}}(x)$.

If $MAJ_{S_{i}}(x)$ is $0$, then
$sum_{S_{i}\backslash\\{y\\}}(x)\geq\frac{\ell-1}{2}$ and
$sum_{S_{i}\backslash\\{y\\}}(x)<\frac{\ell}{2}$, implying $\ell$ is odd and
$sum_{S_{i}\backslash\\{y\\}}(x)=\frac{\ell-1}{2}$. This means that we can
ignore indices from $S_{i}\backslash\\{y\\}$ by claim 1.

If $MAJ_{S_{i}}(x)$ is $1$, then
$sum_{S_{i}\backslash\\{y\\}}(x)<\frac{\ell-1}{2}$ and
$sum_{S_{i}\backslash\\{y\\}}(x)+1\geq\frac{\ell}{2}$, implying $\ell$ is even
and $sum_{S_{i}}(x)=\frac{\ell}{2}$. Again, this means that we can remove
indices from $S_{i}$ from our consideration by claim 1.

So, if there are sets $S_{i}$ with different sizes, we are left with the case
when $MAJ_{S_{i}}(x)=MAJ_{S_{i}\backslash\\{y\\}}(x)$.

Let $M=S_{i}\backslash\\{y\\}$. We know that $|M|=|S_{j}|=\ell-1$ and that
$MAJ_{M}(x)\not=MAJ_{S_{j}}(x)$.

We run find_balanced_set from lemma 4.5 on sets $M$ and $S_{j}$, and get the
set $S$ such that $sum_{S}(x)=|S|/2$. Note that $S$ will be of size $\ell-2$
(in case when $\ell$ is even) or $\ell-1$. Since total number of elements in
$M\cup S_{j}$ is $2\ell-2$, $|M\cup S_{j}\backslash S|\in\\{\ell-1,\ell\\}$.

We then forget about indices in $S$ by claim 1, query $MAJ_{M\cup
S_{j}\backslash S}(x)$ and repeat the case 1 until all $S_{i}$ are of the same
size or have the same answer.

After several iterations we either get sets $S_{i}$ with the same answer
$MAJ_{S_{i}}(x)$ each and we are done, or we get the sets with the same sizes
$|S_{i}|$.

Case 2. All of the sets $S_{i}$ have same sizes.
Take two sets $S_{i}$ and $S_{j}$ with different answers
$MAJ_{S_{i}}\not=MAJ_{S_{j}}$ and run the procedure find_balanced_set(S_i,
S_j, x). Let $S$ be the answer. If $|S_{i}|$ is even, $|S|=|S_{i}|$, so just
exclude $S$ from further considerations, query the oracle $MAJ_{S_{i}\cup
S_{j}\backslash S}(x)$, and repeat the procedure until all the sets left have
the same $MAJ_{S_{i}}(x)$.

The case of odd $|S_{i}|$ requires more careful consideration, since if we
simply leave elements from $S$ (which has size $|S|=|S_{i}|-1$) out,
$|S_{i}\cup S_{j}\backslash S|$ will be $|S_{i}|+1$, and that may become
greater than $k$ after several iterations. If we look at the case of odd $c$
in find_balanced_set, we will see that instead of the set
$\\{i(h+1),\ldots,i(h+c-1)\\}$ of size $c-1$ it can as well return the set
$\\{i(h),\ldots,i(h+c)\\}$ of size $c+1$, since $v(h)\not=v(h+c)$.

So if $|S_{i}\cup S_{j}\backslash S|>k$, then we add elements $i(h),i(h+c)$ to
$S$: $S:=S\cup\\{i(h),i(h+c)\\}$. We can remove indices from $S$ from future
considerations by claim 1, and we query the oracle the value of
$MAJ_{S_{i}\cup S_{j}\backslash S}(x)$. Now we again have sets $S_{i}$ of two
different sizes, so we return to case 1.

The analysis of number of queries this algorithm makes is rather
straightforward: it starts with making $\lceil\frac{n}{\ell-1}\rceil$ queries,
then on each iteration it discards at least $\ell-2$ indices after making at
most $\lceil log(\ell+1)\rceil+2$ queries. So the total number of queries made
is
$\lceil\frac{n}{\ell-1}\rceil+\lceil\frac{n}{\ell-2}\rceil(\lceil\log(\ell+1)\rceil+2)$.
As was stated in the beginning, $\ell\leq k$. On the other hand,
$\ell\geq\frac{n}{\lceil
n/k\rceil}\geq\frac{n}{1+n/k}=\frac{1}{1+k/n}k\geq\frac{k}{2}$. This gives the
following upper bound on the number of queries:

$2\frac{n}{k-2}+1+2(\frac{n}{k-4}+1)(\log k+3)\leq 2(\frac{n}{k-4}+1)(\log
k+4)$

##  5\. Conclusion

We have presented new upper bounds for computing majority functions
$MAJ_{n}(x)$ in two models.

In static model we have constructed the depth-two circuit for $MAJ_{n}$ with
$MAJ_{\frac{2}{3}n+4}$ gates, which gives the first nontrivial upper bound for
this problem.

We have defined an adaptive model with two variations, and presented an
algorithm for computing $MAJ_{n}$ function using $2(\frac{n}{k-4}+1)(\log
k+3)\approx 2\frac{n}{k}\log k$ queries to the oracle in the fixed-threshold
setting, and a simple algorithm for $MAJ_{n}$ with $\lceil
n/k\rceil\lceil\log(k+1)\rceil\approx\frac{n}{k}\log k$ queries in the
adjustable-threshold setting.

There is a large gap between current lower and upper bounds for the static
model: the best lower bound [1] is $k\geq\Omega(n^{4/5})$, and our upper bound
is $\frac{2}{3}n+4$. In the adaptive model there is also a gap between the
bounds, albeit a logarithmic (in $k$) one: our lower bound is $\lceil
n/k\rceil$. So the natural questions for future consideration are those of
narrowing these gaps.

Another interesting possible direction in the adaptive case is the addition of
noise to oracle’s answers — with some probability $\varepsilon$ oracle gives
us a uniformly random answer.

##  6\. Acknowledgements

I would like to thank Vladimir Podolskii for posing the problem and for
fruitful discussions.

## References

  * [1] Christian Engels, Mohit Garg, Kazuhisa Makino, and Anup Rao.  On expressing majority as a majority of majorities.  Electronic Colloquium on Computational Complexity (ECCC), 24:174, 2017.  URL: https://eccc.weizmann.ac.il/report/2017/174.
  * [2] David Eppstein and Daniel S. Hirschberg.  From discrepancy to majority.  In LATIN 2016: Theoretical Informatics - 12th Latin American Symposium, Ensenada, Mexico, April 11-15, 2016, Proceedings, pages 390–402, 2016\.  URL: http://dx.doi.org/10.1007/978-3-662-49529-2_29, doi:10.1007/978-3-662-49529-2_29.
  * [3] Oded Goldreich.  Valiant’s polynomial-size monotone formula for majority, 2001.  Available at http://www.wisdom.weizmann.ac.il/~oded/PDF/mono-maj.pdf.
  * [4] Stasys Jukna.  Boolean Function Complexity - Advances and Frontiers, volume 27 of Algorithms and combinatorics.  Springer, 2012.  URL: http://dx.doi.org/10.1007/978-3-642-24508-4, doi:10.1007/978-3-642-24508-4.
  * [5] Stasys Jukna, Alexander A. Razborov, Petr Savický, and Ingo Wegener.  On P versus NP cap co-NP for decision trees and read-once branching programs.  Computational Complexity, 8(4):357–370, 1999.  URL: http://dx.doi.org/10.1007/s000370050005, doi:10.1007/s000370050005.
  * [6] Jesse Kamp and David Zuckerman.  Deterministic extractors for bit-fixing sources and exposure-resilient cryptography.  SIAM J. Comput., 36(5):1231–1247, 2007.  URL: http://dx.doi.org/10.1137/S0097539705446846, doi:10.1137/S0097539705446846.
  * [7] Alexander S. Kulikov and Vladimir V. Podolskii.  Computing majority by constant depth majority circuits with low fan-in gates.  In 34th Symposium on Theoretical Aspects of Computer Science, STACS 2017, March 8-11, 2017, Hannover, Germany, pages 49:1–49:14, 2017.  URL: http://dx.doi.org/10.4230/LIPIcs.STACS.2017.49, doi:10.4230/LIPIcs.STACS.2017.49.
  * [8] Frédéric Magniez, Ashwin Nayak, Miklos Santha, Jonah Sherman, Gábor Tardos, and David Xiao.  Improved bounds for the randomized decision tree complexity of recursive majority.  Random Struct. Algorithms, 48(3):612–638, 2016.  URL: http://dx.doi.org/10.1002/rsa.20598, doi:10.1002/rsa.20598.
  * [9] Marvin Minsky and Seymour Papert.  Perceptrons - an introduction to computational geometry.  MIT Press, 1987.
  * [10] Elchanan Mossel and Ryan O’Donnell.  On the noise sensitivity of monotone functions.  Random Struct. Algorithms, 23(3):333–350, 2003.  URL: http://dx.doi.org/10.1002/rsa.10097, doi:10.1002/rsa.10097.
  * [11] Ryan O’Donnell.  Analysis of boolean functions.  Cambridge University Press, 2014.