# Magnetoconvection in a horizontal duct flow at very high Hartmann and
Grashof numbers

R. Akhmedagaev1  O. Zikanov1<EMAIL_ADDRESS>Y. Listratov2,3  1University
of Michigan - Dearborn, 4901 Evergreen Road, Dearborn, MI 48128-1491, USA
2Moscow Power Engineering Institute, 14 Krasnokazarmennaya Street, Moscow,
111250, Russia 3Joint Institute for High Temperatures Russian Academy of
Science, Izhorskaya Street 13, Building 2, Moscow, 125412, Russia

###### Abstract

Direct numerical simulations and linear stability analysis are carried out to
study mixed convection in a horizontal duct with constant-rate heating applied
at the bottom and imposed transverse horizontal magnetic field. A two-
dimensional approximation corresponding to the asymptotic limit of very strong
magnetic field effect is validated and applied, together with full three-
dimensional analysis, to investigate the flow’s behaviour in the previously
unexplored range of control parameters corresponding to typical conditions of
a liquid metal blanket of a nuclear fusion reactor. It is found that the
instability to quasi-two-dimensional rolls parallel to the magnetic field
discovered at smaller Hartmann and Grashof numbers in earlier studies also
occurs in this parameter range. Transport of the rolls by the mean flow leads
to magnetoconvective temperature fluctuations of exceptionally high
amplitudes. It is also demonstrated that quasi-two-dimensional structure of
flows at very high Hartmann numbers does not guarantee accuracy of the
classical two-dimensional approximation. The accuracy deteriorates at the
highest Grashof numbers considered in the study.

##  1 Introduction

Combined convection and magnetohydrodynamic (MHD) effects dramatically change
the nature of flows of electrically conducting fluids. The combination appears
in many technological applications such as metallurgy, liquid metal batteries,
and growth of semiconductor crystals (Ozoe, 2005; Davidson, 2016). Another
prominent example is the liquid metal blankets of nuclear fusion reactors
where an electrically conducting fluid (e.g., a PbLi alloy) serves as a
coolant, radiation shield and tritium breeder (Abdou et al., 2015). A
distinctive feature of this system is that the convection and magnetic filed
effects are both exceptionally strong.

Many aspects of the transformation of flows of electrically conducting fluid
under the influence of a strong magnetic field, such as suppression of
turbulent fluctuations, anisotropic or quasi-two-dimensional (Q2D) states with
zero or weak velocity gradients along the field lines, formation of MHD
boundary layers, and delay of laminar-turbulent transition, are relatively
well understood (see, e.g., Branover, 1978; Davidson, 2016; Sommeria & Moreau,
1982; Zikanov et al., 2014). This paper addresses a recently discovered and
still poorly understood phenomenon - the high amplitude fluctuations in flows
in ducts and pipes (see, e.g., Genin et al., 2011; Vetcha et al., 2013;
Zikanov et al., 2013; Belyaev et al., 2021). The term magneto-convective
fluctuations (MCFs) proposed for the phenomenon by Belyaev et al. (2021) will
be used in this paper. As discussed in detail in the review of Zikanov et al.
(2021) and references therein, the fluctuations have been detected in
experimental and computational studies of a large variety of systems: pipes
and ducts of various orientations with respect to gravity, various heating
arrangements, and various configurations of the magnetic field.

The nature of the magnetoconvective fluctuations can be briefly described as
follows. They appear in the conditions of a very strong magnetic field effect,
i.e. in the range of Hartmann numbers, where turbulence is fully suppressed by
magnetic damping. In experiments, the MCFs are manifested by oscillations of
temperature with very high amplitude (up to $50K$ in some cases) and typical
frequencies much lower than the frequencies of turbulence-induced
fluctuations. Specific properties of the MCFs vary with the flow’s
configuration and values of the control parameters (Zikanov et al., 2021). The
effect has potentially serious consequences for design and operation of liquid
metal blankets of future fusion reactors. Should the fluctuations appear in an
actual blanket, they will lead to strong and unsteady thermal stresses in the
walls. Due to their possibly very large amplitude, the stresses will threaten
the structural integrity of a fusion reactor system. Significant effects on
heat transfer, transport of tritium, and wall corrosion are also anticipated.
As we discuss later in this section, it is yet impossible to say how realistic
these expectations are, since no experiments or computations at very high $Ha$
and $Gr$ typical for reactor conditions have been conducted so far.

Flows in a rectangular duct with heating applied at the bottom and imposed
transverse horizontal magnetic field (see figure 1) are considered in this
paper. The configuration is relevant to design of blankets and divertors of
fusion reactors (Kirillov & Muraviev, 1997). It is also important as an
archetypal system, in which the MCFs were first identified (in Genin et al.
(2011) and Zikanov et al. (2013), where they were named anomalous
fluctuations) and explained.

$Z$$Y$$X$${\boldsymbol{U}}$$q$$\boldsymbol{B}$$\boldsymbol{g}$$h$$2d$

Figure 1:  Flow geometry and coordinate system. The arrows marked by letters
$\boldsymbol{g}$, $\boldsymbol{B}$ and $q$ denote, respectively, the
orientations of the gravity acceleration, magnetic field and wall heating.

Similar systems for either ducts or round pipes have been studied
experimentally (Genin et al., 2011; Belyaev et al., 2015; Sahu et al., 2020)
and numerically (Zikanov et al., 2013; Zhang & Zikanov, 2014; Vo et al., 2017;
Listratov et al., 2018). The flow is controlled by four dimensionless
parameters: the Reynolds, Prandtl, Grashof and Hartmann numbers,

${\it Re}=\frac{Ud}{\nu},\ \ \ {\it Pr}=\frac{\nu}{\chi},\ \ \ {\it
Gr}=\frac{g\beta qd^{4}}{\nu^{2}\kappa},\ \ \ {\it
Ha}=Bd\sqrt{{\frac{\sigma}{\rho\nu}}},$ (1)

with the duct half-width $d$, the mean streamwise velocity $U$, the kinematic
viscosity $\nu$, the temperature diffusivity $\chi$, the acceleration due to
gravity $g$, the coefficient of thermal expansion $\beta$, the heat flux of
constant rate $q$, the thermal conductivity $\kappa$, the electrical
conductivity $\sigma$ and the mass density $\rho$. Rectangular duct geometry
adds the aspect ratio $\Gamma={2d}/{h}$ as a parameter, where $h$ is the
height of the duct.

Interesting results were obtained in the linear stability analysis of the
Poiseulle-Rayleigh-Bérnard duct flow with a transverse magnetic field
performed by Vo et al. (2017). Two-dimensional approximation valid in the
limit of strong magnetic field presented later in this paper was used. One
important result of Vo et al. (2017) is relevant to our work even though
different boundary conditions were used. It was demonstrated that the
convection instability occurs at moderate and high Grashof number
(approximately above $10^{6}$) at the Hartmann numbers ($\sim 10^{4}$) typical
for reactor blanket conditions.

The presence of MCFs in a horizontal round pipe with a lower half of the wall
heated was detected in experiments (Genin et al., 2011; Belyaev et al., 2015)
and explained in the linear stability analysis and direct numerical
simulations (DNS) by Zikanov et al. (2013). Flows of mercury with ${\it
Pr}\approx 0.022$, ${\it Re}$ up to $10^{5}$, ${\it Gr}$ up to $10^{8}$, and
${\it Ha}$ up to $500$ were investigated. It was shown that at a strong
magnetic field the suppression of flow structures having large gradients along
the field lines resulted in the most unstable modes in the form of convection
rolls with axes aligned with the field. The instability led to development of
convection structures in the form of Q2D rolls. Transport of the rolls by the
mean flow generated the MCFs.

The analysis was extended in the numerical simulations of Zhang & Zikanov
(2014). Flows in a horizontal duct of aspect ratio $\Gamma=1$ with bottom
heating and transverse magnetic field at ${\it Pr}=0.0321$, ${\it Re}=5000$,
$50\leq{\it Ha}\leq 800$ and $10^{5}\leq{\it Gr}\leq 10^{9}$ were
investigated. The instability leading to the formation of Q2D rolls similar to
those found in the pipe flow was detected at sufficiently high $Gr$ and $Ha$.

Investigations of Zhang & Zikanov (2014) conducted in the broader range of
parameters than for the pipe flow demonstrated existence of two distinct
secondary flow regimes. The realization of the regimes depended on the
relative strength of the convection and MHD effects. The low-${\it Gr}$ type
characterized by Q2D distributions of velocity and temperature dominated by
spanwise rolls appeared at $Gr$ below a certain $Gr^{\ast}(Ha)$. At higher
${\it Gr}$, stronger convection resulted in three-dimensional (3D) flow states
combining the spanwise rolls with streamwise ones (the geometrically preferred
convection structure in pipes and ducts with bottom heating).

Flows of liquid metals in fusion reactor blankets and divertors are subject to
very strong effects of convection ($Gr$ $\sim 10^{10}-10^{12}$) and magnetic
fields ($Ha\sim 10^{4}$) (see, e.g., Smolentsev et al., 2008, 2010). Such
extreme parameters present serious obstacles to analysis, because neither
laboratory experiments nor 3D simulations of unsteady flow regimes in
realistic blanket or divertor geometries can, at this moment, achieve such
values.

In an attempt to reach the typical blanket flow conditions, the data on two
types of the secondary flow regime in a horizontal duct were extrapolated to
high $Gr$ and $Ha$ by Zhang & Zikanov (2014). The extrapolation predicted
existence of MCFs at the typical blanket parameters. It also predicted that
the flow would likely be of the low-${\it Gr}$ type at ${\it Gr}\leq 10^{10}$
and of the high-${\it Gr}$ type at higher ${\it Gr}$. The experiments in the
pipe flow (see the review of recent results in Zikanov et al., 2021), on the
contrary, indicate that MCFs may disappear at high ${\it Ha}$, so the
extrapolation can be wrong. The nature of the convection flow at the
parameters corresponding to ducts in blankets and divertors of an operating
fusion reactor remains unknown, setting up the motivation for the present
study.

The focus of our investigation is on the magnetoconvection in the range of
very high $Gr$ and $Ha$ including the values typical for a reactor blanket and
divertor. To the best of our knowledge, this study is the first to analyze the
MCF effect in this range. Linear stability analysis and DNS of flows in a
horizontal duct with ${\it Pr}=0.025$, $Re=5000$, $10^{8}\leq{\it Gr}\leq
10^{10}$ and $10^{3}\leq{\it Ha}\leq 10^{4}$ are performed. The study follows
the work of Zhang & Zikanov (2014) but differs by much larger values of $Gr$
and $Ha$ and the aspect ratio $\Gamma=3.5$ selected to match the new
experimental facility (see, e.g., Belyaev et al., 2017), on which the same
configuration is to be explored at ${\it Ha}\lesssim 10^{3}$ and ${\it
Gr}\lesssim 10^{8}$ in the near future. Another essential difference between
our work and the work by Zhang & Zikanov (2014) is that we carry out an in-
depth analysis of the accuracy of the two-dimensional approximation applied to
Q2D flows at such high $Ha$.

##  2 Presentation of the problem

The flow of an incompressible, Newtonian, viscous, electrically conducting
fluid (a liquid metal) with constant physical properties is considered. The
fluid moves through a horizontal duct of aspect ratio $\Gamma=3.5$ (see figure
1). Spatially uniform and time-independent magnetic field
$\boldsymbol{B}=B\boldsymbol{e}_{y}$ is imposed in the horizontal transverse
direction. All walls are perfectly electrically insulated. The top and side
walls are perfectly thermally insulated. The bottom wall is subject to uniform
heating with the heat flux of constant rate $q$. The no-slip boundary
conditions for velocity are applied at the walls.

###  2.1 Physical model

The Boussinesq and quasi-static approximations are applied. The quasi-static
approximation is valid at small Reynolds and Prandtl numbers and usually
utilized in numerical, experimental and theoretical studies of MHD flows of
liquid metals (Davidson, 2016). The approximation implies that the imposed
magnetic field $\boldsymbol{B}$ is much stronger than the perturbations of the
magnetic field $\boldsymbol{b}$ induced by the electric currents caused by the
fluid motion. The induced magnetic field can be neglected in the expressions
of the Lorentz force and Ohm’s law. Furthermore, the induced field is assumed
to adjust instantaneously to changes of velocity field.

The governing equations are rendered non-dimensional using the duct half-width
in the magnetic field direction $d$ as the length scale, mean streamwise
velocity $U$ as the velocity scale, wall heating-based group $qd/\kappa$ as
the temperature scale, $B$ as the scale of the magnetic field strength, and
$dUB$ as the scale of electric potential. The equations can be written as:

$\frac{\partial\boldsymbol{u}}{\partial
t}+(\boldsymbol{u}\bcdot\bnabla)\boldsymbol{u}=-\bnabla
p-\bnabla\hat{p}-\bnabla\tilde{p}+\frac{1}{Re}\bnabla^{2}\boldsymbol{u}+{\it\boldsymbol{F}_{b}}+{\it\boldsymbol{F}_{L}},$
(2) $\bnabla\bcdot\boldsymbol{u}=0,$ (3) $\frac{\partial\theta}{\partial
t}+\boldsymbol{u}\bcdot\bnabla\theta={\frac{1}{{\it Re}{\it
Pr}}}\bnabla^{2}\theta-u_{x}\frac{\text{d}{\it T_{m}}}{\text{d}x},$ (4)

where $\boldsymbol{u}$ is the velocity field. The decompositions of the
temperature and pressure fields commonly used in studies of mixed convection
in ducts and pipes (see, e.g., Alboussière et al., 1993; Lyubimova et al.,
2009; Zikanov et al., 2013; Zhang & Zikanov, 2014; Zikanov et al., 2021) are
applied. The decompositions are convenient, since they allow one to recast the
problem in terms of the fluctuation fields, which are statistically uniform in
the streamwise direction and, thus, study the flow in a relatively short
segment of the channel with periodic inlet-exit conditions. The temperature
field is written as a sum

$T(\boldsymbol{x},t)=\it{T_{m}}(x)+\theta(\boldsymbol{x},t)$ (5)

of fluctuations $\theta$ and the mean-mixed temperature
$\it{T_{m}}(x)=A^{-1}\int_{A}u_{x}T\text{d}A$, where $A$ is the cross-section
area of the duct. One can also use the decomposition into fluctuations and
simple mean temperature $\bar{T}(x)=A^{-1}\int_{A}T\text{d}A$. Applying the
energy balance between the wall heating and the streamwise convection heat
transfer, we find that $\it{T_{m}}(x)$ and $\bar{T}(x)$ are linear functions
with the same derivative:

$\frac{\text{d}{\it
T_{m}}}{\text{d}x}=\frac{\text{d}{\it\bar{T}}}{\text{d}x}=\frac{\Pi}{APe}=\frac{\Gamma}{2Pe},$
(6)

where $\Pi$ is the perimeter of the heated portion of the wall.

The total pressure $P$ is presented in (2) as

$P=\hat{p}(x)+\tilde{p}(x,z)+p(\boldsymbol{x},t),$ (7)

where $p(\boldsymbol{x},t)$ is the field of pressure fluctuations
statistically homogeneous in the streamwise direction, and $\hat{p}$ is a
linear function of $x$ corresponding to the spatially uniform streamwise
gradient ${\text{d}{\hat{p}}}/{\text{d}x}$ applied as a flow-driving
mechanism. In the simulations discussed in this paper, the gradient is
adjusted at every time step to maintain constant mean velocity.

The third term of the decomposition becomes necessary in numerical models of
mixed convection in non-vertical channels with periodic inlet–exit conditions.
The component

$\tilde{p}(x,z)=\frac{\text{d}{\it
T_{m}}}{\text{d}x}\frac{Gr}{Re^{2}}xz=\frac{\Pi}{A{\it Re}{\it
Pr}}\frac{Gr}{Re^{2}}xz,$ (8)

arises due to the buoyancy force caused by the mean–mixed temperature ${\it
T_{m}}$:

$\boldsymbol{F}_{b,m}=GrRe^{-2}\boldsymbol{e}_{z}T_{m},$ (9)

where $\boldsymbol{e}_{z}$ is the unit vector opposite to the direction of
gravity (see figure 1). The force has a non-zero curl and, therefore, modifies
the velocity field. Its action on the flow can be described by introducing the
pressure field $\tilde{p}$, such that its vertical gradient balances
$\boldsymbol{F}_{b,m}$. The pressure field is a two-dimensional function
increasing with the streamwise coordinate $x$ and vertical coordinate $z$. Its
$z$-dependent $x$-gradient, which appears in the respective momentum equation,
generates a flow in the positive $x$-direction in the lower part of the
channel and in the negative $x$-direction in the upper part. The result is a
top-bottom asymmetry of the streamwise velocity profile and of the associated
convection heat flux, which can dramatically change the structure of the flow
at high $Gr$ and $Ha$ (see Zikanov et al., 2013; Zhang & Zikanov, 2014, 2017;
Zikanov et al., 2021).

The buoyancy force in (2) is

$\boldsymbol{F}_{b}=GrRe^{-2}\boldsymbol{e}_{z}T.$ (10)

The Lorentz force is computed as

$\boldsymbol{F}_{L}=Ha^{2}Re^{-1}\boldsymbol{j}\times\boldsymbol{e}_{y},$ (11)

where $\boldsymbol{e}_{y}$ is the unit vector along the imposed magnetic field
(see figure 1). The electric current $\boldsymbol{j}$ is determined by the
Ohm’s law

$\boldsymbol{j}=-\nabla\phi+(\boldsymbol{u}\times\boldsymbol{e}_{y}),$ (12)

where the electric potential $\phi$ is a solution of the Poisson equation
expressing the instantaneous electric neutrality of the fluid:

$\nabla^{2}\phi=\nabla\cdot(\boldsymbol{u}\times\boldsymbol{e}_{y}).$ (13)

The inlet–exit conditions are those of periodicity of the velocity
$\boldsymbol{u}$, temperature fluctuations $\theta$, pressure fluctuations
$p$, and potential $\phi$.

###  2.2 Two-dimensional approximation

Flows with a very strong imposed magnetic field are considered, so the
Hartmann number and the Stuart number satisfy ${Ha}\gg 1$ and $N\equiv
Ha^{2}/Re\gg 1$, respectively. The flows are anticipated to have Q2D form with
nearly zero gradients along the magnetic field lines except in the thin
Hartmann layers at the walls perpendicular to the field. The 2D approximation
proposed by Sommeria & Moreau (1982) can be applied in this asymptotic limit.
The problem can be expressed in terms of the variables integrated wall-to-wall
along the direction of the magnetic field, leading to 2D dynamics for
y-averaged quantities. The approximation has been verified and examined by
Pothérat et al. (2000, 2005), and utilized in numerical studies of liquid
metal flows in rectangular ducts (see, e.g., Pothérat, 2007; Smolentsev et
al., 2012; Vetcha et al., 2013; Vo et al., 2017; Zhang & Zikanov, 2018).

The often applied abbreviation SM82 will be used for the model in the
following. The $y-$independent solutions obtained in the framework of the
model will be referred to as 2D solutions, while the full solutions obtained
numerically without resorting to the model will be designated as 3D.

The SM82 model is derived for flows with ${Ha}\gg 1$ and $N\gg 1$, in domains
with electrically insulating walls and constant wall-to-wall distance in the
field direction. It utilizes the fact that the Lorentz force becomes nearly
zero in Q2D flows in such geometries, and the effect of the magnetic field on
the flow is largely reduced to the viscous friction in the thin Hartmann
layers, which can be accurately modeled by the linear friction term in the
momentum equation. In non-isothermal flows, the model can be extended to 2D
approximation of temperature if the imposed heat flux is perpendicular to the
magnetic field (see, e.g., Gelfgat & Molokov, 2011; Vetcha et al., 2013; Zhang
& Zikanov, 2018).

The SM82 version of (2)$-$(4) is

$\frac{\partial\boldsymbol{u}}{\partial
t}+(\boldsymbol{u}\bcdot\bnabla)\boldsymbol{u}=-\bnabla
p-\bnabla\hat{p}-\bnabla\tilde{p}+\frac{1}{Re}\bnabla^{2}\boldsymbol{u}-\frac{Ha}{Re}\boldsymbol{u}+\frac{Gr}{Re^{2}}T\hat{\boldsymbol{e}}_{z},$
(14) $\bnabla\bcdot\boldsymbol{u}=0,$ (15) $\frac{\partial\theta}{\partial
t}+\boldsymbol{u}\bcdot\bnabla\theta={\frac{1}{{\it Re}{\it
Pr}}}\bnabla^{2}\theta-u_{x}\frac{\text{d}{\it T_{m}}}{\text{d}x},$ (16)

where all the flow variables are now 2D fields obtained by wall-to-wall
averaging. The term (${Ha}/{Re}$)$\boldsymbol{u}$ represents the effect of
friction in the Hartmann layers.

###  2.3 Numerical method

The governing equations (2)$-$(4) and (14)$-$(16) are solved numerically using
the finite difference scheme introduced by Krasnov et al. (2011), and later
developed, tested and applied to high-${Ha}$ flows and flows with thermal
convection in numerous works including those by Krasnov et al. (2012), Zhao &
Zikanov (2012), Zikanov et al. (2013), Zhang & Zikanov (2014), and Gelfgat &
Zikanov (2018). The spatial discretization is of the second order and nearly
fully conservative with regards to the mass, momentum, electric charge,
kinetic energy and thermal energy conservation principles (Ni et al., 2007;
Krasnov et al., 2011). The computational grid is clustered towards the walls
according to the coordinate transformation in the horizontal direction
$y=\tanh(A_{y}\eta)/\tanh(A_{y})$ and in the vertical direction
$z=\tanh(A_{z}\xi)/\tanh(A_{z})$. Here $\eta$ and $\xi$ are the transformed
coordinates, in which the grid is uniform, and $A_{y}$ and $A_{z}$ are the
coefficients determining the degrees of clustering. The time discretization is
implicit for the conduction and viscosity terms and based on the
Adams–Bashforth/backward-differentiation method of the second order and the
standard projection algorithm (see, e.g., Zikanov, 2019). The nonlinear
convection and body force terms are treated explicitly. The elliptic equations
for potential, pressure, temperature and velocity components are solved using
the Fourier decomposition in the streamwise coordinate and the direct cyclic
reduction solution of the 2D equations for Fourier components conducted on the
transformed grid (see Krasnov et al., 2011).

The algorithm is parallelized using the hybrid MPI-OpenMP approach. The MPI
memory distribution is along the $y-$coordinate in the physical space and
along the streamwise wavenumber in the Fourier space.

###  2.4 Approach to linear stability analysis

The base flow needs to be selected before conducting the linear stability
analysis. We note that an archetypal structure of a laminar flow with
convection in a horizontal channel heated from below is a superposition of the
streamwise flow $u_{x}(y,z)$ and one or several streamwise-uniform convection
rolls $(u_{y}(y,z),u_{z}(y,z))$. At high $Ha$, the structure can be modified
by the magnetic field and replaced by a 3D structure with the rolls aligned
with the magnetic field and, thus, $x$-dependent velocity and temperature at
high $Ha$. Following Zikanov et al. (2013) and Zhang & Zikanov (2014), we
treat the problem as that of the instability of the laminar steady-state
streamwise-uniform base flow $\boldsymbol{U}(y,z)$, $\Theta(y,z)$, $P(y,z)$ to
$x$-dependent perturbations.

The base flow is calculated by artificially imposing uniformity in the
streamwise direction, i.e. by applying $x$-averaging after every time step. In
order to assure that a fully developed state of the base flow is reached, each
solution is computed for sufficiently long time. Long evolution, in some
regimes up to $1000$ time units, is typically required in order to arrive at
this state. No unsteady base flow solutions have been detected in the studied
range of parameters. The steady-state solutions are discussed in section 3.1.

The linear stability analysis is conducted using a modified version of the
numerical model described in section 2.3. We follow evolution of perturbations
- solutions of the equations linearized around the base flow
$\boldsymbol{U}(y,z)$, $\Theta(y,z)$, $P(y,z)$. The evolution of individual
Fourier modes determined by their streamwise wavelengths $\lambda$ is
computed.

The linear instability is identified by the exponential growth of the
perturbations with the growth rate determined as

$\gamma=\frac{1}{2E^{\prime}}\frac{dE^{\prime}}{dt},$ (17)

where $E^{\prime}=\langle f^{2}\rangle$, $\langle\ldots\rangle$ stands for
volume averaging, and $f$ stands for perturbations of a velocity component or
temperature. The growth rate coefficient is recorded after its values computed
for all three velocity components and temperature coincide with each other and
remain constant within the third digit after the decimal point for at least
$100$ time units. The results of the linear stability analysis are presented
in sections 3.3 and 3.4.

###  2.5 Grid sensitivity study

