© The Authors 2023

1. Introduction

Sources of energy input in the interstellar medium (ISM) are of stellar or galactic origin, and they lead the observed turbulence in the ISM. Supernova explosions (SNe), stellar winds, protostellar outflows, and even ionizing radiation (by decreasing order of importance) are known to contribute to injecting energy into the ISM. It is known that SNe are the main contributors of energy and can solely provide the energy needed to produce the turbulence observed in the ISM (e.g., Gressel et al. 2008). They are also considered one of the main ingredients in determining galactic dynamos (see Brandenburg & Ntormousi 2023, for a recent review). Still, their effect on the amplification of galactic magnetic fields remains poorly understood. It is unclear how they contribute to the seeding and amplification of such fields to the observed values both for Milky Way-like galaxies (e.g., Ntormousi 2018; Ntormousi et al. 2020) and for high-redshift ones (Ntormousi et al. 2022). If they are seen as purely spherical, SNe should not (in principle) produce any vorticity in the surrounding medium.

However, there are many other aspects of the Galactic medium that ought to be taken into account. Stratification, differential rotation, shear, or shocks are examples of physical mechanisms which might help generate vortical flows from the available energy. Three-dimensional (3D) simulations of such environments can be quite demanding, especially considering the supersonic feature of the flow. This field was pioneered by Korpi et al. (1999), who studied how SNe produce turbulence in a multiphase ISM, and then by Mac Low et al. (2005), who analyzed the distribution of pressure in a magnetized ISM. However, in these simulations, the generation of vorticity, which in turn is essential for amplifying the magnetic field, remains somewhat unclear. The question of whether vorticity can be amplified by a purely curl-free (compressible, or irrotational) forcing (e.g., Federrath et al. 2010, 2011) or not (Mee & Brandenburg 2006; Del Sordo & Brandenburg 2011) is still a matter of debate.

Mee & Brandenburg (2006) and Del Sordo & Brandenburg (2011) used localized random expansion waves as forcing terms only in the momentum equation. This potential forcing induces spherically symmetric successive accelerations in random places during a given time interval, reaching locally trans-sonic regimes and leading to a somewhat similar SNe driving. Such an irrotationally forced flow was found not to generate any vorticity up to 512³ mesh points (Mee & Brandenburg 2006). However, when rotation, explicit shear or baroclinicity were added, Del Sordo & Brandenburg (2011) found that vorticity is produced. Magnetic fields were added to the most simple setup (no rotation or shear, nor baroclinicity) by Mee & Brandenburg (2006), who found no evidence of magnetic field amplification. On the contrary, they found that this type of flow has a highly destructive effect on the magnetic field at the smallest scales. In general, the role of SN feedback on regulating star formation rate and the structure of the galactic disk is still not completely understood (e.g., Hennebelle & Iffrig 2014) and it may depend on how SNe are implemented in simulations. The nature of the forcing may have an effect on several features of the ISM, besides the growth of the magnetic field. For instance, Mathew et al. (2023) showed how a purely compressible forcing in magnetohydrodynamic (MHD) simulations of star cluster formation influences the initial mass function.

This work is a continuation of the research done by Mee & Brandenburg (2006) and Del Sordo & Brandenburg (2011). We aim to study how shear, rotation, and baroclinicity, in combination with irrotational forcing in the form of localized random expansion waves, contribute to the decay or growth of an initial magnetic field seed. We want to explore the possible dynamo action in a parameter space defined by the forcing parameters, the rotation and shear rates, the Reynolds and Prandtl numbers, and the initial magnetic seed. Our work concentrates on the subsonic and transonic regime, with the aim of establishing the conditions under which a purely compressible forcing may amplify vorticity and magnetic fields.

Other works have also studied the appearance of small-scale dynamo (SSD), which is dynamo action at length scales equal or smaller than the forcing scale, as a consequence of irrotational forcing. Porter et al. (2015) performed simulations of subsonic turbulence in intra-galactic cluster medium with a purely curl-free (compressible) forcing, finding that dynamo can be excited when divergence-free (solenoidal) modes arise. Achikanath Chirakkara et al. (2021) found turbulent dynamo in a highly subsonic regime by using a turbulent driving force in Fourier space either with solenoidal or with purely compressive contributions, or with a combination of them, following the approach of Federrath et al. (2010) and Federrath et al. (2011).

Other simulations have achieved magnetic field growth with more complex forcing of the same type. Gent et al. (2013a,b) simulated a multi-phase ISM randomly heated and stirred by SNe, in the energy equation, with a stratified, rotating and shearing local domain with a galactic gravitational potential and shock diffusivities. The obtained flow dynamics are robust when three major phases are used (defined with temperature ranges) and when the parameters are adjusted to reproduce individual SN remnants. Käpylä et al. (2018) simulated such events with a combination of mass transfer during SNe and stellar formation, energy deposition in energy equation, and a stably stratified, rotating medium mimicking the galactic plane. They also obtained a SSD in such simulations. This was also confirmed in subsequent studies by Gent et al. (2021, 2023) of a multi-phase ISM. Recently, Seta & Federrath (2022) studied turbulent dynamo in a multi-phase ISM as well, providing a detailed analysis of vorticity sources in their simulations.

These models are fairly complex as they want to reproduce what happens in galactic environment. Our approach has instead been to work with a much simpler model to shed light on the minimum ingredients needed to amplify a magnetic field up to equipartition values in the presence of a purely irrotational forcing. This is in the spirit of studying general aspects of dynamo generated magnetic fields and discuss their applications both in galactic environment, as well as in planetary and stellar ones. The analytical mean field approach by Krause & Raedler (1980) and Rädler & Stepanov (2006) predicts the presence of an electromotive force and, hence, an amplified mean magnetic field, in the presence of a mean flow. While this is surely happening in the case of helical flows (see e.g., Brandenburg & Subramanian 2005; Rincon 2019), a mean electromotive force can (in principle) be observed also with non-helical turbulence in the presence of rotation (Rädler 1968, 1969) or with a large-scale flow with associated vorticity (Rogachevskii & Kleeorin 2003). Large-scale dynamos was found in numerical experiments employing a non-helical forcing to drive small-scale turbulence embedded in a large-scale shear flow (e.g., Yousef et al. 2008; Brandenburg et al. 2008; Singh & Jingade 2015). However, the forcing functions in these calculations did allow the injection of vorticity on small scales. In this work, we want instead to test the possibility of driving a dynamo process with turbulence driven in a purely compressible way, either in the absence or in the presence presence of large-scale flows.

This paper is organized as follows: in Sect. 2 we describe the MHD model and the organization of our numerical experiments. In Sect. 3 we show the results and in Sect. 4 we discuss the implication of these results in the framework of galactic dynamo. Finally, we present our conclusions in Sect. 5.

2. Model and numerical methods

2.1. MHD equations with rotation and shear

