© The Authors 2024

1. Introduction

The occurrence of vortical flows is of high importance in astrophysics due to their close connection with turbulence. Vorticity is a way to characterize turbulent flows, which are ubiquitous in astrophysical contexts. However, the connection between turbulence and the processes responsible for the amplification of astrophysical magnetic fields has not been fully elucidated. Magnetic fields are amplified at all scales in various astrophysical environments, from planets to the extragalactic medium. While the range of scales, densities, and velocities is fairly wide, the instability responsible for the amplification of magnetic fields, known as dynamo, is thought to be active in all of these environments. Numerical models are commonly used to simulate many of these astrophysical contexts, from planets (e.g., Jones 2011) to galaxies (see, e.g., Brandenburg & Ntormousi 2023, for a recent review), and in each of them the turbulence is injected through a forcing mechanism. In this study we focus on two open questions, namely (a) whether the occurrence of dynamo instability depends on the forcing mechanism, and (b) the minimum ingredients needed to trigger a dynamo in a magnetohydrodynamic (MHD) turbulent medium.

In our attempt to address these two points, we expand on the work by Mee & Brandenburg (2006), Del Sordo & Brandenburg (2011), and Elias-López et al. (2023), who made use of a purely irrotational (i.e., curl-free) forcing of the velocity field to create a turbulent medium. The used forcing function mimics the occurrence of spherical expansion waves in an initially homogeneous medium, in a subsonic regime. The aforementioned studies found no vorticity amplification from the forcing alone in an isothermal purely hydrodynamic (HD) environment, and no dynamo amplification unless a background shearing profile is added to the system, independently of the equation of state. A similar result was obtained by Kahniashvili et al. (2012), who employed a curl-free forcing function to model inflationary scenarios in the early Universe. They found that the Lorentz force may produce some vorticity in isothermal models, but the magnetic field eventually dissipates without undergoing any kind of amplification. Dosopoulou et al. (2012) found analogous results in the context of magnetized and rotating cosmological models. Nevertheless, other approaches found vorticity and magnetic fields to be exponentially amplified when a purely irrotational forcing is added in the form of a stochastic function in Fourier space (e.g., Federrath et al. 2010; Achikanath Chirakkara et al. 2021; Seta & Federrath 2022), even in isothermal contexts, in the absence of large-scale contributions to the forcing, such as rotation or shear. As a consequence, it remains unclear which are the minimum ingredients needed to excite a dynamo. It is possible that the difference between these different studies resides either in the used forcing function or in the exploration of different regions of a fairly wide parameter space. This open question in particular motivated the work presented here.

In the present study we took as a starting point the work done by Elias-López et al. (2023), and we investigated the effect of varying the explosion width and the magnetic Prandtl number, as well as the effect of using a Newtonian cooling function, on the amplification of the vorticity and magnetic field. We also studied the difference between (a) forcing the turbulence by imposing an exactly irrotational forcing for the acceleration and (b) a variant that forces a locally fully potential velocity field and hence may include vorticity created by density fluctuations. The article is organized as follows. In Sect. 2 we present the numerical model used to perform the study, in Sect. 3 we describe the main results, and in Sect. 4 we discuss some of the conclusions of our work.

2. Model and numerical methods

2.1. MHD equations with rotation and shear

We employed the same numerical models as in Elias-López et al. (2023), which for completion we briefly explain below. We used the public 3D MHD code Pencil Code¹ (Pencil Code Collaboration 2021), which is a non-conservative, high-order, finite-difference code (sixth-order accurate in space and third-order Runge-Kutta in time). We solved the nonideal fully compressible MHD equations following an approach similar to what was done by Del Sordo & Brandenburg (2011), that is, either in a rigidly rotating frame, with angular velocity ω = Ωe_z, or with a differential velocity (shear) given by u^S = u_y^S(z)e_y, with u_y^S = Acos(kz), similar to Skoutnev et al. (2022), Käpylä et al. (2010), and Käpylä et al. (2009). In our models, z ranges from −π to π, so when we included a shearing profile we set k = 1. This allows simple periodic boundary conditions in the three directions. While in Elias-López et al. (2023) we explored the role played by the shear amplitude, in the present work we used A = 0.2 in all the shearing cases.

The set of equations employed consists of the continuity equation, the momentum equation, the entropy equation (which we only solve for the non-isothermal/baroclinic case), and the induction equation, respectively:

