# Top quark pair production and top quark properties at CDF

Chang-Seong Moon, on behalf of the CDF Collaboration  AstroParticule et
Cosmologie, Université Paris Diderot, CNRS/IN2P3,
10, rue Alice Domon et Léonie Duquet, 75205 Paris, France
and INFN Sezione di Pisa, Largo B. Pontecorvo, 3, 56127 Pisa, Italy

###### Abstract

We present the most recent measurements of top quark pairs production and top
quark properties in proton-antiproton collisions with center-of-mass energy of
1.96 TeV using CDF II detector at the Tevatron. The combination of top pair
production cross section measurements and the direct measurement of top quark
width are reported. The test of Standard Model predictions for top quark
decaying into $b$-quarks, performed by measuring the ratio $R$ between the top
quark branching fraction to $b$-quark and the branching fraction to any type
of down quark is shown. The extraction of the CKM matrix element $|V_{tb}|$
from the ratio $R$ is discussed. We also present the latest measurements on
the forward-backward asymmetry ($A_{FB}$) in top anti-top quark production.
With the full CDF Run II data set, the measurements are performed in top anti-
top decaying to final states that contain one or two charged leptons
(electrons or muons). In addition, we combine the results of the leptonic
forward-backward asymmetry in $t\bar{t}$ system between the two final states.
All the results show deviations from the next-to-leading order (NLO) standard
model (SM) calculation.

###### keywords:

Tevatron, CDF, Standard Model, Top quark, Properties, forward-backward
asymmetry

##  1 Introduction

The top quark, discovered at the Tevatron in 1995 by CDF [1] and D0 [2]
experiments, have completed the third generation of quarks in the SM and been
playing a special role with the heaviest known elementary particle in
electroweak symmetry breaking. This large mass is responsible for the top
quark properties. Any deviation from the SM predictions of the top quark
properties could indicate the presence of beyond SM signal indirectly.

Top quarks are mostly produced in top-antitop pairs via the strong interaction
at the Tevatron. Thousands of top events have been collected by CDF Run II
detector. It allows performing precision measurements of top properties and
many analyses are unique to the Tevatron and complementary to the LHC
measurements.

In the SM, the top quark decays to a $W$-boson and a $b$-quark almost 100 % of
the time. The decay channels of the top quark are labeled as dilepton,
lepton+jets and all hadronic channel classified by whether a leptonical decay
into electron or muon has occurred in both (5 %), one only (30 %) or none of
the two $W$-bosons (45 %) respectively. The dilepton channel is only 5 percent
of top quark pairs decay but it has the cleanest final state to identify
because there is no so many other processes can produce such a striking final
state even before requiring the identification of one of the jets originating
from a $b$-quark ($b$-tagging). The lepton+jets mode is referred to as the
golden channel because of a sizable branching fraction and controllable
background levels using the $b$-tagging. The all hadronic decay channel has an
advantage of the largest cross section but it suffers from a huge background
from the quantum chromodynamics (QCD) multijet production and cannot be easily
simulated with Monte Carlo sample.

##  2 Top quark width and branching fraction ratio

Figure 1:  (top) $m^{reco}_{t}$ distributions for simulated events meeting the
lepton + jets selection are displayed with three values of top and with the
nominal jet-energy scale. (middle) Confidence bands of $\Gamma_{top}$ as a
function of $\Gamma_{meas}$ for 68 % and 95 % CL limits. Results from
simulated experiments assuming 8.7 fb-1 of data at different values of top are
convoluted with a smearing function to account for systematic uncertainties.
The value observed in data is indicated by an arrow. (bottom) Number of events
observed in data and expected for various values of $R$ as a function of
identified $b$-jets.

Top quark has the largest decay width of the known fermions in the SM. The
recent next-to-next-to-leading-order (NNLO) calculation with QCD and
electroweak corrections predicts $\Gamma_{top}$, the top quark decay width to
be 1.32 GeV at top quark mass of $M_{top}$ = 172.5 GeV/c2 [3]. Direct
measurements of the top quark width are performed in fully reconstructed
lepton+jets events using the full CDF Run II data set, corresponding to an
integrated luminosity of 8.7 fb-1. The best-fit value of measured top quark
width ($\Gamma_{meas}$) retrieved from the data is 1.63 GeV and is depicted as
an arrow in the plot of Fig. 1 (middle), This corresponds to an upper limit of
$\Gamma_{top}<$ 6.38 GeV at the 95 % confidence level (CL) We also set a two-
sided limit of 1.10 GeV $<\Gamma_{top}<$4.05 GeV at 68 % CL, which corresponds
to a lifetime of 1.6$\times$10-25 s, this result confirms the prediction that
top quark decays before its hadronization. This is the most precise direct
determination of the top quark width and lifetime and shows no evidence of
non-SM physics in the top quark decay [4].

In the SM, the top quark coupling to a down-type quark $q$, where $q$ = $d$,
$s$, $b$ is proportional to $|V_{tq}|$, the element of the Cabibbo-Kobayashi-
Maskawa (CKM) matrix. Under the assumption that the 3$\times$3 CKM matrix is
unitary and given constraints on $|V_{ts}|$ and $|V_{td}|$, the magnitude of
the top-bottom quark coupling is computed as $|V_{tb}|$ = 0.99915
${}^{+0.00002}_{-0.00005}$ [5]. Using these hypotheses, the ratio of the top
quark branching fraction $R$ is indirectly determined as

$\begin{split}R=\frac{B(t\rightarrow Wb)}{B(t\rightarrow
Wq)}&=\frac{|V_{tb}|^{2}}{|V_{tb}|^{2}+|V_{tb}|^{2}+|V_{tb}|^{2}}\\\
&=0.99830^{+0.00004}_{-0.00009}\end{split}$ (1)

The ratio of the top quark branching fraction $R$ is recently measured in
events with two charged leptons, imbalance in total transverse energy, and at
least two jets (dilepton channel). The sample is divided into 9 subsamples by
dilepton flavor ($ee$, $\mu\mu$, $e\mu$) and 3 $b$-tagging categories
(presence of 0, 1, or 2 tags) in order to better exploit the subsample-
dependent signal-to-background ratio and compare between observed data and
expectations. The number of observed and predicted events are used in the
various subsamples as input to a likelihood function, which is maximized to
extract $R$, since the number of $b$-jets in the event is related to the top
quark branching fraction in the $Wb$ final state. We obtain $R>$ 0.73(0.76) at
95 % (90 %) CL. From Eq. (1) and the assumptions therein we obtain $|V_{tb}|$
= 0.94 $\pm$ 0.04 and $|V_{tb}|>$ 0.85(0.87) at 95 % (90 %) CL. The result,
$R$ = 0.87 $\pm$ 0.07 differs from the SM expectation by $\approx$ 1.8$\sigma$
[6].

##  3 Top quark pair production cross section

|
---|---
Figure 2:  (left) The CDF-only combination result with the four input
$\sigma_{t\bar{t}}$ measurements from the CDF experiment, along with the
D0-only combination for the combination for the Tevatron result. The inner
(red) bars reflect statistical uncertainties while the outer (blue) bars show
the total uncertainties on each measurement. (right) Fraction of cross section
in 10 bins of cos$\theta_{t}$, obtained by integrating the series of Legendre
polynomials over the width of each bin.

CDF collaboration recently combined four measurements of the inclusive top
quark pair ($t\bar{t}$) production cross section ($\sigma_{t\bar{t}}$) from
the full CDF Run II data set corresponding to integrated luminosities of up to
8.8 fb-1 by contribution of a measurement in the dilepton channel [7]. We
obtain a value of $\sigma_{t\bar{t}}$ = 7.63 $\pm$ 0.50 pb for a top quark
mass of $M_{top}$ = 172.5 GeV. The contributions to the uncertainty are 0.31
pb from statistical sources, 0.39 pb from systematic sources including the
uncertainty on the integrated luminosity. The result is in good agreement with
the standard model expectation of 7.35 ${}^{+0.28}_{-0.33}$ pb at NNLO and
next-to-next-to leading logarithms in perturbative QCD [8]. The CDF-only
combination result contributes to 60 % of the combination for the Tevatron
result to be $\sigma_{t\bar{t}}$ = 7.60 $\pm$ 0.41 pb [9].

A measurement of the differential $t\bar{t}$ cross section,
$d\sigma/d(\text{cos}\theta_{t})$ for top quark pair production as a function
of the top quark production angle in proton-antiproton collisions is reported
[10]. This measurement is performed using the lepton+jets events and the full
CDF Run II data set corresponding to an integrated luminosity of 9.4 fb-1. We
employ the Legendre polynomials to characterize the shape of the differential
cross section at the parton-level. We characterize the shape of
$d\sigma/d(\text{cos}\theta_{t})$ by expanding in the Legendre polynomials
where $P_{\ell}$ is the Legendre polynomial of degree $\ell$, and $a_{\ell}$
is the Legendre moment of degree $\ell$ in Eq. (2).

$\frac{d\sigma}{d(\text{cos}\theta_{t})}=\sum^{\infty}_{\ell=0}a_{\ell}P_{\ell}(\text{cos}\theta_{t})$
(2)

The observed Legendre coefficient of the cos$\theta_{t}$ linear-term in the
differential cross section $a_{1}$ = 0.40 $\pm$ 0.12 is excess of the NLO SM
prediction of $a_{1}$ = 0.15 ${}^{+0.07}_{-0:03}$, while the remainder of the
differential cross section is in good agreement within the uncertainties with
the prediction of the NLO SM calculation [10].

##  4 Top quark forward-backward asymmetry

The forward-backward asymmetry ($A_{FB}$) in $p\bar{p}\rightarrow t\bar{t}$
production is defined as

$A_{FB}=\frac{N(\Delta y>0)-N(\Delta y<0)}{N(\Delta y>0)+N(\Delta
y<0)}~{}(\Delta y=y_{t}-y_{\bar{t}}),$ (3)

where $y$ is given by $y=\frac{1}{2}ln\left(\frac{E+p_{z}}{E-pz}\right)$ with
$E$ being the total top quark energy and $p_{z}$ being the component of the
top quark momentum along the beam axis as measured in the detector rest frame.
The $A_{FB}$ is identical to the asymmetry in the top-quark production angle
in the experimentally well-defined $t\bar{t}$ rest frame. An NLO SM
calculation with both electroweak and NLO QCD effects lead to a small $A_{FB}$
of order 8.8 $\pm$ 0.6 % in the $t\bar{t}$ rest frame [11]. Asymmetries in
$t\bar{t}$ production beyond SM prediction could provide the first evidence of
new interactions, such as $t\bar{t}$ production via a heavy axial color octet
or a flavor-changing Z′ boson that might not be easily observed as excesses in
the top-quark production rate or as resonances in the $t\bar{t}$ invariant
mass distribution [12].

CDF measured a reconstruction-level asymmetry of $A_{FB}$ = 6.3 $\pm$ 1.9 %
and an inclusive parton-level $A_{FB}$ = 16.4 $\pm$ 4.7 % in the lepton+jets
decay channel with the full CDF data set corresponding to an integrated
luminosity of 9.4 fb-1. In addition, differential measurements using two bins
each in the top-antitop rapidity difference $|\Delta y|$ and the top-antitop
invariant mass $M_{t\bar{t}}$ are performed. The results agree on a larger
dependence of $A_{FB}$ on $|\Delta y|$ and $M_{t\bar{t}}$ than the NLO SM
expectation. The observed slopes of the asymmetries by linear fit of the data
as a function of $|\Delta y|$ and $M_{t\bar{t}}$ give 2.8 $\sigma$ and 2.4
$\sigma$ excesses from the theoretical prediction respectively. We also study
the dependence of the asymmetry on the $t\bar{t}$ transverse momentum of the
system $p^{t\bar{t}}_{T}$ at the detector level. These results show the excess
asymmetry in the data is consistent with being dependent of $p^{t\bar{t}}_{T}$
and provide additional quantification of the functional dependencies of the
asymmetry [12].

##  5 Leptonic asymmetry in $t\bar{t}$ production

The leptonic forward-backward asymmetry $A^{\ell}_{FB}$ ($\ell$ = $e$ or
$\mu$) kinematically correlated with the top quark pair asymmetry
$A^{t\bar{t}}_{FB}$ as well as the top quark polarizations. If the genuine
value of $A^{t\bar{t}}_{FB}$ would be that measured by CDF collaboration [12],
the predicted value for $A^{\ell}_{FB}$ for top quarks decaying according to
the SM would be 7.0 % $<A^{\ell}_{FB}$ $<$ 7.6 % [13]. The inclusive asymmetry
of the charge-weighted lepton pseudorapidities ($q_{\ell}y_{\ell}$) is defined
as

