# Inertial range statistics of the Entropic Lattice Boltzmann in 3D turbulence

Michele Buzzicotti  Department of Physics and INFN, University of Rome Tor
Vergata, via della Ricerca Scientifica 1, 00133, Rome, Italy.  Guillaume
Tauzin  Department of Physics and INFN, University of Rome Tor Vergata, via
della Ricerca Scientifica 1, 00133, Rome, Italy.  Chair of Applied Mathematics
and Numerical Analysis, Bergische Universität Wuppertal, Gaußstrasse 20, 42119
Wuppertal, Germany

###### Abstract

We present a quantitative analysis of the inertial range statistics produced
by Entropic Lattice Boltzmann Method (ELBM) in the context of 3d homogeneous
and isotropic turbulence. ELBM is a promising mesoscopic model particularly
interesting for the study of fully developed turbulent flows because of its
intrinsic scalability and its unconditional stability. In the hydrodynamic
limit, the ELBM is equivalent to the Navier-Stokes equations with an extra
eddy viscosity term Malaspinas2008 . From this macroscopic formulation, we
have derived a new hydrodynamical model that can be implemented as a Large-
Eddy Simulation (LES) closure. This model is not positive definite, hence,
able to reproduce backscatter events of energy transferred from the subgrid to
the resolved scales. A statistical comparison of both mesoscopic and
macroscopic entropic models based on the ELBM approach is presented and
validated against fully resolved Direct Numerical Simulations (DNS). Besides,
we provide a second comparison of the ELBM with respect to the well known
Smagorinsky closure. We found that ELBM is able to extend the energy spectrum
scaling range preserving at the same time the simulation stability. Concerning
the statistics of higher order, inertial range observables, ELBM accuracy is
shown to be comparable with other approaches such as Smagorinsky model.

##  I Introduction

Turbulence is common in nature and its unpredictable behavior has fundamental
consequences on the understanding and control of various systems, from smaller
engineering devices dumitrescu2004rotational ; Davidson2015 ; Pope2001 , up to
the larger scales geophysical and astrophysical flows gill2016atm ;
alexakis2018 ; barnes2001 . Turbulent flows are described by the Navier-Stokes
equations (NSE),

$\begin{cases}\partial_{t}\bm{u}+\bm{u}\cdot\bm{\nabla}\bm{u}=-\bm{\nabla}p+\nu\nabla^{2}\bm{u}+\bm{f}\\\
\bm{\nabla}\cdot\bm{u}=0\end{cases}$ (1)

which give the evolution of the incompressible velocity field
$\bm{u}(\bm{x},t)$, with kinematic viscosity $\nu$, subject to a pressure
field $p$ and to an external forcing $\bm{f}$. However, even though the
equations of motion are known since almost two hundred years a direct
analytical approach remains elusive Frisch1995 . To overcome mathematical
difficulties, scientists, helped by the exponential growth of the
computational power, have tried to search for approximate solutions using
numerical algorithms fox2003computational ; Pope2001 ; ishihara2009study .
Unfortunately, also in this direction not every effort were successful.
Indeed, the NSE have a very rich non-linear dynamics, where a large range of
scales from the domain size up to the small scales fluctuations, are coupled
together. This results in a very high dimensional problem, with the
dimensionality proportional to the range of active scales Frisch1995 , and
with highly intermittent statistics dominated by the presence of extreme and
rare fluctuations yeung2015extreme ; benzi2010inertial ; iyer2017reynolds . As
a consequence, no matter how powerful new-supercomputers are, numerical
algorithms cannot handle all the degrees of freedom involved in the dynamics
Davidson2015 . The way out is to introduce a scale separation and compute only
the dynamics of degrees of freedom belonging to a subset of scales while
neglecting the other ones Meneveau2000 ; filippova2001multiscale ; Dong2008a ;
Dong2008b . However due to the non-linearity of NSE there is never a real
scale separation in the equations of motion bohr2005dynamical ; Frisch1995 ,
and the small-scale effects on the scales of interest need to be compensated
by the introduction of a model. In other words, the benefit of multi-scale
modeling is to achieve a scale separation, and the main challenge is to find a
“closure”, which guaranties the numerical stability being at the same time the
most accurate as possible in reproducing the coupling of the missing scales on
the resolved ones. This is the principle behind the celebrated Large-Eddy
simulations (LES), which actually solve the flow only on a subset of “large”
scales by filtering each term of the NSE and replacing with a closure the non-
linear coupling term between the resolved and the sub-grid scales (SGS)
Pope2001 ; Lesieur2005 . One of the most important differences between the
real coupling term coming from the filtered NSE and the common LES closures
used in the literature is for the latter to be purely dissipative to ensure
the simulation stability, see Meneveau2000 . As a consequence, it is
impossible for the closures to reproduce the backscatter events of energy
going from the SGS to the large scales, with important consequences on the
statistics of the resolved velocity field. Another possible numerical
approach, who has gained particular popularity, consists in solving the flow’s
macroscopic hydrodynamical properties as an approximation of its mesoscopic
behaviour de2013non . The Lattice Boltzmann Method (LBM) falls into that
category benzi1992lattice ; Succi2001 . In LBM, the flow is simulated by
evolving the Boltzmann equation for the single phase density function,
$f(\bm{x},t)$. The idea is to evolve the streaming and collision of particles
distribution functions, where the possible velocities are restricted on a
subset of discrete lattice directions frisch1986lattice ; Wolf-Gladrow2000 .
It is crucial to choose the collision operator so that macroscopically, (in
the limit of small Knudsen number), the dynamics described by the NSE is
recovered qian1995recent ; chen1998lattice . The most common collision
operator is the Bhatnagar-Gross-Krook model (BGK), see Bhatnagar1954 ,
corresponding to the relaxation towards an equilibrium distribution,
$f^{eq}_{i}(\bm{x},t)$, taken to be a discrete Maxwellian, with a fixed
relaxation time, $\tau_{0}$,

$\displaystyle\begin{split}f_{i}(\bm{x}+\bm{c}_{i}\Delta_{t},t+\Delta_{t})-f_{i}(\bm{x},t)&=\\\
=-\frac{1}{\tau_{\rm
0}}\left[f_{i}(\bm{x},t)-f^{eq}_{i}(\bm{x},t)\right]&+F_{i}.\end{split}$ (2)

Here, $F_{i}$, is a term introduced to model a macroscopic external forcing
Succi2001 , and $i=0,..,q-1$, indexes the different velocity directions on the
lattice. Eq. (2) is obtained discretizing the Boltzmann equation and selecting
the lattice spacing, $\Delta\bm{x}$, such as divided by the time step,
$\Delta_{t}$, they are equal to the lattice velocity
$\bm{c}_{i}=\Delta\bm{x}/\Delta_{t}$. From those mesoscopic quantities, it is
then possible to recover the macroscopic velocity and density fields by
following the perturbative Chapman-Enskog expansion Wolf-Gladrow2000 . It can
be shown that evolving Eq. 2 is equivalent, up to approximations, to evolving
the weakly compressible NSE for a flow with a density
$\rho(\bm{x},t)=\sum_{i=0}^{q-1}f_{i}(\bm{x},t)$, a velocity
$\bm{u}(\bm{x},t)=\sum_{i=0}^{q-1}f_{i}(\bm{x},t)\bm{c}_{i}/\rho(\bm{x},t)$
and a viscosity directly related to the relaxation time Wolf-Gladrow2000 ,

$\nu_{0}=c_{s}^{2}\Delta_{t}\left(\tau_{0}-\frac{1}{2}\right)\text{,}$ (3)

where, $c_{s}$, is the speed of sound, and $\tau_{0}$ is the adimensional
relaxation time. A numerical validation of the hydrodynamic recovery of BGK-
LBM in the context of 2D homogeneous isotropic turbulence (HIT) was performed
in Tauzin2018  showing good agreement with DNS either in decaying that in
forced regimes. Although this method is adapted to describing various physics
of multi-phase and flows with complex boundaries in a highly scalable fashion,
the BGK-LBM model suffers of numerical instabilities when, $\tau_{\rm
0}\rightarrow\frac{1}{2}$, i.e. $\nu_{\rm 0}\rightarrow 0$, which has made the
study of turbulent flows highly prohibitive for this method Tauzin2018 . To
push LBM towards more turbulent regimes a number of collision operators have
been proposed, see eggels1996direct ; yu2006turbulent ; sagaut2010toward ;
malaspinas2012consistent . Here we focus on the Entropic LBM (ELBM) Karlin1999
; ansumali2002single , which tackles the stability issue by equipping LBM with
an H-theorem. To achieve this results the ELBM differs from BGK-LBM by two
major aspects. First, the equilibrium distribution $\bm{f}^{eq}(\bm{x},t)$ is
not anymore a discretization of the Maxwell-Boltzmann distribution, but it is
calculated as the extremum of a discretized H-function defined as;

$H[\bm{f}]=\sum_{i=0}^{q-1}f_{i}\log\left(\frac{f_{i}}{w_{i}}\right),\,\,\bm{f}=\left\\{f_{i}\right\\}_{i=0}^{q-1},$
(4)

where $w_{i}$ are the weights associated to each lattice direction, under the
constraints of mass and momentum conservation, see Ansumali2003 . The second
difference in ELBM, is that the relaxation time is not a constant anymore but
is modified at every time step in order to enforce the non-monotonicity of H
after the collision. This results in an apparent unconditional stability as
$\nu_{0}\rightarrow 0$ Karlin2015 . It follows that ELBM evolution equations
are,

$\displaystyle\begin{split}f_{i}(x+c_{i}\Delta_{t},t+\Delta_{t})-f_{i}({\bm{x}},t)&=\\\
=-\alpha({\bm{x}},t)\beta[f_{i}({\bm{x}},t)&-f_{i}^{eq}({\bm{x}},t)],\end{split}$
(5)

where $\beta=1/(2\tau_{0})$ is constant, while the new relaxation time,
$\tau_{\text{eff}}(\bm{x},t)=1/(\alpha(\bm{x},t)\beta)$, fluctuates in time
and space through the definition of an entropic parameter $\alpha(\bm{x},t)$.
More recently the ELBM method has been extended to a family of multi-
relaxation time (MRT) lattice Boltzmann models dorschner2016entropic ;
dorschner2017transitional ; dorschner2018fluid . Note that $\alpha$ can be
computed as the solution of the entropic equation
$H[\bm{f}]=H[\bm{f}-\alpha\bm{f}^{neq}]$ which represents the maximum
H-function variation due to a collision, with
$\bm{f}^{neq}=\bm{f}-\bm{f}^{eq}$. Following this approach the computation of
$\alpha(\bm{x},t)$ can be performed via an expensive Newton-Raphson algorithm
for every grid and at every time step of ELBM. To alleviate this problem,
after the original ELBM formulation Karlin1999 , a new version has been
proposed where the computation of the entropic parameter is based on an
analytical formulation derived as an first order expansion of the original
model karlin2014gibbs ; bosch2015entropic . However, to our knowledge, a study
of high-order structure functions in the context of forced 3D HIT, has never
been attempted before using the ELBM original model. In this regard, and also
aiming to measure high-order, extremely sensitive statistical observables, we
implemented the ELBM original formulation, relying on the least number of
approximations, even though computationally more expensive. More details about
ELBM will be given in section II. Let us notice that BGK-LBM is recovered from
Eq. (5) with $\alpha=2$ and the specific Maxwellian expression of
$f_{i}^{eq}$. It is important to stress that the bridge relation described in
Eq. (3) connecting viscosity and relaxation time still holds for fluctuating
quantities, hence we can write,