(1)

(2)

(3)

(4)

In Eqs. (1)–(4) ρ is the mass density; u(t) = (u^S + u′(t)) is the total velocity, which can be thought of as the sum of shearing and turbulent velocities; p is the pressure; B the magnetic field; A its vector potential (i.e., B = ∇ × A); J = (∇×B)/μ₀ is the electrical current density (where μ₀ is the vacuum permeability); F_visc = ρ⁻¹∇⋅(2ρνS), where the traceless rate of strain tensor S has components S_ij = (1/2)(u_i, j + u_j, i − (1/3)δ_ij∇⋅u); f the expansion wave forcing (see below for their definitions); χ is the thermal diffusivity; η is the magnetic diffusivity, c_s0 is the initial, uniform sound speed (proportional to the initial temperature) and τ_cool is the cooling term timescale, introduced to avoid an indefinite heating. The advective derivative operator is D/Dt := ∂/∂t + u ⋅ ∇. The differential velocity profile is imposed directly on the y component of the total velocity by an amount proportional to the difference of velocity and the profile itself, so u_y does not deviate much from the differential shearing profile:

where we fix τ_S = 1 (as for τ_cool for the entropy, the smaller the value, the more effective is the keeping u_y close to u_y^S).

In order to close the system of equations, we consider two types of equation of state (EoS): (1) a simple barotropic EoS p(ρ) = c_s²ρ, where we fix the value of the sound speed c_s = 1, or (2) an ideal EoS, also dubbed the 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_s² = (γ − 1)c_pT, 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.

The forcing is applied as an acceleration in the momentum equation and we explore two different ways: either the acceleration itself is irrotational, as in Elias-López et al. (2023), or the force itself is irrotational:

(5)

(6)

where x_f(t) is the randomly changed expansion wave center, ρ₀ is the initial, uniform density, is the normalization factor, c_s0 is the initial sound speed, R is the radius of the Gaussian, Δt is the time interval after which a new expansion wave is forced in a new position (and it can be as short as the time step), and ϕ₀ controls the overall forcing amplitude and has dimensions of velocity squared. In either case, the associated forcing representative wavenumber is thus k_f = 2/R. A spherical symmetric explosion in the interstellar medium (ISM) is arguably better modeled by the second choice, since the expansion takes ISM density fluctuations into account. We note that, in any case, this is a simplification, since important deviation from spherical symmetry are expected in detailed 3D supernova (SN) explosion simulations (Reichert et al. 2023). However, in this work we aim at exploring the difference between the two kinds of forcing, within the spherical symmetric assumption, to understand if the density fluctuations alone can be responsible of vorticity-triggered dynamo action.

The simulation domain consists of a uniform, cubic grid mesh [ − π, π]³, with triply periodic boundary conditions. We also studied resolution convergence varying from 32³ up to 256³ mesh points, and we find that 256³ mesh points are enough to assess our problem. We ran some models at an even higher resolution, 512³, to double-check the validity of the results obtained at lower resolutions. However, simulations at this resolution are computationally expensive, so we limited their use to double check numerical convergence. The chosen velocity profile allows simple periodic boundary conditions in the z direction if k is an integer. We adopted nondimensional variables by measuring speed in units of the initial sound speed, c_s0, and length in units of 1/k₁, where k₁ is the smallest wave number in the periodic domain, implying that the nondimensional size of the domain is (2π)³.

As the initial conditions: pressure and density are set constant and with value of 1 throughout the box (making ρ₀ = 1 for all times), and so are entropy and temperature in the baroclinic case; the fluid is initially at rest so the flow is described by u = 0; the initial magnetic field is a weak seed randomly generated with an amplitude of 10⁻⁶ in code units, uncorrelated at each point for the three components, corresponding to a E_k ∼ k⁴ power law, as reported by Mee & Brandenburg (2006).

2.2. Diagnostics

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 looked at the root mean square of the total velocity, u_rms. We in turn used this to define the fundamental timescale of our problem, which we call turnover time, as

(7)

The turnover time can be understood as the average time for the fluid to cross an explosion width. In the cases where we employed the shearing profile, we also used a similar definition with the shearing wavelength, k, and amplitude, A:

(8)

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

