# Large-scale influence of numerical noises as artificial stochastic
disturbances on a sustained turbulence

Shijie Qin2  Shijun Liao1,2<EMAIL_ADDRESS>1State Key Laboratory of Ocean
Engineering, Shanghai 200240, China 2 Center of Marine Numerical Experiment,
School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong
University, Shanghai 200240, China

###### Abstract

We investigate the large-scale influence of numerical noises as tiny
artificial stochastic disturbances on a sustained turbulence. Using the two-
dimensional (2D) turbulent Rayleigh-Bénard (RB) convection as an example, we
numerically solve the NS equations, separately, by means of a traditional
algorithm with double precision (marked by RKwD) and the so-called clean
numerical simulation (CNS). The numerical simulation given by the RKwD is a
mixture of the “true” physical solution and the “false” numerical noises that
is random and can be regarded as a kind of artificial stochastic disturbances:
unfortunately, the “true” physical solution is mostly at the same level as the
“false” numerical noises. By contrast, the CNS can greatly reduce the
background numerical noise to any a required level so that the “false”
numerical noises are negligible compared with the “true” physical solution and
thus the CNS solution can be used as a “clean” benchmark solution for
comparison. It is found that the numerical noises as tiny artificial
stochastic disturbances could indeed lead to large-scale deviations of
simulations not only in spatio-temporal trajectories but also even in
statistics. Especially, these numerical noises (as artificial stochastic
disturbances) even lead to different types of flows: the shearing convection
occurs for the RKwD simulations, and its corresponding flow field turns to a
kind of zonal flow thereafter, however the CNS benchmark solution always
sustains the non-shearing vortical/roll-like convection during the whole
process of simulation. Thus, we provide a rigorous evidence that numerical
noises as a kind of small-scale artificial stochastic disturbances have
quantitatively and qualitatively large-scale influences on a sustained
turbulence, i.e. the 2D turbulent RB convection considered in this paper.

Keyword Rayleigh-Bénard convection, stochastic disturbances, NS equations,
numerical noises

##  1 Introduction

Turbulent flows of a compressible viscous fluid are mostly described by the
Navier-Stokes (NS) equations in the conservation form

$\frac{\partial\mathbf{U}}{\partial
t}+\nabla\cdot\mathbf{F}_{a}=\nabla\cdot\mathbf{F}_{d},$ (1)

subject to a given proper boundary condition and the initial condition

$\mathbf{U}|_{t=0}=\mathbf{U}_{0}(\mathbf{r}),\hskip
28.45274pt\mathbf{r}\in\Omega\hskip 0.28453pt,$ (2)

where $\mathbf{U}(\mathbf{r},t)=(\rho,\,\rho\hskip 0.85358pt\mathbf{u}\hskip
0.28453pt,\,E)^{\mathbf{T}}$ is a vector of conserved variables, $\nabla$ is
the Hamiltonian/gradient operator,

$\mathbf{F}_{a}=\begin{pmatrix}\rho\hskip 0.85358pt\mathbf{u}\\\ \rho\hskip
0.85358pt\mathbf{u}\otimes\mathbf{u}+P\hskip 0.85358pt\mathbf{I}\\\
(E+P)\hskip 0.85358pt\mathbf{u}\end{pmatrix},\hskip
28.45274pt\mathbf{F}_{b}=\begin{pmatrix}0\\\ \bm{\tau}\\\
\mathbf{u}\cdot\bm{\tau}+\mathbf{q}\end{pmatrix}$ (3)

are the advection/hyperbolic and diffusion fluxes, $\mathbf{r}\in\Omega$ is
the vector of spatial position, $t$ denotes the time,
$\mathbf{U}_{0}(\mathbf{r})$ is a given vector, $\rho$ is the density of mass,
$\mathbf{u}$ is the velocity vector, $\mathbf{u}\otimes\mathbf{u}$ denotes the
tensor product of $\mathbf{u}$ and $\mathbf{u}$, $P$ is the pressure,
$\mathbf{I}$ denotes the unit tensor, $E$ is the total energy per unit mass,
$\bm{\tau}$ denotes the tensor of viscous stress and $\mathbf{q}$ is the
vector of heat diffusion flux, respectively.

Note that tiny stochastic disturbances in fluid flows resulting from either
small-scale thermal fluctuations or environmental noises are unavoidable in
practice and the triggered deviations might exponentially increase (Leith &
Kraichnan, 1972; Boffetta & Musacchio, 2001, 2017). However, these tiny
stochastic disturbances are not considered in the above-mentioned NS equations
(1)-(3) that are deterministic. Landau & Lifshitz (1959) proposed a stochastic
form of the NS equations that models the effect of thermal fluctuations via an
additional stochastic stress tensor (Graham, 1974; Swift & Hohenberg, 1977;
Bell et al., 2007; Donev et al., 2010, 2014). However, the direct numerical
simulation (DNS) of this stochastic model is rather difficult due to the
extremely fine resolution that is required to accurately measure the
velocities at dissipation-range length scales. By contrast, molecular dynamics
(MD) provides molecular-level simulation techniques (Bird, 1998; Donev et al.,
2011; Smith, 2015; McMullen et al., 2022) for directly investigating the role
of thermal fluctuations in turbulent flows.

Currently, McMullen et al. (2022) investigated a decaying turbulent flow and
found that, due to thermal fluctuations, the molecular-gas-dynamics spectra
grow quadratically with wave number in the dissipation range, while the
Navier-Stokes spectra decay exponentially. Furthermore, the transition to
quadratic growth occurs at a length scale much larger than the gas molecular
mean free path, namely in a regime that the NS equations are widely believed
to describe. Thus, they provided the first direct evidence that “the Navier-
Stokes equations do not describe turbulent gas flows in the dissipation range
because they neglect thermal fluctuations” (McMullen et al., 2022), which is
in agreement with the results given by Bandak et al. (2022), Bell et al.
(2022), Eyink & Jafari (2022) and so on. Gallis et al. (2021) separately used
the direct simulation Monte Carlo (DSMC) method (molecular gas dynamics) and
direct numerical simulation (DNS) of the NS equations to simulate a freely
decaying turbulent flow, i.e. the compressible Taylor-Green vortex flow, and
found that both methods produce basically the same energy decay for the Mach
and Reynolds numbers they examined, but the molecular fluctuations in DSMC
(and in experiments) can break symmetries, which in turn can cause the flows
different from but basically similar to those given by the DNS. Note that the
above-mentioned investigations mainly focus on the influence of tiny
stochastic disturbances resulting from thermal fluctuations on small-scale
properties of freely decaying turbulence. Might the micro-level stochastic
disturbances have huge influences on large-scale properties of some sustained
turbulent flows? This is the primary motivation for our investigation of this
paper.

Some turbulent flows might have multiple states according to Frisch (1986) and
his well-known monograph on turbulence (Frisch, 1995) that “it is typical for
dissipative dynamical systems to have more than one attractor” and “each
attractor has an associated basin”, and thus “the statistical properties of
the solution will then depend on to which basin the initial condition
belongs”. The first evidence of this phenomenon in a turbulent flow is given
by Huisman et al. (2014) for their study on the highly turbulent Taylor-
Couette flow. Note that the above-mentioned work mainly focuses on the
influence of initial data on the flow state of turbulence as well as its
statistical properties. Considering that there are some previous
investigations (Knobloch & Weiss, 1989; Kraut et al., 1999; Masoller, 2002; de
Souza et al., 2007) in the effects of noises (in general whose levels are much
larger than the stochastic disturbances considered later in this paper) on the
dynamical systems with multiple attractors, in this paper we aim to
demonstrate that not only the initial conditions but also the weak, small-
scale stochastic noises/disturbances can determine in which basin the solution
of a sustained turbulent flow will reside for a long time.

In theory, the influence of tiny stochastic disturbances on turbulent flows is
fundamental, since this kind of stochastic disturbances is often in a micro-
level of magnitude and is unavoidable in practice. There are two kinds of
stochastic disturbances: natural and artificial. The thermal fluctuation
belongs to the former, whereas the latter can be caused by many sources. Note
that the background numerical noises, i.e. truncation errors and round-off
errors, always exist for all numerical algorithms. Besides, it is widely
believed that turbulence has a close relationship to chaos and thus the
background numerical noises should increase exponentially (and quickly) until
to the same level of “true” physical solution. Therefore, a computer-generated
simulation of the NS equations is a kind of mixture of the “true” physical
solution and “false” numerical noises. Note that the background numerical
noise is tiny and random, dependent on different numerical algorithms and data
accuracy, which itself is a kind of artificial stochastic disturbance.
Therefore, in a natural way, the background numerical noises can be regarded
as the sum of all kinds of tiny artificial stochastic disturbances. In this
paper, we investigate the large-scale influences of numerical noises as tiny
artificial stochastic disturbances on a sustained turbulence by means of
numerically solving the above-mentioned NS equations via a traditional
numerical algorithm (Runge-Kutta’s method with double precision, marked by
RKwD) and the so-called “Clean Numerical Simulation” (CNS) (Liao, 2009, 2013,
2014; Hu & Liao, 2020; Qin & Liao, 2020; Xu et al., 2021; Li et al., 2021;
Liao et al., 2022; Liao & Qin, 2022), separately, with the same
initial/boundary conditions using the same values of physical parameters. The
result given by the traditional numerical algorithm RKwD is a mixture of the
“true” physical solution and “false” numerical noises, which are mostly at the
same order of magnitude, where the background numerical noise is regarded as
the sum of all tiny artificial stochastic disturbances. By contrast, the CNS
can greatly reduce the background numerical noise to any a required level so
that the numerical noises are negligible compared with the “true” physical
solution, and thus the corresponding numerical result is convergent
(reproducible) in an interval of time long enough for statistics, as described
below. In other words, results given by the CNS can be regarded as a “clean”
benchmark solution. Thus, we can investigate the influences of tiny artificial
stochastic disturbances on turbulent flows by means of comparing the RKwD
simulations with the CNS benchmark solution.

Then, let us briefly talk about the motivation and basic ideas of the CNS.
Strictly speaking, all numerical simulations are not “clean”, since background
numerical noises (i.e. truncation errors and round-off errors) always exist
there. Indeed, for non-chaotic systems, the magnitude of numerical noises can
be controlled on a tiny level much smaller than “true” physical solution so
that the influences of the numerical noises can be neglected. However, for
chaotic dynamical systems (Li & Yorke, 1975; Parker & Chua, 1989; Lorenz,
1993; Peter, 1998; Sprott, 2010), numerical noises increase exponentially due
to the “sensitivity dependence on initial condition” (SDIC), which was first
discovered by Poincaré (1890) and later rediscovered by Lorenz (1963) with a
more famous name “butterfly-effect”: a hurricane happening in North America
might be created by a flapping of the wings of a distant butterfly in South
America several weeks earlier. More importantly, it was further found by
Lorenz (1989, 2006) that a chaotic dynamical system has the sensitivity
dependence not only on initial condition (SDIC) but also on numerical
algorithms (SDNA) in single/double precision. Such kind of uncertainty
certainly raises serious doubt about the reliability of numerical simulations
of chaotic systems. For example, Teixeira et al. (2007) carefully investigated
the time-step sensitivity of three nonlinear atmospheric models by means of
some traditional numerical algorithms (in single/double precision), and made a
rather pessimistic conclusion that “for chaotic systems, numerical convergence
cannot be guaranteed forever”.

To overcome the above-mentioned limitations/restrictions of traditional
numerical algorithms (in single/double precision) for chaotic dynamical
systems, Liao (2009) suggested a numerical strategy, namely the “Clean
Numerical Simulation” (CNS). The basic idea of the CNS (Liao, 2013, 2014; Hu &
Liao, 2020; Qin & Liao, 2020; Li et al., 2021; Liao et al., 2022; Liao & Qin,
2022) is to greatly decrease the background numerical noises, i.e. truncation
errors and round-off errors, to such a tiny level that the influence of
numerical noises can be neglected in an interval of time $0\leq t\leq T_{c}$
that is long enough for statistics, where $T_{c}$ is the so-called “critical
predictable time”. The CNS is based on such a well-known phenomenon: for a
numerical/computer-generated simulation of chaotic dynamical system, the level
of simulation-deviation (in an average meaning) from its (“true”) physical
solution increases exponentially to a macroscopic one (at $t=T_{c}$), i.e.

${\cal E}(t)={\cal E}_{0}\exp(K\hskip 0.85358ptt),\hskip
28.45274ptt\in[0,T_{c}],$ (4)

where $K>0$ is the so-called “noise-growing exponent”, ${\cal E}_{0}$ denotes
the level of background numerical noise, ${\cal E}(t)$ is the level of
simulation-deviation (in an average meaning) from the physical solution,
respectively. In theory, the critical predictable time $T_{c}$ is determined
by a critical level ${\cal E}_{c}$ of simulation-deviation from its physical
solution, say,

$T_{c}=\frac{1}{K}\ln\left(\frac{{\cal E}_{c}}{{\cal E}_{0}}\right).$ (5)