The grid sensitivity study has been conducted for the base flow. A detailed
description of the various flow regimes is provided in section 3.1. For the
present discussion, it is sufficient to say that accurate resolution of the
internal flow structure, along with two boundary layers: the Hartmann layers
of thickness $\delta_{Ha}\sim{\it Ha}^{-1}$ at the vertical walls and the
Shercliff layers of thickness $\delta_{Sh}\sim{\it Ha}^{-1/2}$ at the top and
bottom walls, is critically important for accurate representation of the flow
behavior.

As an example, the results obtained at $Ha=1200$, $Gr=10^{8}$ are presented in
table 1. In a fully developed steady-state flow, the integrated Lorentz and
buoyancy forces are zero. The wall friction must be balanced by the driving
pressure gradient according to

$\frac{d\hat{p}}{dx}=A^{-1}\left(\tau_{Ha}+\tau_{Sh}\right),$ (18)

where $\tau_{Ha}$ and $\tau_{Sh}$ are the computed values of the integrated
friction forces at the Hartmann and Shercliff walls of the duct, respectively,
expressed as

$\tau_{Ha}=\tau_{y}=\frac{1}{Re}\sum_{y=\pm
1}\int_{-2/\Gamma}^{2/\Gamma}\frac{\partial U_{x}}{\partial y}dz,\ \ \ \
\tau_{Sh}=\tau_{z}=\frac{1}{Re}\sum_{z=\pm
2/\Gamma}\int_{-1}^{1}\frac{\partial U_{x}}{\partial z}dy.$ (19)

Values of $\tau_{Ha}$, $\tau_{Sh}$ and the error $\epsilon$, with which the
computed solution satisfies (18), found on various grids are compared in table
1. On the basis of these data, we conclude that the grid with $N_{y}\times
N_{z}=192\times 96$, $A_{y}=4.0$ and $A_{z}=2.0$ is sufficient. The maximum
and minimum grid steps of such a grid are $\Delta y_{min}\approx 0.0001$,
$\Delta y_{max}\approx 0.042$, $\Delta z_{min}\approx 0.0009$, $\Delta
z_{max}\approx 0.012$. The Hartmann and Shercliff layers are resolved by,
respectively, $9$ and $32$ grid points.

$N_{y}$ | $N_{z}$ | $A_{y}$ | $A_{z}$ | $-\tau_{Ha}$ | $-\tau_{Sh}$ | $\epsilon$ | $N_{Ha}$ | $N_{Sh}$
---|---|---|---|---|---|---|---|---
128 | 32 | 4.0 | 2.0 | 0.23425 | 0.02914 | 0.00074 | 6 | 10
128 | 64 | 4.0 | 2.0 | 0.23460 | 0.26417 | 0.00070 | 6 | 21
128 | 64 | 4.3 | 2.0 | 0.23477 | 0.26387 | 0.00025 | 8 | 21
128 | 96 | 4.0 | 2.0 | 0.23468 | 0.26407 | 0.00069 | 6 | 32
128 | 96 | 4.3 | 2.0 | 0.23484 | 0.26377 | 0.00024 | 8 | 32
192 | 96 | 4.0 | 2.0 | 0.23459 | 0.26356 | 0.00030 | 9 | 32
192 | 128 | 4.0 | 2.0 | 0.23462 | 0.26347 | 0.00029 | 9 | 43
Table 1: Grid sensitivity study conducted for $Ha=1200$ and $Gr=10^{8}$:
$\tau_{Ha}$ and $\tau_{Sh}$ are the wall friction forces, $\epsilon$ is the
absolute error of the balance (18). The number of grid points inside the
Hartmann and Shercliff boundary layers are $N_{Ha}$ and $N_{Sh}$,
respectively.

The parameters of the grids in the entire studied parameter range of $Ha$ and
$Gr$ have been determined in the same way. It has been found that the value of
$Gr$ does not affect the selection at fixed $Ha$. This effect can be explained
by the presence of very strong magnetic fields which fully suppress transverse
circulation (see section 3.1 for a discussion). The summary of the grids used
in the simulations is presented in table 2.

A grid sensitivity study has also been conducted to determine the minimum
number of grid points $N_{x}$ required in the linear stability analysis. It
has been found that the growth of linear unstable modes is accurately
reproduced at $N_{x}=32$ for modes with $\lambda\lesssim 2$, which needs to be
increased to $N_{x}=64$ for greater $\lambda$.

Computational domains of length $L_{x}=4\pi$ or $L_{x}=2\pi$ are used in DNS.
As we will see below, these lengths are substantially larger than the
streamwise wavelength of the fastest growing instability modes. The flow
structures have been accurately resolved with, respectively, $384$ or $192$
grid points in the $x-$direction.

The time steps adjusted to secure numerical instability and, thus, varying
with $Ha$ and $Gr$, but never exceeding $2.5\times 10^{-3}$, are used in
linear stability and DNS simulations.

$Ha$ | $N_{y}$ | $N_{z}$ | $A_{y}$ | $A_{z}$ | $N_{Ha}$ | $N_{Sh}$
---|---|---|---|---|---|---
1000 | 192 | 96 | 4.0 | 2.0 | 9 | 32
2000 | 256 | 96 | 4.0 | 2.0 | 8 | 28
3000 | 384 | 96 | 4.0 | 2.0 | 9 | 24
4000 | 384 | 96 | 4.0 | 2.0 | 7 | 22
5000 | 512 | 96 | 4.0 | 2.0 | 8 | 20
6000 | 512 | 96 | 4.0 | 2.0 | 7 | 19
7000 | 512 | 96 | 4.3 | 2.0 | 9 | 18
8000 | 512 | 96 | 4.3 | 2.0 | 8 | 17
9000 | 512 | 96 | 4.3 | 2.0 | 7 | 16
10000 | 512 | 96 | 4.3 | 2.0 | 6 | 15
Table 2: Parameters of the computational grids used in the simulations for
$Gr=10^{8}-10^{10}$. The number of grid points inside the Hartmann and
Shercliff boundary layers are $N_{Ha}$ and $N_{Sh}$, respectively.

##  3 Results

###  3.1 Base flow

The flow structure for several typical cases is illustrated in figures 2 and
3. The results for all the completed simulations are summarized in table 3.
The table shows the type of the flow for a particular regime (Q2D or 3D flows
to be discussed shortly), the maximum and minimum values of $U_{x}$, integral
quantities, such the wall friction force $d\hat{p}/dx$, the volume-averaged
kinetic energies of streamwise and transverse velocities

$E_{x}=A^{-1}\int_{-1}^{1}\int_{-1}^{1}U_{x}^{2}dydz,\ \
E_{t}=A^{-1}\int_{-1}^{1}\int_{-1}^{1}\left(U_{y}^{2}+U_{z}^{2}\right)dydz,$
(20)

and the mean square of temperature perturbations

$E_{\theta}=A^{-1}\int_{-1}^{1}\int_{-1}^{1}\Theta^{2}dydz.$ (21)

Similar data for flows with lower values of $Gr$ and $Ha$ in a duct with
$\Gamma=1$ can be found in (Zhang & Zikanov, 2014).

The flow structure is predominantly determined by the effect of
magnetoconvection. Similarly to the findings of Zhang & Zikanov (2014), we
observe two regimes of the base flow depending on whether $Gr$ is smaller or
larger than a certain threshold $Gr^{\ast}(Ha)$. The Q2D regime observed at
$Gr<Gr^{\ast}$ is characterized by the transverse convection-induced
circulation entirely suppressed by the strong magnetic field. The
distributions of the temperature and streamwise velocity are nearly one-
dimensional outside of the Hartmann boundary layers ($\Theta\approx\Theta(z)$
and $U_{x}\approx U_{x}(z)$). In the absence of transverse circulation, the
distributions of temperature are determined by the balance between the heat
conduction and the heat convection by $U_{x}$. Examples of this regime are
shown in figure 2$a$ and figure 3.

The 3D regime observed at $Gr>Gr^{\ast}$($Ha$), when the strength of the
magnetic field is insufficient to suppress convection circulation, is
characterized by significant transverse flow and fully 2D variations of
temperature. (see figure 2$b-e$).

To avoid confusion, it is pertinent to repeat the terminology here. The terms
3D and Q2D are used in this paper to describe the general flow transformation
caused by the magnetic field, i.e. suppression of velocity and temperature
gradients along the magnetic field lines in the core of the duct and formation
of thin Hartmann boundary layers. The base flow, in which streamwise
uniformity is also imposed, becomes, respectively, 2D and Q1D.

We find the Q2D regime in the larger part of the explored range of $Ha$ and
$Gr$ including the most interesting cases of large $Ha$ (see the rightmost
column in table 3). The total friction force increases at stronger magnetic
fields. The visible effect of convection is the asymmetry of the velocity
profile, with $U_{x}$ larger in the bottom than in the top half (see, e.g.,
figure 3$b,d$). The cause of the asymmetry has been explained in section 2.1.

(a)y(b)zyy(c)zyy(d)zyy(e)zyy Figure 2: Base flow at $Ha=1000$, $Gr=10^{8}$
($a$), $Ha=1000$, $Gr=10^{9}$ ($b$), $Ha=1000$, $Gr=10^{10}$ ($c$), $Ha=2000$,
$Gr=10^{10}$ ($d$) and $Ha=3000$, $Gr=10^{10}$ ($e$). Vector fields and
streamlines of transverse circulation ($u_{y}$, $u_{z}$) are shown in the left
column (not in ($a$), since the velocity’s amplitude is virtually zero in this
case). The middle and right columns show distributions of temperature
fluctuations and streamwise velocity $u_{x}$, respectively. Solid and dashed
isolines in the middle column indicate positive and negative values,
respectively. The wall heating is at $z=-0.2857$, and the magnetic field is in
the $y$-direction.

(a)(b)(c)(d)