which are the Reynolds number, magnetic Reynolds number, vorticity Reynolds number, Mach number, magnetic Prandtl number, and a measure of vortical wavelength, respectively.

2.3. Vorticity equation

To study the evolution of the vorticity (ω = ∇ × u), its sources or dissipative terms, we looked at some diagnostic quantities derived from the terms of the vorticity evolution equation:

(9)

Here the first term on the right-hand side is analogous to the advective term ∇ × (u × B) in the induction equation, the second term represents the viscous forces acting on the system, the third is the baroclinic term, related to the EoS, the forth is the effect of the Lorentz force, the fifth appears if the system is rotating, the sixth is due to the effect of the implemented forcing and the last one regards the sinusoidal shearing profile.

Taking the dot product with ω, integrating over the volume, and using the vector identities (∇ × a)⋅b = ∇ ⋅ (a × b)+a ⋅ (∇ × b) and ∇²a = ∇(∇ ⋅ a)−∇ × (∇ × a), we obtain

(10)

where q = ∇ × ω. We did not evolve ω itself; thus, all the quantities are obtained from the evolution of u. We used these diagnostic magnitudes to discriminate which are the most relevant vorticity amplification or destruction terms. Additionally, we can address whether the total sum of terms is similar to the numerical time derivative of ⟨ω²⟩.

3. Results

The main result of our study is that, throughout our fairly wide exploration of the parameter space in terms of forcing scales (see Sects. 3.1 and 3.2), magnetic Prandtl number (Sect. 3.3), cooling times (Sect. 3.4), and Reynolds numbers, we do not obtain an HD or MHD instability unless a background shearing flow is imposed. This result holds both with a forcing acting on the momentum, as in Eq. (6), and on the velocity field alone, as in Eq. (5). The models including a shearing profile instead develop an exponential increase in vorticity, followed by an exponential increase in the magnetic field, unless the scale of the forcing is too small, as explained in Sect. 3.2. In Sect. 2.3 we describe the contribution of vorticity source terms from Eq. (9). In Appendix A we show all tabulated runs and diagnostics described in this work.

3.1. Dependence on the forcing scale without shear

We explored the role played by the forcing scale R, which can be interpreted as the physical scale of the interaction between a SN explosion and the ISM. As a reference in this sense, the typical sizes of a SN remnant are typically 20–50 pc (Franchetti et al. 2012; Asvarov 2014). The size of the box can be converted into physical scales if the shear is included: for the values that we employed here, the box corresponds to about 500 pc (Elias-López et al. 2023). Anyways, we stress the aim of our box simulations is to investigate basic aspects of vorticity dynamo, rather than a detailed model of the effect of SN on the ISM.

We note that when we changed R, we changed ϕ₀ accordingly in order to reach a similar u_rms. This allowed us to compare simulations with different Reynolds numbers (which scale with R), different energy injection scales, and different filling factors within the box. This is especially relevant when the shear is included, since it represents another physical scale and the vorticity production and dynamo depend on both (see Sect. 3.2).

In the left panel of Fig. 1 we plot the temporal evolution of ω_rms in terms of turnover time for different, representative runs without rotation, with isothermal conditions. Independently of the forcing scale R, vorticity reaches a steady state after less than 15 turnover times, with the exception of the smallest value of R = 0.10, which takes less than 5 turnover times. In the model with a forcing scale of R = 0.10, we see the development of local transonic flows as a consequence of the attempt to reach values of the kinetic energy similar to the cases with larger forcing widths. This is accomplished by increasing the parameter ϕ₀ in Eq. (5) and it leads to a different behavior of these models. The mean value of vorticity is observed to decrease with R, while its fluctuations do increase.

Fig. 1. Time evolution for ωrms (left) and time-averaged kinetic spectra at saturation (right) for nonrotating isothermal runs with different explosion widths and similar total energy (see Table A.1). Notice that the peaks in k are close to the corresponding forcing wavenumber kf = 2/R.

In the right panel of Fig. 1, we show the different kinetic spectra obtained from such runs. The corresponding forcing widths change the forcing wavelength k_f, thus changing the inertial range interval, as ν is kept constant. The slope of −2 seems to be independent of R, as already seen by Mee & Brandenburg (2006). No dynamo was observed in these runs despite having reached Rm larger than 200. When rotation is added, the results are the similar to Elias-López et al. (2023): we obtain steeper slopes but a similar inertial range behavior as the nonrotating cases here plotted.