Obviously, for a given critical level ${\cal E}_{c}$, the smaller the level of
the background numerical noise ${\cal E}_{0}$, the larger the critical
predictable time $T_{c}$. This is the reason why in the frame of the CNS we
have to greatly decrease the background numerical noises, i.e. truncation
errors and round-off errors, to a tiny enough level. So, different from the
Taylor series method, the key point of the CNS is the so-called “critical
predictable time” $T_{c}$ that determines a temporal interval $[0,T_{c}]$ in
which the numerical simulations are “reliable” and “clean”, since their
“false” numerical noises are much smaller than the “true” physical solution
and thus are negligible. For more details about the CNS, please refer to Liao
(2009, 2013, 2014) and his coauthors (Hu & Liao, 2020; Qin & Liao, 2020; Xu et
al., 2021; Li et al., 2021; Liao et al., 2022).

The CNS has been successfully applied to many chaotic dynamical systems. For
example, by means of traditional numerical algorithms (in double precision),
one can get convergent (i.e. reproducible) numerical simulations of the famous
Lorenz equations just in a rather short interval of time, say, approximately
$t\in[0,32]$. However, using the CNS, a convergent numerical simulation of the
Lorenz equations was obtained first by Liao (2009) in $t\in[0,1000]$ and then
by Liao & Wang (2014) in a much longer interval of time, i.e. $t\in[0,10000]$.
Besides, since the background numerical noises of the CNS can be much smaller
even than the micro-level physical uncertainty, Lin et al. (2017) successfully
applied the CNS to provide a direct rigorous evidence that the micro-level
thermal fluctuation is the origin of macroscopic randomness of the turbulent
Rayleigh-Bénard convection. Especially, it is worth noting that the CNS has
been successfully applied to find more than $2000$ new families of periodic
orbits of three-body system (Li & Liao, 2017; Li et al., 2018; Li & Liao,
2019; Li et al., 2021; Liao et al., 2022). The discovery of these new periodic
orbits was reported twice in the famous popular magazine New Scientist (Crane,
2017; Whyte, 2018), because only three families of periodic orbits of the
three-body problem had been reported in three hundred years after Newton
mentioned this problem in 1687. All of these illustrate the validity, novelty,
and great potential of the CNS.

Recently, an efficient CNS algorithm has been proposed to solve the spatio-
temporal chaotic systems, i.e. the complex Ginzburg-Landau equation (Hu &
Liao, 2020) and the damped driven sine-Gordon equation (Qin & Liao, 2020),
respectively. Using the CNS result as a benchmark solution, one can
investigate the influence of numerical noises on the computer-generated
simulation of a spatio-temporal chaotic system. It was found (Hu & Liao, 2020;
Qin & Liao, 2020) that numerical noises might lead to huge deviations of
computer-generated simulations of some spatio-temporal chaotic systems not
only in trajectories but also even in statistics.

In this paper, we apply the CNS and a traditional algorithm (based on the 4th-
order Runge-Kutta’s method with double precision, marked by RKwD), separately,
to solve a sustained turbulence, i.e. the two-dimensional (2D) turbulent
Rayleigh-Bénard convection (RBC). Note that, using the CNS, the background
numerical noises can be decreased to such a tiny level that the numerical
noises are negligible in a long enough interval of time, so that the CNS
result can be regarded as a benchmark solution of the “true” physical result
to investigate the influence of tiny artificial stochastic disturbances by
comparing the RKwD simulations with the CNS benchmark solution. In this way,
we provide a rigorous evidence that tiny artificial stochastic disturbances
have huge influences on large-scale properties of the turbulent Rayleigh-
Bénard convection not only in statistics but also even in flow types: the CNS
benchmark solution always keeps the vortical/roll-like turbulent convection,
however, for the RKwD simulations, the shearing convection occurs and its
corresponding flow field turns to zonal flow thereafter. This phenomenon is
reasonable if the boundaries of different attractor basins (mentioned above)
in this multistable system are intricately interwoven as has been observed in
other cases (Shrimali et al., 2008).

##  2 Mathematical model for 2D turbulent RBC

The buoyancy-driven convection in a fluid layer between two horizontal
parallel plates heated from below and cooled from above, known as the
Rayleigh-Bénard convection (RBC) for compressible viscous fluids, is one of
the most fundamental and classic paradigms of nonlinear dynamics in fluid
mechanics. It was first investigated by Rayleigh (1916) and the continuous
efforts devoted to the study of this problem have greatly enriched our
understanding (Chandrasekhar, 1961; Schlüter et al., 1965; Heslot et al.,
1987; Kadanoff, 2001; Niemela & Sreenivasan, 2006; Ahlers et al., 2009; Lohse
& Xia, 2010; Zhou & Xia, 2013; Goluskin et al., 2014; Wang et al., 2020).

As illustrated in Figure 1, a thin layer of fluid is confined between two
horizontal plates that are separated by a distance $H$, where $T_{0}$ and
$T_{0}+\Delta T$ denote the temperatures of the top and bottom boundary
surfaces, $L$ is the horizontal length of computational domain and $g$ is the
gravity acceleration, respectively. The typical vortical/roll-like motions of
two-dimensional (2D) RBC, as illustrated in Figure 1(a), and their
corresponding turbulent states, have been studied extensively by means of the
direct numerical simulation (DNS) (Saltzman, 1962; Fromm, 1965; Veronis, 1968;
Moore & Weiss, 1973; Curry et al., 1984; Zienicke et al., 1998; Johnston &
Doering, 2009; Huang & Zhou, 2013; Zhang et al., 2017; Zhu et al., 2018).
However, there is another type of flow in the 2D turbulent RBC in the case of
the free-slip boundary conditions imposed on two horizontal parallel plates
and the periodic boundary conditions on the left and right sides, namely zonal
flow (Goluskin et al., 2014; van der Pol et al., 2014; von Hardenberg et al.,
2015; Wang et al., 2020), as shown in Figure 1(b). It is worth noting that
such kind of turbulent zonal flow has been widely found in nature and
laboratory such as in the atmosphere of Jupiter (Heimpel et al., 2005; Kaspi
et al., 2018) and some Jovian planets (Sun et al., 1993; Cho & Polvani, 1996;
Yano et al., 2003), in the oceans (Maximenko et al., 2005; Richards et al.,
2006), in the Earth’s outer core (Miyagoshi et al., 2010), in toroidal tokamak
devices (Diamond et al., 2005) and so on. Thus, here we choose the 2D
turbulent RBC with the free-slip boundary conditions at the upper and lower
plates and the periodic boundary conditions in the horizontal direction as our
mathematical model for the sustained turbulence, governed by the NS equations.

---

Figure 1: Schematic drawings of two-dimensional turbulent Rayleigh-Bénard
convections in two totally different flow types: (a) typical vortical/roll-
like flow and (b) zonal flow. The fluid layer between two parallel plates
which are separated by a height $H$ obtains heat from the bottom boundary
surface because of the constant temperature difference $\Delta T>0$, where $L$
is the horizontal length of computational domain and the downward direction of
gravity acceleration $g$ is indicated.

Using the length scale $H$, velocity scale $\sqrt{g\alpha H\Delta T}$ and
temperature scale $\Delta T$ as the characteristic scales, the corresponding
dimensionless NS equations, combined with the Boussinesq approximation
(Saltzman, 1962), read

$\frac{\partial}{\partial
t}\nabla^{2}\psi+\frac{\partial(\psi,\nabla^{2}\psi)}{\partial(x,z)}-\frac{\partial\theta}{\partial
x}-\sqrt{\frac{Pr}{Ra}}\nabla^{4}\psi=0,$ (6) $\frac{\partial\theta}{\partial
t}+\frac{\partial(\psi,\theta)}{\partial(x,z)}-\frac{\partial\psi}{\partial
x}-\frac{1}{\sqrt{PrRa}}\nabla^{2}\theta=0,$ (7)

where $\psi$ is a stream-function with the definition

$u=-\frac{\partial\psi}{\partial z},\hskip
28.45274ptw=\frac{\partial\psi}{\partial x},$ (8)

in which $u$ and $w$ are the horizontal and vertical velocities, $\theta$ is
the temperature departure from a linear variation background, i.e. the
temperature is expressed as $T=\theta-z+1$ in the case of $T_{0}=0$, $t$
denotes the time, $x\in[0,\Gamma]$ and $z\in[0,1]$ are the horizontal and
vertical position coordinates, $\Gamma=L/H$ denotes the aspect ratio,
$\nabla^{2}$ is the Laplace operator and thus
$\nabla^{4}=\nabla^{2}\nabla^{2}$,

$\frac{\partial(a,b)}{\partial(x,z)}=\frac{\partial a}{\partial
x}\frac{\partial b}{\partial z}-\frac{\partial b}{\partial x}\frac{\partial
a}{\partial z}$ (9)

is the Jacobi operator, the Rayleigh number $Ra$ and Prandtl number $Pr$ are
defined by

$Ra=\frac{g\alpha H^{3}\Delta T}{\nu\kappa},\hskip
28.45274ptPr=\frac{\nu}{\kappa},$ (10)

respectively, in which $\alpha$ is the thermal expansion coefficient and
$\nu=\mu/\rho$ is the kinematic viscosity.

Note that the free-slip boundary conditions are adopted at the upper and lower
plates where temperatures are assumed to be constant. Hence, the NS equations
(6)-(7) have the following boundary conditions:

$\psi=\frac{\partial^{2}\psi}{\partial z^{2}}=\theta=0$ (11)

at $z=0$ and $z=1$. On the other hand, since the fluid layer can extend to
infinity in the horizontal direction, we adopt the periodic boundary
conditions for $\psi$ and $\theta$ at the lateral boundaries in the horizontal
direction, i.e. at $x=0$ and $x=\Gamma$.

Without loss of generality, in this paper let us consider the case with the
aspect ratio $\Gamma=L/H=2\sqrt{2}$, which is large enough for the
approximation of heat flux at an infinite aspect ratio (Saltzman, 1962; Curry
et al., 1984; Lin et al., 2017), the Prandtl number $Pr=6.8$ (corresponding to
the water at room temperature 20∘C), and the Rayleigh number $Ra=6.8\times
10^{8}$ corresponding to a turbulent state, respectively. In addition, the
initial temperature and velocity fields are randomly generated by the thermal
fluctuations in Gaussian white noises (Lin et al., 2017), with the temperature
standard deviation $\sigma_{T}=10^{-10}$ and the velocity standard deviation
$\sigma_{u}=10^{-9}$, respectively.

##  3 Influence of numerical noises as artificial stochastic disturbances

The deterministic NS equations (6)-(7) are numerically solved here,
separately, by means of a traditional algorithm (Runge-Kutta’s method with
double precision, marked by RKwD) whose numerical noises are mostly in the
same order of magnitude as the “true” physical solution, and the CNS whose
numerical noises are much less than the physical solution and thus negligible.
By means of comparing the RKwD simulations with the CNS benchmark solution, it
is found that the numerical noises indeed might lead to huge large-scale
differences even in statistics and flow types of the 2D turbulent RBC, as
described below.

First, we apply the clean numerical simulation (CNS) to greatly decrease the
background numerical noises, i.e. the truncation errors and round-off errors,
to such a tiny level that the numerical noises are much smaller than and thus
negligible compared with the “true” physical solution of the 2D turbulent RBC
in an interval of time that is long enough for statistics. In this way, a
convergent (reproducible) solution of the 2D turbulent RBC can be obtained,
which is used here as the “clean” benchmark solution. On the other hand, with
the same initial/boundary conditions and the same physical parameters as
described in the previous section, the NS equations (6)-(7) are also
numerically solved by a traditional algorithm, i.e. the fourth-order Runge-
Kutta’s method with double precision (marked by RKwD) using the time-step
$\Delta t=1\times 10^{-4}$, whose numerical noises increase exponentially
until to the same level of the “true” physical solution and thus are not
negligible. By comparing these RKwD simulations with the CNS benchmark
solution, we can investigate the influence of the numerical noises as tiny
artificial stochastic disturbances on the 2D turbulent RBC in detail. In this
section, we only briefly show some results of comparison. For details about
the CNS algorithms, please refer to Appendix A.

Briefly speaking, to decrease the spatial truncation-error to a small enough
level, we discretize the spatial domain of flow field by a uniform mesh
$N_{x}\times N_{z}=1024\times 1024$, and besides apply the Fourier spectral
method with the $3/2$ rule for dealiasing (Pope, 2001). The corresponding
spatial resolution is high enough for the considered turbulent RBC: the
horizontal (maximum) grid spacing $\Delta_{x}=L/N_{x}=0.00276$ is less than
the minimum Kolmogorov scale (Pope, 2001), which will be shown later in
details. Besides, to decrease the temporal truncation-error to a small enough
level for the CNS, we use the 45th-order (i.e. $M=45$) Taylor expansion with a
time step $\Delta t=10^{-3}$. In addition, to decrease the round-off error to
a small enough level, we use 70 significant digits (i.e. $N_{s}=70$) in
multiple precision for all physical/numerical variables and parameters.
Similarly, we get another CNS result using the Fourier spectral method on the
same uniform mesh $N_{x}\times N_{z}=1024\times 1024$ with the even smaller
background numerical noises by means of a higher order (i.e. $M=47$) Taylor
expansion with the same time step ($\Delta t=10^{-3}$) and the higher multiple
precision with more significant digits (i.e. $N_{s}=72$). Comparing these two
CNS results, it is found that they have no distinct differences in an interval
of time $0\leq t\leq 500$, which is long enough for statistics. This verifies
the convergence and reliability of our CNS result in $t\in[0,500]$ given by
means of $M=45$, $\Delta t=10^{-3}$, and $N_{s}=70$, which is therefore used
in the following parts as the “clean” benchmark solution.