$A^{\ell}_{FB}=\frac{N(q_{\ell}y_{\ell}>0)-N(q_{\ell}y_{\ell}<0)}{N(q_{\ell}y_{\ell}>0)+N(q_{\ell}y_{\ell}<0)}$
(4)

The asymmetric part is decomposed from $q_{\ell}y_{\ell}$ distribution at
parton-level with various physics models [13, 14]. The differential asymmetry
as a function of $q_{\ell}y_{\ell}$, $\mathcal{A}$($q_{\ell}y_{\ell}$) is
parametrized with an functional form of

$\mathcal{A}(q_{\ell}y_{\ell})=a\cdot\mathrm{tanh}\left(\frac{q_{\ell}y_{\ell}}{2}\right),$
(5)

where $a$ is the only free parameter related to $A^{\ell}_{FB}$. A predicted
value in the SM is $A^{\ell}_{FB}$ = 3.8 $\pm$ 0.3 % [11].

A second observable for the leptonic asymmetry is defined in the dilepton
final state analogously to $A^{t\bar{t}}_{FB}$ as

$A^{\ell\ell}_{FB}=\frac{N(\Delta y>0)-N(\Delta y<0)}{N(\Delta y>0)+N(\Delta
y<0)}$ (6)

where $\Delta y$ = $y_{\ell^{+}}$ \- $y_{\ell^{-}}$. An NLO SM prediction
yields $A^{\ell\ell}_{FB}$ = 4.8$\pm$ 0.4 [11].

The first measurement of the production-level inclusive lepton asymmetry in
the lepton+jets events is found to be $A^{\ell}_{FB}$ = 9.4${}^{+3.2}_{-2.9}$
% using the full CDF data set is reported [13] as shown in Fig. 3 (top).
Following measurements of the parton-level inclusive leptonic forward-backward
asymmetry of top-quark pairs decaying into the dilepton final state using the
full CDF data set are presented [14]. The results are $A^{\ell}_{FB}$ = 7.2
$\pm$ 6.0 % and $A^{\ell\ell}_{FB}$ = 7.6 $\pm$ 8.2 % as described in Fig. 3
(bottom), both consistent with the previous determination [13] and the SM
expectations [11].

A combination of CDF $A^{\ell}_{FB}$ measurements with the dilepton and
lepton+jets events yields $A^{\ell}_{FB}$ = 9.0${}^{+2.8}_{-2.6}$, where 80 %
of the measurement weight is due to the lepton + jets result and 20 % is due
to the dilepton result. The difference in the weights is mostly due to the
larger size of the lepton + jets final state sample. This result is about two
standard deviations larger than the NLO SM calculation of $A^{\ell}_{FB}$ =
3.8 $\pm$ 0.3 % [14], but is consistent with the 7.0-7.6 % range expected
assuming unpolarized top quark production and SM top-quark decay, given the
measured value of $A^{\ell}_{FB}$ by the CDF collaboration [13]. These CDF
results of $A^{\ell}_{FB}$ measurements are compared in Fig. 4.

Figure 3:  (top) The binned asymmetry $A^{\ell}_{FB}(q_{\ell}y_{\ell})$ in a
single charged lepton and hadronic jets (lepton+jets) after correcting for
acceptance, compared to the NLO QCD prediction of POWHEG. The best fit to Eq.
(6) for each is shown as the smooth curve of the same color. The dark (light)
gray bands indicate the statistical (total) uncertainty on the fit curve to
the data. The data points for the final state with two charged leptons
(dilepton) in (bottom) are placed at the bin centroids predicted by the POWHEG
simulation. The inner bars on the data points represent the statistical
uncertainties, while the outer bars represent the total uncertainties. The
bands indicate the one standard deviation region for statistical and
statistical + systematic uncertainties.

Figure 4:  The measurements of $A^{lep}_{FB}$ with leption+jets and dilepton
final states from CDF, along with their combination for the CDF result.

##  6 Forward-backward asymmetry in $b\bar{b}$ pairs

The light axigluon models which are viable candidates to explain the forward-
backward asymmetry in $t\bar{t}$ production, predict forward-backward
asymmetries for $b\bar{b}$ and $c\bar{c}$ systems [15, 16, 17, 18]. The
$A_{FB}$ in $b\bar{b}$ pairs at large $b\bar{b}$ mass ($M_{b\bar{b}}>$150 GeV)
is first measured using jet-triggered data and jet charge to identify
$b$-quark from $\bar{b}$-quark. To characterize the posterior and describe the
measured data, we find the highest-probability-density credible intervals at
68 % and 95 % credibility for the $b\bar{b}$ asymmetry in bins of particle-
level mass as displayed in Fig. 5 (top). The observed asymmetry is consistent
with both zero and with the SM prediction as a function of $M_{b\bar{b}}$. The
Axigluon model which has its mass $M_{A}$ = 200 GeV/c2 is excluded, while the
heavier axigluon at $M_{A}$=345 GeV c2 is good agreement with the data as
shown in Fig. 5 (bottom) [19]. CDF is currently working on the measurement of
$A_{FB}$ in $b\bar{b}$ pairs at low mass as well to confirm the results.

Figure 5:  (top) Marginal posterior probability distribution of asymmetry in
each bin of particle-level $b\bar{b}$ mass. The green and yellow bands
represent the 68 % and 95 % credible intervals, respectively. (bottom) Maximum
a posteriori points for the signal asymmetry in each mass bin. The error bars
represent the 68 % credible intervals.

##  7 Conclusion

The recent measurements on the top quark pairs production and properties from
CDF are presented. Most results of the top quark properties at CDF agree with
the SM prediction, while the top quark forward-backward asymmetry and the
leptonic asymmetry in $t\bar{t}$ system show deviations from SM predictions.
The various different experimental checks are performed for the top quark
$A_{FB}$ which is disagreement with the NLO SM calculation with both
electroweak and NLO QCD effects. We found the top quark asymmetry grows faster
than the expectation in the NLO calculation as a function of $t\bar{t}$ mass
and difference of top-antitop rapidities. CDF has been preparing to combine
these $A_{FB}$ results with D0 measurements for the Tevatron combination of
the $A_{FB}$ measurement in $t\bar{t}$ production. The first measurement of
the bottom quark forward-backward asymmetry at large $b\bar{b}$ mass is
reported. The result shows the consistency with both zero and with the SM
predictions and slightly reduces the allowed parameter space for light
axigluon models to explain the top quark forward-backward asymmetry.

##  8 Acknowledgments

I would like to thank my colleagues of the CDF collaboration for their
continuous valuable efforts and dedication to provide these important physics
results. I also thank the organizers of ICHEP 2014 in Valencia for a
productive and successful conference.

## References

  * [1] F. Abe et al. (CDF Collaboration), Phys. Rev. Lett. 74, 2626 (1995).
  * [2] S. Abachi et al. (D0 Collaboration), Phys. Rev. Lett. 74, 2632 (1995).
  * [3] J. Gao, C. S. Li, and H. X. Zhu, Phys. Rev. Lett. 110, 042001 (2013).
  * [4] T. Aaltonen et al. (CDF Collaboration), Phys. Rev. Lett. 111, 202001 (2013).
  * [5] J. Beringer et al. (Particle Data Group), Phys. Rev. D 86, 010001 (2012).
  * [6] T. Aaltonen et al. (CDF Collaboration), Phys. Rev. Lett. 112, 221801 (2014).
  * [7] T. Aaltonen et al. (CDF Collaboration), Phys. Rev. D 88, 091103(R) (2013).
  * [8] M. Czakon, P. Fiedler, and A. Mitov, Phys. Rev. Lett. 110, 252004 (2013).
  * [9] T. Aaltonen et al. (CDF and D0 Collaborations), Phys. Rev. D 89, 072001 (2014).
  * [10] T. Aaltonen et al. (CDF Collaboration), Phys. Rev. Lett. 111, 182002 (2013).
  * [11] Werner Bernreuther and Zong-Guo Si, Phys. Rev. D 86, 034026 (2012).
  * [12] T. Aaltonen et al. (CDF Collaboration), Phys. Rev. D 87, 092002 (2013).
  * [13] T. Aaltonen et al. (CDF Collaboration), Phys. Rev. D 88, 072003 (2013).
  * [14] T. Aaltonen et al. (CDF Collaboration), Phys. Rev. Lett. 113, 042001 (2014).
  * [15] G. M. Tavares and M. Schmaltz, Phys. Rev. D 84, 054008 (2011).
  * [16] C. Gross, G. M. Tavares, M. Schmaltz, and C. Spethmann, Phys. Rev. D 87, 014004 (2013).
  * [17] Seyda Ipek, Phys. Rev. D 87, 116010 (2013).
  * [18] A. Carmona, M. Chala, A. Falkowski, S. Khatibi, M. M. Najafabadi, G., J. Santiago, J. High Energy Phys.07(2014)005.
  * [19] CDF Collaboration, CDF Note Report No. 11092, 2014, http://
www-cdf.fnal.gov/physics/new/top/2014/Afbbpub/index.html.


# Kinetic data-driven approach to turbulence subgrid modeling

G. Ortali  Department of Applied Physics and Science Education, Eindhoven
University of Technology, 5600 MB Eindhoven, The Netherlands  SISSA
(International School for Advanced Studies), Trieste, Italy  A. Gabbana
Department of Applied Physics and Science Education, Eindhoven University of
Technology, 5600 MB Eindhoven, The Netherlands  Eindhoven Artificial
Intelligence Systems Institute, Eindhoven University of Technology, 5600 MB
Eindhoven, The Netherlands  N. Demo  SISSA (International School for Advanced
Studies), Trieste, Italy  G. Rozza  SISSA (International School for Advanced
Studies), Trieste, Italy  F. Toschi  Department of Applied Physics and Science
Education, Eindhoven University of Technology, 5600 MB Eindhoven, The
Netherlands  Eindhoven Artificial Intelligence Systems Institute, Eindhoven
University of Technology, 5600 MB Eindhoven, The Netherlands  CNR-IAC, I-00185
Rome, Italy

###### Abstract

Numerical simulations of turbulent flows are well known to pose extreme
computational challenges due to the huge number of dynamical degrees of
freedom required to correctly describe the complex multi-scale statistical
correlations of the velocity. On the other hand, kinetic mesoscale approaches
based on the Boltzmann equation, have the potential to describe a broad range
of flows, stretching well beyond the special case of gases close to
equilibrium, which results in the ordinary Navier-Stokes dynamics. Here we
demonstrate that, by properly tuning, a kinetic approach can statistically
reproduce the quantitative dynamics of the larger scales in turbulence,
thereby providing an alternative, computationally efficient and physically
rooted approach towards subgrid scale (SGS) modeling in turbulence. More
specifically we show that by leveraging on data from fully resolved Direct
Numerical Simulation (DNS) data we can learn a collision operator for the
discretized Boltzmann equation solver (the lattice Boltzmann method), which
effectively implies a turbulence subgrid closure model. The mesoscopic nature
of our formulation makes the learning problem fully local in both space and
time, leading to reduced computational costs and enhanced generalization
capabilities. We show that the model offers superior performance compared to
traditional methods, such as the Smagorinsky model, being less dissipative
and, therefore, being able to more closely capture the intermittency of
higher-order velocity correlations.

Introduction.— Studying turbulent flows by means of fully resolved numerical
simulations is known to pose outstanding computational challenges. The number
of dynamical degrees of freedom, whose dynamics needs to be accurately
resolved, is typically huge as it rapidly grows as the $9/4$ power of the
Reynolds number ($\rm{Re}$) frisch-book-1995 ; pope-book-2001 ; benzi-pr-2023
. Despite the ever increasing availability of computing resources, Direct
Numerical Simulations (DNS) of turbulent flows are still far out of reach for
most real-world applications, strongly motivating a continuous effort towards
the development of accurate and computationally cheaper reduced order models.
Alongside Reynolds-Averaged Navier-Stokes (RANS) models alfonsi-amr-2009 ,
where the Navier-Stokes equations are averaged over time (separating the flow
into mean and fluctuating components) large-eddy simulations (LES) zhiyin-
cja-2015 are among the most popular choices. In LES only a portion of the
dynamical degrees of freedom, associated to the larger scales, are directly
resolved, whereas the effect of the smaller scales on the large scales is
parametrized by a subgrid-scale (SGS) model. In general, considering the
filtered version of the macroscopic fields of interest, $\overline{\bm{u}}$
and $\overline{p}$ (respectively, velocity and pressure), one can write down
the Navier-Stokes equations filtered to describe only scales larger than
$\ell$:

$\frac{\partial\overline{\bm{u}}}{\partial
t}+\nabla\cdot(\overline{\bm{u}}~{}\overline{\bm{u}}+\overline{\bm{T}})=-\nabla\overline{p}+\nu\nabla^{2}\overline{\bm{u}}.$
(1)

In this filtered equation, we have assumed the density $\rho=1$, $\nu$ is the
kinematic viscosity of the fluid and
$\overline{\bm{T}}=\overline{\bm{u}\bm{u}}-\overline{\bm{u}}~{}\overline{\bm{u}}$
is the Reynolds stress tensor, accounting for the SGS fluctuations, from which
one can define a local subgrid energy flux borue-jfm-1998 . The problem of
defining an SGS closure consists in modelling the unknown term
$\overline{\bm{T}}$ as a function of only the resolved velocity field
$\overline{\bm{u}}$. Over the years, several approaches to SGS turbulence
modeling have been developed, relying on specific assumptions to approximate
the effects of the unresolved scales. The first SGS model was proposed by
Smagorinsky smagorinsky-mwr-1963 and, in view of its simplicity and stability,
it is one of the most widespread adopted approaches to estimate the action of
smaller scales on the larger scales by means of an effective eddy viscosity,
defined on the basis of the local derivatives of the resolved velocity field:

$\overline{\bm{T}}^{\rm smag}=2\nu_{\rm
eff}|\overline{\bm{S}}|-\overline{p}\bm{I},\quad\nu_{\rm
eff}=\nu+(C_{s}\Delta)^{2}|\overline{\bm{S}}|,$ (2)

where
$|\overline{\bm{S}}|=\sqrt{2\overline{S}_{ij}\overline{S}_{ij}},\overline{S}_{ij}=1/2(\partial_{i}\overline{u}_{j}+\partial_{j}\overline{u}_{i})$
represents the magnitude of the filtered strain-rate, $\Delta$ is the filter
width, and $C_{s}$ is a non-dimensional coefficient called the Smagorinsky
constant. The Smagorinsky model can be shown to be a direct consequence of the
Refined Kolmogorov Similarity Hypothesis toschi-prl-2000 and has been
generalized to non homogeneous flows leveque-jfm-2007 . Other approaches, such
as the Spalart-Allmaras (S-A) model spalart-aiaa-1992 and the $k-\epsilon$ and
$k-\omega$ models wilcox-book-1998 , resolve the eddy viscosity in terms of
additional variables, supplementing the equations for mass and momentum with
additional equations. Despite their widespread use, these models often face
limitations in accurately representing the (statistical) properties of
turbulent flows, particularly in scenarios with strong inhomogeneity and
anisotropy.

\begin{overpic}[width=429.28616pt]{slices2.pdf} \put(0.0,30.0){(a)}
\put(32.0,30.0){(b)} \put(70.0,30.0){(c)}
\put(37.0,29.5){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
{FILTERED DNS }}
\put(37.0,0.5){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
{NEURAL LBM}} \put(61.0,29.0){\normalsize$\rm{cg}=2$}
\put(92.0,29.0){\normalsize$\rm{cg}=4$} \end{overpic} Figure 1: Snapshots of
the vorticity magnitude ($|\omega|$) from 3D simulations of HIT at
$\rm{Re}\approx 6000$. In panel a) the upper half of the domain is taken from
filtered DNS data with $\rm{cg}=2$, while the lower half is obtained from a
simulation using NBLM. We show that at a qualitative level the structures
generated by NLBM closely resemble those from filtered DNS. This is further
highlighted in panel b) and c) where we show 2D slices of $|\omega|$ (cf. red
box in panel a) ), respectively at coarse graining factor $\rm{cg}=2$ and
$\rm{cg}=4$.

Recent advances in Machine Learning have opened up new perspectives for
employing Artificial Neural Networks (ANN) goodfellow-book-2016 to enhance
computational fluid dynamic solvers kochkov-pnas-2021 ; vinuesa-ncs-2022 ;
lino-prsa-2023 , and for the development of data-driven turbulence modelling
ling-jfm-2016 ; duraisamy-arfm-2019 ; wang-sigkdd-2020 . Specifically to the
LES context, several attempts have been reported in the literature wang-
pof-2018 ; maulik-jfm-2019 ; xie-prf-2020 ; park-jfm-2021 ; frezat-prf-2021 ;
tian-pnas-2023 to establish SGS closure models from extensive datasets of
fully resolved turbulent flows, leveraging on the capability of ANN to handle
high-dimensional and statistically complex data. In general, two main
approaches have been adopted; the first, which can be regarded as the ”black-
box” approach, consists in making no assumption on the form of the SGS terms
beck-jcp-2019 ; zhou-cf-2019 ; xie-pre-2019 ; kurz-ijhff-2023 , while in the
second one the task of the ANN consists in tuning the free-parameters of an
already established SGS model sarghini-cf-2003 ; xie-pof-2019 ; wang-aip-2021
.

In this Letter we take a fresh look at the problem of establishing a SGS
closure model with ANN, and introduce for the first time, to the best of our
knowledge, a data-driven kinetic-based approach to turbulence modeling. We
employ Physics-Informed Machine Learning (PIML) to enhance the capabilities of
the Lattice Boltzmann Method (LBM) succi-book-2018 , and exploit the extra
degrees of freedom provided by the mesoscopic description to learn a new
collision operator which effectively acts as SGS model. Our framework relies
on the pathway connecting kinetic theory to hydrodynamics, offering the
possibility to learn from data macroscopic equations that extend beyond the
Navier-Stokes level. Thanks to the locality of the collision operator, at no
point this scenario involves the evaluation of spatial derivatives of
macroscopic fields chen-sc-2003 , drastically simplifying the training process
and, in turn, the interpretability of the model. Remarkably, and at variance
with respect to other fully local SGS models, our formalism shows the
potential for capturing the inverse transfer of energy from small to large
scales.

\begin{overpic}[width=429.28616pt]{diagram2.pdf}
\put(19.0,44.0){\large$L^{3}$}
\put(18.0,9.0){\large$\left(\frac{L}{\rm{cg}}\right)^{3}$}
\put(1.5,67.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
{DNS }}
\put(1.5,34.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}
{NLBM }} \put(53.0,64.0){repeat $\times$ cg times}
\put(32.0,38.3){\large\color[rgb]{0.1328125,0.546875,0.1328125}\definecolor[named]{pgfstrokecolor}{rgb}{0.1328125,0.546875,0.1328125}
Filter by $\rm cg$}
\put(78.0,38.3){\large\color[rgb]{0.1328125,0.546875,0.1328125}\definecolor[named]{pgfstrokecolor}{rgb}{0.1328125,0.546875,0.1328125}
Filter by $\rm cg$} \put(59.0,56.0){\huge$\dots$}
\put(28.0,46.5){\normalsize${f}^{\rm pre}_{\rm DNS}(t)$}
\put(38.0,46.5){\normalsize${f}^{\rm post}_{\rm DNS}(t)$}
\put(47.0,46.5){\normalsize${f}^{\rm pre}_{\rm DNS}(t+\Delta t)$}
\put(69.0,46.5){\normalsize$\dots$} \put(80.0,46.5){\normalsize$\dots$}
\put(86.0,46.5){\normalsize${f}^{\rm pre}_{\rm DNS}(t+\rm{cg}\Delta t)$}
\put(32.5,49.5){\tiny COLLIDE} \put(44.0,49.5){\tiny STREAM }
\put(72.0,49.5){\tiny COLLIDE} \put(84.0,49.5){\tiny STREAM }
\put(28.0,16.0){\normalsize${f}^{\rm pre}_{\rm cg}(t)$}
\put(47.0,16.0){\normalsize$\tilde{f}^{\rm post}_{\rm cg}(t)$}
\put(70.0,16.0){\normalsize${f}^{\rm post}_{\rm cg}(t)$}
\put(87.0,16.0){\normalsize${f}^{\rm pre}_{\rm cg}(t+\Delta t)$}
\put(37.0,13.2){\normalsize$\Omega_{\rm NN}$} \put(79.0,13.2){\tiny ANTI-
STREAM}
\put(54.0,4.0){\Large\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}
$\mathcal{L}(\tilde{f}^{\rm post}_{\rm cg},f^{\rm post}_{\rm cg})$}
\end{overpic} Figure 2: Schematic representation of the training process for
the turbulence SGS model. The upper panel corresponds to the DNS simulation on
a $L^{3}$ grid, while in lower panel the Neural LBM (NLBM) operates on a
$(L/\rm{cg})^{3}$ grid, with coarse-graining factor $\rm{cg}$. The mapping
between DNS into coarse-grained data is given by the application of a filter:
the pre-collision state at a generic time step $t$ at the coarse grained level
( $f_{\rm cg}^{\rm pre}(t)$ ) is obtained by filtering the DNS pre-collision
state ( $f_{\rm DNS}^{pre}(t)$ ). Similarly, the post-collision data at the
coarse-grained level ( $f_{\rm cg}^{\rm post}(t)$ ) is obtained by first
filtering the post-collision DNS state at time step $t+\rm{cg}\Delta t$, and
then by applying the inverse of the streaming operator. Following this
procedure it is possible to create an arbitrarily large dataset for training
an ANN to which we assign the task of minimizing the mismatch between
$\tilde{f}_{\rm cg}^{\rm post}(t)=\Omega_{\rm NN}(f_{\rm cg}^{\rm pre}(t))$
and $f_{\rm cg}^{\rm post}(t)$ under a given error-metric $\mathcal{L}$.

Methods.— The Lattice Boltzmann Method (LBM) benzi-pr-1992 ; succi-book-2018 ;
kruger-book-2017 has emerged in the past decades as a popular solver for
computational fluid dynamics. At variance with methods that explicitly
discretize the Navies-Stokes equations, LBM operates at the mesoscopic level,
providing the description of a fluid system in terms of a (small) set of
particle distribution functions (populations) whose dynamic is governed by the
discrete Boltzmann equation:

$f_{i}(\bm{x}+\bm{c}_{i}\Delta t,t+\Delta
t)=f_{i}(\bm{x},t)+\Omega_{i}(\bm{x},t),$ (3)

where at each grid node $\bm{x}$, lattice populations $f_{i}(\bm{x},t)$ are
defined along the discrete components of a velocity stencil
$\\{\bm{c}_{i}\\},i=1,\dots,Q$. A popular choice for the collision operator
$\Omega$ is the single-time relaxation BGK model bhatnagar-pr-1954 :

$\Omega_{i}(\bm{x},t)=-\frac{\Delta
t}{\tau}\left(f_{i}(\bm{x},t)-f_{i}^{\rm{eq}}(\bm{x},t)\right),$ (4)

which models collisions as a relaxation towards an equilibrium distribution,
with $\tau$ the relaxation rate. The macroscopic variables of interest can be
obtained as the lower-order moments of the lattice populations. Moreover, it
can be shown via an asymptotic analysis that Eq. 3 yields a second order
approximation of the Navier-Stokes equations kruger-book-2017 .

It is an expedient for the description of our method to split the time
evolution of Eq. 3 into two steps, following the stream and collide paradigm:

$\displaystyle f_{i}^{\rm post}(\bm{x},t)$ $\displaystyle=f_{i}^{\rm
pre}(\bm{x},t)+\Omega_{i}(\bm{x},t),$ (5) $\displaystyle f_{i}^{\rm
pre}(\bm{x},t+\Delta t)$ $\displaystyle=f^{\rm
post}_{i}(\bm{x}-\bm{c}_{i}\Delta t,t),$ (6)

where here and in what follows we will denote with $f_{i}^{\rm pre}$
($f_{i}^{\rm post}$) the pre(post)-collision populations.

Following the framework introduced in Ref. corbetta-epje-2023 , we replace the
collision operator with a ANN:

$\tilde{f}_{i}^{\rm post}(\bm{x},t)=f_{i}^{\rm pre}(\bm{x},t)+\Omega^{\rm
NN}\left(f_{i}^{\rm pre}(\bm{x},t)\right).$ (7)

The resulting algorithm, to which we refer as Neural Lattice Boltzmann Method
(NLBM), employs a physics-constrained ANN to establish a data-driven SGS
model. Our approach leverages on two key ingredients of the LBM algorithm:

  1. 1.

The LBM mesoscopic representation makes use of a larger number of degrees of
freedom (i.e. the number of discrete lattice populations) than the macroscopic
observable of interest. This observation opens up the possibility of using ANN
to encode extra information in the model.

  2. 2.