This can be seen more clearly in Fig. 2, where we plot the average of various diagnostics calculated during the saturated stage for the isothermal and baroclinic models, with or without rotation, using the forcing f_acc. We note that for R = 0.1 the non-isothermal runs are not shown because they reach locally supersonic flows and remain numerically stable for too few time steps. Despite attempting to select a ϕ₀ that makes u_rms approximately independent of R, we succeeded in doing so only for the isothermal nonrotating case, while for the other non-isothermal or rotating cases u_rms mildly depends on R. Therefore, the more important dimensionless quantity is k_ω/k_f, as it is a measure of vorticity normalized with velocity and forcing wavelength. Figure 2 also illustrates many other features. First, we see how the Reynolds number increases with width for all cases, but no instability is found nevertheless. We notice that the nonrotating cases seem to have a lower overall ω, but it tends to growth with smaller R. In contrast, for the rotating runs there is a maximum at R = 0.4 (i.e., when expansion waves are about a fifth of the simulation domain, k_f = 5). We also see how isothermal rotating case show the highest values for ω_rms, u_rms, and Re, but in k_ω/k_f they are closely matched by the non-isothermal rotating with the addition of a the thermal cooling time. As expected, the cooling term creates an additional source of dissipation, which we had to compensate for by doubling the values of ϕ₀ (for τ_cool = 0.1).

Fig. 2. Different diagnostic quantities, ωrms, urms, kω/kf, and Re (from top to bottom), as a function of explosion width, R, for the forcing facc. These runs have Pm set to 1, i.e., Rm = Re, with a grid size of 2563 and no dynamo present. The data are presented in Tables A.1 and A.2.

In Fig. 3 we plot the same quantities as in the left panel of Fig. 2, comparing runs with different resolutions and keeping R constant with a value of 0.5. Most quantities seem to be resolution independent even when we move to resolution low enough to erase a good part of the inertial ranges. We observe that u_rms is more or less resolution independent for all except two cases: (i) the isothermal nonrotating case, and (ii) the isothermal rotating case. The first case is compatible with the idea of vorticity being created solely by numerical sources.

Fig. 3. Same magnitudes as in Fig. 2 but now for runs with R = 0.5, i.e., kf = 4, different resolutions, and varying the type of forcing.

As expected the contribution of rotation to ω is much greater than the baroclinc one in non-isothermal models, but this difference diminishes at small values of the forcing scale R. No dynamo is found either for the cases with rotation and/or non-isothermality with an ideal gas law and different cooling times. Thus, the results of Elias-López et al. (2023) still hold when k_f grows nearly up to 1, and for Rm values of several hundred.

3.2. Dependence on the forcing scale in the presence of shear

The general behavior with the presence of shear is similar to that found in Elias-López et al. (2023): an HD instability develops first with an exponential growth of ω_rms, which is then closely followed by a magnetic instability leading to an exponential amplification of b_rms. The complete set of runs with the instability is in Table A.3 alongside with their diagnostics. After the linear phase of the dynamo a winding phenomena is seen for all cases, independently of R, Prandtl number and resolution. During this process B_y is further amplified in the shearing direction in a linear way, by winding.

In Fig. 4 we plot the evolution of E_mag and ω_rms for the isothermal runs with different R and using the f_acc forcing. The only model that does not develop a dynamo is that of the smallest forcing scale (i.e., R = 0.1), and this is also true for the baroclinic case. Our interpretation is that, in this case, Re is subcritical, and this allowed us to estimate a critical value of Re ∼ 40 for the vorticity instability to take place. We find instead a critical magnetic Reynolds number of slightly less than 20 (see the tabulated values) for the dynamo instability.

Fig. 4. Time evolution (in units of tturn) of vorticity and magnetic energy of runs with the shearing profile and different R (see Table A.3). If we use the shearing timescales, the magnetic instability is between 70 and 200 tshear with more similar growth rates (see Fig. 5). In both cases, R = 1.00 and R = 1.50 take the least amount of time to reach the instability. Dashed lines represent the exponential fits for Emag during dynamo growth (top panel), ωrms during vorticity growth, and ωrms during dynamo growth (bottom panel).