We ran numerical simulations using the Pencil Code, Pencil Code Collaboration (2021)¹. This is a non-conservative, high-order, finite-difference code (sixth-order accurate in space and third-order Runge-Kutta in time), which we employed to solve the non-ideal compressible MHD equations.

We considered either a rigidly rotating frame, with angular velocity Ω = (0, 0, Ω) in the z-direction or a differential velocity (or shear) in the y-direction given by ${\mathbf{u}}^S=(0,u_y^S(z),0)$. We did not consider stratification. We solved the equations in the co-rotating reference frame, for which the mass and momentum conservation equations are expressed as:

$$\begin{array}{cc}\hfill \frac{\mathrm{D}ln\rho }{\mathrm{D}t}& =-\mathrm{\nabla }·\mathbf{u},\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\mathrm{D}{\mathbf{u}}^{\prime }}{\mathrm{D}t}& =-\frac{\mathrm{\nabla }p}{\rho }+\frac{\mathbf{J}\times \mathbf{B}}{\rho }-2\mathrm{\Omega }\times \mathbf{u}-u_z\frac{\partial u_y^S}{\partial z}\widehat{\mathbf{y}}+{\mathbf{F}}_{\nu }+\mathbf{f}+{\mathbf{f}}_{\mathrm{s}},\hfill \end{array}$$(2)

where: u(t) = (u^S + u′(t)) is the total velocity field combining the fixed shearing velocity with the turbulent velocity u′; ρ is the mass density; p is the pressure; B the magnetic field; and J the electrical current density J := (∇×B)/μ₀ (where μ₀ is the vacuum permeability).

Importantly, the advective derivative operator, D/Dt := ∂/∂t + u ⋅ ∇, applies to the total velocity, that is, to both the evolving turbulent velocity u′ (kept in the left-hand side) and the fixed shearing velocity, u^S. We keep the latter in the right-hand side, in the form of the only term $-\mathbf{u}·\mathrm{\nabla }{\mathbf{u}}^{\mathbf{s}}=-u_z\frac{\partial u_y^S}{\partial z}$. The other source terms are the viscous force, F_ν = ρ⁻¹∇⋅(2ρνS), where the traceless rate of strain tensor S is S_ij = (1/2)(u_i, j + u_j, i − (1/3)δ_ij∇⋅u), the external expansion wave forcing, f, and a shearing forcing, f_s (see below for their definitions).

In order to close the system of equations, we consider two types of equation of state (EoS): 1) a simple barotropic EoS $p(\rho )=c_{\text{s}}^2\rho$, where we fix the value of the sound speed c_s = 1 and 2) an ideal EoS (dubbed also baroclinic case), p(ρ, T) = ρR_gT, with R_g the specific gas constant and T the temperature; in this case, the sound speed squared is $c_{\text{s}}^2=(\gamma -1)c_{\text{p}}T$, where we fix the adiabatic index γ = c_p/c_v = 5/3 (corresponding to a monatomic perfect gas), and c_p and c_v are the specific heats at constant pressure and constant volume, respectively.

In the latter case, the energy equation is expressed in terms of ρ, the specific entropy, s, and T:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill T\frac{\mathrm{D}s}{\mathrm{D}t}=2\nu \mathit{S}\otimes \mathit{S}+{\rho }^{-1}\mathrm{\nabla }(c_{\mathrm{p}}\rho \chi \mathrm{\nabla }T)+{\rho }^{-1}\eta {\mu }_0{\mathit{J}}^2-\frac{1}{{\tau }_{\mathrm{cool}}}(c_{\mathrm{s}}^2-c_{\mathrm{s}0}^2),\end{array}\end{array}$$(3)

where χ is the thermal diffusivity, η is the magnetic diffusivity, c_s0 is the initial, uniform sound speed (proportional to the initial temperature) and τ_cool regulates the timescale of the simple effective cooling term, introduced to avoid an indefinite heating by viscous and resistive dissipation.

The time evolution of the magnetic field is done by uncurling the usual induction equation, using ∇ × A = B, where A is the vector potential. For our system, the coupling between the shearing velocity field and the magnetic field reads (Brandenburg et al. 1995):

$$\begin{array}{c}\hfill \frac{\partial \mathbf{A}}{\partial t}=\mathbf{u}\times (\mathrm{\nabla }\times \mathbf{A})+\eta {\mathrm{\nabla }}^2\mathbf{A}.\end{array}$$(4)

The solenoidal condition, ∇ ⋅ B = 0, is therefore fulfilled at all times, avoiding any possible spurious growth coming from divergence cleaning procedures.

2.2. Forcing

We impose an irrotational forcing, as a gradient of a Gaussian potential function:

$$\begin{array}{c}\hfill \mathbf{f}(\mathbf{x},t)=\mathrm{\nabla }\varphi (\mathbf{x},t)=\mathrm{\nabla }(Kexp\{-{(\mathbf{x}-{\mathbf{x}}_f(t))}^2/R^2\}),\end{array}$$(5)

where x = (x, y, z) is the position vector, x_f(t) is the random position corresponding to the center of the expansion wave, R is the radius of the Gaussian, and N is the normalization factor. We will make use of k_f as the wavenumber corresponding to the scale of the adopted forcing:

$$\begin{array}{c}\hfill k_f=\frac{2}{R}.\end{array}$$

The time dependence of the forcing enters in the coherence time Δt of the expansion waves (i.e., the time interval after which we consider a new random position, x_f). We consider two different cases: in the first one x_f changes at every time-step (which is adaptive), Δt = δt, while in the second case we define an interval forcing time Δt > δt. We chose the normalization factor to be:

$$\begin{array}{c}\hfill K={\varphi }_0\sqrt{c_{\mathrm{s}0}R/\mathrm{\Delta }t},\end{array}$$

where ϕ₀ controls the forcing amplitude and has dimensions of the velocity squared.

After an initial transitory phase, the simulations reach a stationary state, over the course of which, the main average quantities maintain a saturated value. In particular, we will look at the root mean square of the velocity, u_rms. In turn, this is used to define the fundamental timescale of our problem, that we will call turnover time, as

$$\begin{array}{c}\hfill t_{\mathrm{turn}}={(k_f{u}_{\mathrm{rms}})}^{-1}.\end{array}$$(6)

The turnover time can be understood as the average time for the fluid to cross en explosion width.

The root mean square values of velocity u_rms and vorticity ω_rms (see Sect. 2.4 for its definition) are used to define the following dimensionless numbers:

$$\begin{array}{cc}& \text{Re}=\frac{u_{\mathrm{rms}}}{\nu k_f},\text{Rm}=\frac{u_{\mathrm{rms}}}{\eta k_f},{\text{Re}}_{\omega }=\frac{{\omega }_{\mathrm{rms}}}{\nu k_f^2},\hfill \end{array}$$(7)

$$\begin{array}{cc}& \text{Ma}=u_{\mathrm{rms}}/c_{\mathrm{s}},\text{Pm}=\nu /\eta ,\hfill \end{array}$$(8)

