Bifurcation analysis of a two-dimensional magnetic Rayleigh–B´enard problem

We perform a bifurcation analysis of a two-dimensional magnetic Rayleigh– B´enard problem using a numerical technique called deflated continuation. Our aim is to study the influence of the magnetic field on the bifurcation diagram as the Chandrasekhar number Q increases and compare it to the standard (non-magnetic) Rayleigh–B´enard problem. We compute steady states at a high Chandrasekhar number of Q = 10 3 over a range of the Rayleigh number 0 ≤ Ra ≤ 10 5 . These solutions are obtained by combining deflation with a continuation of steady states at low Chandrasekhar number, which allows us to explore the influence of the strength of the magnetic field as Q increases from low coupling, where the magnetic effect is almost negligible, to strong coupling at Q = 10 3 . We discover a large profusion of states with rich dynamics and observe a complex bifurcation structure with several pitchfork, Hopf, and saddle-node bifurcations. Our numerical simulations show that the onset of bifurcations in the problem is delayed when Q increases, while solutions with fluid velocity patterns aligning with the background vertical magnetic field are privileged. Additionally, we report a branch of states that stabilizes at high magnetic coupling, suggesting that one may take advantage of the magnetic field to discriminate solutions.

In this paper, we focus on a two-dimensional (2D) magnetic Rayleigh-Bénard problem, in which an external magnetic field is applied to an electrically conducting fluid, and study how the strength of the magnetic field influences the bifurcation patterns.These so-called magnetohydrodynamics (MHD) convection problems lead to interesting instability [25,26], heat transport [27,28], flow reversal [29] phenomena, and have many applications in astrosciences and geosciences [30][31][32][33].Then, the bifurcation analysis of the magnetic Rayleigh-Bénard problem using either numerical simulations or laboratory experiments is an active field of research [34][35][36][37].As an example, Akhmedagaev et al. performed numerical simulations of turbulent Rayleigh-Bénard convection in a cylinder under a vertical magnetic field at high Rayleigh and Hartmann numbers [37].More recently, Yang et al. performed a laboratory experiment to study the hydromagnetic convection of a liquid metal in a cuboid vessel and analyze the transition from steady to oscillating convection rolls in the regime Ra ∈ [2.3 × 10 4 , 2.6 × 10 5 ] [34].Finally, Naffouti et al. investigated the effect of the direction of the magnetic field on the convection flow patterns [36].
Here, we investigate numerically the effect of the magnetic field on the bifurcation diagram for the Rayleigh number Ra ∈ [0, 10 5 ] as the Chandrasekhar number increases from 1 to 10 3 with Prandtl and magnetic Prandtl numbers fixed to 1.By combining deflation [23] and a continuation technique in the Chandrasekhar number, we are able to relate primary instabilities and secondary branches at Q = 10 3 to the corresponding Rayleigh-Bénard branches in the low magnetic field limit near Q = 0.This allows us to perform qualitative and quantitative comparisons between the convection patterns discovered by deflation at Q = 10 3 and the related branches analyzed in [24], as well as study the effect of the magnetic field strength on the stability of the steady states.
The paper is organized as follows.We first introduce the non-dimensional formulation of the Rayleigh-Bénard problem and describe the numerical methods in Section II.Then, we present the numerical results and analyze the effect of the strength of the magnetic field on the bifurcation diagram in Section III.Finally, in Section IV, we summarize the conclusions of this work and propose potential extensions.