Figure 3: Q2D base flow at $Ha=5000$, $Gr=10^{9}$ ($a$), $Ha=5000$, $Gr=10^{10}$ ($b$), $Ha=10000$, $Gr=10^{9}$ ($c$) and $Ha=10000$, $Gr=10^{10}$ ($d$). Distribution of streamwise velocity $u_{x}$ is shown. The wall heating is at $z=-0.2857$, and the magnetic field is in the $y$-direction. $Ha$ | $Gr$ | $Gr/{\it Ha^{2}}$ | $-d\hat{p}/dx$ | $E_{x}$ | $E_{t}$ | $E_{\theta}$ | $u_{x,min}$ | $u_{x,max}$ | $\it{Regime}$
---|---|---|---|---|---|---|---|---|---
$1000$ | $10^{8}$ | $1.00\times 10^{2}$ | $0.222$ | $1.050$ | $1.62\times 10^{-8}$ | $3.59\times 10^{-3}$ | $0.0$ | $1.139$ | $Q2D$
$1000$ | $10^{9}$ | $1.00\times 10^{3}$ | $0.224$ | $1.056$ | $6.32\times 10^{-3}$ | $9.31\times 10^{-3}$ | $0.0$ | $1.249$ | $3D$
$1000$ | $10^{10}$ | $1.00\times 10^{4}$ | $0.254$ | $1.892$ | $3.11\times 10^{-2}$ | $7.11\times 10^{-4}$ | $-0.402$ | $2.601$ | $3D$
$2000$ | $10^{8}$ | $2.50\times 10^{1}$ | $0.430$ | $1.035$ | $5.42\times 10^{-9}$ | $3.61\times 10^{-3}$ | $0.0$ | $1.079$ | $Q2D$
$2000$ | $10^{9}$ | $2.50\times 10^{2}$ | $0.431$ | $1.037$ | $4.96\times 10^{-7}$ | $3.34\times 10^{-3}$ | $0.0$ | $1.153$ | $Q2D$
$2000$ | $10^{10}$ | $2.50\times 10^{3}$ | $0.447$ | $1.313$ | $8.30\times 10^{-3}$ | $2.19\times 10^{-3}$ | $0.0$ | $1.968$ | $3D$
$3000$ | $10^{8}$ | $1.11\times 10^{1}$ | $0.636$ | $1.029$ | $1.76\times 10^{-9}$ | $3.62\times 10^{-3}$ | $0.0$ | $1.063$ | $Q2D$
$3000$ | $10^{9}$ | $1.11\times 10^{2}$ | $0.637$ | $1.030$ | $1.66\times 10^{-7}$ | $3.45\times 10^{-3}$ | $0.0$ | $1.114$ | $Q2D$
$3000$ | $10^{10}$ | $1.11\times 10^{3}$ | $0.643$ | $1.172$ | $1.58\times 10^{-2}$ | $4.48\times 10^{-3}$ | $0.0$ | $1.683$ | $3D$
$4000$ | $10^{8}$ | $6.25\times 10^{0}$ | $0.842$ | $1.025$ | $3.85\times 10^{-10}$ | $3.63\times 10^{-3}$ | $0.0$ | $1.054$ | $Q2D$
$4000$ | $10^{9}$ | $6.25\times 10^{1}$ | $0.842$ | $1.026$ | $3.70\times 10^{-8}$ | $3.53\times 10^{-3}$ | $0.0$ | $1.093$ | $Q2D$
$4000$ | $10^{10}$ | $6.25\times 10^{2}$ | $0.847$ | $1.109$ | $2.75\times 10^{-6}$ | $2.96\times 10^{-3}$ | $0.0$ | $1.522$ | $Q2D$
$5000$ | $10^{8}$ | $4.00\times 10^{0}$ | $1.046$ | $1.022$ | $1.13\times 10^{-10}$ | $3.63\times 10^{-3}$ | $0.0$ | $1.048$ | $Q2D$
$5000$ | $10^{9}$ | $4.00\times 10^{1}$ | $1.046$ | $1.023$ | $1.10\times 10^{-8}$ | $3.57\times 10^{-3}$ | $0.0$ | $1.080$ | $Q2D$
$5000$ | $10^{10}$ | $4.00\times 10^{2}$ | $1.049$ | $1.078$ | $8.55\times 10^{-7}$ | $3.16\times 10^{-3}$ | $0.0$ | $1.425$ | $Q2D$
$6000$ | $10^{8}$ | $2.78\times 10^{0}$ | $1.251$ | $1.020$ | $1.00\times 10^{-10}$ | $3.63\times 10^{-3}$ | $0.0$ | $1.043$ | $Q2D$
$6000$ | $10^{9}$ | $2.78\times 10^{1}$ | $1.251$ | $1.021$ | $9.79\times 10^{-9}$ | $3.57\times 10^{-3}$ | $0.0$ | $1.071$ | $Q2D$
$6000$ | $10^{10}$ | $2.78\times 10^{2}$ | $1.254$ | $1.060$ | $7.75\times 10^{-7}$ | $3.18\times 10^{-3}$ | $0.0$ | $1.362$ | $Q2D$
$7000$ | $10^{8}$ | $2.04\times 10^{0}$ | $1.454$ | $1.019$ | $9.00\times 10^{-11}$ | $3.63\times 10^{-3}$ | $0.0$ | $1.040$ | $Q2D$
$7000$ | $10^{9}$ | $2.04\times 10^{1}$ | $1.454$ | $1.019$ | $8.80\times 10^{-9}$ | $3.58\times 10^{-3}$ | $0.0$ | $1.064$ | $Q2D$
$7000$ | $10^{10}$ | $2.04\times 10^{2}$ | $1.457$ | $1.048$ | $7.12\times 10^{-7}$ | $3.20\times 10^{-3}$ | $0.0$ | $1.315$ | $Q2D$
$8000$ | $10^{8}$ | $1.56\times 10^{0}$ | $1.658$ | $1.018$ | $8.12\times 10^{-11}$ | $3.63\times 10^{-3}$ | $0.0$ | $1.037$ | $Q2D$
$8000$ | $10^{9}$ | $1.56\times 10^{1}$ | $1.658$ | $1.018$ | $7.96\times 10^{-9}$ | $3.58\times 10^{-3}$ | $0.0$ | $1.058$ | $Q2D$
$8000$ | $10^{10}$ | $1.56\times 10^{2}$ | $1.661$ | $1.040$ | $6.56\times 10^{-7}$ | $3.22\times 10^{-3}$ | $0.0$ | $1.280$ | $Q2D$
$9000$ | $10^{8}$ | $1.23\times 10^{0}$ | $1.862$ | $1.017$ | $2.50\times 10^{-11}$ | $3.63\times 10^{-3}$ | $0.0$ | $1.035$ | $Q2D$
$9000$ | $10^{9}$ | $1.23\times 10^{1}$ | $1.862$ | $1.017$ | $2.46\times 10^{-9}$ | $3.60\times 10^{-3}$ | $0.0$ | $1.054$ | $Q2D$
$9000$ | $10^{10}$ | $1.23\times 10^{2}$ | $1.864$ | $1.035$ | $2.11\times 10^{-7}$ | $3.34\times 10^{-3}$ | $0.0$ | $1.252$ | $Q2D$
$10000$ | $10^{8}$ | $1.00\times 10^{0}$ | $2.067$ | $1.016$ | $2.31\times 10^{-11}$ | $3.63\times 10^{-3}$ | $0.0$ | $1.033$ | $Q2D$
$10000$ | $10^{9}$ | $1.00\times 10^{1}$ | $2.067$ | $1.016$ | $2.28\times 10^{-9}$ | $3.60\times 10^{-3}$ | $0.0$ | $1.050$ | $Q2D$
$10000$ | $10^{10}$ | $1.00\times 10^{2}$ | $2.068$ | $1.031$ | $1.97\times 10^{-7}$ | $3.36\times 10^{-3}$ | $0.0$ | $1.229$ | $Q2D$
Table 3: Integral characteristics and type of the computed base flow states.

The 3D regimes are only found in a limited range of moderate values of $Ha$ at
$Gr=10^{9}$ and $10^{10}$ (see table 3). The transverse circulation consist of
a single roll (see figure 2$b,e$) or two symmetric rolls (see figure 2$c,d$).
The single roll has no preferred circulation direction and may appear in the
solution either as shown in the figure 2$b$ or as a symmetric reflection with
respect to the vertical midplane (see figure 2$e$ for an example). The
circulation causes a visibly 2D distributions of $\Theta$ and $U_{x}$, and, at
the same values of $Gr$, a decrease of $E_{\theta}$ as a result of mixing (see
table 3). We observe stronger top–bottom asymmetries or even formation of
reverse flow in the top portion of the duct as the strength of convection
increases at fixed $Ha$ (see figure 2$a-c$ for an example). As an illustration
of the asymmetry the minimum and maximum values of $u_{x}$ are shown in table
3.

The classification of the flow regimes into Q2D and 3D can also be described
in terms of the values of the average kinetic energy of transverse circulation
$E_{t}$. Our data shown in table 3 and figure 4 are in a good qualitative
agreement with the results of Zhang & Zikanov (2014). The observed differences
can be attributed to the substantially different studied ranges of $Ha$ and
$Gr$ and different aspect ratio.

We see in figure 4 that, at a fixed Reynolds number considered in this work,
the intensity of the transverse circulation is well approximated by a function
of the single control parameter - the combination $Gr/Ha^{2}$. Analyzing the
flow structures at various values of $E_{t}$ we find a clear demarcation
between 3D and Q2D regimes. $E_{t}$ is greater than, approximately, $10^{-4}$
in 3D regimes. Q2D flows all have values of $E_{t}$ less than $10^{-6}$.

Our interest in this study is primarily in the Q2D regimes. The 3D regimes are
not considered in the rest of the paper.

$E_{t}$$Gr/Ha^{2}$Zhang & Zikanov (2014)3DQ2D Figure 4: The average kinetic
energy of transverse circulation in the base flow $E_{t}$ as a function of
$Gr/Ha^{2}$. Circles indicate the numerical results of Zhang & Zikanov (2014)
at $\Gamma=1.0$. Stars and crosses indicate, respectively, Q2D and 3D regimes
found in this work for the flow at $\Gamma=3.5$. Values of $E_{t}$ and
$Gr/Ha^{2}$ for each computed flow can be found in table 3.

###  3.2 Applicability of the SM82 model

In this section, we investigate the applicability of the SM82 model to
analysis of magnetoconvection instability.

####  3.2.1 Base flow

For a streamwise-uniform, unidirectional, steady-state base flow , the SM82
model equations (14)$-$(16) are reduced to a system of linear ordinary
differential equations. The solution satisfying the boundary conditions is:

${U_{x}}(z)=C\frac{Re}{Ha}c_{1}+A\frac{Re}{Ha\Gamma}s_{1}-A\frac{Re}{Ha}z-C\frac{Re}{Ha},$
(22)

$\Theta(z)=\frac{C}{2}\frac{Re\Gamma}{Ha^{2}}c_{1}+\frac{A}{2}\frac{Re}{Ha^{2}}s_{1}-\frac{A}{12}\frac{Re\Gamma}{Ha}z^{3}-\frac{C}{4}\frac{Re\Gamma}{Ha}z^{2}+a_{1}z,$
(23)

where

$a_{1}=-\frac{C}{2}\frac{Re\Gamma}{Ha^{3/2}}t_{1}-\frac{A}{2}\frac{Re\Gamma}{Ha^{3/2}}t_{1}^{-1}+\frac{A}{4}\frac{Re}{Ha\Gamma}+\frac{C}{2}\frac{Re}{Ha},$
(24)

$c_{1}=\frac{cosh(\sqrt{Ha}z)}{cosh(\sqrt{Ha}/\Gamma)},\ \ \
s_{1}=\frac{sinh(\sqrt{Ha}z)}{sinh(\sqrt{Ha}/\Gamma)},\ \ \
t_{1}=tanh(\sqrt{Ha}/\Gamma),$ (25)

$A=\frac{\Gamma}{2RePr}\frac{Gr}{Re^{2}},\ \ \
C=\frac{Ha}{Re\Gamma}\frac{1}{t_{1}/\sqrt{Ha}-1/\Gamma}.$ (26)

The profiles (22) and (23) are shown for $Gr=10^{8},10^{9},10^{10}$ and
several values of $Ha$ in figure 5. For comparison, respective profiles of the
true $U_{x}$ and $\Theta$ obtained for the computed base flow solutions
described in section 3.1 are shown for the midplane $y=0$. The symbols 2D and
3D correspond to, respectively, approximate SM82 and computed solutions.

(a)$Gr=10^{8}$z${\it U_{x}}$${\it Ha}=1000,3D$${\it Ha}=5000,3D$${\it
Ha}=10000,3D$${\it Ha}=1000,2D$${\it Ha}=5000,2D$${\it
Ha}=10000,2D$(d)z${\it\Theta}$(b)$Gr=10^{9}$${\it
U_{x}}$(e)${\it\Theta}$(c)$Gr=10^{10}$${\it U_{x}}$(f)${\it\Theta}$ Figure 5:
Base flow profiles at $Gr=10^{8}$ ($a,d$), $Gr=10^{9}$ ($b,e$), and
$Gr=10^{10}$ ($c,f$). The SM82 model solutions (22), (23) and the
distributions along the midplane $y=0$ obtained in the full numerical
solutions of section 3.1 are denoted, respectively, as 2D and 3D. The curves
for $Ha=10^{3}$ are only shown at $Gr=10^{8}$ ($a,d$), since at higher $Gr$
the base flow is not Q2D (see section 3.1). The top and bottom rows show the
profiles of streamwise velocity and temperature, respectively. The insets in
($d$) and ($e$) show zoomed-in illustrations of the temperature profiles near
the top wall.

Good agreement between the computed solutions and the solutions of the SM82
model is evident at $Gr=10^{8}$ and $10^{9}$. There are some deviations
between the computed and model profiles of $U_{x}$, but they are small and
decrease with increasing $Ha$. The agreement is significantly worse in the
case of the flows at $Gr=10^{10}$. Here we observe large deviations between
the computed and model curves at $Ha=5000$ and smaller, yet still significant
deviations at $Ha=10000$.

####  3.2.2 Linear stability analysis

In order to verify applicability of the SM82 model to linear stability
analysis of Q2D flows, instability of several high$-Ha$ flows was evaluated
twice: once using full 3D model of the base flow and perturbations and once
entirely in the framework of the 2D SM82 model. The results are presented in
figure 6 and table 4. We see good agreement between predictions of 2D and 3D
models at $Gr=10^{8}$ and $10^{9}$. The accuracy improves with growing $Ha$.
As an example, the average relative difference between the values of the
growth rate $\gamma$ for the two models is $33\%$ at $Ha=2000$, $Gr=10^{9}$
and $11\%$ at $Ha=3000$, $Gr=10^{9}$. The situation is less clear for flows at
$Gr=10^{10}$ (see figure 6$c$ and the last six lines of table 4). Here we only
see a qualitative agreement. The shape of the $\gamma$($\lambda$) curves, the
wavelength of the most unstable mode, and the effect of $Ha$ on stability are
similar in the 3D and 2D solution. The quantitative agreement is, however,
poor, with the difference between the values of $\gamma$ found for the two
models being about $50\%$.

