<figure><img src="image_1.png" /></figure>

HAL
open science

# Chaotic edge regimes in magnetohydrodynamic channel flow: An alternative path towards the tipping point

Mattias Brynjell-Rahkola, Yohann Duguet, Thomas Boeck

► To cite this version:

Mattias Brynjell-Rahkola, Yohann Duguet, Thomas Boeck. Chaotic edge regimes in magnetohydrodynamic channel flow: An alternative path towards the tipping point. Physical Review Research, 2024, 6 (3), pp.033066. 10.1103/PhysRevResearch.6.033066. hal-04795955

HAL Id: hal-04795955

https://hal.science/hal-04795955v1

Submitted on 26 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Chaotic edge regimes in magnetohydrodynamic channel flow: An alternative path towards the tipping point

Mattias Brynjell-Rahkola*

Institut für Thermo- und Fluidodynamik, Technische Universität Ilmenau, Postfach 100565, D-98684 Ilmenau, Germany

Yohann Duguet

Laboratoire Interdisciplinaire des Sciences du Numérique – LISN-CNRS, Université Paris-Saclay, F-91400 Orsay, France

Thomas Boeck

Institut für Thermo- und Fluidodynamik, Technische Universität Ilmenau, Postfach 100565, D-98684 Ilmenau, Germany

<figure><img src="image_3.png" /></figure>

(Received 9 February 2024; accepted 5 June 2024; published 15 July 2024)

The effect of an imposed magnetic field on the flow of an electrically conducting fluid in a channel geometry is investigated numerically using high-performance temporal simulations. For the strongest spanwise magnetic fields considered, the turbulent state appears metastable, which makes the determination of the associated tipping point arduous. As an alternative, edge states, i.e., unstable states located on the state space boundary between the laminar and the turbulent basin of attraction, are investigated in detail. Their continuation smoothly leads to the tipping point where they are expected to collide with the turbulent dynamics. As the magnetic field intensity is raised, edge states become, on average, more energetic and more unstable, while their fluctuations become more chaotic, less predictable, and less symmetric. Nevertheless, their continuation allows one to accurately determine the value of the tipping point beyond which the laminar state becomes the only attractor.

DOI: 10.1103/PhysRevResearch.6.033066

## I. INTRODUCTION

Fluid flows are generically found to occur in two different regimes, i.e., laminar or turbulent. A turbulent flow is said to be subcritical if it exists despite the competition with the laminar regime, which, for the same parameter value, is linearly stable to all infinitesimal disturbances. It is, for instance, a well-known property of most wall-bounded shear flows at intermediate Reynolds numbers. Subcriticality makes parametric studies difficult because the identification of the bounding parameter values at which the turbulent state can collapse becomes relatively tedious [1,2]. This property is strongly related to nonlinear mechanisms: unlike for supercritical regimes, these threshold values are not directly connected to linear instabilities of the base flow and cannot be predicted from linear analysis only.

Let us consider as a starting point the case of a subcritical shear flow governed by a unique nondimensional parameter,

namely, the Reynolds number (Re), which is inversely proportional to the fluid viscosity. The phenomenology of the transition process is now partly understood from a dynamical systems perspective [3]. We begin by sketching the state space associated with simple low-order models [4–6] known to capture the desired dynamical behavior in a simplified way. The laminar state is a fixed point of the state space associated with the governing equations, and for low enough Re, it is the only equilibrium of the system. A new pair of fixed points emerges as a single point at a finite value of Re = Re<sub>sn</sub>, at a finite distance in state space from the laminar fixed point. The two twins separate for Re > Re<sub>sn</sub>. This *saddle-node bifurcation* leads to a new stable fixed point and its unstable counterpart, called the *edge state*, as sketched in green in Fig. 1. Conversely, if one is interested in determining the value of Re<sub>sn</sub>, it can be found by smooth continuation of *either* the stable or the unstable branch of the fixed points.

As we move to the more complex Navier-Stokes equations, a similar phenomenology exists, except, rather than fixed points, it involves more complex objects in the state space such as strange chaotic sets [7]. Similarly, the laminar state is the only equilibrium for all values of Re, yet above a given threshold in Re, two new flow regimes emerge, i.e., the turbulent and the edge regimes, that are, respectively, attracting and repelling. The turbulent regime is characterized by complex fluctuations. The repelling edge regime can also support complex dynamics, although generally much less complex than the turbulent one in terms of fractal dimension [8] as well as in energy. Both complex regimes are formally thought to be structured around networks of many finite-amplitude

<footnote>*Present address: Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, United Kingdom.</footnote>

<footnote>†Contact author: thomas.boeck@tu-ilmenau.de</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 1. Qualitative bifurcation diagrams for the HD case with respect to the Reynolds number Re (green), and for the MHD case with respect to the Hartmann number Ha (red) at fixed Re. Dashed, dotted, and solid lines indicate unstable, metastable, and stable regimes, respectively.</figcaption></figure>

solutions (and heteroclinic connections between them) that individually also emerge each at their own saddle-node bifurcation [9,10]. Edge regimes are crucial to the architecture of the whole state space since their stable manifold splits the state space into different basins of attraction [8]. Despite their instability, they can be computed using direct numerical simulation using an iterative stabilization algorithm based on bisection [11–15].

Since neither the edge nor the turbulent state can safely be expected to be fixed points, the concept of saddle-node bifurcation needs to be generalized. The collision between two sets of state space of arbitrary dimension can be topologically highly complex (see, e.g., [16]). Nevertheless, for a parametric study, it is enough to know that such a collision does occur at a given parameter value, labeled, for convenience, also $Re_{sn}$ (see, also, [17]). This value is denoted as the *tipping point* of the system: the turbulent state is no longer defined beyond that point and the system as a whole has no choice but to switch to the only available attractor left, here the laminar state. The notion of tipping point, which originates from catastrophe theory [18], rests on the analogy between Fig. 1 and the bistable dynamics observed in other disciplines, e.g., geophysical and climate sciences [19], ecology [20], or economy [21]. We refer here more specifically to the notion of the *B-tipping point* to distinguish such *bifurcation*-induced tipping from other tipping phenomena due to noise or rate of change in the system parameters [22]. The notion of an edge state also has its direct equivalent in other fields [23], for instance, in the Earth climate problem, where they are referred to as *melancholia states* and mediate between the warm and the cold “snowball” state [24]. Several early warning signs for tipping points have been proposed [25–28], among them an increased variance and a modified skewness in representative time series.

In the hydrodynamic (HD) context, although edge regimes are by definition dynamically unstable, we argue that their

computation can be beneficial to the understanding of the turbulent state itself. In particular, as Re is lowered, turbulent fluctuations become finite-lived [29] and can be missed by simulations, which can have a strong impact on the predicted value of $Re_{sn}$. Edge regimes, on the other hand, are more robust to continuation in Re, and since they are expected to merge with the turbulent regime at $Re \approx Re_{sn}$, their computation can be an alternative way to predict the value of Re at which turbulence can no longer be sustained [30]. This underlying idea has recently led to many edge calculations in various HD flow cases including boundary layers [15,31] where the concepts of edge and saddle-node bifurcation remain to be well defined.

The introduction of an additional parameter likely to dampen turbulence adds yet more complexity to this dynamical systems picture, although this has received comparatively less attention than the ordinary case with one parameter. Common examples of such damping include the Coriolis force in rotating flows [32], the buoyancy force in stably stratified flows [33], or the stress applied by the presence of polymers [34]. In the present work, we focus on another example relevant to electrically conducting fluids, namely, the Lorentz force. Such a flow is considered here between two infinite parallel plates that are electrically insulating. The flow is placed inside an externally imposed transverse magnetic field parallel to the walls. Moreover, we consider the quasistatic magnetohydrodynamic (MHD) approximation (see Sec. II A), in which magnetic self-induction is neglected. This framework applies to most laboratory and industrial flows involving liquid metal [35,36].

Unlike viscosity, which targets the smallest scales of the flow, the damping provided by the Lorentz force typically acts on all length scales [37]. In addition, whereas the viscous dissipation can, in many cases, be regarded as isotropic, the Joule dissipation promotes anisotropy by preserving and damping angular momentum, respectively, parallel and orthogonal to the magnetic field. The coherent structures in turbulent MHD flow have hence a tendency to elongate in the direction of the magnetic field and become two-dimensional [35,38]. This trend is, however, opposed by nonlinear momentum transfer, which strives to restore isotropy. As a result, the flow can be expected to remain three-dimensional for all times unless the magnetic interaction parameter, which measures the ratio between the Lorentz force and inertia, is $O(1)$ or larger. To this end, given the paramount importance of nonlinear and three-dimensional effects in transition to turbulence, we choose a relatively small magnetic interaction parameter (see Sec. II B 1). Despite the similar geometry, this turbulent flow is thus significantly different from the quasi-two-dimensional situation investigated in, e.g., [39,40].

Within the quasistatic MHD framework, the flow dynamics is governed by two nondimensional parameters, the Reynolds number (Re) and the Hartmann number (Ha) proportional to the intensity of the magnetic field. In view of the damping by the Lorentz force mentioned above, every nonzero value of Ha leads to a delay in Re of the transition process compared to the HD case.

For a fixed value of Re capable of supporting HD turbulence, there is now also subcriticality with respect to the parameter Ha: although turbulent flow can be ignited for Ha
<figure><img src="image_5.png" /><figcaption>FIG. 2. Sketch of the channel flow setup with an applied spanwise magnetic field (red arrow) and the laminar Poiseuille velocity profile (blue). The dimensions of the computational domain are indicated in units of the half-gap $h$.</figcaption></figure>

either equal to zero or sufficiently small, increasing Ha beyond a certain (Re-dependent) value leads to partial or even full relaminarization [41–43], as long as the laminar flow itself remains linearly stable. Without any prerequisite, it is hence tempting to hypothesize a globally similar tipping scenario for the transition with respect to changing Ha, at least close to the values of Ha for which full relaminarization occurs. This picture is illustrated in the red plane of Fig. 1 using saddle-node bifurcations. By analogy, the definitive collapse of the turbulent state occurs at a critical value $Ha = Ha_{sn}$. Such a value can be difficult, if not impossible, to determine by following the metastable turbulent branch because sudden relaminarization events make such a determination fully statistical [44] and hence data demanding. We are interested here in approaching that critical value $Ha_{sn}$ from below, i.e., by monitoring the parameters at which the edge regime, rather than the turbulent one, collapses. Our goal is hence to conduct edge state computations at a fixed $Re$ but for increasing values of Ha, and to find where and why edge states cease to be computable. Thereby, we have an opportunity to compute three-dimensional edge states in wall-bounded MHD flows and to perform an educated parameter study of their temporal dynamics, spatial structure, symmetries, and instability. We stress that the main theoretical difficulty compared to previous continuation studies of edge states (see, e.g., [45]) is their chaotic nature.