II. MAGNETIC RAYLEIGH-B ÉNARD PROBLEM IN 2D
A. Formulation and discretisation We consider a two-dimensional magnetic Rayleigh-Bénard problem modeling a conducting fluid heated from below in a unit square cell domain Ω = (0, 1) 2 , with spatial coordinates (x, z).We assume no-slip boundary conditions for the fluid velocity and a background vertical magnetic field pointing upwards.Furthermore, we choose the horizontal walls to be thermally conducting and the vertical walls to be insulating.
We employ the Boussinesq approximation [47,48] to model magnetohydrodynamics convection and assume that the flow is buoyancy-driven and that density differences only appear in the buoyancy term, while other parameters do not depend on density or temperature.Then, one can write down the non-dimensional formulation of the two-dimensional anisothermal MHD equations with Boussinesq approximation on the domain Ω as where u = (u, w) denotes the fluid velocity field, p the fluid pressure, B = (B x , B z ) the magnetic field, E the electric field, T the temperature, ẑ = 0 1 the buoyancy direction, and ε(u) = 1 2 (∇u + ∇u ) the symmetric gradient.For clarity, we denote vector fields with bold letters in Eq. (1).
The dimensionless numbers in Eq. ( 1), i.e. the Prandtl, magnetic Prandtl, Rayleigh, and Chandrasekhar numbers, are respectively defined as Pr Here, ρ 0 denotes the reference density, ν the kinematic viscosity, β is the thermal expansion, η the magnetic resistivity, g the gravitational acceleration, α the thermal diffusivity, µ 0 the magnetic permeability of free space, B a characteristic value for the magnetic field, L a characteristic length scale, and T the difference between two reference temperatures (e.g., the reference temperatures may be defined as the temperature of the hot and cold plates at the top and bottom of our domain).In the numerical experiments performed in Section III, we will study the influence of the magnetic field on the bifurcation diagram with respect to the Rayleigh number by fixing Pr = 1, Pm = 1, and varying Ra and Q in the range Ra ∈ [0, 10 5 ] and Q ∈ [1, 10 3 ].We emphasize that in two dimensions the electric field is a scalar field and is hence denoted as E.Moreover, there exist two curl operators and cross products which are interpreted differently depending on whether the arguments are scalar-or vector-valued as For this reason, we employ the boundary condition E = 0 on ∂Ω, which is the two-dimensional version of the standard condition E × n = 0 on ∂Ω.We then complement Eq. ( 1) with the boundary conditions u = 0, B • n = ẑ, and E = 0 on the boundary of the domain ∂Ω, where n denotes the unit outer normal vector.In addition, we enforce T = 1 at z = 0 and T = 0 at z = 1, as well as ∂ x T = 0 at x = 0 and x = 1.The trivial stationary solution (conducting state) to Eq. (1) with our choice of boundary conditions is given by φ 0 = (u 0 , p 0 , T 0 , B 0 , E 0 ), where Additionally, we remark that the problem has two symmetries which determine the behavior of the arising bifurcations with respect to Ra and Q: the mirror symmetry with respect to and the Boussinesq symmetry, These symmetries are similar to the Rayleigh-Bénard problem (cf.[24, Sec.2]) but also transform the electric and magnetic fields.We reformulate the stationary version of Eq. ( 1) by combining Eqs.(1d) and (1f) to obtain the following augmented Lagrangian formulation: Hu et al. showed that this formulation enforces both ∇ • B = 0 and curl E = 0 and makes the stationary version of Eq. ( 1) well-posed [49].This is also the reason why the electric field E is kept in our formulation and not eliminated via (1c).Moreover, the presented formulation allows for a structure-preserving finite-element discretization that enforces both ∇ • B = 0 and curl E = 0 exactly at the discrete level [49].Indeed, preserving the discrete solenoidality of B is crucial for numerical simulations to produce physically accurate results [50].
We employ a finite element method to solve the stationary magnetic Rayleigh-Bénard problem in two dimensions.The velocity and pressure are discretized using the standard Taylor-Hood elements [51], i.e., piecewise vector valued quadratic polynomials for u and piecewise linear polynomials for p.Then, the temperature and electric fields are respectively approximated by piecewise linear and quadratic polynomials.Finally, we use second order Raviart-Thomas elements [52] for the discretization of the magnetic field, which is a standard choice for a H(div)-conforming finite element that allows for the aforementioned structure-preserving discretization.The domain Ω is discretized using a crossed triangular mesh with 50×50 square cells, where each square cell is divided into four triangles following the diagonals of the squares to preserve the symmetries of the problem.While this is a rather coarse mesh, we emphasize that computing the full bifurcation diagram might require to perform hundreds of thousands of Newton iterations, which is computationally expensive.However, we checked that refining the mesh does not influence the profile of the steady states we are computing in the range of parameters, 0 ≤ Ra ≤ 10 5 and 1 ≤ Q ≤ 10 3 , considered.Moreover, this mesh allows us to employ a sparse LU direct solver to solve the underlying linear systems.We produce our numerical results using the Firedrake finite element software [53], which relies on the PETSc package [54] for solving the discretized equations.