The non-linear terms encountered at the NS level of description (cf. Eq. 1)
are fully embedded in the LBM via the collision operator and thus purely local
(in LBM “non-linearity is local, non-locality is linear” succi-book-2018 ).
This observation drastically reduces the cost of training ANN, since it offers
the possibility to restrict the input and output of the network to local
quantities without explicit evaluation of gradients of macroscopic fields.

We summarize here the main steps required to train the ANN, as sketched in
Fig. 2, leaving full technical details to the Supplementary Material. We
consider simulations of Homogeneous Isotropic Turbulence (HIT), fully resolved
on a $L^{3}$ domain using the standard LBM formulation (Eq. 3) with BGK
collision operator (Eq. 4). We regard this as the target ground truth DNS data
(upper panel in Fig. 2). Next, we define a coarse-graining factor $\rm{cg}$,
and define a filter which projects, at a given time step $t$, data from the
DNS grid to a coarse-grained grid of size $(L/\rm{cg})^{3}$, where the SGS
model will operate (lower panel in Fig. 2). In order to define an arbitrarily
large training dataset consisting of pre and post collision data at the
coarse-grained level, respectively $f_{\rm cg}^{\rm pre}(t)$ and $f_{\rm
cg}^{\rm post}(t)$, we take the following steps: We start from the pre-
collision populations at the DNS level, $f_{\rm DNS}^{\rm pre}(t)$, and
immediately apply the filter in order to obtain $f_{\rm cg}^{\rm pre}(t)$.
Next, we advance the DNS simulation for $\rm{cg}\Delta t$ time-steps, yielding
the pre-collision (and post-streaming) populations $f_{\rm DNS}^{\rm
pre}(t+\rm{cg}\Delta t)$. Since $\Delta t$ steps at the coarse-grained level
correspond to $\rm{cg}\Delta t$ steps at the DNS level, by filtering $f_{\rm
DNS}^{\rm pre}(t+\rm{cg}\Delta t)$ we obtain $f_{\rm cg}^{\rm pre}(t+\Delta
t)$. Finally, we obtain $f_{\rm cg}^{\rm post}(t)$ by reversing the streaming
operation, i.e. we anti-stream populations on the coarse grid with respect to
the corresponding velocity component $-\bm{c}_{i}$.

At this stage, we can train a ANN which, under a given error-metric
$\mathcal{L}$, minimizes the mismatch between $f_{\rm cg}^{\rm post}(t)$ and
the prediction of the network taking $f_{\rm cg}^{\rm pre}(t)$ as input:

$\mathcal{L}\left(\Omega^{\rm NN}(f_{\rm cg}^{\rm pre}(t)),f_{\rm cg}^{\rm
post}(t)\right).$ (8)

We have observed that the two main ingredients allowing for training models
which deliver accurate and stable results throughout an entire simulation are
i) imposing hard-constraints on the preservation of mass and momentum
corbetta-epje-2023 and ii) using unrolled training brandstetter-arxiv-2022 to
compute the loss over consecutive timesteps. Full details on the ANN
architecture and the training process are provided as Supplementary
Information.

Figure 3:  Energy spectrum for simulations of HIT at $\rm{Re}\approx 6000$.
The results from DNS (blue curve) are compared with NLBM (red) and Smagorinsky
(green). For the two SGS models we report the average spectrum from $100$
simulations starting from different initial conditions. The shaded curves
corresponds to one standard deviation from the average value.

Numerical results.— We consider numerical simulations of HIT, with for
reference a DNS on a $L=128^{3}$ grid at Reynolds number $\rm{Re}\approx
6000$. The parameters are selected in such a way that the plain LBM algorithm
would encounter numerical instabilities when working on coarse-grained grids
with $\rm{cg}=\\{2,4\\}$. Fig. 1 qualitatively summarizes our findings.
Starting from the same initial configuration, we present 3D and 2D
representation of the absolute value of the vorticity
$|\bm{\omega}|=|\nabla\times\bm{u}|$ at a late stage of the simulation,
comparing results from the filtered DNS with NLBM. Our model provides stable
simulations, with flow patterns virtually indistinguishable from those of the
filtered DNS.

Figure 4: Structure functions (cf. Eq. 9) of order $p$, ranging between $p=1$
to $p=6$, versus $S^{3}$, with in blue data from filtered DNS with
$\rm{cg}=4$, in red data from simulations using NLBM, and in green data from
simulations using the Smagorinsky model. The inset shows the deviation for the
scaling exponents $\xi_{p}$ from the K41 scaling $\frac{p}{3}$.

On a more quantitative ground, we now turn to the evaluation of the
statistical properties of the turbulent flows produced by NLBM simulations.
For comparison, we hereafter also take into consideration LBM simulations
equipped with the Smagorinsky SGS (Eq. 2), which allows for stable simulations
of the coarse domains. In Fig. 3 we present the kinetic energy spectrum
$E(k)$. We observe that NLBM well compares with the results produced by
Smagorinsky, with slightly superior results at $\rm{cg}=4$, where we can
observe a spectrum less dissipative than the Smagorinsky one, and closer to
the scaling of DNS data. Next, we study the scaling behaviour of high-order
Eulerian Structure function, which for a generic order $p$ can be computed
from numerical data as

$S^{p}(l)=\left\langle\left[\left(\bm{u}(\bm{x}+\bm{l})-\bm{u}(\bm{x})\right)\cdot{\hat{\bm{l}}}\right]^{p}\right\rangle.$
(9)

We use the extended self similarity (ESS) benzi-pre-1993 to determine the
scaling exponents $\xi_{p}$. In Fig. 4 we plot the structure functions of
order $p$, ranging between $p=1$ to $p=6$, versus $S^{3}(l)$, comparing data
from filtered DNS with $\rm{cg}=4$, and data from simulations using NLBM and
the Smagorinsky model. In the inset of Fig. 4 we report the deviation of
$\xi_{p}$ from the K41 scaling ${p}/{3}$. The results highlight that, due to
its lower dissipation at small scales, NLBM provides anomalous scaling within
error bars from filtered DNS. On the other hand, the discrepancies observed in
the Smagorinsky model are due to the fact that under these parameters the
inertial range is shrinked (see the supplementary material for a more detailed
analysis).

A major advantage of our kinetic-based approach lies in the possibility of
having a physical interpretation of the action of the ANN, by projecting the
lattice populations to the velocity moment space (see Supplementary Material).
We have observed that NLBM extends the single-relaxation BGK collision
operator, from where DNS data was generated, to a multi-relaxation collision
operator lallemand-pre-2000 ; dhumieres-ptsa-2002 . The model introduces a
non-linear relation for pre and post collision values of the moments related
to the bulk-viscosity, something which is often used to increase the stability
of numerical simulations. The model preserves a linear-dependency for the pre
and post collision values of the moments related to the kinetic viscosity.
This allows to fit the effective viscosity from the numerical data, which can
be used to establish a direct comparison with the Smagorinsky model (Eq. 2).
In Fig. 5 we compare the probability distribution function (PDF) of the value
of the Smagorinsky constant $C^{2}$ fitted from NLBM data. The average value
of the Smagorinsky constant for NLBM ($\langle C^{2}_{\rm NLBM}\rangle\approx
0.11$) is about a factor two smaller than the one used in simulations with the
Smagorinsky SGS model ($C^{2}=0.2$). Remarkably, in Fig. 5 we observe a tail
taking negative values, suggesting that our model occasionally displays an
inverse transfer of energy from small to large scales, a feature completely
lacking to the Smagorinski model which is, by its own nature, fully
dissipative.

Figure 5:  Probability distribution function (PDF) of the fitted value of the
Smagorinsky constant $C^{2}$ from NLBM data. The inset highlights that the
average value $<C^{2}_{\rm NLBM}>\approx 0.11$ is about a factor two smaller
than the one used in simulations with the Smagorinsky SGS model ($C^{2}=0.2$).
The presence of a tail with negative values highlights the fact that in NLBM
it is possible to capture the inverse transfer of energy from small to large
scales.

Conclusion.— Summarizing, we have introduced a novel kinetic-based approach to
SGS modeling, combining LBM with physics informed ANN. Our model allows for
stable simulations on coarse domains, offering the possibility of reducing the
computational costs of DNS, in turn preserving the statistical properties of
turbulent flows under HIT settings. The model compares well with Smagorinsky,
and our results highlight a better agreement with DNS in terms of energy
spectra and estimation of the anomalous scaling exponents. Moreover, we have
shown that the model in principle supports the possibility to describe the
inverse transfer of energy from small to large scales, which we regard as
promising particularly in vision of future application to more involved
numerical setups. To conclude, this work opens up the possibility of
exploiting the extra degrees of freedom of the mesoscopic representation in
the LBM to develop novel SGS models. In future works, we plan to extend our
analysis to flows subject to anisotropy as well as wall-bounded flows. A
further intriguing question concerns the application of this framework beyond
SGS. Work to explore employing our kinetic data-driven approach in other
contexts is currently underway.

Acknowledgments.— We wish to thank Roberto Benzi for useful discussions. This
work was partially funded by the Dutch Research Council (NWO) through the
UNRAVEL project (with Project No. OCENW.GROOT.2019.044).

## References

  * (1) U. Frisch, Turbulence: The Legacy of A.N. Kolmogorov, (Cambridge University Press, 1995).
  * (2) S. B. Pope, Turbulent Flows, (Cambridge University Press, 2000).
  * (3) R. Benzi and F. Toschi, Physics Reports 1021, 1–106 (2023).
  * (4) G. Alfonsi, Applied Mechanics Reviews 62, 040802 (2009).
  * (5) Y. Zhiyin, Chinese Journal of Aeronautics 28, 11-24 (2015).
  * (6) V. Borue and S. A. Orszag, Journal of Fluid Mechanics 366, 1–31 (1998).
  * (7) J. Smagorinsky, Monthly weather review 91, 99–164 (1963).
  * (8) F. Toschi, E. Lévêque, and G. Ruiz-Chavarria, Phys. Rev. Lett. 85, 1436–1439 (2000).
  * (9) E. Lévêque, F. Toschi, L. Shao, and J. P. Bertoglio, Journal of Fluid Mechanics 570, 491–502 (2007).
  * (10) P. Spalart and S. Allmaras In 30th aerospace sciences meeting and exhibit, AIAA Paper 1992-0439 (1992).
  * (11) D. C. Wilcox, Turbulence Modeling for CFD (DCW Industries, 2006).
  * (12) I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning (MIT Press, 2016).
  * (13) D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer, Proceedings of the National Academy of Sciences 118, e2101784118 (2021).
  * (14) R. Vinuesa and S. L. Brunton, Nature Computational Science 2, 358–366 (2022).
  * (15) M. Lino, S. Fotiadis, A. A. Bharath, and C. D. Cantwell, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 479, 20230058 (2023).
  * (16) J. Ling, A. Kurzawski, and J. Templeton,  Journal of Fluid Mechanics 807, 155–166 (2016).
  * (17) K. Duraisamy, G. Iaccarino, and H. Xiao,  Annual review of fluid mechanics 51, 357–377 (2019).
  * (18) R. Wang, K. Kashinath, M. Mustafa, A. Albert, and R. Yu,  Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 1457–1466 (2020).
  * (19) Z. Wang, K. Luo, D. Li, J. Tan, and J. Fan,  Physics of Fluids 30, 125101 (2018).
  * (20) R. Maulik, O. San, A. Rasheed, and P. Vedula, Journal of Fluid Mechanics 858, 122-144 (2019).
  * (21) C. Xie, J. Wang, and W. E, Phys. Rev. Fluids 5, 054606 (2020).
  * (22) J. Park and H. Choi, Journal of Fluid Mechanics 914, A16 (2021).
  * (23) H. Frezat, G. Balarac, J. Le Sommer, R. Fablet, and R. Lguensat, Phys. Rev. Fluids 6, 024607 (2021).
  * (24) Y. Tian, M. Woodward, M. Stepanov, C. Fryer, C. Hyett, D. Livescu, and M. Chertkov, Proceedings of the National Academy of Sciences 120, e2213638120 (2023).
  * (25) A. Beck, D. Flad, and C.-D. Munz, Journal of Computational Physics 398, 108910 (2019).
  * (26) Z. Zhou, G. He, S. Wang, and G. Jin, Computers & Fluids 195, 104319 (2019).
  * (27) C. Xie, J. Wang, K. Li, and C. Ma, Phys. Rev. E 99, 053113 (2019).
  * (28) M. Kurz, P. Offenhäuser, and A. Beck, International Journal of Heat and Fluid Flow 99, 109094 (2023).
  * (29) F. Sarghini, G. de Felice, and S. Santini, Computers & Fluids 32, 97-108 (2003).
  * (30) C. Xie, J. Wang, H. Li, M. Wan, and S. Chen, Physics of Fluids 31, 085112 (2019).
  * (31) Y. Wang, Z. Yuan, C. Xie, and J. Wang, AIP Advances 11, (2021).
  * (32) H. Chen, S. Kandasamy, S. Orszag, R. Shock, S. Succi, and V. Yakhot, Science 301, 633–636 (2003).
  * (33) R. Benzi, S. Succi, and M. Vergassola,  Physics Reports 222, 145–197 (1992).
  * (34) S. Succi, The Lattice Boltzmann Equation: For Complex States of Flowing Matter (Oxford Scholarship Online, Oxford, 2018).
  * (35) T. Krueger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. Viggen, The Lattice Boltzmann Method: Principles and Practice (Springer, 2016).
  * (36) P.L. Bhatnagar, E.P. Gross, M.Krook, Phys. Rev. 94, 511 (1954).
  * (37) A. Corbetta, A. Gabbana, V. Gyrya, D. Livescu, J. Prins, and F. Toschi, The European Physical Journal E 46, 10 (2023).
  * (38) J. Brandstetter, D. E. Worrall, and M. Welling, arXiv:2202.03376, (2022).
  * (39) R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, F. Massaioli, and S. Succi, Phys. Rev. E 48, R29–R32 (1993).
  * (40) P. Lallemand and L.-S. Luo, Phys. Rev. E 61, 6546–6562 (2000).
  * (41) D. D’Humières, I. Ginzburg, M. Krafczyk, P. Lallemand, and L.-S. Luo, Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 360, 437-51 (2002).