The quantitative disagreement between the base flow profiles and, as an
evident consequence, stability properties found in the 3D and 2D models is
difficult to interpret. The velocity and temperature distributions computed in
the framework of the 3D model clearly show that the base flow is Q2D and
nearly perfectly unidirectional at $Ha$ higher than approximately $4000$ (see
figure 3$b,d$ and values of $E_{t}$ in table 3). As illustrated in section
3.3, fields of growing perturbations also remain Q2D at such high $Ha$.
Additional calculations performed with larger grids and longer times of flow
evolution did not lead to significant changes. The deviations from quasi-two-
dimensionality and inaccuracy of the numerical model are, therefore, excluded
as possible reasons.

$({\it a})$$Gr=10^{8}$${\it\gamma}$${\it\lambda}$${\it Ha}=1000,3D$${\it Ha}=1000,2D$${\it Ha}=2000,3D$${\it Ha}=2000,2D$$({\it b})$$Gr=10^{9}$${\it\gamma}$${\it\lambda}$${\it Ha}=2000,3D$${\it Ha}=2000,2D$${\it Ha}=3000,3D$${\it Ha}=3000,2D$$({\it c})$$Gr=10^{10}$${\it\gamma}$${\it\lambda}$${\it Ha}=4000,3D$${\it Ha}=4000,2D$${\it Ha}=5000,3D$${\it Ha}=5000,2D$${\it Ha}=10000,3D$${\it Ha}=10000,2D$ Figure 6: Rates of exponential growth $\gamma$ shown as functions of the axial wavelength $\lambda$ at $Ha=1000,2000$, $Gr=10^{8}$ ($a$), $Ha=2000,3000$, $Gr=10^{9}$ ($b$), and $Ha=4000,5000$, $Gr=10^{10}$ ($c$). The results of 3D and 2D (SM82) models are denoted as filled and empty circles, respectively. $Ha$ | $Gr$ |  |  |  |  |  | $\lambda$ |  |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---|---
|  | $0.2$ | $0.3$ | $0.4$ | $0.5$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1.0$ | $1.5$ | $2.0$
1000 | $10^{8}$ | $-$ | 0.164 | 0.360 | 0.474 | 0.532 | 0.553 | 0.552 | 0.538 | 0.518 | 0.396 | 0.301
$\underline{1000}$ | $\underline{10^{8}}$ | $-$ | $\underline{0.049}$ | $\underline{0.269}$ | $\underline{0.398}$ | $\underline{0.461}$ | $\underline{0.486}$ | $\underline{0.487}$ | $\underline{0.475}$ | $\underline{0.455}$ | $\underline{0.338}$ | $\underline{0.247}$
2000 | $10^{8}$ | $-$ | $-$ | 0.251 | 0.387 | 0.461 | 0.495 | 0.504 | 0.499 | 0.485 | 0.381 | 0.295
$\underline{2000}$ | $\underline{10^{8}}$ | $-$ | $-$ | $\underline{0.200}$ | $\underline{0.345}$ | $\underline{0.486}$ | $\underline{0.422}$ | $\underline{0.456}$ | $\underline{0.459}$ | $\underline{0.445}$ | $\underline{0.344}$ | $\underline{0.261}$
2000 | $10^{9}$ | 1.563 | 2.602 | 3.010 | 3.136 | 3.134 | 3.069 | 2.972 | 2.859 | 2.741 | 2.175 | $-$
$\underline{2000}$ | $\underline{10^{9}}$ | $\underline{1.028}$ | $\underline{2.029}$ | $\underline{2.459}$ | $\underline{2.626}$ | $\underline{2.671}$ | $\underline{2.651}$ | $\underline{2.593}$ | $\underline{2.513}$ | $\underline{2.421}$ | $\underline{1.941}$ | $-$
3000 | $10^{9}$ | 1.540 | 2.587 | 2.993 | 3.117 | 3.114 | 3.048 | 2.950 | 2.836 | 2.715 | 2.141 | $-$
$\underline{3000}$ | $\underline{10^{9}}$ | $\underline{1.085}$ | $\underline{2.162}$ | $\underline{2.619}$ | $\underline{2.791}$ | $\underline{2.827}$ | $\underline{2.792}$ | $\underline{2.717}$ | $\underline{2.622}$ | $\underline{2.516}$ | $\underline{1.990}$ | $-$
4000 | $10^{10}$ | 10.234 | 11.553 | 11.674 | 11.435 | 11.067 | 10.645 | 10.202 | 9.754 | 9.312 | $-$ | $-$
$\underline{4000}$ | $\underline{10^{10}}$ | $\underline{6.881}$ | $\underline{7.404}$ | $\underline{7.415}$ | $\underline{7.156}$ | $\underline{6.833}$ | $\underline{6.494}$ | $\underline{6.158}$ | $\underline{5.834}$ | $\underline{5.525}$ | $-$ | $-$
5000 | $10^{10}$ | 10.605 | 11.918 | 12.023 | 11.764 | 11.371 | 10.925 | 10.458 | 9.988 | 9.527 | $-$ | $-$
$\underline{5000}$ | $\underline{10^{10}}$ | $\underline{7.454}$ | $\underline{8.269}$ | $\underline{8.273}$ | $\underline{8.018}$ | $\underline{7.673}$ | $\underline{7.303}$ | $\underline{6.934}$ | $\underline{6.577}$ | $\underline{6.237}$ | $-$ | $-$
10000 | $10^{10}$ | 10.807 | 12.064 | 12.106 | 11.791 | 11.351 | 10.865 | 10.366 | 9.871 | 9.389 | $-$ | $-$
$\underline{10000}$ | $\underline{10^{10}}$ | $\underline{8.729}$ | $\underline{9.890}$ | $\underline{9.942}$ | $\underline{9.668}$ | $\underline{9.291}$ | $\underline{8.882}$ | $\underline{8.468}$ | $\underline{8.059}$ | $\underline{7.661}$ | $-$ | $-$
Table 4: Results of the linear stability analysis of 3D and 2D flow solutions.
Rates of exponential growth $\gamma$ are shown as functions of the axial
wavelength $\lambda$. The results of the SM82 model are underlined. The growth
rates are determined as in (17).

A further useful, albeit not fully explaining illustration is provided in
figure 7. Computed base flow distributions of $U_{x}$ and the streamwise
component of the Lorentz force $F_{Lx}$ are shown for $Ha=10000$ and
$Gr=10^{8},10^{9}$ and $10^{10}$ within and near the Hartmann boundary layer.
We see that the strong vertical variation of $U_{x}$ existing at $Gr=10^{10}$
extends toward the Hartmann wall and causes a respective variation of the
Lorentz force. It must be noted that this picture does not contradict to the
identification of the flow as Q2D. The profiles $U_{x}(y,z=const)$ collapse
into one curve corresponding to the Hartmann profile if taken outside the
sidewall layers at the horizontal walls and scaled by the respective maximum
values of $U_{x}$.

(a)z(b)(c)z(d)(e)zy(f)y Figure 7: Base flow at $Gr=10^{8}$ ($a,b$),
$Gr=10^{9}$ ($c,d$), and $Gr=10^{10}$ ($e,f$) for $Ha=10000$. The left and
right columns show, respectively, distributions of streamwise velocity $U_{x}$
and the streamwise component of the Lorentz force $F_{Lx}$ near the Hartmann
wall. The red dashed line shows the boundary of the Hartmann layer of
thickness $\delta_{Ha}={\it Ha}^{-1}$.

####  3.2.3 Nonlinear flows

The results presented so far in this section indicate that the SM82 2D model
may also inaccurately describe the nonlinear flow regimes developing as a
result of the instability at $Gr=10^{10}$. As a test of this possibility,
comparison between the results of 2D and 3D models at $Gr=10^{10}$, $Ha=10000$
is illustrated in figures 8 and 9 and discussed below. The procedure of
computing nonlinear flows is described in section 3.4. Here we only mention
that the same numerical resolution is used in 2D and 3D models. The shorter
wavelength domain length $L_{x}=2\pi$ is used in the 3D model. This rather
small length has no significant effect on the flow evolution as it has been
confirmed in the additional 2D simulations conducted with $L_{x}=4\pi$ (shown
in figure 8) and $L_{x}=2\pi$.

(a)${\it z}$${\it u_{z}}$(b)${\it z}$${\it u_{x}}$(c)${\it
z}$x${\theta}$(d)${\it z}$${\it u_{z}}$(e)${\it z}$${\it u_{x}}$(f)${\it
z}$x${\theta}$(g)$z$$u_{x}$(h)$z$$\theta$$3D$$2D$$3D,base$$2D,base$ Figure 8:
Flow structure in nonlinear regime at $Gr=10^{10}$, $Ha=10000$. The
instantaneous distributions of $u_{z}$ ($a$, $d$), $u_{x}$ ($b$, $e$), and
$\theta$ ($c$, $f$) obtained in 2D and 3D (plotted in the midplane $y=0$)
models are shown in ($a$)-($c$) and ($d$)-($f$), respectively. The profiles of
$u_{x}$ and $\theta$ obtained by averaging over $x$ and time (with the base
flow profiles from figure 5$c,f$) are shown, respectively, in ($g$) and ($h$).
(a)$\theta$$t$(c)$\theta$$t$${\it z}=-0.2857$${\it z}=0$${\it
z}=0.2857$(b)$u_{x}$$t$(d)$u_{x}$$t$${\it z}=-0.2789$${\it z}=0$${\it
z}=0.2789$ Figure 9: Time signals of temperature ($a,c$) and streamwise
velocity ($b,d$) at $y=0$ (in the 3D flow) in the nonlinear regime at
$Ha=10^{4}$, $Gr=10^{10}$ are shown for 2D ($a,b$) and 3D ($c,d$) models.

The simulations show that the 3D flow remains Q2D at these values of $Ha$ and
$Gr$. At the same time, its reproduction by the 2D SM82 model is inaccurate in
some aspects. The instantaneous distributions of velocity components and
temperature shown in figure 8 clearly illustrate the difference. 2D
approximations (see figure 8$a$-$c$) are similar to 3D flows (see figure
8$d$-$f$) in terms of the largest typical streamwise wavelength (about $1.5$)
but demonstrate noticeably less regular pattern and higher energy in shorter
wavelengths. The difference is reflected by the point signals of velocity and
temperature shown in figure 9. Comparing figures 9$b$ and 9$d$, we also see
that the 2D approximation substantially underestimates the typical amplitude
of velocity fluctuations. As an example, the standard deviation for the signal
of $u_{x}$ in the middle of the duct ($z=0$) is $0.33$ for the 2D model and
$0.71$ for the 3D model. Computed values of volume-average kinetic energy of
the fluctuations (not shown) confirm this conclusion. The values of energy
found in the 2D approximation are about an order of magnitude lower than in
the actual flow computed in the framework of the 3D model. Interestingly, the
misrepresentation of the structure of velocity and temperature fluctuations by
the 2D model does not lead to a similarly strong inaccuracy in the prediction
of the effect of mixing by fluctuations. Profiles of average streamwise
velocity and temperature in figures 8($g,h$) show strong change in comparison
with the base flow, but only moderate differences between the 2D and 3D
results.

####  3.2.4 Applicability of the SM82 model: summary

We conclude that the SM82 approximation accurately represents Q2D flows at
moderately large $Gr$ ($10^{8}$ and $10^{9}$ in our system). The accuracy
deteriorates at higher $Gr$ even though the flow remains Q2D. An example of
this is observed at $Ha\gtrsim 4000$ and $Gr=10^{10}$. The base flow profiles
are clearly different between the 2D and 3D models (see figures 5$c,f$). The
2D linear stability analysis is qualitatively correct in the sense that it
correctly predicts the principal type of the unstable perturbations and the
streamwise wavelength of the most unstable mode (see figure 6$c$). The values
of the growth rate $\gamma$ are, however, substantially underestimated by the
2D model (see figure 6$c$ and table 4). The nonlinear flow states resulting
from the instability are predicted incorrectly by the 2D model, which adds
artificial irregularity and short-wave fluctuations and underestimates the
amplitude of velocity fluctuations (see figures 8 and 9).

We do not have a satisfactory explanation to this effect and leave its further
exploration for future studies. It should be mentioned that the quality of the
2D approximation improves with increasing $Ha$. As an example, the values of
the linear instability growth rate $\gamma$ shown in table 4 are
underpredicted by the 2D model by about $35\%$ at $Ha=4000$ and by about only
$18\%$ at $Ha=10000$.

In the remaining part of this paper, the discussion of the instability and
nonlinear states is based on the SM82 2D approximation for flows at
$Gr=10^{8}$ and $10^{9}$ and on the full 3D solutions for flows at $10^{10}$.

###  3.3 Results of linear stability analysis