B. Bifurcation technique
Steady states to the magnetic Rayleigh-Bénard system (1) are computed using a numerical technique called deflation [23].This algorithm aims to discover multiple solutions to nonlinear partial differential equations (PDEs) in an iterative manner by modifying Newton's method in order to avoid its convergence to previously obtained solutions.While deflation has been successfully applied to a wide range of physical systems such as liquid crystals [55,56], quantum mechanics [57][58][59][60], and fluid dynamics [24], to our knowledge this work is its first adaptation to MHD-type problems.To introduce this numerical method, we first rewrite our bifurcation problem as F (Φ, λ) = 0, where Φ = (u, p, T, B, E) is the solution, λ ∈ R denotes the bifurcation parameter, and F is the partial differential operator associated with the stationary version of the MHD system (1).In this work, we employ either the Rayleigh number or the Chandrasekhar number as bifurcation parameters and keep the other dimensionless numbers fixed.For a given λ ∈ R, assuming that Newton's method has converged to a solution Φ 1 satisfying F (Φ 1 , λ) = 0, then deflation constructs a new (deflated) problem of the form The deflation operator M is chosen to penalize solutions that are close to Φ 1 such that lim φ→φ1 M(Φ, Φ 1 ) = ∞.
In this work, we employ the following deflation operator: where the operator N is defined as such that N (φ, φ 1 ) converges to zero as φ approaches φ 1 .Here, • 2 denotes the L 2 -norm and Then, one can apply Newton's method on the deflated problem F 1 from the same initial guess used to discover φ 1 without converging to the same solution, and hopefully discover a different solution to F (φ, λ) = 0.A key benefit of deflation is that solving Eq. ( 3) via Newton's method does not require more computational work than solving the original system since once can express the Newton update of the discrete deflated system as a scaling of the update for the original problem [23].This process can be repeated to discover several steady states to the magnetic Rayleigh-Bénard problem by iteratively composing deflation operators to deflate multiple roots.
Once deflation has been applied at a fixed value of the bifurcation parameter λ, we interlace it with a continuation procedure to reconstruct a bifurcation diagram over a range of parameters λ ∈ [λ min , λ max ].We then use the discovered solutions at λ as initial guesses for a parameter continuation from λ to λ ± ∆λ, where ∆λ is the step-size and the sign determines the sense (forward or backward) of the continuation.The algorithm then iteratively continues by applying a deflation to each discovered solution at λ ± ∆λ until the final values of λ are reached.
Since one initially employs Newton's method on the original problem, the performance of deflation at finding new solutions depends heavily on the available initial guesses for the first deflation step.To address this issue, one technique applied for the Rayleigh-Bénard problem [24] is to construct initial guesses for backward continuation starting from Ra = 10 5 by summing the trivial solution and its normalized unstable eigenmodes.While the same approach works for the magnetic Rayleigh-Bénard problem at the small Chandrasekhar number Q = 1, it fails to provide useful initial guesses for Newton's method when applied in a regime of stronger magnetic fields, such as Q = 10 3 .Therefore, our remedy consists of first performing deflation with respect to the Rayleigh number in the range 0 ≤ Ra ≤ 10 5 at Q = 1, Ra Ra FIG.
1. Bifurcation diagrams of the magnetic Rayleigh-Bénard problem using the kinetic energy u 2 2 (top row), potential energy T 2 2 (middle row), and magnetic energy B 2 2 (bottom row) as diagnostic.The diagrams display the evolution of the branches studied in this paper, i.e., the first four primary instabilities at Q = 10 3 and their secondary bifurcations, in the range 0 ≤ Ra ≤ 10 5 at Q = 1 (a), 1 ≤ Q ≤ 10 3 at Ra = 10 5 (b), and 0 ≤ Ra ≤ 10 5 at Q = 10 3 (c).The first four primary branches for the bifurcation diagrams in (c) are displayed with solid lines while the secondary branches are indicated with dashed lines.The dark line in the diagrams using the potential energy corresponds to the conducting state satisfying T0 2 2 = 1/3.The x-axes of the diagrams in panel (c) are reversed to analyze the evolution of the branches as the strength of the magnetic field increases.
then prolonging this bifurcation diagram at Ra = 10 5 over 1 ≤ Q ≤ 10 3 with deflation.This procedure finally provides good initial guesses at Ra = 10 5 and Q = 10 3 for backward deflated continuation over Ra ∈ [0, 10 5 ].This results in three bifurcation diagrams: over Ra at Q = 1 (∆Ra = 10 3 /3), over Q at Ra = 10 5 (∆Q = 10/3), over Ra at Q = 10 3 (∆Ra = 500), which are respectively displayed in Fig. 1(a-c), and enables the study of the influence of the magnetic field on the bifurcation diagrams by relating the states at Chandrasekhar number Q = 10 3 to the ones at Q = 1.In this work, we are also interested in the evolution of the velocity, temperature, and magnetic profiles of the solutions and therefore used the squared L 2 -norm of these fields (resp.kinetic energy u 2  2 , potential energy T 2 2 , and magnetic energy B 2 2 ) as diagnostic for the bifurcation diagrams.Finally, we note that while a large profusion of states have been discovered by deflation, we focus on a selection arising from primary bifurcations of the trivial state and their second bifurcations to keep their analysis tractable.

