Swirling electrolyte. Part 1. Hierarchy of catastrophic transitions in steady annular flows

Abstract This study extends our previous work (McCloughan & Suslov, J. Fluid Mech., vol. 887, 2020, A23), where the existence of a saddle-node bifurcation of steady axisymmetric electrolyte flows driven by the Lorentz force in a shallow annular domain was first reported. Here we perform further weakly nonlinear analysis over a wider range of the governing parameters to demonstrate that the previously reported saddle-node bifurcation is a local feature of a global fold catastrophe, which, in turn, is a part of cusp catastrophe occurring as the thickness of the fluid layer increases. The amplitude equation characterising multiple flow solutions in the finite vicinity of catastrophe points is derived. The sensitivity of its coefficients and solutions to the distance from the catastrophe points is assessed demonstrating the robustness of the used analytical procedure. The asymptotic flow solution past the catastrophe point is subsequently obtained and its topology is explored confirming the existence of the secondary circulation in the bulk of flow (two-tori background flow structure). The latter is argued to lead to the appearance of experimentally observable vortices on the fluid surface. The rigorous justification of this conjecture is to be given in Part 2 of the study.


Introduction
Electrolytes are weakly electrically conducting fluids that can be easily prepared in laboratory conditions, say, by dissolving salt in water.When the salt concentration is not high, the electric conductivity of such solutions is relatively low, and so are the electric currents occurring in them when an electric potential difference is applied to electrodes submerged in the fluid.However, when such currents interact with an externally applied magnetic field, a Lorentz force arises that is sufficient to set the bulk of electrolyte in motion without any mechanical intervention (Digilov 2007).This makes such set-ups an attractive playground for designing various non-mechanical fluid mixing systems that have a wide range of practical applications (Moffatt 1991;Pérez-Barrera, Ortiz & Cuevas 2016).A large body of literature also exists on using electromagnetically driven flows in physical modelling of atmospheric phenomena in laboratory conditions (Dovzhenko, Novikov & Obukhov 1979;Dovzhenko, Obukhov & Ponomarev 1981;Dovzhenko, Krymov & Ponomarev 1984;Krymov 1989;Manin 1989;Dolzhanskii, Krymov & Manin 1990;Bondarenko, Gak & Gak 2002;Kenjeres 2011).A more detailed review of relevant applications can be found in our previous publication (Suslov, Pérez-Barrera & Cuevas 2017b).
While originally our research reported here was prompted by such applications, our present work is motivated primarily by our interest in fundamental features of such flows that are surprisingly rich and subtle despite their deceptively simple geometry.Here we consider the experimentally realisable (Digilov 2007;Pérez-Barrera et al. 2015, 2016, 2019) flow occurring in a shallow annular container formed by two vertical coaxial copper cylinders that serve as electrodes.The radial electric current flows through an electrolyte filling the annular gap between them, see figure 1 (also refer to figure 2 in Suslov et al. (2017b) for a more detailed view of an experimental set-up and the description of the experimental procedure).If the set-up is placed in a magnetic field created by a vertically polarised magnet, the azimuthal Lorentz force occurs driving the fluid circumferentially.At first glance this motion appears to be unidirectional, but on closer inspection it turns out that even a relatively slow flow has non-zero radial and vertical velocity components (Suslov et al. 2017b).Due to inertia (the perceived centrifugal force) the bulk of fluid tends to move outwards, that is it acquires a radial velocity component.However, the top-bottom symmetry is broken: since the upper surface of the electrolyte layer is free while the bottom is subject to the no-slip condition, the fluid flows towards an outer cylindrical electrode along the free surface.The deceleration of an outward radial flow caused by the outer wall leads to the increase of the bulk pressure there, which, subsequently, forces the fluid near the bottom to flow inwards.In this way, a meridional circulation with fluid flowing towards the outer annulus wall along the free surface and returning along the solid bottom is created that is combined with the circumferential flow so that the overall fluid motion becomes toroidal.Such a steady toroidal flow was termed as the Type 1 solution in Suslov et al. (2017b).A careful comparison of this flow with other similar counterparts such as Ekman (Ekman 1905;Lilly 1966;Greenspan 1968;Faller 1991), Bödewadt (Bödewadt 1940;MacKerrell 2005) or Stewartson (Stewartson 1957;Schaeffer & Cardin 2005) boundary layers or flows between rotating discs (Moisy et al. 2004) that was undertaken in Suslov et al. (2017b) revealed that none of these candidates are equivalent to the electromagnetically driven flow arising in shallow annular layer.
Quite unexpectedly, a different type of azimuthally uniform flow referred to as Type 2 solution was also discovered in Suslov et al. (2017b).The prominent feature distinguishing it from Type 1 flow is the existence of a secondary counter-rotating toroidal flow structure near the corner formed by the free surface and the outer electrode.Finding this solution came as a surprise because the existence of such a spatially complicated flow structure would generally not be expected in thin layers of viscous fluids.Moreover, the basin of attraction of such a solution has been numerically found to be rather small.Because of that, the convergence of Newton-type iterations to it from a randomly chosen initial guess was literally a 'lucky' coincidence.With one such solution found accidentally, the Type 2 solutions for other sets of governing parameters were obtained by using a careful parametric continuation when the solution obtained for one parameter set was used as an initial guess for iterations performed at slightly varied parametric values.
Such a parametric continuation procedure lead to yet another surprise discovery: for a fixed fluid layer depth, Type 2 flow was found to exist only for a finite range of values of the driving electric current or, equivalently, Reynolds numbers (to be defined in § 2) Re * ≤ Re ≤ Re * * .A further careful investigation undertaken in McCloughan & Suslov (2020a) revealed that the ways in which Type 2 solution ceases to exist for Re < Re * and Re > Re * * are completely different.For small Reynolds numbers, Type 2 solution disappears abruptly via a finite-amplitude 'jump' towards Type 1 flow that continues to exist down to Re = 0 when the driving electric current is switched off and the motion stops.In contrast, in McCloughan & Suslov (2020a) we established that both Type 1 and 2 flows cease to exist at Re * * as a result of a local saddle-node bifurcation when these two solutions 'collide' (that is become topologically indistinguishable) and annihilate each other with no other steady azimuthally uniform solution found numerically beyond Re * * .This posed a question: what happens to a physical flow for large driving currents, that is, for Re > Re * * ?
The experimental observations reported in Pérez- Barrera et al. (2015Barrera et al. ( , 2016)), Suslov et al. (2017b) and Pérez-Barrera et al. (2019) provide several hints.First, the circumferential/toroidal flow persists.Second, experiments confirm the existence of a secondary toroidal structure similar to that of Type 2 solution, see figure 17 in Suslov et al. (2017b).Third, and most importantly, laboratory observations demonstrated the existence of a robust system of anticyclonic vortices appearing on the free surface of the electrolyte layer (Suslov, Pérez-Barrera & Cuevas 2017a).Such vortices form an azimuthally propagating periodic structure and thus cannot be attributed to either Type 1 or 2 circumferentially uniform steady solutions.The observations of the flow's temporal evolution also showed that after the driving current is switched on, the circumferential azimuthally invariant flow always establishes first, and only then the vortices appear.This naturally leads to the hypothesis that free-surface vortices result from an instability of an azimuthally uniform steady basic flow, and that they are effectively perturbations superposed on such a basic flow.This led us to perform a linear stability analysis of Type 1 and 2 solutions reported in McCloughan & Suslov (2020a,b) for cases when either the Reynolds number or the electrolyte layer depth were varied.It revealed that only Type 2 flows can be unstable with respect to various azimuthally periodic perturbations whereas Type 1 flows are always linearly stable.Therefore, we concluded that the necessary condition for the existence of free-surface vortices is the presence of the secondary toroidal flow structure distinguishing Type 2 solution from its Type 1 counterpart.Our further investigation of the character of such an instability discovered yet another surprising result: while the observable vortices look very much like Kelvin-Helmholtz 'cat eyes', Kelvin-Helmholtz instability has nothing to do with their origin.They appear primarily due to Rayleigh's centrifugal instability at the border between two counter-rotating toroidal components of the Type 2 flows.
Since free-surface vortices are observed experimentally even for Re > Re * * (Pérez- Barrera et al. 2015Barrera et al. , 2016)), we hypothesise that they continue to arise on a two-tori background even though the numerical procedure reported in Suslov et al. (2017b) and McCloughan & Suslov (2020a) has failed to find it for large Reynolds numbers.Therefore, our goal here is to make an analytic progress to confirm the existence of a two-tori azimuthally uniform steady flow component beyond the saddle-node bifurcation reported in McCloughan & Suslov (2020a).We achieve this goal by employing an amplitude expansion procedure specifically designed for regimes that exist a finite parametric distance away from a bifurcation point (Pham & Suslov 2018).We demonstrate its robustness by using it to develop an asymptotic solution beyond the saddle-node bifurcation and subsequently using it as an initial guess for an iterative search of a full steady solution.We show that the so-initiated iterations converge to a new steady state, which we term Type 3 flow, that is in excellent qualitative and close quantitative agreement with the developed asymptotic approximation.Part 2 of the study then focuses on the stability characteristics of the newly discovered Type 3 flow, compares them with those of the Type 1 and 2 counterparts reported previously and concludes on the roles both Type 2 and Type 3 azimuthally invariant flows play in the formation of experimentally observed free-surface vortices, which are of ultimate interest in various mixing applications.
The structure of Part 1 is as follows.In § 2 we formulate the mathematical problem and specify the appropriate boundary conditions.In § 3, the weakly nonlinear amplitude expansion procedure is developed that aims to approximate solutions at a post-bifurcation point using quantities computed at pre-bifurcation values of the governing parameters.This requires us to deviate from a standard solvability-condition-based derivation of amplitude equations and rigorously justify the variations we introduce to it in order to achieve our goal.In § 4 we discover that the saddle-node bifurcation that leads to the disappearance of Type 1 and 2 solutions is a local feature of a larger-scale fold catastrophe.This enables us to develop an approximate steady-state azimuthally uniform post-bifurcation solution that is different from both Type 1 and 2 solutions but features a prominent two-tori structure satisfying the necessary condition for the existence of experimentally observable free-surface vortices.In § 5 we make further analytical and computational steps and demonstrate that a newly discovered fold catastrophe itself is a part of yet another higher-dimensional (two-parameter) cusp catastrophe.The conclusions are given in § 7, where we also formulate further questions that are answered by a dedicated investigation that is reported in Part 2 of this work.