Supplementary Information for
“Kinetic data-driven approach to turbulence subgrid modeling”

##  I Numerical Method

In this section, we provide a basic introduction to the Lattice Boltzmann
Method (LBM). The reader is referred to Ref. sm-kruger-book-2017 ; sm-succi-
book-2018 for a thorough introduction to the topic. The LBM is a class of
numerical fluid-dynamics solvers. At variance with conventional methods that
explicitly discretize the Navies-Stokes equations, LBM takes root from the
mesoscopic Boltzmann equation:

$\frac{\partial f}{\partial t}+\xi_{\alpha}\frac{\partial f}{\partial
x_{\alpha}}+F_{\alpha}\frac{\partial f}{\partial\xi_{\alpha}}=\Omega(f)$ (S1)

where $f(\bm{x},\bm{\xi},t)$ is the particle distribution function,
representing the average number of particles in a small element of phase-space
centered at position $\bm{x}$ with velocity $\bm{\xi}$ at time $t$, and
$F_{\alpha}$ is the sum of external forces acting on the system. The right
hand side of the equation, $\Omega(f)$, is the collision operator, describing
the changes in $f$ due to particle collisions.

Figure S1: D3Q19 Velocity stencil

Index | $\bm{c}_{i}$ | $w_{i}$
---|---|---
$0$ | $(0,0,0)$ | $1/3$
$1$ | $(-1,0,0)$ | $1/18$
$2$ | $(0,-1,0)$ | $1/18$
$3$ | $(0,0,-1)$ | $1/18$
$4$ | $(0,0,1)$ | $1/18$
$5$ | $(0,1,0)$ | $1/18$
$6$ | $(1,0,0)$ | $1/18$
$7$ | $(-1,-1,0)$ | $1/36$
$8$ | $(-1,0,-1)$ | $1/36$
$9$ | $(-1,0,1)$ | $1/36$
$10$ | $(-1,1,0)$ | $1/36$
$11$ | $(0,-1,-1)$ | $1/36$
$12$ | $(0,-1,1)$ | $1/36$
$13$ | $(0,1,-1)$ | $1/36$
$14$ | $(0,1,1)$ | $1/36$
$15$ | $(1,-1,0)$ | $1/36$
$16$ | $(1,0,-1)$ | $1/36$
$17$ | $(1,0,1)$ | $1/36$
$18$ | $(1,1,0)$ | $1/36$

Table S1: D3Q19 Velocity stencil - discrete velocity vectors $\bm{c}_{i}$ and
corresponding weights $w_{i}$.

The lattice time-discrete counterpart of Eq. S1 reads as

$f_{i}(\bm{x}+\bm{c}_{i}\Delta t,t+\Delta
t)=f_{i}(\bm{x},t)+\Omega_{i}(\bm{x},t)+f_{i}^{\rm ext},$ (S2)

In the above equation, known as lattice Boltzmann equation, the discrete
particle distribution functions (populations) $f_{i}(\bm{x},t)$ move at the
next time step $t+\Delta t$ with velocity $\bm{c}_{i}$ to a neighboring grid
cell $\bm{x}+\bm{c}_{i}\Delta t$. Moreover, $f_{i}^{\rm ext}$ is the discrete
version of the external force term. The discretization of the velocity space
is a key element in the definition of LBM; it consists in replacing the
continuum velocity space with a small set of discrete velocities
$\mathcal{V}=\\{\bm{c}_{i}\in\mathbb{R}^{3}\\}$ corresponding to the abscissa
of a Gauss-Hermite quadrature (with corresponding weights $w_{i}$ ) that
ensure that all moments of the distribution function up to a prescribed order
are correctly recovered sm-philippi-pre-2006 ; sm-shan-jcs-2016 . In this
work, we use the D3Q19 velocity stencil (Fig. S1), defined in Tab. S1, which
allows to express mass and momentum as weighted sums of $f_{i}(\bm{x},t)$:

$\rho=\sum_{i}f_{i},\quad\rho\bm{u}=\sum_{i}f_{i}\bm{c}_{i}.$ (S3)

There are several possible choices to approximate the collision operator
$\Omega_{i}$. Arguably the simplest one that can be used for Navier-Stokes
simulations is the Bhatnagar-Gross-Krook (BGK) sm-bhatnagar-pr-1954 operator:

$\Omega_{i}(f)=\frac{\Delta t}{\tau}(f_{i}^{\rm eq}-f_{i}),$ (S4)

which relaxes the populations at rate $\tau$ towards a local equilibrium
defined by the Maxwell-Boltzmann distribution, for which we use a second order
polynomial expansion:

$f^{\rm{eq}}_{i}(\rho,\bm{u})=\,w_{i}\rho\left(1+\frac{\bm{u}\cdot\bm{c}_{i}}{c_{s}^{2}}+\frac{(\bm{u}\cdot\bm{c}_{i})^{2}-(c_{s}|\bm{u}|)^{2}}{2c_{s}^{4}}\right),$

with $c_{s}=1/\sqrt{3}$ the sound speed in the lattice, and $w_{i}$ a lattice-
dependent set of weighting factors (see Tab. S1).

Several options are available for implementing the body forces in LBM sm-
kruger-book-2017 . In this work we adopt the Exact Difference Method (EDM)
scheme sm-kupershtokh-cma-2009 , where the external forcing term is computed
as

$f_{i}^{\rm ext}=f_{i}^{\rm eq,\rm shift}-f_{i}^{\rm eq},$ (S5)

where the shifted equilibrium $f_{i}^{\rm eq,\rm shift}$ is computed from Eq.
I based on the shifted velocity $\bm{u}^{\text{shift}}=\bm{u}+\frac{F}{\rho}$.

Following an asymptotic analysis, such as the Chapman-Enskog expansion sm-
chapman-book-1952 , it can be shown that Eq. S2 combined with Eq. S4 delivers
in the low-Mach number limit a second order approximation of the Navier-Stokes
equations, with the following expression putting in relationship the
relaxation time parameter $\tau$ with the kinematic viscosity $\nu$ of the
fluid:

$\nu=\left(\tau-\frac{\Delta t}{2}\right)c_{s}^{2}\quad.$ (S6)

###  I.1 Implementing the Smagorinsky Turbulence Model in LBM

In the main text we have reported numerical results of LBM simulations
equipped with the Smagorinsky SGS model. The model involves adding an eddy
viscosity term to the governing equations to account for unresolved turbulent
fluctuations at the grid scale. This eddy viscosity is computed as:

$\nu_{\rm eff}=\nu+(C_{s}\Delta)^{2}|\overline{\bm{S}}|,$ (S7)

where
$|\overline{\bm{S}}|=\sqrt{2\overline{S}_{ij}\overline{S}_{ij}},\overline{S}_{ij}=1/2(\partial_{i}\overline{u}_{j}+\partial_{j}\overline{u}_{i})$
is the magnitude of the filtered strain-rate, $\Delta$ is the filter width,
and $C_{s}$ is a non-dimensional coefficient called the Smagorinsky constant.

In order to incorporate the Smagorinsky model in a LBM scheme one needs to i)
compute the strain rate tensor in each cell to then ii) compute the eddy
viscosity $\nu_{\rm eff}$ to finally iii) locally adjust the relaxation time
to match $\nu_{\rm eff}$. In LBM the strain rate tensor can be computed
locally from the non-equilibrium part of the distribution function:

$\overline{S}_{ij}=-\frac{3\tau}{2}\Pi_{ij}^{\rm neq},\quad\Pi_{ij}^{\rm
neq}=\sum_{q}c_{qi}c_{qj}\left(f_{q}-f_{q}^{\rm eq}\right)\quad.$ (S8)

By combining Eq. S6 with Eq. S7 it is then possible to compute the local
relaxation time due to the Smagorinsky model:

$\tau_{\rm eff}=\frac{\nu_{\rm eff}}{c_{s}^{2}}+\frac{\Delta t}{2}.$ (S9)

##  II Artificial Neural Network Architecture and Training Process

###  II.1 Definition of the Training Dataset

The ground truth data used to train the ANN consists of N pairs of 19-tuples

$\left\\{\left({f_{i,k}^{\rm pre}},{f_{i,k}^{\rm
post}}\right),k=1,2,\dots,N\right\\}.$ (S10)

The data is obtained from Direct Numerical Simulation (DNS) of Homogeneous
Isotropic Turbulence. We conduct LBM simulations in a three-dimensional cubic
domain of side $L=128$, with periodic boundary conditions and $\Delta t=1$.
Simulations are driven by a stationary external force, ensuring a divergence-
free flow, defined for the $x$, $y$ and $z$ components as:

$\begin{cases}F_{x}(y)=F\sin\left(\frac{2\pi}{L}y\right)\\\
F_{y}(z)=F\sin\left(\frac{2\pi}{L}z\right)\\\
F_{z}(x)=F\sin\left(\frac{2\pi}{L}x\right)\end{cases}.$ (S11)

We work at a value of the Reynolds number $\rm{Re}\approx 6000$, achieved
tuning the relaxation time $\tau$ and the magnitude of the external force as
detailed in Tab. S2. We collect data after an initial transient phase of
$T_{0}$ timesteps, ensuring that the dynamics reaches a statistically
stationary state. The full state of the system (in terms of lattice
populations) is sampled every $\delta_{T}$ timesteps, for a total of $N_{T}$
timesteps. Next, we define a coarse graining factor $\rm{cg}$ and a filter
$\phi$, which allows mapping at a generic time step $t$ lattice populations
from the $L^{3}$ DNS grid to the $(L/\rm{cg})^{3}$ coarse-grained grid:

$f^{\rm cg}_{i}(t)\leftarrow\phi(f^{\rm DNS}_{i}(t)).$ (S12)

Several possible choices can be made for the filter $\phi$, which may
substantially impact the quality of the results. In this work we consider a
sub-sampling strategy, and define the coarse-grained grid by sampling one
every $\rm{cg}$ lattice points along the different coordinates. We leave as a
future work a careful evaluation and comparison of different choices for
$\phi$.

At this stage, we can establish the training set (Eq. S10) at the coarse
grained level, following the procedure sketched in Fig. 2 of the main text,
which we here summarize:

  1. 1.

We start from the pre-collision populations at the DNS level, $f_{\rm
DNS}^{\rm pre}(t)$, and apply a filter in order to obtain $f_{\rm cg}^{\rm
pre}(t)$.

  2. 2.

We advance the DNS simulation for $\rm{cg}\Delta t$ time-steps obtaining the
post-streaming state, i.e. $f_{\rm DNS}^{\rm pre}(t+\rm{cg}\Delta t)$.

  3. 3.

By filtering the DNS level we obtain $f_{\rm cg}^{\rm pre}(t+\Delta t)$.

  4. 4.

In order to obtain the pre-streaming state, i.e. the post-collision values at
time $t$, we apply the inverse of the streaming operation at the coarse
grained level by reversing the streaming direction, i.e. with respect to
$-\bm{c}_{i}$.