C. Stability analysis
Once a stationary solution φ 1 = (u 1 , p 1 , T 1 , B 1 , E 1 ) to Eq. ( 1) has been found by deflation, we analyze its stability using the perturbation ansatz where 1 is a small parameter and λ ∈ C is an eigenvalue.Upon inserting Eq. ( 4) into Eq.( 1) to linearize it around the solution (u 1 , p 1 , T 1 , B 1 , E 1 ) and inserting the perturbation ansatz, we obtain the following generalized eigenvalue problem: , where I denotes the identity and the operators F and G are defined as We solve this generalized eigenvalue problem to uncover the eigenvalues with largest real parts using a Krylov Schur solver [61] that is implemented in the SLEPc library [62].
The real and imaginary parts of the computed eigenvalues determine the stability of the corresponding eigenmode.If the real part (growth rate) of all eigenvalues is negative, the solution is stable, whereas if at least one eigenvalue has a positive real part the solution is unstable.Moreover, the type of instability depends on whether the associated imaginary part is zero or non-zero.The change of stability of a solution branch at a bifurcation parameter λ indicates the presence of bifurcations, e.g., a simple bifurcation if the solution has a zero eigenvalue or a Hopf bifurcation if it has a pair of complex conjugate purely imaginary eigenvalues [63].
In this work, we are mostly focusing on primary bifurcations, i.e., bifurcations arising from the trivial state (u 0 , p 0 , T 0 , B 0 , E 0 ) defined by Eq. ( 2).These bifurcations occur at critical values, Ra c , of the bifurcation parameter Ra, which can be computed by re-organizing the eigenvalue problem associated with the trivial state at eigenvalue λ = 0 as a a generalized eigenvalue problem: , where we interpret Ra c as an eigenvalue.Recall that u 0 = E 0 = 0 and hence we dropped these terms above.
We then solve this problem to obtain the first ten critical Rayleigh numbers at Q = 10 3 and associated eigenmodes (see Fig. 2(e)).The first critical Chandrasekhar numbers, Q c , can also be computed at Ra = 10 5 by solving an analogous generalized eigenvalue problem: .

III. RESULTS
In this section, we report the bifurcation diagrams and solutions to the magnetic Rayleigh-Bénard problem discovered by deflation, and analyze their behaviour as the Rayleigh and Chandrasekhar numbers vary.

A. Location of primary instabilities
We begin our study by computing the stability of the conducting state (2) to localize the primary instabilities and resulting bifurcations.To do so, we solve the associated eigenvalue problem over a range of Rayleigh and Chandrasekhar numbers, namely 0 ≤ Ra ≤ 10 5 at Q = 1, 0 ≤ Q ≤ 10 3 at Ra = 10 5 , and 0 ≤ Ra ≤ 10 5 at Q = 10 3 using the numerical method described in Section II C. The positive real parts of the eigenvalues for the different intervals are reported in panels (a-c) of Fig. 2, respectively.
We observe in Fig. 2(a) that eleven supercritical pitch- fork bifurcations emanate from the conducting state at Q = 1 when the Rayleigh number increases to Ra = 10 5 .This plot is nearly identical to the one computed for the standard Rayleigh Bénard problem at Q = 0 [24, Fig. 1] indicating that effect of the magnetic field on the bifurcation structures is negligible when Q = 1.In our computation, the eleventh supercritical bifurcation takes place at Ra (11)   c ≈ 9.95×10 4 , while it starts slightly above Ra = 10 5 in the standard Rayleigh-Bénard problem and is therefore not included in [24].FIG. 3. Evolution of the steady states in the first primary branch (1) with respect to the Rayleigh number, where the Chandrasekhar is fixed to Q = 10 3 , illustrated via the kinetic energy (a), potential energy (b), and magnetic energy (c).The plots in these panel display the magnitude of the velocity, temperature, and magnitude of the magnetic field.The color schemes are defined as follow.For the velocity and magnetic fields, blue corresponds to zero while red stands for high magnitude (the color maps are rescaled to data range).For the temperature, blue corresponds to T = 0 and red to T = 1.The corresponding largest growth rates are reported in (g).(d-f) Same as (a-c) but the evolution is for the steady states in the secondary branch (6), with the largest growth rates plotted in (h).The green dots and red triangles indicate the presence of pitchfork bifurcations (eigenvalue λ = 0) and Hopf bifurcations (pair of complex conjugate purely imaginary eigenvalues).
the third one is a (3, 1)-eigenmode since the velocity field has three rolls in the horizontal direction and one in the vertical direction.
We now analyze the effect of the strength of the magnetic field on the primary bifurcations by fixing Ra = 10 5 and increasing the Chandrasekhar number to Q = 10 3 in Fig. 2(b).We find that the number of unstable eigenmodes decreases to five in this range, leading to the observation that, contrary to the Rayleigh number, increasing the Chandrasekhar number penalizes the birth of primary bifurcations as the growth rates of the conducting state decrease in Fig. 2(b).The delayed onset of instabilities when increasing Q aligns with the findings of Basak et al., while observed in a different parameter regime [42].
Fig. 2(c) displays the positive growth rates of the conducting state for 0 ≤ Ra ≤ 10 5 at Q = 10 3 .The first five primary bifurcations occur at the critical Rayleigh numbers at which there is a zero eigenvalue, i.e., Ra (1)  c ≈ 2.0 × 10 4 , Ra (2) c ≈ 5.9 × 10 4 , and Ra (5)   c ≈ 9.9 × 10 4 .The values of the critical Rayleigh numbers agree with the ones obtained by solving a generalized eigenvalue problem (cf.Section II C), and the locations of the bifurcations from the conducting state in Fig. 1(c).
Then, the velocity magnitude of the first five eigenmode associated with the critical Rayleigh numbers over 0 ≤ Ra ≤ 10 5 at Q = 10 3 is displayed in Fig. 2(e).It is interesting to observe that the velocity profile (in particular the vortices) of the eigenmodes feature a pattern that is oriented in the (upward) direction of the background magnetic field B 0 = ẑ.This indicates that for larger Chandrasekhar number at Q = 10 3 the instabilities aligned with the magnetic field occur at smaller Rayleigh numbers.In general, the order of eigenmodes changes with growing Q.For instance, the (4, 1)-eigenmode is the seventh eigenmode in Fig. 2(a) at Q = 1 but corresponds to the fourth eigenmode at Q = 10 3 .Similarly, we find that the (1, 1) and (2, 1)-eigenmodes are swapped  (d-f) Same as (a-c) but the evolution is displayed over Ra at Q = 10 3 .Panel (g) shows the largest growth rates (real parts of the eigenvalues) corresponding to the states in (a-c), while (h) reports the largest growth rates over Ra at Q = 10 3 , i.e., states in (d-f).
We proceed by displaying the bifurcation diagrams of the magnetic Rayleigh-Bénard problem in Fig. 1 for the kinetic, potential, and magnetic energy.In the bifurcation diagrams, we do not report symmetries of solutions generated by the mirror and Boussinesq symmetries.Additionally, since the effect of the magnetic field at Q = 1 is minor, we just track the evolution of the branches connected to the branches at Q = 10 3 in Fig. 1(a).We then refer to [24] for a detailed bifurcation analysis of the nonmagnetic case.

