Wall-attached convection under strong inclined magnetic fields

Abstract We employ a linear stability analysis and direct numerical simulations to study the characteristics of wall modes in thermal convection in a rectangular box under strong and inclined magnetic fields. The walls of the convection cell are electrically insulated. The stability analysis assumes periodicity in the spanwise direction perpendicular to the plane of a homogeneous magnetic field. Our study shows that for a fixed vertical magnetic field, the imposition of horizontal magnetic fields results in an increase of the critical Rayleigh number along with a decrease in the wavelength of the wall modes. The wall modes become tilted along the direction of the resulting magnetic fields and therefore extend further into the bulk as the horizontal magnetic field is increased. Once the modes localized on the opposite walls interact, the critical Rayleigh number decreases again and eventually drops below the value for onset with a purely vertical field. We find that for sufficiently strong horizontal magnetic fields, the steady wall modes occupy the entire bulk and therefore convection is no longer restricted to the sidewalls. The aforementioned results are confirmed by direct numerical simulations of the nonlinear evolution of magnetoconvection. The direct numerical simulation results also reveal that at least for large values of horizontal magnetic field, the wall-mode structures and the resulting heat transfer are dependent on the initial conditions.


Introduction
Buoyancy-driven flows of electrically conducting fluids under the influence of magnetic fields are a common occurrence in geophysical, astrophysical, as well as in several technological applications.Such flows are called magnetoconvection and their driving mechanism is the temperature dependence of the fluid density, which results in spatial density variations leading to buoyancy forces acting on the fluid.When such a fluid moves under the influence of magnetic fields, electric currents are induced in the fluid due to Faraday's law, which, in turn, induce magnetic fields by the virtue of Ampere's law.These electric currents interact with the applied and induced magnetic fields to generate a Lorentz force distribution that acts on the fluid.Therefore, magnetoconvective flows are governed by equations of conservation of mass, momentum, and thermal energy, along with Maxwell's equations for electromagnetism and Ohm's law (Weiss & Proctor 2014).Magnetoconvection is encountered in the Sun, stars, and planetary dynamos.In industries and technological applications, magnetoconvection is typically encountered in liquid-metal batteries (Kelley & Sadoway 2014;Shen & Zikanov 2016;Kelley & Weier 2018), cooling liquid-metal blankets in fusion reactors (Mistrangelo et al. 2020(Mistrangelo et al. , 2021)), and magnetic stirring and braking of liquid metal melts (Davidson 1999(Davidson , 2017;;Lyubimov et al. 2010).
A simplified paradigm for magnetoconvection consists of a fluid layer that is heated from below and cooled from above (Rayleigh-Bénard convection or RBC) with imposed magnetic fields in different configurations.Typically, the Boussinesq approximation is employed for modelling the above flows; this approximation assumes that the flow is incompressible and the density variations are negligible except in the buoyancy term in the momentum equation (Chandrasekhar 1981;Lohse & Xia 2010;Chillà & Schumacher 2012;Verma 2018).Magnetoconvection is governed by the following nondimensional parameters: i) Rayleigh number Ra -the ratio of buoyancy to dissipative forces, ii) Prandtl number Pr -the ratio of kinematic viscosity to thermal diffusivity, iii) Hartmann number Ha -the ratio of Lorentz to viscous forces, and iv) the magnetic Prandtl number Pm -the ratio of kinematic viscosity to magnetic diffusivity.The important nondimensional output parameters of magnetoconvection are i) the Nusselt number Nu -the ratio of the total heat transport to the diffusive heat transport, ii) the Reynolds number Re -the ratio of inertial to viscous forces, and iii) the magnetic Reynolds number Rm -the ratio of induction to diffusion of the magnetic field.In liquid-metal convection typically encountered in laboratory experiments and most industrial applications, the magnetic Reynolds number is sufficiently small such that the induced magnetic field is negligible compared to the applied magnetic field and is thus neglected in the expressions of the Lorentz force and Ohm's law (Roberts 1967;Davidson 2017;Verma 2019).Such cases are referred to as quasi-static magnetoconvection where the induced magnetic field adjusts instantaneously to the changes in velocity.In the quasi-static approximation, there exists a one-way influence of the magnetic field on the flow only.
Magnetoconvection has been studied theoretically in the past (for example, Chandrasekhar 1981;Houchens et al. 2002;Busse 2008) as well as with the help of experiments (for example, Nakagawa 1957;Fauve et al. 1981;Cioni et al. 2000;Aurnou & Olson 2001;Burr & Müller 2001;King & Aurnou 2015;Vogt et al. 2018Vogt et al. , 2021;;Zürner et al. 2020;Grannan et al. 2022) and numerical simulations (for example, Liu et al. 2018;Yan et al. 2019;Akhmedagaev et al. 2020a,b;Nicoski et al. 2022;Bhattacharya et al. 2023).An application of horizontal magnetic fields causes the large-scale rolls to become quasi two-dimensional and align in the direction of the field (Fauve et al. 1981;Busse & Clever 1983;Burr & Müller 2002;Yanagisawa et al. 2013;Tasaka et al. 2016;Vogt et al. 2018Vogt et al. , 2021)).These self-organized flow structures reach an optimal state wherein the heat transport and convective velocities increase significantly compared to convection without magnetic fields (Vogt et al. 2021).In contrast, strong vertical magnetic fields suppress convection (Chandrasekhar 1981;Cioni et al. 2000;Zürner et al. 2020;Akhmedagaev et al. 2020a,b).It must be noted that in a Rayleigh-Bénard system, convection commences only above a certain critical Rayleigh number, which is Ra c ≈ 1708 for the case with infinite no-slip horizontal walls (Chandrasekhar 1981).For Ra < Ra c , the heat transfer occurs purely by diffusion.The critical Rayleigh number increases when a vertical magnetic field is imposed and scales as Ra c ∼ Ha 2 in the asymptotic limit of large Hartmann numbers.
The dynamics of convection under strong vertical magnetic fields become more intricate in the presence of sidewalls.Houchens et al. (2002) and Busse (2008) analytically showed that magnetoconvection near the sidewalls ceases at Rayleigh numbers below the ones required to completely suppress convection in the bulk.Several numerical and experimental studies on magnetoconvection with sidewalls have also revealed the presence of residual wall-attached convection at Ra < Ra c (Houchens et al. 2002;Liu et al. 2018;Akhmedagaev et al. 2020a,b;Zürner et al. 2020;Teimurazov et al. 2023;McCormack et al. 2023).These so-called wall modes were shown to exhibit a two-layered structure and become more closely attached to the sidewalls with the increase of Hartmann number (Liu et al. 2018).
There are a few studies only on convection with inclined magnetic fields which motivates the present work (Hurlburt et al. 1996;Nicoski et al. 2022).The results of Hurlburt et al. (1996) indicate that the mean flows tend to travel in the direction of the tilt.Nicoski et al. (2022) observed qualitative similarities between convection with inclined magnetic field and that with vertical magnetic field in terms of the behavior of convection patterns, heat transport, and flow speed.However, to the best of our knowledge, there are no studies for the case with inclined magnetic fields where the Rayleigh number is close to but less than the critical Rayleigh number.Therefore, in the present work, we study thermal magnetoconvection in the wall-attached convection regime and explore the effects of additional horizontal magnetic fields on the wall modes.We use a combination of linear stability analysis and direct numerical simulations to study the dependence of the horizontal magnetic field strength, relative to the vertical magnetic field, on the wall-mode structures and their impact on large-scale heat and momentum transport.
The outline of the paper is as follows.In § 2, we discuss the problem setup, linear stability model, and the schemes for direct numerical simulations.The linear stability analysis and the results of direct numerical simulations are described in § 3. We conclude in § 4.