As shown in Figures 2 and 3, the numerical simulation given by the RKwD is
compared with the “clean” benchmark solution given by the CNS, and the
corresponding movie is available at https://github.com/sjtu-
liao/RBC/blob/main/RBC-mv.mp4. Note that these two simulations have exactly
the same initial conditions caused by the micro-level thermal fluctuations.
For both of the CNS and RKwD simulations, the tiny initial disturbances of
velocity and temperature evolve progressively from micro-level to macro-level
until $t\approx 25$ when the transition from laminar flow to turbulence
occurs, then the strong mixing occurs in $t\in[25,36]$ and the typical
vortical/roll-like convection appears at $t\approx 50$, as shown in Figure 2.
Thereafter, as the time increases, the RKwD simulation deviates from the CNS
benchmark solution larger and larger, so that a distinct difference in large-
scale between them can be observed, indicating that the numerical noises (as
artificial stochastic disturbances) could indeed lead to some large-scale
differences of flow fields of velocity and temperature in a macroscopic level,
for example, as shown in Figure 2 at $t=100$ and $t=185$. However, even so,
the flow type of these two simulations keeps the same (i.e. vortical/roll-like
turbulent convection) until $t\approx 188$ when the shearing convection occurs
for the RKwD simulation and its corresponding flow field turns to a kind of
zonal flow thereafter, as shown in Figures 2 and 3. On the contrary, the CNS
benchmark solution thereafter always sustains the non-shearing vortical/roll-
like convection during the whole process of simulation. Therefore, the RKwD
simulation and the CNS benchmark solution have different types of turbulent
convections after $t>188$. It should be emphasized that such kind of
qualitative difference in large-scale is triggered only by the numerical
noises (as artificial stochastic disturbances). All of these highly suggest
that the numerical noises (as artificial stochastic disturbances) have
quantitatively and qualitatively large-scale influences on the sustained
turbulence, i.e. the two-dimensional turbulent RB convection considered in
this paper.

---
Figure 2: Evolutions of the $\theta$ (dimensionless temperature departure from
a linear variation background) field in the case of $Pr=6.8$, $Ra=6.8\times
10^{8}$, and $L/H=2\sqrt{2}$. Left: the CNS benchmark solution; Right: the
RKwD simulation given by the Runge-Kutta’s method with double precision using
$\Delta t=10^{-4}$. The corresponding movie is available at
https://github.com/sjtu-liao/RBC/blob/main/RBC-mv.mp4.
---
Figure 3: Evolutions of the $\theta$ (dimensionless temperature departure from
a linear variation background) field in the case of $Pr=6.8$, $Ra=6.8\times
10^{8}$, and $L/H=2\sqrt{2}$. Left: the CNS benchmark solution; Right: the
RKwD simulation given by the Runge-Kutta’s method with double precision using
$\Delta t=10^{-4}$. The corresponding movie is available at
https://github.com/sjtu-liao/RBC/blob/main/RBC-mv.mp4.

How about the influence of numerical noises (as artificial stochastic
disturbances) on statistical results? The heat transport of 2D turbulent RBC
can be typically quantified by the Nusselt number defined by

$Nu(t)=1-\left.\frac{\partial\langle\theta(x,z,t)\rangle_{x}}{\partial
z}\right|_{z=1},$ (12)

where $\langle a\rangle_{x}=\int_{0}^{\Gamma}adx/\Gamma$ denotes the spatial
average in horizontal direction. As shown in Figure 4(a), the distinct
deviation between the two time histories of $Nu(t)$ given respectively by the
CNS benchmark solution and the RKwD simulation happens at $t\approx 80$ when
the numerical noises of the RKwD simulation have been enlarged to a
macroscopic level because of the butterfly-effect of chaos. Note that the
$Nu(t)$ given by the RKwD simulation drops down greatly at $t\approx 188$
until $Nu(t)\approx$ 30 after $t>300$, one order of magnitude less than that
of the CNS benchmark solution. Such a huge difference is not only quantitative
but also qualitative, which is definitely due to the appearance of the zonal
flow at $t\approx 188$ that is triggered by the numerical noises of the RKwD
simulation (as artificial stochastic disturbances). On the other hand, the
Reynolds number is also calculated to measure the global convection strength,
which is obtained via the root-mean-square (rms) velocity $U_{rms}$ (Sugiyama
et al., 2009; Zhang et al., 2017), i.e.

$Re(t)=\sqrt{\frac{Ra}{Pr}}\hskip 0.85355ptU_{rms}$ (13)

with $U_{rms}=\sqrt{\langle u^{2}+w^{2}\rangle_{A}}$, where $\langle
a\rangle_{A}=\int_{0}^{\Gamma}\int_{0}^{1}adxdz/\Gamma$ denotes the spatial
average. As shown in Figure 4(b), the CNS benchmark solution and RKwD
simulation give almost the same values of $Re(t)$ when t ¡ 188. However, the
departure begins at $t\approx 188$ when the shearing convection occurs, and
thereafter the deviation between the RKwD simulation and the CNS benchmark
solution becomes more and more obvious. Here, it should be emphasized that the
Reynolds number $Re$ of the CNS benchmark solution at $t=500$ is about 3500
larger than that of the RKwD simulation! This is indeed a huge difference.
Besides, the similar phenomena are observed in the comparisons of the
spatially averaged heat flux $\langle wT\rangle_{A}$ and kinetic energy
$\langle E_{V}\rangle_{A}$, given respectively by the CNS benchmark solution
and RKwD simulation, as shown in Figure 5(a) and (b), where the maximum ratio
of kinetic energy reaches about 2.5. Indeed, the numerical noises (as
artificial stochastic disturbances) might lead to qualitatively huge
deviations of the 2D turbulent RBC.

---

Figure 4: Comparisons of the instantaneous Nusselt number $Nu$ and Reynolds
number $Re$ in the case of $Pr=6.8$, $Ra=6.8\times 10^{8}$, and
$L/H=2\sqrt{2}$. (a) The Nusselt number $Nu$; (b) The Reynolds number $Re$.
Solid line in red: the CNS benchmark solution; Dashed line in black: the RKwD
simulation given by the Runge-Kutta’s method with double precision using
$\Delta t=10^{-4}$.

---

Figure 5:  Comparisons of the spatially averaged heat flux $\langle
wT\rangle_{A}$ and kinetic energy $\langle E_{V}\rangle_{A}$ in the case of
$Pr=6.8$, $Ra=6.8\times 10^{8}$, and $L/H=2\sqrt{2}$, where $\langle
a\rangle_{A}=\int_{0}^{\Gamma}\int_{0}^{1}adxdz/\Gamma$ denotes the spatial
average. (a) The heat flux $\langle wT\rangle_{A}$; (b) The kinetic energy
$\langle E_{V}\rangle_{A}$. Solid line in red: the CNS benchmark solution;
Dashed line in black: the RKwD simulation given by the Runge-Kutta’s method
with double precision using $\Delta t=10^{-4}$.

---

Figure 6: Comparisons of the spatially averaged kinetic energy dissipation
rate $\langle\varepsilon_{V}\rangle_{A}$ and thermal dissipation rate
$\langle\varepsilon_{T}\rangle_{A}$ in the case of $Pr=6.8$, $Ra=6.8\times
10^{8}$, and $L/H=2\sqrt{2}$, where $\langle
a\rangle_{A}=\int_{0}^{\Gamma}\int_{0}^{1}adxdz/\Gamma$ denotes the spatial
average. (a) The kinetic energy dissipation rate
$\langle\varepsilon_{V}\rangle_{A}$; (b) The thermal dissipation rate
$\langle\varepsilon_{T}\rangle_{A}$. Solid line in red: the CNS benchmark
solution; Dashed line in black: the RKwD simulation given by the Runge-Kutta’s
method with double precision using $\Delta t=10^{-4}$.

Besides the above-mentioned large-scale quantities, the comparisons of some
small-scale properties of fluid flows such as the kinetic energy dissipation
rate $\langle\varepsilon_{V}\rangle_{A}$ and the thermal dissipation rate
$\langle\varepsilon_{T}\rangle_{A}$ of this 2D turbulent RBC given by the CNS
solution and the RKwD simulation are as shown in Figure 6, where

$\varepsilon_{V}(x,z,t)=\frac{1}{2}\sqrt{\frac{Pr}{Ra}}\sum_{ij}\left[\partial_{i}u_{j}(x,z,t)+\partial_{j}u_{i}(x,z,t)\right]^{2}$
(14)

and

$\varepsilon_{T}(x,z,t)=\frac{1}{\sqrt{Pr\,Ra}}\hskip
0.85355pt\Big{|}\nabla\left[\theta(x,z,t)-z\right]\Big{|}^{2}$ (15)

with $i,\,j=1,\,2$, $u_{1}(x,z,t)=u(x,z,t)$, $u_{2}(x,z,t)=w(x,z,t)$,
$\partial_{1}=\partial/\partial x$, $\partial_{2}=\partial/\partial z$, and
that $\nabla$ is the Hamilton operator and $\langle\;\;\rangle_{A}$ denotes
the spatial average, respectively. Note that, for the RKwD simulation, both of
the kinetic energy dissipation rate $\langle\varepsilon_{V}\rangle_{A}$ and
thermal dissipation rate $\langle\varepsilon_{T}\rangle_{A}$ greatly drop down
at $t\approx 188$ when the shearing convection (i.e. the zonal flow) occurs,
which is triggered only by the numerical noises (as artificial stochastic
disturbances).

By the way, as shown in Figure 6(a), we have the maximum kinetic energy
dissipation rate $(\langle\varepsilon_{V}\rangle_{A})_{max}=0.00646$ at
$t=26.9$ when the transition from the laminar flow to turbulence occurs,
corresponding to the minimum Kolmogorov scale

$\left(\langle\eta\rangle_{A}\right)_{min}\approx(Pr/Ra)^{3/8}\big{[}(\langle\varepsilon_{V}\rangle_{A})_{max}\big{]}^{-1/4}=0.00353.$

Thus, the criteria condition that the maximum horizontal grid spacing
$\Delta_{x}=\Gamma/N_{x}=0.00276<0.8(\langle\eta\rangle_{A})_{min}=0.00282$ is
indeed satisfied, say, the spatial resolution used in this paper is fine
enough for the 2D turbulent RBC under consideration.

Figure 7 shows the comparisons of the kinetic energy dissipation rate
$\langle\varepsilon_{V}\rangle_{x,t}(z)$ and the heat flux $\langle
wT\rangle_{x,t}(z)$, where $\langle
a\rangle_{x,t}=\int_{0}^{\Gamma}\int_{0}^{500}a\;dxdt/\Gamma/500$ denotes the
horizontally spatial and temporal average. Obviously, both of the kinetic
energy dissipation rate $\langle\varepsilon_{V}\rangle_{x,t}$ and the heat
flux $\langle wT\rangle_{x,t}$ of the CNS benchmark solution are significantly
larger than those given by the RKwD simulation. Especially, near the lower and
upper plates, the $\langle\varepsilon_{V}\rangle_{x,t}$ given by the CNS
benchmark solution has a much sharper peak than that given by the RKwD
simulation, as shown in Figure 7(a). These indicate that the numerical noises
(as artificial stochastic disturbances) can lead to large-scale deviations in
statistics of the 2D turbulent RBC under consideration.

---

Figure 7: Comparisons of the horizontally and temporally averaged kinetic
energy dissipation rate $\langle\varepsilon_{V}\rangle_{x,t}(z)$ and the heat
flux $\langle wT\rangle_{x,t}(z)$ in the case of $Pr=6.8$, $Ra=6.8\times
10^{8}$, and $L/H=2\sqrt{2}$, where $\langle
a\rangle_{x,t}=\int_{0}^{\Gamma}\int_{0}^{500}adxdt/\Gamma/500$ denotes the
horizontal and temporal average. (a) The kinetic energy dissipation rate
$\langle\varepsilon_{V}\rangle_{x,t}(z)$; (b) The heat flux $\langle
wT\rangle_{x,t}(z)$. Solid line in red: the CNS benchmark solution; Dashed
line in black: the RKwD simulation given by the Runge-Kutta’s method with
double precision using $\Delta t=10^{-4}$.