The paper is structured as follows. Section II contains the exact mathematical formulation and the numerical procedure. A detailed analysis of the temporal dynamics of edge states, parametrized by the Hartmann number $Ha$, is given in Sec. III. Section IV focuses on the mechanism through which the edge and turbulent dynamics collapse together for increasingly strong magnetic fields, and on a determination of the corresponding tipping point. Conclusions and outlooks are given in Sec. V.

## II. FORMULATION

### A. Governing equations

We consider the incompressible flow of an electrically conducting fluid between two infinite parallel plates. This flow is subject to a spanwise magnetic field (see Fig. 2). We use the classical shear flow convention where $x$ is the streamwise coordinate, $y$ the wall-normal, and $z$ the spanwise coordinate.

The flow is described by the incompressible Navier-Stokes equations subject to the Lorentz force. The magnetic Prandtl number $Pr_m = v/\eta$, the ratio between the kinematic viscosity and the magnetic diffusivity, is chosen to be small as appropriate to liquid metals. The magnetic Reynolds number $Re_m = RePr_m$ is thus also small, which allows the magnetic self-induction to be neglected. By further neglecting displacement currents in Ampère's law and choosing the magnetic field to be steady and irrotational, a perturbation expansion of the electromagnetic quantities in $Re_m$ shows that the Lorentz force to first order may be approximated through Ohm's law subject to a condition of charge conservation (see, e.g., [46,47] for a derivation). The nondimensionalized governing equations hence read

$$ \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -\nabla p + \frac{1}{\mathrm{Re}}\nabla^2 \mathbf{v} + \frac{\mathrm{Ha}^2}{\mathrm{Re}}(\mathbf{j} \times \mathbf{e}_z) + \chi \mathbf{e}_x, \quad (1a) $$

$$ j = -\nabla\phi + \mathbf{v} \times \mathbf{e}_z, \quad (1b) $$

$$ \nabla \cdot \mathbf{v} = 0, \quad (1c) $$

$$ \nabla \cdot \mathbf{j} = 0, \quad (1d) $$

where $\mathbf{v}$ is the velocity, $j$ is the current density, $p$ is the pressure, and $\phi$ is the electric potential. Following previous work on the HD channel flow [48,49], the characteristic velocity and length scale used for the original nondimensionalization are the laminar centerline velocity $V_c$ and the channel half-gap $h$. This defines the (hydrodynamic) Reynolds number $Re = V_c h / v$. The flow is assumed to be driven at a constant mass flow rate such that the Reynolds number is constant in time. The spatially uniform but time-dependent volume force of amplitude $\chi(t)$ is hence applied in the streamwise direction in order to maintain $V_c$ constant. The current density is nondimensionalized using $B_0 V_c / (\mu \eta)$, where $B_0$ is the strength of the magnetic field and $\mu$ the magnetic permeability of free space. This defines the Hartmann number $Ha = B_0 h (\mu \eta \rho v)^{-1/2}$, where $\rho$ is the fluid density. The walls located at $y = \pm 1$ satisfy no-slip and are chosen to be electrically insulating, i.e., $(\mathbf{j} \cdot \mathbf{e}_y)|_{y=\pm 1} = 0$.

### B. Numerical configuration

#### 1. Flow parameters

It is crucial for our interest in subcritical regimes to select well the parameter range in $Ha$ and $Re$ where turbulence (i) can occur and (ii) coexists with a linearly stable laminar regime. We begin by analyzing the first linear instability of the base flow. The HD flow possesses the common Poiseuille base solution $v_x = 1 - y^2$, $v_y = v_z = 0$, which becomes linearly unstable to spanwise-independent Tollmien-Schlichting (TS) waves at $Re \approx 5772$ [50]. As discussed in [51], spanwise-independent solutions to (1) coincide with those of HD channel flow. This can be seen by taking the curl of Fig. 1(b),

$$ \nabla \times \mathbf{j} = -(\mathbf{e}_z \cdot \nabla)\mathbf{v} = \mathbf{0}, \quad (2) $$

whose only solution compatible with the governing equations and the boundary conditions is the trivial $\mathbf{j} = \mathbf{0}$. It implies that the base flow in the MHD case is identical to its HD counterpart, their first instability occurs at the same value $Re \approx 5772$,
<table><caption>TABLE I. Friction Reynolds numbers for a selection of Hartmann numbers. Shear stresses are averaged over the two walls, and for 400 advective time units after the initial transient phase.</caption><thead><tr><th>Ha</th><th>0</th><th>5</th><th>10</th><th>15</th></tr></thead><tbody><tr><th>Re<sub>τ</sub></th><td>207.8</td><td>206.1</td><td>198.3</td><td>182.8</td></tr></tbody></table>

and below that value of Re, the base flow is linearly stable for all Ha. At the other end, the energy Reynolds number $Re_E$, below which all perturbations decay monotonically in energy, is known to vary from $Re_E = 49.6$ for Ha = 0 to $Re_E = 87.6$ for Ha = 10 [52,53]. These numbers reveal the Reynolds number range in which the flow may undergo subcritical transition. We herein set $Re = 5000$ as in Ref. [51]. To study how increasing magnetic field strengths affect the edge dynamics, the Hartmann number is varied between 0 and 20. This corresponds to an interaction parameter of $N \in [0, 0.08]$, defined as $N = Ha^2/Re$. Expressed in terms of the ratio between the Hartmann number and the laminar bulk Reynolds number ($Re_b = 2Re/3$), the case $Ha = 20$ corresponds to $Ha/Re_b = 0.006$. This is close to the critical ratio of 0.0067 reported in laminarization experiments in a small aspect ratio duct approximating the current channel geometry [41]. This suggests that a subcritical transition can be expected for the parameters under consideration, and that $Ha = 20$ will be close to the tipping point beyond which both turbulence and its edge disappear.

## 2. Computational domain

Periodic boundary conditions are assumed in the streamwise (x) and spanwise (z) directions. The two geometric parameters are the streamwise and spanwise periods, $L_x$ and $L_z$, respectively. Borrowing from previous nonlinear studies of HD shear flows in periodic domains, we select $L_x$ and $L_z$ of $O(1)$ such that the computational domain can host at least a pair of modulated streamwise streaks in the turbulent regime [54]. This hypothesis keeps computational costs manageable at the (assumed) risk of missing the spatial localization of edge states reported in larger domains [55,56]. Minimal periodic domains offer, however, the strong advantage that the statistics of the turbulent regime are well reproduced [57] and that the associated self-sustaining process is easier to untangle [58] than in larger computational domains.

With this background, we chose the size of the domain as $L_x = 2\pi$, $L_z = 2\pi/3$ in units of $h$ (see Fig. 2). The friction Reynolds number, defined as

$$ \mathrm{Re}_{\tau} = \left[ \frac{\mathrm{Re}}{2L_x L_z} \int_0^{L_z} \int_0^{L_x} \left. \frac{\partial v_x}{\partial y} \right|_{y=-1} \mathrm{dxdz} - \frac{\mathrm{Re}}{2L_x L_z} \int_0^{L_z} \int_0^{L_x} \left. \frac{\partial v_x}{\partial y} \right|_{y=+1} \mathrm{dxdz} \right]^{\frac{1}{2}}, \quad (3) $$

is listed in Table I for a selection of values of Ha. $Re_τ$ decreases with Ha, which is attributed to a reduced wall shear of the mean flow due to the damping by the Lorentz force.

It is common in wall turbulence studies to use inner units, which are the equivalent of the Kolmogorov scale in isotropic turbulence. Lengths measured in inner units have the typical superscript label “+” and are related to the usual outer units by a simple multiplication by $Re_τ$. Once expressed in inner units, the wall-parallel wavelengths of the computational domain are $L_x^+ = 1157.0-1305.4$ and $L_z^+ = 385.7-435.1$ as Ha goes from 15 to 0. The domain under consideration can thus still be viewed as “large” in comparison to the minimal domains used in Refs. [11,12,57,59], for which $L_z^+ \approx 100$ and $L_x^+ \approx 200-400$. The expectations for yet larger domains will be mentioned in the concluding section.

## 3. Numerical method

All simulations are performed with the freely available open source code NEK5000 [60], celebrated for its high accuracy and parallel performance (see [61] for a validation in the MHD case). Equations (1) are discretized using the staggered $P_N$-$P_{N-2}$ spectral element method (SEM) [62], where the domain is divided into nonoverlapping hexahedral elements, and the solution variables are expanded in high-order Lagrange interpolation polynomials. The velocity and pressure are thus represented on $N+1$ Gauss-Lobatto-Legendre and $N-1$ Gauss-Legendre quadrature nodes, respectively. This scheme has recently been extended to the $j-\phi$ quasistatic MHD formulation in [61], which enables a consistent discretization of Eqs. (1a), (1c) and Eqs. (1b), (1d).

The $x$ and $z$ directions are spanned by 32 and 16 equidistant elements, respectively, whereas the $y$ dimension is split into 32 Gauss-Lobatto-Chebyshev distributed elements. For the direct numerical simulations (DNSs), a polynomial order of $N = 11$ is used. For the HD case, the first grid point is located at a distance $y_1^+ = 0.046$ from the wall, which is comparable to the wall-normal resolution used in [63], while six points are located below $y^+ = 1$. The largest wall-normal grid spacing is $\Delta y_{max}^+ = 2.725$ (within the elements next to the channel center). In the streamwise direction, the smallest and largest grid spacings are $\Delta x_{min}^+ = 1.124$ and $\Delta x_{max}^+ = 5.570$, respectively, while the corresponding values in the spanwise direction are $\Delta z_{min}^+ = 0.748$ and $\Delta z_{max}^+ = 3.713$.

## 4. Edge-tracking algorithm

In order to track the edge state, the classical bisection algorithm [11–13] is implemented based on the values of a scalar observable, here selected as the perturbation cross-flow energy,

$$ E_{cf} = \frac{1}{4L_x L_z} \int_0^{L_z} \int_{-1}^{1} \int_0^{L_x} (v_y^2 + v_z^2) \mathrm{dxdydz}. \quad (4) $$