Numerical model
In this section, we discuss the mathematical model of our problem and the numerics employed for the stability analysis and direct numerical simulations.The study will be conducted under the quasi-static approximation, in which the induced magnetic field is neglected as it is very small compared to the applied magnetic field.This approximation is fairly accurate for magnetoconvection in liquid metals (Davidson 2017).Further, we employ Boussinesq approximation, in which the variations in the density of the fluid are ignored except in the buoyancy term in the momentum equation.Hence, the flow is essentially treated as incompressible.The governing equations of magnetoconvection under the above approximations are as follows: where u, j, p, T , and ϕ are the fields of velocity, current density, pressure, temperature, and electrical potential respectively, and B = (B x , B y , B z ) is the applied magnetic field.
In our work, the magnetic fields are inclined along y-direction only, hence B x = 0. Further, ν is the kinematic viscosity, κ is the thermal diffusivity, ρ is the density, and σ is the electrical conductivity of the fluid.The last term in the momentum equation (2.2) is the Lorentz force density.Equation (2.4) is Ohm's law.The Poisson equation (2.5) for the electric potential is a consequence of the charge conservation condition ∇ • j = 0.
In the following, we will discuss how the above equations have been employed for our linear stability analysis and direct numerical simulations.

Linear stability model
We first discuss the derivation of the perturbation equations for our linear stability analysis.The equations (2.1) to (2.5) are non-dimensionalized using the cell height H as the length scale, κ/H as the velocity scale, the temperature difference ∆ between the two horizontal plates as the temperature scale, and B z , the vertical component of the applied magnetic field.For this part of the analysis, we take the units that are typically chosen for a linear stability analysis to end with a Prandtl number-independent set of equations at the marginal stability threshold.The characteristic units in the subsequent simulation part will differ.The non-dimensionalized governing equations are as follows. (2.9) (2.10) In the above system of equations, Ra is the Rayleigh number, Pr is the Prandtl number, Ha z is the Hartmann number based on the vertical component of the magnetic field, and R is the ratio of the horizontal to the vertical magnetic field strength.These quantities are the governing parameters for our setup and are given by (2.11) The quantity θ is the difference between the temperature and the linear conduction profile, i.e., T (x, t) = θ(x, t) − z . (2.12) Apart from Ra, Pr , Ha z , and R, the dynamics are also governed by the aspect ratio Γ , which is the ratio of the length to the height of the convection cell.We examine a total of 12 cases of magnetoconvection in a bounded, horizontallyextended domain of dimension Γ × 1 in the y-z plane.Two aspect ratios are considered: Γ = 2 and Γ = 4.For Γ = 2, we consider the cases of Ha z = 50 and Ha z = 100, whereas for Γ = 4, we consider only the case of Ha z = 50.For each Ha z analysed in our study, we vary R from 0 (corresponding to a purely vertical magnetic field) to 3 in steps of 1.The convection cell is periodic in x-direction and consists of no-slip horizontal walls at z = ±1/2 and two no-slip sidewalls at y = ±Γ/2.All the walls are electrically insulated.Each horizontal wall is at a constant temperature, with the bottom wall at T = 0.5 and the top wall at T = −0.5.The sidewalls are thermally insulated with ∂T /∂η = 0, where η is the direction normal to the wall.A sketch of our setup is shown in figure 1.
In the present work we assume that the instability is of stationary type.The nonlinear terms as well as the time derivatives in equations (2.6) and (2.7) are therefore neglected.The momentum and continuity equations reduce to the Stokes problem with additional buoyancy and Lorentz force.To avoid complications stemming from the coupling between pressure and velocity we choose a representation for the velocity that satisfies the continuity equation automatically and eliminate the pressure term.The velocity field is written as the curl of a vector streamfunction ψ u = ∇ × ψ, (2.13) and the gauge condition ∇ • ψ = 0 is imposed to determine ψ uniquely as in Priede et al. (2010).The dependence on x is represented by the normal mode ansatz with wavenumber β for all fields, e.g.θ(y, z) exp(iβx) for the temperature perturbation.The gauge condition allows one to express the x-component of ψ by (2.14) The velocity components then read (2.15) (2.16) (2.17) Equations for ψ y and ψ z are obtained by taking the curl of the definition (2.13) and the momentum equation (2.6).They are (2.21) The quantities ω y and ω z are the y-and z-components of the vorticity field ∇ × u.
Equations for the remaining quantities are where N y + 1 and N z + 1 are the number of expansion terms with respect to y and z.
The boundary conditions for the vector stream function and vorticity components are determined with the help of equations (2.15)-(2.17).Zero normal velocity on the horizontal walls requires ψ y = 0 and ∂ z ψ z = 0. On the y=±Γ/2 sidewalls, the corresponding conditions are ψ z = 0 and ∂ y ψ y = 0.The tangential velocity vanishes on the sidewalls if ω y = 0 and ∂ y ψ z = ∂ z ψ y .On the top and bottom walls these conditions are ω z = 0 and ∂ y ψ z = ∂ z ψ y .The remaining boundary conditions for (2.22) are θ = 0 on the top and bottom walls and ∂ y θ = 0 on the y=±Γ/2 sidewalls.The boundary condition for the electric potential supplementing equation (2.23) is the homogeneous Neumann condition.
Since the boundary conditions (zero normal velocity) for (2.18) and (2.19)only involve ψ y and ψ z , respectively, one can represent the expansion coefficients of ψ y and ψ z by linear invertible maps through those of ω y and ω z (assuming the latter are augmented by the zero boundary values to be imposed on either ψ y , ψ z or its normal derivatives).The expansions for ψ y and ψ z therefore contain N y + 3 and N z + 3 terms, respectively.The values of ψ y , ψ z or its derivatives in equations (2.20)-(2.23)(and associated boundary conditions) at the collocation points are represented through expansion coefficients of ω y and ω z via these linear invertible maps.The same can be done for the electric potential, which is the sum of a contribution from ω z and ω y .As a result of the collocation approximation one obtains a vector Y of unknowns containing the expansion coefficients of ω y , ω z , θ with a size of 3(N y + 1)(N z + 1) and a generalized linear eigenvalue problem AY = Ra BY . (2.26) The method was implemented in Matlab (The MathWorks Inc. 2022) using the default double precision.Notice that ω y , ω z and ϕ are real variables.According to equations (2.15-2.17)and (2.20-2.22),u y , u z and θ would be purely imaginary quantities.They are considered as real variables in the code and below.Problem (2.26) was solved with Matlab's eig routine to find all eigenvalues and eigenvectors.The routine also works with a matrix B whose rank is smaller than the rank of A (as it is the case for (2.26)).It associates the spurious solutions that stem from equations not containing the eigenvalue Ra with infinite eigenvalues.
For the cases studied in this paper, the numerical resolution was typically N y = 70, N z = 60 for Ha z = 50 and N y = 90, N z = 70 for Ha z = 100 with aspect ratio Γ = 2.A decrease of N y by 10 typically only resulted in a relative change of the first eigenvalue below 10 −5 .The generation of the matrices and the solution of problem (2.26) for a given wavenumber took about 20 hours for N y = 90, N z = 70 and about 6 hours for N y = 70, N z = 60 on an Intel Xeon E5 CPU.