which are the Reynolds number, magnetic Reynolds number, vorticity Reynolds number, Mach number, and magnetic Prandtl number, respectively. Analogously, we shall consider the maximum values of the velocity, u_max, to define, for example:

$$\begin{array}{c}\hfill {\text{Re}}_M=\frac{u_{max}}{\nu k_f},{\text{Ma}}_M=u_{max}/c_{\mathrm{s}}.\end{array}$$(9)

Regarding the magnetic fields, we shall consider the root-mean square, b_rms, closely related to the magnetic energy density.

2.3. Shear

In order to maintain a given shear velocity background along the y-direction $u_y^S$, we use a forcing term in the y-component of the velocity: ${\mathbf{f}}_{\mathrm{s}}=\tau (u_y^S-u_y)\widehat{\mathbf{y}}$, where τ is the timescale of the forcing (which we keep unity in all cases).

We choose to apply a sinusoidal background flow for the shearing term:

$$\begin{array}{c}\hfill u_y^S=Acos(kz).\end{array}$$(10)

As z ranges from −π to π (see Sect. 2.5), we have set k = 1 which allows simple periodic boundary conditions in the three directions. This profile is similar to what is used by Skoutnev et al. (2022) for studying dynamo in stellar interiors with a non-helical forcing, and was employed by Käpylä et al. (2009, 2010) in the context of stratified convective medium.

We discarded the use of a more standard linear shear term $u_y^S$ = Sz, since the implemented shearing boundary conditions were giving a spurious growth of vorticity at the boundaries (see Appendix C for a more detailed explanation).

2.4. Vorticity

The vorticity is defined as the curl of the velocity field (ω = ∇ × u) and it quantifies the turbulent motions generated in fluid flows. Vortical motions are directly connected to the turbulent dynamo action, which is the mechanism by which our model may grow or maintain the magnetic field. The evolution of vorticity can be obtained by taking the curl of the total velocity evolution equation (∂u/∂t):

$$\begin{array}{cc}& \frac{\partial \mathit{\omega }}{\partial t}=\mathbf{\nabla }\times (\mathbf{u}\times \mathit{\omega })+\mathrm{\nabla }\times {\mathbf{F}}_{\mathrm{visc}}-2\mathrm{\nabla }\times \mathbf{\Omega }\times \mathbf{u}+\frac{\mathrm{\nabla }\rho \times \mathrm{\nabla }p}{{\rho }^2}+\hfill \\ \hfill & +\mathrm{\nabla }\times \frac{\mathbf{J}\times \mathbf{B}}{\rho }+\mathrm{\nabla }\times {\mathbf{f}}_{\mathrm{s}}.\hfill \end{array}$$(11)

This equation is not directly evolved in our calculations and it serves only to understand the production mechanisms of vorticity, as each of the different terms of the RHS can generate or destroy it. We note that the forcing is irrotational by construction (∇ × f = 0), thus it does not appear in Eq. (11) and it cannot directly generate vorticity. However, the forced flow could help generate vorticity by the first two right-hand side terms: the first term (associated to turbulent motions) and the viscous terms. These two terms are quite analogous to the induction equation, hinting that turbulence amplification should be as difficult to achieve as dynamo action, in the absence of any other forcing. Mee & Brandenburg (2006) proved that this forcing produced no measurable vorticity: the amount of vorticity produced showed a decreasing dependence on numerical resolution, hence they concluded its nature was numerical.

The other vorticity sources appearing in the RHS are: rotation, baroclinic (zero for the barotropic case), Lorentz and shear terms, respectively. We expect the amount of vorticity generated by rotation to be proportional to the Coriolis number Co = 2Ωt_Ω, for a fixed viscosity and a rotation that is not too rapid. Also, the t_Ω here is a typical timescale of the system that we can assume t_Ω ≈ t_turn. The contribution of the rotation term to the vorticity, ω_rms, can be then written as:

$$\begin{array}{cc}& \frac{\partial \mathit{\omega }}{\partial t}=\cdots -2\mathrm{\nabla }\times \mathbf{\Omega }\times \mathbf{u}\hfill \\ \hfill & \to {\omega }_{\mathrm{rms}}\approx 2\mathrm{\Omega }{t}_{\mathrm{turn}}{u}_{\mathrm{rms}}{k}_f\to \frac{k_{\omega }}{k_f}\approx 2\mathrm{\Omega }{t}_{\mathrm{turn}},\hfill \end{array}$$(12)

where we have defined:

$$\begin{array}{c}\hfill k_{\omega }=\frac{{\omega }_{\mathrm{rms}}}{u_{\mathrm{rms}}},\end{array}$$(13)

as a proxy of how much of the available kinetic energy is in the form of vortical flows and provides an idea of the typical wavenumber of vortical structures.

The shearing term contribution to Eq. (11) is:

$$\begin{array}{cc}& \frac{\partial \mathit{\omega }}{\partial t}=\cdots +\mathrm{\nabla }\times (-u_z\frac{\partial u_y^S}{\partial z}\widehat{}\}y)\to \hfill \\ \hfill & \to \frac{{\omega }_{\mathrm{rms}}}{t_{\mathrm{turn}}}\approx k_f{u}_{\mathrm{rms}}{k}_f{u}^S\to \frac{k_{\omega }}{k_f}\approx A{t}_{\mathrm{turn}}{k}_f.\hfill \end{array}$$(14)

Finally, the baroclinic contribution ∇ρ × ∇p (misalignment between the pressure and density gradients) is proportional to ∇T × ∇s (see Eqs. (12) and (13) in Del Sordo & Brandenburg 2011). To study the effect of such baroclinic term we can monitor the mean angle between these two gradients:

$$\begin{array}{c}\hfill \text{sin}\theta =\frac{⟨\mathrm{\nabla }T\times \mathrm{\nabla }s⟩}{⟨\mathrm{\nabla }T⟩⟨\mathrm{\nabla }s⟩},\end{array}$$(15)

so that the evolution of vorticity can be seen as:

$$\begin{array}{c}\hfill \frac{\partial \mathit{\omega }}{\partial t}=\cdots -\mathrm{\nabla }T\times \mathrm{\nabla }s\to {\omega }_{\mathrm{rms}}\approx t_{\mathrm{turn}}{(\mathrm{\nabla }T)}_{\mathrm{rms}}{(\mathrm{\nabla }s)}_{\mathrm{rms}}\text{sin}\theta .\end{array}$$(16)

2.5. Boundary and initial conditions

The simulation domain consists of a uniform, cubic grid mesh [ − π, π]³, with triply periodic boundary conditions. We consider resolutions varying from 128³ up to 512³, which are enough to assess our problem.

We adopted non-dimensional variables by measuring speed in units of the initial sound speed, c_s0, length in units of 1/k₁ where k₁ is the smallest wave number in the periodic domain, implying that the non-dimensional size of the domain is (2π)³.