$E_{cf}$ vanishes exactly for laminar flow, whereas it reaches values of $O(10^{-3})$ in all turbulent regimes. The bisection algorithm recursively adjusts the amplitude of a given initial perturbation at $t = 0$, so that the corresponding time series of $E_{cf}$ stays away from both its laminar and turbulent levels for the longest possible time. When finite precision limits the refinement of the amplitude, the process can be restarted from the state of the perturbation at a later time.
<figure><img src="image_6.png" /><figcaption>FIG. 3. Edge dynamics (dark red) for different Hartmann numbers: (a) Ha = 0, (b) Ha = 5, (c) Ha = 10, (d) Ha = 12, (e) Ha = 13, and (f) Ha = 15. Trajectories that approach turbulent and laminar flow are shown with blue solid and dashed lines, respectively. (Note that the laminar fixed point corresponds to $E_{cf}$ = 0.) Green markers in (c) and (f) indicate the snapshots visualized in Figs. 7 and 8, respectively.</figcaption></figure>

This generates an arbitrary long trajectory converging to an edge state.

It is useful in practice to define threshold values of $E_{cf}$, denoted $E^{\text{lam}}$ and $E^{\text{turb}}$, where $O(10^{-3}) > E^{\text{turb}} > E^{\text{lam}} > 0$. If they are well chosen, a given trajectory experiencing values of $E_{cf}$ below $E^{\text{lam}}$ has no choice but to eventually reach the laminar point. Similarly, a trajectory experiencing values of $E_{cf}$ above $E^{\text{turb}}$ is already turbulent or reaches the turbulent state in finite time. Well-tuned threshold values lead to an important gain in computational time, while the duration of the learning phase remains negligible. Here, $E^{\text{lam}} = 10^{-6}$ for all Ha, whereas $E^{\text{turb}} = 10^{-3}$ for all Ha except Ha = 18 where a slightly lower value of $E^{\text{turb}} = 6 \times 10^{-4}$ was chosen due to the reduced energy levels of the turbulent state. In practice, for every pair of initial conditions bracketing the edge, bisection is performed until the distance between them (measured using the $L^2$-norm) falls below $10^{-4}$. A new bisection is initiated once the separation between the trajectories has diverged to $10^{-2}$. These numbers are a compromise between the computational cost and the required accuracy of the edge trajectories. Since edge states from the literature feature much less energy at small scales than their turbulent counterparts (see, e.g., [64]), the polynomial order used during the bisections is lowered to $N = 7$.

In general, the results of the bisection do not depend on the choice for an initial condition and the edge state is considered, by default, to be unique. For reproducibility purposes, the initial velocity field used in the present edge-tracking method is described in Appendix B.

### III. CHARACTERIZATION OF CHAOTIC EDGE STATES

#### A. A first exploration using energy-based observables

We begin our description of the edge states of the MHD channel with a qualitative description based on time series associated with physically meaningful observables. We first consider the cross-flow perturbation energy $E_{cf}(t)$, i.e., the observable used in the edge-tracking method.

The first parameter configuration that is investigated corresponds to the ordinary HD channel, for which $Ha = 0$. The time series associated with the converged edge state, shown in Fig. 3(a), feature long recurrent quiescent phases interrupted by rapid bursts of $E_{cf}$. Such time series, as well as the timescale of $O(10^3)$ for the bursting period, are directly reminiscent of periodic edge states computed in the same system for lower Re values [12,48,59] or in other fluid systems [65]; however, strict periodicity is lost here in favor of a relatively gentle chaotic state characterized by fluctuating periods. Exact periodic orbits and traveling wave solutions are expected to exist for nearby parameter values of $L_x$ or $L_z$, as shown recently in Ref. [66]. An interburst time is evaluated from the energy time series of Fig. 3 using adequate thresholding. The mean interburst time $t_{\text{interb}}$ is shown in Fig. 4 to decrease exponentially fast with Ha. For Ha above 13, it becomes increasingly difficult to separate rapid fluctuations in the bursts from quiescent phases and hence this region is excluded from Fig. 4.

Perhaps the most direct way to assess what these bursts in energy correspond to in physical space is to choose a given
<figure><img src="image_7.png" /><figcaption>FIG. 4. Mean interburst time ($t_{\text{interb}}$) as a function of Ha (logarithmic scale). The error bars represent 95% confidence intervals derived using Student's t-distribution and the dashed line is an exponential fit.</figcaption></figure>

line parallel to the walls (with y constant), average the flow in the homogeneous x direction, and visualize the velocity field in a (z, t) space-time diagram. We selected the streamwise velocity perturbation $v_x'$ because it is the quantity whose modulation along z is associated with the coherent structures called streaks. The edge state is distributed in space (both x and z directions) and consists, as shown in the space-time diagram of Fig. 5(a), of two high- and low-speed streaks. These streaks shift by a distance $L_z/4$ in the spanwise direction after every burst.

Similar dynamics were reported for other realizations of the HD channel flow [12,48] and for the asymptotic suction boundary layer (ASBL) [31,65]. In ASBL, the quiescent and bursting phases of this recurrence cycle were linked to counter-rotating vortices that wrap around a low-speed streak. This streak becomes increasingly wavy until it breaks up into smaller segments, before it forms again and the cycle continues. Strict periodicity of the corresponding edge state was reported in Ref. [65]. The corresponding high-amplitude periodic orbit was shown to originate from the collision of a stable and an unstable fixed point in a saddle-node infinite-period (SNIPER) bifurcation [67] as the suction velocity exceeds a critical value. The bursting dynamics should thus not be interpreted as a transient visit to an unstable fixed point, but rather as a visit to its ghost left behind by the bifurcation. The only notable difference between that case and the present channel simulation is the absence of exactly recurrent dynamics, consistent with weak chaos.

As Ha departs from zero and a weak magnetic field is applied, the above picture is initially unchanged, but the duration of the quiescent phases shortens, while the frequency of the bursts increases and their duration extends [cf. Figs. 3(a)–3(d)]. It is useful at that stage to introduce a second global observable, i.e., the longitudinal perturbation energy $E_{\text{stk}}$ defined by

$$E_{\text{stk}} = \frac{1}{4L_x L_z} \int_{0}^{L_z} \int_{-1}^{1} \int_{0}^{L_x} [v_x - (1 - y^2)]^2 \, dxdydz, \quad (5)$$

where the suffix “stk” stands for streaks. An increase in $E_{\text{cf}}$ is accompanied by a decrease in the longitudinal perturbation energy $E_{\text{stk}}$ [see Fig. 6(a)]. As Ha continues to increase, the formerly short bursts turn into extended episodes of chaotic, higher-energy fluctuations. As shown in Fig. 5(b)

<figure><img src="image_8.png" /><figcaption>FIG. 5. Space-time diagrams showing the temporal evolution of the x-averaged streamwise velocity perturbation $\langle v_x' \rangle_x(z, t)$ at $y = -0.7$, for (a) $Ha = 0$, (b) $Ha = 10$, and (c) $Ha = 15$.</figcaption></figure>

for $Ha = 10$, this change is associated with a change in the spanwise shifts of the streaks: the shifts no longer occur over a well-defined fraction of $L_z$ and they instead become erratic. Three-dimensional visualizations of the streamwise velocity contours in Fig. 7 show that coherent streaks are still present both in the quiescent and the active parts of the cycle. For Ha above approximately 10, the quiescent phase becomes more clearly evaded, and above $Ha = 15$, there is no longer any proper quiescent episode, only rapid fluctuations. The motion of the streaks can no longer be described as recurrent, but as spatially and temporally chaotic with brief quiescent episodes [cf. Figs. 3(d)–3(f)]. The repercussions on the state portraits are clear in Fig. 6(a): trajectories have progressively switched towards another part of the reduced state space characterized by higher values of $E_{\text{cf}}$ (including a higher ground level) but less variations in $E_{\text{stk}}$. For $Ha = 15$, the space-time diagram
<figure><img src="image_9.png" /><figcaption>FIG. 6. State portraits of the edge state for different Ha. (a) Cross-flow energy ($E_{cf}$) plotted against streak energy ($E_{stk}$). (b) Energy input ($I$) plotted against total dissipation ($D_v + D_μ$). (c) Joule dissipation ($D_μ$) plotted against viscous dissipation ($D_v$). The color coding is the same for all state portraits.</figcaption></figure>

of Fig. 5(c) shows a clearly erratic path for the streaks with no clear bursting phase, despite the fact that the fluctuations of the edge state still span a decade in $E_{cf}$. The three-dimensional visualizations of Fig. 8 show that the lowest-$E_{cf}$ episodes still maintain spatial coherence in the form of a solitary meandering streak localized near one wall only. In contrast, even though the fluctuations remain located close to one wall, the more energetic snapshots show less spatial coherence and a striking lack of spatial symmetry compared to lower Ha.

As discussed in Sec. I, the presence of the Lorentz force introduces a second dissipation mechanism besides the common viscous dissipation. Starting from (1), the energy balance now reads

$$ \frac{\partial E}{\partial t} = I - D_v - D_μ, \qquad (6) $$

where the different terms

$$ E = \frac{1}{2} \int_{0}^{L_z} \int_{-1}^{1} \int_{0}^{L_x} \mathbf{v} \cdot \mathbf{v} \, dxdydz, \qquad (7a) $$

$$ I = \int_{0}^{L_z} \int_{-1}^{1} \int_{0}^{L_x} \chi(\mathbf{v} \cdot \mathbf{e}_x) \, dxdydz, \qquad (7b) $$

$$ D_v = \frac{1}{\text{Re}} \int_{0}^{L_z} \int_{-1}^{1} \int_{0}^{L_x} \nabla \mathbf{v} : \nabla \mathbf{v} \, dxdydz, \qquad (7c) $$

$$ D_μ = \frac{Ha^2}{\text{Re}} \int_{0}^{L_z} \int_{-1}^{1} \int_{0}^{L_x} \mathbf{j} \cdot \mathbf{j} \, dxdydz \qquad (7d) $$

are, respectively, the kinetic energy, the energy input by the streamwise driving force, the viscous dissipation, and the Joule dissipation linked to MHD effects. State portraits of the input $I$ vs the total dissipation $D_v + D_μ$ of the edge states are shown in Fig. 6(b). The trajectories spend most of their time in the vicinity of the lower-left corner of the figure and only explore the central portions during the bursts. During these events, trajectories always stay above the diagonal (where equilibria lie), indicating that the energy input always dominates over the dissipation during these events and that the bursts are always out of equilibrium. As Ha continues to increase beyond 13, excursions away from the diagonal are less and less marked.

In Fig. 6(c), the Joule dissipation is directly plotted versus the viscous dissipation in another state portrait. The Joule dissipation always appears smaller than the viscous one by one order of magnitude. Both quantities burst in synchrony with $E_{cf}$ and the energy input $I$. Except for these excursions, the flow becomes, on average, more dissipative as the edge state turns more chaotic.