Direct numerical simulations
We conduct direct numerical simulations (DNS) of our magnetoconvection setup using a second-order finite difference code developed by Krasnov et al. (2011).The governing equations are made dimensionless by using the cell height H, the imposed temperature difference ∆, and the free-fall velocity U = √ αg∆H (where g and α are respectively the gravitational acceleration and the volumetric coefficient of thermal expansion of the fluid).The following non-dimensional equations are employed for our DNS: (2.30) We numerically solve equations (2.27) to (2.31).The Prandtl number Pr is chosen to be 0.025, which is the same as that of mercury.We choose two Hartmann numbers based on vertical magnetic field: Ha z = 100 and Ha z = 200.For these Hartmann numbers, the corresponding critical Rayleigh numbers (Ra c ) for the case without horizontal walls obtained using the linear stability analysis of Chandrasekhar (1981) is: Since we are interested in the wall-attached convection regime, we choose the Rayleigh number to be slightly below the critical Rayleigh numbers given by (2.32).Thus, the chosen Rayleigh numbers are Ra = 10 5 for the case of Ha z = 100 and Ra = 4 × 10 5 for the case of Ha z = 400.For each Ha z , we vary R from 0 to 3. We choose the domain-size to be Γ × Γ × H = 4 × 4 × 1 and employ a grid-resolution ranging from 1200 × 1200 × 300 points to 1800 × 1800 × 400.The horizontal walls are at z = ±1/2 and the sidewalls are at x = ±Γ/2 and y = ±Γ/2.The mesh is nonuniform with stronger clustering of the grid points near the top and bottom boundaries.The elliptic equations for pressure, electric potential, and the temperature are solved based on applying cosine transforms in x-and y-directions and using a tridiagonal solver in the z-direction.The diffusive term in the temperature transport equation is treated implicitly.The time discretization of the momentum equation uses the fully explicit Adams-Bashforth/Backward-Differentiation method of second order (Peyret 2002).A constant time step size ranging from 5 × 10 −5 to 1 × 10 −4 free fall time unit was chosen for our simulations, which satisfied the Courant-Friedrichs-Lewy (CFL) condition for our runs.
All the walls are rigid and electrically insulated such that the electric current density j forms closed field lines inside the cell.The top and bottom walls are held fixed at T = −0.5 and T = 0.5 respectively, and the sidewalls are adiabatic with ∂T /∂η = 0 (where η is the component normal to sidewall).All the simulations are initialized with the linear conduction profile for temperature (which is a function of the z-coordinate only) and a random noise of amplitude A = 0.001 along the z-direction for velocity.We run the simulations initially on a coarse grid of 120 × 120 × 30 points for 100 free-fall time units in which they converge to a statistically steady state.Following this, we successively refine the mesh to the required resolutions and allow the simulations to converge after each refinement.Once the simulations reach the statistically steady state at the highest resolution, they are run for another 20 to 21 free-fall time units and a snapshot of the flow field is saved after every free-fall time unit.
Since all the walls are no-slip, thin velocity boundary layers are formed adjacent to the walls.For our simulations to be well-resolved, an adequate number of gridpoints need to  1. Parameters of the simulations: the Hartmann number (Haz) based on the vertical magnetic field, the Rayleigh number (Ra), the ratio R of the magnetic field strength along the horizontal to the vertical directions, the grid-size, the number of points in the velocity boundary layers along the horizontal walls (nz), the y=±Γ/2 sidewalls (ny), and the x=±Γ/2 sidewalls (nx).The Prandtl number is constant in all cases at a value of Pr = 0.025.
be present in these boundary layers.It must be noted that the boundary layer profiles are strongly influenced by the magnetic fields.For a purely vertical magnetic field, these boundary layers are categorized into Hartmann layers adjacent to the top and bottom walls and Shercliff layers adjacent to the sidewalls.The thickness of Hartmann layers is given by δ H = 1/Ha z and that for the Shercliff layers is given by δ S = 1/ √ Ha z .However, in our case, the magnetic field is inclined with respect to the vertical direction; therefore, both the horizontal walls and y=±Γ/2 sidewalls will have a mix of Hartmann and Shercliff layers.On the other hand, the x=±Γ/2 sidewalls will have purely Shercliff layers.Thus, we use δ only for representing the boundary-layer thickness.For a conservative analysis, we estimate the thicknesses of these boundary layers as follows. (2.33) Table 1 lists the important parameters of our simulation runs.In this table, we also report the number of points in the different velocity boundary layers.We ensure that a minimum of 10 points is present in the boundary layers so as to adequately resolve our simulation runs.
In the next section, we will discuss the results obtained from our stability analysis and numerical simulations.