Problem formulation and governing equations
We consider an annular layer of an electrolyte of depth h laterally confined by vertical coaxial cylindrical electrodes located at r * = R 1 and r * = R 2 (asterisk denotes  dimensional quantities), see figure 1.The electrically non-conducting stationary bottom of a container is placed on top of a vertically polarised disc magnet of radius R 2 that produces magnetic field B 0 at the reference location marked by the black dot in figure 1.The top of the layer is assumed to be stress-free.
The magnetic field created by a permanent magnet is assumed to be axisymmetric.In general, it has both vertical and radial components that vary with horizontal and vertical coordinates: The magnetic field configuration and details of its computation were discussed in Suslov et al. (2017b), see figure 4 in  § 4 there.When the electric potential difference Δφ 0 is applied between the electrodes, the total current flows mostly radially (see figure 2c) through the layer between them.The interaction of the current and the applied magnetic field creates a Lorentz force F * L = j * × B * , where j * is the current density.This force drives the flow predominantly circumferentially.
Because of a low conductivity of electrolytes (see Pérez-Barrera et al. (2016) and Suslov et al. (2017b) for typical physical characteristics of an electrolyte), hydrodynamic flows arising in the considered set-up practically do not influence the magnetic field of a magnet and the so-called small magnetic Reynolds number approximation (Davidson 2001) to the equations describing fluid motion can be used.Non-dimensional axisymmetric Poisson's equation for the electric potential, the momentum equations and the continuity equation for an incompressible fluid are then written in cylindrical coordinates as (2.4) where φ is the electric potential, (u, v, w) are the velocity components in the radial (r), azimuthal (θ) and vertical (z) directions, respectively, and p is the pressure including the hydrostatic component.In the above equations, the non-dimensional magnetic field created by a disc magnet is B = (B r , 0, B z ).The electric current density components are (2.7a-c) where star denotes dimensional quantities.Here the aspect ratio of the electrolyte layer , the square of Hartmann number Ha 2 characterising electromagnetic effects and Reynolds number Re quantifying viscous effects are where the velocity scale is defined as The typical experimental values of these parameters are ∼ 10 −1 , Ha 2 10 −5 and Re ∼ 10 3 (Suslov et al. 2017b).
The no-slip/no-penetration velocity boundary conditions at the cylindrical electrodes and non-conducting bottom are and Suslov et al. (2017b) and in our current computations).Neglecting the deformation of the free surface (Satijn et al. 2001;Suslov et al. 2017b) the tangential stress-free condition is written as The boundary conditions for the electric potential at the electrode surfaces, the non-conducting bottom and the free surface are respectively.
We also note that a steady-state θ -independent pressure field component satisfies the following Poisson equation (unfortunately, a similar equation given in Suslov et al. (2017b) contained a number of unnoticed typesetting errors; it should be replaced with (2.14), with the Froude number multiplying its right-hand side due to a different pressure scaling) where z are squares of the non-dimensional flow vorticity and velocity and of the applied magnetic field.The pressure must satisfy the following boundary conditions that are consistent with the momentum equations in the vicinity of the physical boundaries and take into account the velocity boundary conditions (2.10) and (2.11).These equations with the specified boundary conditions admit steady θ-independent solutions (referred to as Type 1 and 2 solutions) that have been discussed in detail in Suslov et al. (2017b) and briefly outlined in the introduction here.Linear stability of such solutions with respect to the m-periodic infinitesimal perturbations in the form w (r, z) exp(σ t + imθ) has been subsequently investigated in McCloughan & Suslov (2020a).Note that the solution of Poisson's equation ( 2.2) for electric potential with boundary conditions (2.12a,b) and (2.13a,b) can be written as where This potential induces the electric current j = j + j, where The example of such steady azimuthally invariant electric potential and current distributions is shown in figure 2. It demonstrates that the maximum variation of the electric potential across the flow domain induced by a steady fluid motion field ū = [ū(r, z), v(r, z), w(r, z)] T does not exceed 1.5 × 10 −4 % of the applied electric potential difference so that the associated variation of the Lorentz force remains negligible.The induced current tends to reduce the total (predominantly radial) current, but quantitatively such a reduction is negligible compared with the current applied externally.
In view of this, in what follows we neglect any time-dependent perturbations φ (r, z, t) of the steady-state potential φ + φ that could be induced by an unsteady velocity field perturbations |u (r, z, t)| |ū(r, z)|.This reduces the system of perturbation equations and the corresponding computational cost.A systematic relative error |φ |/|u | introduced by neglecting φ in the linearised equations discussed in the following is of the order of 2 Ha 2 ∼ 10 −6 that we deem acceptable compared with the error of the order of 10 −4 introduced by truncating the asymptotic series developed in § 3. The resulting linearised perturbation equations have been given explicitly in McCloughan & Suslov (2020a) and will not be repeated here.We just mention for the future reference that they can be written in an operator form as where σ = σ R + iσ I is the complex amplification rate and L σ ;Π , A Π and B are matrix-differential and matrix operators, respectively, arising from the linearised perturbation equations and boundary conditions, is the vector of meridional perturbation quantities, Π = {Re, , m} is the set of the governing parameters.Hartmann number is fixed to zero in stability equations.In the subsequent analysis we also take m = 0 so that only θ-independent component of the flow is considered.