B. First and second primary branches
We begin our study of the solutions computed via deflation by analyzing the evolution of the steady states in branch (1) at Rayleigh number Ra = 10 5 , illustrated in Fig. 3(a-c).In panel (g), we report the largest growth rates (real part of eigenvalues) associated with the solutions, and highlight the location of simple bifurcations in green and Hopf bifurcations by a red dot in the different diagrams.As the Chandrasekhar number increases to Q = 10 3 in Fig. 3(a-c), the deflation algorithm found a rather surprising evolution of the states in this branch.Hence, while the branch emanates from the (2, 1)-eigenmode at Ra ≈ 2 × 10 4 and Q = 10 3 by breaking the Boussinesq symmetry but keeping the mirror symmetry, with two convection rolls on the magnitude velocity field, the states in the branch evolve to follow a pattern with four rolls for Rayleigh number larger than Ra ≈ 8.4 × 10 4 .The latter profile near Ra ≈ 10 5 in Fig. 3(a) is reminiscent to the (4, 1)-eigenmode depicted in Fig. 2(e).We have investigated this branch by refining the grid resolution and the continuation step size but could not discover additional branches connecting to the fourth eigenmode.To deepen the understanding of the evolution of the states in the branch, we have provided a more detailed diagram in the inlet plot of Fig. 3(a) for the magnitude of the velocity field of branch (1) in the range 8.4 × 10 4 ≤ Ra ≤ 9.4 × 10 4 .This additional diagram indicates that the correct primary branch was discovered and reported and that there is no jump to secondary branches.There exist two saddle-node bifurcations located at Ra ≈ 8.87 × 10  (d-f) Same as (a-c) but the evolution is for the steady states in the fourth primary branch (4), with the largest growth rates plotted in (h).The secondary branch (7) emanates at the two bifurcation points and connect the two primary branches.
leading to an S-shaped curve.Moreover, three pitchfork bifurcations take place in this region at respectively Ra ≈ 8.47 × 10 4 , Ra ≈ 8.89 × 10 4 , and Ra ≈ 9.22 × 10 4 , highlighted by green dots in Fig. 3(g).We have indicated these secondary branches in dashed black curves in Fig. 3(a) to emphasize that the branches have been discovered by deflation but will not be analyzed further in this work.Finally, we note the presence of a Hopf bifurcation at Ra ≈ 8.9×10 4 and a subcritical pitchfork bifurcation at Ra ≈ 6.2 × 10 4 .The secondary branch (6) arising from this bifurcation is depicted in Fig. 3(d-f) and will be analyzed later.Based on the plot of the growth rates in Fig. 3(g), we see that this branch is stable for Rayleigh numbers smaller than 6.2 × 10 4 and believe that all the bifurcations from branch (1) have been obtained.Hence, the bifurcation diagram is complete.Then, as we decrease the Chandrasekhar number towards Q = 1, we find in Fig. 1(b) that branch (1) undergoes two saddle-node bifurcations around Q ≈ 610 and Q ≈ 415.Finally, after fixing Q = 1 and decreasing the Rayleigh number, we observe in the bifurcation diagrams depicted in Fig. 1(a) that the branch originates from the (4, 1)-eigenmode at the seventh critical Rayleigh number Ra ≈ 4.7 × 10 4 , as expected by the profile of the velocity magnitude for Ra ≥ 9 × 10 4 in Fig. 3(a).
We proceed with the secondary branch (6), which emanates from branch (1) at Ra ≈ 6.2 × 10 4 and connects to branch (2) at Ra ≈ 4.47 × 10 4 .As we observe in Fig. 3(d-f), this branch has lost the mirror symmetries of its parent branches and evolves from a velocity field with three convection rolls to two convection rolls.Additionally, the branch is weakly unstable over its range of existence with a maximum growth rate of R(λ) ≈ 17, as confirmed by Fig. 3(h).
We conclude our analysis of the first and second primary branches and their secondary bifurcations by focusing on branch (2) in Fig. 4(a-c) and (g).We find that the states in branch (2) originally emanate from the (3, 1)eigenmode around Ra = 0 and Q = 1 in Fig. 1(a) and have the Boussinesq symmetry.The profile of the velocity, temperature, and magnetic fields does not change much in Fig. 4(a-c) when increasing the Chandrasekhar number as the three convection rolls of the magnitude velocity field, i.e., two smaller rolls on the sides and a FIG. 6. Evolution of the steady states in the secondary branch (7), bifurcating from branches (3) and ( 4), with respect to the Rayleigh number, where the Chandrasekhar is fixed to Q = 10 3 , and illustrated via the kinetic energy (a), potential energy (b), and magnetic energy (c).This branch features a fold bifurcation at Ra ≈ 9.17 × 10 4 , with an upper branch in the range 4.83 × 10 4 ≤ Ra ≤ 9.17 × 10 4 and a lower branch in the range 7.09 × 10 4 ≤ Ra ≤ 9.17 × 10 4 .The largest growth rates for the upper and lower branches are respectively reported in (d) and (e).
stronger in the middle, are preserved from Q = 1 to Q = 10 3 .Interestingly, we observe in Fig. 4(g) that this branch goes through a Hopf bifurcation at Q ≈ 700 and becomes stable afterwards.This indicates that increasing the strength of the magnetic field may help stabilizing solutions to the magnetic Rayleigh-Bénard problem.We then focus on Fig. 4(d-f) to study the behavior of the steady states in the branch at Q = 10 3 as the Rayleigh number decreases to find the linear limit from which the branch arises.One of the eigenvalue (depicted by a green dot in Fig. 4(h)) crosses zero as Ra ≈ 4.47 × 10 4 , which gives birth to a subcritical pitchfork bifurcation, leading to branch (6).Surprisingly, branch (2) undergoes a sequence of two saddle-node bifurcations afterwards around Ra ≈ 3.7 × 10 4 , as confirmed by the inlet plots of Fig. 4(d-f) and finally originate from the (1, 1) eigenmode by breaking the initial mirror symmetry of the eigenmode.This behavior is counter intuitive as one would have expected branch (2) to bifurcate from the (3, 1)eigenmode instead.