For the results reported in the main text we have used $\rm{cg}=2$ and
$\rm{cg}=4$, resulting in effective lattice sizes of, respectively,
$L^{\prime}=64$ and $L^{\prime}=32$. Note that to achieve the same $\rm{Re}$
on the coarse grids one could use the same numerical parameters provided in
Tab. S2, alongside a rescaling of the relaxation time parameter

$\tau^{\prime}=\frac{1}{\rm
cg}\left(\tau-\frac{1}{2}\right)+\frac{1}{2}\quad.$ (S13)

However, we shall remark that the choice of the parameters is such that plain
LBM simulations at the coarse grid level would lead to numerical
instabilities.

Parameter |  Description |  Value
---|---|---
$\rm{Re}$ |  Reynolds Number |  $6\times 10^{3}$
$\tau$ |  Relaxation time Parameter |  $0.5032$
$F$ |  Force term Amplitude |  $5\times 10^{-6}$
$L$ |  Lattice Size for DNS |  $128$
$\Delta_{x}$ |  Grid size |  $1$
$\Delta_{t}$ |  Timestep |  $1$
$T_{0}$ |  Transient Phase Duration |  $5\times 10^{4}$
$\delta_{T}$ |  Sampling Interval |  $10^{3}$
$N_{T}$ |  Total Number of Timesteps |  $10^{3}$
$M$ |  Lattice Points Sampled per Timestep |  $512$
$\lambda$ |  Learning rate |  $1\times 10^{-4}$
$\omega$ |  Weight decay rate |  $1\times 10^{-6}$
Table S2:  The table reports the numerical parameters used in DNS numerical
simulations (top), for the definition of the training dataset (middle) and for
the hyper-parameters used by the ADAM optimizer to train the ANN.

###  II.2 Neural Network Architecture and Optimization

In this section we detail our ANN model, which consists in learning a
correction term adding to the BGK collision operator. We define the post-
collision state as

$\tilde{f}_{i}^{\rm post}=f_{i}^{\rm pre}+\Omega(f_{i}^{\rm
pre})+\frac{1}{\tau}\Omega^{\rm NN}(f_{i}^{\rm pre})$ (S14)

where $\Omega$ is the BGK collision operator in Eq. S4, while $\Omega^{\rm
NN}$ is the correction term due to the ANN.

We consider an ANN consisting of a Multi-layer Perceptron (MLP) with input and
output layers of size $19$ (equal to the number of lattice populations sitting
on one single grid point), and three fully connected hidden layers of size
$[300,300,100]$, in combination with LeakyReLU activation function. Following
Ref. sm-corbetta-epje-2023 , before the output layer we also include a non-
trainable layer posing hard constraints on the conservation of mass and
momentum. In our case we require that the correction term given by the ANN
does not introduce additional mass and momentum in the output:

$\displaystyle\sum_{i}\Omega^{\rm NN}(f^{\rm pre}_{i})$ $\displaystyle=0,$
(S15) $\displaystyle\sum_{i}\Omega^{\rm NN}(f^{\rm pre}_{i})\bm{c}_{i}$
$\displaystyle=\bm{0}.$ (S16)

The above equations are satisfied using

$\displaystyle\hat{f}_{i}$ $\displaystyle=\Omega^{\rm NN}(f^{\rm
pre}_{i})+\kappa_{1}+\kappa_{2}c_{i,x}+\kappa_{3}c_{i,y}+\kappa_{4}c_{i,z},$
(S17)

with the following coefficients specific to the D3Q19 model:

$\displaystyle\kappa_{1}$ $\displaystyle=-\frac{1}{19}\sum_{i}\Omega^{\rm
NN}(f^{\rm pre}_{i}),$ $\displaystyle\kappa_{2}$
$\displaystyle=-\frac{1}{10}\sum_{i}\Omega^{\rm NN}(f^{\rm pre}_{i})c_{i,x},$
$\displaystyle\kappa_{3}$ $\displaystyle=-\frac{1}{10}\sum_{i}\Omega^{\rm
NN}(f^{\rm pre}_{i})c_{i,y},$ $\displaystyle\kappa_{4}$
$\displaystyle=-\frac{1}{10}\sum_{i}\Omega^{\rm NN}(f^{\rm pre}_{i})c_{i,z}.$

We have observed that the model defined in Eq. S14, albeit performing well in
one-step predictions, once plugged in a simulation typically leads to
numerical instabilities after just few iterations. This behavior is commonly
observed in autoregressive models, i.e. models that feed in the output as an
input for the next step sm-brandstetter-arxiv-2022 . This is due to
accumulations of errors that eventually deviate the dynamics of the system
towards configurations not covered in the training dataset, a phenomenon known
in literature as distribution shift sm-goodfellow-book-2016 . In the next
section we discuss a procedure which allows to mitigate this problem.

###  II.3 Unrolled training

In this section we discuss the Unrolled training technique, which we employ to
enhance the stability of the model. In an auto-regressive model, i.e. a model
that recursively uses its output as an input for the successive timesteps, one
must ensure the stability of the model with respect to the accumulation of
errors. This is a requirement, for example, for any time integration scheme,
where the stability hinges on both the system’s dynamics and the
discretization scheme employed. Moreover, in the domain of Machine Learning,
another phenomenon may arise, since the accumulation of errors may lead the
dynamics of the system into unexplored regions of the optimization space,
deviating from the training conditions. This problem is typically framed in
the context of distribution shift, and deals with scenarios in which the data
distribution during training differs from that encountered during testing.
Fields such as language modeling and robotics commonly encounter such
challenges. For instance, maintaining coherence in generated text over
multiple iterations is a complex task for language models sm-holtzman-
arxiv-2019 .

Several solutions have been proposed to deal with this problem. One
possibility is to modify the architecture of the Neural Network, for example
forcing the spectra of the linear matrices, stabilizing the dynamics. Examples
of this approach are Antisymmetric Neural Networks sm-goles-dam-1986 . Another
possibility is to use regularization strategies such as noise injection in the
input and/or the output of the Neural Network sm-labach-arxiv-2019 .

In this work, we employ a technique known in literature as unrolled training
or trajectory learning sm-brandstetter-arxiv-2022 . In short, it consists in
applying the model to make predictions over $L_{T}$ consecutive timesteps at
training time, backpropagating the gradients through some or all such
timesteps. This is similar to the Backpropagation Through Time method, used
for training Recurrent Neural Networks, with the main difference being that in
unrolled training no memory term is present in the Neural Network.

We train the model using the following error metric:

$\mathcal{L}=\sum_{x_{i}\in S}\left(\tilde{f}^{\rm
post}(x_{i},t_{0}+L_{t})-f^{\rm post}(x_{i},t_{0}+L_{t})\right)^{2},$ (S18)

where for each training epoch $S$ represents a randomly selected subset of the
lattice points. One single training epoch consists of the following steps:

  1. 1.

We select a snapshot from the ground truth dataset, i.e. the pre-collision
populations on the entire filtered grid of size $(L/\rm{cg})^{3}$ at a fixed
time step. This serves as the initial condition.

  2. 2.

We randomly select a subset of lattice points, denoted as
$S=\\{x_{i}\\}_{i=0}^{batch\\_size}$, which will be used compute the loss.

  3. 3.

We evolve the system for $L_{T}$ timesteps, alternating the collision step via
Eq. S14, and the streaming step.

  4. 4.

We compute the loss via Eq. S18 comparing against ground-truth values at time
$t+L_{T}$.

  5. 5.

Finally, we backpropagate the loss through all the timesteps to compute the
gradients and update the weights of the Neural Network.

We remark that, while the full domain is evolved using the ANN, the loss term
is computed only on randomly selected subset of points. The reason for this is
twofold. First, it allows reducing computational costs, both in terms of
memory requirements as well as overall training time. Second, it allows to
increase stochasticity in the computation of gradients for updating the
weights of the network.

This strategy introduces a few extra hyperparameters in the training of the
network. The results reported in the main text have been obtained using a
batch size $|S|$ of $512$ elements, with a trajectory lenght of $L_{T}=12$.
The choice of these parameters impacts the time and memory requirements of the
training process (longer trajectories and bigger batch sizes require more time
and memory), and on the resulting stability of the model (too short
trajectories give less stable models). The optimization is done using Adam
Optimizer, training for a total of $2\times 10^{4}$ epochs with early stopping
based on test loss. The optimization parameters are reported in Table S2.

##  III Structure Functions and Extended Self Similarity

In Fig. 4 in the main text we have reported results with of high order
statistical properties of turbulence, analyzing data from simulations using
our data-driven SGS model. In this section we provide further details; for the
convenience of the reader we state once again the definition of the Eulerian
Structure function of order $p$:

$S^{p}(l)=\left\langle\left[\left(\bm{u}(\bm{x}+\bm{l})-\bm{u}(\bm{x})\right)\cdot{\hat{\bm{l}}}\right]^{p}\right\rangle.$
(S19)

\begin{overpic}[width=429.28616pt]{structure_functions_extended.pdf}
\put(0.0,63.0){(a)} \put(34.0,63.0){(b)} \put(68.0,63.0){(c)}
\put(0.0,30.5){(d)} \put(34.0,30.5){(e)} \put(68.0,30.5){(f)}
\put(25.0,61.0){$\rm{cg}=2$} \put(59.0,61.0){$\rm{cg}=2$}
\put(93.0,61.0){$\rm{cg}=2$} \put(25.0,28.5){$\rm{cg}=4$}
\put(59.0,28.5){$\rm{cg}=4$} \put(93.0,28.5){$\rm{cg}=4$} \end{overpic} Figure
S2: Results for structure functions (cf. Eq. 9) of order $p$, ranging between
$p=1$ to $p=6$, with in blue data from filtered DNS, in red data from
simulations using NLBM, and in green data from simulations using the
Smagorinsky model. On the top row (panels (a)-(c)) results for $cg=2$, on the
bottom row (panels (d)-(f)) results for $cg=4$. In the first column we report
the Structure function $S^{p}(l)$ versus $l$, in the second column $S^{p}$
versus $S^{3}$, and in the third column the logarithmic derivative
$\frac{\partial log(S^{p})}{\partial log(S^{3})}$.

In Fig. S2, the first row refers to results for a grid with coarse graining
factor $\rm{cg}=2$, whereas the second row shows results for $\rm{cg}=4$. For
all cases we show in blue color filtered DNS data, in red our NLBM model, and
for comparison in green results using Smagorinsky. In the first column we plot
the structure functions $S^{p}$ vs $l$, with $p$ ranging between $1$ to $6$.
The central column shows the results of the extended self similarity analysis,
where this time we plot the structure functions $S^{p}$ versus versus $S^{3}$.
The insets shows the deviation of the scaling exponents $\xi_{p}$ from K41
scaling $p/3$. The scaling exponents have been computed fitting the slopes,
which are better highlighted in the third column, where we show the
logarithmic derivative $\frac{\partial log(S^{p})}{\partial log(S^{3})}$. The
fitting ranges are defined considering intervals where the logarithmic
derivative stays approximately constant, avoiding the dissipative interval by
excluding the first point, resulting in $S^{3}$ values of $[4,20]$ for $cg=2$
and $[8,40]$ for $cg=4$.

The results here presented strengthen the message and the considerations made
in the main text, showcasing that NLBM better captures the kinetic energy at
the large scales with respect than Smagorinsky, which exhibits instead larger
dissipation.

##  IV Generalization

A primary concern in data-driven modelling is over-fitting and poor
generalization goodfellow-book-2016 , connected to the challenge of generating
a comprehensive, high-quality training dataset that can cover all potential
uses. Additionally, ANN models often lack interpretability and physical
consistency, rendering them unreliable in application contexts. In this
section, we discuss the generalization capabilities of our model.

We consider the ANN trained with the parameters from Tab. S2, with
$\rm{cg}=4$, and consider a few scenarios slightly departing from the training
conditions. In Fig. S3 we report the energy spectrum for four different cases,
comparing the results of DNS simulations, NLBM and Smagorinsky, changing the
following parameters. Panels from (a) to (d) cover the following set of
parameters:

  1. (a)

Different magnitude of the forcing term (for reference, value used during
training: $F=5\times 10^{-6}$, $Re\approx 6\times 10^{3}$)

     * •

$F=1\times 10^{-5}$ ( $Re\approx 1.2\times 10^{4}$) ,

     * •

$F=1\times 10^{-6}$ ( $Re\approx 1.2\times 10^{3}$) .

  2. (b)

A non-homogeneous forcing term:

$\begin{cases}F_{x}(y)=2F\sin\left(\frac{2\pi}{L}y\right)\\\ F_{y}(z)=0\\\
F_{z}(x)=0\end{cases}.$

  3. (c)