$\displaystyle\begin{split}\nu_{\text{eff}}\left(\alpha\right)&=c_{s}^{2}\Delta_{t}\left(\frac{1}{\alpha\beta}-\frac{1}{2}\right)=\\\
&=c_{s}^{2}\Delta_{t}\left(\frac{1}{2\beta}-\frac{1}{2}\right)+c_{s}^{2}\Delta_{t}\frac{2-\alpha}{2\alpha\beta}=\\\
&=\nu_{0}+\delta\nu_{\alpha},\end{split}$ (6)

where $\nu_{0}$ represents the constant kinematic viscosity and
$\delta\nu_{\alpha}$ is the fluctuating term. Following this idea, it has been
shown that ELBM is implicitly enforcing a SGS model of an eddy-viscosity type
Malaspinas2008 ; Karlin2015 . In particular, as initially done in
Malaspinas2008 , and then rederived in chapter 4 of tauzin2019implicit ,
performing a third-order Chapman-Enskog perturbative expansion in the limit of
small Knudsen number (Kn), it is possible to obtain a macroscopic
approximation of $\delta\nu_{\alpha}$, which can be written as,

$\delta\nu_{\alpha}^{M}=-c_{s}^{2}\Delta t^{2}\frac{S_{\ell
j}S_{ij}S_{i\ell}}{S_{ij}S_{ij}},$ (7)

where $S_{ij}=1/2(\partial_{j}u_{i}+\partial_{i}u_{j})$ is the strain rate
tensor. The entropic eddy viscosity in Eq. (7) is particularly interesting
because it is not positive definite and can reproduce events of energy
backscatter, which is generally not the case among the other LES closures,
i.e. see the Smagorinsky eddy viscosity Smagorinsky1963 . After the
introduction of the new LES model, see section II, we compare it with standard
Smagorinsky closure and with fully resolved data obtained from DNS. At the
same time, we also present for the first time a quantitative investigation of
the inertial range statistics provided by the mesoscopic ELBM approach in the
context of 3d turbulence. Results provide evidence that ELBM is a good
approximation of the 3d flows up to turbulent regimes never reachable with the
standard BGK-LBM. We found that ELBM guaranties the simulation stability
producing a considerable extension of the inertial range scaling, i.e. the
extension of the energy spectrum power law. Measuring statistical properties
of higher order observables such as the structure functions, we found that
ELBM is also able to reproduce qualitatively the intermittent features of real
3d turbulent flows with an accuracy comparable to the Smagorinsky model. To
get a more accurate estimation of the anomalous exponents more refined models
are required biferale2019self .

The paper is organized as follows. In section II we introduce the details of
the ELBM and the LES models considered in this work. In section III we present
the set of simulations performed. In sect. IV we evaluate the quality of the
two LES closures by comparing them with the real SGS energy transfer measured
from fully resolved DNS. In sect. V we focus on the intermittent properties by
analysing high-order inertial range statistics Arneodo2008 ; benzi2010inertial
. In sect. VI we discuss our conclusions.

##  II Turbulence Modelling

In this section we give a description of the SGS modelling approaches
considered in this work. We start discussing the mesoscopic ELBM approach
highlighting the differences with respect to the standard BGK-LBM. Following
we discuss the new hydrodynamic LES closure inspired by the ELBM macroscopic
approximation first derived in Malaspinas2008 . In the end of this section we
briefly recall the well known Smagorinsky model.

###  II.1 Entropic Lattice Boltzmann Method

Using the same formalism as in Karlin1999 , the ELBM, Eq. (5), can be
rewritten as,

$\displaystyle\begin{split}f_{i}(x+c_{i}&\Delta_{t},t+\Delta_{t})=\\\
=&f_{i}({\bm{x}},t)-\alpha({\bm{x}},t)\beta\left(f_{i}({\bm{x}},t)-f_{i}^{eq}({\bm{x}},t)\right)\\\
=&\left(1-\beta\right)\,f^{pre}_{i}({\bm{x}},t)+\beta\,f_{i}^{mir}({\bm{x}},t)\\\
=&f_{i}^{post}(\bm{x},t),\end{split}$ (8)

where the fluctuating relaxation time is
$\tau_{\text{eff}}({\bm{x}},t)=1/(\alpha({\bm{x}},t)\beta)$, with
$\beta=1/(2\tau_{0})$ and where $\alpha(\bm{x},t)$ is the time and space
dependent, locally-calculated, entropic parameter. The post-collision
distribution, $\bm{f}^{post}(\beta)$, can be understood as a convex
combination between the pre-collision distribution, $\bm{f}^{pre}=\bm{f}$, and
the so-called mirror distribution,
$\bm{f}^{mir}(\alpha)=\bm{f}^{pre}-\alpha\,\bm{f}^{neq}$, with
$\bm{f}^{neq}=\bm{f}^{pre}-\bm{f}^{eq}$, the non-equilibrium part of
$\bm{f}^{pre}$. This convex combination is parametrized by the parameter
$\beta$ in the range $0<\beta<1$ for which we have $0.5<\tau_{0}<+\infty$.
From the definition of the H-functional given in Eq. (4) the discrete
H-theorem can then be expressed as a the local decrease of the H-functional
between the pre-collision and post-collision distributions,

$\displaystyle\begin{split}\Delta H&=H[\bm{f}^{post}]-H[\bm{f}^{pre}]\\\
&=H[(1-\beta)\bm{f}^{pre}+\beta\bm{f}^{mir}(\alpha)]-H[\bm{f}]\leq
0,\end{split}$ (9)

The equilibrium distribution function $\bm{f}^{eq}$ can be calculated as the
extremum of the convex H-functional introduced in Eq. (4), which has an
analytical solution for the D1Q3 lattice, whose tensorial product is solution
for three D2Q9 and the D3Q27 lattice,

$\displaystyle\begin{split}f_{i}^{eq}&(\bm{x},t)=\\\
w_{i}\rho&\prod_{j=1}^{d}\left\\{\left(2-\sqrt{1+\frac{u_{j}^{2}}{c_{s}^{2}}}\right)\left[\frac{\frac{2u_{j}}{\sqrt{3}c_{s}}+\sqrt{1+\frac{u_{j}^{2}}{c_{s}^{2}}}}{1-\frac{u_{j}}{\sqrt{3}c_{s}}}\right]^{\frac{c_{ij}}{\sqrt{3}c_{s}}}\right\\},\end{split}$
(10)

where $d$ is the dimension of the DdQq lattice, $c_{s}$ is the speed of sound,
$w_{i}$ are the weights associated with each lattice direction and
$u_{j}(\bm{x},t)=\sum_{i=0}^{q-1}f_{i}(\bm{x},t)c_{ij}/\rho(\bm{x},t)$ is the
flow macroscopic velocity. It is important to remark that the first three
moments of the entropic equilibrium distribution Eq. (10) are exactly the same
as the one coming from the 3rd order Hermite polynomial expansion of the
Maxwell-Boltzmann equilibrium distribution, namely,

$\displaystyle\begin{split}f^{\rm{eq}}_{i}(\rho,{\bm{u}})=w_{i}\rho&\bigg{(}1+\frac{{\bm{u}}\cdot{\bm{c}}_{i}}{c_{s}^{2}}+\frac{{\bm{u}}{\bm{u}}:{\bm{c}}_{i}{\bm{c}}_{i}-c_{s}^{2}|{\bm{u}}|^{2}}{2c_{s}^{4}}\\\
&+\frac{{\bm{u}}{\bm{u}}{\bm{u}}:\\!\cdot\,{\bm{c}}_{i}{\bm{c}}_{i}{\bm{c}}_{i}-3c_{s}^{2}|{\bm{u}}|^{2}{\bm{u}}\cdot{\bm{c}}_{i}}{6c_{s}^{6}}\bigg{)},\end{split}$
(11)

therefore, it allows the recovery of the same athermal weakly compressible NSE
as in the case of BGK-LBM Ansumali2003 . This recovery, obtained by performing
a Chapman-Enskog expansion at the second order in Kn, is also valid for ELBM
as fluctuation of $\alpha$ around $2$ leads to higher-order terms in Kn,
absorbed in $\mathcal{O}(\text{Kn}^{2})$ tauzin2019implicit . In this work,
following the approach used in Karlin1999 , we calculate $\alpha(\bm{x},t)$ as
the solution of,

$H[\bm{f}^{pre}]=H[\bm{f}^{mir}],$ (12)

which can be estimated via the popular Newton-Raphson algorithm. In this way,
$f_{i}^{post}$ being a convex combination between two distributions,
$f_{i}^{pre}$ and $f_{i}^{mir}$ of equal H-value, and being H at the same time
a convex functional, the monotonic decrease of the H is ensured. Let us stress
that, as it was shown in tauzin2019implicit , the ELBM equation cannot be
considered macroscopically as an approximation to the weakly compressible NSE
with the addiction of a sole eddy viscosity term of the form of Eq. (7).
Indeed, this term appears in a macroscopic equation of motion that requires a
Chapman-Enskog expansion of third order in the Knudsen number, while the NSE
are recovered at the second order. As a consequence a number of extra third-
order terms are part of the implicit ELBM SGS model. This makes the actual
ELBM closure even more complex than a simple eddy viscosity, and in principle,
able to outperform standard methods. On top of this, as already discussed in
the introduction, the macroscopic approximation of the ELBM eddy viscosity,
Eq. (7), has itself a very interesting formulation, being similar to a
Smagorinsky eddy viscosity Smagorinsky1963 , but being not positive-definite
and therefore allows events of energy backscatter, i.e. energy transfer from
the unresolved to the resolved scales. Indeed, while energy in 3d turbulence
is on average cascading from the large towards the small scales, in real flows
there are local events of energy going backward with non-trivial implications
on the statistical properties of the resolved scales.

###  II.2 Large-Eddy Simulations

The LES governing equations can be directly derived by filtering each term of
the incompressible NSE, see Lesieur2005 ; Meneveau2000 . The filtering
operation consist in a convolution between the full velocity field and a
filter kernel. There are several choices that can be made for the filter
kernel, see Pope2001 , in this work, we consider a “sharp spectral cutoff” in
Fourier space. This choice is convenient for two reasons, first, the sharp
cutoff produces a clear separation between resolved and sub-grid scales,
defined respectively as all scales above and below the cutoff wavenumber,
$k_{c}$. Second, it is a Galerkin projector that produces the same results
when operating multiple times on a field, which allows to have a clear scale
separation also in a dynamical sense, namely it projects on the same support
all terms (non-linear ones included) of the equations of motion, see
buzzicotti2018effect . In the following we briefly sketch the main operations
required to arrive at the LES equations. Given a filter kernel
$G_{\Delta}(\bm{x})$, the filtered velocity
$\overline{\boldsymbol{u}}(\bm{x},t)$ can be defined by the following
convolution operation,

$\displaystyle\begin{split}\overline{\boldsymbol{u}}(\bm{x},t)\equiv\int_{\Omega}d\bm{y}\
G_{\Delta}(|\bm{x}-\bm{y}|)\,\bm{u}(\bm{y},t)=\\\
=\sum_{\bm{k}\in\mathbb{Z}^{3}}\hat{G}_{\Delta}(|\bm{k}|)\boldsymbol{\hat{u}}(\bm{k},t)e^{i\bm{k}\bm{x}}\,\end{split}$
(13)

