# Large-scale convective flow sustained by thermally active Lagrangian tracers

Lokahith Agasthya1 2 3<EMAIL_ADDRESS>Andreas Bartel2  Luca Biferale1
Matthias Ehrhardt2  Federico Toschi4  1Department of Physics, University of
Rome “Tor Vergata” and INFN, Via della Ricerca Scientifica 1, 00133, Rome RM,
Italy 2Angewandte Mathematik und Numerische Analysis, Bergische Universität
Wuppertal, Gaußstrasse 20 D-42119 Wuppertal, Germany 3Computation-based
Science and Technology Research Center, The Cyprus Institute, 20 Kavafi Str.,
Nicosia 2121, Cyprus 4Department of Applied Physics, Eindhoven University of
Technology, The Netherlands

###### Abstract

Non-isothermal particles suspended in a fluid lead to complex interactions –
the particles respond to changes in the fluid flow, which in turn is modified
by their temperature anomaly. Here, we perform a novel proof-of-concept
numerical study based on tracer particles that are thermally coupled to the
fluid. We imagine that particles can adjust their internal temperature
reacting to some local fluid properties and follow simple, hard-wired active
control protocols. We study the case where instabilities are induced by
switching the particle temperature from hot to cold depending on whether it is
ascending or descending in the flow. A macroscopic transition from a stable to
unstable convective flow is achieved, depending on the number of active
particles and their excess negative/positive temperature. The stable state is
characterized by a flow with low turbulent kinetic energy, strongly stable
temperature gradient, and no large-scale features. The convective state is
characterized by higher turbulent kinetic energy, self-sustaining large-scale
convection, and weakly stable temperature gradients. The particles
individually promote the formation of stable temperature gradients, while
their aggregated effect induces large-scale convection. When the Lagrangian
temperature scale is small, a weakly convective laminar system forms. The
Lagrangian approach is also compared to a uniform Eulerian bulk heating with
the same mean injection profile and no such transition is observed. Our
empirical approach shows that thermal convection can be controlled by pure
Lagrangian forcing and opens the way for other data-driven particle-based
protocols to enhance or deplete large-scale motion in thermal flows.

###### keywords:

Flow control, Bénard convection, Turbulent convection

MSC Codes  76F70.

##  1 Introduction

Thermally driven flows play an important role in both nature and industry.
They are notoriously hard to predict and control. In the presence of gravity,
temperature fluctuations cause density fluctuations, which in turn drive
convective motions through buoyancy in the atmosphere (Salesky & Anderson,
2018; Markowski, 2007), in oceans (Marshall & Schott, 1999), and especially in
idealized systems such as Rayleigh-Bénard convection (Ahlers et al., 2009;
Lohse & Xia, 2010), and horizontal convection (Gayen et al., 2014).

It is well known that the two-way interactions between particles suspended in
a fluid and the fluid phase itself are complex and highly nonlinear. They
exhibit behaviour such as preferential concentration due to ejection from
vortical regions (Squires & Eaton, 1991; Cencini et al., 2006) and
modification of turbulence (Yang & Shy, 2005). The dynamics of particles
suspended in turbulence plays an important role in several natural as well as
industrial processes, for example in the dispersal of pollutants (Fernando et
al., 2010), clouds (Falkovich et al., 2002; Mazin, 1999), planet formation
(Bec et al., 2014), combustion of jet sprays (Irannejad et al., 2015).

When suspended particles are thermally coupled to the fluid and are non-
isothermal, the particles cause local temperature fluctuations in the fluid,
which in turn can further modify a turbulent flow, either purely by thermal
action or also in conjunction with the momentum-coupling (Carbone et al.,
2019) while momentum coupling alone can also alter the heat-transfer dynamics
of a thermal flow (Elperin et al., 1996). Modification of specific thermal
flows due to suspended, thermal particles has also been studied, for example
in the Rayleigh-Bénard convection (Park et al., 2018), where heavy particles
with fixed initial temperatures are introduced into a Rayleigh-Bénard
convection system. In this case, the particles are found to enhance vertical
heat transfer, an effect that is most pronounced when the particle
concentration is greatest due to turbulence (preferential concentration),
while attenuating turbulent kinetic energy due to momentum-coupling.
Furthermore, the feasibility of achieving control of Rayleigh-Bénard
convection solely by applying small temperature or velocity fluctuations has
been studied (Tang & Bau, 1994). Here, deviations from the stable profiles
near the thermal boundaries are detected and compensated, leading to stable
Rayleigh-Bénard flows well above the critical Rayleigh number and also the
possibility of control of flow patterns is given. Increasing the critical
Rayleigh number and delaying the onset of convection can further be improved
by applying reinforcement learning techniques to apply the temperature
fluctuations near the boundary, as shown by (Beintema et al., 2020).

External radiation acting solely by heating particles suspended in a flow have
shown to modify the global motion and to induce turbulent thermal convection.
The work of Zamansky et al. (2014, 2016) considered a transparent fluid with
suspended inertial particles subject to a constant radiation and at local
thermal equilibrium with the fluid. Convection induced in such a system was
found to be driven by individual plumes rising out of each particle with
turbulent kinetic energy being the largest in the presence of a strong
particle preferential concentration where the plumes of individual particles
are reinforced by one-another due to their spatial proximity. This eventually
led to a sustained turbulent thermal convection, albeit with the temperature
of the system constantly increasing due to the permanently applied incident
radiation.

Internally heated convection (IHC) – induced and sustained by the application
of a bulk heating term in a fluid – has also been studied as an idealised
theoretical model by Wang et al. (2021). They consider a uniformly heated
domain with the top and bottom walls kept at the same constant temperature. In
this scenario, the bulk attains a stationary temperature depending on the
strength of the heating and other parameters such as gravity or the height of
the domain, while the fixed temperature boundaries works as a sink of heat,
ensuring that the temperature does not increase indefinitely.

The study of fluid systems where the heating in the bulk rather than boundary
forcing is the dominant mode of thermal forcing has important implications for
several natural systems. For example, in the mantle of the earth, the
radiogenic heating from the decay of radioactive elements plays a significant
role in addition to the heat transfer from the hotter inner core (Lay et al.,
2008). The atmosphere of Venus which contains a high amount of sulphurous
gases absorbs a large part of the incoming solar radiation, making this the
dominant mode of heat transfer (Tritton, 1975) in contrast to the earth where
the majority of the radiation is absorbed by the land surface and in-turn
forces the atmosphere. The mantle of Venus is driven in large part by internal
heating (Limare et al., 2015). Finally, in industrial applications, chiefly in
the interior of liquid-metal batteries, convection due to internal heating is
of crucial importance (Kim et al., 2013).

In this study we set up a ”theoretical experiment” to study the possibility of
controlling the global properties of a thermal flow by applying temperature
fluctuations locally along particle trajectories. In our proposed
idealisation, the particles are equipped with a hard-wired active protocol
capable of releasing or absorbing heat by setting the temperature of each
Lagrangian, thermally active tracer as a function of the local velocity field
of the underlying fluid background. Our system is internally heated/cooled by
these virtual particles so that the average heating term $\Phi$ is
statistically zero and hence the average temperature attained by the fluid is
unchanged by the forcing. The heat injection by the particles is the only
energy source for the system, since the horizontal boundaries are periodic and
the top and bottom walls are adiabatic. The aims of the set-up are multifold.
First, as a proof of concept, we wish to demonstrate that it is possible to
invent hard-wired Lagrangian protocols that can cause global flow transitions.
Second, we hope to trigger more studies using phenomenological or data-driven
approaches to achieve control of complex systems. Finally, by acting on
thermal plumes, we hope one can better understand their role in determining
the organisation of the global flow.

The remainder of the article is organised in the following manner. In Section
2, we introduce the model equations for the system, the particle temperature
protocol and describe the numerical experiments conducted. In Section 3, we
present and discuss our main findings from the numerical experiments and
finally in Section 4, we present our conclusions as well as possible future
directions for further investigation.

##  2 Methods

The protocol for particle forcing is as follows: virtual tracer particles are
initially randomly placed in a 2D region of length $L_{x}$ and height $L_{z}$
with a fluid at rest. The initial temperature of the fluid is set to an
unstable configuration with warmer temperatures at the bottom of the domain
and colder temperatures at the top of the domain. The particles are idealised
to have an infinite heat capacity and a temperature determined by an imposed
protocol in which rising particles moving vertically upward are warm with a
positive temperature $T_{+}$, while the temperature of falling particles is
set to $-T_{+}$ (see figure 1) so the average temperature of the fluid remains
constant. The temperature of the fluid near the particle relaxes to the
temperature of the particle at a rate proportional to the difference between
the local fluid temperature $T$ and the particle temperature $T_{p}$, with a
relaxation time $\tau=1/\alpha$.

Figure 1: An overview of the methods applied in the study. The domain consists
of adiabatic walls at the top and bottom while the lateral boundaries are
periodic. In inset (i) and (ii), we show a rising hot particle with
temperature $T_{+}$ and a falling cold particle with temperature $-T_{+}$
respectively.

###  2.1 Fluid Equations of Motion

The fluid velocity $\boldsymbol{u}=(u,v)$ and temperature $T$ follow the
equations

$\displaystyle\boldsymbol{\nabla}\cdot\boldsymbol{u}=0,$ (1)
$\displaystyle\partial_{t}\boldsymbol{u}+(\boldsymbol{u}\cdot\boldsymbol{\nabla})\boldsymbol{u}=-\boldsymbol{\nabla}p+\nu\nabla^{2}\boldsymbol{u}-\beta
T\boldsymbol{g},$ (2) $\displaystyle\frac{\partial T}{\partial
t}+\boldsymbol{u}\cdot\boldsymbol{\nabla}T=\kappa\nabla^{2}T-\sum_{i=1}^{N_{p}}\Bigl{(}\alpha_{i}(\boldsymbol{r},t)\bigl{[}T(\boldsymbol{r},t)-T_{i}(t)\bigr{]}\Bigr{)}.$
(3)