As the initial conditions, pressure (and entropy and temperature in the baroclinic case), and density are set constant and with value 1 throughout the box. The flow is initially still, v = 0.

Finally, we set a weak initial magnetic field. Even though the initial flow is static, it is soon shaken by the expansion waves, so that the small initial magnetic seed does not influence the flow. For the initial magnetic topology, we consider either a uniform field in one given direction, or a spatially random field. The latter is assigned by picking, at each point, uncorrelated random values to the three components of the initial magnetic potential. The corresponding initial magnetic energy spectrum follows a k⁴ power law, as reported by Mee & Brandenburg (2006).

3. Results

3.1. Explored parameters

We ran a series of different numerical models to explore the role of each term in the RHS of Eq. (2). The forcing and viscous terms are always included in all of our models, while we separately consider the activation of either rotation or shear. The forcing amplitudes are high enough to create transonic flows, but we do not investigate supersonic flow. We have restrained ourselves from the usage of shock viscosities to avoid the introduction of numerical artifacts and to concentrate on the effect of a uniform viscosity. Most simulations have Pm = 1 (and therefore Rm = Re) up to values of a few hundred. In some runs, we went to Pm up to 100.

We start from some basic hydrodynamic (HD) simulations, similar to Del Sordo & Brandenburg (2011). We then expand them to magnetized cases, exploring broadly the parameter space. Unless otherwise specified, the expansion width has been chosen to be R = 0.2 (thus, k_f is 10).

In Tables 1–6, we list the input parameters of the performed simulations. The tagging names starting with “H” indicate HD runs, those starting with “M” indicate MHD runs. By default, we employed the barotropic EoS, while “B” identifies the baroclinic cases. The first number stands for the value of Ω. By default, we employ Δt = 0.02, s in the end of the name identifies Δt = 1, while c means Δt = δt. Unless stated otherwise (by the suffix 128, 512), the resolution is 256³. The width (R) and constant (ϕ₀) of the Gaussian are marked with the number following W and F in the run names, and, if not stated, they are assumed to be 0.2 and 1, respectively. The cooling term of the baroclininc cases is used only in the runs listed in Table 4. The initial magnetic field by default takes the form of random values of the potential vector component, unless a background uniform field in a given i direction is assigned (Bi). Finally, we use “high” or “low” to indicate the changes in the viscous terms.

We monitor the volume-integrated energies and the above mentioned root mean square values of quantities. Additionally, we look at the spectral energy distribution of velocity (total or turbulent), vorticity, and magnetic field, as a function of the wavenumber, k (which includes the factor 2π, see Appendix A). Spectra are also averaged over different times to filter out fluctuations, which are quite noticeable on the smallest scales. To evaluate the rotationality of the flow we also make use the Helmholtz decomposition for the velocity field, as defined in Appendix B.

3.2. Effects of forcing and rigid rotation on the flow

Table 1 lists some examples of HD runs, without magnetic field. To visualize the general behavior at early times, Fig. 1 shows the values of u² in the xy-plane for one of them, the representative model H_0cW1. Rapidly, the imprints of the most recent expansion waves are visible on top of the turbulent state. Such an evolution leads to a homogeneous turbulent flow, which becomes stationary when the dissipating forces counterbalance the forced energy injection.

Fig. 1. Slice in the xy-plane showing the value of u2 at t/tturn = 0, 0.01 (top), 0.1, and 1 (bottom) for the non-rotating model H_0cW1 of Table 1.

Such features are seen also in the MHD cases without shear, for which Tables 2 and 3 list the barotropic and baroclinic configurations, respectively, considering or not rigid rotation. In these same tables, we provide the diagnostic magnitudes. As a first comparison, in Fig. 2 we can see how different representative simulations (with or without rotation, and for the two EoS) saturate their Re and Ma values within ∼1000 t_turn. A faster rotation increases the final saturation value for the velocities, but the growth takes places at a similar rate, so simulations with faster rotation take longer to saturate. The ideal EoS (i.e., baroclinic) cases show a lower saturation value than the barotropic one. We interpret this result as a consequence of the presence of the extra dissipating terms in the energy equation. The maximum values for Re and Ma for these simulations, Re_M and Ma_M, oscillate much more (but with values that are similar for each), namely, ranging from 150 to 300, and 0.3 to 0.6, respectively. These values slightly increase with the rotation rate, and again are smaller in the baroclinic case.

Fig. 2. Comparison of Re and Ma evolution in the HD cases, considering or not the rigid rotation (Ω = 0 or 2) and for different EoS (H vs. HB).

Figure 3 shows the comparison between kinetic spectra at saturation for different EoS and Ω. In general, for all cases kinetic spectra peak around the value k_f (simulations M_0s and M_0 having R = 0.2, show the bump for k ∼ 5 − 10, while H_0cW1 accumulates energy at the largest scales, since R = 1). This is in agreement with the characteristic forcing wavenumber k_f of the adopted forcing, which in the Fourier space is $\propto \exp (-k^2/k_f^2)$. The inertial part of the spectra shows a ∼k⁻² slope. This irrotational forcing does not lead to any vorticity production, which might help explain the difference with respect to the isotropic turbulent slope of k^−5/3. These general spectral features (peak and slope) are compatible with those found by Mee & Brandenburg (2006) and Del Sordo & Brandenburg (2011).

Fig. 3. Time-averaged kinetic spectra at saturation of some representative simulations of Tables 1–3. We notice how the results on large scales (small k) are very similar for all the models, excluding H_0cW1 (dash-dotted green line), which has a forcing with kf = 2. The spectrum with lowest power is H_0s (dashed green line), which has Δt = 1. Models marked with other colors explore variation of Ω and EoS. In all these models the magnetic field decays rapidly and kinetic spectra coincide with the corresponding HD runs of Fig. 2.

Rotation tends to inhibit smaller scales and promotes the accumulation of kinetic energy at the largest scales, also slightly displacing the spectra peak to the left. This results in a steeper slope in the inertial range. A similar (but less noticeable) effect is seen in the baroclinicity, as compared to the barotropic cases. This is in agreement with the lower total kinetic energy, due to the extra diffusive terms.

To quantitatively measure the dependence of vorticity with rotation, we compared the saturation values of Re_ω (see Table 2). The saturation value increases with Ω, but reaches a maximum around Ω = 5, from which rising rotation does not lead to more vorticity. This translates into a saturation of Rossby number, Ro := u_rmsk_f/2Ω, of the order of u_rms. The low amount of vorticity appearing with Ω = 0 appears to be spurious, since it decreases when resolution is increased, as already reported by Mee & Brandenburg (2006). Such spurious contribution is anyway much smaller than the physical one when the rotation is considered.