Results
In this section, we make a detailed analysis of the structure of the wall-modes and their impact on the heat and momentum transport.We will first present the linear stability analysis which is followed by a discussion of the results obtained from our direct numerical simulations.

Linear stability analysis
In this subsection, we discuss the results from our linear stability model.We assume that the set of finite positive eigenvalues of problem (2.26) is sorted in ascending order, i.e. the smallest one will be referred to as the first eigenvalue etc.
In figures 2(a-c), we plot the first eigenvalue, denoted as Rayleigh number Ra, at which the instability sets in due to disturbances at wavenumber β for different values of R. The minimum value of Ra is the critical Rayleigh number Ra c , and the wavenumber corresponding to Ra c is the most unstable wavenumber.We also plot Ra for the infinite plane layer for the corresponding Ha z , where the flow is assumed to be uniform in the y-direction.The figures show that Ra c is smaller than that corresponding to the plane infinite layer (Ra c,∞ ) for all R, Ha z , and aspect ratios considered in this study.This clearly implies that sidewalls destabilize the magnetoconvection system.Figure 3(a) exhibits the dependence of Ra c on R for different aspect ratios and strengths of the vertical magnetic field.It is evident from the above figure that for R < 1, the critical Rayleigh number increases for all aspect ratios.For the larger aspect ratio box (Γ = 4), Ra c continues to increase with R beyond R = 1 and saturates at R = 2. On the other hand, for the smaller aspect ratio box, the critical Rayleigh number starts to decrease with R for R > 1, a trend that does not seem to depend on Ha z .
In figure 3(b), we exhibit the plots of the most unstable wavelength λ versus R. The most unstable wavelength is calculated as where β c is the most unstable wavenumber.The figure shows that λ, like Ra c , exhibits a non-monotonic variation with R and lies between 1.1 and 1.6 for the entire range of parameters considered in our study.For small values of R, λ decreases with R and for the larger aspect-ratio box, it saturates at R = 2.For the smaller aspect-ratio box, λ increases with increasing R beyond R = 2.We now examine the behavior of the first three eigenvalues (Ra) for our cases.In figures 4(a)-(l), we plot the variations of these eigenvalues with β.The neutral stability curve for the infinite plane convection layer for the corresponding Ha z is also shown in dashed curves.For R = 0, it can be seen in figures 4(a,e,i) that the minima of the third eigenvalue (represented by green triangles) overlaps with that of the plane convection layer.This indicates that the third eigenvalue corresponds to the onset of bulk convection for different wavenumbers.The figures also indicate that minima of the third eigenvalue, which corresponds to the critical Rayleigh number for convection in the bulk, increases with R.This indicates that the bulk convection is further suppressed as the horizontal magnetic field increases.
The first and second eigenvalues (denoted by black squares and red circles respectively) correspond to wall-attached convection, and their minima are less than that for the third eigenvalue (corresponding to bulk convection).These eigenvalues nearly overlap for small horizontal magnetic fields but begin to diverge at large horizontal magnetic fields.These variations are shown more explicitly in figures 5(a,b) in which we exhibit the variations of the first two eigenvalues with R for β = 4 (near the most unstable wavenumber).It is evident from these figures that the eigenvalues diverge visibly for R > R c ≈ 1.5 for Γ = 2 and for R > R c ≈ 2.5 for Γ = 4.It is interesting to note that the value of R for which the eigenvalues begin to diverge does not seem to depend on Ha z .
We will now examine the spatial structure of the wall-modes, i.e. the eigenfunctions that correspond to the first and second eigenvalues.The eigenfunctions are normalized such that the maximum of the vertical velocity becomes equal to unity.
In figures 6(a)-(p), we exhibit the contour plots of the temperature perturbation θ on the vertical y-z plane corresponding to the first and second eigenvalues for β = 4, Γ = 2, and different Ha z and R. The first and third rows of the figure correspond to the first eigensolutions for Ha z = 50 and Ha z = 100, respectively, whereas the second and fourth rows correspond to the second eigensolutions.The figures show that as the horizontal magnetic field strength is increased, the wall modes get elongated and tilt along the direction of the resultant magnetic field.In fact, for R ⩾ 2, the wall modes are no longer confined adjacent to the sidewalls; instead, they occupy almost the entire bulk and interact with each other as explained later in this section.
Figures 6(a)-(p) also show that two there are two types of spatial structures corresponding to the first two eigensolutions.The first type consists of a hot plume (θ > 0) adjacent to the y = −Γ/2 sidewall and a cold plume (θ < 0) adjacent to y = Γ/2 sidewall; the corresponding eigensolution will be referred to as antisymmetric solution, see panel (a).The second type consists of cold (or hot) plumes adjacent to both the sidewalls; the corresponding eigensolution will be referred to as symmetric solution, see panel (e).For R < 1.5, there is only a marginal difference between the symmetric (Ra s ) and antisymmetric (Ra a ) eigenvalues.As exhibited in figure 7, this difference is less than 0.5% for R < 0.8 and is approximately 10% at R = 1.5 for the case of Ha z = 50 and Γ = 2. Thus, for R < 1.5, there is nearly an equal preference for symmetric and antisymmetric structures to develop at the onset of convection.However, these eigenvalues deviate significantly once R exceeds the threshold R c ≈ 1.5, above which the eigenvalue corresponding to the symmetric solution becomes significantly smaller.This implies that there is a stronger preference for the symmetric structures to develop at the onset of convection for R > R c .In this regime of R, the symmetric eigensolution comprises of a large plume developed by the merging of two cold (or hot) plumes adjacent to the opposite walls as visible in figures 6(c), (d), (k), and (l).On the other hand, as seen in figures 6(g), (h), (o), and (p), the antisymmetric eigensolutions for R > R c comprise of a cold plume extending on the top of a hot plume (or vice-versa) extended from the opposite wall, resulting in two convection rolls on top of each other.Our observations imply that a clear separation of the symmetric and antisymmetric solutions occurs when the wall modes adjacent to the opposite walls for symmetric solutions start interacting with each other.The merging of the plumes and the resultant formation of a merged roll results in an increased heat and momentum transport and hence in the decrease of the critical Rayleigh number.
Figures 8(a-h) exhibit the contour plots of θ on vertical y-z midplane corresponding to the first and second eigensolutions for Γ = 4, β = 4, and Ha z = 50.These figures again show that the wall modes tend to elongate along the direction of the resultant magnetic field.Similar to the Γ = 2 cases, the two eigenvalues correspond to symmetric and antisymmetric solutions respectively.The wall modes are symmetric for the first eigensolution and antisymmetric for the second.It can be recalled from figure 5(a) that the first and second eigenvalues for Γ = 4 box begin to diverge at a higher value of R compared to the Γ = 2 box.A clear separation occurs near R c ≈ 2.5; this is because owing to the larger aspect ratio of the box, the plumes adjacent to opposite walls are able to interact and merge only at higher tilts, and hence at larger R. The merged plume is exhibited in figure 8(g) which displays the contours of θ corresponding to the symmetric solution for R = 3.
Our analysis suggests that the nonmonotonic behaviour of Ra c and wavelength λ with respect to R is due to the increasing interaction of the plumes on opposite sidewalls.As long as the opposite plumes do not merge, an increase in the horizontal magnetic field component further stabilizes the magnetoconvection system with Ra c increasing and the wavelength λ decreasing with R. The system gets destabilized due to the merging of the opposite plumes and the variations of Ra c and λ with R get reversed.In the next subsection, we will discuss the results of direct numerical simulations of the nonlinear evolution of magnetoconvection.