where (1) and (2) are the incompressible Navier-Stokes equations for a fluid
with unit density and average temperature $T_{0}=0$ with a buoyancy-force term
according to the Boussinesq approximation, where the density variations are
small and enter the equations only via the gravity-force term. Here $p$ is the
fluid pressure, $\nu$ is the kinematic viscosity, and $\beta$ is the thermal
expansion coefficient. Temperature is advected and diffused by Equation (3)
where $\kappa$ the thermal conductivity and the last term on the rhs is a heat
source term (i.e., a thermal forcing) that depends on the particles (see
later).

The domain is periodic in the horizontal $x$-direction while the top and
bottom walls at $z=0$ and $z=L_{z}$ are adiabatic with $\boldsymbol{u}=0$,
that is

$\displaystyle\partial_{z}T|_{z=0}=\partial_{z}T|_{z=L_{z}}=0,$ (4)
$\displaystyle\boldsymbol{u}(z=0)=\boldsymbol{u}(z=L_{z})=0.$ (5)

Note that the only source of energy injected into the system is the heat
supplied by the particles.

###  2.2 Equations of Particle Motion

Each particle is assumed to be a point-like, thermally active tracer. The
$N_{p}$ particles with positions
$\\{\boldsymbol{r}_{1},\boldsymbol{r}_{2},\dots,\boldsymbol{r}_{N_{p}}\\}$ and
temperatures $\\{T_{1},T_{2},\dots,T_{N_{p}}\\}$ follow the local fluid
velocity

$\frac{d\boldsymbol{r}_{i}}{dt}=\boldsymbol{u}(\boldsymbol{r}_{i}(t),t).$ (6)

To mimic an effective particle size, concerning its thermal properties, we
imagine that each particle exerts a thermal forcing on the fluid in its
immediate vicinity up to a cut-off distance $\eta$. The feedback of the
particle is defined as a local heat injection term proportional to the
difference between the underlying fluid temperature, at the location of the
particle, and the instantaneous particle temperature. Furthermore, to have a
smooth thermal forcing, we assume that the strength of the coupling
$\alpha_{i}$ (with dimension inverse of time) between the $i$-th particle and
the fluid at time $t$ and position $\boldsymbol{r}$ has the form of a Gaussian
with a peak at the particle location $\boldsymbol{r}_{i}(t)$ (see inset (ii)
of Fig. 1), given by

$\alpha_{i}(\boldsymbol{r},t)=\begin{cases}\alpha_{0}\,\exp{\Bigl{(}-\frac{|\boldsymbol{r}-\boldsymbol{r}_{i}(t)|^{2}}{2c^{2}}\Bigr{)}},&\text{
if }|\boldsymbol{r}-\boldsymbol{r}_{i}(t)|\leq\eta,\\\ 0,&\text{ if
}|\boldsymbol{r}-\boldsymbol{r}_{i}(t)|>\eta.\end{cases}$ (7)

Here, $\alpha_{0}$ is the coupling strength at the particle location and $c$
is the size of the virtual particle (referred to as particle size). In fact,
$c$ determines the sharpness of the peak of the Gaussian function
$\alpha_{i}$: the Gaussian peaks more sharply and falls off more quickly for
smaller $c$. On the other hand, $\eta$ is simply a cut-off length for the
thermal forcing by the particle. While setting a cut-off length introduces a
discontinuity in the system, we note that for the typical parameters chosen in
the study, the value of $\alpha$ at the cut-off length is $\approx
0.011\alpha_{0}$. This value is only $1$% of the value at
$\boldsymbol{r}=\boldsymbol{r}_{i}$. The cut-off is used to save computational
resources and to avoid calculating vanishingly small values of the Gaussian
curve for points far away from a particle (large values of
$|\boldsymbol{r}-\boldsymbol{r}_{i}|)$.

The thermal forcing due to the $i$-th particle $\Phi_{i}$ at location
$\boldsymbol{r}$ is

$\Phi_{i}(\boldsymbol{r},t)=-\alpha_{i}(\boldsymbol{r},t)\bigl{[}T(\boldsymbol{r},t)-T_{i}(t)\bigr{]},$
(8)

and the total thermal forcing at a given location $\boldsymbol{r}$ due to all
$N_{p}$ particles reads

$\Phi(\boldsymbol{r},t)=-\sum_{i=1}^{N_{p}}\Bigl{(}\alpha_{i}(\boldsymbol{r},t)\bigl{[}T(\boldsymbol{r},t)-T_{i}(t)\bigr{]}\Bigr{)}.$
(9)

Here $\Phi(\boldsymbol{r},t)$ is of exactly the same form as eq. (3). It
defines the relationship between the particle positions and the thermal
forcing at a given location in the fluid, and thus can be directly substituted
into the heat-equation, which is an Eulerian description of the temperature.

To summarise, each particle influences a fixed region surrounding itself and
when two particles are within distance $2\eta$, their thermal effects are
additive in the overlapping region. The coupling between the particle and the
fluid is two-way, since the particle actively heats or cools the fluid, while
the fluid acceleration due to this thermal forcing also accelerates the
particle.

###  2.3 Particle Temperature Policy

The temperatures of the particles are determined by a binary policy where the
$i$-th particle has either a positive value $T_{+}$ or a negative value
$-T_{+}$ depending on the sign of the vertical velocity of the particle
$v_{i}(t)$:

$T_{i}=\begin{cases}T_{+},&\text{if }v_{i}>0,\\\ -T_{+},&\text{if
}v_{i}<0.\end{cases}$ (10)

Since the particle is a tracer, $v_{i}$ is the same as the vertical velocity
of the fluid at the particle location $v(\boldsymbol{r}_{i},t)$. $T_{+}$ is a
parameter that sets the temperature scale of the system. By heating the upward
moving fluid regions and conversely, cooling the downward moving regions, this
policy should enhance thermal convection by amplifying any updrafts or
downdrafts if they exist. Particles are coupled to each other via their
effects on the fluid and because of the flow thermal diffusivity.

Our policy leads to a sharp discontinuity in the particle temperature when the
particle changes direction. Furthermore, the temperature would rapidly
fluctuate between $T_{+}$ and $-T_{+}$ at the top and bottom walls where the
velocity is very small and the flow is mainly horizontal. To ensure that this
doesn’t affect our results, we verified that setting $T_{i}=0$ for particles
within one grid length from the top and bottom walls, where the vertical
velocity of the particle fluctuates rapidly from small positive values to
small negative values, leads to (statistically) the same flows. We have also
verified that all results reported below are robust against small change of
the above protocol, e.g. by setting a threshold velocity $v_{0}$ such that
$T_{i}=0$ when $|v|<v_{0}$.

Parameter |  $L_{x}$ |  $L_{y}$ |  $\nu$ |  $\kappa$ |  $\alpha_{0}$ |  $T_{+}$ |  $g$ |  $c$ |  $\beta$ |  $N_{p}$
---|---|---|---|---|---|---|---|---|---|---
Range |  864 |  432 |  $\frac{1}{1500}$ |  $\frac{1}{1500}$ |  $[0.0001-0.005]$ |  $[2.5\times 10^{-8}-0.05]$ |  $8\times 10^{-6}$ |  $[0.5-\sqrt{2}]$ |  1 |  $[48-960]$
Table 1: List of parameters used in the study along with the range of values
in simulation units.

###  2.4 Numerical Experiments

The fluid equations (1)–(3) are solved by the Lattice-Boltzmann method (see
Appendix A for details), together with the particle evolution as a tracer
given by equation (6). The particle evolution is solved by the two-step Adams-
Bashforth method. We start from an initially unstable vertical temperature
profile of

$T(z)=T_{+}\tanh\Bigl{(}\frac{L_{z}}{2}-z\Bigr{)}.$ (11)

The two-way coupled particle-fluid system is evolved until the flow reaches a
statistically stationary kinetic energy independent of the initial conditions
for the flow velocity, temperature and particles positions. All measurements
and analyses are performed at this steady state for different sets of
parameters, varying $T_{+}$, $N_{p}$, $\alpha_{0}$, and $c$. The cut-off
distance for the particles $\eta$ is kept constant throughout the study.

All results presented in this study are for a 2D fluid domain resolved with
864 grid points in the horizontal direction and 432 grid points in the
vertical direction. With the Lattice Boltzmann grid spacing $\Delta x=1$, we
have $L_{x}=864$ and $L_{z}=432$. The particles have a fixed cut-off distance
$\eta=3$ in computational units, while their size $c$, is varied. $\alpha_{0}$
is varied from $10^{-4}$ to $5\times 10^{-3}$ in simulation units. The
temperature $T_{+}$ is varied over several orders of magnitude. All
temperatures in this study are reported in units of $T_{s}/0.025$ where
$T_{s}$ is the temperature in simulation units. Thus, $T=0.1$ corresponds to a
temperature of $T_{s}=0.0025$ in simulation units. This convention is chosen
solely to make it easier to compare the scales of the various $T_{+}$ and make
the manuscript more readable. The values of the parameters are summarised in
table 1.

In order to have dimensionless quantities, we define a typical velocity
$u_{0}$, given by

$u_{0}=\sqrt{cg\beta\frac{\alpha_{0}}{\alpha_{0}+\frac{\kappa}{2c^{2}}}T_{+}},$
(12)

where $c$ is the size of the particle as defined in equation (7). The form
(12) was suggested by studying the evolution of single particles experiments
at varying $\alpha_{0}$ and $c$, where the rms value of the vertical particle
velocity was found to scale as in (12). In particular, we find that the
particle velocity statistics remain independent of the domain height $L_{z}$,
justifying the choice of $c$ as the length scale of the system. The fluid near
the particle relaxes to the temperature of the particle, and this relatively
hotter/cooler local plume rises/falls. The tracer particle in turn responds to
the fluid and accelerates at a rate that depends on the temperature anomaly,
gravity $g$ and $\beta$. This is similar to other thermal flows such as
Rayleigh-Bénard convection. The local heating is high when $c$ is large
because a wider region around each particle is thermally forced.

In the given protocol, the temperature scale alone is set by the particles and
this is in turn sets a natural velocity field for the system. Given that the
particles are active and the fluid-particle coupling is highly non-linear, it
would be a hard task to estimate other quantities associated with thermal
systems usually defined in the literature such as a Richardson number or a
Brunt-Valsala frequency as this would require a-priori knowledge of the
average state of the system.

The quantity

$T_{a}=\frac{\alpha_{0}}{\alpha_{0}+\frac{\kappa}{2c^{2}}}T_{+}$ (13)