Weakly nonlinear formulation
Solutions of linearised equations (2.23) for various values of the governing parameters Π are eigenfunctions corresponding to a discrete eigenvalue spectrum σ i .The eigenvalues are sorted in the order of a decreasing real part and only eigenfunctions corresponding to σ R i > 0, i = 0, 1, 2, . . .(temporally growing perturbations) are deemed to represent physically observable flow patterns.For Type 1 and 2 basic flows, they have been discussed in detail in McCloughan & Suslov (2020a), where it was shown computationally that at most one unstable eigenfunction can be identified for each set of the governing parameters.Furthermore, it was also established in McCloughan & Suslov (2020a) that Type 1 basic flow remained linearly stable whereas Type 2 solutions were unstable with respect to perturbations with wavenumbers m belonging to a finite range.Importantly, it was demonstrated that azimuthally invariant perturbations with m = 0 grew most rapidly.This prompted their weakly nonlinear consideration that was carried out up to the second order in perturbation amplitude, which enabled the authors to conclude that Type 1 and 2 solutions undergo a local saddle-node bifurcation at Re = Re * * .
Here we extend that study to include third-order terms so that we can understand global dynamics of the considered physical system beyond the bifurcation point and discover a hierarchy of catastrophes experienced by the flow as its governing parameters change.The same Chebyshev pseudo-spectral collocation method was used as that described in Suslov et al. (2017b) and McCloughan & Suslov (2020a).Further details of the numerical implementation will be presented in Part 2 of this paper.
Let w represent the steady θ-invariant basic flow solution vector computed at Re 0 Re * * (Re 0 cannot be exactly equal to or greater than Re * * because both the Type 1 and 2 solutions cease to exist and the computational procedure fails to converge there).As discussed in McCloughan & Suslov (2020a) the generalised eigenvalue problem (2.23) derived for m = 0 from equations linearised about w results in a real growth rate σ 0 = σ R 0 , which tends to zero at Re 0 → Re * * .The corresponding eigenfunction is defined up to an arbitrary multiplicative amplitude A. Since the maximum perturbation growth rate is close to zero, the magnitude of the amplitude must be small.Therefore, we rescale it by introducing a formal order parameter ζ : A = ζ Ã, where Ã = O(1).We also introduce the scaled parametric distance from point Re 0 : We look for the asymptotic solution in the vicinity of Re 0 in the form that is informed by the fact that the governing equations have a quadratic nonlinearity (e.g.Pham & Suslov 2018): where the vector of flow quantities is w = [u, v, w, p] T .Subscripts 21, 31 and 22, 32 denote flow deviations from w computed for Re 0 due to the variation of the governing parameter (Re / = Re 0 ) and contributions due to the nonlinearity of the problem, respectively.Substituting (3.2) into the governing equations (2.3)-(2.7a-c)and boundary conditions (2.10) and (2.11) we regain the basic flow equations from the ζ 0 terms and the linearised perturbation equations at the order of ζ 1 , that is, and that L σ 0 ;Π 0 is necessarily a singular operator, which is required to obtain a non-trivial perturbation solution w 0 .Equations arising at the second order of ζ can be written as where 21 f (2) 22 f (2) 22 f (3) 22 0] T are defined in Appendix A. Equation (3.5) is equivalent to a system of two-component equations The operators in the left-hand sides of these equations are not singular and, thus, (3.7a,b) can be solved uniquely.However, since the computational parametric point Π 0 is arbitrarily chosen in the vicinity of Π * * = {Re * * , , 0}, it can approach it asymptotically closely.In this limit, σ 0 → 0 and the operators in the left-hand sides of (3.7a,b) become singular.Therefore, the existence of solutions of these equations with non-zero right-hand sides cannot be guaranteed unless the right-hand-side vectors are put in the range of the operators.The additional degree of freedom required for adjusting the right-hand sides is obtained by noting that (3.4) describing a fast exponential evolution is only valid for an infinitesimal amplitude and it must be modified whenever the amplitude acquires a finite size.Algebraically, this is done by adding terms to the evolution equation (3.4) that have functional forms identical to those of terms appearing in the amplitude expansion (3.2) at any particular order.Therefore, we write which leads to the appearance of additional terms in (3.7a,b): Taking the limit of Π 0 → Π * * or, equivalently, σ 0 → 0, w 0 → w * * , which makes the operators in the left-hand side of (3.9a,b) identical and singular, and considering the inner product of these equations with the adjoint eigenvector w † * * of L 0;Π * * normalised as w † * * , Bw * * = 1 we obtain the classical solvability conditions that define constants K 21 and K 22 (Landau 1944;Stuart 1960;Watson 1960)  This means that all computations necessarily have to be performed at Π 0 / = Π * * so that σ 0 / = 0 and the operators in the left-hand sides of (3.9a,b) remain non-singular, thus, making the choice of K 21 and K 22 ambiguous (given that for non-singular left-hand sides these equations can be solved for any values of the constants).As was explained in Pham & Suslov (2018) such an ambiguity is not an algebraic artefact but rather is a reflection of physical reality: away from the critical point one must chose explicitly the angle of the desired projection of the full solution of the physical problem onto a space spanned by the eigenfunctions of the linearised perturbation problem.In other words, this means that one has to provide an additional condition specifying the main features of the physical problem that are to be retained as closely as possible under the solution projection given by (3.2).As algorithmically proven in Pham & Suslov (2018) and discussed previously in Suslov & Paolucci (1997, 1999) and references therein this additional requirement is formulated as the weighted orthogonality conditions where M is an appropriate weight matrix.Here we choose so that the higher-order terms in (3.2) are weight-orthogonal to the eigenfunctions of the linearised perturbation problem and the low-order truncation of this asymptotic series retains most of the kinetic energy of the full nonlinear solution of the original problem.Subsequently, as shown in Pham & Suslov (2018) equations (3.9a,b) are recast as extended systems that are proven to be non-singular and, thus, can be solved using any standard computational routine for systems of linear equations to simultaneously define (w 21 , K 21 ) and (w 22 , K 22 ).
Proceeding similarly with equations arising at the third order of ζ and taking into account that and where 31 f (2) 32 f (2) 32 f (3) 32 0] T are defined in Appendix A. Then (w 31 , K 31 ) and (w 32 , K 32 ) are found from respectively.Finally, rewriting Ã and δ in terms of A and δ in (3.15) leads to the amplitude equation governing the temporal evolution of the θ-independent flow component near Π 0 .