A higher wave number in the sine function:

$\begin{cases}F_{x}(y)=F\sin\left(2\frac{2\pi}{L}y\right)\\\
F_{y}(z)=F\sin\left(2\frac{2\pi}{L}z\right)\\\
F_{z}(x)=F\sin\left(2\frac{2\pi}{L}x\right)\end{cases}.$

  4. (d)

Higher and lower values of the kinematic viscosity, changing the relaxation
time (for reference, value used in training $\tau=0.5032$, $Re\approx 6\times
10^{3}$)

  * •

$\tau=0.5048$ ( $Re\approx 4\times 10^{3}$) ,

  * •

$\tau=0.5024$ ( $Re\approx 8\times 10^{3}$) .

The results show that the NLBM yields stable simulations for all the cases
covered, considering both higher and lower value for $\rm{Re}$ with respect to
the one considered for training the ANN. NLBM provides results systematically
in better agreement with DNS with respect to Smagorinsky, with the exception
of one case in Fig. S3(a) where the absolute value of the external force was
taken to be 5 times smaller than the one used during training. We attribute
this to a possible overfitting of the small scales, which, as can be seen in
the plot, overlap for NLBM regardless of the magnitude of the external force.

We should also remark that the generalization capabilities of the model do not
extend to different grid sizes and different coarse graining factor; such
cases currently require the training of a new model from scratch.

\begin{overpic}[width=429.28616pt]{energy_spectrum_generalization.pdf}
\put(0.0,61.0){(a)} \put(0.0,30.0){(c)} \put(52.0,30.0){(d)}
\put(52.0,61.0){(b)} \end{overpic} Figure S3:  Energy spectrum for simulations
of HIT with different parameters to evaluate the capability of NLBM to
generalize outside of the training dataset. Blue curves represent DNS results,
red for NLBM and green for the Smagorinsky SGS. In panel a) we present results
increasing or decreasing the magnitude of the forcing term, giving
respectively $Re\approx 1.2\times 10^{4}$ and $Re\approx 1.2\times 10^{3}$ (vs
$Re\approx 6\times 10^{3}$ used at training time). In panel b) we show results
using a non-homogeneous forcing term. In panel c) we show the effect of
forcing with a higher wave number. In panel d) we show the effect of
increasing and decreasing the relaxation time parameter (hence the kinematic
viscosity).

##  V Physical Interpretation of model action

In this section we provide a physical interpretation of the action of the ANN.
While our model works at the kinetic level, it is simple to map the lattice
populations $f$ to the moments space, and observe the action of the model in
terms of physical quantities. In order to outline the procedure, we start from
Eq. S2 with the single-relaxation time BGK operator in Eq. S4, introduce an
invertible matrix $\bm{M}$, and recast the equation in the following form:

$\begin{split}f^{\rm post}-f^{\rm
pre}=\bm{M}^{-1}\frac{1}{\tau}\bm{M}\left(f^{\rm eq}-f^{\rm pre}\right)+f^{\rm
ext}\\\ =\bm{M}^{-1}\bm{S}\bm{M}\left(f^{\rm eq}-f^{\rm pre}\right)+f^{\rm
ext}\end{split},$ (S20)

where $\bm{M}$ defines a transformation from lattice populations $f(\bm{x},t)$
to macroscopic moments $\bm{m}=\bm{M}f(\bm{x},t)$, and the relaxation matrix
$\bm{S}=\text{diag}\left(\frac{1}{\tau},\dots,\frac{1}{\tau}\right)$ acts on
macroscopic moments space.

A generalization of the BGK collision operator in given by the Multi-
Relaxation Time (MRT) collision operator sm-dhumieres-ptsa-2002 , which allows
for individual relaxation rates for the different macroscopic moments.

\begin{overpic}[width=216.81pt]{{moments_dist_13.pdf}} \put(2.0,87.0){(a)}
\end{overpic}\begin{overpic}[width=212.47617pt]{moments_dist_1.pdf}
\put(2.0,90.0){(b)} \end{overpic}

Figure S4:  Scatterplots of $\bm{m^{\rm neq}}$ versus $\bm{m^{\rm diff}}$ (see
text for their definition) from NLBM data for the components corresponding to
the component $p_{xy}$ of the shear stress tensor in (a) and to the kinetic
energy $e$ in (b). The data shown here comes from a model trained using a
coarse graining factor $cg=4$. In both panels, the red dotted line correspond
to the value used for the BGK model at training time $\tau^{\prime}=0.5008$
(rescaled with respect to $\rm{cg}$ as per Eq. S13). A clear linear dependence
can be observed for case $p_{xy}$, which allows to define an effective
relaxation time $\tau_{\rm eff}\approx 0.5023$, slightly larger than
$\tau^{\prime}$, hence leading to a larger effective value of the kinetic
viscosity. In panel b) we see instead strong non-linear effects which does not
allow for fitting a value of the bulk viscosity.

In this framework, several possible choices can be operated for the matrix
$\bm{M}$, corresponding to a map of the lattice populations into different
sets of macroscopic moments $\bm{m}$, in which the individual relaxation times
$\tau_{1}$, …, $\tau_{19}$ can have different physical interpretations. In
what follows we consider

$\bm{M}=\setcounter{MaxMatrixCols}{19}\begin{pmatrix}1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1\\\
-30&-11&-11&-11&-11&-11&-11&8&8&8&8&8&8&8&8&8&8&8&8\\\
12&-4&-4&-4&-4&-4&-4&1&1&1&1&1&1&1&1&1&1&1&1\\\
0&-1&0&0&0&0&1&-1&-1&-1&-1&0&0&0&0&1&1&1&1\\\
0&4&0&0&0&0&-4&-1&-1&-1&-1&0&0&0&0&1&1&1&1\\\
0&0&-1&0&0&1&0&-1&0&0&1&-1&-1&1&1&-1&0&0&1\\\
0&0&4&0&0&-4&0&-1&0&0&1&-1&-1&1&1&-1&0&0&1\\\
0&0&0&-1&1&0&0&0&-1&1&0&-1&1&-1&1&0&-1&1&0\\\
0&0&0&4&-4&0&0&0&-1&1&0&-1&1&-1&1&0&-1&1&0\\\
0&2&-1&-1&-1&-1&2&1&1&1&1&-2&-2&-2&-2&1&1&1&1\\\
0&-4&2&2&2&2&-4&1&1&1&1&-2&-2&-2&-2&1&1&1&1\\\
0&0&1&-1&-1&1&0&1&-1&-1&1&0&0&0&0&1&-1&-1&1\\\
0&0&-2&2&2&-2&0&1&-1&-1&1&0&0&0&0&1&-1&-1&1\\\
0&0&0&0&0&0&0&1&0&0&-1&0&0&0&0&-1&0&0&1\\\
0&0&0&0&0&0&0&0&0&0&0&1&-1&-1&1&0&0&0&0\\\
0&0&0&0&0&0&0&0&1&-1&0&0&0&0&0&0&-1&1&0\\\
0&0&0&0&0&0&0&-1&1&1&-1&0&0&0&0&1&-1&-1&1\\\
0&0&0&0&0&0&0&1&0&0&-1&-1&-1&1&1&1&0&0&-1\\\
0&0&0&0&0&0&0&0&-1&1&0&1&-1&1&-1&0&-1&1&0\\\ \end{pmatrix},$ (S21)

which maps to the following macroscopic moments:

$\bm{m}=\setcounter{MaxMatrixCols}{19}\begin{pmatrix}\rho,&e,&\epsilon,&j_{x},&q_{x},&j_{y},&q_{y},&j_{z},&q_{z},&p_{xx},&\pi_{xx},&p_{ww},&\pi_{ww},&p_{xy},&p_{yz},&p_{xz},&m_{x},&m_{y},&m_{z}\end{pmatrix},$
(S22)

where:

  * •

$\rho$, $j_{x}$, $j_{y}$ and $j_{z}$ are the mass and the components of
momentum along $x,y,z$.

  * •

$e$ is the kinetic energy;

  * •

$\epsilon$ is the kinetic energy squared;

  * •

$q_{x}$, $q_{y}$ and $q_{z}$ are the heat fluxes along $x,y,z$.

  * •

$p_{xx}$, $p_{ww}$, $p_{xy}$, $p_{yz}$ and $p_{xz}$ are the components of the
symmetric traceless viscous stress tensor;

  * •

$\pi_{xx}$, $\pi_{ww}$, $m_{x}$, $m_{y}$ and $m_{z}$ correspond to higher
order moments, with no obvious physical interpretation (in the LBM jargon they
go under the name of ghost modes).

In turn, the macroscopic moments are associated to the following relaxation
parameters:

$\bm{S}=\rm{diag}\setcounter{MaxMatrixCols}{19}\begin{pmatrix}0,&\tau_{e},&\tau_{\epsilon},&0,&\tau_{q},&0,&\tau_{q},&0,&\tau_{q},&\tau_{\nu},&\tau_{\pi},&\tau_{\nu},&\tau_{\pi},&\tau_{\nu},&\tau_{\nu},&\tau_{\nu},&\tau_{m},&\tau_{m},&\tau_{m}\end{pmatrix},$
(S23)

where the relaxation times associated to conserved quantities have been set to
zero.

\begin{overpic}[width=429.28616pt]{c2_est_sm.pdf} \put(1.0,25.0){(a)}
\put(34.0,25.0){(b)} \put(68.0,25.0){(c)} \end{overpic} Figure S5:
Probability distribution function (PDF) of the fitted value of the Smagorinsky
constant $C^{2}$ from NLBM data, for coarse graining factor $cg=4$. In panel
a) data fitted considering the $p_{xy}$ component of the shear stress tensor,
in b) using $p_{yz}$ and in c) using $p_{xx}$. The inset highlights that the
average value, marked with vertical bars, comparing the results of NLBM (red)
against the value used in the Smagorinsky SGS model in green ($C^{2}=0.2$).

Likewise for the BGK, also in MRT sm-dhumieres-ptsa-2002 it is possible to
establish a connection between the relaxation time parameters and the
macroscopic tranport coefficients. For example, the kinematic kinematic
viscosity $\nu$ and bulk viscosity $\zeta$ of the model are given by:

$\nu=c_{s}^{2}\left(\tau_{\nu}-\frac{1}{2}\right);\>\>\>\zeta=\frac{5-9c_{s}^{2}}{9}\left(\tau_{e}-\frac{1}{2}\right),$
(S24)

While for NLBM an asymptotic analysis for defining the transport coefficients
is not viable from an analytic point of view, it is still possible to
numerically evaluate how the ANN relaxes the different moments $\bm{m}$. To
this aim we multiply the LHS and the RHS of Eq. S20 by $\bm{M}$:

$\begin{split}\underbrace{\bm{M}\left(f^{\rm post}-f^{\rm pre}-f^{\rm
ext}\right)}_{\bm{m^{\rm diff}}}=\bm{S}\underbrace{\bm{M}\left(f^{\rm
eq}-f^{\rm pre}\right)}_{\bm{m^{\rm neq}}}\end{split},$ (S25)

and define $\bm{m^{\rm diff}}$ and $\bm{m^{\rm neq}}$, two column vectors
defined in the moment space similarly to $\bm{m}$.

From numerical data, we can plot the different components of $\bm{m^{\rm
diff}}$ vs $\bm{m^{\rm neq}}$. In Fig. S4 we report and example of this
analysis, for the case $\rm{cg}=4$, where the individual components are
centered with respect to their average value. We remark that in a MRT
collision operator we would observe a linear dependence, with a slope of the
regression line $1/\tau_{i}$ depending on the choice of the relaxation
parameters in $\bm{S}$. For NLBM we observe, to good approximation, a linear
relationship for the moments related to the shear stress tensor $p_{ij}$. In
Fig. S4(a) we provide an example for the component $p_{xy}$. From a linear fit
we can obtain an effective value for the relaxation time $\tau_{\rm eff}$. At
variance with Eq. S23, we observe two different relaxation times, one for the
components $p_{xx}$ and $p_{ww}$ ($\tau_{\rm eff}\approx 0.5125$), and a
second one for $p_{xy}$, $p_{yz}$ and $p_{xz}$, ($\tau_{\rm eff}\approx
0.5023$). Both values are larger than the corresponding value of the
relaxation time rescaled as for Eq. S13, $\tau=0.5008$, thus implying a larger
effective shear viscosity.

For all the other moments we generally observe a non linear dependence. In
Fig. S4(b) we show an example for the kinetic energy coefficient which relates
to the bulk viscosity. It is common practice in literature to adjust the bulk
viscosity to enhance the stability of numerical methods, and recently an ANN-
based approach has been reported also for LBM sm-horstmann-cf-2024 .