is interpreted as an effective temperature reached in the vicinity of each
particle. The empirical prefactor
$\alpha_{0}/(\alpha_{0}+\frac{\kappa}{2c^{2}})$ by which $T_{+}$ is multiplied
is a constant that gives the rate of relaxation of the fluid temperature to
the particle temperature compared with the rate at which heat is diffused away
from the particle by conduction, which is proportional to $\kappa/c^{2}$. When
$\alpha_{0}\to 0$, then $T_{a}\to 0$, because the fluid is no longer coupled
to the particle and there is no energy input to the system. When
$\alpha_{0}\gg\kappa/c^{2}$, then $T_{a}\to T_{+}$, meaning the fluid attains
the local particle temperature. For large $\kappa$, the heat is rapidly
conducted away from the particle so that effective temperature is lower, where
again $T_{a}\to 0$ when $\kappa\to\infty$ while the case of small $\kappa$ is
similar to that of large $\alpha_{0}$. In our study, $\alpha_{0}$ and
$\kappa/c^{2}$ are of comparable magnitude.

Furthermore, we define the normalized turbulent kinetic energy $E_{k}(t)$ of
the system as

$E_{k}(t)=\frac{1}{2}\frac{\bigl{\langle}|\boldsymbol{u}(t)|^{2}\bigr{\rangle}_{V}}{u_{0}^{2}N_{p}},$
(14)

where $\langle\cdot\rangle_{V}$ represents the average over the entire domain
at a given time. We also define with an overline $\overline{E}_{k}$ as the
average normalized turbulent kinetic energy (TKE), i.e.

$\overline{E}_{k}=\bigl{\langle}E_{k}(t)\bigr{\rangle}_{t},$ (15)

where $\langle\cdot\rangle_{t}$ denotes the time average after the flow
reaches a statistically stationary regime. If the particles are sparse and
their motion is independent of each other, the kinetic energy of the system
would simply be a sum of the motion of the individual particles and we would
expect $\overline{E}_{k}$ of the flow to attain a constant value. However, if
the motions of the particle are not independent of each other, the variation
of $\overline{E}_{k}$ as a function of $N_{p}$ is not possible to predict
a-priori.

##  3 Results

###  3.1 Stable and Convective Configurations

First, we vary the number of virtual particles $N_{p}$. Figure 2 shows four
cases, where we visualize snapshots of the temperature and velocity fields.
Thereby the rising particle temperature $T_{+}$, particle-fluid coupling
strength $\alpha_{0}$ and particle size $c$ are fixed. The figure indicates
that there are two distinct stationary typical configurations. The first,
which we term stable, is shown in the top panels (a) and (b) of figure 2. In
this state, kinetic energy is low and large scale circulation is absent.
Particles are either nearly still and close to the top and bottom walls or
they execute a slow vertical motion independently one from the others,
propelled by their higher or lower temperature compared to the bulk. When the
particle concentration reaches beyond a certain threshold, the individual
thermal effect of the particles aggregates and triggers a transition to a
second state shown in the bottom panels (c) and (d) of figure 2. This
convective state enjoys a large scale circulation, the presence of rising and
falling plumes with the particles trajectories synchronized with the large-
scale recirculation regions. In figure 2, this transition occurs for
$N_{p}\sim 150$.

Figure 2: Snapshots of the temperature field $T(\boldsymbol{r},t)/T_{+}$ at a
given instant of time for $T_{+}=0.1$, $\alpha_{0}=0.005$, $c=1$ and at
changing $N_{p}=120,140,160,180$ in panels (a), (b), (c) and (d),
respectively. The colour palette varies from red to blue where red indicates
$T=T_{+}$ and blue indicates $T=-T_{+}$. The black arrows show the velocity
field with the length of the arrow representing the relative magnitude of the
velocity with identical scaling for all four panels. The top panels show a
stable configuration while the bottom panels show a convective configuration.
Figure 3: Time evolution of $E_{k}(t)$ for flows with (a) $T_{+}=0.02$, (b)
$T_{+}0.1$ and (c) $T_{+}=1.0$ with $\alpha_{0}=0.005$ and $c=1$ kept fixed.
Stable configurations are plotted in blue while convective configurations are
plotted in red. The time is in simulation time units.

In figure 3 we show the time evolution of the TKE for parameters before and
after the transitions. Panels (a), (b) and (c) corresponds to $T_{+}=0.02$,
$T_{+}=0.1$ and $T_{+}=1.0$, respectively, with $\alpha_{0}$ and $c$ remaining
fixed. The blue curves represent stable configurations while the red curves
represent convective configurations. The kinetic energy first increases due to
the unstable temperature gradient imposed on the initial condition. At later
times, the thermal forcing by the tracers is dominant and the flow attains a
statistically stationary kinetic energy where $E_{k}(t)$ either shows a large
value (red curves), corresponding to a convective flow shown qualitatively in
figure 2 or a low value (blue curves) corresponding to a quasi stable flow.

Two further points are note-worthy about the transition from figure 3.
Firstly, the transition is abrupt:it is enough to add very few particles to
hvae a jump $\gtrsim 5$ in the normalised kinetic energy. It should be noted
that the expression of $E_{k}(t)$ is normalized by $N_{p}$ in the denominator,
so the absolute increase in kinetic energy is even greater. Secondly, the
critical $N_{p}$ depends slightly on $T_{+}$, where for larger $T_{+}$, the
transition occurs at a slightly larger $N_{p}$. We see that in panel (a) with
$T_{+}=0.02$, the transition lies between $N_{p}=120$ and $N_{p}=140$ while in
panel (c) with $T_{+}=1.0$, the transition lies between $N_{p}=160$ and
$N_{p}=180$, with the case of $T_{+}=0.1$ in panel (b) showing an intermediate
behavior. This weak dependence on $T_{+}$ which will be further commented upon
later. It has been verified that the transitions are robust by replacing the
initial unstable profile with an initial temperature field of $T=0$ everywhere
with particles being either hot or cold with probability $0.5$ each.

Figure 4: Time-averaged vertical temperature profile divided by $T_{+}$
plotted against the vertical height for various $N_{p}$ close to the
transition $N_{p}$ for $T_{+}=0.02$ (a), $T_{+}=0.1$ (b) and $T_{+}=1.0$ (c).
Stable configurations are plotted in blue while convective configurations are
shown in red.

In figure 4, we show a comparison of the normalised time-averaged vertical
temperature profiles $\overline{T}(z)$ for the same set of flows given by

$\overline{T}(z)=\frac{\bigl{\langle}T(\boldsymbol{r},t)\bigr{\rangle}_{x,t}}{T_{+}},$
(16)

where $\langle\cdot\rangle_{x,t}$ represents the time-average at a given
height $z$. Notice that the temperature gradients for the stable flows (blue)
show a strongly stable profile ($\partial_{z}T>0$) while the convective flows
still show an overall stable temperature profile but with weaker gradients so
that the temperature difference between the top and the bottom adiabatic walls
are much smaller. In the presence of a large-scale circulation, the
temperature field is more effectively transported and mixed throughout the
domain. We also see that with increase in $T_{+}$, the convective
configurations show a flatter temperature profile for the corresponding
$N_{p}$ of lower $T_{+}$ flows, i.e., for example, the red curves in panel (c)
are much flatter than those in panel (a).

The dual-nature of the effect of the virtual particles is observed here – the
particles tend to make the flow more stable by carrying heat away from the
lower half of the domain while carrying heat towards the upper half of the
domain. Thus, the larger $T_{+}$ is, the more stable the system becomes.
However, when a certain threshold of particles is reached, the situation
changes – the virtual particles together create a persistent large-scale flow
and now the convection is strong enough to overcome the stable temperature
gradient.

Figure 5: Time averaged normalized TKE $\overline{E}_{k}$ as a function of
$N_{p}$ for various $T_{+}$ (shown in legend) for $\alpha_{0}=0.005$ and
$c=1$. The inset (in the lower right corner) shows the behavior of a flow with
$T_{+}=0.1$ very close to the transition $N_{p}$. Also shown are instantaneous
snapshots of the temperature field for a stable configuration (bottom) and a
convective configuration (top right). The error bars indicate the standard
deviation of the temporal fluctuations of $E_{k}(t)$ around the average
kinetic energy in the stationary regime.

In figure 5, we take a closer look at the transition by plotting the average
normalised TKE of the flows as defined in equation (15) against $N_{p}$, for
the same $\alpha_{0}$ as above, for various $T_{+}$. Here we see that
$\overline{E_{k}}$ initially remains constant for small values of $N_{p}$ and
then increases rapidly during the transition from a stable flow to a
convective flow. Once the convective flow has been set up, the value of
$\overline{E_{k}}$ saturates to a higher constant value for large $N_{p}$. The
sharp increase of TKE at a transition $N_{p}$ is once again clearly visible.
We empirically define a value of $E^{0}_{k}=0.225$ indicated by the horizontal
red line as the transition point where for stable end states,
$\overline{E_{k}}<E^{0}_{k}$ and vice-versa for the convective end state. The
sharpness of the transition is examined more closely in the inset of the
figure for a given $T_{+}$. It is seen that the transition occurs for an
increase of just one single particle. The dependence of the transition on
$T_{+}$ is weak, for $T_{+}$ varying over 2 orders of magnitude the transition
occurs at nearly the same $N_{p}$.

###  3.2 Large-scale Circulation and Heat Transfer

While the existence of the large-scale circulation is apparent from the
visualisations of the temperature and velocity fields, it is possible to infer
its presence quantitatively from the fluid energy spectrum. In particular, we
consider the spectrum in the horizontal direction taken at the mid-plane
$z_{0}=L_{z}/2$, given by

$E_{\boldsymbol{u}}(k_{x})=\frac{1}{2}\Bigl{\langle}\big{|}\hat{\boldsymbol{u}}(k_{x},z_{0},t)\bigr{|}^{2}\Bigr{\rangle}_{t},$
(17)

and $\hat{\boldsymbol{u}}(k_{x},z_{0},t)$ are the Fourier coefficients of the
field $\boldsymbol{u}$ and $\langle\cdot\rangle_{t}$ denotes the time
averaging. We denote by $E_{1}$ the energy contained in the first Fourier mode
with wavenumber $k_{x}=2\pi/L_{x}$, $E_{2}$ is used for energy of the second
mode ($k_{x}=4\pi/L_{x}$), and so on. Moreover, we define $E_{\rm tot}$ as the
sum of the energy contained in all the Fourier modes,