Determining fold points
The cubic amplitude equation (3.18) with real coefficients can have either one or three real fixed points, see figure 3(a).Of interest here is the parameter set, where the number of such points changes, that is, the fold catastrophe.In terms of parameters defined in the previous sections, the task is to estimate Re * * at which the Type 1 and 2 steady θ-independent solutions cease to exist given (3.18) with coefficients estimated at Re 0 < Re * * .This is required since direct computations at Re * * are impossible, see figure 3(b), where a closeup of a fold is shown.The circles show the closest location to the catastrophe at which computing the two distinct steady θ-independent basic flow solutions is still possible.The phase diagram in figure 3(b) is topologically equivalent to that discussed in detail in McCloughan & Suslov (2020a).Therefore, we conclude that the saddle-node bifurcation analysed there is just a local feature of a larger-scale fold catastrophe that we discovered here by considering a higher-order amplitude expansion of the flow fields.
Setting dA/dt = 0 in (3.18), differentiating it with respect to A, taking into account that at the catastrophe point dδ/dA = 0 and solving for δ we find the location of the fold catastrophe in terms of A: Substituting δ = δ * * from (4.1) to (3.18) leads to a cubic equation for A at the fold catastrophe point: which for σ 0 → 0 has an asymptotically small solution Finally, from (4.1) and (4.3) we obtain the approximate location of the fold catastrophe point 4.2.Parameter sensitivity considerations While (4.4a,b) enables one to determine the fold catastrophe point by computing coefficients of the amplitude equation (3.18) at a parametric point Π 0 some distance away from it, the question arises how sensitive such a prediction is to the choice of Π 0 .As seen from figure 4(a-e) the coefficients of (3.18) do not remain constant as Re deviates from Re * * , which is expected since they depend on the structure of both basic flow and perturbation fields that vary with the strength of an electric current that determines the value of Re.However, the estimations of Re * * remain sufficiently robust, see figure 4( f ).For example, for Re 0 = 1493.88(the largest value of Reynolds number for which the steady θ -independent basic flow solutions can be computed, see black circles in figure 4) we obtain Re * * = 1493.92,whereas for Re 0 = 1470.00we compute Re * * = 1491.49.Therefore, shifting the computational point by 1493.88/1470.00− 1 ≈ 0.016 leads to the Re * * estimation error that is at least an order of magnitude smaller than this shift: 1493.92/1491.49− 1 ≈ 0.0016.Thus, the suggested weakly nonlinear expansion procedure enables one to obtain accurate solutions without approaching the catastrophe point asymptotically closely.This is a practically important conclusion: while the solution at Re 0 = 1470.00can be obtained starting with a relatively crude initial guess, a very careful parametric continuation procedure with at least 50 intermediate computations, each providing an initial guess for the subsequent run, was required to approach the value of Re 0 = 1493.88.This is so because such an approach has to be strictly one-sided with progressively decreasing parametric increments since overshooting Re * * completely breaks the catastrophe search process given that for Re > Re * * solutions cease to exist and Newton-type iterations cannot converge in principle.

