inductionless approximation holds in the present case of small $Pm$, the code updates the induction Eq. (2) in time in a general manner without invoking this approximation. Further details of the code and its validation tests on HD and MHD problems in TC flows can be found in [@Guseva_etal15].

The cylindrical flow domain in nondimensional units spans the radial range $(R_i,R_o)=(1,2)$ divided into $N_r=100$ grid points, azimuthal range $(0, 2\pi)$ and axial range $(0,L_z)=(0,10)$ with, respectively, $N_{\phi}=80$ and $N_z=400$ spatial Fourier modes. The radial grid points are placed at collocation points of Chebyshev polynomials (despite the code is based on finite differences in $r$) that ensures higher resolution near the cylinder surfaces. Thus, although the HMRI is thought to be predominantly axisymmetric at least in the linear regime, we still allow for non-axisymmetric modes to naturally develop in the nonlinear regime by encompassing azimuthal wavenumbers, $m$, in the range $(-N_{\phi}/2,
N_{\phi}/2)$. A resolution test we performed (not shown here) showed that the adopted resolution is well sufficient for the present problem, capturing a whole range of wavenumbers from dynamically important intermediate ones, corresponding to the most unstable HMRI modes, down to wavenumbers dominated by viscous dissipation, as indicated by the 