$E_{\rm tot}=\sum_{i=1}^{N_{k}}E_{i},$ (18)

where $N_{k}$ is the Fourier mode corresponding to the smallest resolved
length-scale. The strength of the large-scale circulation with a rising plume
and a falling plume can be measured by the value $E_{1}/E_{\rm tot}$ (Xi et
al., 2016), which measures the fraction of energy contained in the first mode,
that is the smallest wave number. This corresponds to a cosine mode for the
velocity field in the bulk, which is a close approximation when there exist
two counter-rotating vortices. When such a large-scale flow is present, we
would have $E_{1}/E_{\rm tot}\gg 0$, while if the flow lacks large-scale
convection, we would have a flatter energy spectrum with $E_{\rm tot}\gg
E_{1}$ and $E_{1}\sim E_{2}$.

(a)

(b)

(c)

Figure 6: (a) $E_{1}/E_{\rm tot}$ for varying $N_{p}$ for various values of
$T_{+}$. Error bars show the temporal fluctuations of $E_{1}/E_{\rm tot}$ (b)
$\operatorname{Nu}$ for varying $N_{p}$ for various values of $T_{+}$. The
black solid line shows a linear scaling with $N_{p}$. (c) Plot of the average
normalised Nusselt number $\operatorname{Nu}_{\Phi}$ against the average
normalised thermal energy injection $\overline{\epsilon}_{T}$ for flows with
varying parameters. Stable flows are marked with blue filled circles,
convective with red filled hexagons and the two black lines scale as
$(\overline{\epsilon}_{T})^{1.2}$. Inset of panel (c) shows the scaling of
$\operatorname{Nu}_{\Phi}$ as a function of $N_{p}$ for 3 different $T_{+}$.

In figure 6(a), we plot the strength of the large-scale circulation
$E_{1}/E_{\rm tot}$ for varying $N_{p}$. We see clearly here that
corresponding to a jump in the magnitude of the TKE seen in figure 5, there is
also a similar large increase in the ratio of kinetic energy contained in the
largest-scale. Given that $\overline{E}_{k}$ takes into account the typical
velocity of a single particle as well as the number of particles, the excess
kinetic energy clearly comes from the large-scale circulation that arises
after the transition, a cumulative particle effect.

Figure 6(b) shows the dimensionless Nusselt number, $\operatorname{Nu}$,
defined as

$\operatorname{Nu}=\frac{\Bigl{\langle}vT-\kappa\frac{\partial T}{\partial
z}\Bigr{\rangle}_{V,t}}{\frac{\kappa\overline{\Delta T}}{L_{z}}},$ (19)

where $\langle\cdot\rangle_{V,t}$ represents average over the entire domain
and time, $v$ is the vertical fluid velocity and $\overline{\Delta T}$ is the
time-averaged temperature difference between the top and bottom walls given by

$\overline{\Delta T}=\left\langle T(x,L_{z})\right\rangle_{x,t}-\left\langle
T(x,0)\right\rangle_{x,t}.$ (20)

Here, the Nusselt number is defined in analogy with Rayleigh-Bénard
convection: it is the ratio of heat transfer due to convection and the heat
transfer by conduction with the difference that the temperature jump is taken
in the opposite sense because of the presence of a stable mean profile. Due to
the adiabatic boundary conditions imposed at the top and bottom walls
($\partial_{z}T=0$) along with the no-slip boundary condition for the velocity
($\boldsymbol{u}=0$), the value of the Nusselt number is $0$ at the top and
bottom walls. Thus, the boundary walls do not contribute to the heat transfer.
The Nusselt number naturally increases proportionally with the number of
particles.

We see in figure 6(b) that the value of $\operatorname{Nu}$ increases
gradually with increase in $N_{p}$, followed by a large increase around the
transition $N_{p}$ and then settling to a roughly linear increase with $N_{p}$
in the convective regime. The reason for the large increase of
$\operatorname{Nu}$ at the transition is two-fold. Firstly, the increase in
TKE overall leads to an increase in the convective heat transfer which further
increases $vT$. Secondly, with more effective mixing of the temperature and a
weakly stable temperature gradient, $\overline{\Delta T}$ in the denominator
also has a smaller magnitude.

Another way to quantify the heat transfer is to divide it by the typical
forcing $\Phi$ multiplied by the length-scale of the system. This is similar
to the normalisation procedure of (Wang et al., 2021) applied to internally
heated convection (with volume forcing). The effective temperature $T_{a}$
defined in equation (13) was introduced as a typical value of the temperature
attained by the fluid in the vicinity of the particle with an associated
length-scale $c$ for each particle. In a similar vein,
$\alpha_{0}(T_{+}-T_{a})$ can be considered the typical thermal forcing acting
on the fluid. We use two dimensionless response parameters of the system.
First, we define the normalised Nusselt number $\operatorname{Nu}_{\Phi}$
given by

$\operatorname{Nu}_{\Phi}=\frac{\Bigl{\langle}vT-\kappa\frac{\partial
T}{\partial z}\Bigr{\rangle}_{V,t}}{c\alpha_{0}(T_{+}-T_{a})}.$ (21)

$\operatorname{Nu}_{\Phi}$ measures the heat transfer by convection relative
to the input typical thermal forcing multiplied by the length scale of the
system.

In the stationary regime, the thermal dissipation rate $\epsilon_{T}$ is given
by $\langle\Phi T\rangle_{V,t}$ (see Appendix C) and is normalised as

$\overline{\epsilon}_{T}=\frac{\langle\Phi
T\rangle_{V,t}}{\alpha_{0}(T_{+}-T_{a})T_{+}}.$ (22)

The normalisation factor is once again the typical forcing multiplied by the
temperature scale.

In figure 6(c) we plot the normalised Nusselt number
$\operatorname{Nu}_{\Phi}$ against the normalised thermal dissipation
$\overline{\epsilon}_{T}$, quantifying the measured convective response of the
fluid to the measured input thermal forcing for varying $T_{+}$, $c$,
$\alpha_{0}$ and $N_{p}$. It is seen that there exists a global scaling of
these two quantities for both the flow regimes, stable and convective with a
rough scaling of
$\operatorname{Nu}_{\Phi}\propto(\overline{\epsilon}_{T})^{1.2}$. However, the
higher magnitude of the normalised Nusselt number in the convective case
differentiates it from the stable flows. In the inset of panel (c) we can see
that NuΦ also shows a transition and a sharp increase in magnitude from stable
to convective.

The above findings are consistent with a situation that can be briefly
described as such – individual particles thermally coupled to the fluid have a
small zone of influence and release or absorb heat in their immediate
vicinity. Thus, each particle contributes to the thermal injection into the
domain as well as the vertical heat transfer across the domain both of which
increase with the increase in number of particles. In the stable regime, the
main effect of the particles is to maintain the strongly stable temperature
gradient. However, since there exists a maximum and minimum temperature that
the fluid can attain ($T_{+}$ and $-T_{+}$ respectively), there is
consequently a limit on the strength of the positive temperature gradient that
can be set-up in the system, independent of the number of particles. Hence,
increasing the number of particles eventually leads to a saturation of the
stable temperature gradient while it continues to further destabilise the
system. Thus, eventually at a critical $N_{p}$, the instability dominates the
flow and leads to a transition to a fully convective flow.

At the transition to the convective regime, the development of the large-scale
convective flow patterns and more turbulent flow leads to a large increase in
the heat transfer relative even to the thermal energy injection, while also
seeing a weaker stable vertical temperature gradient across the domain.

Figure 7: Average particle kinetic energy $\overline{E}_{k}^{(L)}$ plotted
against $N_{p}$ for 3 values of $T_{+}$. All measurements are made with $c=1$
and $\alpha_{0}=0.005$.

###  3.3 Particle Kinetic Energy

Given that the particles are not passive and are actively interacting with the
fluid, their behaviour is of particular interest. We define the Lagrangian
kinetic $\overline{E}_{k}^{(L)}$ as

$\overline{E}_{k}^{(L)}=\frac{1}{u_{0}^{2}N_{p}}\left\langle\sum_{i=1}^{N_{p}}\frac{1}{2}\left(u_{i}^{2}+v_{i}^{2}\right)\right\rangle_{t},$
(23)

where $u_{i}$ and $v_{i}$ are the horizontal and vertical velocity of the
$i$-th particle respectively. $\overline{E}_{k}^{(L)}$ measures the average
kinetic energy of a particle in the system. In figure 7, we plot this value
against $N_{p}$ for various $T_{+}$, which shows the characteristic jump at
the transition to a convective flow. The absolute value of
$\overline{E}_{k}^{(L)}$ is far larger than the corresponding fluid kinetic
energy $\overline{E}_{k}$ due to the fact that the particle locations are the
active regions of the flow where the energy is input into the system. Thus, in
spite of being tracers, these particles selectively sample higher kinetic
energy regions of the flow. $\overline{E}_{k}^{(L)}$ also shows an increase as
a function of $N_{p}$. The reason for this behaviour is two-fold – first, with
increasing $N_{p}$, the fluid becomes more energetic and turbulent as
discussed before. Secondly, the particles themselves are better mixed and are
found more often in the more energetic regions of the fluid away from the
stationary walls.

###  3.4 Comparison with Eulerian imposed thermal forcing

We consider a thermal fluid system with a thermal forcing $\phi(z)$ uniformly
applied at all times. The forcing is a close approximation of the vertical
profile of the measured forcing $\Phi$ in the Lagrangian system as shown in
figure 8. Defining $Q(z)$, the numerator of the Nusselt number, as the average
net heat transfer in the positive $z$ direction at height $z$ given by

$Q(z)=\Bigl{\langle}v(\boldsymbol{r},t)T(\boldsymbol{r},t)-\kappa\partial_{z}T|_{(\boldsymbol{r},t)}\Bigr{\rangle}_{x,t},$
(24)

where $\langle\cdot\rangle_{x,t}$ indicates the time and spatial averages at a
given height $z$, notice that averaging equation (3) over time gives