In Fig. 5 we show the growth rates r as a function of R for both isothermal and baroclinic cases. We show the values of r calculated using different time units: either t_shear or t_turn. When choosing t_shear as the time unit, r approximately displays a constant behavior, while when choosing t_turn the growth rates get vary considerably with R, as smaller expansion waves lead to much slower growths.

Fig. 5. Vorticity and magnetic energy growth rates as a function of explosion width, R. The blue lines represent the isothermal models and the green ones the baroclinic cases. The dash-dotted lines correspond to growth rates in time units of tshear and solid lines in forcing turnover times, tturn. Both time units are tabulated in Appendix A. The points marked with crosses are the corresponding dynamo runs with fmom.

We observe that both magnetic and kinetic helicities grow in the dynamo cases, and start oscillating when b_rms and ω_rms saturate. These oscillation resemble those seen in other instabilities such as the Tayler instability (e.g., Guerrero et al. 2019; Stefani et al. 2021; Monteiro et al. 2023).

3.3. Magnetic Prandtl number dependence

The magnetic Prandtl number (Pm) controls the scales at which the kinetic and magnetic energy cascades are truncated respectively by viscosity and resistivity. We varied Pm in our simulations so as to observe whether a difference in these truncation scales leads to a different behavior for the amplification of vorticity and magnetic field. In the ISM, Pm is usually much larger than 1 (Ferrière 2020), which is different from the range explored in our models. We also explored a range of Pm slightly smaller than the unity, which is more typical of planetary environments.

We observe only a weak dependence on Pm for the models that do not develop a dynamo instability. In the isothermal case, we varied Pm from 0.25 up to 4 and see that Re_ω either increases with Pm in the absence of rotation or slightly decrease with Pr when rotation is added (for more details see Table 1). We recall here that Re_ω is a dimensionless number introduced in Elias-López et al. (2023) with the aim of quantifying the vorticity. In the baroclinic case, Re_ω slightly decreases when Pm = 4 compared to the case of Pm = 0.25, independently on the presence of rotation. However, in the range of the explored values of Pm, we always report a decrease in the initial magnetic field.

Models with shear, conversely, develop a dynamo instability unless Pm = 4 or above (see the tabulated values). We interpret this as a consequence of the lack of vorticity instability that does not develop when the physical viscosity increases above a certain value. The growth rates for the magnetic field increase with Pm, even in the cases where η is kept constant and ν increased, of course up until the point where viscosity does not allow the hydro instability. Conversely, the growth of vorticity is approximately constant in the explored range.

We find that for Pm = 0.1 the vorticity is amplified, but, differently from the other cases, this instability is not followed by an exponential amplification of the magnetic field. This can be seen as a consequence of having a Rm below the critical value, since, in general, it is also possible to excite a dynamo at a Pm value of less than 0.1 (e.g., Warnecke et al. 2023). Our results are therefore consistent with what found in literature.

3.4. Dependence on the cooling time

When the isothermal condition is relaxed, we let the temperature evolve according to Eq. (3), where we use a Newtonian cooling term, regulated by the timescale τ_cool. The use of the cooling function leads to spectra cut at short wavelengths in models that does not develop instabilities. This is shown in Fig. 6 (compare the two purple lines). Differently, in the presence of shear and hence after a dynamo is excited, all the spectra recover their small-scale contribution and show again a wide dynamical range decreasing in a k⁻² fashion down to a dissipation scale, regardless of the cooling term (compare the orange curves). The use of the cooling term slightly diminishes the total amount of energy, but does not change notably the shape of the spectral distribution. However, we observe that in the presence of this cooling term the instability kicks in at much earlier times, independently of the value of τ_cool, at least within the explored range. This can be seen as a quicker injection of vorticity in the system due to the cooling function. We also observe that the average angle between ∇T and ∇s slightly increases when a cooling function is used, hence leading to a larger contribution of the baroclinic term in seeding the vorticity.

Fig. 6. Kinetic spectra for models with R = 0.5 and an optional shear and cooling term. The models with shear are shown after dynamo growth.