The results of the linear stability analysis are summarized in figures 10$-$12
and table 5. We need to mention that the wavelength $\lambda$ is varied with
step $0.1$ in the simulations. The computed values of the exponential growth
rate $\gamma$ as a function of the wavelength $\lambda$ for various
combinations of $Ha$ and $Gr$ are shown in figure 10. Two trends of the linear
stability behaviour were proposed by Zhang & Zikanov (2014): (1) a higher
growth rate and shorter wavelength appear at higher $Gr$; (2) an increase of
$Ha$ leads to a higher growth rate. The second, apparently counterintuitive
effect was attributed by Zhang & Zikanov (2014) to modification of the base
flow, namely to suppression of the transverse circulation resulting in weaker
mixing and stronger unstable temperature stratification.

In order to investigate these trends for our system and parameter range, we
present the exponential growth rate $\gamma_{max}$ of the fastest growing
modes and the corresponding wavelengths $\lambda_{max}$ in figure 11 and table
5. Our results are clearly consistent with the first trend, $\gamma_{max}$
increases and $\lambda_{max}$ decreases with growing $Gr$. However, we find a
different behaviour in regard of the second trend. The increase of $Ha$ leads
to noticeable or slight decrease of the exponential growth rates at
$Gr=10^{8}$ or $Gr=10^{9}$, respectively (see figure 10$c$). For $Gr=10^{10}$,
$\gamma$ is nearly insensitive to the values of $Ha$.

We conclude that the counterintuitive behaviour of stronger instability at
higher $Ha$ is not observed in Q2D flows at high $Ha$ considered in our study.
It is, nevertheless, noteworthy that $\gamma_{max}$ does not decrease with
$Ha$ at $Gr=10^{10}$. Increase slightly at $Ha\lesssim 6000$ and remains
practically constant at higher $Ha$ The enhanced friction in the Hartmann
boundary layers is compensated by another effect, the only plausible candidate
for which is the strong modification of the streamwise velocity profile
visible in figure 5$c$. $U_{x}$ strongly grows with $Ha$ in the bottom half of
the duct, i.e., its part with the strongest unstable temperature gradient.

The findings have critical implications for design and operation of the fusion
reactor systems, since they indicate that strength of the convection
instability is not diminished by strong magnetic fields at $Gr\geq 10^{10}$
typical for reactor blankets.

(a)${\it\gamma}$${\it\lambda}$${\it Ha}=1000$${\it Ha}=2000$${\it Ha}=3000$${\it Ha}=4000$${\it Ha}=5000$${\it Ha}=6000$${\it Ha}=7000$${\it Ha}=8000$${\it Ha}=9000$${\it Ha}=10000$(b)${\it\gamma}$${\it\lambda}$(c)${\it\gamma}$${\it\lambda}$ Figure 10: Rates of exponential growth $\gamma$ shown as functions of axial wavelength $\lambda$ at $Gr=10^{8}$ ($a$), $Gr=10^{9}$ ($b$), and $Gr=10^{10}$ ($c$) for various values of $Ha$. Results of 2D approximation and 3D computational analysis are shown, respectively, for $Gr=10^{8},10^{9}$ and $Gr=10^{10}$. Additional plots in ($b$) and ($c$) show zoomed-in areas around ($\lambda_{max},\gamma_{max}$).  $Gr$ | $Ha$ | $\lambda_{max}$ | $\gamma_{max}$ |  | $Gr$ | $Ha$ | $\lambda_{max}$ | $\gamma_{max}$ |  | $Gr$ | $Ha$ | $\lambda_{max}$ | $\gamma_{max}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$10^{8}$ | $1000$ | $0.8$ | 0.487 |  | $10^{9}$ | $1000$ | $-$ | $-$ |  | $10^{10}$ | $1000$ | $-$ | $-$
$10^{8}$ | $2000$ | $0.8$ | 0.465 |  | $10^{9}$ | $2000$ | $0.6$ | 2.671 |  | $10^{10}$ | $2000$ | $-$ | $-$
$10^{8}$ | $3000$ | $0.9$ | 0.400 |  | $10^{9}$ | $3000$ | $0.6$ | 2.827 |  | $10^{10}$ | $3000$ | $-$ | $-$
$10^{8}$ | $4000$ | $0.9$ | 0.336 |  | $10^{9}$ | $4000$ | $0.6$ | 2.856 |  | $10^{10}$ | $4000$ | $0.4$ | 11.674
$10^{8}$ | $5000$ | $1.0$ | 0.277 |  | $10^{9}$ | $5000$ | $0.6$ | 2.833 |  | $10^{10}$ | $5000$ | $0.4$ | 12.023
$10^{8}$ | $6000$ | $1.1$ | 0.226 |  | $10^{9}$ | $6000$ | $0.6$ | 2.786 |  | $10^{10}$ | $6000$ | $0.4$ | 12.094
$10^{8}$ | $7000$ | $1.1$ | 0.184 |  | $10^{9}$ | $7000$ | $0.6$ | 2.726 |  | $10^{10}$ | $7000$ | $0.4$ | 12.112
$10^{8}$ | $8000$ | $1.2$ | 0.149 |  | $10^{9}$ | $8000$ | $0.6$ | 2.659 |  | $10^{10}$ | $8000$ | $0.4$ | 12.106
$10^{8}$ | $9000$ | $1.3$ | 0.120 |  | $10^{9}$ | $9000$ | $0.6$ | 2.589 |  | $10^{10}$ | $9000$ | $0.4$ | 12.161
$10^{8}$ | $10000$ | $1.4$ | 0.096 |  | $10^{9}$ | $10000$ | $0.6$ | 2.517 |  | $10^{10}$ | $10000$ | $0.4$ | 12.106
Table 5: Results of linear stability analysis. Wavelengths $\lambda_{max}$ and
exponential growth rates $\gamma_{max}$ of the fastest growing modes are shown
as functions of $Ha$ and $Gr$. Only the data for flow regimes identified as
Q2D in the analysis of the base flow are shown.
(a)${\lambda_{max}}$$Ha$(b)${\gamma_{max}}$$Ha$$Gr=10^{8}$$Gr=10^{9}$$Gr=10^{10}$$Gr=10^{6}$$Gr=10^{7}$$Gr=10^{8}$$Gr=10^{9}$
Figure 11: The wavelength $\lambda_{max}$ ($a$) and the exponential growth
rate $\gamma_{max}$ ($b$) of the fastest growing perturbations as a function
of $Ha$. The data are taken from table 5. The results of the current work and
that of Zhang & Zikanov (2014) are denoted as filled and unfilled elements,
respectively. Results of 2D approximation and 3D computational analysis are
shown, respectively, for $Gr=10^{8},10^{9}$ and $Gr=10^{10}$.

The spatial structure of the unstable modes is illustrated in figure 12 for
$Gr=10^{8},10^{9},10^{10}$ and $Ha=10^{4}$. The structures are qualitatively
similar to those found for Q2D instabilities by Zhang & Zikanov (2014).
Consistent with the first trend mentioned above and with the base flow
modification illustrated in figure 5 is the fact that the energy of growing
perturbations becomes contained to the lower part of the duct at higher values
of $Gr$.

(a)zx${\Theta}$(b)zx${\Theta}$(c)zx${\Theta}$(d)${\it U_{z}}$yx Figure 12:
Spatial structure of the fastest growing instability modes during the stage of
exponential growth at $Gr=10^{8}$ ($a$), $Gr=10^{9}$ ($b$), and $Gr=10^{10}$
($c$),($d$) for $Ha=10000$. Results of 2D approximation and 3D computational
analysis are shown, respectively, for $Gr=10^{8},10^{9}$ and $Gr=10^{10}$.
Perturbations of temperature and vector fields of velocity perturbations
($u^{\prime}$, $w^{\prime}$) in the vertical midplane ($y=0$) are shown in
($a$)-($c$). Perturbations of vertical velocity in the horizontal midplane
section $z=0$ are shown in ($d$). Solid and dashed isolines indicate positive
and negative values, respectively.

###  3.4 Results of DNS of nonlinear flows

The results concerning the nonlinear flow regimes are illustrated in figures
8($d,e,f$), 13, 14, and 15. DNS approach based on direct solution of the
nonlinear governing equations is utilized. The 2D SM82 model (14)$-$(16) and
the computational domain of length $L_{x}=4\pi$ are used for flows at
$Gr=10^{8}$ and $10^{9}$. Full 3D equations (2)$-$(4) are solved and the
domain is reduced to $L_{x}=2\pi$ at $Gr=10^{10}$. Other parameters of
computational model are described in section 2.5. Each simulation starts with
the streamwise-independent base flow (see section 3.2) computed at the same
$Gr$ and $Ha$, to which random small-amplitude ($\sim 10^{-3}$) random
perturbations of velocity and temperature are added.

$({\it a})$ $Gr=10^{8}$${\it E_{x}}$${\it t}$${\it Ha}=5\times 10^{3}$${\it
Ha}=10^{4}$$({\it b})$ $Gr=10^{9}$${\it t}$$({\it c})$ $Gr=10^{10}$${\it t}$
Figure 13: Time signals of the kinetic energy of the streamwise velocity
obtained in the DNS of flows at $Gr=10^{8}$ ($a$), $Gr=10^{9}$ ($b$), and
$Gr=10^{10}$ ($c$) for $Ha=5000$ and $Ha=10000$. Results of 2D approximation
and 3D computational analysis are shown, respectively, for $Gr=10^{8},10^{9}$
and $Gr=10^{10}$. (a)${\it z}$${\it u_{z}}$(b)${\it z}$${\it u_{x}}$(c)${\it
z}$${\theta}$(d)${\it z}$${\it u_{z}}$(e)${\it z}$${\it u_{x}}$(f)${\it
z}$x${\theta}$(g)$z$$u_{x}$(h)$z$$\theta$$Gr=10^{8}$$Gr=10^{9}$$Gr=10^{10}$
Figure 14: Flow structure in nonlinear regime at $Gr=10^{8}$ ($a$-$c$) and
$Gr=10^{9}$ ($d$-$f$) for $Ha=10000$. The instantaneous distributions of
$u_{z}$ ($a$, $d$), $u_{x}$ ($b$, $e$), and $\theta$ ($c$, $f$) obtained in
the 2D model are shown. The profiles of $u_{x}$ and $\theta$ obtained by
averaging over $x$ and time are shown, respectively, in ($g$) and ($h$) $({\it
a})$ $Ha=5\times 10^{3},Gr=10^{8}$$\theta$$t$$({\it b})$
$Ha=10^{4},Gr=10^{8}$$\theta$$t$$({\it c})$ $Ha=5\times
10^{3},Gr=10^{9}$$\theta$$t$$({\it d})$ $Ha=10^{4},Gr=10^{9}$$\theta$$t$$({\it
e})$ $Ha=5\times 10^{3},Gr=10^{10}$$\theta$$t$$({\it f})$
$Ha=10^{4},Gr=10^{10}$$\theta$$t$${\it z}=-0.2857$${\it z}=0$${\it z}=0.2857$
Figure 15: Time signals of temperature measured at top and bottom walls and in
the middle of the duct in fully developed flows at $Gr=10^{8},10^{9},10^{10}$
are shown for $Ha=5\times 10^{3}$ in ($a,c,e$) and for $Ha=10^{4}$ in
($b,d,f$). Results of 2D approximation and 3D computational analysis are
shown, respectively, for $Gr=10^{8},10^{9}$ ($a-d$) and $Gr=10^{10}$ ($e-f$).

The typical flow evolution is illustrated by the curves of average kinetic
energy shown in figure 13. The flow reaches a fully developed state after the
instability and initial development. The evolution of the fully developed flow
is computed for at least $500$ time units for the 2D model at
$Gr=10^{8},10^{9}$ and at least $100$ time units for the 3D model at
$Gr=10^{10}$. At this stage, the integral parameters fluctuate around steady
means (at $Gr=10^{9}$ and $10^{10}$) or remain steady (at $Gr=10^{8}$). The
amplitudes of the fluctuations are small at $Gr=10^{9}$ and large, but still
moderate at $Gr=10^{10}$.

Structure of fully developed flows is illustrated in figures 8 and 14. The
velocity field shows finite-amplitude roll-like structures (hereafter referred
to as rolls) resulting from the instability, which are superimposed on a
streamwise-independent mean flow (see figures 8$d$ and 14$a,d$). The rolls
cause variations of temperature (see figures 8$f$, 14$c$ and 14$f$). Transport
of the rolls by the mean flow is a known reason of MCFs in horizontal channels
(Zikanov et al., 2013; Zhang & Zikanov, 2014; Zikanov et al., 2021).