From all global observables investigated until here ($E_{cf}$, $E_{stk}$, energy budget terms), although there is a trend towards higher energies and less fluctuations as Ha increases, no sharp transition has been noted. The apparent chaoticity does not make a bifurcation study trivial. A costly computation of the Lyapunov exponents would be required, both to locate bifurcations and to confirm the chaotic property itself. In the next section, we adopt another strategy by moving to physically relevant symmetry-based observables.

## B. Symmetry indicator

The physical space visualizations of Figs. 7 and 8 suggest, at least during the higher-energy intervals, a loss of spatial coherence for stronger magnetic fields. In this section, a symmetry-based indicator is constructed in order to quantify this trend in the hope of identifying (or not) a bifurcation at a finite value of Ha. The symmetry of interest apparently persists during the lower-energy episodes, at least in both mentioned figures for Ha = 10 and 15, respectively. A careful look at these figures suggest to consider the *shift-and-reflect*
<figure><img src="image_10.png" /><figcaption>FIG. 7. Snapshots of the edge state for Ha = 10 at times marked with green circles in Fig. 3(c). (a) t = 4494 quiescent phase and (b) t = 4680 bursting phase. Isosurface of $v_x = 0.3$ (blue) at the lower wall plotted together with contours (white and red) of $\omega_x = \pm 0.1$. For guidance, isolines of $v_x$ separated by 0.1 units ranging from 0.0 to 0.9 are shown on two-dimensional planes at $x = 2\pi$ and $z = 0$.</figcaption></figure>

symmetry defined for a three-dimensional vector field $v$ by

$$ \|v - S(b)v\| = 0, \qquad (8) $$

where the shift-and-reflect operator is defined as

$$ S(b) : [v_x, v_y, v_z](x, y, b+z) \\ \rightarrow [v_x, v_y, -v_z] \left( x + \frac{L_x}{2}, y, b-z \right), \qquad (9) $$

and $z = b$ is the plane of symmetry. We chose the $L^2$-norm scaled with the volume of the domain such that