The comparison of the probability density functions (PDFs) of the stream
function $\psi(x,z,t)$ in $0\leq x<\Gamma$, $0\leq z\leq 1$ and $0\leq t\leq
500$ given by the CNS benchmark solution and the RKwD simulation is as shown
in Figure 9. Unlike the PDF of the CNS benchmark solution that has a kind of
asymmetry about $\psi=0$, the PDF of the RKwD simulation has no such kind of
asymmetry but two peaks at $\psi\approx 0$ and $\psi\approx 0.25$.
Furthermore, the comparison of the PDFs of $\theta(x,z,t)$ given by the CNS
benchmark solution and RKwD simulation is as shown in Figure 9. As shown in
Figure 9(a), except at $\theta\approx 0$, the PDF of $\theta(x,z,t)$ given by
the CNS benchmark solution remains almost the same value. In contrast, the PDF
of $\theta$ given by the RKwD simulation is relatively more typical, as shown
in Figure 9(b). Thus, the numerical noises (as artificial stochastic
disturbances) indeed lead to large-scale deviations even in the PDFs of the 2D
turbulent RBC under consideration.

---

Figure 8: Probability density functions (PDFs) of $\psi(x,z,t)$ in $0\leq
x<\Gamma$, $0\leq z\leq 1$ and $0\leq t\leq 500$ in the case of $Pr=6.8$,
$Ra=6.8\times 10^{8}$, and $L/H=2\sqrt{2}$. (a) The PDF given by the CNS
benchmark solution; (b) The PDF given by the RKwD simulation obtained via the
Runge-Kutta’s method with double precision using $\Delta t=10^{-4}$.

---

Figure 9: Probability density functions (PDFs) of $\theta(x,z,t)$ in $0\leq
x<\Gamma$, $0\leq z\leq 1$ and $0\leq t\leq 500$ in the case of $Pr=6.8$,
$Ra=6.8\times 10^{8}$, and $L/H=2\sqrt{2}$. (a) The PDF given by the CNS
benchmark solution; (b) The PDF given by the RKwD simulation obtained via the
Runge-Kutta’s method with double precision using $\Delta t=10^{-4}$.

All of the above-mentioned comparisons indicate that the micro-level
background numerical noises as a kind of artificial stochastic disturbances
can lead to large-scale differences not only in spatio-temporal trajectories
but also even in flow types of the 2D turbulent RBC, which further affects the
statistics of the Nusselt number, the Reynolds number, the kinetic energy, the
kinetic energy dissipation rate, the thermal dissipation rate and so on. Note
that it is currently reported by McMullen et al. (2022) that “the Navier-
Stokes equations do not describe turbulent gas flows in the dissipation range
because they neglect thermal fluctuations”, say, tiny stochastic disturbances
resulting from thermal fluctuations might influence the small-scale properties
of the freely decaying turbulent flows under their consideration. In this
paper, the detailed comparisons between the CNS benchmark solution and the
RKwD simulation provide us the rigorous evidence that numerical noises as a
kind of small-scale artificial stochastic disturbances might influence the
large-scale properties of a sustained turbulence, i.e. the 2D turbulent RB
convection considered in this paper.

##  4 Concluding remarks and discussions

All numerical algorithms have the background numerical noises, i.e. the
truncation errors and round-off errors, which are tiny and random. It was
reported that, for a chaotic dynamic system, the random numerical noises
increase exponentially due to the butterfly-effect of chaos, until up to the
same order of magnitude as its “true” physical solution (Hu & Liao, 2020; Qin
& Liao, 2020). Therefore, numerical simulations of a deterministic chaotic
system are mostly a mixture of the “true” physical solution, which is
deterministic in physics, and the “false” numerical noises, which is however
stochastic. This is the reason why numerical simulations of a deterministic
chaotic system given by traditional algorithms in single/double precision
often look like stochastic, and why it has been wrongly believed that a
deterministic chaotic system can lead to randomness. This is also the reason
why Teixeira et al. (2007) made a rather pessimistic conclusion that “for
chaotic systems, numerical convergence cannot be guaranteed forever”.
Obviously, even given an accurate initial condition, this kind of randomness
of chaotic system comes from the randomness of artificial background numerical
noises. Thus, the background numerical noises can be naturally regarded as a
kind of artificial stochastic disturbances to a chaotic system when it is
solved numerically.

In this paper, the so-called clean numerical simulation (CNS) is adopted to
accurately investigate the influence of numerical noises as a kind of tiny
artificial stochastic disturbances on the two-dimensional turbulent RB
convection under consideration. This is mainly because the CNS can reduce the
background numerical noises (i.e. round-off errors and truncation errors) to
any a required level, which can be so small that the “false” numerical noises
are negligible compared with its “true” physical solution and thus the
numerical simulations of turbulence are convergent/reproducible in an interval
of time long enough for statistics, as illustrated in this paper. Strictly
speaking, the CNS solution is also a mixture of the “true” physical solution
and the “false” numerical noises. However, unlike the simulation given by a
traditional algorithm in single/double precision, the “false” numerical noises
of a CNS solution are often several orders of magnitude smaller than its
“true” physical solution in a region $t\in[0,T_{c}]$ long enough for
statistics, so that its numerical noises are negligible and the CNS result is
“convergent” and “reproducible”. Such kind of CNS solution in $t\in[0,T_{c}]$
can be used as a “clean” benchmark solution for comparison with those given by
traditional algorithms in single/double precision so as to investigate the
influence of numerical noises as tiny artificial stochastic disturbances on
the two-dimensional turbulent RB convection under consideration. So, unlike
the Taylor series method, the key point of the CNS is the so-called “critical
predictable time” $T_{c}$ that determines a temporal interval $[0,T_{c}]$ in
which the numerical simulations are “reliable” and “clean”, since their
“false” numerical noises are much smaller than the “true” physical solution
and thus are negligible.

It was currently reported by McMullen et al. (2022) that tiny stochastic
disturbances resulting from thermal fluctuations might influence the small-
scale properties of the freely decaying turbulence under their consideration,
which is in agreement with the conclusions given by Gallis et al. (2021),
Bandak et al. (2022), Bell et al. (2022), Eyink & Jafari (2022) and so on. In
this paper, we investigate the large-scale influence of numerical noises as a
kind of tiny artificial stochastic disturbances on a sustained turbulence.
Using the two-dimensional (2D) turbulent Rayleigh-Bénard convection (RBC) as
an example, we illustrate that the numerical noises as a kind of micro-level
artificial stochastic disturbances could indeed lead to large-scale deviations
not only in spatio-temporal trajectories but also even in statistics of the
sustained turbulence considered in this paper. Especially, such kind of tiny
artificial stochastic disturbances even leads to different types of flows: the
shearing convection occurs for the RKwD simulations, and its corresponding
flow field turns to a kind of zonal flow thereafter, however the CNS benchmark
solution always sustains the non-shearing vortical/roll-like convection during
the whole process of simulation, as shown in Figures 2 and 3. Thus, we provide
a rigorous evidence that numerical noises as a kind of tiny artificial
stochastic disturbances have not only quantitatively but also qualitatively
large-scale influences on a sustained turbulence, i.e. the 2D turbulent RB
convection considered in this paper. Of course, for various types of turbulent
flows governed by the NS equations, more investigations are needed in the
future.

---

Figure 10: Comparison between the CNS benchmark solution $\theta_{CNS}$ and
the numerical noises $\varepsilon_{n}=\theta_{CNS}-\theta_{RKwD}$ at the point
$x=1$ and $y=1/2$, where $\theta_{RKwD}$ is the RKwD simulation. Solid line:
the CNS benchmark solution; Dashed line: the numerical noises. (a) $1\leq
t\leq 500$ ; (b) $75\leq t\leq 100$.

---

Figure 11: Temporal evolution of the Nusselt number $Nu$ and the Reynolds
number $Re$ in the case of $Pr=6.8$, $Ra=6.8\times 10^{8}$, and
$L/H=2\sqrt{2}$. (a) Nusselt number $Nu$; (b) Reynolds number $Re$; Solid line
in red: the CNS benchmark solution; Dashed line in black: the RKwD simulation
(with $\Delta t=10^{-4}$); Dashed line in blue: the RKwD′ simulation (with
$\Delta t=2\times 10^{-4}$).

Why? We try to give an explanation. This 2D turbulent RBC system might have
two possible final states: the vortical/roll-like flow and the zonal flow,
which can be seen as two minima of a double-well potential. Once the CNS
benchmark solution falls in one of these two minima (e.g. the roll-like flow),
it remains there forever, since the corresponding “false” numerical noises are
much less than the “true” physical solution and thus can not trigger off the
transition to another minimum. On the contrary, the RKwD simulation can
perform the transition from one state to the other, because its “false”
numerical noises as artificial stochastic disturbances might be at the same
order of magnitude as the “true” physical solution, as shown in Figure 10, so
that the RKwD simulation might depart very far away from its “true” physical
solution and thus fall in another minimum. If so, the transition of the RKwD
simulation from the vortical/roll-like flow to the zonal flow of the 2D
turbulent RBC should occur at random times for different numerical noises as a
kind of artificial stochastic disturbances. In order to check this hypothesis,
we perform a new RKwD simulation of the 2D turbulent RBC with the identical
initial/boundary conditions and physical parameters of the previous one, but
using a different time-step $\Delta t=2\times 10^{-4}$, marked by RKwD′. The
corresponding numerical noises as a different kind of artificial stochastic
disturbances are verified to be at the same tiny level as the previous ones.
As shown in Figure 11 (that describes the temporal evolutions of the Nusselt
number $N{u}$ and the Reynolds number $Re$), the transition from the
vortical/roll-like flow to the zonal flow of the RKwD′ simulation with $\Delta
t=2\times 10^{-4}$ occurs indeed at a different time $t\approx 230$, compared
to $t\approx 188$ of the previous RKwD simulation with $\Delta t=10^{-4}$.
Thus, our above-mentioned explanation in terms of a stochastic dynamical
system in the double-well potential should be reasonable. This test provides
useful information to better understand the origin of the phenomenon reported
in this paper.

Thus, generally speaking, if a chaotic system has $N$-well potential, where
$N\geq 2$, its numerical simulations given by traditional algorithms in
single/double precision should fall randomly in one of the $N$ minimums, and
besides the transition between different minimums should occur frequently and
randomly. Similarly, if a turbulent flow has multiple states in a large scale,
small disturbances might lead to the transition between different states. Such
kind of small disturbances can be either natural (such as the environmental
perturbations) or artificial (such as the background numerical noises). In
this case, it is impossible to make a correct prediction of flow states by
means of traditional algorithms in single/double precision. By contrast, if a
turbulent flow has an unique state, the background numerical noises should not
have a large-scale influence on the flow. This should be a piece of good news
for researchers in the field of CFD.

In theory, even if a turbulent flow has multiple states in large-scale, we can
still give a deterministic prediction of its large-scale states by means of
the CNS, since the artificial background numerical noises are negligible in a
long enough interval of time for CNS solution. However, a natural disturbance,
say, thermal fluctuation, always exists in fluid flows, which is random,
small-scale, but unavoidable in practice. It is an open question whether or
not, like tiny artificial numerical noises, the micro-scale thermal
fluctuation might have a large-scale influence on some turbulent flows with
multiple states. Note that the influence of thermal fluctuation is not
considered in the NS equations. So, it is highly suggested to investigate the
large-scale influence of thermal fluctuation on turbulent flows with multiple
states by means of the Landau-Lifshitz-Navier-Stokes (LLNS) equations (Landau
& Lifshitz, 1959). It is also worthwhile seriously discussing which one is
better for turbulent flows (especially those with large-scale multiple
states), either the deterministic NS equations or the stochastic LLNS
equations. Obviously, more investigations are necessary in the future.

It should be emphasized that the external noises considered in most of related
articles are mostly about ten orders of magnitude larger than the numerical
noises considered in this paper, which might greatly change the
characteristics of the chaotic/turbulent systems under their consideration.
Besides, even when their external noises are zero, their numerical simulations
are in fact a mixture of the “true” physical solution and the “false”
numerical noises: both of them might be at the same level, as shown in this
paper. Thus, strictly speaking, the conclusions based on this kind of
“mixtures” should be not rigorous in theory. Thus, the CNS provides us, for
the first time, a rigorous way to investigate the influence of external
disturbances and very tiny artificial numerical noises on chaotic dynamical
systems and turbulent flows.

In summary, numerical noises as weak, small-scale stochastic perturbations
increase exponentially to a macro level of numerical simulations, and besides
might have a large influence on the macroscopic statistics of turbulent flows.
Therefore, we should pay more attentions to the influences of small-scale
stochastic perturbations on turbulence. Finally, this work also illustrates
the validity, novelty, and great potential of the CNS as a reliable and
accurate tool in theoretical studies of turbulence. Hopefully, the CNS might
provide a brand-new, extremely accurate numerical tool to study turbulent
flows.

[Acknowledgements]Thanks to the anonymous reviewers for their valuable
suggestions and constructive comments.

[Funding]This work is partly supported by the National Natural Science
Foundation of China (No. 91752104) and Shanghai Pilot Program for Basic
Research - Shanghai Jiao Tong University (No. 21TQ1400202).

[Declaration of Interests]The authors report no conflict of interest.

[Data availability statement]The data that support the findings of this study
are available on request from the corresponding author. The movie related to
Figures 2 and 3 is available at https://github.com/sjtu-
liao/RBC/blob/main/RBC-mv.mp4