Another possible interpretation is to invoke an effect similar to what observed by Rädler et al. (2011). Irrotationally forced flows present peculiarities such as the possibility of having a negative magnetic diffusivity contribution from turbulent flows, especially at low Reynolds numbers, as shown analytically by Krause & Raedler (1980), Rädler & Rheinhardt (2007) and Rädler et al. (2011). The contributions to the diffusivity coming from the turbulence may affect the occurrence of dynamo instability even when one moves toward higher values of the magnetic Reynolds number, a regime that is closer to that of astrophysical bodies. Rädler et al. (2011) found, with mean-field approaches based on the second-order correlation approximation (SOCA), that a negative contribution to the magnetic diffusivity can come from the presence of turbulence in irrotational flows in the case of small Péclet numbers Pe = uL/k, where u, L, and k are the typical velocity, length scale, and diffusivity of a system. In our case, the presence of a cooling time introduces an additional diffusion term for thermodynamical quantities, which results in Rm being smaller than it would be in the absence of cooling. However, the SOCA approximation is valid only for small magnetic Reynolds numbers, which is a condition that is not satisfied in our models.

Although our cooling function is not meant to model any specific astrophysical environment, we can attempt a comparison with typical values of the cooling in the ISM. Using the hydrogen cooling function Λ (e.g., Sutherland & Dopita 1993) and assuming a temperature of 10⁴ K,

If, as done by Elias-López et al. (2023), we considered a time unit of 8 Myr, then the cooling time τ_cool should be on the order of 0.001. In our models we can reach down to τ ∼ 0.01 for numerical stability reasons. This is one order of magnitude lower than typical ISM values (i.e., lower cooling rates and higher temperature differences than a one-phase ISM).

3.5. Vorticity source terms

From Eq. (10) we can evaluate which are the most relevant terms in the vorticity equation for the different runs. We find that in the isothermal nonrotating runs there are no positive terms comparable to the viscous ones. This fact along side the decrease in vorticity growth with resolution (see Fig. 3), makes us think that in such cases only numerical diffusive sources are in action, in agreement with what we already observed in Mee & Brandenburg (2006) and Elias-López et al. (2023).

In Fig. 7 we can see the time series of the different source terms as a function of time, for a representative non-isothermal case. All terms are very small. The baroclinic term dominates as the most positive contribution. As vorticity grows, the turbulent contribution ⟨(u × ω)⋅q⟩ gains importance, and the sum of both are counteracted by viscous forces, so that ω_rms saturates. From the viscous contributions, only ν⟨q²⟩ is relevant, while 2ν⟨S∇lnρ ⋅ q⟩ is more than one order of magnitude lower. This last statement holds true for all the isothermal, non-isothermal and rotating, nonrotating cases. Obviously, the Lorentz term is irrelevant in the cases without dynamo, orders of magnitude below by comparison.

Fig. 7. Time evolution for ωrms2/2 (left axis and gray line) and vorticity growth terms (right axis and various colored lines for each term) for a non-isothermal run with R = 0.5 and forcing in acceleration form (and by construction its contribution to ⟨∂tω2/2⟩ is zero). This plot zooms in on the beginning of the temporal evolution of all the terms until t ≃ 130 tturn. However, the saturation regime does not present any relevant changes in time for more than 2000 tturn, and no instability is reached. Both the 2ν⟨S∇lnρ ⋅ q⟩ and ⟨(j × b/ρ)⋅q⟩ terms are negligible and very close to 0, so this makes the total contribution of viscous forces, i.e., ⟨Fvisc ⋅ q⟩, overlap with ν⟨q2⟩.

When the forcing is exactly irrotational, its corresponding term does not contribute to vorticity generation. But when it is applied in its second form (not exactly irrotational due to density fluctuations), there is indeed a small vorticity growth that leads to a similar behavior in strength and shape compared to the baroclinic source term. When applying this type of forcing in the non-isothermal runs, the forcing vorticity growth overtakes the baroclinic term to such an extent that this latter becomes negative.

When the rotation is included, the vorticity generation is more relevant. In this case, the Coriolis source creates an amount of vorticity that is later counteracted by viscous terms, so that the steady state is reached. In Fig. 8 we show the same plot as in Fig. 7 but having added rotation. We can see that in the beginning the rotation term, 2Ω⟨(e_z × u)⋅q⟩, has a big positive spike that leads to the initial growth of vorticity. We note that it is larger than the baroclinic term in Fig. 7 by more than one order of magnitude and oscillates substantially even becoming negative at times, leading to an overall a noisier ω_rms. In the beginning of the run, the rotation contribution is mostly positive, and, when ω_rms saturates, the term is slightly positive in average, and of the same order of the viscous main contribution. Thus, in the presence of rotation, all other source terms become much less important.