Post-fold catastrophe solution
Once the location of the catastrophe point Re * * is determined, the natural question arises: what happens to the physical flow at larger Reynolds numbers?Since Type 1 and 2 steady θ-independent basic flow solutions cease to exist, multiple scenarios may be considered.One of them is that another steady θ-invariant flow state may exist.However, as reported in Suslov et al. (2017b) and McCloughan & Suslov (2020a), attempts to find the third solution directly using Newton-type iterations, similarly to how Type 1 and 2 basic flows were found, fail.This indicates that even if such a solution does exist, its basin of attraction is sufficiently far from solutions found previously.Thus, we resort to the analysis of the third fixed point of (3.18) that belongs to the upper branch of the S-shaped curve in figure 3(a), see the top circle there.It exists in the finite neighbourhood of Π * * for Re ≷ Re * * .We denote the value of the amplitude at that point by A 3 and find an explicit asymptotic expression for it in the limit σ 0 → 0 (or, equivalently, Π → Π * * ): where γ = (Re − Re 0 )/(Re * * − Re 0 ), Re is the Reynolds number of interest, Re 0 is the Reynolds number at which the coefficients of (3.18) are evaluated and Re * * is the Reynolds number of the catastrophe point (approximately given by (4.4a,b)).This enables us to reconstruct the corresponding steady azimuthally invariant component of the flow field using (3.2).To do that, we choose Re 0 = 1493.88< Re * * for which we still can find two distinct steady-state solutions by a standard iterative process.It has been shown in McCloughan & Suslov (2020a) that both Type 1 and 2 basic flows can be used to determine equally good estimations of the catastrophe point Re * * .However, since our goal here is to approximate the third steady-state solution that is topologically 'closer' to the Type 2 basic flow (corresponding to the middle segment of the S-shaped curve in figure 3a), we use the latter as w in (4.4a,b).Subsequently, we compute the coefficients and component fields w , w 21 , w 22 , w 31 and w 32 .The coefficient values that we obtain using N r × N z = 51 × 40 spectral collocation modes (rather than N r × N z = 45 × 31 used in Suslov et al. (2017b) and McCloughan & Suslov (2020a), where the details of numerical approximation have been given) are listed in table 1.Note that the σ 0 , K 21 and K 22 values reported here are somewhat different from those reported previously in McCloughan & Suslov (2020a).This is partially because of a higher spatial resolution used to produce the current results as well as because of the properties of the Matlab (The MathWorks 2020) routine eigs used to compute eigenvectors.It has an internal eigenvector scaling procedure returning results that generally depend on a particular computer architecture so that the same code can return differently scaled eigenvectors when it runs on different computers.Of course, this has no effect on the reported physical results: different scaling of eigenvectors leads to the corresponding change of the computed coefficients and amplitudes so that the product Aw remains invariant and so does the location Re * * ≈ 1493.92 of the catastrophe point.This enables us to estimate the value of A 3 ≈ 1.04 × 10 −1 at a representative value of Re = 1500 > Re * * (γ = 153), see figure 3(a), where we now can reconstruct the post-catastrophe steady θ-independent flow component that cannot be computed directly.Figure 5 shows the converged Type 2 basic flow field.It corresponds to a toroidal flow, which, for lower Reynolds numbers, possess a secondary recirculation in the top outer corner of the meridional cross-section.However, for Re 0 in the close vicinity of Re * * the secondary circulation is completely suppressed so that visually the Type 2 solution depicted in figure 5 is indistinguishable from a typical Type 1 flow.However, such a basic flow is linearly weakly unstable (the growth rate σ 0 is slightly positive, see table 1) with respect to a perturbation that contains a prominent counter-rotating toroidal structure near the outer cylinder with the flow in the rest of the meridional plane almost not affected by it, see figure 6.It is also seen from figure 6(a) that the secondary circulation occurs along with the formation of an azimuthal jet near the corner between the outer cylinder and the free surface.It is directed against the main flow and, consequently, decelerates the overall azimuthal fluid motion there.Therefore, the secondary circulation occurs at the expense of energy taken from the local azimuthal flow.
The second-order perturbation field induced by the problem's nonlinearity enhances both this secondary toroidal structure and the primary meridional circulation of the basic flow and further retards azimuthal flow near the upper right corner, see figure 7. The secondary circulation near the right top corner is also enhanced by the third-order nonlinear perturbation shown in figure 8.However, the second-and third-order perturbations have opposite influences on the main flow structure in the meridional plane: while the second-order perturbations are characterised by clockwise main circulation, its third-order counterpart counteracts it with anti-clockwise motion in the bulk of fluid.The ultimate balance between these two influences leads to the existence of an equilibrium amplitude A 3 defining the fixed point of (3.18), see the top circle in figure 3(a).
Figures 9 and 10 show the variation of the basic flow and the linearised perturbation fields with the parametric shift δ away from the computational point Π 0 , respectively.Even though the contribution of these variations to the overall balance is negligible given a very small parametric shift δ = 1500/1493.88− 1 ≈ 4.1 × 10 −3 considered here, these figures demonstrate a noteworthy but somewhat counterintuitive tendency: the increase of Reynolds number, that is strengthening of electromagnetic flow forcing, leads to the decreasing intensity of both primary and secondary circulations (fluid circulation in figures 5 and 9 and in figures 6 and 10 occurs in opposite directions).From the energy point of view this cannot continue indefinitely and eventually as Re increases the flow is forced to transition to a more energetic state, which provides a physical explanation to the detected fold catastrophe.
Figure 11 depicts the overall asymptotic solution (3.2) at Re = 1500.Despite the original basic flow shown in figure 5 having no visible secondary circulation, the reconstituted solution has it very strongly pronounced.Having constructed this asymptotic solution we now can use it as an initial guess for Newton iterations hoping to find the expected steady state that would differ from the previously reported Type 1 and 2 flows and exist for Re > Re * * .Reporting numerical implementation details of such an attempt are left for Part 2 of this study whereas here we state that indeed the derived asymptotic solution is found to be sufficiently close to the previously undetected steady state that we call the Type 3 flow so that Newton iterations do converge.The resulting field that represents the solution of nonlinear steady θ-independent governing equations (2.2)-(2.6)obtained for Re > Re * * is illustrated in figure 12.It is remarkable that while the approximate asymptotic solution shown in figure 11 differs from the full Type 3 solution illustrated in figure 12 quantitatively, it accurately captures all qualitative features of the full solution including the existence of the secondary circulation and of the submerged jet, where the magnitude of the azimuthal velocity component v exceeds that at the free surface.This demonstrates once again the robustness of the asymptotic procedure developed here.
Even though the analysis here says nothing about the stability of the steady θ-independent Type 3 flow that is different from both Type 1 and 2 basic flows investigated previously, our past studies (McCloughan & Suslov 2020a) demonstrated that the experimentally observable vortices arise on the free surface due to instabilities at the border between the two toroidal flow structures.Therefore, the new θ-independent flow state obtained here has the necessary physical features to serve as a background for to the formation of such vortices in post-catastrophe regimes at Re > Re * * as indeed will be confirmed in Part 2 of this study.