[Author ORCID]Shijie Qin, https://orcid.org/0000-0002-0809-1766; Shijun Liao,
https://orcid.org/0000-0002-2372-9502

[Author contributions]Liao conceived and designed the analysis. Qin performed
the analysis. Both wrote the manuscript.

##  Appendix A The CNS algorithm for 2D turbulent RBC

Lin et al. (2017) combined the clean numerical simulation (CNS) with the
traditional Fourier-Galerkin spectral method (in spectral space) to solve a
two-dimensional (2D) turbulent Rayleigh-Bénard convection (RBC) with free-slip
boundary conditions. However, their approach is rather time-consuming.
Currently, Hu & Liao (2020) and Qin & Liao (2020) proposed an efficient CNS
algorithm in physical space for spatio-temporal chaos to overcome the
shortcoming of the CNS algorithm in spectral space. Here, the basic idea of
this kind of efficient CNS algorithm in physical space is briefly described by
using the 2D turbulent RBC as an example.

###  A.1 The CNS algorithm in physical space

Applying the coordinate transformations $\tilde{x}=\lambda x$ and
$\tilde{z}=\mu z$ to the NS equations (6) and (7), where $\lambda=2\pi/\Gamma$
and $\mu=\pi$, we obtain the following governing equations:

$\displaystyle\frac{\partial}{\partial
t}(\lambda^{2}\psi_{xx}+\mu^{2}\psi_{zz})=\lambda\mu\hskip
0.85358pt\psi_{z}\left(\lambda^{2}\psi_{xxx}+\mu^{2}\psi_{xzz}\right)-\lambda\mu\hskip
0.85358pt\psi_{x}\left(\lambda^{2}\psi_{xxz}+\mu^{2}\psi_{zzz}\right)$
$\displaystyle+\lambda\hskip
0.85358pt\theta_{x}+\sqrt{\frac{Pr}{Ra}}\left(\lambda^{4}\psi_{xxxx}+2\lambda^{2}\mu^{2}\psi_{xxzz}+\mu^{4}\psi_{zzzz}\right),$
(16) $\frac{\partial\theta}{\partial t}=\lambda\mu\hskip
0.85358pt(\psi_{z}\theta_{x}-\psi_{x}\theta_{z})+\lambda\hskip
0.85358pt\psi_{x}+\frac{1}{\sqrt{PrRa}}(\lambda^{2}\theta_{xx}+\mu^{2}\theta_{zz}),$
(17)

with $t\geq 0$, $x\in[0,2\pi]$ and $z\in[0,\pi]$, where $x$ and $z$ as
subscripts denote the spatial derivatives and the overhead tildes are omitted
in the remainder of this article.

We prolongate the computational domain from $z\in[0,\pi]$ to $z\in[0,2\pi]$ so
as to easily satisfy the free-slip boundary conditions at lower ($z=0$) and
upper ($z=\pi$) plates by means of the Fourier series. We use $N_{x}\times
N_{z}$ equidistant points, i.e.

$x_{j}=\frac{2\pi}{N_{x}}j,\hskip 28.45274ptz_{k}=\frac{2\pi}{N_{z}}k,$ (18)

where $j=0,\,1,\,2,\,...,\,N_{x}-1$ and $k=0,\,1,\,2,\,...,\,N_{z}-1$, to
discretize $\psi$ and $\theta$, respectively.

To reduce truncation errors in the temporal dimension, the high-order Taylor
expansions are adopted, i.e.

$\psi(x_{j},z_{k},t+\Delta
t)\approx\sum^{M}_{m=0}\psi^{[m]}(x_{j},z_{k},t)(\Delta t)^{m},$ (19)
$\theta(x_{j},z_{k},t+\Delta
t)\approx\sum^{M}_{m=0}\theta^{[m]}(x_{j},z_{k},t)(\Delta t)^{m},$ (20)

where $\Delta t$ is the time step, $M$ is the order of Taylor expansion, with
the definitions

$\psi^{[m]}(x_{j},z_{k},t)=\frac{1}{m!}\frac{\partial^{m}\psi(x_{j},z_{k},t)}{\partial
t^{m}},\hskip
28.45274pt\theta^{[m]}(x_{j},z_{k},t)=\frac{1}{m!}\frac{\partial^{m}\theta(x_{j},z_{k},t)}{\partial
t^{m}}.$ (21)

Here, the order $M$ should be large enough so as to reduce the truncation
errors (in the temporal dimension) to a required tiny level.

Differentiating $(m-1)$ times both sides of (16) and (17) with respect to $t$
and then dividing them by $m!$, we obtain the governing equations of
$\psi^{[m]}$ and $\theta^{[m]}$:

$\displaystyle\lambda^{2}\psi^{[m]}_{xx}(x_{j},z_{k},t)+\mu^{2}\psi^{[m]}_{zz}(x_{j},z_{k},t)=\frac{1}{m}\Big{\\{}\sqrt{\frac{Pr}{Ra}}\big{[}2\lambda^{2}\mu^{2}\psi^{[m-1]}_{xxzz}(x_{j},z_{k},t)$
$\displaystyle+\lambda^{4}\psi^{[m-1]}_{xxxx}(x_{j},z_{k},t)+\mu^{4}\psi^{[m-1]}_{zzzz}(x_{j},z_{k},t)\big{]}$
$\displaystyle+\sum^{m-1}_{r=0}\lambda\mu\hskip
0.85358pt\psi^{[r]}_{z}(x_{j},z_{k},t)\big{[}\lambda^{2}\psi^{[m-1-r]}_{xxx}(x_{j},z_{k},t)+\mu^{2}\psi^{[m-1-r]}_{xzz}(x_{j},z_{k},t)\big{]}$
$\displaystyle-\sum^{m-1}_{r=0}\lambda\mu\hskip
0.85358pt\psi^{[r]}_{x}(x_{j},z_{k},t)\big{[}\lambda^{2}\psi^{[m-1-r]}_{xxz}(x_{j},z_{k},t)+\mu^{2}\psi^{[m-1-r]}_{zzz}(x_{j},z_{k},t)\big{]}$
$\displaystyle+\lambda\hskip
0.85358pt\theta^{[m-1]}_{x}(x_{j},z_{k},t)\Big{\\}},$ (22)
$\displaystyle\theta^{[m]}(x_{j},z_{k},t)=\frac{1}{m}\Big{\\{}\frac{1}{\sqrt{PrRa}}\big{[}\lambda^{2}\theta^{[m-1]}_{xx}(x_{j},z_{k},t)+\mu^{2}\theta^{[m-1]}_{zz}(x_{j},z_{k},t)\big{]}$
$\displaystyle+\lambda\mu\sum^{m-1}_{r=0}\psi^{[r]}_{z}(x_{j},z_{k},t)\theta^{[m-1-r]}_{x}(x_{j},z_{k},t)$
$\displaystyle-\lambda\mu\sum^{m-1}_{r=0}\psi^{[r]}_{x}(x_{j},z_{k},t)\theta^{[m-1-r]}_{z}(x_{j},z_{k},t)$
$\displaystyle+\lambda\hskip
0.85358pt\psi^{[m-1]}_{x}(x_{j},z_{k},t)\Big{\\}},$ (23)

where $m\geq 1$.

Note that there exist some spatial partial derivatives (denoted by subscripts)
in (22) and (23), such as $\partial^{s_{1}+s_{2}}\psi^{[r]}/(\partial
x^{s_{1}}\partial z^{s_{2}})$ and
$\partial^{s_{1}+s_{2}}\theta^{[r]}/(\partial x^{s_{1}}\partial z^{s_{2}})$
with $r,\,s_{1},\,s_{2}\geq 0$. In order to approximate these spatial partial
derivative terms with high computational efficiency and precision from the
known discrete variables $\psi^{[r]}(x_{j},z_{k},t)$ and
$\theta^{[r]}(x_{j},z_{k},t)$, we adopt the spatial Fourier series

$\psi^{[r]}(x,z,t)\approx\sum^{\frac{N_{x}}{2}-1}_{n_{x}=-\frac{N_{x}}{2}+1}\,\sum^{\frac{N_{z}}{2}-1}_{n_{z}=-\frac{N_{z}}{2}+1}\Psi^{[r]}(n_{x},n_{z},t)e^{i\hskip
0.85355ptn_{x}x}e^{i\hskip 0.85355ptn_{z}z},$ (24)
$\theta^{[r]}(x,z,t)\approx\sum^{\frac{N_{x}}{2}-1}_{n_{x}=-\frac{N_{x}}{2}+1}\,\sum^{\frac{N_{z}}{2}-1}_{n_{z}=-\frac{N_{z}}{2}+1}\Theta^{[r]}(n_{x},n_{z},t)e^{i\hskip
0.85355ptn_{x}x}e^{i\hskip 0.85355ptn_{z}z},$ (25)

where $i=\sqrt{-1}$ is the imaginary unit and

$\Psi^{[r]}(n_{x},n_{z},t)=\frac{1}{N_{x}\hskip
0.85355ptN_{z}}\sum^{N_{x}-1}_{j=0}\sum^{N_{z}-1}_{k=0}\psi^{[r]}(x_{j},z_{k},t)e^{-i\hskip
0.85355ptn_{x}x_{j}}e^{-i\hskip 0.85355ptn_{z}z_{k}},$ (26)
$\Theta^{[r]}(n_{x},n_{z},t)=\frac{1}{N_{x}\hskip
0.85355ptN_{z}}\sum^{N_{x}-1}_{j=0}\sum^{N_{z}-1}_{k=0}\theta^{[r]}(x_{j},z_{k},t)e^{-i\hskip
0.85355ptn_{x}x_{j}}e^{-i\hskip 0.85355ptn_{z}z_{k}},$ (27)

are determined by the known $\psi^{[r]}(x_{j},z_{k},t)$ and
$\theta^{[r]}(x_{j},z_{k},t)$, respectively, at discrete points
($x_{j},z_{k}$) with $j=0,\,...,\,N_{x}-1$ and $k=0,\,...,\,N_{z}-1$. Then, we
can obtain the rather accurate approximations of the spatial partial
derivative terms in (22) and (23), say,

$\displaystyle\frac{\partial^{s_{1}+s_{2}}\psi^{[r]}}{\partial
x^{s_{1}}\partial z^{s_{2}}}(x_{j},z_{k},t)$ (28) $\displaystyle\approx$
$\displaystyle
i^{s_{1}+s_{2}}\sum^{\frac{N_{x}}{2}-1}_{n_{x}=-\frac{N_{x}}{2}+1}\,\sum^{\frac{N_{z}}{2}-1}_{n_{z}=-\frac{N_{z}}{2}+1}(n_{x})^{s_{1}}(n_{z})^{s_{2}}\hskip
0.85355pt\Psi^{[r]}(n_{x},n_{z},t)e^{i\hskip 0.85355ptn_{x}x_{j}}e^{i\hskip
0.85355ptn_{z}z_{k}},$
$\displaystyle\frac{\partial^{s_{1}+s_{2}}\theta^{[r]}}{\partial
x^{s_{1}}\partial z^{s_{2}}}(x_{j},z_{k},t)$ (29) $\displaystyle\approx$
$\displaystyle
i^{s_{1}+s_{2}}\sum^{\frac{N_{x}}{2}-1}_{n_{x}=-\frac{N_{x}}{2}+1}\,\sum^{\frac{N_{z}}{2}-1}_{n_{z}=-\frac{N_{z}}{2}+1}(n_{x})^{s_{1}}(n_{z})^{s_{2}}\hskip
0.85355pt\Theta^{[r]}(n_{x},n_{z},t)e^{i\hskip 0.85355ptn_{x}x_{j}}e^{i\hskip
0.85355ptn_{z}z_{k}}.$

Here, the fast Fourier transform (FFT) algorithm is used to increase
computational efficiency. In this way, the spatial truncation error can be
decreased to a required tiny level, as long as the mode numbers $N_{x}$ and
$N_{z}$ are large enough.

Note that, if the order $M$ of the Taylor expansions (19) and (20) is large
enough, the temporal truncation errors can be controlled below a required tiny
level. Besides, if spatial discretizations are fine enough, say, the mode
numbers $N_{x}$ and $N_{z}$ are large enough, the spatial truncation errors in
Fourier expressions (24) and (25) and besides the corresponding spatial
derivative terms in (22) and (23) can be reduced to a required tiny level.
However, this is not enough, since there always exist some round-off errors,
which might be larger than the truncation errors. So, in addition, all
physical/numerical variables and parameters are expressed in multiple
precision (MP) with a large enough number $N_{s}$ of significant digits so
that the round-off errors can also be controlled below a required tiny level.
In this way, the background numerical noises, i.e. both of the
spatial/temporal truncation errors and the round-off error as a whole, can be
controlled below a required tiny level. This is different from other numerical
algorithms including the DNS. In theory, the CNS results are more accurate
than those given by the DNS, since the CNS adopts the multiple precision,
however the DNS mostly uses the double precision. In addition, note that the
CNS results are useful only in a limited interval of time $t\in[0,T_{c}]$, in
which numerical noises can be neglected.

###  A.2 Realization of free-slip boundary conditions