Fig. 8. Time evolution for ωrms2/2 and vorticity growth terms for a non-isothermal rotating run with R = 0.5 and forcing in acceleration form. The notation is the same as in Fig. 7, and similarly 2ν⟨S∇lnρ ⋅ q⟩ and ⟨(j × b/ρ)⋅q⟩ overlap near 0, also making ⟨Fvisc ⋅ q⟩ overlap with ν⟨q2⟩.

For the cases with the shearing profile, both the HD and MHD instabilities make the contributions change substantially. In Fig. 9 we show all the relevant terms for an isothermal shearing run. As before, within the viscous forces, ν⟨q²⟩ still dominates over 2ν⟨S∇lnρ ⋅ q⟩, but the latter becomes more relevant in comparison to the cases described above.

Fig. 9. Time evolution for ωrms2/2 and vorticity growth terms for a isothermal run with the presence of shear, hence leading to instability, and with R = 0.5 and the exactly irrotational forcing. The notation is the same as in Fig. 7.

The background shearing profiles increases ω_rms up to a certain values in a very few time steps. This value is kept approximatively until the vorticity instability kicks in (for example, until t ≃ 400t_turn in Fig. 9). Afterward, when vorticity is amplified, the advective term (which includes both shear and turbulence), ⟨(u × ω)⋅q⟩, brings the main positive contribution, as expected. The viscous forces are not enough to counteract this completely, thus leading to growth in ω_rms.

The Lorentz term is negligible up to the point where dynamo starts. Then there is a brief time when it becomes slightly negative exactly when vorticity starts growing exponentially but still the dynamo has not kicked in. This behavior was also observed by Seta & Federrath (2022). When the kinetic phase of the dynamo starts, the Lorentz term increases but as the Lorentz forces act against the flow, ⟨(u × ω)⋅q⟩ decreases more than ⟨(j × b/ρ)⋅q⟩ increases. This leads to a negative overall contribution and a decrease in vorticity, which later stabilizes at the end of the kinetic phase, and to the amplification of B_y by winding. In all dynamo runs the Lorentz term always ends up surpassing the advective term ⟨(u × ω)⋅q⟩.

4. Conclusion

This work, a continuation of Elias-López et al. (2023), studied the vorticity and dynamo instability in the presence of irrotationally forced turbulence. We first explored the role played by the scale on which the irrotational forcing is acting. We find that, independent of the scale of this forcing, no dynamo instability develops for systems that do not include any shear. Also, while the root mean square values of vorticity weakly depend on the forcing scale, in no case do we observe an exponential amplification of vorticity. When shear is added to the picture, the vorticity is always exponentially amplified after a transient time if the kinematic viscosity is low enough, with the exception of the case of a very small forcing scale, which does not lead to any growth over thousands of turnover times. Then it is followed by a dynamo instability if the magnetic diffusivity is not too high. By analyzing kinematic spectra, we see how the typical scale of the system is provided by the forcing scale of the turbulence before the vorticity is amplified and, conversely, by the scale of the shear, in the saturation phase. Based on that, we observe that the growth rate of the dynamo depends on the scale of the expansion waves if the growth rate is calculated using the turnover time as the time unit, but it exhibits no R dependence if the time unit is the shear timescale. The scale of the forcing also sets the time needed for the instability to develop. Models with a larger forcing scale amplify vorticity and magnetic fields after shorter times. The only models that immediately develop the instability are those that include both the baroclinic term and a cooling function. We observe an increase of one order of magnitude for the growth rate with the magnetic Prandtl number, which we varied between 0.1 and 10. We extrapolated a critical value for the magnetic Reynolds number of slightly less than 20.

With these results in hand, we conclude that the presence of shear remains the basic ingredient for triggering a dynamo instability when subsonic turbulence is driven by spherical expansion waves. Future work will have to take into account turbulence forced on more than one scale at the same time as well as the role played by plane waves, before moving toward more complex models that include density stratification and also take shocks and supersonic flows into account.