$$ \|f\| = \sqrt{\frac{1}{2L_x L_z} \int_0^{L_z} \int_{-1}^1 \int_0^{L_x} f \cdot f \, dxdydz, \qquad (10) $$

for some vector $f$. This same symmetry was identified in the primary instability of $x$-independent streaks in several HD shear flows, including channel flow and ASBL [65,68]. The definition of $S$ given in [68] is generalized here to reflections around some arbitrary spanwise coordinate $z = b$. We

introduce the subsets

$$ \mathcal{X}_b \equiv \{v : \|v - S(b)v\| = 0\} \qquad (11) $$

of velocity fields, and their union

$$ \mathcal{X}_u = \bigcup_{0 \le b < L_z} \mathcal{X}_b. \qquad (12) $$

Our goal is to numerically evaluate the instantaneous distance of a given edge trajectory to all the symmetry subspaces $\mathcal{X}_b$, and use the time-averaged shortest distance to $\mathcal{X}_u$ as a new scalar quantity likely to quantify the degree of symmetry of an edge state.

We note that the shift-reflect operator in (9) is independent of the $y$ coordinate. Given the size of the three-dimensional data sets, the workload associated with probing the $S$-symmetry property in time and Ha is reduced if the flow is averaged across the channel height. Let us consider the $y$-averaging operation $\langle \cdot \rangle_y = (1/2) \int_{-1}^{1} (\cdot) dy$ and define the associated two-dimensional (2D) measure of a

<figure><img src="image_11.png" /><figcaption>FIG. 8. Snapshots of the edge state for Ha = 15 at the times marked with green circles in Fig. 3(f). (a) $t = 3288$ visit to a traveling wavelike state and (b) $t = 3408$ chaotic state. Isosurface of $v_x = 0.3$ (blue) at the lower wall plotted together with contours (white and red) of $\omega_x = \pm 0.3$. For guidance, isolines of $v_x$ separated by 0.1 units ranging from 0.0 to 0.9 are shown on two-dimensional planes at $x = 2\pi$ and $z = 0$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 9. Time series of the pseudodistance $d(t)$ from the edge state to $X_u$ (blue, left ordinate), compared to $E_{cf}(t)$ (dark red, right ordinate) for (a) Ha = 10 and (b) Ha = 15. Local minima of $d$ are marked using gray vertical dashed lines. (c) Time-averaged value of $d$ as a function of Ha (full circles) together with an exponential fit (dashed line). Error bars represent ± one standard deviation.</figcaption></figure>

three-dimensional vector field $f$,

$$ \|f\|_{2D} = \sqrt{\frac{1}{L_x L_z} \int_0^{L_z} \int_0^{L_x} \langle f \rangle_y \cdot \langle f \rangle_y \, dxdz.} \quad (13) $$

Note that $\| \cdot \|_{2D}$ is not strictly a norm. We then define the two-dimensional residual as

$$ \zeta(b) = \|v - S(b)v\|_{2D}. \quad (14) $$

In order to approximate the shortest distance from the edge state to the set $X_u$ [i.e., the union of all $S(b)$-symmetric subspaces], the one-dimensional minimization problem

$$ d = \min_{0 \le b < L_z} \zeta(b) \quad (15) $$

is solved at each discrete time using Brent's method of parabolic interpolation. Due to the wall-normal integration, $d$ does not correspond exactly to the shortest state space distance from the edge state to $X_u$, unless the components of the residual vector $r = v - S(b)v$ satisfy $\langle r_i^2 \rangle_y = \langle r_i \rangle_y^2$ for all $i \in \{x, y, z\}$. Nevertheless, since $d$ vanishes on all $X_b$'s, the solution to (15) can be regarded as a practical approximation of the desired state space distance.

Time series of $d(t)$ have been computed for all Ha. Two values, Ha = 10 and Ha = 15, are shown in Figs. 9(a) and 9(b), plotted together with the time series of $E_{cf}(t)$. The two signals look visually strongly anticorrelated with a local minimum of $d(t)$ immediately preceding every burst [see the dashed vertical lines of Figs. 9(a) and 9(b)]. In spite of this, here we will analyze $d(t)$ because it admits a more direct interpretation from a state space perspective.

From these time series, a temporal average can be directly computed using the time-averaging operator defined by $\langle \cdot \rangle_t = T^{-1} \int_{t_0}^{t_0+T} (\cdot) dt$. The mean quantity $\langle d \rangle_t$ is plotted against Ha in Fig. 9(c), where $t_0$ is an arbitrary time instance after the initial transient phase. The values of $\langle d \rangle_t$ are accompanied by error bars that correspond to one standard deviation in either direction. The edge state drifts away from $X_u$ as Ha increases. An exponential fit $\langle d \rangle_t \sim \exp(0.167Ha)$ can be applied over the whole range of values of Ha under study.

There are several teachings from this fit. The quantity $\langle d \rangle_t$ is never strictly zero even in the HD case Ha = 0 and it increases by almost two decades over the interval $0 \le Ha \le 20$.

The various $X_b$'s cease to be the relevant invariant subspaces for sufficiently large Ha. For Ha $\gtrsim$ 15, since the distributions shown in Fig. 10 are approximately Gaussian, the error bars corresponding to the standard mean deviation can be interpreted as the approximate distance from $X_u$ reached for approximately 68% of the time. In particular, this continuous and monotonic trend leaves no room for a symmetry-breaking bifurcation. The loss of symmetry of the edge state is hence gradual with increasing Ha over the whole range of values investigated, starting from Ha = 0 which is itself not perfectly symmetric.

Despite the shortness of the edge signals due to the cost of carrying out bisection coupled with a DNS, a statistical analysis of the time series of $d(t)$ is possible beyond the analysis of the mean values. It is likely that convergence of the relevant probability density functions (PDFs) would improve if the edge tracking was continued. PDFs of $d$ are plotted in Fig. 10 for various values of Ha, normalized by their maximum for easier reading. The PDFs are highly skewed towards the origin for vanishing to small Ha, but tend to a more symmetric, approximately normal distribution as Ha increases. This is interpreted as an indication that the edge dynamics become increasingly random with Ha.

<figure><img src="image_13.png" /><figcaption>FIG. 10. Probability density function (PDF) of the distance measure $d$ from the edge state to $X_u$, i.e., the union of symmetric subspaces $X_b$, for different Hartmann numbers Ha. To facilitate the comparison between the distributions, the PDFs are normalized to have $\max(PDF) = 1$.</figcaption></figure>
Eventually, the probability that the edge trajectory visits $X_u$ becomes vanishingly small for $Ha \ge 16$, whereby the edge state can be regarded as nonsymmetric at all times. The role of the Lorentz force in this gradual loss of symmetry can be justified based on physical arguments. The occurrence of streamwise rolls and streaks that are invariant under a shift-and-reflect symmetry is central to the self-sustaining process (SSP) explaining turbulent [4,58] as well as edge [64,68] solutions. The Lorentz force modifies the force balance at play in the SSP. From Eqs. (1a) and (1b), it is proportional to $-(\partial_y\phi + v_x)e_x + (\partial_x\phi - v_y)e_y$. This implies that the streamwise and the wall-normal velocity components will be increasingly affected as Ha becomes larger. As a result, edge states are forced out of the shift-and-reflect symmetric subspace, which becomes less frequently visited. This behavior relies on three-dimensional effects and is in sharp contrast to the quasi-two-dimensional MHD channel flow, where three-dimensional mechanisms are precluded and the edge state corresponds to a two-dimensional traveling wave [40] unrelated to the present chaotic solutions.

### C. Instability and predictability horizon

Edge states are, by construction, linearly unstable, yet their rate of instability has seldom been reported in the literature. In the simplest cases where they are fixed points of the original system, the governing equations are routinely linearized around such a fixed point and standard algebra methods are used to compute the corresponding eigenvalues. Only one of these eigenvalues is positive [64,69]. If edge states are traveling waves, they can be transformed into fixed points by symmetry reduction, while in the case of periodic orbits, Floquet analysis is the adequate tool. In the generic case where edge states are chaotic, eigenvalues and Floquet exponents are generalized into Lyapunov exponents, which, however, are notoriously difficult to compute and to converge for large-scale systems of partial differential equations [70]. Chaoticity along the edge trajectory implies that several Lyapunov exponents must be positive, although we are typically only interested in the rate of instability away from the edge manifold [71]. We decided instead to quantify here the instability of the chaotic edge state using the related concept of a finite-size Lyapunov exponent (FSLE). It was originally introduced to determine the predictability time for finite-amplitude perturbations in chaotic or turbulent systems [72,73] and is sometimes also referred to as the finite-amplitude Lyapunov exponent [74]. Although the FSLE concerns the growth of perturbations at different scales, it coincides with the largest Lyapunov exponent for perturbation amplitudes (and short trajectories) that are sufficiently small so that the influence of nonlinearities can be neglected [72–75]. While the FSLE measures the average time needed for nearby state space trajectories to separate by a given amount, the (largest) Lyapunov exponent measures the separation of nearby trajectories over a predefined time interval. For a detailed account of the numerical calculation of the FSLE, the reader is referred to Appendix A.

The average FSLE of the edge states for different Hartmann numbers, denoted as $\lambda$, is plotted in Fig. 11 together with its fluctuations. The exponent $\lambda$ is observed to grow with Ha. The data suggest that $\lambda$ is well approximated by a fit

<figure><img src="image_14.png" /><figcaption>FIG. 11. Finite-size Lyapunov exponent $\lambda$ as a function of Ha (open circles), and its fluctuations (colored dots). The dots are colored according to the instantaneous cross-flow energy $E_{cf}$ of the edge states. The dashed line is a power-law fit of the form $c + bHa^a$ with $a = 2.31$, $b = 2.86 \times 10^{-5}$ and $c = 2.67 \times 10^{-2}$.</figcaption></figure>

of the form $c + bHa^a$, with $a = 2.31$ ($b = 2.86 \times 10^{-5}$ and $c = 2.67 \times 10^{-2}$). The temporal fluctuations of the FSLE are also shown in Fig. 11, using dots colored by the value of $E_{cf}$. For $Ha \le 5$, the quiescent phase (characterized by low $E_{cf}$) appears more unstable than the bursting phase, which may seem counterintuitive. We also computed the predictability time $T_p$ resulting from an integration of $\lambda^{-1}$ over all relevant scales, and deduced from the above fit for $\lambda$ that $T_p \sim Ha^{-2.31}$ (see, again, Appendix A). This implies that the predictability time decreases with Ha and, correspondingly, that the complexity of the edge dynamics can be expected to increase with Ha. This matches well with the conclusions from the time series of Sec. III A and the associated visualizations.

Note, for comparison, that in former HD studies of plane Couette flow, the instability rate of the edge states was computed as an eigenvalue in the case where the edge state is a fixed point [64]. An asymptotic dependence $Re^{-0.48}$ of the unstable eigenvalue $\lambda$ on the Reynolds number $Re$ was reported. Although the comparison with the present case rests on different governing parameters, respectively, Ha and $Re$, a common trend emerges: the edge state is most unstable close to the saddle-node bifurcation in which it is created.

### IV. APPROACHING THE TIPPING POINT

In this section, we focus both on the turbulent and the edge dynamics as the global tipping point of the system is approached. This parameter range corresponds to the largest values of Ha where the turbulent and edge regimes still exist. The main hypothesis, as exposed in Sec. I, is that parameter continuation along the turbulent branch is fraught with difficulties, whereas continuation along the edge branch is smooth. This allows for an alternative way to determine the location
<figure><img src="image_15.png" /><figcaption>FIG. 12. Time series of $E_{cf}$ in the turbulent regime. Presence of supertransients [76] for $Ha \ge 17$.</figcaption></figure>

of the tipping point. This section is devoted to a quantitative confirmation of these hypotheses.

### A. Before the tipping point: Metastability of the turbulent state

The main focus shifts here towards the turbulent dynamics. For an account of turbulence in the Hartmann channel and the present configuration with a spanwise magnetic field, the reader is referred to Refs. [42,77]. The concept of tipping point is entirely linked with the expectation that the turbulent flow has to laminarize beyond some magnetic field intensity. Channel flow laminarization under the influence of a streamwise magnetic field was investigated in [78]. It was reported that the critical Ha for laminarization increases with both Re and domain size.

In the present work, a turbulent parameter study is performed for different Ha starting from the synthetic perturbation in Appendix B, while keeping all other parameters fixed. The time series for a set of Hartmann numbers are plotted in Fig. 12. For $0 \le Ha \le 16$, turbulence appears as sustained over the time (necessarily finite) of the simulations, which was pushed for some parameter values up to $O(10^4)$ time units. In this Ha range, no laminarization event is observed after turbulence is initiated in the domain.

By a minute increase to $Ha = 17$ (and $Ha = 18$), sudden relaminarization events are observed, without any clear warning. The time series shown in Fig. 12 suggests that these transient regimes are type-II supertransients [76]. Such supertransients are common in HD studies [2] and have also been reported in either simulations or models of finite magnetic Reynolds number MHD turbulence [79–81]. Only a detailed statistical analysis of the relaminarization times as a function of Ha can tell us for which value of Ha mean lifetimes are finite or diverge [82]. Nevertheless, such a task would be computationally extremely costly in terms of number of simulations and it is clearly outside the scope of this article. From a dynamical systems point of view, the appearance of supertransients suggests the possibility for a *boundary crisis*, in which an attractor turns into a chaotic saddle that allows trajectories to escape to another attracting state [83,84]. The present data suggest a critical value of $Ha_{bc} \le 17$ for such a bifurcation. We stress that this prediction, as well as the validity of the whole boundary

<figure><img src="image_16.png" /><figcaption>FIG. 13. Time-averaged square root of the cross-flow energy, $\langle\sqrt{E_{cf}}\rangle_t$, plotted vs Ha for the turbulent states (squares) and the edge states (circles). The error bars represent $\pm$ one standard deviation. The axis $\langle\sqrt{E_{cf}}\rangle_t = 0$ corresponds to laminar flow. With reference to Fig. 12 and the definition of the edge state [8], stable and unstable/metastable states are plotted with full and open symbols, respectively. To estimate the value of Ha for the saddle-node bifurcation, a quadratic function is fitted to the data for $Ha \ge 16$. The critical value is estimated as $Ha_{sn} \approx 19$. Furthermore, power laws are fitted to the data as a guide to the eye. The fits are plotted using a solid line (stable states), a dotted line (metastable states), and a dashed line (unstable edge states).</figcaption></figure>

crisis scenario, are potentially influenced by the finite observation time and, more importantly, by the computational domain.

Continuing further to $Ha = 19$, only brief episodes of chaotic motion are observed. Attempts to initiate turbulent flow from amplitudes $10^{-3} \le A \le 10^{-1}$ [see Eq. (B2)] were numerous and many failed. The turbulent fluctuation level still seems to be reached, a sign that the chaotic saddle still exists but that its basin of attraction has shrunken considerably.

### B. Estimation of the tipping point

The evidence for supertransients in Sec. IV A is a first obstacle to the continuation of the turbulent set with Ha. Besides there is further evidence from the data that the mean lifetimes shrink with increasing Ha. Therefore, even if a rigorous statistical analysis of the lifetimes is, in principle, achievable, we do not, in general, expect the variations of the lifetime with Ha to be a faithful indicator of where the expected tipping point lies.

#### 1. Expectations from normal form theory

We instead base our prediction of the tipping point on the graphical expectation that the bifurcation diagram for the mean of a well-chosen observable has to feature a saddle-node bifurcation. The canonical normal form for a saddle-node bifurcation is typically quadratic in the amplitude of the main variable [67]. We therefore choose to use, rather than the observable $E_{cf}$ which is a quadratic in the velocity field, its square root $\sqrt{E_{cf}}$. The time-averaged values $\langle\sqrt{E_{cf}}\rangle_t$ are reported versus Ha in Fig. 13, together with error bars based on the standard mean deviation in the time
series. Both the edge branch and the turbulent branch are reported, although the data for the turbulent branch is, consistently with the presence of supertransients, not exploitable close enough to the tipping point. A least-squares fit for $Ha \ge 16$ shows that the data conform locally to the quadratic normal form,

$$Ha = a_2 \langle \sqrt{E_{cf}} \rangle_t^2 + a_1 \langle \sqrt{E_{cf}} \rangle_t + a_0, \quad (16)$$

with $a_2 = -32572.82$, $a_1 = 902.43$, and $a_0 = 12.79$. From this fit, the critical Hartmann number $Ha_{sn} \approx 19.04 \pm 0.01$ is deduced. Above this threshold, only the laminar solution exists as an attracting state. Away from the critical Ha, both the upper and the lower branches are well described by power laws in ($Ha_{sn} - Ha$).

Note that the demonstration above is merely a proof of concept: the present method for determining the tipping point might seem superfluous in the present case because the metastable turbulent state is still within reach. However, in the presence of shorter turbulent lifetimes, only the edge branch might be quantitatively exploitable. Clearer evidence may be obtained by attempting predictions based on early warning signals.

## 2. Early warnings

We now consider the ability of early warnings [25–28] to predict the same value $Ha_{sn}$ of the tipping point. We consider, in particular, the turbulent time series of the scalar $L^2$-norm scaled by the domain volume,

$$|f| = \sqrt{\frac{1}{2L_x L_z} \int_0^{L_z} \int_{-1}^{1} \int_0^{L_x} f^2 \, dxdydz}, \quad (17)$$

applied to both the velocity and the vorticity components (with the base flow subtracted). The change in the standard deviation of these quantities relative to HD is plotted in Fig. 14(a) versus $Ha$. All variables show an increase in their standard deviation as $Ha \to Ha_{sn}$. (This trend is also visible in the error bars of Fig. 13.) The relative increase is most prominent for the streamwise velocity perturbation $|v_x'| = |v_x - (1 - y^2)|$, which shows a monotonically growing standard deviation as the tipping point is approached. The slight reduction in the other quantities at $Ha = 18$ may possibly stem from an insufficient averaging time due to the short turbulence lifetime (see Fig. 12). While the data is a *posteriori* consistent with the observed phenomenology and the previous prediction, it unfortunately does not allow for a precise estimation of $Ha_{sn}$.

Next, the change in the skewness of the same quantities is shown in Fig. 14(b) versus $Ha$. This quantity fluctuates considerably with respect to the bifurcation parameter, which makes it difficult to identify any clear trend. The fluctuations possibly intensify as the system approaches the bifurcation point, but the results are uncertain and here there is an apparent risk for *false warnings* or missed tipping. Again, it is not excluded that this failure is linked to the limited amount of data.

We conclude by noting that the traditional early warning methods based on *turbulent* statistics, at most, may signal the approach of a tipping point, but not how far away from it the system is. In that respect, the predictive use of these indicators

<figure><img src="image_17.png" /><figcaption>FIG. 14. Early warning statistics for the tipping point based on the statistics of various state variables (components of the velocity and vorticity perturbation). (a) Normalized standard deviation vs Ha; (b) skewness vs Ha.</figcaption></figure>

is limited. By contrast, the previous method documented in Fig. 13, based on a fit of mean observable values in the edge regime according to the standard normal form of a saddle-node bifurcation, yields a better quantitative estimation of that tipping point.

## C. Around the tipping point: Physical space visualizations

This section is devoted to a comparison of typical velocity fields of both the turbulent and edge regimes around the tipping point. The aim here is to illustrate, using three-dimensional rendering of the velocity field, that the edge and the turbulent regimes continuously connect at the tipping point. The edge and the turbulent regimes are depicted in Fig. 15, respectively, in the lower and upper rows, with the value of $Ha$ increasing from left to right. This organization parallels that of the bifurcation diagram in Fig. 13.

In Figs. 15(a)–15(c), instantaneous snapshots of the turbulent state are shown for $Ha = \{14, 16, 18\}$. For $Ha = 14$, turbulence is space filling; it occupies both walls as well as large portions of the bulk during the entire simulation. This is qualitatively indiscernible from HD turbulence at high $Re$. Upon increasing the Hartmann number to $Ha = 16$, the turbulent regions occupy a smaller portion of the domain and the turbulent activity lies only intermittently near both walls, with frequent switches to one wall only, a situation already described at lower $Re$ in HD channel turbulence [11]. Numerous wall switchings are also witnessed. When moving to $Ha = 18$, turbulence stays near one given wall for all times after the initial transient phase and it also displays a trend towards spanwise localization. No transition
<figure><img src="image_18.png" /><figcaption>FIG. 15. Snapshots of the turbulent state (upper row) and the edge state (lower row) at arbitrary times. (a),(d) Ha = 14, (b),(e) Ha = 16, (c),(f) Ha = 18. Isosurface of $v_x = 0.3$ (blue) plotted together with contours (white and red) of (a)–(c) $\omega_x = \pm 1.2$ and (d)–(f) $\omega_x = \pm 0.4$. In similarity to Figs. 7 and 8, isolines of $v_x$ separated by 0.1 units ranging from 0.0 to 0.9 are shown on two-dimensional planes at $x = 2\pi$ and $z = 0$.</figcaption></figure>

between the two walls was observed. The flow eventually laminarizes.

In Figs. 15(d)–15(f), instantaneous snapshots of the edge state are shown for Ha = {14, 16, 18}. The velocity field is streaky and always located on one wall (note that both walls are equivalent by symmetry). In addition, it shows strong signs of spanwise localization, especially for lower Ha. Strict localization was reported for all edge states of planar shear flows for sufficiently large $L_z$ [31,48].

Comparison between the turbulent and the edge snapshots in Fig. 15 is delicate because the energy levels differ between the two branches, and thus the values of the isolevels of, e.g., streamwise vorticity that are displayed need to be adjusted depending on the flow state. Starting from the HD edge state, moving along the lower branch of Fig. 13 around the fold, back along the upper branch towards the HD turbulent state, the flow progressively loses all its localization properties: it evolves from two low-speed streaks near one wall to full-fledged turbulent activity near both walls as well as in the channel bulk. Visual comparison between the upper and the lower rows fully supports the concept that the two flow regimes approach each other with increasing Ha.

### D. Beyond the tipping point

Beyond the tipping point, no turbulent or edge dynamics is possible, and the only attracting solution left in the state space is the fixed point corresponding to the laminar Poiseuille flow. In order to probe the validity of the tipping point predicted at $Ha_{sn} \approx 19$, a bisection is initiated at $Ha = 20 > Ha_{sn}$ in the same conditions as in Sec. III. The corresponding time series of $E_{cf}$ are shown in Fig. 16(a) where, unlike the curves in Fig. 3, no “transient phase” is excluded. Initially, there appears

to be a turbulent and laminar state to bisect between, with levels of $E_{cf}$ consistent with those found earlier before $Ha_{sn}$. However, as the bisection progresses, the value of the upper threshold $E^{turb}$ (defined in Sec. II B 4) needs to be repeatedly readjusted. Figure 16(b) reveals that *all* trajectories that are generated eventually laminarize, no matter how large the initial $E_{cf}$ value and the initial growth. The time series never settle to a constant level as with the supertransients found for $Ha_{sn} \lesssim 19$. The situation becomes even more pronounced for $Ha \approx 30$ (not shown), where all perturbation amplitudes plotted in Fig. 16(b) have decayed to $O(10^{-7})$ around $t \approx 150$. Consequently, no such thing as an equilibrium turbulent set (whether an attractor or a chaotic saddle) can be claimed here. In these conditions, the bisection algorithm cannot select a trajectory fulfilling the properties of an edge trajectory. This failure is interpreted as a sign that the edge state and its turbulent counterpart no longer exist, which confirms that $Ha = 20$ is beyond $Ha_{sn}$.

### V. CONCLUSIONS AND FUTURE DIRECTIONS

The present study has addressed the dynamical competition between the laminar and the turbulent regimes in MHD channel flow when the magnetic field is oriented spanwise. For a fixed Reynolds number and a given periodic computational domain, the Hartmann number $Ha$ assumes the role of a bifurcation parameter. The two competing regimes being stable (or metastable), another regime of unstable edge states is known to exist at the boundary between the two attraction basins. Edge states have been computed here using high-performance DNS, coupled together with the standard bisection algorithm.
<figure><img src="image_19.png" /><figcaption>FIG. 16. Time series of $E_{cf}(t)$ beyond the tipping point for $Ha = 20$. (a) Trajectories reaching the predefined (but repeatedly readjusted) threshold levels $E^{turb}$ and $E^{lam}$. (b) Trajectories with initial perturbation amplitudes $A$ in (B2) logarithmically distributed between $10^{-3}$ (dark blue) and $10^{-1}$ (yellow). All trajectories eventually laminarize without reaching any equilibrium level, making it impossible to define a turbulent level and hence to select an edge trajectory in (a).</figcaption></figure>

The edge state evolves continuously, without any marked bifurcation, from an almost periodic state with recurrent quiescent and bursting phases to a more random and more unstable chaotic state as $Ha$ is increased from 0 to 18. Beyond a tipping point located at $Ha_{sn} \approx 19$, this edge state ceases to exist, while the turbulent state shows signs of metastability already from $Ha \approx 17$ on. It was hypothesized that the edge state and the turbulent state, independently of their stability, disappear together in what constitutes a generalization of a saddle-node bifurcation to chaotic regimes. All of the presented results are consistent with this hypothesis. A large part of this paper is devoted to a characterization of the dynamics close to the corresponding tipping point and to the methodology for estimating the tipping parameter value. Note that the saddle-node bifurcation hypothesis, although common in HD transition studies, is only the simplest scenario compatible with a subcritical transition to turbulence. More complex scenarios exist; see, e.g., Ref. [30] where the branch of preturbulent states bifurcates from the edge branch before its own saddle-node bifurcation.

The present investigation was performed for one single value of the Reynolds number $Re$. It would be interesting to know how the tipping value $Ha_{sn}$ depends on $Re$. The duct experiments summarized by Branover [41] have suggested a scaling of the form $Ha_{sn} \sim Re$ over a wide range of values of $Re$, which includes the present case. Numerical simulations in the channel geometry with a spanwise field and no spanwise walls point towards a different scaling, instead involving the

interaction parameter $N$ [42,85]. An accurate numerical verification would be a first step towards a better characterization of tipping phenomena in MHD flows.

Another line of research, suggested by Figs. 1 and 13, would be to determine the exact coherent structures transiently visited by the edge state near the tipping point [see Fig. 9(b)] and continue these to the upper branch, possibly towards the $Ha = 0$ limit. Such homotopy methods, making use of a second governing parameter (here, $Ha$) to uncover hidden solutions of the original Navier-Stokes equations, were popular in the 1990–2000s after the discovery of the first unstable steady solution of plane Couette flow [86]. To our knowledge, homotopy based on the introduction of MHD effects has not been attempted so far.

Another important generalization of the present work should deal with larger computational domains. It is well known that in many wall-bounded shear flows, large-scale intermittency and laminar-turbulent patterning is possible, characterized by large length scales requiring computational domains or experimental setups of size at least 100h. It was demonstrated in Ref. [43] that the addition of a damping force, such as the Lorentz force in this work, shifts the $Re$ regime of intermittency and patterns upwards. What becomes of the tipping point in this context is so far unknown. In the limit of infinite domains (the so-called thermodynamic limit [1]), the dynamical systems picture has to be replaced by a phase-transition formalism. Recent work suggests that the phase transition, governed by $Re$ as an order parameter, could switch from a continuous to a discontinuous type as a second governing parameter is added [87]. However it is yet unknown whether fluid flows influenced by Lorentz forces enter this category.

Finally, the approach suggested in this paper to tackle the tipping point of the system, focusing on the edge branch rather than on the more energetic one, can be generalized to many physical systems with bistability properties. This approach remains more computational than experimental since unstable states become the main protagonists. The main prerequisite is to have deterministic governing equations, such that the concept of an edge state stays well defined. Extensions of these concepts to stochastic systems require more formal progress as to how unstable states can be defined and determined.

## ACKNOWLEDGMENTS

M.B.R. and T.B. acknowledge financial support from The German Research Foundation (DFG) under Project No. 470628784. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. [88] for funding this project by providing computing time on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre [89]. Post-processing and data visualization were performed at the local computing center of TU Ilmenau. M.B.R. also thanks LISN-CNRS for the hospitality during his research visit in 2023.

## APPENDIX A: COMPUTATION OF FINITE-SIZE LYAPUNOV EXPONENTS (FSLEs) FOR EDGE STATES

We describe here how the mean FSLE and its fluctuations are determined for the edge states directly as a by-product of
the bisection algorithm described in Sec. II B 4. During the bisection, two sets of trajectories with $M_{bis}$ entries are obtained that straddle the edge manifold and stay arbitrarily close to each other for a sufficiently long time. We consider the converged edge trajectory, reported earlier, as our trajectory of reference and exploit the fact that the other off-the-edge trajectories can be considered as trajectories initially linearized around that reference edge trajectory [13].

Every ith pair of trajectories corresponds to a restart of the bisection algorithm, and we measure the initial state space distance between them at each restart, denoted by $\tilde{\Delta}_i$, using the $L^2$-norm. We assume that the two trajectories in each pair diverge exponentially with time, i.e., the evolution of the distance for small times is $\Delta_i(t) = \tilde{\Delta}_i \exp(\gamma_i t)$. Two threshold values for $\Delta_i(t)$, $\delta_{1,i}$, and $\delta_{0,i}$ are specified such that $\delta_{1,i} > \delta_{0,i} > 0$ and $\delta_{1,i} = \varrho_i \delta_{0,i}$. The divergence rate associated with the ith pair of trajectories (called the doubling rate) then reads

$$ \gamma_i = \frac{\ln(\varrho_i)}{\tau_i}, \quad i = 1, \dots, M_{bis}, \qquad (A1) $$

where $\tau_i$ is the time needed for $\Delta_i$ to grow from $\delta_{0,i}$ to $\delta_{1,i}$, and $M_{bis}$ is the number of times the bisection algorithm has been restarted. The time $\tau_i$ is frequently called the doubling time (even if $\varrho_i$ is different from two) [74,75].

In practice, due to the size of the system, trajectories can only be sampled at a number $M_{snap,i}$ of points, from which $\gamma_i$ is determined through a linear least-squares fit between $\{\ln(\Delta_i(t_k)/\tilde{\Delta}_i)\}_{k=0}^{M_{snap,i}}$ and $\{t_k\}_{k=0}^{M_{snap,i}}$. With the tolerances given in Sec. II B 4, we take $\delta_{0,i} = \tilde{\Delta}_i \le 10^{-4}$ (i.e., the tolerance of the bisection) and $\delta_{1,i} = 10^{-3}$ (arbitrarily chosen). While these thresholds are far from infinitesimal, they are small compared to the size of the attractor. Since the lower threshold varies between the $M_{bis}$ pair of trajectories, the factor $\varrho_i$ varies as well. The recorded doubling time $\tau_i$ is taken to be the largest discrete time for which $\Delta_i(\tau_i) < \delta_{1,i}$. Following [75], the FSLE is computed as the time-averaged doubling rate, i.e.,

$$ \lambda = \langle \gamma \rangle_t = \lim_{T \to \infty} \frac{1}{T} \int_0^T \gamma dt = \lim_{M_{bis} \to \infty} \frac{\sum_{i=1}^{M_{bis}} \gamma_i \tau_i}{\sum_{i=1}^{M_{bis}} \tau_i}. \qquad (A2) $$

The resulting value of $\lambda$ is plotted in Fig. 11 using open circles, while the individual values of $\gamma_i$ are plotted with colored dots.

The predictability time from an initial error $\delta_0$ to a given tolerance $\delta_M$ is obtained by integrating the inverse growth rate $\lambda(\delta)$ over all scales [75],

$$ T_p(\delta_0, \delta_M) = \int_{\delta_0}^{\delta_M} \frac{d \ln(\delta')}{\lambda(\delta')}. \qquad (A3) $$

In the infinitesimal limit, $T_p \approx \ln(\delta_M/\delta_0)/\lambda_1$, where $\lambda_1$ is the largest Lyapunov exponent. This limit is also applicable in the

present case with only two defined thresholds, and yields the Ha scaling proposed in Sec. III C.

## APPENDIX B: INITIAL PERTURBATION

A synthetic velocity field is devised in order to trigger turbulence in the system. In the application of the bisection algorithm, the associated perturbation velocity (defined by subtracting the Poiseuille base flow) is rescaled in amplitude until it lies exactly on the edge manifold. This Appendix is devoted to an analytical expression for the chosen velocity field. Although it is known that most three-dimensional velocity fields can be used to converge to the same attractors, we detail the one used in this study for reproducibility purposes.

The perturbation velocity field $v'$ has to be divergence free, which is ensured by considering a potential $\psi$ from which it derives via

$$ v'_{x} = \frac{\partial}{\partial z} \left( \frac{\partial \psi}{\partial y} \right), \qquad (B1a) $$

$$ v'_{y} = \frac{\partial}{\partial z} \left( \frac{\partial \psi}{\partial x} \right), \qquad (B1b) $$

$$ v'_{z} = -2 \frac{\partial^2 \psi}{\partial x \partial y}. \qquad (B1c) $$

The function $\psi$ is assumed separable in the form

$$ \psi(x, y, z) = A f(x, y) g(z). \qquad (B2) $$

Following the work on linear optimal disturbances in the flow at hand [51], we expect three-dimensional structures that reside close to the wall and are inclined against the flow direction to undergo large energy amplification, and hence to be more likely to trigger transition to turbulence. Moreover, the sought disturbance must satisfy the imposed boundary conditions (no slip at the walls and periodicity in $x$ and $z$). To comply with these requirements, we define

$$ f(x, \hat{y}) = \cos[\alpha(x + \sigma\hat{y})] \exp[-(\xi\hat{y})^2](\xi\hat{y})^2, \qquad (B3a) $$

$$ g(z) = \cos(\beta z), \qquad (B3b) $$

where $\hat{y} = 1 - |y|$. The variables $\alpha$ and $\beta$ are, respectively, the streamwise and spanwise wave numbers of the perturbation, respectively, $\sigma$ its planar inclination with respect to the streamwise direction, and $\xi$ a wall-normal decay rate. We set $\alpha = 1$ and $\beta = 3$, which corresponds to near-optimal values at time $t \approx 15-20$ for Re = 5000 and Ha = 10 [51] (this motivated the size of the numerical domain mentioned in Sec. II B). The other parameters are set to $\xi = 4$ and $\sigma = 5$. These values are kept constant throughout the study, regardless of Ha. The perturbation is normalized such that A corresponds to a measure of the disturbance amplitude relative to the laminar centerline velocity $V_c$, i.e., $A = \|v'\|$.

<footnote>[1] P. Manneville, On the transition to turbulence of wall-bounded flows in general, and plane Couette flow in particular, Eur. J. Mech. B-Fluids 49, 345 (2015).</footnote>

<footnote>[2] M. Avila, D. Barkley, and B. Hof, Transition to turbulence in pipe flow, Annu. Rev. Fluid Mech. 55, 575 (2023).</footnote>
[3] B. Eckhardt, Transition to turbulence in shear flows, *Physica A* **504**, 121 (2018).

[4] F. Waleffe, On a self-sustaining process in shear flows, *Phys. Fluids* **9**, 883 (1997).

[5] O. Dauchot and P. Manneville, Local versus global concepts in hydrodynamic stability theory, *J. Phys. II France* **7**, 371 (1997).

[6] D. Barkley, Simplifying the complexity of pipe flow, *Phys. Rev. E* **84**, 016309 (2011).

[7] T. M. Schneider and B. Eckhardt, Edge states intermediate between laminar and turbulent dynamics in pipe flow, *Philos. Trans. R. Soc. A* **367**, 577 (2009).

[8] J. D. Skufca, J. A. Yorke, and B. Eckhardt, Edge of chaos in a parallel shear flow, *Phys. Rev. Lett.* **96**, 174101 (2006).

[9] B. Eckhardt, H. Faisst, A. Schmiegel, and T. M. Schneider, Dynamical systems and the transition to turbulence in linearly stable shear flows, *Philos. Trans. R. Soc. A* **366**, 1297 (2008).

[10] Y. Duguet, A. P. Willis, and R. R. Kerswell, Transition in pipe flow: the saddle structure on the boundary of turbulence, *J. Fluid Mech.* **613**, 255 (2008).

[11] T. Itano and S. Toh, The dynamics of bursting process in wall turbulence, *J. Phys. Soc. Jpn.* **70**, 703 (2001).

[12] S. Toh and T. Itano, A periodic-like solution in channel flow, *J. Fluid Mech.* **481**, 67 (2003).

[13] T. M. Schneider, B. Eckhardt, and J. A. Yorke, Turbulence transition and the edge of chaos in pipe flow, *Phys. Rev. Lett.* **99**, 034502 (2007).

[14] C. S. Paranjape, Y. Duguet, and B. Hof, Oblique stripe solutions of channel flow, *J. Fluid Mech.* **897**, A7 (2020).

[15] M. Beneitez, Y. Duguet, P. Schlatter, and D. S. Henningson, Edge manifold as a Lagrangian coherent structure in a high-dimensional state space, *Phys. Rev. Res.* **2**, 033258 (2020).

[16] A. S. Pikovsky, U. Feudel, and S. P. Kuznetsov, *Strange Non-chaotic Attractors: Dynamics Between Order and Chaos in Quasiperiodically Forced Systems*, World Scientific Series on Nonlinear Science Series A: Volume 56 (World Scientific, Singapore, 2006).

[17] A. Gesla, L. Martin Witkoswki, Y. Duguet, and P. Le Quéré, Subcritical axisymmetric solutions in rotor-stator flows, *arXiv:2312.09874*.

[18] V. I. Arnol'd, *Catastrophe Theory*, 3rd ed. (Springer-Verlag, Berlin, Heidelberg, 2004).

[19] S. Rahmstorf, M. Crucifix, A. Ganopolski, H. Goosse, I. Kamenkovich, R. Knutti, G. Lohmann, R. Marsh, L. A. Mysak, Z. Wang et al., Thermohaline circulation hysteresis: A model intercomparison, *Geophys. Res. Lett.* **32**, L23605 (2005).

[20] M. Scheffer, S. Carpenter, J. A. Foley, C. Folke, and B. Walker, Catastrophic shifts in ecosystems, *Nature (London)* **413**, 591 (2001).

[21] R. M. May, S. A. Levin, and G. Sugihara, Ecology for bankers, *Nature (London)* **451**, 893 (2008).

[22] P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox, Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system, *Philos. Trans. R. Soc. A* **370**, 1166 (2012).

[23] S. Wieczorek, C. Xie, and P. Ashwin, Rate-induced tipping: thresholds, edge states and connecting orbits, *Nonlinearity* **36**, 3238 (2023).

[24] V. Lucarini and T. Bódai, Transitions across melancholia states in a climate model: Reconciling the deterministic and

stochastic points of view, *Phys. Rev. Lett.* **122**, 158701 (2019).

[25] J. M. T. Thompson and J. Sieber, Predicting climate tipping as a noisy bifurcation: A review, *Intl. J. Bifurcat. Chaos* **21**, 399 (2011).

[26] T. M. Lenton, Early warning of climate tipping points, *Nat. Clim. Change* **1**, 201 (2011).

[27] S. R. Carpenter and W. A. Brock, Rising variance: a leading indicator of ecological transition, *Ecol. Lett.* **9**, 311 (2006).

[28] V. Guttal and C. Jayaprakash, Changing skewness: an early warning signal of regime shifts in ecosystems, *Ecol. Lett.* **11**, 450 (2008).

[29] K. Avila, D. Moxey, A. De Lozar, M. Avila, D. Barkley, and B. Hof, The onset of turbulence in pipe flow, *Science* **333**, 192 (2011).

[30] C. S. Paranjape, G. Yalniz, Y. Duguet, N. B. Budanur, and B. Hof, Direct path from turbulence to time-periodic solutions, *Phys. Rev. Lett.* **131**, 034002 (2023).

[31] T. Khapko, T. Kreilos, P. Schlatter, Y. Duguet, B. Eckhardt, and D. S. Henningson, Localized edge states in the asymptotic suction boundary layer, *J. Fluid Mech.* **717**, R6 (2013).

[32] T. Tsukahara, N. Tillmark, and P. Alfredsson, Flow regimes in a plane Couette flow with system rotation, *J. Fluid Mech.* **648**, 5 (2010).

[33] J. W. Miles, On the stability of heterogeneous shear flows, *J. Fluid Mech.* **10**, 496 (1961).

[34] R. Benzi and E. S. Ching, Polymers in fluid flows, *Annu. Rev. Condens. Matter Phys.* **9**, 163 (2018).

[35] P. A. Davidson, *Introduction to Magnetohydrodynamics*, Cambridge Texts in Applied Mathematics (Cambridge University Press, Cambridge, 2017).

[36] B. Knaepen and R. Moreau, Magnetohydrodynamic turbulence at low magnetic Reynolds number, *Annu. Rev. Fluid Mech.* **40**, 25 (2008).

[37] O. Zikanov and A. Thess, Direct numerical simulation of forced MHD turbulence at low magnetic Reynolds number, *J. Fluid Mech.* **358**, 299 (1998).

[38] A. Pothérat and K. Kornet, The decay of wall-bounded MHD turbulence at low *Rm*, *J. Fluid Mech.* **783**, 605 (2015).

[39] A. Pothérat, Quasi-two-dimensional perturbations in duct flows under transverse magnetic field, *Phys. Fluids* **19**, 074104 (2007).

[40] C. J. Camobreco, A. Pothérat, and G. J. Sheard, Subcritical transition to turbulence in quasi-two-dimensional shear flows, *J. Fluid Mech.* **963**, R3 (2023).

[41] H. Branover, *Magnetohydrodynamic Flow in Ducts*, A Halsted Press Book (Wiley, New York; Israel Universities Press, Jerusalem, 1978).

[42] D. Krasnov, O. Zikanov, J. Schumacher, and T. Boeck, Magnetohydrodynamic turbulence in a channel with spanwise magnetic field, *Phys. Fluids* **20**, 095105 (2008).

[43] G. Brethouwer, Y. Duguet, and P. Schlatter, Turbulent–laminar coexistence in wall flows with Coriolis, buoyancy or Lorentz forces, *J. Fluid Mech.* **704**, 137 (2012).

[44] B. Hof, J. Westerweel, T. M. Schneider, and B. Eckhardt, Finite lifetime of turbulence in shear flows, *Nature (London)* **443**, 59 (2006).

[45] M. Avila, F. Mellibovsky, N. Roland, and B. Hof, Streamwise-localized solutions at the onset of turbulence in pipe flow, *Phys. Rev. Lett.* **110**, 224502 (2013).
[46] T. Boeck, Coherent structures, instabilities, and turbulence in interfacial and magnetohydrodynamic flows, Habil. thesis, Department of Mechanical Engineering, Technische Universität Ilmenau, 2010.

[47] V. Bandaru, T. Boeck, D. Krasnov, and J. Schumacher, A hybrid finite difference-boundary element procedure for the simulation of turbulent MHD duct flow at finite magnetic Reynolds number, J. Comput. Phys. **304**, 320 (2016).

[48] S. Zammert and B. Eckhardt, Periodically bursting edge states in plane Poiseuille flow, Fluid Dyn. Res. **46**, 041419 (2014).

[49] S. Zammert and B. Eckhardt, Crisis bifurcations in plane Poiseuille flow, Phys. Rev. E **91**, 041003(R) (2015).

[50] S. A. Orszag, Accurate solution of the Orr-Sommerfeld stability equation, J. Fluid Mech. **50**, 689 (1971).

[51] D. Krasnov, M. Rossi, O. Zikanov, and T. Boeck, Optimal growth and transition to turbulence in channel flow with spanwise magnetic field, J. Fluid Mech. **596**, 73 (2008).

[52] F. H. Busse, Bounds on the transport of mass and momentum by turbulent flow between parallel plates, Z. Angew. Math. Phys. **20**, 1 (1969).

[53] T. Boeck, M. Brynjell-Rahkola, and Y. Duguet, Energy stability of magnetohydrodynamic flow in channels and ducts, J. Fluid Mech. **987**, A33 (2024).

[54] F. Waleffe, Three-dimensional coherent states in plane shear flows, Phys. Rev. Lett. **81**, 4140 (1998).

[55] Y. Duguet, P. Schlatter, and D. S. Henningson, Localized edge states in plane Couette flow, Phys. Fluids **21**, 111701 (2009).

[56] T. Khapko, T. Kreilos, P. Schlatter, Y. Duguet, B. Eckhardt, and D. S. Henningson, Edge states as mediators of bypass transition in boundary-layer flows, J. Fluid Mech. **801**, R2 (2016).

[57] J. Jiménez and P. Moin, The minimal flow unit in near-wall turbulence, J. Fluid Mech. **225**, 213 (1991).

[58] J. M. Hamilton, J. Kim, and F. Waleffe, Regeneration mechanisms of near-wall turbulence structures, J. Fluid Mech. **287**, 317 (1995).

[59] S. A. Neelavara, Y. Duguet, and F. Lusseyran, State space analysis of minimal channel flow, Fluid Dyn. Res. **49**, 035511 (2017).

[60] NEK5000 Version 20.0-dev (2020). Argonne National Laboratory, Illinois, https://nek5000.mcs.anl.gov.

[61] M. Brynjell-Rahkola, A spectral element discretization for quasi-static magnetohydrodynamic flows, Int. J. Numer. Meth. Fl. (2024), doi: 10.1002/fld.5321.

[62] Y. Maday and A. T. Patera, Spectral element methods for the incompressible Navier-Stokes equations, in *State-of-the-Art Surveys in Computational Mechanics*, edited by A. K. Noor and J. T. Oden (ASME, New York, 1989), Chap. 3, pp. 71–143.

[63] J. Kim, P. Moin, and R. Moser, Turbulence statistics in fully developed channel flow at low Reynolds number, J. Fluid Mech. **177**, 133 (1987).

[64] J. Wang, J. Gibson, and F. Waleffe, Lower branch coherent states in shear flows: Transition and control, Phys. Rev. Lett. **98**, 204501 (2007).

[65] T. Kreilos, G. Veble, T. M. Schneider, and B. Eckhardt, Edge states for the turbulence transition in the asymptotic suction boundary layer, J. Fluid Mech. **726**, 100 (2013).

[66] S. Li, X. Shi, L. Wang, and S. Dong, Effects of a spanwise magnetic field on the exact coherent states in a channel flow, Phys. Fluids **36**, 047119 (2024).

[67] S. H. Strogatz, *Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering* (Westview Press, Boulder, CO, 2015).

[68] F. Waleffe, Exact coherent structures in channel flow, J. Fluid Mech. **435**, 93 (2001).

[69] T. M. Schneider, J. F. Gibson, M. Lagha, F. De Lillo, and B. Eckhardt, Laminar-turbulent boundary in plane Couette flow, Phys. Rev. E **78**, 037301 (2008).

[70] M. Beneitez, Y. Duguet, P. Schlatter, and D. S. Henningson, Instability of the optimal edge trajectory in the Blasius boundary layer, J. Fluid Mech. **971**, A42 (2023).

[71] Y. Duguet, A. Monokrousos, L. Brandt, and D. S. Henningson, Minimal transition thresholds in plane Couette flow, Phys. Fluids **25**, 084103 (2013).

[72] E. Aurell, G. Boffetta, A. Crisanti, G. Paladin, and A. Vulpiani, Growth of noninfinitesimal perturbations in turbulence, Phys. Rev. Lett. **77**, 1262 (1996).

[73] E. Aurell, G. Boffetta, A. Crisanti, G. Paladin, and A. Vulpiani, Predictability in the large: An extension of the concept of Lyapunov exponent, J. Phys. A: Math. Gen. **30**, 1 (1997).

[74] A. Pikovsky and A. Politi, *Lyapunov Exponents: A Tool to Explore Complex Dynamics* (Cambridge University Press, Cambridge, 2016).

[75] M. Cencini and A. Vulpiani, Finite size Lyapunov exponent: review on applications, J. Phys. A: Math. Theor. **46**, 254019 (2013).

[76] Y.-C. Lai and T. Tél, *Transient Chaos: Complex Dynamics on Finite Time Scales* (Springer Science & Business Media, New York, 2011), Vol. 173.

[77] D. Lee and H. Choi, Magnetohydrodynamic turbulent flow in a channel at low magnetic Reynolds number, J. Fluid Mech. **439**, 367 (2001).

[78] S. Dong, D. Krasnov, and T. Boeck, Secondary energy growth and turbulence suppression in conducting channel flow with streamwise magnetic field, Phys. Fluids **24**, 074101 (2012).

[79] E. L. Rempel, G. Lesur, and M. R. E. Proctor, Supertransient magnetohydrodynamic turbulence in Keplerian shear flows, Phys. Rev. Lett. **105**, 044501 (2010).

[80] A. Riols, F. Rincon, C. Cossu, G. Lesur, P.-Y. Longaretti, G. I. Ogilvie, and J. Herault, Global bifurcations to subcritical magnetorotational dynamo action in Keplerian shear flow, J. Fluid Mech. **731**, 1 (2013).

[81] F. Carbone, D. Telloni, G. Zank, and L. Sorriso-Valvo, Transition to turbulence in a five-mode Galerkin truncation of two-dimensional magnetohydrodynamics, Phys. Rev. E **104**, 025201 (2021).

[82] B. Hof, A. de Lozar, D. J. Kuik, and J. Westerweel, Repeller or attractor? Selecting the dynamical model for the onset of turbulence in pipe flow, Phys. Rev. Lett. **101**, 214501 (2008).

[83] C. Grebogi, E. Ott, and J. A. Yorke, Fractal basin boundaries, long-lived chaotic transients, and unstable-unstable pair bifurcation, Phys. Rev. Lett. **50**, 935 (1983).

[84] T. Kreilos, B. Eckhardt, and T. M. Schneider, Increasing lifetimes and the growing saddles of shear flow turbulence, Phys. Rev. Lett. **112**, 044503 (2014).

[85] O. Zikanov, D. Krasnov, T. Boeck, A. Thess, and M. Rossi, Laminar-turbulent transition in magnetohydrodynamic duct,
pipe, and channel flows, Appl. Mech. Rev. 66, 030802 (2014).

[86] M. Nagata, Three-dimensional finite-amplitude solutions in plane Couette flow: bifurcation from infinity, J. Fluid Mech. 217, 519 (1990).

[87] Y. Zhuang, B. Yang, V. Mukund, E. Marensi, and B. Hoł, Discontinuous transition to shear flow turbulence, arXiv:2311.11474.

[88] www.gauss-centre.eu.

[89] www.lrz.de.