where $\hat{G}_{\Delta}$ is the Fourier transform of $G_{\Delta}$, and
$\Delta\sim\pi/k_{c}$, is the filter cutoff scale, see Pope2001 . Applying the
filtering operation to all terms in the Navier-Stokes equations we get,

$\partial_{t}\overline{\boldsymbol{u}}+\nabla\cdot(\overline{\overline{\boldsymbol{u}}\otimes\overline{\boldsymbol{u}}})=-\nabla\overline{{p}}-\nabla\cdot\overline{{\tau}}^{\Delta}(\bm{u},\bm{u})+\nu\Delta\overline{\boldsymbol{u}}\
.$ (14)

Here, we have introduced the SGS tensor,
$\overline{{\tau}}^{\Delta}(\bm{u},\bm{u})$, defined as,

$\overline{{\tau}}^{\Delta}_{ij}(\bm{u},\bm{u})=\overline{u_{i}u_{j}}-\overline{\overline{u}_{i}\overline{u}_{j}},$
(15)

which is the only term of Eq. (14) that depends on SGS scales. Hence, it is
the only term that needs to be replaced by a model to close the equations in
terms of the resolved-scales dynamics. From
$\overline{{\tau}}^{\Delta}(\bm{u},\bm{u})$, we can easily get the exact
formulation of the SGS energy transfer $\overline{{\Pi}}^{\Delta}$, namely,
the energy transfer across the filter scale produced by the real non-linear
coupling in the NSE. To do so we need to multiply with a scalar product each
term of Eq. (14) and the velocity field to obtain the filtered energy balance
equations;

$\frac{1}{2}\partial_{t}(\overline{u}_{i}\overline{u}_{i})+\partial_{j}A_{ij}+{\overline{\Pi}}^{\Delta}_{L}=-\overline{{\Pi}}^{\Delta}.$
(16)

The terms on the lhs of Eq. (16) are defined respectively as
$\partial_{j}A_{j}=\partial_{j}\overline{u}_{i}(\overline{\overline{u}_{i}\overline{u}_{j}}+\overline{p}\delta_{ij}+\overline{{\tau}}^{\Delta}_{ij}-\frac{1}{2}\overline{u}_{i}\overline{u}_{j})$
and
${\overline{\Pi}}^{\Delta}_{L}=-\partial_{j}\overline{u}_{i}\left(\overline{\overline{u}_{i}\overline{u}_{j}}-\overline{u}_{i}\overline{u}_{j}\right)$.

ID |  | Eddy viscosities |  | Stress tensors |  |  | Energy transfers
---|---|---|---|---|---|---|---
DNS |  | — |  | $\overline{{\tau}}^{\Delta}_{ij}=\overline{u_{i}u_{j}}-\overline{\overline{u}_{i}\overline{u}_{j}}$ |  |  | $\overline{{\Pi}}^{\Delta}=-\partial_{j}\overline{u}_{i}\,\,\overline{{\tau}}^{\Delta}_{ij}(\bm{u},\bm{u})$
S-LES |  | $\nu^{S}_{e}=(C_{S}\Delta)^{2}\sqrt{2S_{ij}S_{ij}}$ |  | $\overline{{\tau}}^{\text{S}}_{ij}=-2\nu_{e}^{S}\bar{S}_{ij}$ |  |  | $\overline{{\Pi}}^{\text{S}}_{\text{LES}}=-\overline{{\tau}}^{\text{S}}_{ij}\overline{S}_{ij}$
E-LES |  | $\nu^{E}_{e}=(C_{E}\Delta)^{2}\frac{S_{\ell j}S_{ji}S_{i\ell}}{S_{ij}S_{ij}}$ |  | $\overline{{\tau}}^{\text{E}}_{ij}=-2\nu_{e}^{E}\bar{S}_{ij}$ |  |  | $\overline{{\Pi}}^{\text{E}}_{\text{LES}}=-\overline{{\tau}}^{\text{E}}_{ij}\overline{S}_{ij}$
ELBM |  | $\delta\nu_{\alpha}=c_{s}^{2}\frac{2-\alpha}{2\alpha\beta}\Delta_{t}$ |  | $\overline{{\tau}}^{\alpha}_{ij}=-2\delta\nu_{\alpha}\bar{S}_{ij}$ |  |  | $\overline{{\Pi}}^{\text{E}}_{\text{LBM}}=-\overline{{\tau}}^{\alpha}_{ij}\overline{S}_{ij}$
Table 1: Summary of definitions of eddy viscosities, sub-grid scales stresses
and energy transfers. The ID column indicates the names of the corresponding
set of simulations, see sect. III, in particular S-LES and E-LES correspond to
the hydrodynamical LES respectively with Smagorinsky and macroscopic ELBM
model, while with ELBM we refer to the SGS energy transfer measured from the
mesoscopic quantities.

As shown in leonard1974energy ; buzzicotti2018effect , to get the correct
contribution to the energy transfer, it is important to distinguish between
${\overline{\Pi}}^{\Delta}_{L}$ and the actual SGS energy transfer
$\overline{{\Pi}}^{\Delta}$ because the former depends only on resolved-scales
quantities and does not contribute to the mean energy flux across the cutoff
scale. On the other hand,

$\overline{{\Pi}}^{\Delta}=-\partial_{j}\overline{u}_{i}\,\,\overline{{\tau}}^{\Delta}_{ij}(\bm{u},\bm{u})=-\partial_{j}\overline{u}_{i}\left(\overline{u_{i}u_{j}}-\overline{\overline{u}_{i}\overline{u}_{j}}\right),$
(17)

is the flux which depends on both the SGS and the resolved scales. In this
work, as already mentioned, we consider as possible closure the Smagorinsky
LES model (referred to as S-LES in the sequel),
$\overline{{\tau}}^{\text{S}}_{ij}(\overline{\boldsymbol{u}},\overline{\boldsymbol{u}})$,
which aims to model the deviatoric part of the stress tensor,
$\overline{{\tau}}^{\Delta}_{ij}-\frac{1}{3}\overline{{\tau}}^{\Delta}_{kk}\delta_{ij}\,\,\rightarrow\,\,\overline{{\tau}}^{\text{S}}_{ij}$,
as follows,

$\overline{{\tau}}^{\text{S}}_{ij}=-2\nu_{e}^{S}(\bm{x},t)\bar{S}_{ij},\quad\nu^{S}_{e}=(C_{S}\Delta)^{2}\sqrt{2S_{ij}S_{ij}}$
(18)

where,
$\bar{S}_{ij}=1/2(\partial_{j}\overline{u}_{i}+\partial_{i}\overline{u}_{j})$,
is the resolved scales strain-rate tensor, $\nu_{e}^{S}$ is the Smagorinsky
eddy viscosity depending on the filter cutoff scale $\Delta$ and the non-
dimensional factor $C_{S}$. From the definition of the macroscopic
approximation of the ELBM eddy viscosity in Eq. (7), we can now define the
hydrodynamic ELBM-LES model (called E-LES in the sequel),

$\overline{{\tau}}^{\text{E}}_{ij}=-2\nu_{e}^{E}(\bm{x},t)\bar{S}_{ij},\quad\nu^{E}_{e}=(C_{E}\Delta)^{2}\frac{S_{\ell
j}S_{ji}S_{i\ell}}{S_{ij}S_{ij}},$ (19)

where $C_{E}$ is the entropic dimensionless coefficient. Comparing the
definition of $\nu^{E}_{e}$ with $\delta\nu_{\alpha}^{M}$, in Eq. (7), we can
see that they both have the same functional dependency on the strain-rate
tensor, but different signs and multiplicative constants. In particular, the
minus sign of the E-LES closure has been absorbed in the definition of
$\overline{{\tau}}^{\text{E}}_{ij}$ to align with the Smagorinsky closure
formulation. Let us stress that the E-LES model has the same scaling as the
Smagorinsky model, proportional to the strain rate, but it is not positive
definite. From the above definitions of the S-LES and E-LES models the two
corresponding SGS energy transfers can be written as,

$\overline{{\Pi}}^{\text{S}}_{\text{LES}}=-\overline{{\tau}}^{\text{S}}_{ij}\overline{S}_{ij};\qquad\overline{{\Pi}}^{\text{E}}_{\text{LES}}=-\overline{{\tau}}^{\text{E}}_{ij}\overline{S}_{ij}.$
(20)

To compare the behaviour of the mesoscopic ELBM model with respect to the two
hydrodynamical approaches just introduced, we can approximate the SGS energy
transfer from the ELBM as,

$\overline{{\Pi}}^{\text{E}}_{\text{LBM}}=-2\delta\nu_{\alpha}\overline{S}_{ij}\overline{S}_{ij},$
(21)

where $\overline{{\Pi}}^{\text{E}}_{\text{LBM}}$ stands for ELBM SGS energy
transfer and
$\delta\nu_{\alpha}=c_{s}^{2}\Delta_{t}\frac{2-\alpha}{2\alpha\beta}$ is the
mesoscopic fluctuating viscosity depending on $\alpha(\bm{x},t)$. The strain
rate tensor can be measured from the ELBM data after the calculation of the
macroscopic velocity in terms of the mesoscopic ones, namely,
$\bm{u}(\bm{x},t)=\sum_{i=0}^{q-1}f_{i}(\bm{x},t)\bm{c}_{i}/\rho(\bm{x},t)$. A
summary of these SGS energy transfer definitions with their respective SGS
tensors, eddy viscosities is given in table 1. It is worth noting that the
ELBM in the limit of small Knudsen numbers is not equivalent to the entropic
LES. Indeed, as previously mentioned, the eddy viscosity term appears in the
Chapman-Enskog expansion at the third order in Kn along with various extra
terms that are not contained in the entropic LES formulation.

##  III Numerical Simulations

All simulations performed in this work are intended to model HIT on a three
dimensional domain with periodic boundary conditions. In the following we
provide some details about the sets of simulations performed with the
different modelling techniques. Concerning the lattice Boltzmann simulation
with entropic formulation of the relaxation time, ELBM, we have conducted a
set of 3d simulations with a number of $512$ collocation points along each
spatial direction. To reach stationarity the flow is forced at large scales,
$1\leq|\bm{k}|\leq 2$ with a constant and isotropic forcing. More precisely we
have used in all simulations the same forcing, defined in Fourier space with
constant phases and amplitudes, added isotropically to all wave-vectors at
large-scales. To ensure incompressibility the forcing is projected on its
solenoidal component. The ELBM simulation uses a lattice with 27 discrete
velocities (see Fig. 1), the D3Q27 Succi2001 ; Wolf-Gladrow2000 ; Kruger2017 .
The spectral forcing is implemented using the exact-difference method forcing
scheme Kuperstokh2004  for a relaxation time $\tau_{0}=0.5003$ corresponding
to $\beta\approx 0.9994$.

Figure 1: Schematic representation of the D3Q27 lattice stencil used for the
ELBM simulation.