In order to see if the estimate in Eq. (12) holds, we show in Table 2 some physical quantities of these runs. We find a linear trend valid up to Ω ∼ 3, which we quantify by a trend between some simulations (marked in blue in Table 2) as: k_ω/k_f = (0.180  ±  0.004)Ω + (0.010  ±  0.009). The slope seems much smaller for Ω ≳ 4 (but further simulations would be needed to quantify it).

Finally, the saturated values increase as we decrease Δt from 1 (H_0s) to 0.02 (H_0). There is not much of a difference between Δt being 0.02 or δt, since the latter is indeed on the order 𝒪(10⁻²). Overall, these results are compatible with the idea that this irrotational forcing alone does not produce any significant amount of vorticity.

3.3. Absence of dynamo with rigid rotation only

In all the models of Tables 2 (barotropic) and 3 (baroclinic), the initial magnetic field quickly decays, with or without rotation. In other words, in spite of the vorticity growth induced by rotation, the system does not experience SSD. Indeed, there is no significant difference in the flow between these MHD models the purely HD ones that correspond to them.

An example of a decay of the magnetic spectrum is shown in Fig. 4, where it can be observed how the smallest scales decay faster and the k⁴ slope rapidly changes. Washing away first the smaller scales is natural, as both numerical and physical diffusivities have more influence on the very small scales. After some turnover times the inertial range of the magnetic spectra resemble (in some way) the kinetic ones. Rotation also tends to favor larger scales and increase the slope, although both slopes are steeper than in the kinetic spectra. There is a minimal interaction between the flow and the magnetic field through the first term of the induction, as expressed by Eq. (4).

Fig. 4. Magnetic spectra decay for the M_0 and M_2 models. The spectra evolution is only shown for M_0, and the first corresponds to less than a turnover time when it still shows the initial slope of k4. The arrow indicates the direction of temporal evolution. The last spectrum is approximately at turnover time 20 and it is shown for both models.

The decaying evolution of b_rms, for the initial random values, were found to follow power laws. The parameter fits do not depend on Ω for our runs. With no forcing we obtain very similar decays with the same exponent, hence, finding no substantial difference between M_0 and the corresponding simulation with no forcing. This leads us to think that the forced flow does not have a major destructive influence on the field, contrary to what Mee & Brandenburg (2006) noted, although it does not lead to any growth of the magnetic field either.

We have explored this for a large range of parameters, especially for the baroclinic cases, with rotation only. There was no SSD seen for any of the cases, despite thousands of runs at various turnover times. In particular, this conclusion is reached in each of these cases: (i) changing R, Δt (Tables 2 and 3); (ii) using a simple uniform initial magnetic field, which experiences an early fast (exponential) decay, followed by a power law, with or without rotation; (iii) changing resolution; and (iv) changing τ_cool in the baroclinic cases (Table 4).

We found that by not including the cooling term in the entropy equation, the sine of the baroclinic angle, Eq. (15) takes peak values that are higher than one. The mean angle takes similar values of the ones found in Korpi et al. (1999). These results are in contrast with Achikanath Chirakkara et al. (2021) who found SSD in similar turbulent setups. In particular, even with a completely irrotational driving (implemented in the Fourier space) and without any rotation, they see that after ∼10³ t_turn there is an increase of magnetic energy, which saturates at about 1/1000 of the kinetic one. To compare them, simulations MB_0_128 and MB_2_128 were run to approximately 10⁴t_turn. We also ran simulations with high values for the forcing amplitude (up to ϕ₀ = 500, see Table 5 for non-rotating cases), which required higher values of viscosity and thermal diffusivities, leading to overall lower values of Re. However, we still could not observe any growth in the magnetic field.

As we wanted to investigate as much as possible the parameter space looking for possible dynamo action, we tried to push the code capability to its limits by increasing the values of Ω or decreasing diffusivites. In both cases, we could see that a growth of the field could be found only in cases that quickly became numerical unstable. In these cases, we were using either Ω = 10 or diffusivities of the order of ν = 10⁻⁶. As in these cases it was not possible to clearly establish growth and saturation phases, we did not consider them reliable and, hence, we did not include them in our analysis.

Summarizing, the irrotational forcing in combination with solid body rotation can produce vorticity. The chosen ideal EoS favors it more than the barotropic one. However, SSD is never activated, regardless of the rotation, EoS, and correlation length of the seed magnetic field, at least in the explored range of Rm up to a few hundred units.

3.4. Dynamo in the presence of a shear

A further effect to be investigated is the shearing flow in combination with this irrotational forcing. In Table 6 we report the numerical experiment performed with the sinusoidal shearing flow u^S introduced above in Sect. 2.3. Dynamo growth is always found, unless the shearing profile is rather weak (first three models in Table 6), or the forcing acts on a rather small length-scale (case M_0A020_W0.10_128), or if the forcing is not acting, meaning there is only the background flow (case M_0A020_F0_128). In terms of Reynolds numbers, we tentatively indicate a threshold value of Rm_crit ∼ 50, for the chosen setups.

The system evolves as illustrated in Fig. 5. There is an initial growth of the vorticity that after a few turnover times seems to saturate, along with the value of u_rms. The small-scale magnetic field initially decays. For about 1000 turnover times, the system slowly evolves by slightly increasing the vorticity and keeping both u_rms and b_rms almost constant. We then observe a sudden growth of vorticity followed by that of the kinetic energy, and then also of the magnetic energy. This process occurs with a strong enough shear amplitude, for a given a set of diffusivities and forcing parameters. We understand it as a vorticity dynamo produced by a Kelvin-Helmholtz instability, which develops in the system after it is perturbed beyond a certain threshold.

Fig. 5. Time evolution for the average turbulent kinetic energy density (solid blue and red lines), the average magnetic energy density (yellow line) and the root mean square of the turbulent vorticity ${\omega }_{\text{rms}}^{\prime }$ (dashed blue and red lines) for one HD setup and for an MHD set up with the same physical parameters, M_0A020_128.

In Fig. 5, we show the evolution of the average magnetic energy density, E_mag = ⟨B²⟩/2μ₀, the average turbulent kinetic energy density, $E_{\text{kin}}^{\prime }=\langle \rho u^{\prime 2}\rangle /2$, (i.e., neglecting the shearing contributions) and the root mean square of the turbulent vorticity ω′ = ∇ × u′, for one MHD and one HD representative runs. All averages are taken over the entire domain. We notice that $E_{\text{kin}}^{\prime }$ and ${\omega }_{\text{rms}}^{\prime }$ grow almost exactly the same way in the HD and MHD cases during the initial instability. Then the MHD run shows SSD shortly after the vorticity dynamo. When the magnetic field increases to a dynamical regime, the MHD run sees a decrease in $E_{\text{kin}}^{\prime }$ and ${\omega }_{\text{rms}}^{\prime }$, until we almost reach equipartition between the magnetic and the turbulent kinetic energies. If, instead, the shear velocity is included, it dominates the total kinetic energy by a factor of a few, the value depending on the shear amplitude.