To further test the importance of non-linearities introduced by the ANN we
have tested a MRT collision operator in which the matrix $\bm{S}$ has been
obtained performing a linear fit of $\tau_{\rm eff}$ for all the different
moments independently (hence discarding the non-linear contributions):

$\bm{S}=\rm{diag}\setcounter{MaxMatrixCols}{19}\begin{pmatrix}0,&\tau_{1},&\tau_{2},&0,&\tau_{3},&0,&\tau_{4},&0,&\tau_{5},&\tau_{6},&\tau_{7},&\tau_{8},&\tau_{9},&\tau_{10},&\tau_{11},&\tau_{12},&\tau_{13},&\tau_{14},&\tau_{15},\end{pmatrix}$
(S26)

with

$\begin{split}\tau_{1}&=0.1554,\quad\tau_{2}=0.0947,\quad\tau_{3}=0.5220,\quad\tau_{4}=0.6921,\quad\tau_{5}=0.5614,\quad\\\
\tau_{6}&=0.5122,\quad\tau_{7}=0.5564,\quad\tau_{8}=0.5121,\quad\tau_{9}=0.5483,\quad\tau_{10}=0.5020,\quad\\\
\tau_{11}&=0.5028,\quad\tau_{12}=0.5023,\quad\tau_{13}=0.4597,\quad\tau_{14}=0.4910,\tau_{15}=0.4872\end{split}$

With this choice of $\bm{S}$ numerical simulations became unstable after just
few iterations, similarly to the BGK case. We leave as a future work further
investigation on the role of the heat flux and of the ghost modes. We focus
here instead on the effective viscosity, which can be computed from the shear
stress tensor akin to the Smagorinsky turbulence model. From the linear fit of
the components of the shear stress tensor it is possible to compute $\nu_{\rm
eff}$, which in combination with Eq. S7 allows to define what would be the
equivalent of the Smagorinsky constant for NLBM. In Fig. S5 we report the PDF
of the fitted value of the Smagorinsky constant from NLBM data, with in panel
(a) and (b) the results for $p_{xy}$ and $p_{yz}$ respectively, and in $(c)$
for $p_{xx}$. The average values are comparable with Smagorinsky for cases (a)
and (b), about a factor two smaller than the one used in simulations with the
Smagorinsky SGS model ($C^{2}=0.2$), with the largest discrepancies observed
for case (c). For all cases, the presence of a tail with negative values
highlights the fact that in NLBM it is possible to capture the inverse
transfer of energy from small to large scales.

## References

  * (1) T. Krueger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. Viggen, The Lattice Boltzmann Method: Principles and Practice (Springer, 2016).
  * (2) S. Succi, The Lattice Boltzmann Equation: For Complex States of Flowing Matter (Oxford Scholarship Online, Oxford, 2018).
  * (3) P.C. Philippi, L.A. Hegele, L.O.E. dos Santos, R. Surmas,  Phys. Rev. E 73, 056702 (2006).
  * (4) X. Shan,  Journal of Computational Science 17, 475–481 (2016).
  * (5) P.L. Bhatnagar, E.P. Gross, M.Krook, Phys. Rev. 94, 511 (1954).
  * (6) A. L. Kupershtokh, D.A. Medvedev, and D. Karpov, Computers & Mathematics with Applications 58(5), 965–974 (2009).
  * (7) S. Chapman and T.G. Cowling, The Mathematical Theory of Non-Uniform Gases (American Journal of Physics, 1952).
  * (8) A. Corbetta, A. Gabbana, V. Gyrya, D. Livescu, J. Prins, and F. Toschi, The European Physical Journal E 46, 10 (2023).
  * (9) J. Brandstetter, D. E. Worrall, and M. Welling, arXiv:2202.03376, (2022).
  * (10) I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning (MIT Press, 2016).
  * (11) A. Holtzman, J. Buys, L. Du, M. Forbes, and Y. Choi,  arXiv:1904.09751, (2019).
  * (12) E. Goles, Discrete Applied Mathematics 13, 97-100 (1986).
  * (13) A. Labach, H. Salehinejad, and S. Valaee, arXiv:1904.13310, (2019).
  * (14) D. D’Humières, I. Ginzburg, M. Krafczyk, P. Lallemand, and L.-S. Luo, Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 360, 437-51 (2002).
  * (15) J. T. Horstmann, M. C. Bedrunka, and H. Foysi, Computers & Fluids 272, 106191 (2024).


# Full LST-1 data reconstruction with the use of convolutional neural networks

Jakub Juryšek1, Etienne Lyard1 and Roland Walter1 for the CTA LST project

###### Abstract

The Cherenkov Telescope Array (CTA) will be the world’s largest and most
sensitive ground-based gamma-ray observatory in the energy range from a few
tens of GeV to tens of TeV. The LST-1 prototype, currently in its
commissioning phase, is the first of the four largest CTA telescopes, that
will be built in the northern site of CTA in La Palma, Canary Islands, Spain.
In this contribution, we present a full-image reconstruction method using a
modified InceptionV3 deep convolutional neural network applied on non-
parametrized shower images. We evaluate the performance of optimized networks
on Monte Carlo simulations of LST-1 shower images, and compare the results
with the performance of the standard reconstruction method. We also show how
both methods work on real-data reconstruction.

1University of Geneva, Astronomy Department, Geneva, Switzerland;
<EMAIL_ADDRESS>
##  1 Introduction

When primary gamma-ray photons or cosmic-ray particles enter the atmosphere,
they interact with atomic nuclei and produce showers of secondary particles
emitting Cherenkov radiation. In ground-based gamma-ray astronomy, this
radiation can be observed by Imaging Atmospheric Cherenkov Telescopes (IACTs),
forming roughly elliptical images of the shower in the camera plane of a
focusing telescope. In typical IACT observations, the trigger rate from
diffuse cosmic-ray background is about 1000x higher than the trigger rate from
gamma-ray photons even for the brightest gamma-ray sources in the sky, and
therefore a strong background suppression is necessary. Besides reducing this
background, the goal of the image analysis is to reconstruct the energy and
incoming direction of the primary gamma-ray. In this study, we compare the
reconstruction performance of standard Random Forests (RFs) (Breiman 2001)
with a novel method using a deep convolutional neural network (CNN) on Monte
Carlo (MC) simulations (Heck et al. 1998; Bernlöhr 2008) and first data from
the LST-1 telescope prototype for CTA (CTA Consortium 2019).

##  2 Methods

We used a modified InceptionV3 CNN architecture (Fig. 1) (Szegedy et al. 2016;
Lyard et al. 2020) consisting of a variable number of inception modules
followed by five sequential convolutional layers, whose capability of
gamma/hadron separation was previously proven by Lyard et al. (2020). In this
contribution, we study the performance of this network for full event
reconstruction. Our networks were trained on MC gamma and MC proton shower
images from a standard CTA MC production data set (Prod5) for LST-1. We used
adadelta and rmsprop optimizers for gamma/hadron separation and
energy/direction reconstruction, respectively. Unlike RFs, which need image
parameters on input leading to a loss of information, we trained CNN directly
on shower images and times of maxima of waveform in each pixel which follow a
temporal gradient, calibrated in lstchain v0.6.3 111https://github.com/cta-
observatory/cta-lstchain/releases/tag/v0.6.3. In order to get a standard
rectangular input image for the Keras Python library, hexagonal images from
the telescope had to be reshaped following the method introduced in Lyard et
al. (2020). Our results suggest that such reshaping doesn’t introduce any
position reconstruction bias if more than 3 inception modules are used. We
also applied relatively weak selection cuts on MC events (Intensity222The
total integrated signal in pixels containing Cherenkov light emitted by the
shower particles, after removal of the background pixels. > 50 p.e.,
leakage333The fraction of the shower signal contained in the two outermost
pixel rings of the camera. < 0.2), to create a MC sample that was used for
training/testing of CNNs and RFs consistently. The number of inception modules
and dropout were optimised for each reconstruction task on a subset of
$10^{5}$ training events (left Fig. 2). Parameters providing the best results
in terms of Area Under the Curve (AUC)444Area under the ROC curve, usually
being used as a measure of the ability of a classifier to distinguish between
classes., energy and angular resolution were used for final training of each
CNN on $5\times 10^{5}$ events.

Figure 1.: Architecture of modified InceptionV3 used in this study. Number of
inception modules was optimized for each reconstruction task. Figure from
Lyard et al. (2020).

Figure 2.: _Left:_ An example of angular resolution evaluated on testing point
gammas for different parameters of CNN. _Right:_ Angular resolution of
reconstruction by optimized CNN compared with Random Forests.

##  3 Reconstruction performance

Figure 3 (left) shows receiver operating characteristic curve (ROC) for
gamma/hadron separation trained on simulated diffuse gammas and diffuse
protons. A comparison with the separation power of RFs trained/tested on
Hillas parameters extracted from the same MC sample is also shown, and it can
be clearly seen that CNN significantly outperform RFs. Angular and energy
resolution curves are shown in Figure 2 (right) and Figure 3 (right),
respectively. For the mid-energy range the angular resolution improvement
gained from CNN is not significant, but it’s very prominent for low and high
energies, which can be a hint that CNN are particularly good in the
reconstruction of noise-dominated low energy shower images, or showers cut by
the camera edges. In energy reconstruction, RFs outperform CNNs by 5%. The
worse reconstruction power of CNN in this case could be caused by image
normalisation in the preprocessing step and needs to be addressed in a
following study.

Figure 3.: _Left:_ ROC for gamma/hadron separation evaluated on on-axis point
gammas and diffuse protons. _Right:_ Energy resolution of CNN reconstruction
compared with Random Forests.

##  4 Analysis of real data

Despite its good performance on simulations, CNN often fails when applied to
real data due to subtle differences between the two (e.g. different night sky
background level, or stars present in the real field of view), and thus the
evaluation of the performance on the real data is essential. We used CNN and
RF trained on the same MC events to reconstruct events from Crab Nebula
observations from November 2020 (19 mins on-source and 15 mins off-source),
and reached a similar detection significance of 6.2$\sigma$ for both methods
after carefully matching background levels using an adaptive gammaness cut on
the CNN reconstructed data (Gammaness cut for RFs processed data was 0.8), and
common $\Theta^{2}$ < 0.05 $\mathrm{deg}^{2}$ cut on the source region. A
significance map for the CNN reconstructed data is shown in Figure 4 (left).
The $\Theta^{2}$ distribution for both methods (right Fig. 4) shows that
contrary to the performance evaluated on MC, the angular resolution of the CNN
on data is worse than we expected, which is probably caused by MC-data
differences, which will be addressed in a following study.

Figure 4.: _Left:_ Significance map of events reconstructed with CNN. _Right:_
$\Theta^{2}$ distribution of events in energy range 0.01- 10 TeV.

##  5 Conclusions

In this contribution, we demonstrated a potential of CNNs to be used not only
for gamma/hadron separation, but also for full event reconstruction with a
competitive performance to the standard RFs reconstruction. Even though the
CNN performance on real data seems to be slightly worse than expected from the
full MC study, there is much room for improvement, particularly in MC-data
tuning.

### Acknowledgments

We gratefully acknowledge financial support from the agencies and
organizations listed here: www.cta-observatory.org/consortium_acknowledgments.
We would like also to thank CSCS in Lugano for the computer hours on PizDaint
computer.

## References

  * Bernlöhr (2008) Bernlöhr, K. 2008, Astroparticle Physics, 30, 149. 0808.2253
  * Breiman (2001) Breiman, L. 2001, Machine Learning, 45, 5. URL https://doi.org/10.1023/A:1010933404324
  * CTA Consortium (2019) CTA Consortium 2019, Science with the Cherenkov Telescope Array (WORLD SCIENTIFIC). https://www.worldscientific.com/doi/pdf/10.1142/10986, URL https://www.worldscientific.com/doi/abs/10.1142/10986
  * Heck et al. (1998) Heck, D., Knapp, J., Capdevielle, J. N., Schatz, G., & Thouw, T. 1998, CORSIKA: a Monte Carlo code to simulate extensive air showers. (TIB Hannover)
  * Lyard et al. (2020) Lyard, E., Walter, R., Sliusar, V., & Produit, N. 2020, in Journal of Physics Conference Series, vol. 1525 of Journal of Physics Conference Series, 012084\. 1907.02428
  * Szegedy et al. (2016) Szegedy, C., Vanhoucke, V., Ioffe, S., Shlens, J., & Wojna, Z. 2016, in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)