Cusp catastrophe
All numerical results presented so far have been obtained for a fixed value of .To complete the discussion we now turn to investigating what happens as the thickness of an electrolyte layer changes.A significant effort was invested into drawing a parametric existence map of the Type 1 and 2 steady θ-independent basic flow solutions locating fold (saddle-node bifurcation) points by direct computations with parametric continuation in McCloughan & Suslov (2020a).The weakly nonlinear consideration based on the third-order amplitude expansion developed in the present work offers a computationally much cheaper procedure for achieving the same goal and providing a much more insightful view of the arising flow structures.
In figure 13 we present the collection of fixed-point diagrams for (3.18) computed for various non-dimensional electrolyte layer depths .The individual S-shaped amplitude curves have two turning points.The one corresponding to a larger value of Re approximates Re * * , where Type 1 and 2 basic flows annihilate each other in a local saddle-node bifurcation, see the blue solid line in figure 13(a,c).The red solid lines in the same panels connect the second fold points.Qualitatively, they correspond to Re * < Re * * at which the Type 2 solution ceases to exist as found in McCloughan & Suslov (2020a).No conclusive evidence that at Re * the Type 2 solution is replaced with yet another θ-independent solution different from both Type 1 and 2 flows was given in McCloughan & Suslov (2020a) (a numerical steady-state solver used there failed to produce it starting from a random initial guess or using parametric continuation from the Type 2 solution).Yet, the S-shaped curves in figure 13 strongly suggest that such a flow topologically similar to the asymptotic solution shown in figure 11 should exist as indeed will be confirmed in Part 2 of this study.
As the depth of the electrolyte layer increases, the difference between the two turning points of the S-shaped curves approaches zero at ≈ 0.649 (which correspond to the depth of 19.6 mm in the experiments of Pérez-Barrera et al. 2015, 2016), where the fold disappears.Algebraically, this can be seen from ( 4.2) that in the limit σ 0 → 0 factorises to which has three solutions: A 1 = 0 and The behaviour of the discriminant D is shown in figure 14.When it becomes negative, the real amplitude solutions A 2,3 cease to exist.Therefore, D = 0 corresponds to the two-parameter cusp catastrophe.Indeed, the projection of the lines connecting the two fold points in figure 13(a,c) onto the (Re, )-plane (the solid blue and red lines) form a cusp shown in figure 13(b,d) with the fold existing between them.
The location of folds predicted by the analysis of (3.18) remains very close to the values of Re * and Re * * determined from very expensive direct numerical computations completed in McCloughan & Suslov (2020a), compare the dashed and solid lines in figure 13(b,d).
In figure 13(b) the coefficients of (3.18) were computed for Re 0 chosen close to Re * * so that the dashed and solid blue lines are indistinguishable within the plot resolution.The difference between the parametric location of the second fold and Re * (the red lines) is noticeable, yet topologically the two curves are identical and lead to the same cusp location.
Of course, from a practical point of view one would wish to avoid expensive direct computations of the catastrophe points Re * and Re * * replacing them with a much less time-consuming analysis of (3.18).Thus, one generally cannot rely on choosing Π 0 at which its coefficients are to be computed to be close to Π * or Π * * as they may not be known.To demonstrate the robustness of the developed analysis in figure 13(c,d) we present the cusp catastrophe diagrams computed using coefficients of (3.18) evaluated far away from both Re * and Re * * , namely, at a half of parametric distance between them, see the black dashed line in figure 13(d).Despite such a remote computational position, the predicted locations of Re * and Re * * (the solid red and blue lines) remain close to the true values (the dashed lines) with the position of the cusp point determined accurately regardless of where the coefficients were computed.
Finally, we note that once the fold disappears in the cusp catastrophe the topologically different solutions obtained for the past-catastrophe parametric values are expected to morph into each other in a continuous way.This is indeed demonstrated in figure 15.It shows that as the strength of electric current (the magnitude of the Reynolds number) increases near the cusp catastrophe the single-torus flow field very quickly develops the secondary circulation.However, unlike in the fold catastrophe existing in thinner electrolyte layers, this transition occurs in a continuous way.6.Is there a similarity with vortex breakdown or its variants?
Prior to concluding Part 1 of this study we touch upon a question that a reader familiar with other types of swirling flows in a cylindrical geometry may naturally ask: do the azimuthally invariant steady flows discussed in this paper bear any similarity with those described previously in the literature?Our past paper (Suslov et al. 2017b) has addressed it in length with the main conclusion being that despite a superficial similarity the electromagnetically bulk-driven flow considered here has no physically equivalent counterparts among mechanically induced flows.However, consideration in Suslov et al. (2017b) was primarily given to boundary/shear layer-type flows as mentioned in the introduction in this paper.With the full topology of various steady azimuthally invariant flows bulk-driven by the Lorentz force in a shallow annular channel now being established, it is of interest to re-examine this conclusion in a specific context of vortex breakdown, the phenomenon that received much attention starting from the 1960s (Benjamin 1962;Vogel 1968;Escudier 1984).
A classical vortex breakdown occurs in a closed cylinder filled with a fluid the swirling motion of which is driven by the rotation of one of the solid end walls (see, for example, Gelfgat, Bar-Yoseph & Solan (1996) and Blackburn & Lopez (2002) and references therein).The resulting Ekman layer is centrifugally driven towards the side wall next to it.The end-to-end symmetry is broken when the opposite wall of the container is either stationary or rotates with a different speed.The overall flow then consists of the fluid moving along the side wall from the quickly rotating end to the opposite and spiralling back near the axis of the vessel.Superficially, this may appear to be similar to the Type 1 flow in the context of the present study (with the obvious difference that the fluid has to move up along the inner solid cylinder in our case rather than to spiral down freely).
As the rotation speed of the end wall increases the recirculation bubble near the axis of the cylinder forms and this is referred to as the vortex breakdown.Again, this may appear similar to the formation of a secondary circulation as in the Type 2 and 3 solutions reported here (with the difference that in our set-up the secondary circulation occurs near the outer stationary solid wall rather than at the axis).This is where the similarities stop.
The physical reasons for that may be explained following the discussion in Shtern, del Mar Torregrosa & Herrada (2012).These authors argue that the classical vortex breakdown is closely related to the so-called swirl decay, that is, to the fact that the influence of the end wall rotation driving the swirling motion of fluid decreases with the distance cylindrical flow between two counter-rotating discs.Subsequently, the author hypothesised that there should be a third, unstable and thus unaccessible by DNS, state sandwiched between the two stable flows in a parametric space.It was then qualitatively suggested that these three states could form a cusp/fold similar to that we computed and presented in figure 13 here with the Type 1 and 3 being stable and Type 2 unstable solutions (as shown in McCloughan & Suslov (2020a) for the Type 1 and 2 flows and in Part 2 of this manuscript for the Type 3 state).Of course, the physical mechanisms responsible for the existence of multiple steady flow states, which are associated with topologically different boundary layer separation patterns, between two counter-rotating discs are completely different from those acting in a bulk-driven electrolyte layer, where boundary layers either do not exist (at the free surface) or remains safely attached (at the solid bottom).
Finally, we mention a set-up of a free-surface flow in a cylindrical tank with a rotating bottom (e.g.Spohn, Mory & Hoppinger 1993, and references therein).In an idealised formulation the free surface is assumed to be flat, and this approximation has also been used in our computations, see Suslov et al. (2017b) for justification.It might be tempting then to refer to such a similarity and argue the potential equivalence of the bottom-driven cylindrical flow set-up with bulk-driven flows investigated here.Following symmetry arguments Lopez et al. (2004) calculated azimuthally invariant steady flow states by doubling the height of the computational domain and treating the free surface as a symmetry plane of such an extended cylinder with co-rotating end walls.They found that steady flows computed for small values of aspect ratio (Λ = 0.25), which are of relevance here, contained a prominent solid-body-rotation core, the conclusion also confirmed by Yang et al. (2019) for computations performed accounting for a realistic free-surface deformation.No such structure has been found in flows discussed here.Therefore, reinforcing the arguments of Suslov et al. (2017b) we are led to conclude that despite a superficial resemblance none of the well-studied wall-driven swirling flows share their features with flows caused by the action of the bulk Lorentz force.