C. Third and fourth primary branches
This section focuses on the third and fourth primary branches arising from the conducting states at Q = 10 3 over Ra ∈ [0, 10 5 ].First, the third primary branch (3) plotted in Fig. 5(a-c) emanates from the (3, 1)-eigenmode of the conducting state through a supercritical pitchfork bifurcation at the critical Rayleigh number Ra ≈ 3.2×10 4 (cf.Fig. 2(e)).As we see in the magnitude velocity profile of Fig. 5(a), the branch has lost the mirror symmetry of the third eigenmode but kept the Boussinesq symmetry.
As the Rayleigh number increases, we observe that the velocity field evolves to feature two convection rolls located symmetrically at the bottom left and top right corners of the domain.Interestingly, this behavior is similar to the evolution of the states in the branch arising from the third eigenmode in the non-magnetic case (see [24,Fig. 7]) and does not seem to be influenced by the presence of an external magnetic field.Then, looking at the growth rates reported in Fig. 5(g), we note the presence of a supercritical pitchfork bifurcation at Ra ≈ 4.7 × 10 4 , leading to the secondary branch (7) illustrated in Fig. 6, and a Hopf bifurcation located at Ra ≈ 5.96 × 10 4 .After reaching Rayleigh number Ra = 10 5 , we study the influence of the magnetic field on steady states in this branch by decreasing the Chandrasekhar number in Fig. 1(b).We observe that the branch bifurcates from branch (5) at Q ≈ 567 after going through a saddle-node bifurcation around Q ≈ 540.
Here, branch (5), depicted in green in the bifurcation diagrams of Fig. 1(a,b) is the branch bifurcating from the (2, 2)-eigenmode at Ra (4)   c ≈ 2.3 × 10 4 and Q = 1, i.e., the magnitude velocity profiles of the steady states feature four symmetric convection rolls, with the mirror and Boussinesq symmetries preserved throughout the branch.As the strength of the magnetic field increases at Ra = 10 5 , the convection rolls become thinner in the x direction and longer in the buoyancy direction (see [64,Fig. 4.7(b)]).However, we find that this branch dies around Q ≈ 965 and is therefore not present in the bifurcation diagrams at Q = 10 3 depicted in Fig. 1(c).
Moving to the fourth primary branch (4) illustrated in Fig. 5(d-f), we see that the states in this branch emanates from the (4, 1)-eigenmode of the conducting state through a supercritical pitchfork bifurcation at Ra ≈ 5.9 × 10 4 and Q = 10 3 .Contrary to branch (3), we observe in the magnitude velocity plots of Fig. 5(d) that the Boussinesq symmetry of the eigenmode is broken at the bifurcation but the mirror symmetry is preserved.Despite the birth of the state from the (4, 1)eigenmode of the conducting state, the velocity and temperature fields of the states in the branch rapidly evolve to a pattern reminiscent to the one displayed in [24, Fig. 13(a,b)].As the Chandrasekhar number decreases in Fig. 1, branch (4) bifurcates from branch (5) through a subcritical pitchfork bifurcation at Q ≈ 330 after a saddle-node bifurcation at Q ≈ 270, similarly to branch (3).The connection between branch (4) and branch (5), which itself emanates from the (2, 2)-eigenmode at Q = 1, might explain the similarities between the velocity and temperature fields in branch (4) and the branch depicted in [24, Fig. 13(a,b)].Hence, the latter branch also bifurcates from a primary branch emanating from the (2, 2)eigenmode.Moving back to the study of branch (4) at Q = 10 3 , we observe in Fig. 5(h) that the branch is unstable throughout its range of existence (in terms of the Rayleigh number) but gives rise to the secondary branch (7) at Ra ≈ 7.1 × 10 4 through a subcritical pitchfork bifurcation and undergoes two Hopf bifurcations later on at Ra ≈ 8.97 × 10 4 and Ra ≈ 9 × 10 4 .
We now analyze the secondary branch (7) depicted in Fig. 6, which bifurcates from branch (3) at Ra ≈ 4.7×10 4 through a supercritical pitchfork bifurcation and branch (4) at Ra ≈ 7.1 × 10 4 through a subcritical bifurcation.This branch breaks the symmetry of its parent branches and features a large central convection roll with two diagonally opposite smaller rolls in the magnitude velocity field illustrated in Fig. 6, similarly to the states in [24,Fig. 14] for the non magnetic case.While the branch is unstable, it goes through a saddle-node bifurcation at Ra ≈ 9.2 × 10 4 .

