mainly determines the linear dynamics of HMRI, it is natural to explore how the nonlinear dynamics and saturation properties of HMRI depend on this parameter, which is our main goal. This is the main reason why we decided to follow the dynamics as a function of $\Lambda$. We anticipate that, like in the linear theory, this parameter will be decisive in the nonlinear outcome of the HMRI. Previous related nonlinear studies [@Szklarski_Ruediger06; @Szklarski07] focused mostly on the variation of the saturated state’s characteristics with the Reynolds number at fixed Hartmann number. Besides, for these parameters, HMRI was in the weakly nonlinear regime. In this connection, we mention that the interaction parameter also plays an important role in the dynamics of forced MHD turbulence with mean field in the small-$Pm$ regime [@Zikanov_Thess98].

Numerical method
----------------

We solve the basic Eqs. (1)-(3) using the pseudo-spectral code from [@Guseva_etal15]. It is based on the Fourier expansion in the axial $z$- and azimuthal $\phi$-directions and the finite-difference method in the radial direction. The nonlinear terms are calculated using the pseudospectral method and are de-aliased using the 3/2-rule. Integration in time is done with a second-order scheme based on the implicit Crank-Nicolson method. Although the 