$\Phi(z)\coloneqq\bigl{\langle}\Phi(\boldsymbol{r},t)\bigr{\rangle}_{x,t}=\partial_{z}Q(z).$
(25)

Figure 8: The measured value of the average vertical profile of thermal
forcing $\Phi(z)$ for (a) a stable Lagrangian flow and for (b) a convective
Lagrangian flow shown along with the profile of the corresponding applied
Eulerian thermal forcing.

(a)

(b)

Figure 9: (a) The normalised temperature profile for a stable Lagrangian flow
(blue) compared with the measured temperature profile of a flow with an
identical imposed profile of thermal forcing. (b) The normalised temperature
profile for a convective Lagrangian flow (red) compared with the measured
temperature profile of a flow with an identical imposed profile of thermal
forcing Figure 10: Snapshots of the temperature fields: (a) a stable
Lagrangian flow (upper left), (c) a convective Lagrangian flow (lower left),
and the two uniformly forced flows to mimic the stable (b) and convective
flows (d) in right column. The temperature fields $T$ are divided by the
respective $T_{+}$. The black arrows show the velocity field. The length of
the arrows indicate the magnitude of the velocity field within each panel –
the arrow lengths are scaled differently for different flows to allow for the
clearest viewing of the flow structure.

The comparison is made for one stable and one convective flow. Given identical
vertical profiles of thermal forcing (see figure 8), one would expect that the
resulting temperature profile and hence the nature of the flows would remain
identical. However, as shown in figure 9, the temperature profiles show a
dramatic difference, with the Eulerian flows showing an unstable temperature
profile similar to the Rayleigh-Bénard Convection. Further, as shown in figure
10, even when the thermal forcing matches the measured value from a stable
configuration, the Eulerian flow with uniform thermal forcing shows a
convective behavior with clear, well-defined hot and cold plumes and an
unstable temperature gradient. Even in the convective Lagrangian case, the
corresponding Eulerian flow is convective. Thus, the presence of the stable
temperature gradients and the two typical configurations outlined previously
is not a result of the net thermal forcing applied on the system but of the
particular Lagrangian nature of the thermal tracers and the two-way coupling
with the fluid.

While it is possible to obtain a stable flow by imposing a comparatively very
small value of thermal forcing in the Eulerian case, it should be emphasised
here that the Lagrangian forcing is a measured quantity and not a quantity
that can similarly be set independently. The thermal forcing in the Lagrangian
case is not known a-priori due to the non-linear, 2-way coupling between the
particles and the fluid.

###  3.5 Anomalous Behavior for Small $T_{+}$

We have already noted in previous sections that there is weak dependence of
the transition of the system on the value of $T_{+}$. In particular, it was
observed that for larger $T_{+}$, the transition occurs at a larger $N_{p}$
and the stable configurations for larger $T_{+}$ have relatively flatter
temperature gradients. One would conclude then that for any given $N_{p}$,
there exists a $T_{+}$ small enough such that the system is convective.
However, at very small $T_{+}$, the system attains a third columnar state
where the temperature profile is still stable ($\partial_{z}T>0$) and the
system has a weak convective flow (see snapshot in figure 11 (b)).

(a)

(b)

Figure 11: (a) The ratio of kinetic energy contained in the first Fourier mode
$E_{1}/E_{\rm tot}$ (dashed lines) and $E_{2}/E_{\rm tot}$ (solid lines) to
the total energy contained in all modes for $N_{p}=140$ and $N_{p}=240$
plotted against $T_{+}$. Inset shows the averaged normalised TKE
$\overline{E}_{k}$ for the same parameters and the horizontal line represents
$\overline{E}_{k}=E_{k}^{0}$. (b) A snapshot of the temperature field for a
columnar flow with $N_{p}=240$ and $T_{+}=10^{-5}$. The colour palette varies
from red to blue where red indicates $T=T_{+}$ and blue indicates $T=-T_{+}$.

In figure 11 (a), we plot the fraction of energy contained in the first
Fourier mode ($E_{1}/E_{\rm tot}$) as well as the second Fourier mode
($E_{2}/E_{\rm tot}$) to understand the large-scale behavior of the flow. We
can see clearly that for smaller $T_{+}$, the second mode dominates the
kinetic energy while the energy contained in the first mode approaches $0$.
This is the case until a transition $T_{+}$, where now the flow turns
convective from columnar, with a dominance of $E_{1}$. At larger $T_{+}$ for
$N_{p}=240$ (orange, filled symbols), we see that while $E_{2}/E_{\rm tot}$
remains small, the value of $E_{1}/E_{\rm tot}$ shows a decreasing trend. This
is because as $T_{+}$ is increased, the flow becomes more turbulent and small-
scale velocity features begin to appear, increasing the energy contained at
higher modes. For $N_{p}=140$ (cyan, empty symbols), the flow is columnar for
$T_{+}\lesssim 10^{-3}$ and transitions to convective at $T_{+}\sim 0.02$, as
evidenced by the values of $E_{1}/E_{\rm tot}$ and $E_{2}/E_{\rm tot}$.
However, on increasing $T_{+}$ further, the flow again moves to a stable
configuration, as evidenced by the fact that $E_{1}\sim E_{2}$ which indicates
the lack of any large-scale velocity flow. This transition is due to the
effect already observed, that for increasing $T_{+}$, the $N_{p}$ of
transition from stable to convective is greater.

The inset of figure 11(a) shows the normalised TKE plotted for the two given
$N_{p}$ and varying $T_{+}$. Notice that at small $T_{+}$, when the flow is
columnar, it is characterised by a smaller normalised TKE and kinetic energy
smoothly approaches $0$ as $T_{+}\to 0$.

##  4 Conclusions and Discussion

We have performed numerical simulations of an idealized non-isothermal 2D
fluid system under the Boussinesq approximation with suspended tracer
particles. The particles act as heat sources or sinks depending on their
vertical velocity. The particles are coupled to the fluid only thermally, the
fluid is forced only by the action of the particles. Individually, each
particle aids in the transport of heat away from the bottom of the domain
towards the top of the domain, thus working to create a thermally more stable
system. However, under certain conditions, the cumulative effect of the
particles overpowers the tendency towards stability and the result is a system
with a large-scale convective flow pattern with increased turbulent kinetic
energy, larger heat transfer across the domain, maximum energy in the largest
Fourier modes and a (weakly) stable vertical temperature gradient. The main
parameters of the system are the temperature of the hot, rising particles
$T_{+}$, the number of particles $N_{p}$, the strength of the thermal coupling
between the fluid phase and the particles $\alpha_{0}$ and the size of the
particle $c$. Increasing $N_{p}$, $c$ and $\alpha_{0}$ makes the flow
increasingly convective while increasing $T_{+}$ weakly contributes to making
the flow more stable.

This Lagrangian protocol is compared with a system with a uniform thermal
forcing identical to the measured Lagrangian forcing and it is found that the
temperature profiles of the Eulerian system is unstable rather than stable and
a convective flow always develops.

Extension to 3D set-ups and to cases with larger domain and/or a larger number
of particles to study whether the intensity of turbulence can be increased
indefinitely would also be interesting.

Independently of the possibility to realize a protocol like the one we studied
here in a realistic experimental set-up, out study is meant to gain a new
insight about the impact of Lagrangian control on turbulent convection. A
real-world example would be a cloud of droplets moving along with an updraft –
the droplet remains uniformly hotter than the surroundings due to condensation
of water onto its surface and similarly, falling cloud droplets constantly
lose water to the atmosphere thus remaining cooler while moving downward.

Our study also opens several further interesting avenues for investigation
including -but not limited to- the formulation of similar protocols where the
properties of the suspended particles is optimized by a data-driven approach
to attain complex controls and modulation of fluid flows. It would also be
very interesting to explore other protocols where the particle properties are
coupled with the underlying dynamics of the fluid local to the particle, for
example where the temperature of the particle is proportional to its vertical
velocity as well as cases where the momentum exchange of the particle with the
fluid is also considered. Whether other protocols could also exhibit a
transition like the one observed in this study and the physical reasons for
the presence or absence of a transition may also be a question for future
researchers.

Other avenues for further research include a dynamic comparison between the
convective system described in this study driven by Lagrangian forcing with
the Rayleigh-Bénard convection as well as a statistical comparison between the
trajectories of the active tracers with the trajectories of tracers suspended
in a RB convective system. These could provide new insights into the behaviour
of thermal convection.

## Funding

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodowska-Curie grant
agreement No 765048. This work was supported by the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (Grant Agreement No. 882340).

Declaration of interests The authors report no conflict of interest

## Data Availability Statement

Data available on request from the authors – The data that support the
findings of this study are available from the corresponding author upon
reasonable request.

##  Appendix A Numerical Methods

###  A.1 Lattice Boltzmann Method

The fluid equations are solved by the Lattice Boltzmann method with two sets
of populations using a standard D2Q9 grid.

$\displaystyle f_{i}(\boldsymbol{r}+\boldsymbol{c}_{i}\Delta t,t+\Delta t)$
$\displaystyle=f_{i}(\boldsymbol{r},t)-\frac{f_{i}-f^{\rm eq}}{\tau_{f}}\Delta
t+S_{i}\Delta t,$ (26) $\displaystyle
g_{i}(\boldsymbol{r}+\boldsymbol{c}_{i}\Delta t,t+\Delta t)$
$\displaystyle=g_{i}(\boldsymbol{r},t)-\frac{g_{i}-g^{\rm eq}}{\tau_{g}}\Delta
t+q_{i}\Delta t.$ (27)

The evolution of the two sets of populations $f$ and $g$, representing the
fluid and the thermal phase respectively, follow the Lattice Boltzmann
equations with a Bhatnagar–Gross–Krook (BGK) collision operator. The vectors
$\boldsymbol{c}_{i}$ for $i=1,\dots,9$ are the discrete particle velocities,
$\Delta t$ is the lattice time-step, so that $\boldsymbol{c}_{i}\Delta t$ go
from each lattice point to the 8 nearest neighbouring lattice points in the
uniform 2D grid and $\boldsymbol{c}_{0}=0$. $S_{i}$ and $q_{i}$ represent the
momentum forcing (buoyancy) and the thermal forcing respectively. The time-
step and the grid spacing respectively $\Delta t=\Delta r=1$, as is the
standard practice. $f^{\rm eq}$ and $g^{\rm eq}$ are the equilibrium
population distributions as defined in He et al. (1998) given by