D. Disconnected branch
We conclude the study of the branches of solutions to the magnetic Rayleigh-Bénard problem at Q = 10 3 and Ra ∈ [0, 10 5 ] with branch (8), illustrated in Fig. 7.This branch is highly unstable, non symmetric, and has a saddle-node bifurcation at Ra ≈ 8 × 10 4 .Observing Fig. 7(a), we see that the magnitude velocity profile of the lower branch features two large convection rolls and a smaller roll in the bottom left corner while the upper branch has one dominant convection roll and a smaller one located in the lower part of the domain, similarly to the branch reported in [24,Fig. 6].We emphasize that branch (8) is a disconnected branch in the bifurcation diagram for the Rayleigh number at Q = 10 3 and would be difficult to compute with standard bifurcation tools.This branch was obtained without the use of branch-switching techniques by following our strategy to continue branches at low Chandrasekhar numbers first in Fig. 1(b).
We first analyze the effect of the strength of the magnetic field on upper part of this branch (in Fig. 7(a)) by tracking the steady states as the Chandrasekhar number decreases to Q = 1 in Fig. 1(b).We find that the upper part of branch (8) persists throughout the range of Chandrasekhar numbers considered by going through two saddle-node bifurcations at Q ≈ 700 and Q ≈ 310.Then, after fixing the Chandrasekhar number to Q = 1 and decreasing the Rayleigh number, we observe in Fig. 1(a) that the branch finally emanates from branch (5) through a supercritical pitchfork bifurcation at Ra ≈ 3.4 × 10 3 .
Focusing on the evolution of the lower part of branch (8) as the Rayleigh number is fixed to Ra = 10 5 and the Chandrasekhar number decreases in Fig. 1(b), we observe the presence of two saddle-node bifurcations at Q ≈ 610 and Q ≈ 310.Finally, and similarly to branch (3), the lower part of branch (8) bifurcates from branch (5) around Q ≈ 567.However, the discretization size employed in this work does not allow us to characterize the nature of these bifurcations; either one transcritical bifurcation giving birth to the two branches or two distinct pitchfork bifurcations located very closely.

