2. Formulation of problem

Consider the flow of an incompressible viscous electrically conducting liquid with density ${\it\rho}$ , kinematic viscosity ${\it\nu}$ and electrical conductivity ${\it\sigma}$ driven by a constant gradient of pressure $p$ in a channel of width $2h$ between two parallel walls in the presence of a transverse homogeneous magnetic field $\boldsymbol{B}$ . The velocity distribution of the flow $\boldsymbol{v}(\boldsymbol{r},t)$ is governed by the Navier–Stokes equation

(2.1) $$\begin{eqnarray}\partial _{t}\boldsymbol{v}+(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{v}=-{\it\rho}^{-1}\boldsymbol{{\rm\nabla}}p+{\it\nu}{\rm\nabla}^{2}\boldsymbol{v}+{\it\rho}^{-1}\boldsymbol{f},\end{eqnarray}$$

where $\boldsymbol{f}=\boldsymbol{j}\times \boldsymbol{B}$ is the electromagnetic body force containing the induced electric current $\boldsymbol{j}$ , which is governed by Ohm’s law for a moving medium

(2.2) $$\begin{eqnarray}\boldsymbol{j}={\it\sigma}(\boldsymbol{E}+\boldsymbol{v}\times \boldsymbol{B}),\end{eqnarray}$$

where $\boldsymbol{E}$ is the electric field in the stationary frame of reference. The flow is assumed to be sufficiently slow that the induced magnetic field is negligible relative to the imposed one. This supposes a small magnetic Reynolds number $\mathit{Rm}={\it\mu}_{0}{\it\sigma}v_{0}\,L\ll 1$ , where ${\it\mu}_{0}$ is the permeability of vacuum and $v_{0}$ and $L$ are the characteristic velocity and length scale of the flow, respectively. Using the thickness of the Hartmann layer $L={\it\delta}=B^{-1}\sqrt{{\it\rho}{\it\nu}/{\it\sigma}}$ as the relevant length scale in strong magnetic field, we have $\mathit{Rm}\sim R\mathit{Pm}$ , where $R=\mathit{Re}/\mathit{Ha}$ is a local Reynolds number based on ${\it\delta}$ and $\mathit{Pm}$ is the magnetic Prandtl number. Taking into account that for typical liquid metals $\mathit{Pm}\sim 10^{-5}$ , the constraint on $\mathit{Rm}$ translates into $R\ll \mathit{Pm}^{-1}\sim 10^{5}$ . In addition, we assume that the characteristic time of velocity variation ${\it\tau}\sim L/v_{0}$ is much longer than the magnetic diffusion time ${\it\tau}_{m}={\it\mu}_{0}{\it\sigma}L^{2}$ . This is also satisfied by the constraint on $\mathit{Rm}$ and thus allows us to use the quasi-stationary approximation leading to $\boldsymbol{E}=-\boldsymbol{{\rm\nabla}}{\it\phi}$ , where ${\it\phi}$ is the electrostatic potential (Roberts Reference Roberts1967).

The velocity and current satisfy mass and charge conservation $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{j}=0$ . Applying the latter to Ohm’s law (2.2) yields

(2.3) $$\begin{eqnarray}{\rm\nabla}^{2}{\it\phi}=\boldsymbol{B}\boldsymbol{\cdot }{\bf\omega},\end{eqnarray}$$

where ${\bf\omega}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{v}$ is vorticity. At the channel walls $S$ , the normal $(n)$ and tangential $({\it\tau})$ velocity components satisfy the impermeability and no-slip boundary conditions $v_{n}|_{s}=0$ and $v_{{\it\tau}}|_{s}=0$ .

We employ right-handed Cartesian coordinates with the origin set at the mid-height of the channel, and the $x$ - and the $z$ -axes directed, respectively, against the applied pressure gradient $\boldsymbol{{\rm\nabla}}p_{0}=P\boldsymbol{e}_{x}$ and along the magnetic field $\boldsymbol{B}=B\boldsymbol{e}_{z}$ so that the channel walls are located at $z=\pm h$ , as shown in figure 1, and the velocity is defined as $\boldsymbol{v}=(u,v,w)$ . Subsequently, all variables are non-dimensionalized by using $h$ , $h^{2}/{\it\nu}$ and $Bh{\it\nu}$ as the length, time and electric potential scales, respectively. The velocity is scaled by the viscous diffusion speed ${\it\nu}/h$ , which we employ as the characteristic velocity instead of the commonly used centreline velocity.

Figure 1.

Sketch of the problem showing velocity profiles of Poiseuille and Hartmann flows.

The problem admits a rectilinear base flow

(2.4) $$\begin{eqnarray}\boldsymbol{v}_{0}(z)={\bar{u}}_{0}(z)\boldsymbol{e}_{x}=\mathit{Re}\,{\bar{u}}(z)\boldsymbol{e}_{x}\end{eqnarray}$$

for which (2.1) reduces to

(2.5) $$\begin{eqnarray}{\bar{u}}^{\prime \prime }-\mathit{Ha}^{2}{\bar{u}}=\bar{P},\end{eqnarray}$$

where $\mathit{Re}=Uh/{\it\nu}$ is the Reynolds number based on the centreline velocity $U$ of unperturbed flow, $\mathit{Ha}=Bh\sqrt{{\it\sigma}/{\it\rho}{\it\nu}}$ is the Hartmann number, and $\bar{P}$ is a dimensionless coefficient which relates $U$ with the applied pressure gradient $P=\bar{P}U{\it\nu}{\it\rho}/h^{2}$ by satisfying the normalization condition ${\bar{u}}(0)=1$ . This Reynolds number is convenient for characterizing a flow driven by a fixed pressure gradient. Besides this so-called pressure Reynolds number, one can also use a Reynolds number based on the flow rate, $\mathit{Re}_{q}$ , which is more suitable for the case of fixed flow rate and will be introduced in § 5. Note that the use of either of these two Reynolds numbers is a matter of convenience as long as one remains a time-independent function of the other. Equation (2.5) defines the well-known Hartmann flow profile

(2.6) $$\begin{eqnarray}{\bar{u}}(z)=\frac{\cosh (\mathit{Ha})-\cosh (z\mathit{Ha})}{\cosh (\mathit{Ha})-1}\end{eqnarray}$$

with $\bar{P}=-(\mathit{Ha}^{2}\cosh (\mathit{Ha})/\!\cosh (\mathit{Ha})-1)$ , which corresponds to a channel with perfectly conducting walls. In a weak magnetic field $(\mathit{Ha}\ll 1)$ , the Hartmann flow reduces to the classic plane Poiseuille flow ${\bar{u}}(z)=1-z^{2}$ .

Note that electrical conductivity of the walls affects only the relationship between the applied pressure gradient and the centreline velocity but not the profile of the Hartmann flow. Thus, the stability of Hartmann flow with respect to transverse 2D disturbances, which are not affected by the conductivity of the walls, is determined entirely by the centreline velocity. This is the case considered in the present study.

3. Theoretical background

3.1. Linear stability of the base flow

The two-dimensional travelling waves considered in this study are expected to emerge as the result of linear instability of the Hartmann flow (2.4) with respect to infinitesimal perturbations $\boldsymbol{v}_{1}(\boldsymbol{x},t)$ . Owing to the invariance of the base flow in both $t$ and $\boldsymbol{x}=(x,y)$ , perturbations are sought as Fourier modes

(3.1) $$\begin{eqnarray}\boldsymbol{v}_{1}(\boldsymbol{r},t)=\hat{\boldsymbol{v}}(z)\text{e}^{{\it\lambda}t+\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}+\text{c.c.}\end{eqnarray}$$

defined by the complex amplitude distribution $\hat{\boldsymbol{v}}(z)$ , temporal growth rate ${\it\lambda}$ and the wave vector $\boldsymbol{k}=({\it\alpha},{\it\beta})$ . The incompressibility constraint, which takes the form $\boldsymbol{D}_{k}\boldsymbol{\cdot }\hat{\boldsymbol{v}}=0$ , where $\boldsymbol{D}_{k}\equiv \boldsymbol{e}_{z}(\text{d}/\text{d}z)+\text{i}\boldsymbol{k}$ is a spectral counterpart of the nabla operator, is satisfied by expressing the component of the velocity perturbation in the direction of the wave vector as ${\hat{u}}_{\Vert }=\boldsymbol{e}_{\Vert }\boldsymbol{\cdot }\hat{\boldsymbol{v}}=\text{i}k^{-1}{\hat{w}}^{\prime }$ , where $\boldsymbol{e}_{\Vert }=\boldsymbol{k}/k$ and $k=|\boldsymbol{k}|$ . Taking the curl of the linearized counterpart of (2.1) to eliminate the pressure gradient and then projecting it onto $\boldsymbol{e}_{z}\times \boldsymbol{e}_{\Vert }$ , after some transformations we obtain a modified Orr–Sommerfeld-type equation which includes a magnetic term

(3.2) $$\begin{eqnarray}{\it\lambda}\boldsymbol{D}_{k}^{2}{\hat{w}}=\left[\boldsymbol{D}_{k}^{4}-\mathit{Ha}^{2}(\boldsymbol{e}_{z}\boldsymbol{\cdot }\boldsymbol{D}_{k})^{2}+\text{i}k\mathit{Re}({\bar{u}}^{\prime \prime }-{\bar{u}}\boldsymbol{D}_{k}^{2})\right]{\hat{w}}.\end{eqnarray}$$

The no-slip and impermeability boundary conditions require

(3.3) $$\begin{eqnarray}{\hat{w}}={\hat{w}}^{\prime }=0\quad \text{at}\ z=\pm 1.\end{eqnarray}$$

The equation above is written in a non-standard form corresponding to our choice of the characteristic velocity. Note that the Reynolds number appears in this equation as a factor in the convective term rather than its reciprocal in the viscous term as in the standard form. As a result, the growth rate ${\it\lambda}$ differs by a factor $\mathit{Re}$ from its standard definition.

Since (3.2) like its non-magnetic counterpart admits Squire’s transformation, in the following we consider only two-dimensional perturbations $(k={\it\alpha})$ , which are the most unstable (Lock Reference Lock1955). The problem is solved numerically using the Chebyshev collocation method which is described in detail by Hagan & Priede (Reference Hagan and Priede2013a ).

3.2. 2D nonlinear travelling waves

Two-dimensional travelling waves emerge as follows. First, the neutrally stable mode (3.1) with a purely real frequency ${\it\omega}=-\text{i}{\it\lambda}$ interacting with itself through the quadratically nonlinear term in (2.1) produces a steady streamwise-invariant perturbation of the mean flow as well as a second harmonic ${\sim}\text{e}^{2\text{i}({\it\omega}t+{\it\alpha}x)}$ . Further nonlinear interactions produce higher harmonics, which similarly to the fundamental and second harmonic travel with the same phase speed $c=-{\it\omega}/{\it\alpha}$ . Thus, the solution can be sought in the form

(3.4) $$\begin{eqnarray}\boldsymbol{v}(\boldsymbol{r},t)=\mathop{\sum }\limits_{n=-\infty }^{\infty }E^{n}\hat{\boldsymbol{v}}_{n}(z),\end{eqnarray}$$

where $E=\text{e}^{\text{i}({\it\omega}t+{\it\alpha}x)}$ contains ${\it\omega}$ , which needs to determined together $\hat{\boldsymbol{v}}_{n}$ by solving a nonlinear eigenvalue problem (Huerre & Rossi Reference Huerre, Rossi, Godreche and Manneville1998). The reality of solution requires $\hat{\boldsymbol{v}}_{-n}=\hat{\boldsymbol{v}}_{n}^{\ast }$ , where the asterisk stands for the complex conjugate. The incompressibility constraint applied to the $n$ th velocity harmonic results in $\boldsymbol{D}_{{\it\alpha}_{n}}\boldsymbol{\cdot }\hat{\boldsymbol{v}}_{n}=0$ , where $\boldsymbol{D}_{{\it\alpha}_{n}}\equiv \boldsymbol{e}_{z}(\text{d}/\text{d}z)+\text{i}\boldsymbol{e}_{x}{\it\alpha}_{n}$ with ${\it\alpha}_{n}={\it\alpha}n$ stands for the spectral counterpart of the nabla operator. This constraint can be satisfied by expressing the streamwise velocity component

(3.5) $$\begin{eqnarray}{\hat{u}}_{n}=\boldsymbol{e}_{x}\boldsymbol{\cdot }\hat{\boldsymbol{v}}_{n}=\text{i}{\it\alpha}_{n}^{-1}{\hat{w}}_{n}^{\prime },\end{eqnarray}$$

in terms of the transverse component ${\hat{w}}_{n}=\boldsymbol{e}_{z}\boldsymbol{\cdot }\hat{\boldsymbol{v}}_{n}$ , which we employ instead of the commonly used stream function. Henceforth, the prime is used as a shorthand for $\text{d}/\text{d}z$ . Note that (3.5) is not applicable to the zeroth harmonic, for which it yields ${\hat{w}}_{0}\equiv 0$ . Thus, ${\hat{u}}_{0}$ needs to be considered separately in this velocity-based formulation.

Taking the curl of (2.1) to eliminate the pressure gradient and then projecting it onto $\boldsymbol{e}_{y}$ , we obtain

(3.6) $$\begin{eqnarray}[\boldsymbol{D}_{{\it\alpha}_{n}}^{2}-\text{i}{\it\omega}n]\hat{{\it\zeta}}_{n}-\mathit{Ha}^{2}{\hat{u}}_{n}^{\prime }={\hat{h}}_{n},\end{eqnarray}$$

where

(3.7) $$\begin{eqnarray}\hat{{\it\zeta}}_{n}=\boldsymbol{e}_{y}\boldsymbol{\cdot }\boldsymbol{D}_{{\it\alpha}_{n}}\times \hat{\boldsymbol{v}}_{n}=\left\{\begin{array}{@{}ll@{}}\text{i}{\it\alpha}_{n}^{-1}\boldsymbol{D}_{{\it\alpha}_{n}}^{2}{\hat{w}}_{n},\quad & n\neq 0,\\ {\hat{u}}_{0}^{\prime },\quad & n=0,\end{array}\right.\end{eqnarray}$$

and

(3.8) $$\begin{eqnarray}{\hat{h}}_{n}=\displaystyle \mathop{\sum }\limits_{m}\hat{\boldsymbol{v}}_{n-m}\boldsymbol{\cdot }\boldsymbol{D}_{{\it\alpha}_{m}}\hat{{\it\zeta}}_{m}\end{eqnarray}$$

are the $y$ -components of the $n$ th harmonic of the vorticity ${\bf\zeta}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{v}$ and that of the curl of the nonlinear term $\boldsymbol{h}=\boldsymbol{{\rm\nabla}}\times (\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{v}$ . Henceforth, the omitted summation limits are assumed to be infinite. Separating the terms involving ${\hat{u}}_{0}$ in (3.8), it can be rewritten as ${\hat{h}}_{n}=\text{i}{\it\alpha}_{n}^{-1}({\hat{h}}_{n}^{w}+{\hat{h}}_{n}^{u})$ , where

(3.9) $$\begin{eqnarray}\displaystyle & {\hat{h}}_{n}^{w}=n\displaystyle \mathop{\sum }\limits_{m\neq 0}m^{-1}({\hat{w}}_{n-m}\boldsymbol{D}_{{\it\alpha}_{m}}^{2}{\hat{w}}_{m}^{\prime }-{\hat{w}}_{m}^{\prime }\boldsymbol{D}_{{\it\alpha}_{n-m}}^{2}{\hat{w}}_{n-m}), & \displaystyle\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle & {\hat{h}}_{n}^{u}=\text{i}{\it\alpha}_{n}[{\hat{u}}_{0}-{\hat{u}}_{0}^{\prime \prime }\boldsymbol{D}_{{\it\alpha}_{n}}^{2}]{\hat{w}}_{n}\equiv \mathscr{N}_{n}({\hat{u}}_{0}){\hat{w}}_{n}. & \displaystyle\end{eqnarray}$$

Eventually, using the expressions above, (3.6) can be written as

(3.11) $$\begin{eqnarray}\mathscr{L}_{n}(\text{i}{\it\omega},{\hat{u}}_{0}){\hat{w}}_{n}={\hat{h}}_{n}^{w},\end{eqnarray}$$

with the operator

(3.12) $$\begin{eqnarray}\mathscr{L}_{n}(\text{i}{\it\omega},{\hat{u}}_{0})=[\boldsymbol{D}_{{\it\alpha}_{n}}^{2}-\text{i}{\it\omega}n]\boldsymbol{D}_{{\it\alpha}_{n}}^{2}-\mathit{Ha}^{2}(\boldsymbol{e}_{z}\boldsymbol{\cdot }\boldsymbol{D}_{{\it\alpha}_{n}})^{2}-\mathscr{N}_{n}({\hat{u}}_{0}).\end{eqnarray}$$

This equation governs all harmonics except the zeroth one, for which, in accordance with the incompressibility constraint (3.5), it implies ${\hat{w}}_{0}\equiv 0$ . The zeroth velocity harmonic, which has only the streamwise component ${\hat{u}}_{0}$ , is governed directly by the $x$ -component of the Navier–Stokes equation (2.1):

(3.13) $$\begin{eqnarray}{\hat{u}}_{0}^{\prime \prime }-\mathit{Ha}^{2}{\hat{u}}_{0}=\hat{P}_{0}+{\hat{g}}_{0},\end{eqnarray}$$

where $\hat{P}_{0}=\bar{P}\mathit{Re}$ is a dimensionless mean pressure gradient and

(3.14) $$\begin{eqnarray}{\hat{g}}_{0}=\text{i}\mathop{\sum }\limits_{m\neq 0}{\it\alpha}_{m}^{-1}({\hat{w}}_{m}^{\ast }{\hat{w}}_{m}^{\prime })^{\prime }\end{eqnarray}$$

is the $x$ -component of the zeroth harmonic of the nonlinear term $\boldsymbol{g}=(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{v}$ . Velocity harmonics are subject to the usual no-slip and impermeability boundary conditions

(3.15) $$\begin{eqnarray}{\hat{w}}_{n}={\hat{w}}_{n}^{\prime }={\hat{u}}_{0}=0\quad \text{at}\ z=\pm 1.\end{eqnarray}$$

3.3. Linear stability of 2D travelling waves

Weakly subcritical equilibrium states, which exist in this case (Hagan & Priede Reference Hagan and Priede2013b ), are unconditionally unstable (Schmid & Henningson Reference Schmid and Henningson2001). This is because the growth rate of subcritical disturbances increases with their amplitude (Hagan & Priede Reference Hagan and Priede2013b ). Thus, a disturbance with an amplitude slightly lower or higher than the equilibrium one will respectively decay or grow, so diverging from the equilibrium state. The stability of strongly subcritical equilibrium states is not obvious. Orszag & Patera (Reference Orszag and Patera1983) originally suggested that the subcritical equilibrium state appearing at the linear stability threshold remains linearly unstable down the lowest possible Reynolds number admitting such states. At this limiting Reynolds number, which is the main concern of the present study, the unstable subcritical state disappears by merging with another equilibrium state of a higher amplitude. The latter was thought by Orszag & Patera (Reference Orszag and Patera1983) to be linearly stable as in the saddle–node bifurcation. This simple picture was amended by Pugh & Saffman (Reference Pugh and Saffman1988) who showed that this is the case when the flow is driven by a fixed flow rate but not by a fixed pressure gradient. Although travelling-wave states at fixed flow rate are physically equivalent to those at fixed pressure gradient (Soibelman & Meiron Reference Soibelman and Meiron1991), it is not the case in general when the flow rate and the mean pressure gradient depend not only on each other but also on time. First of all, the distinction between flows driven by fixed pressure gradient and fixed flow rate becomes important when the stability of travelling-wave states is considered.

Linear stability of the travelling-wave states, which in contrast to the rectilinear base state are periodic rather than invariant in both the time and the streamwise direction, is described by Floquet theory (Bender & Orsag Reference Bender and Orsag1978; Herbert Reference Herbert1988) according to which a small-amplitude velocity disturbance can be sought similarly to (3.4) as

(3.16) $$\begin{eqnarray}\boldsymbol{v}_{1}(\boldsymbol{r},t)=\text{e}^{\tilde{{\it\lambda}}t}\mathop{\sum }\limits_{n=-\infty }^{\infty }E^{n+{\it\epsilon}}\tilde{\boldsymbol{v}}_{n}(z)+\text{c.c.},\end{eqnarray}$$

where $\tilde{{\it\lambda}}$ is generally a complex growth rate and ${\it\epsilon}$ is a real detuning parameter which defines the sideband wavenumber $\tilde{{\it\alpha}}={\it\epsilon}{\it\alpha}$ (Huerre & Rossi Reference Huerre, Rossi, Godreche and Manneville1998), also called the subharmonic wavenumber by Soibelman & Meiron (Reference Soibelman and Meiron1991). The case ${\it\epsilon}=0$ corresponds to the fundamental mode, which is also called superharmonic (Pugh & Saffman Reference Pugh and Saffman1988; Soibelman & Meiron Reference Soibelman and Meiron1991), whereas ${\it\epsilon}=\pm \frac{1}{2}$ correspond to the so-called subharmonic mode. The disturbances with other values of ${\it\epsilon}$ are referred to as combination or detuned modes. Note that (3.16) with ${\it\epsilon}\pm 1$ is equivalent to the mode with ${\it\epsilon}$ and the index $n$ shifted by one. In addition, as seen from (3.16) the modes with $\pm {\it\epsilon}$ are complex conjugate. Thus, it suffices to consider $0\leqslant {\it\epsilon}<1$ or any other interval of ${\it\epsilon}$ of length one. Separating the two modes with opposite $z$ -symmetries, as discussed in the next section, this interval can be reduced by half.

Adding disturbance (3.16) to the travelling-wave base state (3.4), we obtain a linearized counterpart of (3.6) for the transverse velocity harmonic $\tilde{w}_{n}$ :

(3.17) $$\begin{eqnarray}\left[\boldsymbol{D}_{\tilde{{\it\alpha}}_{n}}^{2}-\tilde{{\it\lambda}}-\text{i}{\it\omega}n\right]\boldsymbol{D}_{\tilde{{\it\alpha}}_{n}}^{2}\tilde{w}_{n}-\mathit{Ha}^{2}\tilde{w}_{n}^{\prime \prime }=-\text{i}\tilde{{\it\alpha}}_{n}\mathop{\sum }\limits_{m}\boldsymbol{D}_{\tilde{{\it\alpha}}_{n}}\boldsymbol{\cdot }(\hat{\boldsymbol{v}}_{n-m}\tilde{{\it\zeta}}_{m}+\tilde{\boldsymbol{v}}_{m}\hat{{\it\zeta}}_{n-m}),\end{eqnarray}$$

where $\tilde{{\it\alpha}}_{n}={\it\alpha}_{n}+\tilde{{\it\alpha}}=(n+{\it\epsilon}){\it\alpha}$ is a detuned wavenumber for the $n$ th harmonic, ${\tilde{u}}_{n}=\boldsymbol{e}_{x}\boldsymbol{\cdot }\tilde{\boldsymbol{v}}_{n}=\text{i}\tilde{{\it\alpha}}_{n}^{-1}\tilde{w}_{n}^{\prime }$ is the streamwise velocity and $\tilde{{\it\zeta}}_{n}=\text{i}\tilde{{\it\alpha}}_{n}^{-1}\boldsymbol{D}_{\tilde{{\it\alpha}}_{n}}^{2}\tilde{w}_{n}$ is the spanwise vorticity of the disturbance. The boundary conditions, as usual, are $\tilde{w}_{n}(1)=\tilde{w}_{n}^{\prime }(1)=0$ .

Cutting the series (3.16) off at $n=\pm N$ , we obtain a linear eigenvalue problem represented by $2N+1$ complex equations (3.17) for the eigenvector consisting of the same number of harmonics $\tilde{w}_{n}$ with the eigenvalue $\tilde{{\it\lambda}}$ , which depends on the subharmonic wavenumber $\tilde{{\it\alpha}}$ . This eigenvalue problem is solved in the same way as that for the linear stability of the rectilinear base flow. To avoid the division by zero in the expressions for ${\tilde{u}}_{0}$ and $\tilde{{\it\zeta}}_{0}$ above, which occurs for $\tilde{{\it\alpha}}=0$ , we use the substitution $\tilde{w}_{n}=-\text{i}\tilde{{\it\alpha}}_{n}\tilde{{\it\psi}}_{n}$ , where $\tilde{{\it\psi}}_{n}$ is the stream function. This makes superharmonic disturbances treatable in the same way as detuned ones. For the superharmonic disturbances corresponding to $\tilde{{\it\alpha}}\rightarrow 0$ , the eigenvalue problem becomes self-adjoint with $\tilde{w}_{-n}=\tilde{w}_{n}^{\ast }$ , which means that eigenvalues are either real or complex conjugate.

Note that the detuned disturbances with $\tilde{{\it\alpha}}\neq 0$ affect neither the mean pressure gradient nor the flow rate. Both of these quantities are associated with a zero-wavenumber mode, which occurs only for the fundamental disturbances ${\it\epsilon}=0$ when the zeroth harmonic of the streamwise velocity perturbation is an odd function of  $z$ , i.e. ${\tilde{u}}_{0}(-z)=-{\tilde{u}}_{0}(z)$ . This is the case for the fundamental modes with $z$ -parities opposite to those of the travelling-wave state which are discussed in the next section.

The problem posed by (3.17) corresponds to the case of fixed flow rate. For $\tilde{{\it\alpha}}_{0}=0$ we have $\boldsymbol{D}_{\tilde{{\it\alpha}}_{0}}\equiv \boldsymbol{e}_{z}(\text{d}/\text{d}z)$ and thus for $n=0$ (3.17) can be integrated once. This results in a linearized counterpart of (3.13) for $\tilde{{\it\psi}}_{0}^{\prime }={\tilde{u}}_{0}$ containing a constant of integration $\tilde{P}_{0}$ which represents perturbation of the mean pressure gradient:

(3.18) $$\begin{eqnarray}{\tilde{u}}_{0}^{\prime \prime }-(\tilde{{\it\lambda}}+\mathit{Ha}^{2}){\tilde{u}}_{0}=\mathop{\sum }\limits_{m}({\hat{w}}_{m}^{\ast }{\tilde{u}}_{m}+{\hat{u}}_{m}^{\ast }\tilde{w}_{m})^{\prime }+\tilde{P}_{0}.\end{eqnarray}$$

Using this equation with $\tilde{P}_{0}=0$ instead of the original one we obtain a linear stability problem for the travelling waves driven by a fixed mean pressure gradient. Equation (3.18) can integrated once more leading to an equivalent equation in terms of $\tilde{{\it\psi}}_{0}$ which can be used instead (3.18). Alternatively, the integrated equation can be used as an effective but rather complicated boundary condition replacing the fixed flow rate condition $\tilde{{\it\psi}}_{0}(1)=0$ in the original equation (3.17) for $\tilde{{\it\psi}}_{0}$ (Soibelman & Meiron Reference Soibelman and Meiron1991).

4. Numerical method

The problem is solved numerically using a Chebyshev collocation method with the Chebyshev–Lobatto nodes

(4.1) $$\begin{eqnarray}z_{i}=\cos \left(i{\rm\pi}/2M\right),\quad i=0,\dots ,M,\end{eqnarray}$$

at which the discretized solution $({\hat{w}}_{n},{\hat{u}}_{0})(z_{i})=(\tilde{w}_{n},{\tilde{u}}_{0})_{i}$ is sought in the upper half of the channel. The reduction to half the channel is due to the following symmetries implied by (3.9), which is the nonlinear term of the vorticity (3.6). The $z$ -symmetry of this term is the same of as that of ${\hat{w}}_{n}$ governed by (3.6), and it is determined by the first derivative of the terms quadratic in ${\hat{w}}_{n}$ . Note that the first derivative inverts $z$ -symmetry and the second derivative conserves it. Since the product of two harmonics with indices $m$ and $n$ produces two harmonics with the same $z$ -symmetry and indices $n\pm m$ , which are separated by an even number $2m$ , all harmonics with the same parity of indices have the same $z$ -symmetry. Harmonics with even indices are produced by the products of harmonics with indices of the same parity, which have the same $z$ -symmetry. Thus, these products are always even functions of $z$ . The inversion of the $z$ -symmetry by the first derivative in the nonlinear term results in ${\hat{w}}_{n}$ with even $n$ being an odd function of $z$ . All harmonics with odd indices can have either odd or even $z$ -symmetry depending on the symmetry of the fundamental mode $(n=1)$ , which is usually determined by the linear stability analysis. This gives rise to two types of possible 2D solutions satisfying ${\hat{w}}_{n}(-z)=(-1)^{n+1}{\hat{w}}_{n}(z)$ and ${\hat{w}}_{n}(-z)=-{\hat{w}}_{n}(z)$ , which are referred to as even or odd depending on the $z$ -symmetry of the transverse velocity harmonics with odd indices. In secondary stability analysis, small-amplitude perturbations of each of these two solution branches splits into two further symmetry types. The first type of possible perturbation consists of the harmonics with the same $z$ -parities as the travelling-wave solution itself, i.e. $\tilde{w}_{n}(-z)=(-1)^{n+1}\tilde{w}_{n}(z)$ and $\tilde{w}_{n}(-z)=-\tilde{w}_{n}(z)$ for even and odd branches, respectively. The second type of perturbation admitted by (3.17) consists of harmonics with $z$ -symmetries opposite to those of the travelling-wave solution, i.e. $\tilde{w}_{n}(-z)=(-1)^{n}\tilde{w}_{n}(z)$ and $\tilde{w}_{n}(-z)=\tilde{w}_{n}(z)$ for even and odd branches, respectively (Pugh & Saffman Reference Pugh and Saffman1988). Note that for the even 2D solution branch, which consists of harmonics with alternating $z$ -parities, one type of perturbation is changed into the other by the transformation ${\it\epsilon}\rightarrow {\it\epsilon}\pm 1$ which effectively shifts the index of harmonics in (3.16) by one.

Equations are approximated at the internal collocation points $1\leqslant i\leqslant M$ by using differentiation matrices which express the derivatives in terms of the collocation variables $(\unicode[STIX]{x1D66C}_{n},\unicode[STIX]{x1D66A}_{0})_{i}$ . Odd or even $z$ -symmetry of each particular harmonic is taken into account by the reduced-size differentiation matrices based on the collocation points only in one half of the channel. The boundary conditions (3.3) are imposed at the boundary point $i=0$ (Peyret Reference Peyret2002). When series (3.4) are truncated at $n=\pm N$ , (3.11) reduce to $N\times M$ complex algebraic equations with respect to the same number of complex unknowns $\unicode[STIX]{x1D66C}_{n}$ for $n=1,2,\dots ,N$ . Note that $\unicode[STIX]{x1D66C}_{0}\equiv 0$ and $\unicode[STIX]{x1D66C}_{-n}=\unicode[STIX]{x1D66C}_{n}^{\ast }$ . The equations also contain $M$ real unknowns ${\tilde{u}}_{0}$ , which are governed by the same number of real equations resulting from the collocation approximation of (3.13). The equations are linear and, thus, directly solvable for ${\tilde{u}}_{0}$ in terms of $\unicode[STIX]{x1D66C}_{n}$ by inverting the respective matrices. This leads to a system of $N\times M$ nonlinear complex equations for the same number of complex unknowns $\unicode[STIX]{x1D66C}_{n}$ . Since the equations also contain $\unicode[STIX]{x1D66C}_{n}^{\ast }$ , the actual unknowns are the real and imaginary parts of $\unicode[STIX]{x1D66C}_{n}$ , which need to be determined by solving the same number, i.e. $2N\times M$ , of real equations represented by the real and imaginary parts of the original complex equations.

However, there is one more unknown: the frequency ${\it\omega}$ , which is the eigenvalue of this nonlinear problem, and needs to be determined along with $\unicode[STIX]{x1D66C}_{n}$ . There are two ways to balance the number of unknowns and equations. First, as the problem is homogeneous, non-trivial solution requires a solvability condition to be satisfied, which provides another equation analogous to the characteristic determinant in the case of a linear eigenvalue problem. The second possibility, which is used here, follows from the fact that owing to the translational invariance of the problem $\unicode[STIX]{x1D66C}_{n}$ is defined up an arbitrary phase. The phase, which determines the $x$ -offset of the wave, can be fixed by imposing the condition $\text{Im}[\unicode[STIX]{x1D66C}_{1,i}]=0$ at some collocation point $i$ where the solution is not already fixed by boundary or symmetry conditions. This is equivalent to setting $\unicode[STIX]{x1D66C}_{1,i}=A$ , where $A$ is a real parameter defining the amplitude of the transverse velocity. Thus, the number of unknowns is reduced by one and the system of $2N\times M$ nonlinear algebraic equations can be written in the general form

(4.2) $$\begin{eqnarray}\unicode[STIX]{x1D641}(\unicode[STIX]{x1D66C},A,{\it\omega};{\it\alpha},\mathit{Re})\unicode[STIX]{x1D66C}=\mathbf{0},\end{eqnarray}$$

where $\unicode[STIX]{x1D66C}$ are the real and imaginary parts of $\unicode[STIX]{x1D66C}_{n}$ normalized with the real amplitude $A$ and $\unicode[STIX]{x1D641}$ is a real square matrix of size $2N\times M$ depending on the listed parameters. For given ${\it\alpha}$ and $\mathit{Re}$ , this problem can be solved by the Newton–Raphson method with respect to $A$ , ${\it\omega}$ and $2N\times M-2$ unknown $\unicode[STIX]{x1D66C}$ . In some cases, instead of $\mathit{Re}$ , it is more convenient to fix $A$ and then to solve for $\mathit{Re}$ depending on $A$ and ${\it\alpha}$ .

The solution is traced using a quadratic extrapolation along the arclength in logarithmic coordinates. For a general function $f(p)$ of an argument $p$ the scale-independent arclength element is defined as ${\it\delta}s^{2}=\ln ^{2}(1+({\it\delta}f/f))+\ln ^{2}(1+({\it\delta}p/p))$ . Starting with a reference arclength based on the solution at the first three chosen parameter values, a subsequent parameter value and an initial guess for the Newton–Raphson method are extrapolated from the previous three values. When Newton–Raphson iterations fail, the step size along the arclength is reduced until a solution is recovered, and then gradually increased to its original value as the solution is successfully traced.

6. Summary and conclusions

The present study was concerned with 2D nonlinear travelling-wave states in MHD channel flow subject to a transverse magnetic field. Such states are thought to mediate transition to turbulence, which is known to take place in this flow at Reynolds numbers more than two orders of magnitude below the linear stability threshold. Using the Newton–Raphson method, we determined the extension of such 2D nonlinear travelling waves below the linear stability threshold, at which they appear as the result of a subcritical bifurcation. Starting from the non-magnetic plane Poiseuille flow, where these travelling-wave states extend from the linear stability threshold $\mathit{Re}_{l}=5772$ down to $\mathit{Re}_{n}=2939$ , and gradually increasing the magnetic field, we found that in a sufficiently strong magnetic field $(Ha\gtrsim 10)$ such states extend down to $\mathit{Re}_{n}\sim 6.50\times 10^{3}\mathit{Ha.}$ Although this Reynolds number is almost an order of magnitude lower than the linear stability threshold of the Hartmann flow, it is still more than an order of magnitude greater than that at which turbulence is observed in this type of flow.

Besides the solution which evolves from plane Poiseuille flow as the magnetic field is increased, we found another inherently magnetohydrodynamic solution which bifurcates subcritically from the Hartmann flow at $Ha\approx 6.5$ . The two solutions differ by their $z$ -symmetries: the harmonics of the transverse velocity are odd functions of $z$ for the latter and they have alternating parities starting with even one for the former. The odd-symmetry solution was found to have two branches. The first branch, which exists in a limited range of the magnetic field strength with $2.8\lesssim Ha\lesssim 20$ , has a relatively long wavelength and consists of a pair of intense localized vortices. The second branch, which emerges at $Ha\approx 10$ with a much shorter wavelength, has very similar characteristics to the original hydrodynamic branch, and becomes practically indistinguishable from the latter when $Ha\gtrsim 20$ . We also showed that the lower-amplitude 2D travelling waves, as in the non-magnetic case, are invariably unstable with respect to 2D infinitesimal superharmonic disturbances down to the limiting value of the Reynolds number based on the flow rate.

Subcritical 2D travelling waves in the Hartmann flow are also likely to be unstable to more general three-dimensional disturbances similar to those considered by Orszag & Patera (Reference Orszag and Patera1983) in the non-magnetic case. Three-dimensional equilibrium states bifurcating either from 2D travelling waves (Ehrenstein & Koch Reference Ehrenstein and Koch1991) or infinity (Waleffe Reference Waleffe2001), as in plane Poiseuille flow, may extend to significantly lower Reynolds numbers and, thus, provide a more adequate threshold for the onset of turbulence in the Hartmann flow. Such a possibility is supported by the interaction of only two mirror-symmetric oblique waves with the resulting 2D second harmonic considered by Zinovév & Shtern (Reference Zinov’ev and Shtern1987), who found, for the Hartmann layer, $R_{n}\approx 4670$ . As for the 2D waves considered in this study, it is likely that a more adequate three-dimensional model including a sufficient number of higher harmonics could result in a substantially lower $R_{n}$ .