Considering the LES we have performed two sets of pseudo-spectral fully
dealiased simulations on a domain $\Omega=[0,2\pi]^{3}$ with periodic boundary
conditions both at the resolution of $512^{3}$ grid points. The first LES is
equipped with the Smagorinsky model (S-LES), see Eq. (18), and a second with
the macroscopic formulation of the entropic eddy viscosity, see Eq. (19),
(E-LES). As discussed above, all simulations are forced with the same
isotropic constant forcing mechanisms acting only on the larger system scales
($1\leq|\bm{k}|\leq 2$). In the expression of Smagorinsky eddy viscosity Eq.
(18), we use the standard value of $C_{S}=0.16$, while for the entropic eddy
viscosity Eq. (19), we use $C_{E}=0.45$, found to be optimal values for the
best compromise between the maximization of inertial range extension and the
minimization of spurious effects produced by the model buzzicotti2020synch .
For both we have $\Delta=\pi/k_{max}\approx 0.0184$ with $k_{max}=171$, which
comes from the $2/3$ rule for the dealiasing projection Patterson71 .
Additionally, as a reference, we have run a pseudo-spectral fully resolved DNS
of the NSE with the same forcing scheme on the same 3d domain
$\Omega=[0,2\pi]^{3}$, using a number of $512^{3}$ (DNSx1) and $1024^{3}$
(DNSx2) collocation points. The resolution in both DNS is kept such as
$\eta_{\alpha}/dx\simeq 0.7$, where $dx$ is the grid spacing and
$\eta_{\alpha}=(\nu^{3}/\epsilon)^{1/4}$ is the Kolmogorov microscale Borue95
with $\varepsilon$ denoting the mean energy dissipation rate. In order to
create ensembles of statistically independent data all the ELBM, LES and DNS
are sampled on time intervals of one large-scale eddy turnover time after
reaching a statistically stationary state. It is worth mentioning that both
the ELBM and E-LES simulations remain stable even though their modeling terms
are not purely dissipative being their eddy viscosities not positive definite.
In Fig. 2 we show the time-averaged energy spectra, $E(k)$, for all
simulations. It is visible that all modelled simulations have an extended
inertial range with respect to DNS at the same resolution (DNSx1) with an
inertial range slope close to the Kolmogorov prediction of $k^{-5/3}$
Frisch1995 . Let us stress that the ELBM spectrum reaches a maximum wavenumber
higher than the pseudo-spectral data, this is because in the ELBM simulations
there is not any dealiasing operation. Anyway also in the ELBM case the
spectrum loses the Kolmogorov inertial range scaling at wavenumber larger that
the pseudo-spectral dealiasing cutoff.

Figure 2: Time-averaged spectra for the conducted simulations at $512^{3}$
grid points, measured from the mesoscopic ELBM simulation (empty squares, red
color), the hydrodynamical LES with entropic model (E-LES, empty circles) and
with Smagorinsky model (S-LES, empty triangles). The energy spectra from fully
resolved DNS at $512^{3}$ (DNSx1) and $1024^{3}$ (DNSx2) are presented
respectively with full triangles and full circles. The curves are shifted
vertically for the sake of data presentation. The Kolmogorov predicted slope
of $k^{-5/3}$ is given as a reference Frisch1995 .

##  IV SGS energy transfer analysis

In this section we provide a statistical comparison of the SGS energy
transfers measured in the modelled simulations together with the original SGS
transfer measured a priori from higher resolution DNS, see Table 1. Before
staring the comparison between macroscopic and mesoscopic simulations, we
analyse the quality of the approximation made in the Chapman-Enskog expansion
to obtain the macroscopic formulation of the ELBM eddy viscosity
Malaspinas2008 . In this direction, we have computed the SGS energy transfer
defined in Eq. (21) using the two different definitions of the fluctuating
viscosity. Namely, either the correct definition, $\delta\nu_{\alpha}$,
depending on the entropic parameter or its third order expansion in the limit
of small Kn, $\delta\nu_{\alpha}^{M}$ see Eq. (7). Their statistical
comparison is shown in Fig. 3, where on the left panel we show the probability
density functions (PDF) measured from the ELBM SGS energy transfer using the
two different formulations. Here we can see that the PDFs once re-scaled by
their standard deviations have almost an identical shape. From the center and
right panels of the same figure, we can qualitatively see two visualizations
of the SGS energy transfers measured by selecting the same plane of the
velocity field. From these visualizations we can appreciate that there is a
very high spatial correlation between them. These results suggest that the
approximation of neglecting the extra third order terms coming from Chapman-
Enskog expansion is a good approximation of the ELBM eddy viscosity.

Figure 3: Standardized PDFs of the ELBM SGS energy transfers measured from the
correct fluctuating viscosity, $\delta\nu_{\alpha}$, and its approximated
formulation, $\delta\nu_{\alpha}^{M}$ (left panel). Visualization of a plane
of ELBM SGS energy transfer measured using the approximated (center panels)
and correct definitions of eddy viscosity (right panel).

Let us now analyse the statistics of the SGS energy transfers, comparing them
also with the statistics of the real SGS energy transfer, see Eq. (17),
measured a priori from fully resolved DNSx2. To obtain the a priori
$\overline{{\Pi}}^{\Delta}$ we filter the velocity field with a sharp
projector in Fourier space with a cutoff at the maximum wavenumber allowed in
the modelled simulations, which corresponds to the dealiasing cutoff
($k_{max}=171$). As known, the presence of a forward energy cascade, as in 3d
turbulence, reflects in a skewed PDF of the a priori
$\overline{{\Pi}}^{\Delta}$, see gotoh2005statistics ; buzzicotti2018energy ;
buzzicotti2018effect , see Fig. 4. Instead, the negative tail describes the
presence of intense backscattering events with fluctuations up to two orders
of magnitude larger than the standard deviation. The main remarkable
difference between the different models considered here is that, as expected,
the ELBM and E-LES produce backscatter events, while the Smagorinsky model is
positive definite in its energy formulation and it produces a zero tail in the
negative region. Let us notice that ELBM mesoscopic model shows qualitatively
a better overlap with respect to the DNS data.

Figure 4: Standardized PDF of the SGS energy transfer measured from the a
posteriori data obtained via the ELBM simulations (empty squares, red color),
the LES with hydrodynamical entropic closure (E-LES, empty circles) and from
LES with Smagorinsky model (S-LES, empty triangles). For comparison the PDF
measured for the real SGS energy transfer measured a priori by filtering data
from higher resolution simulations is presented (DNSx2, full circles).

##  V Inertial range statistics

In this last section we analyse the inertial range statistics by measuring the
longitudinal velocity increments defined as
$\delta_{r}u=(\bm{u}(\bm{x}+\bm{r})-\bm{u}(\bm{x}))\cdot\bm{r}/r$. In this way
we can quantify the effects produced by the different models at different
scales, $r=|\bm{r}|$, hence, we can get an accurate estimation of the quality
of the models in capturing the correct intermittent properties of the NSE. We
study the scaling properties of the longitudinal structure functions (SF)
defined as,

$S_{p}(r)\equiv\langle[\delta_{r}u]^{p}\rangle$ (22)

where the angular brackets indicate the ensemble average, that assuming
spatio-temporal ergodicity can be evaluated averaging over space and time,
$\langle(...)\rangle=\frac{1}{V}\frac{1}{T}\int_{V}\int_{t_{0}}^{t_{0}+T}(...)\,d\bm{x}dt$.
In the limit of large Reynolds number, where $r$ can be taken arbitrarily
small the structures function follows a powerlaw scaling behavior,
$S_{p}(r)\sim r^{\xi_{p}}$ Frisch1995 , where a p-th order scaling exponent
that, according to the phenomenological theory of Kolmogorov (K41)
Kolmogorov1941 , is $\xi_{p}=p/3$. Nevertheless, both experimental and
numerical studies have highlighted as bi-product of intermittency the presence
of anomalous exponents in turbulent data, with important deviations from the
K41 predicted values Gotoh2002 ; sinhuber2017dissipative ; benzi2010inertial ;
biferale2019self . On the other hand, to get an accurate measurement of these
exponents is extremely difficult. The reason is twofold, first it is required
to have large scaling range (very well resolved simulations) and second it is
simultaneously required to have large statistical ensemble. The first question
we ask here is connected to the first of the aforementioned problems, namely
whether those models are able or not to extend the length of inertial range of
scales in our simulations. To answer this in the left column of Fig. 5 we show
the 2nd and 4th order structure functions measured from all modelled and fully
resolved simulations. In the right column on the same figure we show the local
scaling exponents,

$\xi_{p}(r)=\frac{d\,\log S_{n}(r)}{d\,\log(r)}\,,$

measured form the structure functions shown on the left side. From these plots
it is evident that all modelled simulations present an extension of the
inertial range with respect to the same resolution, similar to the inertial
range observed in the simulations with a number of grid points two times
larger along each spatial direction, DNSx2.

Figure 5: Second-order longitudinal structure functions (left panels) and
corresponding local slopes (right panels) for the conducted simulations at
$512^{3}$ grid points, using ELBM (empty squares, red color), E-LES with
entropic inspired model (empty circles), S-LES with Smagorinsky model (empty
squares), DNSx1 at $512^{3}$ grid points (full triangles) and DNSx2 at
$1024^{3}$ grid points (full circles). The straight line corresponds to the
K41 prediction in the inertial range ($\xi_{p}=p/3$), while the dashed line
corresponds to the intermittent measure as reported from the literature
Gotoh2002 ; biferale2019self .

Which means that the models allows to save a factor $8$ in terms of the number
of degrees of freedom required in the simulations. Considering that also the
time-step needs to be changed accordingly in the higher resolution DNS to
resolve the smaller time-scales, it means that the modelled simulations are
more than an order of magnitude cheaper than a DNS with the same inertial
range. On the other hand, by measuring the local scaling exponents as reported
on the right column of Fig. 5, we can see that the inertial range extension
produced by the model is not as accurate as the fully resolved DNS. This is a
problem if we want the model to correctly describe intermittency. Let us
stress that the correction to the K41 prediction at the level of the second
and fourth order exponents is very small, hence a model needs to be extremely
accurate to correctly capture the intermittent scaling biferale2019self . In
both right panels of Fig. 5 we report with solid line the $\xi_{p}$ value of
the K41 prediction, and with dotted lines the values measured from DNS as
reported in Gotoh2002 ; sinhuber2017dissipative . To highlight intermittency
we can look at the ratio between SF at different orders. In particular, any
systematic non-linear dependency of $\xi_{p}$ vs $p$, will introduce a scale-
dependency in the Kurtosis, defined by the dimensionless ratio among fourth
and second order SF,

$K(r)=\frac{S_{4}(r)}{\left[S_{2}(r)\right]^{2}}.$ (23)

In Fig. 6 we see that in all simulations, at large scales the increments are
Gaussian ($K\sim 3$), while the Kurtosis quickly increases, decreasing the
scale. This observation shows non-self similarity in the statistics of all
data. It is interesting to observe that at this level the inertial range of
scale observed in the DNSx2 simulation are well captured by all closures up to
the dissipative scales $r\approx 0.1$ where deviations from the DNSs and
models arise.

Figure 6: Kurtosis of the velocity increment for the simulations at $512^{3}$
grid points, using ELBM (empty squares, red color), E-LES with the entropic
inspired model (empty circles), S-LES with Smagorinsky model (ampty
triangles), DNSx1 at $512^{3}$ grid points (full triangles) and DNSx2 at
resolution of $1024^{3}$ (full circles). The dashed horizontal line at 3
corresponds to the value of a Gaussian distribution.

Going further, we measure the most refined quantity we can observe to quantify
intermittency, namely the local scaling exponent in Extended Self-Similarity
Benzi1993 ,

$\zeta(r)=\frac{\xi_{4}}{\xi_{2}}.$ (24)