$\displaystyle f^{\rm eq}$
$\displaystyle=w_{i}\rho\Bigl{(}1+\frac{\boldsymbol{u}\cdot\boldsymbol{c}_{i}}{c^{2}_{s}}+\frac{(\boldsymbol{u}\cdot\boldsymbol{c}_{i})^{2}}{2c^{4}_{s}}-\frac{\boldsymbol{u}\cdot\boldsymbol{u}}{2c^{2}_{s}}\Bigr{)},$
(28) $\displaystyle g^{\rm eq}$
$\displaystyle=w_{i}T\Bigl{(}1+\frac{\boldsymbol{u}\cdot\boldsymbol{c}_{i}}{c^{2}_{s}}+\frac{(\boldsymbol{u}\cdot\boldsymbol{c}_{i})^{2}}{2c^{4}_{s}}-\frac{\boldsymbol{u}\cdot\boldsymbol{u}}{2c^{2}_{s}}\Bigr{)},$
(29)

where $w_{i}$ are the weights for each population set by the grid used, D2Q9
in this study. $c_{s}$ is the lattice speed of sound set by the choice of
$\boldsymbol{c}_{i}$.

$\tau_{f}$ and $\tau_{g}$ are respectively the fluid and the thermal
relaxation times which set the values for kinematic viscosity $\nu$ and
thermal conductivity $\kappa$ as

$\displaystyle\nu$ $\displaystyle=c_{s}^{2}(\tau_{f}-0.5),$ (30)
$\displaystyle\kappa$ $\displaystyle=c_{s}^{2}(\tau_{g}-0.5).$ (31)

To account for the buoyancy force term, the Guo-forcing scheme (Guo et al.,
2002) is employed with

$S_{i}=\Bigl{(}1-\frac{\Delta
t}{2\tau_{f}}\Bigr{)}w_{i}\Bigl{(}\frac{\boldsymbol{c_{i}}-\boldsymbol{u}}{c^{2}_{s}}+\frac{(\boldsymbol{c_{i}}\cdot\boldsymbol{u})\boldsymbol{c_{i}}}{c^{4}_{s}}\Bigr{)}\cdot\boldsymbol{F},$
(32)

where $\boldsymbol{F}$ is the physical force vector.

The fluid hydrodynamic quantities at each point in space and time are obtained
from the various moments of the populations as

$\displaystyle\rho$ $\displaystyle=\sum_{i}f_{i},$ (33)
$\displaystyle\boldsymbol{u}$
$\displaystyle=\frac{1}{\rho}\sum_{i}f_{i}\boldsymbol{c}_{i}+\frac{\boldsymbol{F}}{2\rho}.$
(34)

The ease of implementation of the Guo-forcing scheme is from the fact that the
velocity $\boldsymbol{u}$ that enters the expression for $f^{\rm eq}$ in
equation (28) is the same as the hydrodynamic velocity obtained in equation
(34). This isn’t the case for other forcing schemes.

The addition of a heat source term (thermal forcing term) is performed
according to (Seta, 2013) with

$q_{i}=\Bigl{(}1-\frac{1}{2\tau_{g}}\Bigr{)}w_{i}\Phi\Delta t,$ (35)

where $\Phi=-\alpha(T-T_{p})$ is the required source term. The temperature is
then obtained at each lattice grid point from the thermal populations $g_{i}$
as

$T=\sum_{i}g_{i}+\Bigl{(}1-\frac{1}{2\tau_{g}}\Bigr{)}\Phi.$ (36)

The no-slip boundary condition for the velocity at the top and bottom walls
are imposed using the bounce-back method (Ladd, 1994). The adiabatic boundary
condition for the top and bottom walls are imposed using the Inamuro method
for setting the normal flux at a boundary for an advected scalar in a fluid
(Yoshino & Inamuro, 2003) by setting the flux equal to $0$.

##  Appendix B Effects of varying $\alpha_{0}$ and $c$

It is clear from the main text that an increase in the number of particles
$N_{p}$ strongly pushes the system towards the convective configuration while
increasing $T_{+}$ weakly causes the system to tend towards stability. The
other ways a phase change from a stable configuration to a convective
configuration can be triggered is by increasing the fluid-particle coupling
strength $\alpha_{0}$ or the size of the particle $c$, both of which serve to
increase the typical velocity $u_{0}$.

(a)

(b)

Figure 12: (a) Normalised vertical temperature profiles for $T_{+}=0.01$,
$N_{p}=180$ for different $\alpha$. The red curves correspond to convective
flows while the blue curves represent the stable flows. (b) Normalised TKE for
$T_{+}=0.02$ plotted against $N_{p}$ for 3 values of $\alpha_{0}$. Horizontal
red line represents $E^{0}_{k}=0.225$.  Figure 13: Normalised TKE for varying
virtual particle size $c$ for two different $\alpha_{0}$.

The former effect can be gauged in figure 12. In panel (a), we see the
behavior of the temperature profile for varying $\alpha_{0}$. It has already
been seen that the stable regime is characterised by a strongly stable
temperature profile while the convective regime is characterised by a weakly
stable temperature gradient. The temperature profile remains nearly identical
for changing values of $\alpha_{0}$ except when the flow changes from stable
(blue curves) to convective (red curves). We also note that the time taken for
the flow to relax from the initial unstable configuration (see equation (11))
to the eventual stationary state is larger for smaller $\alpha_{0}$. It
indicates that for a given temperature scale $T_{+}$ and $N_{p}$, there exists
a temperature difference $\overline{\Delta T}$ for which the flow is stable
independent of $\alpha_{0}$. Panel (b) of the same figure where we plot the
average normalised TKE $\overline{E}_{k}$ shows the transition from stable to
convective for 3 different $\alpha_{0}$. That the increase in TKE corresponds
to the transition from stable to convective was verified from visualisations
of the flow field as well as the strength of the large-scale circulation as
already discussed in Section 3.2. We see that decreasing $\alpha_{0}$
increases the $N_{p}$ of the transition and still note that the empirical
value of $E_{k}^{0}$ for the transition holds.

Increasing $c$ too shows a similar effect, as clear in figure 13 where keeping
the other parameters fixed, a transition to convective configuration is
triggered by enlarging the size of the individual virtual particle.

##  Appendix C Thermal Dissipation

We define the thermal dissipation rate as standard in the turbulence
literature as

$\epsilon_{T}\equiv\kappa\bigl{\langle}(\partial_{i}T(\boldsymbol{x},t))^{2}\bigr{\rangle}_{V},$
(37)

and note that in the statistically stationary regime, the thermal dissipation
is equal to the thermal injection. We have the heat equation given by

$\partial_{t}T+\boldsymbol{u}\cdot\boldsymbol{\nabla}T=\kappa\nabla^{2}T+\Phi.$
(38)

Following (Siggia, 1994) and as shown explicitly by Ching (2014, pp. 5-7) for
the Rayleigh-Bénard convection, we multiply equation (38) with $T$ and average
over the entire domain and time to give

$\frac{1}{2}\frac{d\langle
T^{2}\rangle_{V,t}}{dt}+\frac{1}{2}\bigl{\langle}\boldsymbol{u}\cdot\boldsymbol{\nabla}(T^{2})\bigr{\rangle}_{V,t}-\bigl{\langle}\Phi
T\bigr{\rangle}_{V,t}\\\
=\kappa\bigl{\langle}T\nabla^{2}T\bigr{\rangle}_{V,t}=\kappa\bigl{\langle}\boldsymbol{\nabla}\cdot(T\boldsymbol{\nabla}T)\bigr{\rangle}_{V,t}-\kappa\bigl{\langle}|\boldsymbol{\nabla}T|^{2}\bigr{\rangle}_{V,t},$
(39)

and then use the stationary condition
($\partial_{t}\langle\cdot\rangle_{V,t}=0$) and the incompressibility
($\boldsymbol{\nabla}\cdot\boldsymbol{u}=0$) condition to give

$\bigl{\langle}\boldsymbol{u}\cdot\boldsymbol{\nabla}(T^{2})\bigr{\rangle}_{V}=\bigl{\langle}\boldsymbol{\nabla}\cdot(\boldsymbol{u}T^{2})\bigr{\rangle}_{V}=0.$
(40)

Then, equation (39) becomes

$\kappa\bigl{\langle}|\boldsymbol{\nabla}T|^{2}\bigr{\rangle}_{V,t}=\kappa\bigl{\langle}\boldsymbol{\nabla}\cdot(T\boldsymbol{\nabla}T)\bigr{\rangle}_{V,t}+\bigl{\langle}\Phi
T\bigr{\rangle}_{V,t},$ (41)

or

$\epsilon_{T}=\kappa\bigl{\langle}\boldsymbol{\nabla}\cdot(T\boldsymbol{\nabla}T)\bigr{\rangle}_{V,t}+\bigl{\langle}\Phi
T\bigr{\rangle}_{V,t}$ (42)

The first term of $\epsilon_{T}$ can further be simplified using the Gauss
theorem and writing it in terms of a surface integral

$\kappa\bigl{\langle}\boldsymbol{\nabla}\cdot(T\boldsymbol{\nabla}T)\bigr{\rangle}_{V,t}=\frac{\kappa}{L_{z}}\biggl{[}\Bigl{\langle}T\partial_{z}T\Bigr{\rangle}_{z=L_{z}}-\Bigl{\langle}T\partial_{z}T\Bigr{\rangle}_{z=0}\biggr{]}.$
(43)

In this study, we set $\partial_{z}T=0$ at $z=0$ and $z=L_{z}$. Thus, finally
we get simply

$\epsilon_{T}=\langle\Phi T\rangle_{V,t}.$ (44)