Results of direct numerical simulations
In this subsection, we analyze our steady-state DNS results and examine the structures of the wall modes, their role in the heat and momentum transport, and the effects of initial conditions on the formation of the wall modes.

Structure of the wall modes
We use our numerical data to study the spatial convection structures for Ha z = 100 and Ha z = 200 with R ranging from 0 to 3.
In figures 9(a)-(j), we exhibit the isosurfaces of u z = ±0.01 for all our runs.In figures 10(a)-(h), we display the contours of u z ⩾ 0.01 (red) and u z ⩽ 0.01 (blue); these contours give a visualization of the upwelling and downwelling plumes respectively.These figures show the presence of wall-attached convection with suppressed fluid flow in the bulk.Consistent with the results of the stability analysis in the previous subsection, these wall modes become elongated and align themselves along the direction of the resultant magnetic field as the horizontal magnetic field is increased.Figures 9 and 10 also show that with the exception of the case with Ha z = 100 and R = 3, the wall modes are antisymmetric, that is, the plumes are upwelling (or downwelling) adjacent to y = −Γ/2 wall and downwelling (or upwelling) adjacent to y = Γ/2 wall.The wall modes occupy the entire bulk for R = 3, consistent with the stability analysis of the Γ = 4 box.For the above field ratio, the wall modes for the Ha z = 100 case are symmetric and the plumes adjacent to the opposite walls merge to form a large plume as displayed in figures 9(e) and 10(d).
It can also be observed in figures 9(a-j) that there are slight irregularities in the wall-mode structures.These irregularities seem to be associated with small temporal fluctuations (∼ 1%) in integral quantities that remain after our simulations reach an apparently stationary state.These fluctuations indicate that the wall modes are still evolving, albeit slowly.However, it is important to note that we could not observe any noticeable change in the structures of the wall modes over the time-span of our simulations.Any change in the spatial arrangement of the modes is likely to be visible only after the simulations are run for several times the diffusion time scale, which is around 150 to 200 free-fall time units.Hence our solutions can be considered to be effectively steady.
We estimate the wavelength of the wall modes by visual inspection of the vertical velocity isosurfaces in figures 9(a)-(j).We consider only those wall modes that are adjacent to y = ±Γ/2 sidewalls; these walls are not parallel to the resultant magnetic field.The estimation of the wavelength is done as follows.We count the number of upwelling (or downwelling) plumes adjacent to both y = −Γ/2 and y = Γ/2 sidewalls and calculate their average as N p .The wavelength λ wm of the wall modes is given by We plot the estimated λ wm versus R in figure 11.The figure shows that λ wm tends to decrease as the horizontal magnetic field increases and saturates at large R. Further, the wavelength corresponding to Ha z = 200 is smaller than that for Ha z = 100 except for R = 3.This is in line with the fact that the threshold wavelength decreases with increasing vertical magnetic field strength (Chandrasekhar 1981).It can be recalled that a similar trend was observed in the variation of the most unstable wavelength with R in our stability analysis in § 3.1.However, it must be kept in mind that since the Rayleigh numbers considered in our DNS are higher than the threshold Rayleigh number above which the wall modes appear, there is always a chance for secondary modes with smaller wavenumbers to develop.