Comparison of the flow structures in figures 8$d$-$f$ and 14 reveal the effect
of the value of $Gr$ on convection rolls. As anticipated, increase of $Gr$
leads to higher non-dimensional amplitude of the velocity fluctuations. This
results in stronger vertical mixing as illustrated by the streamwise- and
time-averaged profiles in figure 14$g,h$. In particular, a nearly uniform
vertical distribution of average temperature with a narrow boundary layer at
the bottom is observed at $Gr=10^{10}$.

As we discussed earlier, MCFs caused by the instability have potentially
critical implications for design and operation of liquid-metal components of
nuclear fusion reactors. The DNS results allow us, for the first time, to
evaluate the properties of the MCFs at the high values of $Gr$ and $Ha$
corresponding to the actual reactor conditions.

In addition to the instantaneous temperature distributions in figures 8$f$,
14$c$ and 14$f$, the discussion will be based on the point-signals of
temperature measured at the top and bottom walls and in the middle of the duct
(see figure 15). As discussed, e.g., by Zikanov et al. (2021), measuring such
signals is the most reliable and commonly used tool for studying MCFs in
experiments.

The evident conclusion from the DNS data is that MCFs are fully present in
flows with $Gr=10^{8},10^{9},10^{10}$ and the highest values $Ha=5\times
10^{3}$ and $10^{4}$ considered in this study. The fluctuations are observed
in the entire duct. The temperature signals are regular and dominated by one
or several low frequencies (the typical period is $2$-$3$ non-dimensional time
units at $Gr=10^{8}$ and $10^{9}$). The signal is less regular and
characterized by higher typical frequencies at $Gr=10^{10}$.

Interestingly, the non-dimensional amplitude of the temperature fluctuations
decreases noticeably with growing $Gr$. Comparison of the signals in the two
columns of figure 15 demonstrates that value of $Ha$ has practically no effect
on the MCFs.

Considering the practical implications, it is interesting to evaluate the
parameters of the MCFs in dimensional units. We will do that for the
temperature signals at the bottom of the duct ($z=-0.2857$) assuming the duct
half-width $d=5$ $cm$ and using the physical properties of PbLi at $573$ $K$
(Zikanov et al., 2021). The wall heat rate is $q=10.56$ $kWm^{-2}$ at
$Gr=10^{8}$, $q=105.6$ $kWm^{-2}$ at $Gr=10^{9}$ and $q=1056$ $kWm^{-2}$ at
$Gr=10^{10}$. We find, by applying the temperature scale $qd/\kappa$, that the
largest amplitude of fluctuations of wall temperature is in the range $5-6$
$K$ at $Gr=10^{8}$, $44-62$ $K$ at $Gr=10^{9}$, and somewhat unrealistic
$180-300$ $K$ at $Gr=10^{10}$. The typical time period of the fluctuations is
$5.64$ s at $Gr=10^{8}$, $6.45$ $s$ at $Gr=10^{9}$, and $4.51$ $s$ at
$Gr=10^{10}$ for $Ha=10^{4}$.

Similar evaluations have been done for the future experiments on the recently
built experimental facility (see, e.g., Belyaev et al., 2017), in which liquid
mercury flows in the duct with the half-width of $d=2.8$ $cm$. The physical
properties of Hg are taken at $303$ $K$ (Zikanov et al., 2021). The wall heat
rate is $q=9.59$ $kWm^{-2}$ at $Gr=10^{8}$ and $q=95.9$ $kWm^{-2}$ at
$Gr=10^{9}$. The results of nonlinear simulations allow us to predict the
largest amplitude of fluctuations of temperature in the middle of the duct.
The amplitudes are in the range of $2-3$ $K$ at $Gr=10^{8}$ and $12-18$ $K$ at
$Gr=10^{9}$. The typical time period of the fluctuations is $1.1$ s at
$Gr=10^{8}$ and $1.58$ s at $Gr=10^{9}$ for $Ha=1000$.

##  4 Concluding remarks

We have analysed mixed convection in a liquid metal flow in a duct with bottom
heating and transverse magnetic field. The analysis is extended to much higher
values of $Ha$ and $Gr$ than the previous analysis of similar effects by Zhang
& Zikanov (2014).

The main conclusion of our work is that magnetoconvective fluctuations appear
at the parameters anticipated for operational regimes of blankets and
divertors of future fusion rectors. The fluctuations are not suppressed or
even significantly reduced in amplitude by the very strong magnetic field. The
amplitude remains high, reaching tens or hundreds (depending on the value of
$Gr$) $K$ in a typical duct geometry. This has significant far-reaching
implications for mixing, heat and mass transfer, and structural integrity of
reactor components. The most dangerous modes of the instability have the form
of rolls localized in the lower half of the duct and having the streamwise
wavelength measured in horizontal half-widths of the duct, approximately,
between $0.8$ and $1.4$ at $Gr=10^{8}$, $0.6$ at $Gr=10^{9}$, and $0.4$ at
$Gr=10^{10}$.

Another conclusion concerns applicability of the two-dimensional approximation
by Sommeria & Moreau (1982) to flows with thermal convection. We have found
that the approximation may become inaccurate at high values of $Gr$ even
though the flow remains quasi-two-dimensional. Reasons of this phenomenon
remain to be understood. The conclusion itself, however, is important as a
warning against application of the model without a proper verification.

Further study of MCFs at high Hartmann and Grashof numbers is warranted by
their practical importance and theoretical significance. Many interesting
possible directions of future work can be suggested. We mention two of them.
One is the exploration of the phenomenon for other geometries, where strong
MCFs are known to exist, for example for downward flow in a vertical duct.
Another particularly interesting direction is the analysis of the effects of
finite thermal and electrical conductivities of the walls.

Authors are thankful to D. Krasnov for continuing assistance with the
numerical model and S. Molokov for interesting and useful discussions. Work of
R. Akhmedagaev and O. Zikanov is supported by the US NSF (Grant CBET 1803730
"Extreme magnetoconvection"). Work of Y. Listratov is supported by the
Ministry of Science and Education of the Russian Federation (Grant
14.Z50.31.0042) and by the Russian Foundation for Basic Research (Grant NNIO
18-508-12005).

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

## References

  * Abdou et al. (2015) Abdou, M., Morley, N. B., Smolentsev, S., Ying, A., Malang, S., Rowcliffe, A. & Ulrickson, M. 2015 Blanket/first wall challenges and required R&D on the pathway to DEMO. Fusion Eng. Des. 100, 2–43.
  * Alboussière et al. (1993) Alboussière, T., Garandet, J. P. & Msoreau, R. 1993 Buoyancy-driven convection with a uniform magnetic field. Part 1. asymptotic analysis. J. Comp. Phys. 253, 545–563.
  * Belyaev et al. (2015) Belyaev, I., Ivochkin, Y. P., Listratov, Y. I., Razuvanov, N. & Sviridov, V. 2015 Temperature fluctuations in a liquid metal MHD-flow in a horizontal inhomogeneously heated tube. High Temp. 53 (5), 734–741.
  * Belyaev et al. (2021) Belyaev, I., Sardov, P., Melnikov, I. & Frick, P. 2021 Limits of strong magneto-convective fluctuations in liquid metal flow in a heated vertical pipe affected by a transverse magnetic field. Int. J. Therm. Sci. 101, 106773.
  * Belyaev et al. (2017) Belyaev, I. A., Sviridov, V. G., Batenin, V. M., Biryukov, D. A., Nikitina, I. S., Manchkha, S. P., Pyatnitskaya, N. Yu., Razuvanov, N. G. & Sviridov, E. V. 2017 Test facility for investigation of heat transfer of promising coolants for the nuclear power industry. Therm. Eng. 64 (11), 841–848.
  * Branover (1978) Branover, H. 1978 Magnetohydrodynamic Flow in Ducts. New York: Wiley.
  * Davidson (2016) Davidson, P. A. 2016 Introduction to Magnetohydrodynamics, 2nd Ed. Cambridge University Press.
  * Gelfgat & Molokov (2011) Gelfgat, A. Yu. & Molokov, S. 2011 Quasi-two-dimensional convection in a three-dimensional laterally heated box in a strong magnetic field normal to main circulation. Phys. Fluids 23, 034101.
  * Gelfgat & Zikanov (2018) Gelfgat, A. Yu. & Zikanov, O. 2018 Computational modeling of magnetoconvection: effects of discretization method, grid refinement and grid stretching. Comp. Fluids 175, 66–82.
  * Genin et al. (2011) Genin, L. G., Zhilin, V. G., Ivochkin, Y. P., Razuvanov, N. G., Belyaev, I. A., Listratov, Y. I. & Sviridov, V. G. 2011 Temperature fluctuations in a heated horizontal tube affected by transverse magnetic field. In Proc. $8^{th}$ PAMIR Conf. Fund. Appl. MHD, pp. 37–41. Borgo, Corsica, France.
  * Kirillov & Muraviev (1997) Kirillov, I.R. & Muraviev, E.V. 1997 Review of liquid metal divertor concepts for tokomak reactors. Fusion Technology 1996 pp. 251–254.
  * Krasnov et al. (2011) Krasnov, D., Zikanov, O. & Boeck, T. 2011 Comparative study of finite difference approaches to simulation of magnetohydrodynamic turbulence at low magnetic Reynolds number. Comp. Fluids 50, 46–59.
  * Krasnov et al. (2012) Krasnov, D., Zikanov, O. & Boeck, T. 2012 Numerical study of magnetohydrodynamic duct flow at high Reynolds and Hartmann numbers. J. Fluid Mech. 704, 421–446.
  * Listratov et al. (2018) Listratov, Y. I., Ognerubov, D., Zikanov, O. & Sviridov, V. 2018 Numerical simulations of mixed convection in liquid metal flow within a horizontal pipe with transverse magnetic field. Fluid Dyn. Res. 50 (5), 051407.
  * Lyubimova et al. (2009) Lyubimova, T. P., Lyubimov, D. V., Morozov, V. A., Scuridin, R. V., Hadid, H. B. & Henry, D. 2009 Stability of convection in a horizontal channel subjected to a longitudinal temperature gradient. Part 1. Effect of aspect ratio and Prandtl number. J. Fluid Mech. 635, 275–295.
  * Ni et al. (2007) Ni, M.-J., Munipalli, R., Huang, P., Morley, N. B. & Abdou, M. A. 2007 A current density conservative scheme for incompressible MHD flows at a low magnetic Reynolds number. Part I: On a rectangular collocated grid system. J. Comp. Phys. 227, 174–204.
  * Ozoe (2005) Ozoe, H. 2005 Magnetic Convection. Imperial College Press.
  * Pothérat (2007) Pothérat, A. 2007 Numerical simulations of an effective two-dimensional model for flows with a transverse magnetic field. Phys. Fluids 19 (7), 074104.
  * Pothérat et al. (2000) Pothérat, A., Sommeria, J. & Moreau, R. 2000 An effective two-dimensional model for MHD flows with transverse magnetic field. J. Fluid Mech. 424, 75–100.
  * Pothérat et al. (2005) Pothérat, A., Sommeria, J. & Moreau, R. 2005 Numerical simulations of an effective two-dimensional model for flows with a transverse magnetic field. J. Fluid Mech. 534, 155–143.
  * Sahu et al. (2020) Sahu, S., Courtessole, C., Ranjan, A., Bhattacharyay, R., Sketchley, T. & Smolentsev, S. 2020 Thermal convection studies in liquid metal flow inside a horizontal duct under the influence of transverse magnetic field. Phys. Fluids 32 (6), 067107\.
  * Smolentsev et al. (2008) Smolentsev, S., Moreau, R. & Abdou, M. 2008 Characterization of key magnetohydrodynamic phenomena for PbLi flows for the US DCLL blanket. Fusion Eng. Des. 83, 771–783.
  * Smolentsev et al. (2010) Smolentsev, S., Moreau, R., Bühler, L. & Mistrangelo, C. 2010 MHD thermofluid issues of liquid-metal blankets: Phenomena and advances. Fusion Eng. Des. 85 (7-9), 1196–1205.
  * Smolentsev et al. (2012) Smolentsev, S., Vetcha, N. & Moreau, R. 2012 Study of instabilities and transitions for a family of quasi-two-dimensional magnetohydrodynamic flows based on a parametrical model. Phys. Fluids 24, 024101.
  * Sommeria & Moreau (1982) Sommeria, J. & Moreau, R. 1982 Why, how and when MHD-turbulence becomes two-dimensional. J. Fluid Mech. 118, 507–518.
  * de les Valls et al. (2008) de les Valls, E. Mas, Sedano, L., Batet, L., Ricapito, I., Aiello, A., Gastaldi, O. & Gabriel., F. 2008 Lead–lithium eutectic material database for nuclear fusion technology. J. Nucl. Mater. 376, 353–357.
  * Vetcha et al. (2013) Vetcha, N., Smolentsev, S., Abdou, M. & Moreau, R. 2013 Study of instabilities and quasi-two-dimensional turbulence in volumetrically heated magnetohydrodynamic flows in a vertical rectangular duct. Phys. Fluids 25 (2), 024102.
  * Vo et al. (2017) Vo, T., Pothérat, A. & Sheard, G. J. 2017 Linear stability of horizontal, laminar fully developed, quasi-two-dimensional liquid metal duct flow under a transverse magnetic field and heated from below. Phys. Rev. Fluids 2, 033902\.
  * Zhang & Zikanov (2014) Zhang, X. & Zikanov, O. 2014 Mixed convection in a horizontal duct with bottom heating and strong transverse magnetic field. J. Fluid Mech. 757, 33–56.
  * Zhang & Zikanov (2017) Zhang, X. & Zikanov, O. 2017 Thermal convection in a toroidal duct of a liquid metal blanket. part ii. effect of axial mean flow. Fusion Eng. Des. 116, 40–46.
  * Zhang & Zikanov (2018) Zhang, X. & Zikanov, O. 2018 Convection instability in a downward flow in a vertical duct with strong transverse magnetic field. Phys. Fluids 30, 117101.
  * Zhao & Zikanov (2012) Zhao, Y. & Zikanov, O. 2012 Instabilities and turbulence in magnetohydrodynamic flow in a toroidal duct prior to transition in Hartmann layers. J. Fluid Mech. 692, 288–316.
  * Zikanov (2019) Zikanov, O. 2019 Essential Computational Fluid Dynamics, 2nd Ed. Wiley.
  * Zikanov et al. (2014) Zikanov, O., Krasnov, D., Boeck, T., Thess, A. & Rossi, M. 2014 Laminar-turbulent transition in magnetohydrodynamic duct, pipe, and channel flows. Appl. Mech. Rev. 66 (3), 030802.
  * Zikanov et al. (2021) Zikanov, O., Listratov, Y., Razuvanov, N., Belyaev, I., Frick, P. & Sviridov, V. 2021 Mixed convection in pipe and duct flows with strong magnetic fields. Appl. Mech. Rev. 73 (1), 010801.
  * Zikanov et al. (2013) Zikanov, O., Listratov, Ya. & Sviridov, V. G. 2013 Natural convection in horizontal pipe flow with strong transverse magnetic field. J. Fluid Mech. 720, 486–516.