A linear K41 behavior would recover in the inertial range a plateau value of
$\zeta$ equal to $2$. The correction, accounting for intermittency, measured
in both experimental and DNS data gives the plateau for $\zeta$ at the value
of $1.86$ Gotoh2002 ; sinhuber2017dissipative ; biferale2019self . As we can
see in Fig. 7, all models show deviations from the K41 self-similar prediction
meaning that they all capture the non-self-similarity of the turbulent
inertial range. However none of them is really accurate enough to extend the
length of the plateau displayed, hence to improve the prediction obtainable
from the DNS at the same resolution. Indeed, if we compare the modelled data
with the fully resolved simulations we can see that former are not showing any
flat plateau in the inertial range and we cannot estimate precisely the
correction to K41 of the structure function scaling exponents.

Figure 7: Extended Self-Similarity, $\zeta(r)$, for the simulations at
$512^{3}$ grid points, using ELBM (empty squares, red color), the LES with
entropic inspired model E-LES (empty circles), the LES with Smagorinsky model
S-LES (empty triangles), the fully resolved DNSx1 (full triangles) and the
fully resolved DNS at $1024^{3}$ grid points DNSx2 (full circles). The
straight line corresponds to the K41 prediction in the inertial range equal to
$2$, while the dashed line corresponds to the intermittent measure as reported
from numerical Gotoh2002  and experimental data sinhuber2017dissipative .

It is interesting to point out that the models show a very similar accuracy up
to this last analysis. This suggest the backscatter events of energy
introduced by the entropic closures are not accurate enough to improve
quantitatively the statistics of the Smagorinsky model. This results supports
the observation that intermittency in turbulent flows comes as a result of
highly non-trivial correlations among all degrees of freedom at different
scales buzzicotti2016phase ; buzzicotti2016lagrangian ; lanotte2015turbulence
. The observation that for all models we have a very similar inertial range
statistics goes in agreement with the common property of these models to have
the same scaling, proportional to the strain rate tensor.

##  VI Conclusions

In this paper, we performed a quantitative assessment of the ELBM capabilities
in the modelling of 3d homogeneous isotropic turbulent flows by comparing the
inertial range statistics of ELBM data with the one of high resolution DNS of
the NSE. We also compared the quality of ELBM with respect to the
hydrodynamical Smagorinsky model, popular in the realm of LES. Furthermore, in
this work we have proposed and investigated for the first time, a new
hydrodynamical closure for LES simulation inspired from the macroscopic
approximation of the ELBM model introduce by Malaspinas2008 . We found that
ELBM extends the length of the inertial range with respect to the fully
resolved DNS, allowing to reduce the computational cost by an order of
magnitude and at the same time preserving the simulation stability. Results
showed that, in both the macroscopic and mesoscopic formulations, ELBM is able
to reproduce an inertial range with a non self-similar dynamics. ELBM captures
the correct deviations from the large scales Gaussian statistics as observed
in the fully resolved DNS, with an accuracy comparable to the one produced by
the Smagorinsky model. From the measure of the structure functions scaling
exponents in ESS, we have highlighted the limitations of these models to get
with high accuracy the turbulence corrections to the Kolmogorov scaling. In
this context we found that the modelled data are not producing the same
inertial range plateau as observed in the fully resolved DNS and experiments.
To conclude, we found that ELBM suffers in the modeling of extreme and rare
intermittent fluctuations, while on the other hand it is very efficient in
modeling the mean properties of 3d turbulence. Which makes ELBM a good
candidate for the modeling of 3d turbulent flows in complex geometries.

## Acknowledgement

The authors thank Prof. Luca Biferale for inspiration and many useful
discussions. This work was supported also by the European Unions Framework
Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie
Skłodowska-Curie grant [grant number 642069] and by the European Research
Council under the ERC grant [grant number 339032]. The authors would like to
thanks Prof. Dirk Pleiter as well as the Juelich Supercomputing Center for
providing access to the JURON cluster.

## References

  * [1] O. Malaspinas, M. Deville, and B. Chopard.  Towards a physical interpretation of the entropic Lattice Boltzmann method.  Physical Review E, 78(6):066705, 2008.
  * [2] Horia Dumitrescu and Vladimir Cardos.  Rotational effects on the boundary-layer flow in wind turbines.  AIAA journal, 42(2):408–411, 2004.
  * [3] P. A. Davidson.  Turbulence: An Introduction for Scientists and Engineers.  Oxford University Press, 2015.
  * [4] S. B. Pope.  Turbulent Flows.  Cambridge University press, 2000.
  * [5] Adrian E Gill.  Atmosphere—ocean dynamics.  Elsevier, 2016.
  * [6] Alexandros Alexakis and Luca Biferale.  Cascades and transitions in turbulent flows.  Physics Reports, 767:1–101, 2018.
  * [7] Sydney A Barnes.  An assessment of the rotation rates of the host stars of extrasolar planets.  The Astrophysical Journal, 561(2):1095, 2001.
  * [8] U. Frisch.  Turbulence : the legacy of A.N. Kolmogorov.  Cambridge University Press, 1995.
  * [9] Rodney O Fox.  Computational models for turbulent reacting flows.  Cambridge university press, 2003.
  * [10] Takashi Ishihara, Toshiyuki Gotoh, and Yukio Kaneda.  Study of high–reynolds number isotropic turbulence by direct numerical simulation.  Annual Review of Fluid Mechanics, 41:165–180, 2009.
  * [11] PK Yeung, XM Zhai, and Katepalli R Sreenivasan.  Extreme events in computational turbulence.  Proceedings of the National Academy of Sciences, 112(41):12633–12638, 2015.
  * [12] R Benzi, L Biferale, R Fisher, D Lamb, and F Toschi.  Inertial range eulerian and lagrangian statistics from numerical simulations of isotropic turbulence.  2010\.
  * [13] Kartik P Iyer, Katepalli R Sreenivasan, and PK Yeung.  Reynolds number scaling of velocity increments in isotropic turbulence.  Physical Review E, 95(2):021101, 2017.
  * [14] C. Meneveau and J. Katz.  Scale-invariance and turbulence models for large-eddy simulation.  Annual Review of Fluid Mechanics, 32(1):1–32, 2000.
  * [15] Olga Filippova, Sauro Succi, Francesco Mazzocco, Cinzio Arrighetti, Gino Bella, and Dieter Hänel.  Multiscale lattice boltzmann schemes with turbulence modeling.  Journal of Computational Physics, 170(2):812–829, 2001.
  * [16] Y.-H. Dong and P. Sagaut.  A study of time correlations in Lattice Boltzmann-based large-eddy simulation of isotropic turbulence.  Physics of Fluids, 20(3):035105, 2008.
  * [17] Y.-H. Dong, P. Sagaut, and S. Marie.  Inertial consistent subgrid model for large-eddy simulation based on the Lattice Boltzmann method.  Physics of Fluids, 20(3):035104, 2008.
  * [18] Tomas Bohr, Mogens H Jensen, Giovanni Paladin, and Angelo Vulpiani.  Dynamical systems approach to turbulence.  Cambridge University Press, 2005.
  * [19] M. Lesieur, O. Métais, and P. Comte.  Large-eddy simulations of turbulence.  Cambridge University Press, 2005.
  * [20] Sybren Ruurds De Groot and Peter Mazur.  Non-equilibrium thermodynamics.  Courier Corporation, 2013.
  * [21] Roberto Benzi, Sauro Succi, and Massimo Vergassola.  The lattice boltzmann equation: theory and applications.  Physics Reports, 222(3):145–197, 1992.
  * [22] S. Succi.  The Lattice Boltzmann Equation for Fluid Dynamics and Beyond.  Oxford University Press, 2001.
  * [23] Uriel Frisch, Brosl Hasslacher, and Yves Pomeau.  Lattice-gas automata for the navier-stokes equation.  Physical review letters, 56(14):1505, 1986.
  * [24] D. A. Wolf-Gladrow.  Lattice-gas cellular automata and Lattice Boltzmann models : an introduction.  Springer, 2000.
  * [25] Yue-Hong Qian, S Succi, and SA Orszag.  Recent advances in lattice boltzmann computing.  In Annual reviews of computational physics III, pages 195–242. World Scientific, 1995.
  * [26] Shiyi Chen and Gary D Doolen.  Lattice boltzmann method for fluid flows.  Annual review of fluid mechanics, 30(1):329–364, 1998.
  * [27] P. L. Bhatnagar, E. P. Gross, and M. Krook.  A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems.  Phys. Rev., 94:511–525, 1954.
  * [28] G. Tauzin, L. Biferale, M. Sbragaglia, A. Gupta, F. Toschi, A. Bartel, and M. Ehrhardt.  A numerical tool for the study of the hydrodynamic recovery of the Lattice Boltzmann method.  Computers & Fluids, 172:241 – 250, 2018.
  * [29] Jack GM Eggels.  Direct and large-eddy simulation of turbulent fluid flow using the lattice-boltzmann scheme.  International journal of heat and fluid flow, 17(3):307–323, 1996\.
  * [30] Huidan Yu, Li-Shi Luo, and Sharath S Girimaji.  Les of turbulent square jet flow using an mrt lattice boltzmann model.  Computers & Fluids, 35(8-9):957–965, 2006.
  * [31] Pierre Sagaut.  Toward advanced subgrid models for lattice-boltzmann-based large-eddy simulation: Theoretical formulations.  Computers & Mathematics with Applications, 59(7):2194–2199, 2010\.
  * [32] Orestis Malaspinas and Pierre Sagaut.  Consistent subgrid scale modelling for lattice boltzmann methods.  Journal of Fluid Mechanics, 700:514–542, 2012.
  * [33] I. V. Karlin, A. Ferrante, and H. C. Öttinger.  Perfect entropy functions of the Lattice Boltzmann method.  Europhysics Letters (EPL), 47(2):182–188, 1999.
  * [34] Santosh Ansumali and Iliya V Karlin.  Single relaxation time model for entropic lattice boltzmann methods.  Physical Review E, 65(5):056312, 2002.
  * [35] S. Ansumali, I. V. Karlin, and H. C. Öttinger.  Minimal entropic kinetic models for hydrodynamics.  Europhysics Letters (EPL), 63(6):798–804, 2003.
  * [36] I. V. Karlin, F. Bösch, S. S. Chikatamarla, and S. Succi.  Entropy-assisted computing of low-dissipative systems.  Entropy, 17(12):8099–8110, 2015.
  * [37] Benedikt Dorschner, Fabian Bösch, Shyam S Chikatamarla, Konstantinos Boulouchos, and Ilya V Karlin.  Entropic multi-relaxation time lattice boltzmann model for complex flows.  Journal of Fluid Mechanics, 801:623, 2016.
  * [38] Benedikt Dorschner, Shyam S Chikatamarla, and Iliya V Karlin.  Transitional flows with the entropic lattice boltzmann method.  Journal of Fluid Mechanics, 824:388, 2017.
  * [39] Benedikt Dorschner, Shyam S Chikatamarla, and Ilya V Karlin.  Fluid-structure interaction with the entropic lattice boltzmann method.  Physical Review E, 97(2):023305, 2018.
  * [40] Ilya V Karlin, Fabian Bösch, and SS Chikatamarla.  Gibbs’ principle for the lattice-kinetic theory of fluid dynamics.  Physical Review E, 90(3):031302, 2014.
  * [41] Fabian Bösch, Shyam S Chikatamarla, and Ilya V Karlin.  Entropic multirelaxation lattice boltzmann models for turbulent flows.  Physical Review E, 92(4):043309, 2015.
  * [42] Guillaume Tauzin.  Implicit Sub-Grid Scale Modeling within the Entropic Lattice Boltzmann Method in Homogeneous Isotropic Turbulence.  PhD thesis, Universität Wuppertal, Fakultät für Mathematik und Naturwissenschaften …, 2019.
  * [43] J. Smagorinsky.  General circulation experiments with the primitive equations.  91(3):99–194, 1963.
  * [44] Luca Biferale, Fabio Bonaccorso, Michele Buzzicotti, and Kartik P Iyer.  Self-similar subgrid-scale models for inertial range turbulence and accurate measurements of intermittency.  Physical review letters, 123(1):014503, 2019.
  * [45] A. Arnèodo, R. Benzi, J. Berg, L. Biferale, E. Bodenschatz, A. Busse, E. Calzavarini, B. Castaing, M. Cencini, L. Chevillard, R. T. Fisher, R. Grauer, H. Homann, D. Lamb, A. S. Lanotte, E. Lévèque, B. Lüthi, J. Mann, N. Mordant, W.-C. Müller, S. Ott, N. T. Ouellette, J.-F. Pinton, S. B. Pope, S. G. Roux, F. Toschi, H. Xu, and P. K. Yeung.  Universal intermittent properties of particle trajectories in highly turbulent flows.  Phys. Rev. Lett., 100:254504, 2008.
  * [46] M Buzzicotti, M Linkmann, H Aluie, L Biferale, J Brasseur, and C Meneveau.  Effect of filter type on the statistics of energy transfer between resolved and subfilter scales from a-priori analysis of direct numerical simulations of isotropic turbulence.  Journal of Turbulence, 19(2):167–197, 2018.
  * [47] A Leonard et al.  Energy cascade in large-eddy simulations of turbulent fluid flows.  Adv. Geophys. A, 18(A):237–248, 1974.
  * [48] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen.  The Lattice Boltzmann Method.  Graduate Texts in Physics. Springer International Publishing, Cham, 2017\.
  * [49] A. L. Kuperstokh.  New method of incorporating a body force term into the Lattice Boltzmann equation.  In Proceedings of the 5th International EDH Workshop, pages 241–246, Poitiers, France, 2004.
  * [50] Michele Buzzicotti and Patricio Clark Di Leoni.  Synchronizing subgrid scale models of turbulence to data.  Physics of Fluids, 32(12):125116, 2020.
  * [51] G. S. Patterson and S. A. Orszag.  Spectral Calculations of Isotropic Turbulence: Efficient Removal of Aliasing Interactions.  Phys. Fluids, 14:2538–2541, 1971.
  * [52] V. Borue and S. A. Orszag.  Self-similar decay of three-dimensional homogeneous turbulence with hyperviscosity.  Phys. Rev. E, 51:2859(R), 1995.
  * [53] Toshiyuki Gotoh and Takeshi Watanabe.  Statistics of transfer fluxes of the kinetic energy and scalar variance.  Journal of Turbulence, (6):N33, 2005.
  * [54] Michele Buzzicotti, Hussein Aluie, Luca Biferale, and Moritz Linkmann.  Energy transfer in turbulence under rotation.  Physical Review Fluids, 3(3):034802, 2018.
  * [55] A. N. Kolmogorov.  The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers.  In Dokl. Akad. Nauk SSSR, volume 30, pages 301–305. JSTOR, 1941\.
  * [56] T. Gotoh, D. Fukayama, and T. Nakano.  Velocity field statistics in homogeneous steady turbulence obtained using a high-resolution direct numerical simulation.  Physics of Fluids (1994-present), 14(3):1065–1081, 2002.
  * [57] Michael Sinhuber, Gregory P Bewley, and Eberhard Bodenschatz.  Dissipative effects on inertial-range statistics at high reynolds numbers.  Physical review letters, 119(13):134502, 2017.
  * [58] R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, F. Massaioli, and S. Succi.  Extended self-similarity in turbulent flows.  Phys. Rev. E, 48:R29–R32, 1993.
  * [59] Michele Buzzicotti, Brendan P Murray, Luca Biferale, and Miguel D Bustamante.  Phase and precession evolution in the burgers equation.  The European Physical Journal E, 39(3):1–9, 2016.
  * [60] Michele Buzzicotti, Akshay Bhatnagar, Luca Biferale, Alessandra S Lanotte, and Samriddhi Sankar Ray.  Lagrangian statistics for navier–stokes turbulence under fourier-mode reduction: fractal and homogeneous decimations.  New Journal of Physics, 18(11):113047, 2016.
  * [61] Alessandra S Lanotte, Roberto Benzi, Shiva K Malapaka, Federico Toschi, and Luca Biferale.  Turbulence on a fractal fourier set.  Physical review letters, 115(26):264502, 2015.