Conclusions
Despite the geometric simplicity of an annular flow domain and of the electromagnetic forcing driving the primary circumferential fluid motion, the investigated flow reveals a surprisingly rich range of experimentally observable flow patterns that is shown to differ from other swirling flow set-ups such as vortex breakdown in a cylinder with a rotating end wall.The current study has focused primarily on the azimuthally invariant flow component.It continues the computational and analytic effort made in our previous study that discovered a saddle-node bifurcation in which two different steady flow patterns termed Type 1 and 2 in Suslov et al. (2017b) and McCloughan & Suslov (2020a) merge into one and then cease to exist as the Lorentz force driving the flow increases.The current study undertook a higher order (cubic) asymptotic analysis based on perturbation amplitude (A) expansion.It enabled us to strengthen the hypothesis proposed in McCloughan & Suslov (2020a) that another azimuthally invariant flow pattern exists for sufficiently large values of the radial current (and, thus, the driving Lorentz force parameterised by the Reynolds number Re).Such a pattern is different from Type 1 and 2 flows but has a secondary circulation that is qualitatively similar to that of Type 2 flow reported previously.In the amplitude-forcing (A, Re) phase space, along with Type 1 and 2 flows such a pattern forms a classical fold catastrophe.The previously reported saddle-node bifurcation is found to be its local feature.