## References

  * Ahlers et al. (2009) Ahlers, Guenter, Grossmann, Siegfried & Lohse, Detlef 2009 Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection. Reviews of Modern Physics 81 (2), 503.
  * Bec et al. (2014) Bec, Jérémie, Homann, Holger & Ray, Samriddhi Sankar 2014 Gravity-driven enhancement of heavy particle clustering in turbulent flow. Physical Review Letters 112 (18), 184501.
  * Beintema et al. (2020) Beintema, Gerben, Corbetta, Alessandro, Biferale, Luca & Toschi, Federico 2020 Controlling Rayleigh-Bénard convection via reinforcement learning. Journal of Turbulence 21 (9-10), 585–605.
  * Carbone et al. (2019) Carbone, Maurizio, Bragg, Andrew D & Iovieno, Michele 2019 Multiscale fluid-particle thermal interaction in isotropic turbulence. Journal of Fluid Mechanics 881, 679–721.
  * Cencini et al. (2006) Cencini, M, Bec, J, Biferale, L, Boffetta, Guido, Celani, A, Lanotte, AS, Musacchio, S & Toschi, F 2006 Dynamics and statistics of heavy particles in turbulent flows. Journal of Turbulence 7, N36.
  * Ching (2014) Ching, Emily SC 2014 Statistics and scaling in turbulent Rayleigh-Bénard convection. Springer.
  * Elperin et al. (1996) Elperin, Tov, Kleeorin, Nathan & Rogachevskii, Igor 1996 Turbulent thermal diffusion of small inertial particles. Physical Review Letters 76 (2), 224.
  * Falkovich et al. (2002) Falkovich, Gregory, Fouxon, A & Stepanov, MG 2002 Acceleration of rain initiation by cloud turbulence. Nature 419 (6903), 151–154.
  * Fernando et al. (2010) Fernando, Harindra JS, Zajic, Dragan, Di Sabatino, Silvana, Dimitrova, Reneta, Hedquist, Brent & Dallman, Ann 2010 Flow, turbulence, and pollutant dispersion in urban atmospheres. Physics of Fluids 22 (5), 051301.
  * Gayen et al. (2014) Gayen, Bishakhdatta, Griffiths, Ross W & Hughes, Graham O 2014 Stability transitions and turbulence in horizontal convection. Journal of Fluid Mechanics 751, 698–724.
  * Guo et al. (2002) Guo, Zhaoli, Zheng, Chuguang & Shi, Baochang 2002 Discrete lattice effects on the forcing term in the lattice Boltzmann method. Physical Review E 65, 046308.
  * He et al. (1998) He, Xiaoyi, Chen, Shiyi & Doolen, Gary D. 1998 A novel thermal model for the lattice Boltzmann method in incompressible limit. Journal of Computational Physics 146 (1), 282 – 300.
  * Irannejad et al. (2015) Irannejad, Abolfazl, Banaeizadeh, Araz & Jaberi, Farhad 2015 Large eddy simulation of turbulent spray combustion. Combustion and Flame 162 (2), 431–450.
  * Kim et al. (2013) Kim, Hojong, Boysen, Dane A, Newhouse, Jocelyn M, Spatocco, Brian L, Chung, Brice, Burke, Paul J, Bradwell, David J, Jiang, Kai, Tomaszowska, Alina A, Wang, Kangli & others 2013 Liquid metal batteries: past, present, and future. Chemical Reviews 113 (3), 2075–2099.
  * Ladd (1994) Ladd, Anthony JC 1994 Numerical simulations of particulate suspensions via a discretized boltzmann equation. part 1. theoretical foundation. Journal of fluid mechanics 271, 285–309.
  * Lay et al. (2008) Lay, Thorne, Hernlund, John & Buffett, Bruce A 2008 Core-mantle boundary heat flow. Nature Geoscience 1 (1), 25–32.
  * Limare et al. (2015) Limare, Angela, Vilella, Kenny, Di Giuseppe, Erika, Farnetani, CG, Kaminski, E, Surducan, Emanoil, Surducan, Vasile, Neamtu, Camelia, Fourel, Loïc & Jaupart, Claude 2015 Microwave-heating laboratory experiments for planetary mantle convection. Journal of Fluid Mechanics 777, 50–67.
  * Lohse & Xia (2010) Lohse, Detlef & Xia, Ke-Qing 2010 Small-scale properties of turbulent rayleigh-bénard convection. Annual Review of Fluid Mechanics 42, 335–364.
  * Markowski (2007) Markowski, Paul 2007 An overview of atmospheric convection. Atmospheric Convection: Research and Operational Forecasting Aspects pp. 1–6.
  * Marshall & Schott (1999) Marshall, John & Schott, Friedrich 1999 Open-ocean convection: Observations, theory, and models. Reviews of Geophysics 37 (1), 1–64.
  * Mazin (1999) Mazin, Ilia 1999 The effect of condensation and evaporation on turbulence in clouds. Atmospheric Research 51 (2), 171–174.
  * Park et al. (2018) Park, Hyungwon John, O’Keefe, Kevin & Richter, David H 2018 Rayleigh-Bénard turbulence modified by two-way coupled inertial, nonisothermal particles. Physical Review Fluids 3 (3), 034307.
  * Salesky & Anderson (2018) Salesky, ST & Anderson, W 2018 Buoyancy effects on large-scale motions in convective atmospheric boundary layers: implications for modulation of near-wall processes. Journal of Fluid Mechanics 856, 135–168.
  * Seta (2013) Seta, Takeshi 2013 Implicit temperature-correction-based immersed-boundary thermal lattice Boltzmann method for the simulation of natural convection. Physical Review E 87, 063304.
  * Siggia (1994) Siggia, Eric D 1994 High Rayleigh number convection. Annual Review of Fluid Mechanics 26 (1), 137–168.
  * Squires & Eaton (1991) Squires, Kyle D & Eaton, John K 1991 Preferential concentration of particles by turbulence. Physics of Fluids A: Fluid Dynamics 3 (5), 1169–1178.
  * Tang & Bau (1994) Tang, Jie & Bau, Haim H 1994 Stabilization of the no-motion state in the Rayleigh-Bénard problem. Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 447 (1931), 587–607.
  * Tritton (1975) Tritton, DJ 1975 Internally heated convection in the atmosphere of venus and in the laboratory. Nature 257 (5522), 110–112.
  * Wang et al. (2021) Wang, Qi, Lohse, Detlef & Shishkina, Olga 2021 Scaling in internally heated convection: a unifying theory. Geophysical Research Letters 48 (4), e2020GL091198.
  * Xi et al. (2016) Xi, Heng-Dong, Zhang, Yi-Bao, Hao, Jian-Tao & Xia, Ke-Qing 2016 Higher-order flow modes in turbulent Rayleigh-Bénard convection. Journal of Fluid Mechanics 805, 31–51.
  * Yang & Shy (2005) Yang, TS & Shy, SS 2005 Two-way interaction between solid particles and homogeneous air turbulence: particle settling rate and turbulence modification measurements. Journal of Fluid Mechanics 526, 171–216.
  * Yoshino & Inamuro (2003) Yoshino, Masato & Inamuro, Takaji 2003 Lattice boltzmann simulations for flow and heat/mass transfer problems in a three-dimensional porous structure. International journal for numerical methods in fluids 43 (2), 183–198.
  * Zamansky et al. (2014) Zamansky, Rémi, Coletti, Filippo, Massot, Marc & Mani, Ali 2014 Radiation induces turbulence in particle-laden fluids. Physics of Fluids 26 (7), 071701.
  * Zamansky et al. (2016) Zamansky, Rémi, Coletti, Filippo, Massot, Marc & Mani, Ali 2016 Turbulent thermal convection driven by heated inertial particles. Journal of Fluid Mechanics 809, 390–437.


# A New Kind of Adversarial Example

Ali Borji
Quintic AI, San Francisco, CA
<EMAIL_ADDRESS>
###### Abstract

Almost all adversarial attacks are formulated to add an imperceptible
perturbation to an image in order to fool a model. Here, we consider the
opposite which is adversarial examples that can fool a human but not a model.
A large enough and perceptible perturbation is added to an image such that a
model maintains its original decision, whereas a human will most likely make a
mistake if forced to decide (or opt not to decide at all). Existing targeted
attacks can be reformulated to synthesize such adversarial examples. Our
proposed attack, dubbed NKE, is similar in essence to the fooling images, but
is more efficient since it uses gradient descent instead of evolutionary
algorithms. It also offers a new and unified perspective into the problem of
adversarial vulnerability. Experimental results over MNIST and CIFAR-10
datasets show that our attack is quite efficient in fooling deep neural
networks. Code is available at https://github.com/aliborji/NKE.

##  1 Introduction

Figure 1: Sample adversarial examples generated by our attack for MNIST and
CIFAR-10 images. The original images (top row) are highly perturbed (bottom
row) such that the model still maintains its decision, whereas humans can no
longer make a decision (or will make an error if forced to decide).

Deep neural networks are very brittle [9, 4, 1]. A small imperceptible
perturbation is enough to fool them. A sample $x^{\prime}$ is said to be an
adversarial example for $x$ when $x^{\prime}$ is close to $x$ under a specific
distance metric, while $f(x^{\prime})\neq y$. Formally:

$x^{\prime}:D(x,x^{\prime})<\epsilon_{s}\ \ \text{s.t}\ \ f(x^{\prime})\neq y$
(1)

where $D(.,.)$ is a distance metric, $\epsilon_{s}$ is a predefined distance
constraint (_a.k.a._ allowed perturbation), $f(.)$ is a neural network, and
$y$ is the true label of sample $x$.

The bulk of research on adversarial attack has focused on crafting adversarial
examples by adding imperceptible perturbations to images. Some works have
shown that perceptible perturbations can fool neural networks in the physical
work (_e.g._ [3, 8, 7]). Here, we introduce a new attack which is in essence
similar to the existing ones but differs in an important aspect. Instead of
fooling a model to change its decision on the perturbed image (small
imperceptible perturbation; images still recognizable by humans), it aims to
fool the model to keep its original decision when there is almost no signal in
the image (large perceptible perturbation; image unrecognizable by humans).
Our attack is similar to the fooling images [7], but instead of evolutionary
techniques here we employ gradient descent which is much more efficient (See
Fig. 2). In a sense, our work gives a new perspective to unify adversarial
examples and fooling images under one umbrella. Gradient descent can be used
to generate both types of perturbations to fool a model111Except images with
patterns and textures similar to the images in the target class. See the left
images in Fig. 2.. Adversarial images synthesized by our attack for some
images from MNIST and CIFAR-10 datasets are shown in Fig. 1. A model
confidently maintains its original decision on the perturbed images that look
just like noise to humans.