Corresponding author: J. Hubmayr<EMAIL_ADDRESS>
# Tolerance Analysis of Octave Bandwidth Millimeter-Wave Planar Orthomode
Transducer

Johannes Hubmayr  NIST Quantum Sensors Group, 325 Broadway Ave, Boulder, CO,
80305, USA  Jason E. Austermann  NIST Quantum Sensors Group, 325 Broadway Ave,
Boulder, CO, 80305, USA  James A. Beall  NIST Quantum Sensors Group, 325
Broadway Ave, Boulder, CO, 80305, USA  Jake A. Connors  NIST Quantum Sensors
Group, 325 Broadway Ave, Boulder, CO, 80305, USA  Shannon M. Duff  NIST
Quantum Sensors Group, 325 Broadway Ave, Boulder, CO, 80305, USA  Jeffrey J.
McMahon  Department of Astronomy and Astrophysics, University of Chicago, 5640
S. Ellis Ave, Chicago, IL 60637, USA

###### Abstract

Planar Orthomode Transducers (OMTs) are commonly used for polarization
measurements at millimeter wavelengths. We present an optical coupling study
of an octave bandwidth planar OMT in circular waveguide based on 3D
electromagnetic simulations. We quantify results through metrics such as co-
and cross- polar coupling, reflection, and waveguide leakage as a function of
the OMT construction geometry. We evaluate the tolerance of these metrics to
the waveguide backshort distance, probe impedance, waveguide gap size, and
waveguide-to-probe misalignment. Two probe geometries are studied: the
‘classic’ shape used in several previous experiments, and a new ‘wineglass’
geometry. The bandwidth ratio of both optimized OMTs is 2.0:1, defined where
co-polar coupling exceeds 80%. The average co-polar coupling, cross-polar
coupling, reflection, and waveguide leakage of the classic probe is
approximately 93%, $<$ -50 dB, 5% and 2%, respectively and depends slightly on
the exact frequency range. The wineglass probe co-polar coupling is $\sim$ 2%
larger. Radial waveguide misalignment at the level of 4% of the waveguide
radius can result in up to a 10% reduction in co-polar coupling and -20 dB
cross-polar coupling in one polarization. These results may be used to guide
the detector module designs of future Cosmic Microwave Background experiments
and beyond.

###### keywords:

planar OMT, orthomode transducer, CMB

##  1 INTRODUCTION

A traditional orthomode transducer (OMT) is a three-dimensional waveguide
component, which either separates or combines two orthogonally polarized
microwave signals. In a planar OMT of the type described herein, membrane
supported polarization-sensitive probes typically made of a superconducting
thin film are inserted into a waveguide cavity. Once on a planar structure,
subsequent on-chip circuitry conditions the signal for power sensing. Planar
OMTs of this type have found application in multiple Cosmic Microwave
Background (CMB) polarization experiments [1, 2, 3, 4, 5, 6] because they can
separately measure orthogonal linear polarization states within one spatial
pixel and can be fabricated into large arrays.

The optimization and geometric tolerance analysis of 30% fractional bandwidth
planar OMTs within circular waveguide operating near 150 GHz has previously
been studied [7]. Octave bandwidth planar OMTs demonstrated at microwave
frequencies [8] have been frequency-scaled and integrated to form multichroic
transition-edge-sensor (TES) polarimeters [9]. Detectors of this type have
been deployed in two generations of Atacama Cosmology Telescope (ACT)
receivers, spanning the center frequency range 27–230 GHz[10, 11]. Multiple
variations of this technology are either in production or in development for
upcoming instruments [12, 13, 14, 15].

Given their wide-spread adoption, here we present an optimization and
tolerancing study for octave bandwidth planar OMTs based on electromagnetic
simulations. As such this paper serves as an update to McMahon et al. 2009 [7]
for these broader bandwidth OMTs. Section 2 describes the physical model,
simulation parameters, and defines the optical coupling metrics. Section 3
presents the results of the optimized geometry and the variation of the
metrics to the most relevant geometric parameters. We include band-averaged
results within two bands since in most applications the bandwidth is
partitioned to form a dichroic, dual-polarization sensitive pixel. In Sec.
3.4, we present a new OMT probe geometry, which has been discussed in the
context of CMB-S4, and compare the performance to the ‘classic’ probe geometry
of Grimes et al. 2007 [8]. Discussion and conclusions are in Sec. 4 and Sec.
5, respectively.

##  2 MODEL DESCRIPTION

---
Figure 1:  Left: Simulation volume with TE11 mode superimposed over the input waveguide. Center: Expanded view of the lumped port used to terminate each of the four OMT probes. Right: Cross-section of the simulation model with labeled components (not to scale).  Parameter | Description | Value | Value/$a$
---|---|---|---
$a$ | waveguide radius | 1200 $\mu$m | 1
$a_{\mathrm{c}}$ | ground-plane radial width | 725 $\mu$m | 0.604
$z_{\mathrm{bs}}$ | backshort distance | 675 $\mu$m | 0.5625
$g_{1}$ | gap to input (top) waveguide | 15 $\mu$m | 0.0125
$g_{2}$ | gap to backshort (bottom) waveguide | 10 $\mu$m | 0.0083
$h_{\mathrm{p}}$ | probe height | 390 $\mu$m | 0.325
$\ell_{\mathrm{p1}}$ | probe length, rectangular section | 612 $\mu$m | 0.51
$\ell_{\mathrm{p2}}$ | probe length, triangular section | 268.8 $\mu$m | 0.224
$w_{\mathrm{cpw}}$ | CPW width | 2 $\mu$m |
$g_{\mathrm{cpw}}$ | CPW gap | 4.5 $\mu$m |
$t_{\mathrm{SiN_{x}}}$ | SiN${}_{\text{x}}$ thickness | 2.0 $\mu$m |
$t_{\mathrm{SiO_{2}}}$ | SiO2 thickness | 0.45 $\mu$m |
$t_{\mathrm{p}}$ | thickness of probe and ground plane | 0.2 $\mu$m |
$\lambda$ | Nb London penetration depth | 0.085 $\mu$m |
$\epsilon_{\mathrm{r,SiN_{x}}}$ | SiN${}_{\text{x}}$ relative dielectric constant | 6.8 |
$\epsilon_{\mathrm{r,SiO_{2}}}$ | SiO2 relative dielectric constant | 4.0 |
Table 1: Planar OMT model parameters separated into three sections:
geometrical parameters that scale with frequency proportional to the waveguide
radius, geometrical parameters that do not scale with frequency, and material
parameters.
---
Figure 2:  Geometry of ‘classic’ (left) and ‘wineglass’ (right) planar OMT
probes. Both red and blue areas are thin film Nb. The port numbers are defined
in the image on the right.