Considering the Fourier expansions (24) and (25), the complex coefficients
$\Psi^{[r]}(n_{x},n_{z},t)$ and $\Theta^{[r]}(n_{x},n_{z},t)$ can be expressed
in terms of their real and imaginary parts, respectively, as follows:

$\Psi^{[r]}(n_{x},n_{z},t)=\Psi_{1}^{[r]}(n_{x},n_{z},t)+i\hskip
0.85355pt\Psi_{2}^{[r]}(n_{x},n_{z},t),$ (30)
$\Theta^{[r]}(n_{x},n_{z},t)=\Theta_{1}^{[r]}(n_{x},n_{z},t)+i\hskip
0.85355pt\Theta_{2}^{[r]}(n_{x},n_{z},t),$ (31)

where $r\geq 0$, $-N_{x}/2+1\leq n_{x}\leq N_{x}/2-1$ and $-N_{z}/2+1\leq
n_{z}\leq N_{z}/2-1$. Using the conjugate symmetry of
$\Psi^{[r]}(n_{x},n_{z},t)$ and $\Theta^{[r]}(n_{x},n_{z},t)$, we enforce

$\Psi_{1}^{[r]}(n_{x},n_{z},t)=\Psi_{1}^{[r]}(-n_{x},-n_{z},t)=-\Psi_{1}^{[r]}(n_{x},-n_{z},t)=-\Psi_{1}^{[r]}(-n_{x},n_{z},t),$
(32)
$\Psi_{2}^{[r]}(n_{x},n_{z},t)=-\Psi_{2}^{[r]}(-n_{x},-n_{z},t)=-\Psi_{2}^{[r]}(n_{x},-n_{z},t)=\Psi_{2}^{[r]}(-n_{x},n_{z},t),$
(33)
$\Theta_{1}^{[r]}(n_{x},n_{z},t)=\Theta_{1}^{[r]}(-n_{x},-n_{z},t)=-\Theta_{1}^{[r]}(n_{x},-n_{z},t)=-\Theta_{1}^{[r]}(-n_{x},n_{z},t),$
(34)
$\Theta_{2}^{[r]}(n_{x},n_{z},t)=-\Theta_{2}^{[r]}(-n_{x},-n_{z},t)=-\Theta_{2}^{[r]}(n_{x},-n_{z},t)=\Theta_{2}^{[r]}(-n_{x},n_{z},t),$
(35)

so as to automatically satisfy the free-slip boundary conditions at $z=0$ and
$z=\pi$, corresponding to the upper and lower boundaries, respectively.
Besides, there should exist

$\Psi_{1}^{[r]}(0,n_{z},t)=\Psi_{1}^{[r]}(n_{x},0,t)=\Psi_{2}^{[r]}(n_{x},0,t)=0,$
(36)
$\Theta_{1}^{[r]}(0,n_{z},t)=\Theta_{1}^{[r]}(n_{x},0,t)=\Theta_{2}^{[r]}(n_{x},0,t)=0.$
(37)

For more details, please refer to Saltzman (1962). Actually, the equations
(32)-(37) imply that $\psi^{[r]}(x,z,t)$ and $\theta^{[r]}(x,z,t)$ in (24) and
(25) are expanded as the sine series in the vertical direction, which
automatically satisfy the free-slip boundary conditions

$\psi^{[r]}(x,z,t)=\frac{\partial^{2}\psi^{[r]}(x,z,t)}{\partial
z^{2}}=\theta^{[r]}(x,z,t)=0$ (38)

at $z=0$ and $z=\pi$, say,

$\psi^{[r]}(x_{j},z_{k},t)=\psi_{zz}^{[r]}(x_{j},z_{k},t)=\theta^{[r]}(x_{j},z_{k},t)=0$
(39)

at $k=0$ and $k=N_{z}/2$.

###  A.3 Computational efficiency

In this section, we illustrate that the above-mentioned CNS algorithm in
physical space needs much less amount of calculation than the previous CNS
algorithm combined with the Fourier-Galerkin spectral method (Lin et al.,
2017). Table 1 illustrates the required CPU times of the CNS algorithm
combined with the Fourier-Galerkin spectral method (marked by $T_{1}$) and the
CNS algorithm in physical space described above (marked by $T_{2}$) by means
of different numbers $np$ of threads in parallel computing. Note that $T_{1}$
is much larger than $T_{2}$ in the case of $np=16$, $32$, $64$ or $128$: the
corresponding time ratio $T_{1}/T_{2}$ is more than $1500$ in all cases. In
consequence, the CNS algorithm in physical space described above is much more
efficient (about three level of magnitude higher) than the previous CNS
algorithm combined with the Fourier-Galerkin spectral method (Lin et al.,
2017), and thus is used in this paper.

$~{}~{}np$ | $T_{1}$ | $T_{2}$ | $T_{1}/T_{2}~{}~{}~{}$
---|---|---|---
(in seconds) | (in seconds)
$~{}~{}16$ | $475348.2$ | $282.0$ | $1686~{}~{}~{}$
$~{}~{}32$ | $237953.8$ | $138.4$ | $1719~{}~{}~{}$
$~{}~{}64$ | $123709.2$ | $77.5$ | $1596~{}~{}~{}$
$~{}~{}128$ | $65149.4$ | $43.2$ | $1508~{}~{}~{}$
Table 1: Required CPU times of the CNS algorithm combined with the Fourier-
Galerkin spectral method (marked by $T_{1}$) and the CNS algorithm in physical
space (marked by $T_{2}$) for solving the governing equations (6) and (7) in
the temporal interval $t\in[0,0.1]$ with the same physical and numerical
parameters: the aspect ratio $\Gamma=L/H=2\sqrt{2}$, Prandtl number $Pr=6.8$,
Rayleigh number $Ra=10^{7}$, the mode numbers $|n_{x}|,|n_{z}|\leq 127$ of
spatial discretizations, the $10$-th order ($M=10$) of the truncated Taylor
series in temporal dimension with a fixed time step $\Delta t=10^{-3}$, and
the number $N_{s}=100$ of significant digits in the multiple precision scheme,
where $np$ denotes the number of threads in parallel computing.