11institutetext: Department of Physical Sciences,
Indian Institute of Science Education and Research (IISER) - Mohali,
Knowledge City, Sector 81, SAS Nagar, Punjab 140306, India

# Lattice Supersymmetry and Holography

Anosh Joseph 11<EMAIL_ADDRESS>
###### Abstract

Over the last twenty years, work based on lattice supersymmetry has generated
many new results and insights into the non-perturbative nature of string
theory, quantum black holes, and gravity. This endeavor is a broad research
program encompassing lattice field theory, supersymmetry, string theory, and
quantum gravity. In this volume, we look at a selected subset of the topics
covering recent progress in lattice supersymmetry and holography.

##  1 Introduction

Supersymmetric theories play a prominent role in our efforts to understand
string theory, strong dynamics, quantum gravity, and various extensions to the
Standard Model of particle physics. The holographic duality conjecture reveals
the connections between theories of quantum gravity living on curved
spacetimes and quantum field theories without gravity living on the boundaries
of such spacetimes. This conjecture also provides promising directions for
studying the nature of quantum gravity and black holes.

It is possible to describe certain black hole geometries in terms of the
world-volume theories of the D-branes that compose them. These are the
maximally supersymmetric Yang-Mills theories in various spacetime dimensions,
with many colors, taken in the ‘t Hooft limit and at finite temperatures.
These theories can be strongly coupled in the regimes in which they describe
string theory backgrounds, including black holes. Solving these field theories
at strong coupling and finite temperature would allow us to directly study the
quantum properties of the dual black holes, including their thermodynamic
features. As the ‘t Hooft coupling and the number of colors are reduced,
classical and quantum string corrections become more prominent. These less
understood limits can be investigated using a lattice formulation - a first-
principles definition of quantum field theories. Simulations of lattice
discretized field theories can validate the holographic duality conjecture and
provide new insight into the non-perturbative structure of string theory and
quantum gravity.

We hope that this special issue would be beneficial to beginning researchers
and practitioners in string theory, quantum field theory, and lattice field
theory who would want to contribute to this exciting interdisciplinary topic
of lattice supersymmetry and its usefulness in testing and validating the
holographic duality conjecture.

##  2 A Brief Overview of Contributions

Below we provide a brief overview of the contributions within this special
issue.

In Ref. Schaich:2022xgy , David Schaich reviews the recent progress and near
future prospects in lattice investigations of supersymmetric field theories
and some of the challenges that remain to be overcome. He focuses on the
progress in three areas: supersymmetric Yang-Mills (SYM) theories in fewer
than four spacetime dimensions, four-dimensional ${\cal N}=1$ SYM theory and
maximally supersymmetric Yang-Mills theory in four dimensions. He also
highlights supersymmetric QCD (SQCD) and the sign problem as significant
challenges that will be important to address in future work.

In Ref. Asano:2022jag , Yuhma Asano outlines some basics of the matrix model
conjecture and the gauge/gravity duality conjecture for the matrix models. He
reviews various numerical evidence provided for the gauge/gravity duality
conjecture for the BFSS and BMN matrix models and their flavored cousins, the
Berkooz-Douglas (BD) and Kim-Lee-Yi (KLY) matrix models.

Masanori Hanada and Hiromasa Watanabe review the basic properties of partial
deconfinement in Ref. Hanada:2022wcq  and discuss its applications. The
confinement deconfinement transition in gauge theory plays an important role
in physics, including describing thermal phase transitions in the dual
gravitational theory. In the scenario of partial deconfinement, there is an
intermediate phase in which the color degrees of freedom split into the
confined and deconfined sectors. The partially deconfined phase is dual to the
small black hole that lies between the large black hole and graviton gas. A
better understanding of partial deconfinement may explain how gravity emerges
from the degrees of freedom of the field theory.

In Ref. Bliard:2022 , Gabriel Bliard, Ilaria Costa, and Valentina Forini
review the lattice study of the Green-Schwarz gauge-fixed string action
describing the worldsheet fluctuations about the minimal surface
holographically dual to the null cusp Wilson loop. A numerical study of this
system using the Monte Carlo method helps evaluate the cusp anomaly of ${\cal
N}=4$ super Yang-Mills. They comment on the discretization, numerical
explorations, and challenges for the non-perturbative study of this benchmark
model of gauge-fixed worldsheet actions.

In Ref. Jha:2022 , Raghav Jha proposes additional tests of holography by
studying supersymmetric Wilson loops in $p+1$-dimensional maximally
supersymmetric Yang-Mills (SYM) theories on a lattice, taken in the large-$N$
limit. In the dual gravity description, this computation involves calculating
the area of a fundamental string worldsheet in certain Type II supergravity
backgrounds. Even though thermodynamic observables have been computed on the
lattice using Monte Carlo methods, and the results agree with the supergravity
results in various dimensions, more needs to be done for the gauge-invariant
operators, such as the Wilson loop. Jha provides analytical predictions for
these loops for various non-conformal D$p$-brane background cases, with $p<3$,
in the large $N$ limit. He also comments on how these can be computed on non-
orthogonal lattices for various supersymmetric models.

Daisuke Kadoh and Naoya Ukita in Ref. Kadoh:2022  propose a supersymmetric
gradient flow for four-dimensional ${\cal N}=1$ SQCD. They gave expressions
for the flow equation in the superfield formalism and the component fields
formalism in the Wess-Zumino gauge. They also discuss a simplified flow using
the gradient of supersymmetric Yang-Mills (SYM) action instead of SQCD action
to define a gauge multiplet flow.

##  3 Future Research Directions

In this section, we briefly outline the various research directions the effort
of lattice supersymmetry and holography can take in the near future.

  * •

For the case of the ${\cal N}=4$ Yang-Mills in four dimensions, there is a
famous holographic prediction for the Coulomb coefficient $C(\lambda)$ that it
is proportional to $\sqrt{\lambda}$ up to ${\cal O}(1/\sqrt{\lambda})$
corrections. For the $N=\infty$ planar limit more general analytic results
have also been obtained. It would be interesting to search for this behavior
in more detail, although some promising preliminary results have already been
reported.

  * •

A more detailed understanding of the non-trivial scaling dimension of the
simplest conformal primary operator of the four-dimensional ${\cal N}=4$ SYM
theory, the Konishi operator is much needed. Preliminary lattice results have
already been obtained and are consistent with existing perturbation theory
results.

  * •

It would be interesting to study the behavior of the four-dimensional ${\cal
N}=4$ SYM theory around the $S$-dual point, where the ‘t Hooft coupling takes
the form $\lambda_{\rm sd}=4\pi N$.

  * •

Another direction is to adjust the scalar potential to study the four-
dimensional ${\cal N}=4$ SYM on the Coulomb branch of the moduli space. In
this context, the $S$-duality connects the masses of the U(1)-charged
elementary ‘$W$ bosons’ and the magnetically charged topological ‘t Hooft-
Polyakov monopoles. They can be accessed from lattice calculations with
appropriate boundary conditions.

  * •

Non-perturbative lattice calculations can be used to study the free energy of
four-dimensional ${\cal N}=4$ SYM theory. The weak-coupling perturbative
prediction and the strong coupling holographic calculation differ by a famous
factor of $3/4$. A lattice setup can be used to interpolate between these two
coupling regimes.

  * •

In the case of maximally supersymmetric Yang-Mills in three dimensions, one
can explore, through lattice investigations, the phase transition between the
‘D2 phase’ and the spatially deconfined ‘D0 phase’ dual to a localized black
hole geometry.

  * •

A more detailed and careful non-perturbative investigation of the finite
temperature phase diagrams in two-dimensional SYM theories is still needed.
These theories also possess rich zero-temperature dynamics, such as the
‘meson’ spectrum and spontaneous symmetry breaking, that are important to
explore non-perturbatively.

  * •

In the context of the BMN model, the critical temperature of the confinement
transition can be predicted by perturbative calculations in the weak-coupling
regime and by a dual supergravity calculation for strong coupling. An open
research direction would be non-perturbatively connecting the transition
temperatures by mapping out the intermediate coupling regime, where
perturbative and holographic approaches are unreliable.

  * •

There has been excellent progress in validating and testing holography in the
context of BFSS and BMN matrix models. It would be interesting to extend these
investigations into more exotic models such as the KLY matrix model.

  * •

Partial deconfinement can be understood as the coexisting phenomenon in the
space of the color degrees of freedom. Historically, it was proposed for the
SYM theory as the dual of the small black hole phase. An open research
direction would be to generalize the idea of partial deconfinement to systems
with finite $N$ by utilizing their chiral symmetry.

We hope this volume will contribute to constructive discussions on advances in
lattice supersymmetry and holography and stimulate further studies.

### Acknowledgements

The author was supported in part by the Start-up Research Grant (No. SRG /
2019 / 002035) from the Science and Engineering Research Board (SERB),
Government of India, and in part by the Indian Institute of Science Education
and Research (IISER) - Mohali.

### Data Availability Statement

No Data associated in the manuscript.

## References

  * (1) D. Schaich, Lattice studies of supersymmetric gauge theories, Eur. Phys. J. Spec. Top. (2022). doi:10.1140/epjs/s11734-022-00708-1
  * (2) Y. Asano, Numerical development of the matrix model for non-perturbative formulation of M-theory, Eur. Phys. J. Spec. Top. (2022). doi:10.1140/epjs/s11734-022-00718-z
  * (3) M. Hanada and H. Watanabe, Partial deconfinement: a brief overview, Eur. Phys. J. Spec. Top. (2022). doi:10.1140/epjs/s11734-022-00709-0
  * (4) G. Bliard, I. Costa and V. Forini, Holography on the lattice: the string worldsheet perspective, Eur. Phys. J. Spec. Top. (2023). doi: 10.1140/epjs/s11734-023-00769-w
  * (5) R. G. Jha, Supersymmetric Wilson loops on the lattice in the large-$N$ limit, Eur. Phys. J. Spec. Top. (2023). doi: 10.1140/epjs/s11734-023-00768-x
  * (6) D. Kadoh and N. Ukita, Supersymmetric gradient flow in $4d$ ${\cal N}=1$ SQCD, Eur. Phys. J. Spec. Top. (2022). doi:10.1140/epjs/s11734-022-00707-2