We use ANSYS HFSS111The use of this software is not an endorsement by the US
government. We list by name to be comprehensive. to optimize the planar OMT
and determine the tolerance to geometric variations. The physical model is
shown in Fig. 1 with labeled parameters that are defined and listed in Tab. 1.
An input waveguide of radius $a$ feeds four planar OMT probes of geometry
defined in Fig. 2 that lie a distance $z_{bs}$ away from a reflective
backshort. The waveguide walls are modeled as perfect conductors. Vertical
gaps $g_{1}$ and $g_{2}$ from the probes to the top and bottom waveguide
sections are necessary to insert the probes within the waveguide. For the
summed height of the gaps and at a radial distance $a+a_{c}$, we implement a
radiation boundary that allows the determination of the waveguide leakage, the
power which escapes the waveguide. The physical realization of this geometry
uses three separate silicon parts and is described in Ward et al. 2016[16].

We model the OMT probes as 2D sheets with an impedance boundary
$Z_{probe}=j\omega L$, where

$L=\mu_{0}\frac{\lambda}{\tanh{(t_{p}/\lambda)}}.$

$\omega$ is the angular frequency, $L$ is the kinetic inductance modeled with
London penetration depth $\lambda~{}=~{}85~{}$nm (suitable for thin film Nb),
and $t_{\mathrm{p}}$ is the thickness of the metal. Each probe output is a
coplanar waveguide (CPW) of width $w_{\mathrm{cpw}}$ and gap
$g_{\mathrm{cpw}}$, which is terminated in a lumped port of impedance
$Z_{\mathrm{term}}=90~{}\Omega$ that matches the CPW impedance. As compared to
waveports, the internal lumped ports provide a better representation of the
true geometry and avoid possible low-level reflection artifacts at the expense
of supporting only a single mode. The probes sit on a membrane comprised of a
layer of SiN${}_{\text{x}}$ and SiO2 of thickness $t_{\mathrm{SiN_{x}}}$ = 2
$\mu$m and $t_{\mathrm{SiO_{2}}}$ = 0.45 $\mu$m, respectively and relative
dielectric constant $\epsilon_{\mathrm{r,SiN_{x}}}$=6.8 and
$\epsilon_{\mathrm{r,SiO_{2}}}$=4.0, respectively. These choices match the
planar OMT fabrication process described in Duff 2016 [17]. We do not include
loss in the dielectrics to aid our determination of waveguide leakage, as we
define below.

We excite the input waveguide with eight modes, and align the first mode
(TE11) with the x-axis. For each parameter combination, we perform frequency
sweeps that span a typical 90/150 GHz dichroic OMT. However, since results may
be frequency-scaled as 1/$a$, we report both geometric parameters normalize to
the waveguide radius $a$ and the frequency axis normalized by the waveguide
cut-off frequency $f_{\text{c}}$. Thus, the frequency sweep is in the range
0.9 $<f/f_{\mathrm{c}}<~{}$2.38 with stepsize 0.054 $f/f_{c}$. The solution
frequency is 2.05 $f/f_{\mathrm{c}}$.

We largely use default HFSS parameters. The convergence criteria is standard,
that successive passes change the S-parameters by $<0.02$. The maximum number
of passes is set to 20. We use the direct solver and ‘first order’ basis
functions. All of our simulations achieve the defined convergence criteria.
Using a modest machine (4 GHz processor with 4 cores and 256 GB RAM), the
average simulation time is 30s per frequency point.

For simulations that require no waveguide misalignment, we make use of a
quarter-sized model with appropriate E- and H- symmetry planes to
substantially decrease the simulation time. Radial misalignment of the
waveguide is incompatible with these boundary conditions, and thus a full
model is used in these situations. We find good agreement between the full and
quarter model. For example the maximum difference in co-polarization is $<$
2%.

From the S-parameters of the eight modes of the input waveguide port (port 1)
to the four lumped ports which terminate each OMT probe (ports 2-5 defined
defined in Fig. 2), we define the co- and cross- polar power coupling (Co-pol
and X-pol), power reflection (Reflection), and waveguide leakage (Leakage) as
follows:

$\displaystyle\mathrm{Co-
pol}\overset{\mathrm{def}}{=}|S(2,1:1)|^{2}+|S(4,1:1)|^{2}$ (1)
$\displaystyle\mathrm{X-pol}\overset{\mathrm{def}}{=}|S(2,1:2)|^{2}+|S(4,1:2)|^{2}$
(2)
$\displaystyle\mathrm{Reflection}\overset{\mathrm{def}}{=}|S(1:1,1:1)|^{2}$
(3)
$\displaystyle\mathrm{Leakage}\overset{\mathrm{def}}{=}1-\sum_{i=1}^{8}|S(1:i,1:1)|^{2}-\sum_{j=2}^{5}|S(j,1:1)|^{2}$
(4)

Here the convention $S(i:n,j:m)$ denotes the voltage at port $i$ in mode $n$
divided by the voltage produced at port $j$ in mode $m$. Since the lumped
ports contain only a single mode, we drop the mode index for ports 2-5.

Coupling to the TE11 mode (mode $n$=1 in Eqns. 1-4) is desirable, whereas
coupling to higher order modes is not. This goal is represented by these
definitions. Simply put, the design goal is to construct a planar OMT that
maximizes Co-pol and minimizes X-pol, Reflection, and Leakage over the
broadest bandwidth achievable. As can be seen by these definitions, we sum the
power in opposite probes incoherently, which is accomplished in real devices
by ancillary on-chip circuitry. In the following section, we report these
metrics as a function of model geometry.

##  3 RESULTS

---
Figure 3:  Linear (left) and logarithmic (right) OMT coupling summary for the
‘classic’ probes. Co- (cross-) polar coupling peaks at 96% (-47 dB), and the
leakage out of the waveguide is 2%. The full bandwith ratio of the OMT is
2.0:1, defined where co-pol $>$ 0.8. For guidance, the two shaded regions
indicate fiducial bands in a dichroic OMT-coupled pixel, commonly used in
ground-based CMB instruments.

Figure 3 presents the coupling summary of the optimized planar OMT geometry.
Since results scale with the waveguide radius $a$, and this has been
explicitly confirmed in simulation, we report the frequency axis normalized by
$f_{c}$. In this manner, results are applicable to any specific octave of
bandwidth. The bandwidth ratio, defined where Co-pol$~{}>0.8$, is 2.0:1. For
guidance, the two shaded regions in Fig. 3 indicate fiducial bands (low band:
1.05$<f/f_{\mathrm{c}}<$1.45 and high band: 1.6$<f/f_{\mathrm{c}}<$2.2)
commonly used in dichroic OMT-coupled pixels. Band averaged metrics defined in
Eqns. 1-4 for the optimized geometry are reported in Tab. 2 along with these
values for the ‘wineglass’ probes discussed in Sec. 3.4.

###  3.1 Higher-Order Mode Coupling

---
Figure 4:  Coupling of higher order modes to the ‘classic’-probe OMT relative
to the TE11 mode. The TM01 (mode 3) and TE21 (mode 4) are the only parasitic
modes that couple stronger than -20 dB across the useable bandwidth.
---
Figure 5:  Effect of distance to reflective backshort $z_{\mathrm{bs}}$. Left:
Co-pol (solid) and Reflection (dashed) as a function of frequency for several
backshort distances listed in the legend as a fraction of the waveguide radius
$a$. Right: band-averaged co-polar coupling as a function of normalized
backshort distance for three different frequency ranges listed in the legend.
The vertical dashed-line indicates the fiducial value. A 25% change in
$z_{\mathrm{bs}}$ leads to only a 10% decrease in co-polar coupling,
demonstrating a broad optimum.

Figure 4 shows power coupling of the seven higher order modes relative to the
TE11 mode. The legend lists mode indices $n$, which correspond to the TE11
mode that is anti-aligned with the probes, TM01, 2 polarizations of TE21,
TE01, and 2 polarizations of TM11. For nearly all frequencies of interest
($f/f_{\mathrm{c}}<$ 2.2), only two modes couple at any appreciable level. The
TM01 (mode 3) coupling is $>$-10 dB for $f/f_{c}~{}>$1.4, plateaus to -3.5 dB
for 1.6 $<~{}f/f_{\mathrm{c}}<$ 2, and reaches 0 dB at $f/f_{\mathrm{c}}$ =
2.26. The TE21 polarization best aligned with the probes has $>-10$ dB
coupling when $f/f_{\mathrm{c}}~{}>$1.7, plateaus to -6 dB for 1.8
$<~{}f/f_{\mathrm{c}}<$ 2.1, and reaches 0 dB at $f/f_{\mathrm{c}}$ = 2.34.
All other modes couple a negligible amount. To reject the TM01 and TE21 modes,
on-chip microwave components have been implemented [18, 9], which exploit the
even/odd mode difference between higher order modes and TE11.

###  3.2 Optimal Parameter Values

In this subsection we discuss the change in metrics as a function of the most
relevant geometric parameters. We begin with the backshort distance
$z_{\mathrm{bs}}$, which is summarized in Fig. 5. The backshort distance
affects the bandwidth. As expected, larger backshort distances increase
coupling at lower frequencies at the expense of higher frequencies and vice
versa (Fig. 5 left). However, as illustrated in Fig. 5 right, a broad optimum
exists. For example in the low band, a 25% decrease from the optimum
$z_{\mathrm{bs}}$ results in only a 10% decrease in co-polar coupling.

---
Figure 6:  Effect of waveguide gaps $g_{1,2}$. The average co-polar coupling
(left) and Leakage (right) is plotted as a function of normalized total
waveguide gap, $(g_{1}+g_{2})/a$. The colors shown in the legend denote the
different parameter sweeps. The solid (dashed) lines show the average in the
low (high) frequency band.

As the waveguide gaps $g_{1,2}$ increase, the leakage radiation out of the
waveguide increases and in turn the co-polar coupling decreases. We determine
the co-polar and leakage radiation in three parameter sweeps. First, we vary
the top (bottom) waveguide while the bottom (top) waveguide has fixed
$g_{2}$=10 $\mu$m ($g_{1}$=15 $\mu$m). Modeling non-equal gaps is intentional,
reflects the actual assemblies, and is the result of using wafers with
standard thickness tolerance that are readily available. In the third case, we
set $g_{1}$ = $g_{2}$. We show the average Co-pol and Leakage within the low
and high bands in Fig. 6 as a function of normalized total waveguide gap,
$(g_{1}+g_{2})/a$. The Leakage and decrement in Co-pol are roughly linear with
the waveguide gap. We observe that leakage is more sensitive to gaps of the
top waveguide rather than the backshort waveguide, particularly at lower
frequencies. For the symmetric gap case, we see that Leakage $\sim$ 10% when
$g_{1}+g_{2}$ is 10% of the waveguide radius.

The membrane radius $a_{c}$ has no effect on the metrics to within simulation
uncertainty. Similarly the CPW impedance has wide tolerance. We vary
$g_{\mathrm{cpw}}$ from 3$\mu$m to 6$\mu$m. A 20% change in $g_{\mathrm{cpw}}$
leads to only a 1% change in Co-pol in both bands. Our fiducial value produces
a co-polar coupling within 1% of the maximum in either band. Lastly we note
that, other than affecting coupling to evanescent modes, the results are
independent of the input waveguide length, as expected. From these parameter
sweeps, we determine the optimal parameter values listed in Tab. 1.