Figure 2 .
Figure 2. Steady azimuthally invariant distributions of the externally applied (a) and fluid-motion-induced (b) components of the electric potential given by (2.18) and the current components given by (2.21) (c) and (2.22) (b), respectively, for the Type 2 flow at (Re, ) = (1493.88,0.248) and 2 Ha 2 = 1.59 × 10 −6 .The colour in the bottom right panel represents the value of the azimuthal electric current component induced by fluid motion.
.10a,b) However, implementing the above in practice in the current problem is impossible because neither steady-state basic flow solutions, nor their linearised perturbation, nor the adjoint eigenvector can be found at Π * * (see also the discussion in McCloughan & Suslov 2020a).

Figure 3 .
Figure 3. (a) Fold catastrophe and (b) saddle-node bifurcation phase diagrams computed for the Type 2 basic flow at Π 0 = (Re 0 , , Ha, m) =(1493.88,0.248, 5.08 × 10 −3 , 0).The curves show the location of fixed points of (3.18) and the vectors correspond to the value of dA/dt.The circles indicate values of the fixed-point amplitudes at δ = 0 corresponding to the computational parametric point Π 0 and the square marks the position of the predicted fold catastrophe point that is not accessible by direct computations.

Figure 4 .
Figure 4. Coefficients of (3.18) evaluated for Type 2 basic flow and estimations of Re * * as functions of Re 0 for ( , Ha) = (0.248, 5.08 × 10 −3 ).Black circles show the values computed at the parametric point that is closest to the catastrophe point, where convergence of iterations can still be achieved.

Figures 5 -
Figures 5-8 depict the flow fields for each component of the asymptotic solution and figure 11 shows the approximated full flow field given by (3.2).Panels (a) and (b) in each figure illustrate the meridional velocity and the azimuthal vorticity of the flow field, respectively.Figure5shows the converged Type 2 basic flow field.It corresponds to a toroidal flow, which, for lower Reynolds numbers, possess a secondary recirculation in the top outer corner of the meridional cross-section.However, for Re 0 in the close vicinity of Re * * the secondary circulation is completely suppressed so that visually the Type 2 solution depicted in figure5is indistinguishable from a typical Type 1 flow.However, such a basic flow is linearly weakly unstable (the growth rate σ 0 is slightly positive, see table1) with respect to a perturbation that contains a prominent counter-rotating toroidal structure near the outer cylinder with the flow in the rest of the meridional plane almost not affected by

Figure 6 .
Figure 6.Meridional θ-independent linearised perturbation velocity (a) and azimuthal vorticity (b) fields for the same parameters as in figure 5.

Figure 7 .
Figure 7. Second-order meridional θ -independent perturbation velocity (a) and azimuthal vorticity (b) fields for the same parameter values as in figure 5.

Figure 8 .
Figure 8. Third-order meridional θ-independent perturbation velocity (a) and azimuthal vorticity (b) fields for the same parameter values as in figure 5.

Figure 9 .
Figure 9. Meridional variation of θ -independent Type 2 basic flow velocity (a) and azimuthal vorticity (b) fields near the computational point Π 0 for the same parameters as in figure 5.

Figure 10 .
Figure 10.Third-order variation of meridional θ-independent perturbation velocity (a) and azimuthal vorticity (b) fields for the same parameters as in figure 5.

Figure 11 .
Figure 11.Approximation of the steady θ -independent component of the flow for Re = 1500 > Re * * : meridional velocity (a) and azimuthal vorticity (b) fields for the same values of and Ha as in figure 5.

Figure 12 .
Figure 12.Full numerical steady θ -independent (Type 3) flow solution: meridional velocity (a) and azimuthal vorticity (b) fields for the same parameters as in figure 11.

Figure 13 .Figure 14 .
Figure 13.Cusp catastrophe diagrams computed for (a,b) Re 0 = Re * + 0.99(Re * * − Re * ) and (c,d) Re 0 = (Re * + Re * * )/2 in the (Re, ) parametric space.The red and blue dashed lines in panels (b) and (d) correspond to the values Re * and Re * * , respectively, reported in McCloughan & Suslov (2020a), the black dashed lines show the values Re 0 at which coefficients of (3.18) have been computed in each panel.The solid lines illustrate the positions of folds.