Such an exponential amplification of vorticity, or vorticity dynamo, was initially predicted by Blackman & Chou (1997) who, however, showed how an helical forcing was needed. The effect we see in our simulations was predicted by Elperin et al. (2003) for non-helical homogeneous turbulence with a mean velocity shear. It was then observed in the HD case by Käpylä et al. (2009) in a setup with linear shear and plane wave forcing. Since the Kelvin-Helmholtz is a kinetic instability, the MHD case leads to SSD just after the vortical structures are formed: the magnetic fields are twisted and stretched by advection, until the magnetic energy is strong enough to provide feedback on the fluid, slightly reducing the mean kinetic energy density. At a later stage, also a large-scale component of the magnetic field is amplified: this can be interpreted as the shear-current effect described by Rogachevskii & Kleeorin (2003), where an electromotive force proportional to vorticity is produced.

In the presence of a differential velocity profile and once a sufficiently large magnetic field exists, we observe a relevant effect known as winding (e.g., Fujisawa 2015). Given such a shear profile as ours, $u_y^S$(z), the advective term of the induction equation provides an increase of the B_y component:

$$\begin{array}{c}\hfill \frac{\partial \mathbf{B}}{\partial t}=\mathrm{\nabla }\times (\mathbf{u}\times \mathbf{B})\to \frac{\partial B_y}{\partial t}=-B_z\frac{\partial u_y^S}{\partial z}.\end{array}$$(17)

This winding behavior leads to a linear growth of the magnetic field in the shearing direction which will stop when magnetic resistivity and/or dragging become important; namely: when the Lorentz force becomes relevant and the system non-linear. This can be seen for the model M_0A020 in the upper plot of Fig. 6, where during the time of dynamo growth each component of the magnetic field is equally important, until the non-linear problem sets in. Although our setup is more complex, as the potential forcing produces flow in all directions, this results in different growth of the field in the y direction during the latest stage of the field growth. In the middle and bottom panels of Fig. 6, we also show the root mean square of turbulent vorticity and turbulent velocity. They show a mild preferred contribution only in the y-component of the turbulent velocity (probably for the feedback of the strong magnetic field).

Fig. 6. Vector component evolution for the magnetic field, turbulent velocity and turbulent vorticity of the M_0A020 run. The y direction of the magnetic field shows an evolution different from the other components before and after, but not during, the dynamo phase. This is only slightly observed in the turbulent velocity field, while vorticity shows no preferential direction at all.

In Table 7, we show the different growth rates (in units $t_{\mathrm{turn}}^{-1}$) obtained for all dynamo runs. In Fig. 7 we also show the volume-integrated magnetic energy E_mag and the root mean square of the vorticity ω_rms evolution for some of the 256³ runs. We have found that the growth rate of the magnetic field is roughly proportional to the shear rate, similar to what was found by Käpylä et al. (2010). The dynamo growth seems to have saturated already in the 256³ run, but the vorticity growth seems to vary more. The 512³ run was stopped after the dynamo growth for computational reasons, thus we do not show the last slope. Baroclinicity does not change much the vorticity instability growth, but allows for a more rapid dynamo growth. For the runs with constant amplitude (A = 0.2) and increasing resolution there is little difference between the 256³ run and the 512³, which indicates numerical convergence. Comparing the 128³ dynamo runs with different forcing, we may notice that the dynamo and vorticity growth rates increase with both ϕ₀ and R, but when they exceed a certain value, both rates decrease. This is probably caused by the flow reaching a transonic regime and decreasing the growth efficiency (Haugen et al. 2004; Federrath et al. 2011; Schleicher et al. 2013).

Fig. 7. Evolution of the volume-integrated magnetic energy and the root mean square of the vorticity, ωrms, for some dynamo runs (as in the legend). The dynamo and vorticity growth rates increases with the shearing amplitude. The baroclinicity slightly enhances the dynamo growth rate even though it takes longer for the instability to start. The vorticity keeps growing after its exponential growth, during the magnetic dynamo phase. It then decays to lower values. The straight lines indicate the fits we obtained for the growth rates of the dynamo, r, of the vorticity dynamo, rω, and of the growth of the vorticity during the first part of its saturation phase, rω,sat. Their values are reported in Table 7.

We did not reach dynamo for the cases with a shear amplitude A ≲ 0.1, which can be translated as a critical Reynolds value of Re_crit ∼ 50, and a minimal rotational component of the flow of about 50% which in our models corresponds to k_ω/k_f, crit of 0.073. Although these specific values might depend on the resolution, this transition from non-dynamo to dynamo generating flows was not observed in the rigidly rotating cases.

Finally, in Fig. 8 we represent the spectra and snapshots of M_0A020 to give a more visual and detailed description of the instability. This figure summarizes the evolution of the main quantities during the three main evolutionary stages of our simulations.

Fig. 8. Spectra and snapshots over three faces of the domain, for the squares of velocity, u2, magnetic field, b2, and vorticity, ω2, of the run M_0A020. We display three different phases of the evolution of the system, evolving over time from left to right. Color bars have different ranges to allow a better visualization: for u2 the range is (0,0.2); for b2 the range is (0,1e−13) for the first snapshot and (0,0.015) for the others; for ω2 the range is (0,0.2) for the first snapshot and (0, 2.0) for the others, all in code units. The x, y and z components of the spectra correspond to the red, green and blue lines, respectively. Slope lines in the kinetic energy spectra are proportional to k−2 before the dynamo and get a bit steeper afterwards. Magnetic energy spectra take a typical Kazantsev k3/2 slope during the kinematic phase of the dynamo (Kazantsev 1968). We note that the scale of the kinetic energy and vorticity spectra (or enstrophy spectra) are kept constant, but the magnetic ones are not, in order to follow the magnetic field growth.

4. Possible astrophysical applications

The presented model is rather general and can therefore be applied to several astrophysical environments. In particular, the implemented forcing can be interpreted as the expansion waves coming from SNe, the main forcing in ISM. Although SNe tend to be highly supersonic and spherically asymmetric, the expansion waves under consideration here can be thought as the long-time expanding wave of the SN remnant.

We can interpret our box as a small cube with side 500 pc inside a galaxy (i.e., a resolution of about 2 pc). The forcing width of 1/10 of the box corresponds to 50 pc, slightly larger but of the same order of magnitude of some SN remnants. We can take the reference value for density of 10⁻²³ g cm⁻³, and the speed of sound of 10 km s⁻¹. Galactic rotation curves lead to a shearing amplitude of about 5 km s⁻¹ for a 500 pc radial distance, thus, it is on the same order of magnitude, A = 0.2, for the sinusoidal shear in our model.