Heat and momentum transport
We will now explore the influence of the wall modes on the heat and momentum transport.We compute the Nusselt (Nu) and the Reynolds numbers (Re) using our 0.0 0.5 1.0 1.5 2.0 2.5 3.0 numerical data as follows: where U rms = ⟨u 2 x + u 2 y + u 2 z ⟩ V is the root mean square velocity and ⟨•⟩ V denotes volume averaging.The second term of the right hand side of (3.2) is the normalized convective heat flux.We plot the normalized heat flux and the Reynolds number versus R for all our runs in figures 12(a) and (b) respectively.The figures show that for R < 1, both Nu and Re decrease with R.This is consistent with our conclusion from the stability analysis in § 3.1 that for small values of R, an increase in the horizontal magnetic field results in the stabilization of the magnetoconvection system.There is, however, no clear trend in the variations of Nu and Re for R > 1.
Having studied the trends of the global heat and momentum transport with R, we will now explore the spatial variation of heat and momentum transport for different regimes of R. In figures 13(a)-(e), we plot twice the kinetic energy E = 0.5(u 2 x +u 2 y +u 2 z ), averaged over the x-z plane, versus y, which is the direction parallel to the horizontal magnetic field.For R = 0, there are two sharp peaks close to y = −2 and y = 2 for both Ha z = 100 and Ha z = 200; these peaks correspond to the wall modes.The bulk region lies in between Figure 13.distribution of kinetic energy E and the convective heat flux uzT .Plots of 2⟨E⟩x,z = ⟨u , and (e) R = 3. Plots of ⟨uzT ⟩x,z for (f) R = 0, (g) R = 0.3, (h) R = 1, (i) R = 2, and (j) R = 3.In the above, ⟨•⟩x,z represents averaging over x-z plane.
these two peaks.The bulk consists of several smaller peaks and the kinetic energy in this region is much less compared to the near-wall regions.As R is increased, the height of the near-wall peaks decreases and of those in the bulk increases.Thus, the distribution of kinetic energy along y becomes more uniform as R is increased.It is a known fact that in magnetohydrodynamic flows, the Lorentz force tends to modify the flow-field so as to minimize the gradients of velocity along the direction of the magnetic field (Davidson 2017).Thus, in our case, as R is increased, the Lorentz force generated due to the ycomponent of the magnetic field becomes strong and suppresses the gradients of velocity along y.
In figures 13(f)-(j), we plot the local convective heat flux u z T , averaged over the x-z plane, along y.Although the volume-averaged heat flux is positive in thermal convection, the local heat flux for R = 0 fluctuates between positive and negative values as one proceeds along y.These fluctuations even out as R is increased, and for R = 3, the local heat flux remains positive throughout with very small gradients along y.Again, this is due to the strong Lorentz forces generated by the horizontal component of the magnetic field which suppresses the gradients of velocity along the magnetic field's direction.
In figures 14(a)-(e), we plot twice the kinetic energy, averaged over the y-z plane, versus x, which is the direction perpendicular to the horizontal magnetic field.In the absence of horizontal magnetic field (R = 0), the variation of 2⟨E⟩ y,z versus x is similar to that of 2⟨E⟩ x,z versus y due to the symmetry of the problem.Again, there are two sharp peaks corresponding to the wall modes close to x = −2 and x = 2 for both Ha z = 100 and Ha z = 200.However, unlike ⟨E⟩ x,z , the values of ⟨E⟩ y,z at the peaks recede only marginally as R is increased from 0 to 2. As the wall modes extend fully into the bulk at R = 3, the above peaks recede sharply with the value of 2⟨E⟩ y,z at these peaks being close to its values at the peaks in the bulk.It must be noted that ⟨E⟩ y,z continues to fluctuate as it is varied with x at R = 3 and does not smoothen out unlike ⟨E⟩ x,z .
Figures 14(f-j) exhibit the variations of the local convective heat flux u z T , averaged over the y-z plane, along x.The figures show that ⟨u z T ⟩ y,z fluctuates between positive and negative values.It is clear from the figures that for the Ha z = 100 case, the spatial fluctuations increase with increasing R.This is because as discussed before, the gradients along the y-direction are reduced as R is increased.Therefore, for large R, if the gradients along z are also small, any quantity averaged over y-z plane will fluctuate between the quantity's two extremities.The case of R = 3 for Ha z = 100 consists of plumes from the opposite y=±Γ walls merged with each other; therefore, the gradients of the velocity, temperature, and heat flux along the z-direction are small.Thus, ⟨u z T ⟩ y,z exhibits strong fluctuations along the x-direction at R = 3.
For Ha z = 200 and R ⩽ 2, ⟨u z T ⟩ y,z follows a similar trend in that its fluctuations increase with increasing R.However, for R = 3, its fluctuations get suppressed and ⟨u z T ⟩ y,z is mostly positive.Now, let us recall that unlike in the case of Ha z = 100, the structures for the case of Ha z = 200 are antisymmetric.Thus, although the gradients along the y-direction are small, there are fluctuations along the z-direction for Ha z = 200 due to the presence of upwelling and downwelling plumes on top of each other.An averaging over the y-z plane cancels out the opposing effects of the upwelling and downwelling plumes, thus resulting in the suppression of fluctuations of ⟨u z T ⟩ y,z .
Finally, we analyze the contributions by the bulk and near-wall regions to the total kinetic energy and heat flux of the system.Towards this objective, we define the near-wall region as follows.For R = 0 and for both the Hartmann numbers, we determine Nu S , which is the Nusselt number averaged over successively smaller concentric volumes S = Γ ×[r y , Γ −r y ]×1.In this definition, r y is the normal distance from the y=±Γ/2 sidewalls.We plot Nu S versus r y in figure 15.The figure shows that Nu S initially decreases with r y upto a point of local minima and then begins to increase with r y .The distance between the point of the local minima and the sidewall is taken as the width δ w of the near-wall region.These widths are computed to be δ w = 0.37 for Ha z = 100 and δ w = 0.33 for Ha z = 200.Although the above values of δ w were computed for R = 0, these values will be assumed to hold for all R.
The total kinetic energy E and the convective heat flux H can be expressed as the sum of their bulk and near-wall contributions.Therefore, where E bulk and E nw are respectively the bulk and near-wall contributions to the total kinetic energy, and H bulk and H nw are respectively the bulk and near-wall contributions to the total heat flux.The total and bulk contributions to these quantities are defined as We compute the relative strengths of the bulk and near-wall kinetic energies and heat fluxes using our numerical data and plot them versus R in figures 15(a)-(d).For small values of R, the bulk contribution to the total kinetic energy and heat flux is very small (less than 10%).This is expected because convection is completely suppressed in the bulk at small R. It is clear from the figures that as R increases, the bulk contributions to the total kinetic energy and heat flux increase.This is due to the fact that the wall modes extend further into bulk as R increases.In fact, for R > 2, the bulk and sidewall contributions become comparable to each other because the wall modes at such high values of R extend fully into the bulk.It is also interesting to note that for Ha z = 100, the bulk contribution to the heat flux decreases as R is increased from 2 to 3. The reason for this anomalous behaviour is yet to be understood.