###  3.3 Waveguide Misalignment

---
Figure 7:  Effect of waveguide misalignment. Results depend on the direction
of misalignment with respect to the TE11 polarization axis. The planar OMT is
more (less) affected by misalignment parallel (perpendicular) to the
polarization axis. The top row of images illustrate the misalignment geometry.
The left (right) column of plots show the optical coupling metrics as a
function of frequency for the parallel (perpendicular) misalignment case. The
legend, common to all plots, lists the radial misalignment normalized to the
waveguide radius $a$. The notation X-pol⟂ indicates that the cross-
polarization belongs to the OMT probes perpendicular to the TE11 polarization
axis. The results shown correspond to simulations that vary the bottom
waveguide section while the top section is aligned and static. The magnitude
of the results vary slightly for the other parameter sweep combinations
explored, but the general phenomenology holds.

The standard construction of the OMT uses three separate silicon parts: the
top waveguide, the detector wafer that contains the OMT probes, and the bottom
backshort waveguide. As such, misalignment of these parts is expected. In this
section we determine the tolerance to radial misalignment. We perform
parameter sweeps shifting the top and bottom waveguide sections independently,
together in the same direction, and in opposite directions.

The general effects are illustrated Fig. 7, which show that the OMT is more
sensitive to misalignment in the direction parallel to the TE11 polarization
axis. The performance degradation is an order of magnitude less for equal
radial shifts in the perpendicular direction. Waveguide misalignment results
in decreased sensitivity, differential gain, and increased cross-polarization.
An observed effect of lesser importance is the large reflection in a narrow
band near 2.08$f/f_{\mathrm{c}}$, associated with the turn-on of TE01 and TM11
modes, for shifts perpendicular to the polarization axis.

The band-averaged co- and cross- polar couplings are shown in Fig. 8.
Performance degradation is more sensitive to shifts of the backshort waveguide
than the top waveguide when the sections are displaced independently. Most
notably, the cross-polarization in the low frequency band is a factor of two
larger. However, the OMT is most sensitive to shifts of both sections in the
same direction. Shifts in opposite directions up to $dr/a$ = 0.04 fare better,
and as such we do not include these results in Fig. 8.

---
Figure 8:  Waveguide misalignment band-averaged co- and cross- polar coupling.
Left shows the band-averaged change in co-polar coupling from perfect
alignment as a function of the normalized waveguide displacement $dr/a$ in the
direction parallel to the polarization axis. Right shows the band-averaged
cross-polarization in probes orthogonal to the displacement direction as a
function of $dr/a$. Colors denote the different parameter sweep variations.
The legend indicates the displaced waveguide part. The ‘both’ parameter sweep
shifts both the top and bottom waveguide sections in the same direction. The
solid (dashed) line are the results for the low (high) frequency band.

Misalignment sweeps were run with larger waveguide gaps $g_{1,2}$, and the
fractional changes in the optical coupling metrics were the same as those for
the nominal gap sizes to within simulation accuracy.

###  3.4 ‘Wineglass’ Probes

An alternative to the classic probe geometry has been discussed in the context
of CMB-S4, and we refer to this design as the ‘wineglass probe’ due to its
shape (see Fig. 2 right). The wineglass probe shape is parameterized by the
following equation:

$\displaystyle y(x)=0.3a(bx-1)\tanh{(-\alpha x)}.$ (5)

Through independent parameters sweeps, we find that the values which maximize
co-polar coupling are $\alpha$ = 3.6 and $b$ = 0.72. However, a fairly broad
optimum exists. We also find that the optimal $z_{\mathrm{bs}}/a$ = 0.521,
slightly less than the distance in the classic probes. Figure 9 summarizes the
wineglass probe coupling, and Fig. 10 directly compares the co-polar coupling
to the classic probes. Table 2 lists the band averaged metrics for the
wineglass and classic probes within three frequency ranges.

The wineglass probes provide $\sim$ 2% higher average co-polar coupling. The
same waveguide misalignment study was carried out, and no appreciable
differences were discovered relative to the results of Sec. 3.3.

---
Figure 9:  Linear (left) and logarithmic (right) ‘wineglass’ probe OMT optical
coupling summary. The bandwidth ratio is identical to the ‘classic’-probe OMT.
---
Figure 10:  Direct comparison of the classic and wineglass planar OMT. The wineglass probe OMT Co-pol is on-average 2% higher than the classic probe OMT in the range 1.2$<$f/f${}_{c}<$2.25. | CLASSIC PROBE | WINEGLASS PROBE
---|---|---
| [1.05, 1.45] | [1.6, 2.2] | [1.2, 2.0] | [1.05, 1.45] | [1.6, 2.2] | [1.2, 2.0]
co-pol | 0.907 | 0.926 | 0.947 | 0.918 | 0.944 | 0.970
reflection | 0.072 | 0.045 | 0.033 | 0.058 | 0.028 | 0.008
leakage | 0.021 | 0.024 | 0.021 | 0.024 | 0.023 | 0.022
Table 2: Comparison of ‘classic’ to ‘wineglass’ probe band-averaged metrics
for several frequency bands. The frequency range is listed in the column
heading as a fraction of the waveguide cut-off frequency.

##  4 DISCUSSION

To vastly decrease the simulation volume, the excitation begins in circular
waveguide. The model does not include a feedhorn attached to the waveguide. A
limited set of simulations including a feedhorn profile have been carried out,
and there is good agreement for the co-polar coupling. The main difference is
that cross-polarization reaches $\sim$ -20 dB, which is expected from the horn
design [19].

On average 2% of the power in one polarization leaks outside of the waveguide.
While small, this power is uncontrolled and may lead to non-idealities. In our
experience from ground-based CMB observations with planar OMTs, waveguide
leakage has not been problematic. Nonetheless, opportunities to improve in
this aspect of the design exist. Approaches to mitigate waveguide leakage
include the addition of microwave absorber, inclusion of a waveguide choke,
metallic bump bonds to form a seal around the waveguide [20], and passing the
CPW probes through minimum-sized waveguide wall holes [21].

The wineglass OMT probes provide modest gain in co-polar coupling, but we note
that the reduced reflection may also help limit image ‘ghosting’, which can be
of concern for large, flat focal planes. Both the classic and wineglass planar
OMTs have the same tolerance to waveguide misalignment. A good rule of thumb
is to keep any radial waveguide misalignment to the probes less than 4% of the
waveguide radius. This criteria will ensure $<$10% co-polar coupling reduction
and that cross-polarization remains $<$-20 dB. The top-waveguide misalignment
study presented may be a useful reference when mating the planar OMT array to
metal feedhorns where differential thermal contraction is of concern.

##  5 CONCLUSION

We have presented a detailed simulation study of the optical coupling
performance of octave bandwidth planar OMTs. Prior to this work, relevant
results were limited to OMTs of 30% fractional bandwidth. For the designs
presented, the bandwidth ratio defined where co-polar coupling exceeds 80% is
2.0:1. The average co-polar coupling, cross-polar coupling, reflection, and
waveguide leakage is approximately 93%, $<$ 50 dB, 5% and 2%, respectively and
depends slightly on the exact frequency range. The ‘wineglass’ probes increase
co-polar coupling and decrease reflection by $\sim$ 2% as compared to the
‘classic’ probes. A good rule of thumb to avoid strong performance degradation
is to align the top and bottom waveguides to the OMT probes to better than 4%
of the waveguide radius. We hope these results provide general guidance in the
development of future focal planes for CMB measurements and beyond.

###### Acknowledgements.

Support for this work was provided in part by the NASA APRA program, proposal
#17-APRA17-0045 .

## References

  * [1] Essinger-Hileman, T., et al. “The Atacama B-mode Search: CMB polarimetry with transition-edge-sensor bolometers,” AIP conference proceedings 1185(1), 494–497, (2009).
  * [2] Stanchfield, S., et al., “Development of a microwave SQUID-multiplexed TES array for MUSTANG-2,” JLTP 184(1), 460–465 (2016).
  * [3] Austermann, J. E., et al., “SPTpol: an instrument for CMB polarization measurements with the South Pole Telescope,” SPIE proceedings 8452 (Sept. 2012).
  * [4] Niemack, M. D., et al., “ACTPol: a polarization-sensitive receiver for the Atacama Cosmology Telescope,” SPIE proceedings 7741 (July 2010).
  * [5] Hubmayr, J., et al., “Design of 280 GHz feedhorn-coupled TES arrays for the balloon-borne polarimeter SPIDER,” SPIE proceedings 9914, 99140V, (2016).
  * [6] Essinger-Hileman,T., et al., “CLASS: the cosmology large angular scale surveyor,” SPIE proceedings 9153, 91531I, (2014).
  * [7] McMahon, J., AIP Conference Proceedings, 1185(1), 490–493, American Institute of Physics (2009).
  * [8] Grimes, P., King, O., Yassin, G., and Jones, M., “Compact broadband planar orthomode transducer,” Electronics Letters 43(21), 1146–1148 (2007).
  * [9] McMahon, J., et al., “Multi-chroic feed-horn coupled TES polarimeters,” Journal of Low Temperature Physics 167(5-6), 879–884 (2012).
  * [10] Thornton, R., et al., “The Atacama Cosmology Telescope: the polarization-sensitive ACTpol instrument,” The Astrophysical Journal Supplement Series 227(2), 21 (2016).
  * [11] Henderson, S., et al., “Advanced ACTPol cryogenic detector arrays and readout,” JLTP 184(3), 772–779 (2016).
  * [12] Ade, P., et al., “The Simons Observatory: science goals and forecasts,” JCAP 2019(02), 056 (2019).
  * [13] Li, H., et al., “Probing primordial gravitational waves: Ali cmb polarization telescope,” National Science Review 6(1), 145–154 (2019).
  * [14] Abazajian, K., et al., “CMB-S4 science case, reference design, and project plan,” arXiv preprint arXiv:1907.04473 (2019).
  * [15] Hazumi, M., et al., “LiteBIRD satellite: JAXA’s new strategic L-class mission for all-sky surveys of cosmic microwave background polarization,” SPIE proceedings 11443, 431 – 450, (2020).
  * [16] Ward, J. T., et al., “Mechanical designs and development of tes bolometer detector arrays for the advanced actpol experiment,” SPIE proceedings 9914, 794–802, (2016).
  * [17] Duff, S. M., et al., “Advanced ACTPol multichroic polarimeter array fabrication process for 150 mm wafers,” JLTP 184(3), 634–641 (2016).
  * [18] U-yen, K., Wollack, E. J., Papapolymerou, J., and Laskar, J., “A broadband planar magic-t using microstrip slotline transitions,” IEEE Transactions on Microwave Theory and Techniques 56(1), 172–177 (2008).
  * [19] Simon, S. M., et al., “Feedhorn development and scalability for Simons Observatory and beyond,” SPIE proceedings 10708, 787–798, (2018).
  * [20] Rostem, K., et al., “Silicon-based antenna-coupled polarization-sensitive millimeter-wave bolometer arrays for cosmic microwave background instruments,” SPIE proceedings 9914, 54 – 63, (2016).
  * [21] Gualtieri, R., Barry, P. S., Cecil, T., Bender, A. N., Chang, C. L., Hood, J. C., Lisovenko, M., and Yefremenko, V. G., “Optical leakage mitigation in ortho-mode transducer detectors for microwave applications,” JLTP, 1–8 (2022).