The estimated SN rate of 2-3 SN per century per galaxy (Murphey et al. 2021) can be then translated to a 3–4 SN/Myr per (500 pc)³-box, assuming a galactic volume of the order of 10¹² pc³. The value Δt of 0.02 leads to a rate of ∼6 SN/My per (500 pc)³-box. As an example, for the M_0A020 run, these units lead to a u_rms of the order of 2 km s⁻² and a mean magnetic field of 10 μG, which is of the order of the estimated galactic one (e.g., Jansson & Farrar 2012; Brandenburg & Ntormousi 2023).

Therefore, the results can be indeed of direct application for ISM. More realistic simulations could consider density variations and non spherically symmetric expansion waves (which might lead easier to SSD). Finally, the shear here considered could be refined in order to replicate for instance the rotational curve of the galaxy or the shear corresponding to the spiral arms.

5. Conclusions

In this work, we investigate the possibility of having vorticity production and dynamo action driven by a purely curl-free forcing of the velocity field. We ran numerical models that employed an irrotational forcing, which does not lead to any growth of the vorticity neither in the HD nor in the MHD scenarios, if it is acting independently. However, a vortical flow is produced when the forcing interacts with a solid body rotation, in the presence of baroclinicity, or of a background shear. In the case of a rotating and baroclinic system we have not found any dynamo action within the explored parameter space, neither with an initial random seed magnetic field, nor with an initial uniform field. This result, therefore, seems not to be dependent on the topology of the initial field.

Instead, in the presence of a background, sinusoidal shearing velocity, we could observe an amplification of the initial seed magnetic field as a consequence of a dynamo process, both in the barotropic and in the baroclinic cases. The main novelty in the presented work is to show how such a dynamo phase occurs after the onset of a hydrodynamical instability driving an exponential growth of the vorticity. After the vorticity has grown, the magnetic energy is amplified too, approaching equipartition with the turbulent kinetic energy. This exponential amplification of vorticity occurs both in the purely HD case and in the MHD one, starting from small scales and then spreading up to the scale of the box. Similarly, the magnetic field is amplified first on small scales, but the inverse cascade makes the large scales dominate during the saturation phase of the simulations. Such large-scale field is related to the winding process and made visible after the equipartition. This winding enhances the field in the shearing direction, but the turbulent vorticity and velocity remain isotropic.

If the baroclinic term is at work, we noticed how the vorticity growth takes place at a slightly later time, with a growth rate similar to the barotropic case. Instead, the magnetic field grows at a bit faster rate.

Our results cannot firmly exclude that a dynamo is produced with irrotational forcing, rotation, and no shear, as, for example, seen by Seta & Federrath (2022) using a different Fourier-based forcing in a multi-phase MHD setup with the Flash code. We can, however, affirm that in our setup, such a putative dynamo could only be produced in a regime of higher Re than what tested in this work (up to a few hundreds). In any case, the critical threshold would therefore be much larger than what seen in the shear case (Re_crit ∼ 50).

Possible follow-up works in the context of dynamo triggered by an irrotational forcing could be focused on the exploration of higher Re and resolutions, as well as the inclusion of a stratified shearing medium. Such results are ultimately applicable to study galactic environments and other astrophysical scenarios, such as accretion disks and planetary or stellar atmospheres.

Acknowledgments

This work has been carried out within the framework of the doctoral program in Physics of the Universitat Autònoma de Barcelona and it is partially supported by the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. D.V. and A.E. are supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (ERC Starting Grant “IMAGINE” No. 948582, PI: DV). F.D.S. acknowledges support from a Marie Curie Action of the European Union (Grant agreement 101030103). The authors acknowledge support from “María de Maeztu” award to the Institut de Ciències de l’Espai (CEX2020-001058-M). We are grateful to Axel Brandenburg, Maarit Korpi-Lagg, Frederick Gent, Matthias Rheinhardt, Eva Ntormousi, Amit Seta and Claudia Soriano-Guerrero for fruitful discussions. We also want to acknowledge the whole Pencil Code community for support. A.E. and F.D.S. gratefully acknowledge NORDITA for hospitality during the Program “Magnetic field evolution in low density or strongly stratified plasmas” in May 2022, where part of this work was performed.

References

Appendix A: Spectra

Spectra of a continuous variable ($\stackrel{\sim }{u}$) can be calculated by numerically obtaining the 3D Fourier transform:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \stackrel{\sim }{\mathit{u}}(\mathit{k})=\frac{1}{{(2\pi )}^3}{\displaystyle {\int }_0^L{\int }_0^L{\int }_0^L\mathit{u}(\mathit{r})e^{-i\mathit{k}·\mathit{r}}dxdydz\approx }\\ \hfill \approx \frac{1}{N^3}{\displaystyle \sum _{p=0}^{N-1}\sum _{q=0}^{N-1}\sum _{r=0}^{N-1}u(x_p,y_q,z_r)e^{-i{k}_x{x}_p}{e}^{-i{k}_y{y}_q}{e}^{-i{k}_z{z}_r}}\end{array},\end{array}$$(A.1)

where N is the numbers of points along each box direction, L = 2π is the length of each box side, (x_p, y_q, z_r) indicate the discrete points of the mesh, and k_i are the set of integer wavenumbers related to each directions, each defined in the range [1, N/2]. The three dimensional (3D) spectrum is defined as:

$$\begin{array}{c}\hfill P(k)=\frac{1}{2}\stackrel{\sim }{\mathit{u}}(\mathit{k}){\stackrel{\sim }{\mathit{u}}}^{\ast }(\mathit{k})\text{where}k=\sqrt{k_x^2+k_y^2+k_z^2},\end{array}$$(A.2)

where the three-dimensional dependence of the spectrum on k is reduced to a function of the modulus k alone, by rebinning in the Fourier space.

This has been used for the velocity (kinetic energy spectra, neglecting the small density deviations from unity), magnetic field (magnetic energy spectra), and vorticity (enstrophy spectra).

Appendix B: Helmholtz decomposition

Helmholtz showed that any vector field F which vanishes suitably quickly at infinity can be decomposed into irrotational (longitudinal, purely divergent) and solenoidal (transverse) components. Therefore it can be expressed as the gradient and curl of a scalar and gradient potentials, respectively:

$$\begin{array}{c}\hfill \mathbf{F}={\mathbf{F}}_{\mathit{pot}}+{\mathbf{F}}_{\mathit{rot}}=-\mathrm{\nabla }\varphi +\mathrm{\nabla }\times \mathbf{A}.\end{array}$$(B.1)

Aided by the Pencil code available tools, we have used this decomposition to evaluate how much of the flow is irrotational, with no associated vorticity, and how much is rotational. To numerically obtain such decomposition we made use of the 3D Fourier transform to obtain the flow components in Fourier space: ${\stackrel{\sim }{\mathit{u}}}_{\mathit{pot}}(\mathit{k})=\stackrel{\sim }{\mathit{u}}(\mathit{k})·\mathit{k}$, ${\stackrel{\sim }{\mathit{u}}}_{\mathit{rot}}(\mathit{k})=\stackrel{\sim }{\mathit{u}}(\mathit{k})\times \mathit{k}$. If we inverse transform these functions, we can obtain the decomposition. We can check both the percentage of the rotational flow and the error of the numerical procedure by comparing the squared volume quantities with the relation:

$$\begin{array}{c}\hfill \mathit{u}(\mathit{r})={\mathit{u}}_{\mathit{pot}}(\mathit{r})+{\mathit{u}}_{\mathit{rot}}(\mathit{r})\to {⟨|\mathit{u}|}^2⟩=⟨|{\mathit{u}}_{\mathit{pot}}{|}^2⟩+⟨|{\mathit{u}}_{\mathit{rot}}{|}^2⟩.\end{array}$$(B.2)

In the analysis, we have chosen only one single snapshot (once saturation is reached) to perform the decomposition, for computational practical reasons. We verified that averaging among many snapshots does indeed lead to the percentage of rotational flow oscillating between values that vary by no more than 5 %.

Appendix C: Linear shear

We could opt for the use of a linear shear of the form u^S = (0, Sx, 0). With this choice, the yz-faces of the box cannot satisfy periodic boundary conditions, thus shearing-periodic boundary conditions are needed: the x direction is periodic with respect to positions in y that shift in time:

$$\begin{array}{c}\hfill f(-\frac{1}{2}{L}_x,y,z,t)=f(\frac{1}{2}{L}_x,y+L_{x}St,z,t),\end{array}$$

where f is each of the evolving fields. This condition is routinely used in numerical studies of shear flows in Cartesian geometry since it was introduced by Wisdom & Tremaine (1988) and Hawley et al. (1995). These boundary conditions are known to potentially produce spurious vorticity at the boundary, which, also due to the irrotational nature of our forcing, we clearly notice, as we discuss below.

We show the tested cases in Table C.1 (see Table C.2 for the diagnostics of some cases). For reasons of numerical stability, we set higher values of viscosity for the 128³ tests with the linear shear term. This leads us to a lower Re compared to the most of the other shear-less simulations. Regarding the magnetic field evolution (last column of Table C.1), in most cases we find a decay of magnetic energy, which depends on the shear parameter S. Only in cases with a low magnetic diffusivity (2 ⋅ 10⁻⁵), we observe some magnetic field exponential growth. However, we cannot draw any reliable conclusion about dynamo action, for the following reasons.

Looking at the trend between the vorticity indicator (at saturation) and the shear amplitude, S, k_ω/k_f = (0.108±0.001) + (0.151±0.005) S, we notice an nonphysical non-zero (and non-negligible) vorticity for vanishing shear, S, and we checked that it is not cured by increasing the resolution. This was already observed by Del Sordo & Brandenburg (2011) and attributed to numerical artifacts due to the interaction between the expansion waves laying on both sides of the shearing boundary.

A visual inspection of the results reveals indeed that almost the entire vorticity and magnetic field production were near the shearing boundaries. As an attempt to mitigate this artifact, we spatially constrained the expansion waves in the middle half of the box (i.e., as far away as possible from the shearing boundaries). This indeed resulted in a more spatially uniform and slower growth of vorticity and magnetic field, but we still could not keep the system from being affected by the aforementioned spurious contributions at the shearing boundaries after a certain time. Higher values of S or the addition of rotation seem to amplify some magnetic field in the middle of the box, but the boundary contribution was still dominant. This numerical issue has led us to use the more complex sinusoidal shearing profile which allows simple periodic boundary conditions, with no spurious effects.

All Tables

All Figures

	Fig. 1. Slice in the xy-plane showing the value of u2 at t/tturn = 0, 0.01 (top), 0.1, and 1 (bottom) for the non-rotating model H_0cW1 of Table 1.
In the text

	Fig. 2. Comparison of Re and Ma evolution in the HD cases, considering or not the rigid rotation (Ω = 0 or 2) and for different EoS (H vs. HB).
In the text

Fig. 3. Time-averaged kinetic spectra at saturation of some representative simulations of Tables 1–3. We notice how the results on large scales (small k) are very similar for all the models, excluding H_0cW1 (dash-dotted green line), which has a forcing with kf = 2. The spectrum with lowest power is H_0s (dashed green line), which has Δt = 1. Models marked with other colors explore variation of Ω and EoS. In all these models the magnetic field decays rapidly and kinetic spectra coincide with the corresponding HD runs of Fig. 2.

In the text

	Fig. 4. Magnetic spectra decay for the M_0 and M_2 models. The spectra evolution is only shown for M_0, and the first corresponds to less than a turnover time when it still shows the initial slope of k4. The arrow indicates the direction of temporal evolution. The last spectrum is approximately at turnover time 20 and it is shown for both models.
In the text

	Fig. 5. Time evolution for the average turbulent kinetic energy density (solid blue and red lines), the average magnetic energy density (yellow line) and the root mean square of the turbulent vorticity ${\omega }_{\text{rms}}^{\prime }$ (dashed blue and red lines) for one HD setup and for an MHD set up with the same physical parameters, M_0A020_128.
In the text

	Fig. 6. Vector component evolution for the magnetic field, turbulent velocity and turbulent vorticity of the M_0A020 run. The y direction of the magnetic field shows an evolution different from the other components before and after, but not during, the dynamo phase. This is only slightly observed in the turbulent velocity field, while vorticity shows no preferential direction at all.
In the text

Fig. 7. Evolution of the volume-integrated magnetic energy and the root mean square of the vorticity, ωrms, for some dynamo runs (as in the legend). The dynamo and vorticity growth rates increases with the shearing amplitude. The baroclinicity slightly enhances the dynamo growth rate even though it takes longer for the instability to start. The vorticity keeps growing after its exponential growth, during the magnetic dynamo phase. It then decays to lower values. The straight lines indicate the fits we obtained for the growth rates of the dynamo, r, of the vorticity dynamo, rω, and of the growth of the vorticity during the first part of its saturation phase, rω,sat. Their values are reported in Table 7.

In the text

Fig. 8. Spectra and snapshots over three faces of the domain, for the squares of velocity, u2, magnetic field, b2, and vorticity, ω2, of the run M_0A020. We display three different phases of the evolution of the system, evolving over time from left to right. Color bars have different ranges to allow a better visualization: for u2 the range is (0,0.2); for b2 the range is (0,1e−13) for the first snapshot and (0,0.015) for the others; for ω2 the range is (0,0.2) for the first snapshot and (0, 2.0) for the others, all in code units. The x, y and z components of the spectra correspond to the red, green and blue lines, respectively. Slope lines in the kinetic energy spectra are proportional to k−2 before the dynamo and get a bit steeper afterwards. Magnetic energy spectra take a typical Kazantsev k3/2 slope during the kinematic phase of the dynamo (Kazantsev 1968). We note that the scale of the kinetic energy and vorticity spectra (or enstrophy spectra) are kept constant, but the magnetic ones are not, in order to follow the magnetic field growth.

In the text