Effects of initial conditions
The direct numerical simulations of magnetoconvection for R = 3 resulted in different structural arrangements of the convection plumes for Ha z = 100 and Ha z = 200.On one hand, we obtained a symmetric arrangement of the plumes adjacent to the opposite walls merging with each other for Ha z = 100.On the other hand, an antisymmetric arrange- ment of the plumes was obtained for Ha z = 200 with the upwelling and downwelling plumes on top of each other.In this subsection, we will explore the sensitivity of the above results to initial conditions.We conduct two more direct numerical simulations of magnetoconvection for R = 3; one with Ha z = 100, Ra = 10 5 , and the other with Ha z = 200, Ra = 4 × 10 5 .However, we change the initial conditions in this set of simulations as follows.The results of the old simulation of Ha z = 200 are taken as the initial condition for the new simulation of Ha z = 100.In the same way, the results of the old simulation of Ha z = 100 are taken as the initial condition for the new simulation of Ha z = 200.Both these simulations were allowed to run for 20 free-fall time units after reaching the steady state.Henceforth, we refer to the old set of simulations as IC1 and the new set of simulations of R = 3 as IC2.
We examine the structure of the wall-modes for the new set of simulations by plotting the vertical velocity isosurfaces for u z = ±0.01 in figures 17(a) and (c).The figures show that the simulation of Ha z = 100 with the new initial conditions yields antisymmetric arrangement structures with upwelling and downwelling plumes on top of each other.This is unlike the result of the simulation with old initial conditions which resulted in a symmetric arrangement of structures with merged plumes.Similarly, the simulation of Ha z = 200 with the new initial conditions yield symmetric structures with merged plumes, unlike the case with old initial conditions which resulted in antisymmetric arrangement of structures.Our results indicate that that the solution of magnetoconvection with R = 3 for Ha z = 100 and 200 is non-unique, and the arrangement of the wall modes depends on the initial conditions.
We will now explore the impact of the arrangement of the plumes on the global heat and momentum transport.In figures 18(a,b), we replot the computed values of Nu − 1 and Re for our runs with the old initial conditions and also add the plots of these quantities computed using the results of our simulations with new initial conditions.The figures show that for Ha z = 200, the Reynolds and Nusselt numbers computed using the solutions of simulations with new initial conditions are significantly higher than those corresponding to the old initial conditions.On the other hand, for Ha z = 100, the Reynolds and Nusselt numbers computed using the solutions of simulations with new initial conditions are lower, albeit marginally, than those corresponding to the old initial conditions.The above observations clearly indicate that the solutions that are symmetric consisting of merged plumes have higher heat and momentum transport.This is because merged plumes give rise to larger convection rolls which, in turn, result in more efficient transport of heat and momentum.It is worth noting that for Ha z = 200, the increase of Reynolds number at R = 3 due to the merged plumes is so significant that its value is larger than for R = 0.
Our results indicate that the non-unique nature of the solutions for R = 3 has a profound effect on the heat and momentum transport, especially for Ha z = 200.We conclude in the next section.