## References

  * Ahlers et al. (2009) Ahlers, G., Grossmann, S. & Lohse, D. 2009 Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection. Rev. Mod. Phys. 81 (2), 503–503.
  * Bandak et al. (2022) Bandak, D., Goldenfeld, N., Mailybaev, A.A. & Eyink, G. 2022 Dissipation-range fluid turbulence and thermal noise. Phys. Rev. E 105 (6), 065113.
  * Bell et al. (2007) Bell, J.B., Garcia, A.L. & Williams, S.A. 2007 Numerical methods for the stochastic Landau-Lifshitz Navier-Stokes equations. Phys. Rev. E 76 (1), 016708.
  * Bell et al. (2022) Bell, J.B., Nonaka, A., Garcia, A.L. & Eyink, G. 2022 Thermal fluctuations in the dissipation range of homogeneous isotropic turbulence. J. Fluid Mech. 939.
  * Bird (1998) Bird, G.A. 1998 Molecular gas dynamics and the direct simulation of gas flows. Clarendon Press, Oxford.
  * Boffetta & Musacchio (2001) Boffetta, G. & Musacchio, S. 2001 Predictability of the inverse energy cascade in 2D turbulence. Phys. Fluids 13 (4), 1060–1062.
  * Boffetta & Musacchio (2017) Boffetta, G. & Musacchio, S. 2017 Chaos and predictability of homogeneous-isotropic turbulence. Phys. Rev. Lett. 119 (5), 054102.
  * Chandrasekhar (1961) Chandrasekhar, S. 1961 Hydrodynamic and Hydromagnetic Stability. Oxford University Press, London and New York.
  * Cho & Polvani (1996) Cho, J. Y-K. & Polvani, L.M. 1996 The morphogenesis of bands and zonal winds in the atmospheres on the giant outer planets. Science 273 (5273), 335–337.
  * Crane (2017) Crane, L. 2017 Infamous three-body problem has over a thousand new solutions. New Scientist (https://www.newscientist.com/article/2148074-infamous-three-body-problem-has-over-a-thousand-new-solutions/).
  * Curry et al. (1984) Curry, J.H., Herring, J.R., Loncaric, J. & Orszag, S.A. 1984 Order and disorder in two- and three-dimensional Bénard convection. J. Fluid Mech. 147, 1–38.
  * Diamond et al. (2005) Diamond, P.H., Itoh, S.I., Itoh, K. & Hahm, T.S. 2005 Zonal flows in plasma-a review. Plasma Phys. Control. Fusion 47 (5).
  * Donev et al. (2011) Donev, A., Bell, J.B., De la Fuente, A. & Garcia, A.L. 2011 Enhancement of diffusive transport by non-equilibrium thermal fluctuations. J. Stat. Mech.-Theory Exp. 2011 (06), P06014.
  * Donev et al. (2014) Donev, A., Nonaka, A., Sun, Y., Fai, T.G., Garcia, A.L. & Bell, J.B. 2014 Low Mach number fluctuating hydrodynamics of diffusively mixing fluids. Commun. Appl. Math. Comput. Sci. 9 (1), 47–105.
  * Donev et al. (2010) Donev, A., Vanden-Eijnden, E., Garcia, A. & Bell, J. 2010 On the accuracy of finite-volume schemes for fluctuating hydrodynamics. Commun. Appl. Math. Comput. Sci. 5 (2), 149–197.
  * Eyink & Jafari (2022) Eyink, G. & Jafari, A. 2022 High Schmidt-number turbulent advection and giant concentration fluctuations. Phys. Rev. Research 4 (2), 023246.
  * Frisch (1986) Frisch, U. 1986 Fully developed turbulence: where do we stand? In Dynamical Systems—A Renewal of Mechanism: Centennial of Georges David Birkhoff, pp. 13–28. World Scientific.
  * Frisch (1995) Frisch, U. 1995 Turbulence: The Legacy of AN Kolmogorov. Cambridge University Press.
  * Fromm (1965) Fromm, J.E. 1965 Numerical solutions of the nonlinear equations for a heated fluid layer. Phys. Fluids 8 (10), 828–828.
  * Gallis et al. (2021) Gallis, M.A., Torczynski, J.R., Krygier, M.C., Bitter, N.P. & Plimpton, S.J. 2021 Turbulence at the edge of continuum. Phys. Rev. Fluids 6, 013401.
  * Goluskin et al. (2014) Goluskin, D., Johnston, H., Flierl, G.R. & Spiegel, E.A. 2014 Convectively driven shear and decreased heat flux. J. Fluid Mech. 759 (6), 360–385.
  * Graham (1974) Graham, R. 1974 Hydrodynamic fluctuations near the convection instability. Phys. Rev. A 10 (5), 1762.
  * von Hardenberg et al. (2015) von Hardenberg, J., Goluskin, D., Provenzale, A. & Spiegel, E.A. 2015 Generation of large-scale winds in horizontally anisotropic convection. Phys. Rev. Lett. 115 (13), 134501.
  * Heimpel et al. (2005) Heimpel, M., Aurnou, J. & Wicht, J. 2005 Simulation of equatorial and high-latitude jets on Jupiter in a deep convection model. Nature 438 (7065), 193–6.
  * Heslot et al. (1987) Heslot, F., Castaing, B. & Libchaber, A. 1987 Transitions to turbulence in helium gas. Phys. Rev. A 36 (12), 5870.
  * Hu & Liao (2020) Hu, T.L. & Liao, S.J. 2020 On the risks of using double precision in numerical simulations of spatio-temporal chaos. J. Comput. Phys. 418, 109629.
  * Huang & Zhou (2013) Huang, Y.X & Zhou, Q. 2013 Counter-gradient heat transport in two-dimensional turbulent Rayleigh-Bénard convection. J. Fluid Mech. 737 (R3), 1–12.
  * Huisman et al. (2014) Huisman, S.G., van der Veen, R.C.A., Sun, C. & Lohse, D. 2014 Multiple states in highly turbulent Taylor-Couette flow. Nat. Commun. 5 (1), 1–5.
  * Johnston & Doering (2009) Johnston, H. & Doering, C.R. 2009 Comparison of turbulent thermal convection between conditions of constant temperature and constant flux. Phys. Rev. Lett. 102 (6), 064501.
  * Kadanoff (2001) Kadanoff, L.P. 2001 Turbulent heat flow: structures and scaling. Phys. Today 54 (8), 34–39.
  * Kaspi et al. (2018) Kaspi, Y., Galanti, E., Hubbard, W.B., Stevenson, D.J., Bolton, S.J., Iess, L., Guillot, T., Bloxham, J., Connerney, J.E.P. & Cao, H. 2018 Jupiter’s atmospheric jet streams extend thousands of kilometres deep. Nature 555 (7695), 223–226.
  * Knobloch & Weiss (1989) Knobloch, E. & Weiss, J.B. 1989 Effect of noise on discrete dynamical systems with multiple attractors. In Noise in Nonlinear Dynamical Systems, , vol. 2, pp. 65–86. Cambridge University Press.
  * Kraut et al. (1999) Kraut, S., Feudel, U. & Grebogi, C. 1999 Preference of attractors in noisy multistable systems. Phys. Rev. E 59 (5), 5253.
  * Landau & Lifshitz (1959) Landau, L.D. & Lifshitz, E.M. 1959 Fluid Mechanics, volume 6 of Course of Theoretical Physics. Pergamon.
  * Leith & Kraichnan (1972) Leith, C.E. & Kraichnan, R.H. 1972 Predictability of turbulent flows. J. Atmos. Sci. 29 (6), 1041–1058.
  * Li & Yorke (1975) Li, T.Y. & Yorke, J.A. 1975 Period three implies chaos. Am. Math. Mon. 82 (10), 985–992.
  * Li et al. (2018) Li, X.M., Jing, Y.P. & Liao, S.J. 2018 Over a thousand new periodic orbits of a planar three-body system with unequal masses. Publ. Astron. Soc. Jpn. 70 (4), 64.
  * Li et al. (2021) Li, X.M., Li, X.C. & Liao, S.J. 2021 One family of 13315 stable periodic orbits of non-hierarchical unequal-mass triple systems. Sci. China Phys. Mech. Astron. 64, 219511.
  * Li & Liao (2017) Li, X.M. & Liao, S.J. 2017 More than six hundred new families of Newtonian periodic planar collisionless three-body orbits. Sci. China-Phys. Mech. Astron. 60 (12), 129511.
  * Li & Liao (2019) Li, X.M. & Liao, S.J. 2019 Collisionless periodic orbits in the free-fall three-body problem. New Astron. 70, 22–26.
  * Liao (2009) Liao, S.J. 2009 On the reliability of computed chaotic solutions of non-linear differential equations. Tellus Ser. A-Dyn. Meteorol. Oceanol. 61 (4), 550–564.
  * Liao (2013) Liao, S.J. 2013 On the numerical simulation of propagation of micro-level inherent uncertainty for chaotic dynamic systems. Chaos Solitons Fractals 47, 1–12.
  * Liao (2014) Liao, S.J. 2014 Physical limit of prediction for chaotic motion of three-body problem. Commun. Nonlinear Sci. Numer. Simul. 19 (3), 601–616.
  * Liao et al. (2022) Liao, S.J., Li, X.M. & Yang, Y. 2022 Three-body problem: from Newton to supercomputer plus machine learning. New Astronomy 96, 101850.
  * Liao & Qin (2022) Liao, S.J. & Qin, S.J. 2022 Ultra-chaos: an insurmountable objective obstacle of reproducibility and replication. Adv. Appl. Math. Mech. 14 (4), 799–815.
  * Liao & Wang (2014) Liao, S.J. & Wang, P.F. 2014 On the mathematically reliable long-term simulation of chaotic solutions of lorenz equation in the interval [0, 10000]. Sci. China-Phys. Mech. Astron. 57 (2), 330–335.
  * Lin et al. (2017) Lin, Z.L., Wang, L.P. & Liao, S.J. 2017 On the origin of intrinsic randomness of Rayleigh-Bénard turbulence. Sci. China-Phys. Mech. Astron. 60 (1), 1–13.
  * Lohse & Xia (2010) Lohse, D. & Xia, K.Q. 2010 Small-scale properties of turbulent Rayleigh-Bénard convection. Annu. Rev. Fluid Mech. 42, 335–364.
  * Lorenz (1963) Lorenz, E.N. 1963 Deterministic nonperiodic flow. J. Atmos. Sci. 20 (2), 130–141.
  * Lorenz (1989) Lorenz, E.N. 1989 Computational chaos-a prelude to computational instability. Physica D 35 (3), 299–317.
  * Lorenz (1993) Lorenz, E.N 1993 The Essence of Chaos. University of Washington Press.
  * Lorenz (2006) Lorenz, E.N. 2006 Computational periodicity as observed in a simple system. Tellus Ser. A – Dyn. Meteorol. Oceanol. 58 (5), 549–557.
  * Masoller (2002) Masoller, C. 2002 Noise-induced resonance in delayed feedback systems. Phys. Rev. Lett. 88 (3), 034102.
  * Maximenko et al. (2005) Maximenko, N.A., Bang, B. & Sasaki, H. 2005 Observational evidence of alternating zonal jets in the world ocean. Geophys. Res. Lett. 32 (12).
  * McMullen et al. (2022) McMullen, R.M., Krygier, M.C., Torczynski, J.R. & Gallis, M.A. 2022 Navier-Stokes equations do not describe the smallest scales of turbulence in gases. Phys. Rev. Lett. 128 (11), 114501.
  * Miyagoshi et al. (2010) Miyagoshi, T., Kageyama, A. & Sato, T. 2010 Zonal flow formation in the Earth’s core. Nature 463 (7282), 793.
  * Moore & Weiss (1973) Moore, D.R. & Weiss, N.O. 1973 Two-dimensional Rayleigh-Bénard convection. J. Fluid Mech. 58 (2), 289–312.
  * Niemela & Sreenivasan (2006) Niemela, J.J. & Sreenivasan, K.R. 2006 Turbulent convection at high Rayleigh numbers and aspect ratio 4. J. Fluid Mech. 557, 411–422.
  * Parker & Chua (1989) Parker, T.S. & Chua, L.O. 1989 Practical numerical algorithms for chaotic systems. Math. Comput. 66 (1), 125–8.
  * Peter (1998) Peter, S. 1998 Explaining Chaos. Cambridge University Press.
  * Poincaré (1890) Poincaré, H. 1890 Sur le problème des trois corps et les équations de la dynamique. Acta Math. 13 (1), A3–A270.
  * van der Pol et al. (2014) van der Pol, E.P., Ostilla-Mónico, R., Verzicco, R. & Lohse, D. 2014 Effect of velocity boundary conditions on the heat transfer and flow topology in two-dimensional Rayleigh-Bénard convection. Phys. Rev. E 90 (1), 13017–13017.
  * Pope (2001) Pope, S.B. 2001 Turbulent Flows. IOP Publishing.
  * Qin & Liao (2020) Qin, S.J. & Liao, S.J. 2020 Influence of numerical noises on computer-generated simulation of spatio-temporal chaos. Chaos Solitons Fractals 136, 109790.
  * Rayleigh (1916) Rayleigh, L. 1916 On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side. Philos. Mag. 32 (192), 529–546.
  * Richards et al. (2006) Richards, K.J., Maximenko, N.A., Bryan, F.O. & Sasaki, H. 2006 Zonal jets in the Pacific Ocean. Geophys. Res. Lett. 33 (3), 155–170.
  * Saltzman (1962) Saltzman, B. 1962 Finite amplitude free convection as an initial value problem I. J. Atmos. Sci. 19 (4), 329–341.
  * Schlüter et al. (1965) Schlüter, A., Lortz, D. & Busse, F. 1965 On the stability of steady finite amplitude convection. J. Fluid Mech. 23 (01), 129.
  * Shrimali et al. (2008) Shrimali, M.D., Prasad, A., Ramaswamy, R. & Feudel, U. 2008 The nature of attractor basins in multistable systems. Int. J. Bifurcation Chaos 18 (06), 1675–1688.
  * Smith (2015) Smith, E.R. 2015 A molecular dynamics simulation of the turbulent Couette minimal flow unit. Phys. Fluids 27 (11), 115105\.
  * de Souza et al. (2007) de Souza, S.L.T., Batista, A.M., Caldas, I.L., Viana, R.L. & Kapitaniak, T. 2007 Noise-induced basin hopping in a vibro-impact system. Chaos Solitons Fractals 32 (2), 758–767.
  * Sprott (2010) Sprott, J.C. 2010 Elegant Chaos: Algebraically Simple Chaotic Flows. World Scientific, Singapore.
  * Sugiyama et al. (2009) Sugiyama, K., Calzavarini, E., Grossmann, S. & Lohse, D. 2009 Flow organization in two-dimensional non-Oberbeck-Boussinesq Rayleigh-Bénard convection in water. J. Fluid Mech. 637, 105–135.
  * Sun et al. (1993) Sun, Z.P., Schubert, G. & Glatzmaier, G.A. 1993 Banded surface flow maintained by convection in a model of the rapidly rotating giant planets. Science 260 (5108), 661–664.
  * Swift & Hohenberg (1977) Swift, J. & Hohenberg, P.C. 1977 Hydrodynamic fluctuations at the convective instability. Phys. Rev. A 15 (1), 319.
  * Teixeira et al. (2007) Teixeira, J., Reynolds, C.A. & Judd, K. 2007 Time step sensitivity of nonlinear atmospheric models: numerical convergence, truncation error growth, and ensemble design. J. Atmos. Sci. 64 (1), 175–189.
  * Veronis (1968) Veronis, G. 1968 Large-amplitude Bénard convection. J. Fluid Mech. 31 (1), 113–139.
  * Wang et al. (2020) Wang, Q., Chong, K.L., Stevens, R.J.A.M., Verzicco, R. & Lohse, D. 2020 From zonal flow to convection rolls in Rayleigh-Bénard convection with free-slip plates. J. Fluid Mech. 905, A21.
  * Whyte (2018) Whyte, C. 2018 Watch the weird new solutions to the baffling three-body problem. New Scientist (https://www.newscientist.com/article/2170161-watch-the-weird-new-solutions-to-the-baffling-three-body-problem/).
  * Xu et al. (2021) Xu, T.Z., Li, J.and Li, Z.H. & Liao, S.J. 2021 Accurate predictions of chaotic motion of a free fall disk. Phys. Fluids 33 (3), 037111.
  * Yano et al. (2003) Yano, J.I., Talagrand, O. & Drossart, P. 2003 Outer planets: origins of atmospheric zonal winds. Nature 421 (6918), 36.
  * Zhang et al. (2017) Zhang, Y., Zhou, Q. & Sun, C. 2017 Statistics of kinetic and thermal energy dissipation rates in two-dimensional turbulent Rayleigh-Bénard convection. J. Fluid Mech. 814, 165–184.
  * Zhou & Xia (2013) Zhou, Q. & Xia, K.Q. 2013 Thermal boundary layer structure in turbulent Rayleigh-Bénard convection in a rectangular cell. J. Fluid Mech. 721, 199–224.
  * Zhu et al. (2018) Zhu, X., Mathai, V., Stevens, R., Verzicco, R. & Lohse, D. 2018 Transition to the ultimate regime in two-dimensional Rayleigh-Bénard convection. Phys. Rev. Lett. 120 (14), 144502.1–144502.6.
  * Zienicke et al. (1998) Zienicke, E., Seehafer, N. & Feudel, F. 1998 Bifurcations in two-dimensional Rayleigh-Bénard convection. Phys. Rev. E 57 (1), 428.


[a,b]K. K. Yadav

# Status update of the MACE Gamma-ray telescope

###### Abstract

MACE (Major Atmospheric Cherenkov Experiment), an imaging atmospheric
Cherenkov telescope, has recently been installed by the HiGRO (Himalayan
Gamma-Ray Observatory) collaboration at Hanle (32.8∘N, 78.9∘E, 4270m asl) in
Ladakh region of North India. The telescope has a 21m diameter large light
collector consisting of indigenously developed 1424 square-shaped diamond
turned spherical aluminum mirror facets of size $\sim$ 0.5m$\times$0.5m. MACE
is the second largest Cherenkov telescope at the highest altitude in the
northern hemisphere. The imaging camera of the telescope consists of 1088
photo-multiplier tubes with a uniform pixel resolution of $\sim 0.125^{\circ}$
covering a field of view of $\sim$ 4.0∘ $\times$ 4.0∘. The main objective of
the MACE telescope is to study gamma-ray sources mainly in the unexplored
energy region 20 -100 GeV and beyond with high sensitivity. In this paper, we
describe the key design features and current status of MACE including results
from the trial observations of the telescope.

##  1 Introduction

MACE (Major Atmospheric Cherenkov Experiment), a very high energy (VHE)
$\gamma$-ray telescope based on imaging atmospheric Cherenkov technique, is
currently under commissioning at Hanle, India [1, 2]. The Indian astronomical
site at Hanle (32.8∘N, 78.9∘E, 4270 m asl) in the Himalayan range is at the
highest altitude in the world for any existing imaging atmospheric Cherenkov
telescope (IACT) at present. The main objective of the MACE telescope is the
observation of the VHE $\gamma$-ray photons with energy 30 GeV and above,
which is traditionally considered as the window for space-based telescopes.
Observations with the MACE telescope will help in understanding the production
and propagation of energetic $\gamma$-ray photons in the non-thermal Universe
and to also explore some phenomena which are not accessible to terrestrial
accelerators. Some of the important astrophysical motivations for the MACE
telescope include: discovery of new $\gamma$-ray sources, identification of
unassociated sources in the high energy _Fermi_ catalog, detection of the
$\gamma$-ray emission from various types of known sources (active galactic
nuclei, gamma-ray bursts, pulsars and so on) at different ages of the
Universe, providing evidence for the enigmatic jet formation process etc.
Beyond probing the non-thermal Universe and cosmic accelerators, the VHE
observations using the MACE telescope are also expected to address a range of
cosmological topics such as measurement of the intensity of extragalactic
background light, cosmic ray electron spectrum, search for nature of dark
matter candidates like weakly interacting massive particles [3, 4, 5].

##  2 The MACE Telescope: Salient Features

The MACE telescope is equipped with a 21m diameter quasi-parabolic reflector
with $f/d$ ratio of $\sim$ 1.2. The quasi-parabolic design of such a large
reflector helps in reducing the optical aberrations of the telescope. To
achieve such a large light collector area, reflector of the MACE telescope is
segmented into 1424 small, square-shaped spherical mirror facets of size
0.488m $\times$ 0.488m each with varying focal lengths (Figure 1(Left)). Four
such mirrors with similar focal length are mounted on a single panel of size
0.986m $\times$ 0.986m each and a total of 356 mirror panels with varying
focal length in the range 25m to 26.2m are obtained. All the 356 mirror panels
mounted on the telescope basket are manually aligned in such a way that the
resulting surface behaves like a single quasi-parabolic reflector with focal
length gradually increasing from the center of the basket towards the
periphery. An active mirror alignment control system is used for orienting
each mirror panel to achieve the desired optical quality of the reflector. The
metallic mirror facets are made up of Aluminium alloy supported by a honeycomb
structure. Diamond turning technology [6] has been used to fabricate all the
mirror facets in order to achieve desired mirror quality in terms of surface
finish and accuracy within the country (Figure 1(Left)). The reflectance of
these mirror facets is greater than 85$\%$ in the wavelength range 280-700 nm.
Due to its large size, the MACE telescope is not protected by a dome and the
mirrors are continuously exposed to the environment. The reflecting surface of
each mirror facet is coated with a thin (100-150 nm) layer of SiO2 for
protecting the reflecting surface and ensuring its longevity.

Figure 1: Left: A diamond turned mirror facet for the MACE light collector.
Right: 16-channel camera integrated module for the camera of MACE telescope.

The telescope deploys a 1088-pixel imaging camera with resolution of 0.125∘ at
its focal plane, covering a field of view of $4.36^{\circ}\times
4.03^{\circ}$. Six stage photomultiplier tubes (pixel) with 38 mm diameter
arranged at a triangular pitch of 55 mm are used in the MACE-camera for
detection of the Cherenkov photons. Each photomultiplier tube is provided with
a hexagoanl front-coated light concentrator in order to enhance the light
collection efficiency by collecting the Cherenkov photons incident in the dead
space between adjacent pixels. The light collection efficiency of these
compound parabolic concentrators is more than 85$\%$ in the wavelength range
of 230-600 nm. The MACE-camera is modular in design and consists of 68 camera
integrated modules (CIMs) of 16 channels each. Each module channel has its
signal processing electronics (photomultiplier tubes with programmable high
voltage generators, pre-amplifiers, amplifiers, discriminators, first level
trigger generation logic, and signal digitization circuitry) housed in it
(Figure 1(Right)). An analogue switched capacitor array DRS-4 is used as ring
sampler at 1 GHz speed for continuous digitization of the signal from the
photomultiplier tubes. The signal from the photomultiplier tubes is
simultaneously amplified at low and high gain of 14 and 140 respectively to
ensure a large dynamic range. The discriminator output amplitude of each
channel is used for monitoring its single channel rate and also for generating
the first-level trigger from an individual CIM. The first-level triggers from
all the modules are collated in a second-level trigger generator where
proximity of the triggered pixels in adjacent modules is checked. The
innermost 576 pixels (24$\times$24) are used for generating the trigger
according to predefined logic for nearest pairs, triplets, quadruplets, etc of
the pixels. After the generation of the second-level trigger, the data from
all the 68 modules are collected by the data concentrator, which in turn sends
them to the data acquisition computer in the control room through optical
fibres. It is expected that about 50GB of data will be stored during every
hour of observation with the MACE telescope [7, 8].

An altitude-azimuth mount drive system of the MACE telescope provides the
stability of 21m diameter large mechanical structure with 180 ton weight by
using a track and wheel design. Two azimuth drive wheels are coupled to three
phase, permanent magnet brushless AC servo motors through multi-stage
gearboxes for providing the azimuth motion. The elevation movement is provided
through a gearbox coupled to a 13-section bull-gear assembly of 11.6m radius.
All the drives are provided with counter torque capability to avoid gear
backlash error. The motors are driven by pulse width modulated drive
amplifiers powered by 480 Volt DC from a solar power station. The positions of
the two axes are monitored by a 25-bit absolute optical encoders with $\sim$
20 arcsec accuracy. Both the azimuth and elevation gear boxes also have high-
speed options to move the telescope at speed of 3${}^{\circ}s^{-1}$ to quickly
point the telescope in the direction of transient events like gamma-ray bursts
in the sky. The MACE drive system can provide tracking accuracy of better than
1 arcmin in wind speeds of up to 30 km h-1. The telescope is automatically
brought to the parking position if the sustained wind speed is more than 40 km
h-1. Moreover, the telescope structure is designed for survival at a wind
speed of 150 km h-1.

Figure 2: The MACE telescope installed at Hanle-Ladakh, India.

Figure 3: Flowchart for the MACE data Analysis Package (MAP).

##  3 Present Status of the MACE Telescope

The installation of the MACE telescope at Hanle site was successfully
completed in October 2020. A recent picture of the telescope at the site is
depicted in Figure 2. The trial runs for testing the performance of different
components of the fully assembled telescope started in November 2020. Monte
Carlo simulation studies carried out using the CORSIKA package suggest that
the MACE telescope is expected to have a $\gamma$-ray trigger energy threshold
of $\sim$ 20 GeV in the low zenith angle range below 40∘ [9]. At high zenith
angles above 40∘, the threshold energy increases to $\sim$ 150 GeV. The
dynamic energy range of the MACE telescope is expected to be 20 GeV-5 TeV with
a total trigger rate of $\sim$ 1 kHz. The integral sensitivity for a point
source with Crab Nebula like spectrum above 30 GeV is 2.7$\%$ at 5$\sigma$
statistical significance level in 50 hours of observation [10]. The angular
resolution of MACE telescope is estimated as $\sim$ 0.21∘ in the energy range
of 30-50 GeV and this improves to $\sim$ 0.06∘ in the energy band 1.8-3 TeV
[11]. The telescope is expected to have energy resolutions of $\sim$ 40$\%$
and $\sim$ 19$\%$ in the energy bands 30-50 GeV and 1.8-3 TeV respectively.
For an effective analysis of the real time data collected with the MACE
telescope, a data analysis software called MAP (MACE Analysis Package) has
been developed. The MAP represents collection of ROOT-based programs/routines
written in C-language under LINUX platform. Different analysis steps included
in the MAP are outlined in Figure 3. Another version of the MAP based on
Python libraries is also under development.

Figure 4: Left: The point spread function of the MACE telescope measured for
Polaris. Right: Hit pattern uniformity in the trigger region of the camera
during the Crab Nebula observations for $\sim$ 1 hour on April 1, 2021.

##  4 Preliminary Results

The MACE telescope, currently under its commissioning phase, has been deployed
for special trail runs since October 2020. Alignment of the 356 mirror panels
on telescope basket is carried out with the help of an active mirror alignment
control system while tracking the pole star. A CCD camera along with 180mm
focal length lens mounted on the center of the telescope basket is used to
capture the reflected image of an optical source on telescope focal plane. The
image of pole star captured by the CCD camera on the focal plane after
background subtraction and image cleaning is shown in Figure 4(Left). D80
(defined as the diameter of circle, concentric with the centroid of the image,
which contains 80$\%$ intensity of the image) of the image is found to be 46.5
mm after performing a Gaussian fit to the pole star image data. This
characterizes the point sperad function of the telescope.

In order to record the Cherenov event data, the information of hit channels in
the camera of the telescope is also stored. The number of hits for each
triggered pixel in the trigger region (inner 576 pixels) of the camera during
the trial run in the direction of the Crab Nebula on the night of April 1,
2021 for about 1 hour of observation is reported in Figure 4(Right). This
indicates the homogeneity of the camera of the MACE telescope in generating
the trigger in response to the Cherenkov events.

Preliminary analysis of the data collected on the Crab Nebula on April 1, 2021
has been performed using the standard Hillas Parameterization method [12]
implemented in the MAP. The distribution of the Hillas orientation parameter
($\alpha$) for the recorded Cherenkov events is depicted in Figure 5. The
$\alpha$-plot shown here only represents the initial performance of the MACE
telescope without standard data quality checks for IACT observations. Although
there is an indication of the presence of $\gamma$-ray like events in the
source region at $\sim$ 5.4$\sigma$ statistical significance level in Figure
5, yet these preliminary results in the present form cannot be referred to as
the first light of the MACE telescope. A detailed analysis of the data
collected on the Crab Nebula with MACE in the month of April 2021 is currently
underway and the final results of the first light are expected very soon.

Figure 5: Preliminary results from the Crab Nebula observations on April 1,
2021 for a live time of $\sim$ 1 hour.

##  5 Summary

A group of very high energy $\gamma$-ray astronomers are actively involved in
setting up the MACE telescope in India under HiGRO (Himalayan Gamma Ray
Observatory) Collaboration. The construction and installation of the telescope
at Hanle-site has been successfully completed. Currently, the MACE telescope
is under its commissioning phase and its first light is expected very soon.
Subsequently, the telescope will enter into the science operatione phase. Once
fully operational, the MACE telescope, equipped with a 21m diameter light
collector and 1088-pixel camera, will be the second largest $\gamma$-ray
telescope (after the Large Size Telescope prototype (LST-1) at the Cherenkov
Telescope Array) in the northern hemisphere. However, altitude of the MACE-
site is the highest in the world. This geographical advantage along with the
large light collector will help in exploring the $\gamma$-ray sky above 15GeV
and also provide an excellent energy overlap with the space-based $\gamma$-ray
observatories.

## Acknowledgements

The HiGRO collaboration thanks all the colleagues from various Divisions at
Bhabha Atomic Research Centre, who have been involved at different stages of
the construction and installation of the MACE telescope at Hanle, India.

## References

  * [1] R. Koul, _MACE gamma-ray telescope-a status update_ , _CURRENT SCIENCE_ 113 (2017) 691
  * [2] K. K. Singh, & K. K. Yadav, _20 Years of Indian Gamma Ray Astronomy Using Imaging Cherenkov Telescopes and Road Ahead_ , _Universe_ 7 (2021) 96
  * [3] F. Aharonian et al., _High energy astrophysics with ground-based gamma ray detectors_ , _Reports on Progress in Physics_ 71 (2008) 096901
  * [4] A. M. Hillas, _Evolution of ground-based gamma-ray astronomy from the early days to the Cherenkov Telescope Arrays_ , _Astroparticle Physics_ 43 (2013) 19
  * [5] G. Di Sciascio, _Ground-based Gamma-Ray Astronomy: an Introduction_ , _Journal of Physics Conference Series_ 1263 (2019) 012003
  * [6] T. Sugano et al., _Diamond Turning of an Aluminum Alloy for Mirror_ , _CIRP Annals_ 36 (1987) 17
  * [7] S. Srivastava et al., _MACE camera controller embedded software: Redesign for robustness and maintainability_ , _Astronomy and Computing_ 30 (2020) 100358
  * [8] D. Sarkar et al., _Comparison of data storage and analysis throughput in the light of high energy physics experiment MACE_ , _Astronomy and Computing_ 33 (2020) 100409
  * [9] C. Borwankar et al., _Simulation studies of MACE-I: Trigger rates and energy thresholds_ , _Astroparticle Physics_ 84 (2016) 97
  * [10] M. Sharma et al., _Sensitivity estimate of the MACE gamma ray telescope_ , _Nuclear Instruments and Methods in Physics Research A_ 851 (2017) 125
  * [11] C. Borwankar et al., _Estimation of expected performance for the MACE $\gamma$-ray telescope in low zenith angle range_, _Nuclear Instruments and Methods in Physics Research A_ 953 (2020) 163182
  * [12] A. M. Hillas, _Cerenkov Light Images of EAS Produced by Primary Gamma Rays and by Nuclei_ , _19th International Cosmic Ray Conference_ 3 (1985) 445

## Full Authors List: HiGRO Collaboration

N. Bhatt1, S. Bhattacharyya1,2, C. Borwankar1, K. Chanchalani1, P. Chandra1,
V. R. Chitnis3, N. Chouhan1, M. P. Das1, V. K. Dhar1, B. Ghosal1, S. Godambe1,
S. Godiyal1, K. K. Gour1, J. Hariharan1, M. Khurana1, M. Kothari1, S. Kotwal1,
M. K. Koul1, N. Kumar1, N. Kumar1, C. P. Kushwaha1, N. Mankuzhiyil1, P.
Marandi1, K. Nand1, S. Norlha1, D. Sarkar1, M. Sharma1, K. K. Singh1, R.
Thubstan1, A. Tolamatti1, K. Venugopal1, K. K. Yadav1,2

1Astrophysical Sciences Division, Bhabha Atomic Research Centre, Mumbai
400094, India.
2Homi Bhabha National Institute, Anushakti Nagar,Mumbai 400094, India.
3Department of High Energy Physics, Tata Institute of Fundamental Research,
Homi Bhabha Road, Colaba, Mumbai 400005, India.