Highly perturbed images with the predicted label as the original label are
dangerous and equally important as adversarial examples with imperceptible
perturbations. As a real-world scenario, imagine you see a green light from a
distance while driving. As you get closer to the light, it starts to rain such
that you are no longer certain about the state of the light. In this
situation, you may opt to slow down or pull over. Imagine a self-driving
vehicle in the same situation, except now an adversary is manipulating the
image of the light. It would be very dangerous if the car maintains its
_green_ decision with high confidence, while the light may actually be red.
Here, the adversary aims to maximize the chance of fooling the system by
adding a highly perceptible perturbation.

Figure 2: Example meaningless patterns (_a.k.a._ fooling images) that are
classified as familiar objects by a deep neural network [7]. Figure 3: A
general view of adversarial robustness. Hypothetically, a system is fully
robust if it can function perfectly regardless of the perturbation magnitude.
In reality, however, the performance curves look like a reverse sigmoid
function (_a.k.a._ psychometric functions for humans) indicating that
performance drops as the perturbation grows (panel b). a) Current adversarial
attacks are formulated to lower the left side of the curve (_i.e._ zero
accuracy for all values of $\epsilon$) and emphasize on imperceptible
perturbations within the yellow $L_{p}$ ball (Eq. 1). In contrast, our attack
aims to lift the right side of the curve (_i.e._ 100% accuracy for all values
of $\epsilon$). These two types of attacks can also be combined, _i.e._
lowering the left side and lifting the right side like a sigmoid. We emphasize
on large perturbations (gray region) for which humans are more likely to fail.
b) Hypothetical performance curves of humans vs. two models, one under
performing humans (model A), and the other outperforming humans (model B).
Here, model B is less robust than model A since its psychometric function
falls below the psychometric function of model A.

##  2 Attack formulation

We name our attack NKE after the “new kind of adversarial example” in the
paper title. Adversarial examples generated by NKE fulfill the following
condition:

$x^{\prime}:\Bigl{[}g(x^{\prime})\neq g(x)\ \ \land\ \
f(x^{\prime})=g(x)\Bigr{]}$ (2)

where $g(x)$ is the oracle ($y$ in Eq. 1). Here, we assume that the initial
prediction for sample $x$ is correct (_i.e._ $f(x)=g(x)$). According to this
formulation, an input is adversarial if the model predicts the same label as
the original label while humans222human expert or the best possible classifier
in a Bayesian sense. will make a mistake if forced to make a decision (or it
may opt not to make a decision at all). Since the function $g(x)$ is usually
not available, a distance metric can be used as a proxy. We revise the Eq. 1
to reflect the above definition as follows:

$x^{\prime}:\Bigl{[}D(x,x^{\prime})>\epsilon_{l}\ \ \land\ \
f(x^{\prime})=y\Bigr{]}$ (3)

This new definition mirrors the classic definition. Instead of small
imperceptible perturbations, it emphasizes highly perceptible perturbations
(beyond a certain threshold $\epsilon_{l}$).

Fig.3 depicts a general view on adversarial robustness. A robust system should
be insensitive to small perturbations while being sensitive to large
perturbations. It should become more sensitive as the perturbation grows. Two
perturbation bounds are particularly interesting. The first one,
$\epsilon_{s}$, is considered in the classic formulation in Eq. 1. It is used
to constrain the amount of the added perturbation (_i.e.__maximum allowed
perturbation_ _e.g._ 8/255)333The intuition here is that humans will not make
a mistake at small perturbations, thus $g(x^{\prime})=g(x)$.. The second one,
$\epsilon_{l}$, determines the _minimum allowed perturbation_ , and is chosen
such that the perturbations beyond it will enforce a human error. The
perturbed image would be considered adversarial if the model maintains its
prior decision. The adversarial examples generated this way are reminiscent of
fooling images (Fig. 2). Our attack also resembles deep
dream444https://en.wikipedia.org/wiki/DeepDream. Instead of altering the image
to be classified in another class, we want it to be classified in the correct
class, even with higher confidence.

The formulation in Eq. 3 is basically a targeted attack in which the desired
target class is actually the true class $y$. To increase the confidence of the
model in choosing the label $y$ for image $x$, we can reduce the loss $J(x,y)$
by performing gradient descent on the image. Any existing attack can be
reformulated to serve this purpose. Here, we use the iterative FGSM attack [5,
4] (_a.k.a._ IFGSM; commonly known as BIM) with $L_{\infty}$-norm:

$\displaystyle x^{\prime}_{0}$ $\displaystyle=$ $\displaystyle x,$
$\displaystyle x^{\prime}_{N+1}$ $\displaystyle=$
$\displaystyle\text{clip}_{x,\epsilon}\Bigl{\\{}x^{\prime}_{N}-\epsilon\operatorname{sign}\bigl{(}\nabla_{x}J(x^{\prime}_{N},y)\bigr{)}\Bigr{\\}}$
(4)

##  3 Experiments and results

Figure 4: Attack performance. The lower accuracy means higher efficiency for
the BIM attack. The opposite happens for the NKE attack. The $x$ axis shows
the magnitude of perturbation. Please see also Fig. 5 for sample images and
their perturbed versions.

On MNIST, we train a CNN with two convolutional layers for 10 epochs. On
CIFAR-10, we train a ResNet50 for 30 epochs. We do not use any data
augmentation. Images are scaled to [0,1] range. The classifiers score about
95% and 61% test accuracy, respectively.

To analyze the attack performance, we only consider the images that were
classified correctly. We then measure the fraction of them that retain their
labels after perturbation (_i.e._ classified correctly again). Results are
shown in Fig. 4 and Fig. 5. The NKE attack is less effective as the
perturbation grows. This is because there is more signal at small $\epsilon$
values, making it easy to alter the input to be classified correctly.
Conversely, at higher perturbations when the signal fades, it becomes
increasingly harder to alter the images such that its new label matchs the
initial prediction. The NKE attack becomes more effective with more steps.
With five steps, it is almost always possible to perturb the image such that
it can be predicted correctly regardless of the magnitude of the signal in the
image. Results are consistent across both datasets.

At $\epsilon=0.4$, the BIM attack with only one step completely punctures the
MNIST classifier. Up to around $\epsilon_{s}=0.1$, perturbations generated by
this attack are imperceptible. $\epsilon_{s}$ decreases with more steps of
BIM. A similar observation can be made for the NKE attack. The image becomes
almost completely unrecognizable at $\epsilon_{l}=0.9$ with one step of NKE.
The classifier at this $\epsilon$ achieves around 75% accuracy (_i.e._ 75% vs.
0% by humans!). The NKE attack with five steps scores above 99% at
$\epsilon_{l}=0.9$ perturbation. Similar qualitative results are obtained on
the CIFAR-10 dataset. Here, $\epsilon_{l}$ is about $0.4$ and $0.1$
corresponding to one and five steps of NKE, respectively.

##  4 Conclusion and future work

We introduced a new attack to generate a new kind of adversarial example.
These adversarial images are similar to the fooling images but are much easier
and faster to synthesize. Our work offers a new perspective in viewing
adversarial attacks and unifies methods that generate adversarial examples and
fooling images. Theoretically, any attack can be revised to generate
adversarial examples that can fool humans, but not models. Some may be more
suitable than others for this purpose. For example, the Carlini-Wagner attack
[2] can be revised as follows:

$\displaystyle\text{Maximize}\ \ ||\epsilon||_{2}\ \ s.t$ (5) $\displaystyle
1.\ \ f(x+\epsilon)=f(x)$ $\displaystyle 2.\ \ \ \ x+\epsilon\in[0,1]^{m}$

which can be written as $J(x+\epsilon,f(x))-c\cdot||\epsilon||_{2}\ \ s.t\ \
x+\epsilon\in[0,1]^{m}$ using Lagrange multiplies, and can be optimized via
the Quasi Newton method.

Some questions that should be answered in future research include: a) how
transferable are the adversarial examples generated by NKE attack, b) is it
possible to craft universal perceptible patterns that can make a classifier to
maintain its initial decision for images from different classes (akin to [6]),
and c) how effective is adversarial training (or other defenses) against the
NKE attack. Finally, we hope that our work will spark follow-up research and
help gain deeper insights for solving the problem of adversarial
vulnerability.

## References

  * [1] Ali Borji.  Is current research on adversarial robustness addressing the right problem?  arXiv, 2022.
  * [2] Nicholas Carlini and David A. Wagner.  Towards evaluating the robustness of neural networks.  In IEEE Symposium on Security and Privacy, pages 39–57, 2017.
  * [3] Kevin Eykholt, Ivan Evtimov, Earlence Fernandes, Bo Li, Amir Rahmati, Chaowei Xiao, Atul Prakash, Tadayoshi Kohno, and Dawn Song.  Robust physical-world attacks on deep learning visual classification.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1625–1634, 2018.
  * [4] Ian Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and harnessing adversarial examples.  In Proc. ICLR, 2015.
  * [5] Alexey Kurakin, Ian J. Goodfellow, and Samy Bengio.  Adversarial examples in the physical world.  CoRR, abs/1607.02533, 2016.
  * [6] Seyed-Mohsen Moosavi-Dezfooli, Alhussein Fawzi, Omar Fawzi, and Pascal Frossard.  Universal adversarial perturbations.  In Proc. CVPR, pages 86–94, 2017.
  * [7] Anh Nguyen, Jason Yosinski, and Jeff Clune.  Deep neural networks are easily fooled: High confidence predictions for unrecognizable images.  In Proc. CVPR, 2015.
  * [8] Mahmood Sharif, Sruti Bhagavatula, Lujo Bauer, and Michael K Reiter.  A general framework for adversarial examples with objectives.  ACM Transactions on Privacy and Security (TOPS), 22(3):1–30, 2019\.
  * [9] Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus.  Intriguing properties of neural networks.  In In Proc. ICLR, 2014.

Figure 5: Sample images from MNIST and CIFAR-10 datasets along with the
adversarial examples generated by NKE attack for different values of
$\epsilon$. Columns correspond to 1, 3, and 5 steps of NKE attack. Notice that
$\epsilon_{l}$ drops as the number of steps grows which means that
perturbations become visible sooner with more step of the attack.