Summary and conclusions
In this paper, we systematically examined the effects of strong inclined magnetic fields on wall-attached magnetoconvection using a combination of linear stability analysis and direct numerical simulations.The linear stability analysis was conducted for lower Hartmann numbers and aspect-ratio box whereas the DNS were conducted for higher Hartmann numbers and aspect ratio.The ratio R = B y /B z of the outer imposed horizontal to vertical magnetic field was varied from 0 to 3.
The linear stability analysis revealed that the critical Rayleigh number varies nonmonotonically with the relative strength of the horizontal magnetic field.The plumes at the onset of instability get elongated along the direction the resultant magnetic field and extend more into the bulk as R is increased.At sufficiently high R, the plumes extend fully into the bulk.The linear stability analysis further revealed the existence of nonunique solutions at the onset at low R.One eigensolution corresponds to a symmetric arrangement of hot or cold plumes adjacent to the opposite sidewalls, whereas the other eigensolution corresponds to an antisymmetric arrangement of hot or cold plumes.For the symmetric solution, when the opposite plumes extend sufficiently into the bulk at a large R, they merge into a single large plume.Below this critical value of R = R c , the symmetric and antisymmetric eigensolutions overlap; however, at R > R c , these solutions start to diverge with the symmetric eigensolution being more unstable than the antisymmetric solution.The critical Rayleigh number increases with R for R < R c and decreases with R for R > R c .
The direct numerical simulations of the fully nonlinear regime reinforced the observation from the linear stability analysis that the wall modes get elongated along the direction of the resultant magnetic field and extend further into the bulk as R is increased.The heat and momentum transport was observed to decrease with R for R < 1 but not show any visible trend for R > 1.The fluctuations of the local heat flux and kinetic energy along the direction of the horizontal magnetic field get suppressed as R is increased.Since the wall modes extend more into the bulk as R is increased, the relative contribution to the total heat transport and kinetic energy by the near-wall regions decrease with an increase of R.
The analysis from direct numerical simulations further revealed that at least for R = 3, the solutions are dependent on the initial conditions.The solutions for R = 3 can either comprise of an antisymmetric arrangement of upwelling and downwelling plumes on top of each other, or a symmetric arrangement of merged upwelling or downwelling plumes.The solution with merged plumes corresponds to higher heat and momentum transport because the merged plumes give rise to larger convection rolls and hence more efficient heat transfer.
Our present work provides important insights into the dynamics of wall-attached convection under inclined magnetic fields which will be a typical situation in view to most applications.Our work may find applications in several industrial flows such as cooling blankets in fusion reactors.Although we worked on a small set of parameters, we expect our results to hold for higher Rayleigh numbers as well.In the future, we plan to conduct a similar analysis for fluids at different Prandtl numbers.

Figure 1 .
Figure 1.A sketch of the Rayleigh-Bénard convection setup with inclined magnetic field employed for our linear stability analysis.

Figure 11 .Figure 12 .
Figure 11.Variations of the wavelength λwm of the wall modes with R. The wavelength tends to decrease with an increase in the horizontal magnetic field.(a) (b)

Figure 15 .
Figure 15.Determination of the near-wall region for our analysis from the results of our DNS for R = 0. Variation of the Nusselt number NuS averaged over successively smaller concentric volumes with the sidewall-normal distance ry.The width of the near-wall region is given by the point of the first minima (represented as black dashed vertical line for Haz = 100 and blue dashed vertical line for Haz = 200) in the NuS profile

Figure 16 .
Figure 16.Results of direct numerical simulations: Relative strengths of the bulk and boundary layer contributions to the total kinetic energy E for (a) Haz = 100 and (b) Haz = 200, and to the total heat flux F for (c) Haz = 100 and (d) Haz = 200.The bulk contributions to kinetic energy and heat flux increase with increasing R.

Figure 17 .Figure 18 .
Figure 17.Results of direct numerical simulations with new sets of initial conditions (IC2): Isosurfaces of uz = 0.01 (red) and uz = −0.01(blue) for R = 3 and (a) Haz = 100 and (c) Haz = 200.For comparison, the vertical velocity isosurfaces using results of the simulations of R = 3 with the old initial conditions are shown for (b) Haz = 100 and (d) Haz = 200