This work both confirms previous results and extends them to a very wide, previously unexplored range of parameters and using two different kinds of spherically symmetric forcing: one that is curl-free and one that is not and takes density fluctuations into account. With the exploration of parameters, we definitely confirm that the shear is a key ingredient for a vorticity dynamo, while rigid rotation alone is not able to trigger it, regardless of the magnitude of the Reynolds number, the forcing, and the EoS.

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. DV and AE 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). FDS 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, Matthias Rheinhardt, Eva Ntormousi, Federico Stasyszyn, Andrea Pallottini, Amit Seta and Claudia Soriano-Guerrero for fruitful discussions. We also want to acknowledge the whole Pencil Code community for support. AE gratefully acknowledges Scuola Normale Superiore for hospitality during October 2023. We acknowledge the use the SCAYLE supercomputer of the Spanish Supercomputing Network, via project RES/BSC Call AECT-2023-2-0034 (PI FDS).

References

Appendix A: Simulations with sinusoidal shear

The following tables contain all the simulations included in this work, as well as the most relevant runs in Elias-López et al. (2023) (marked in violet). We note that the runs from the past article do not have all diagnostics. For all models we used B₀ = 10⁻⁶ as the seed field, Δt = 0.02 as interval between two different explosions, c_s0² = 1 in the case of non-isothermal runs and an amplitude of the shearing profile A= 0.2 in the shearing cases.

All Tables

All Figures

	Fig. 1. Time evolution for ωrms (left) and time-averaged kinetic spectra at saturation (right) for nonrotating isothermal runs with different explosion widths and similar total energy (see Table A.1). Notice that the peaks in k are close to the corresponding forcing wavenumber kf = 2/R.
In the text

	Fig. 2. Different diagnostic quantities, ωrms, urms, kω/kf, and Re (from top to bottom), as a function of explosion width, R, for the forcing facc. These runs have Pm set to 1, i.e., Rm = Re, with a grid size of 2563 and no dynamo present. The data are presented in Tables A.1 and A.2.
In the text

	Fig. 3. Same magnitudes as in Fig. 2 but now for runs with R = 0.5, i.e., kf = 4, different resolutions, and varying the type of forcing.
In the text

Fig. 4. Time evolution (in units of tturn) of vorticity and magnetic energy of runs with the shearing profile and different R (see Table A.3). If we use the shearing timescales, the magnetic instability is between 70 and 200 tshear with more similar growth rates (see Fig. 5). In both cases, R = 1.00 and R = 1.50 take the least amount of time to reach the instability. Dashed lines represent the exponential fits for Emag during dynamo growth (top panel), ωrms during vorticity growth, and ωrms during dynamo growth (bottom panel).

In the text

Fig. 5. Vorticity and magnetic energy growth rates as a function of explosion width, R. The blue lines represent the isothermal models and the green ones the baroclinic cases. The dash-dotted lines correspond to growth rates in time units of tshear and solid lines in forcing turnover times, tturn. Both time units are tabulated in Appendix A. The points marked with crosses are the corresponding dynamo runs with fmom.

In the text

	Fig. 6. Kinetic spectra for models with R = 0.5 and an optional shear and cooling term. The models with shear are shown after dynamo growth.
In the text

Fig. 7. Time evolution for ωrms2/2 (left axis and gray line) and vorticity growth terms (right axis and various colored lines for each term) for a non-isothermal run with R = 0.5 and forcing in acceleration form (and by construction its contribution to ⟨∂tω2/2⟩ is zero). This plot zooms in on the beginning of the temporal evolution of all the terms until t ≃ 130 tturn. However, the saturation regime does not present any relevant changes in time for more than 2000 tturn, and no instability is reached. Both the 2ν⟨S∇lnρ ⋅ q⟩ and ⟨(j × b/ρ)⋅q⟩ terms are negligible and very close to 0, so this makes the total contribution of viscous forces, i.e., ⟨Fvisc ⋅ q⟩, overlap with ν⟨q2⟩.

In the text

	Fig. 8. Time evolution for ωrms2/2 and vorticity growth terms for a non-isothermal rotating run with R = 0.5 and forcing in acceleration form. The notation is the same as in Fig. 7, and similarly 2ν⟨S∇lnρ ⋅ q⟩ and ⟨(j × b/ρ)⋅q⟩ overlap near 0, also making ⟨Fvisc ⋅ q⟩ overlap with ν⟨q2⟩.
In the text

	Fig. 9. Time evolution for ωrms2/2 and vorticity growth terms for a isothermal run with the presence of shear, hence leading to instability, and with R = 0.5 and the exactly irrotational forcing. The notation is the same as in Fig. 7.
In the text