IV. CONCLUSION
In this work, we performed a bifurcation analysis of a two-dimensional magnetic Rayleigh-Bénard problem.By combining deflation with a continuation of steady states at low Chandrasekhar number Q, we were able to compute the entire bifurcation diagram at Q = 10 3 up to Rayleigh number Ra = 10 5 .The profusion of steady states obtained reveals the rich dynamics of the problem, similar to the classical Rayleigh-Bénard problem, with the occurrence of several pitchfork, Hopf, and saddlenode bifurcations.
Additionally, our numerical scheme allowed us to study the influence of the strength of the magnetic field on the branches and the profile of their velocity, temperature, and magnetic fields by continuing them to relate them to steady states at low Chandrasekhar number.Our findings show that the presence of a vertical magnetic field tends to privilege primary branches emanating from eigenmodes with convection rolls aligned with the direction of the magnetic field.This leads to a swapping of the eigenmodes of the conducting states as the Chandrasekhar number increases.Moreover, one of the states computed by deflation and depicted in Fig. 4 becomes stable as the Chandrasekhar number increases beyond Q ≈ 700.This demonstrates that increasing the strength of the magnetic field may help stabilizing solutions.We also conjecture that one could select certain solution profiles by modifying the direction of the background magnetic field.
A potential extension of this work would be to investigate the dependence of the steady states to the magnetic Rayleigh-Bénard problem on the magnetic Prandtl number.Finally, one could also exploit parameter-robust preconditioners [64,Sec. 4.3] to build fast numerical solvers and perform bifurcation analysis at large values of the Chandrasekhar number on finer meshes or in three dimensions.

4 FIG. 2 .
FIG. 2. (a-c) Evolution of the positive growth rates of the conducting state (2) as the Rayleigh and Chandrasekhar numbers vary in the intervals 0 ≤ Ra ≤ 10 5 and 1 ≤ Q ≤ 10 3 .(d) First ten eigenmodes of the primary bifurcations that emanate from the conducting state in the range 0 ≤ Ra ≤ 10 5 at Q = 1.The plots display the magnitude of the velocity, where blue colors indicate the zero-velocity magnitude and red colors a high magnitude.(e) Eigenmodes associated with the first five primary bifurcations that emanate from the conducting state in the range 0 ≤ Ra ≤ 10 5 at Q = 10 3 .

Fig. 2 (
Fig. 2(d) displays the velocity magnitude of the first

FIG. 4 .
FIG.4.Evolution of the steady states in the second primary branch(2) as Q increases at Ra = 10 5 , illustrated via the kinetic energy (a), potential energy (b), and magnetic energy (c).(d-f) Same as (a-c) but the evolution is displayed over Ra at Q = 10 3 .Panel (g) shows the largest growth rates (real parts of the eigenvalues) corresponding to the states in (a-c), while (h) reports the largest growth rates over Ra at Q = 10 3 , i.e., states in (d-f).

FIG. 5 .
FIG.5.Evolution of the steady states in the third branch (3) with respect to the Rayleigh number, where the Chandrasekhar number is fixed to Q = 10 3 , illustrated via the kinetic energy (a), potential energy (b), and magnetic energy (c).The corresponding largest growth rates are reported in (g).(d-f) Same as (a-c) but the evolution is for the steady states in the fourth primary branch (4), with the largest growth rates plotted in (h).The secondary branch(7) emanates at the two bifurcation points and connect the two primary branches.

FIG. 7 .
FIG.7.Evolution of the steady states in the disconnected branch(8) with respect to the Rayleigh number, where the Chandrasekhar is fixed to Q = 10 3 , illustrated via the kinetic energy (a), potential energy (b), and magnetic energy (c).This branch features a fold bifurcation at Ra ≈ 8.04 × 10 4 , with an upper branch in (a) featuring a simple bifurcation at Ra ≈ 10 5 , as indicated by the green dot.The largest growth rates for the upper and lower branches are respectively reported in (d) and (e).