Issue 
A&A
Volume 690, October 2024



Article Number  A77  
Number of page(s)  13  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202450398  
Published online  01 October 2024 
Vorticity and magnetic dynamo from subsonic expansion waves
II. Dependence on the magnetic Prandtl number, forcing scale, and cooling time
^{1}
Institut de Ciències de I’Espai (ICECSIC), Campus UAB, Carrer de Can Magrans s/n, 08193 Cerdanyola del Vallès, Barcelona, Catalonia, Spain
^{2}
Institut d’Estudis Espacials de Catalunya (IEEC), 08860 Castelldefels, Barcelona, Catalonia, Spain
^{3}
INAF, Osservatorio Astrofisico di Catania, via Santa Sofia, 78, Catania, Italy
^{4}
Institute of Applied Computing & Community Code (IAC3), University of the Balearic Islands, Palma 07122, Spain
Received:
16
April
2024
Accepted:
12
July
2024
Context. The amplification of astrophysical magnetic fields takes place via dynamo instability in turbulent environments. Vorticity is usually present in any dynamo, but its role is not yet fully understood.
Aims. This work is an extension of previous research on the effect of an irrotational subsonic forcing on a magnetized medium in the presence of rotation or a differential velocity profile. We aim to explore a wider parameter space in terms of Reynolds numbers, the magnetic Prandtl number, the forcing scale, and the cooling timescale in a Newtonian cooling. We studied the effect of imposing that either the acceleration or the velocity forcing function be curlfree and evaluated the terms responsible for the evolution vorticity.
Methods. We used direct numerical simulations to solve the fully compressible, resistive magnetohydrodynamic equations with the Pencil Code. We studied both isothermal and nonisothermal regimes and addressed the relative importance of different vorticity source terms.
Results. We report no smallscale dynamo for the models that do not include shear. We find a hydro instability, followed by a magnetic one, when a shearing velocity profile is applied. The vorticity production is found to be numerical in the purely irrotational acceleration case. Nonisothermality, rotation, shear, and densitydependent forcing, when included, contribute to increasing the vorticity.
Conclusions. As in our previous study, we find that turbulence driven by subsonic expansion waves can amplify the vorticity and magnetic field only in the presence of a background shearing profile. The presence of a cooling function makes the instability occur on a shorter timescale. We estimate critical Reynolds and magnetic Reynolds numbers of 40 and 20, respectively.
Key words: ISM: magnetic fields / galaxies: evolution / intergalactic medium / galaxies: magnetic fields
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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 EliasLópez et al. (2023), who made use of a purely irrotational (i.e., curlfree) 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 curlfree 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 largescale 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 EliasLó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 EliasLópez et al. (2023), which for completion we briefly explain below. We used the public 3D MHD code Pencil Code^{1} (Pencil Code Collaboration 2021), which is a nonconservative, highorder, finitedifference code (sixthorder accurate in space and thirdorder RungeKutta 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 EliasLó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 nonisothermal/baroclinic case), and the induction equation, respectively:
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)/μ_{0} is the electrical current density (where μ_{0} is the vacuum permeability); F_{visc} = ρ^{−1}∇⋅(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}^{2}ρ, 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_{g}T, with R_{g} the specific gas constant and T the temperature; in this case, the sound speed squared is c_{s}^{2} = (γ − 1)c_{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.
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 EliasLópez et al. (2023), or the force itself is irrotational:
where x_{f}(t) is the randomly changed expansion wave center, ρ_{0} 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 ϕ_{0} 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 vorticitytriggered dynamo action.
The simulation domain consists of a uniform, cubic grid mesh [ − π, π]^{3}, with triply periodic boundary conditions. We also studied resolution convergence varying from 32^{3} up to 256^{3} mesh points, and we find that 256^{3} mesh points are enough to assess our problem. We ran some models at an even higher resolution, 512^{3}, to doublecheck 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_{1}, where k_{1} is the smallest wave number in the periodic domain, implying that the nondimensional size of the domain is (2π)^{3}.
As the initial conditions: pressure and density are set constant and with value of 1 throughout the box (making ρ_{0} = 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^{−6} in code units, uncorrelated at each point for the three components, corresponding to a E_{k} ∼ k^{4} 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
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:
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:
Here the first term on the righthand 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 ∇^{2}a = ∇(∇ ⋅ a)−∇ × (∇ × a), we obtain
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 ⟨ω^{2}⟩.
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 (EliasLó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 ϕ_{0} 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 ϕ_{0} 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 timeaveraged 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 k_{f} = 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 EliasLó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 nonisothermal runs are not shown because they reach locally supersonic flows and remain numerically stable for too few time steps. Despite attempting to select a ϕ_{0} that makes u_{rms} approximately independent of R, we succeeded in doing so only for the isothermal nonrotating case, while for the other nonisothermal 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 nonisothermal 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 ϕ_{0} (for τ_{cool} = 0.1).
Fig. 2. Different diagnostic quantities, ω_{rms}, u_{rms}, k_{ω}/k_{f}, and Re (from top to bottom), as a function of explosion width, R, for the forcing f_{acc}. These runs have Pm set to 1, i.e., Rm = Re, with a grid size of 256^{3} 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., k_{f} = 4, different resolutions, and varying the type of forcing. 
As expected the contribution of rotation to ω is much greater than the baroclinc one in nonisothermal 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 nonisothermality with an ideal gas law and different cooling times. Thus, the results of EliasLó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 EliasLó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 t_{turn}) 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 t_{shear} 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 E_{mag} 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 dashdotted lines correspond to growth rates in time units of t_{shear} and solid lines in forcing turnover times, t_{turn}. Both time units are tabulated in Appendix A. The points marked with crosses are the corresponding dynamo runs with f_{mom}. 
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 EliasLó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.
Vorticity and magnetic energy growth rates for different values of Pm.
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 smallscale contribution and show again a wide dynamical range decreasing in a k^{−2} 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 meanfield approaches based on the secondorder 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^{4} K,
If, as done by EliasLó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 onephase 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 EliasLó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 nonisothermal 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^{2}⟩ is relevant, while 2ν⟨S∇lnρ ⋅ q⟩ is more than one order of magnitude lower. This last statement holds true for all the isothermal, nonisothermal 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 ω_{rms}^{2}/2 (left axis and gray line) and vorticity growth terms (right axis and various colored lines for each term) for a nonisothermal 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 t_{turn}. However, the saturation regime does not present any relevant changes in time for more than 2000 t_{turn}, 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., ⟨F_{visc} ⋅ q⟩, overlap with ν⟨q^{2}⟩. 
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 nonisothermal 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 ω_{rms}^{2}/2 and vorticity growth terms for a nonisothermal 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 ⟨F_{visc} ⋅ q⟩ overlap with ν⟨q^{2}⟩. 
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^{2}⟩ 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 ω_{rms}^{2}/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 EliasLó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 curlfree 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 CEX2020001058M. 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 (CEX2020001058M). We are grateful to Axel Brandenburg, Matthias Rheinhardt, Eva Ntormousi, Federico Stasyszyn, Andrea Pallottini, Amit Seta and Claudia SorianoGuerrero 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 AECT202320034 (PI FDS).
References
 Achikanath Chirakkara, R., Federrath, C., Trivedi, P., & Banerjee, R. 2021, Phys. Rev. Lett., 126, 091103 [Google Scholar]
 Asvarov, A. I. 2014, A&A, 561, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brandenburg, A., & Ntormousi, E. 2023, ARA&A, 61, 561 [Google Scholar]
 Del Sordo, F., & Brandenburg, A. 2011, A&A, 528, A145 [Google Scholar]
 Dosopoulou, F., Del Sordo, F., Tsagas, C. G., & Brandenburg, A. 2012, Phys. Rev. D, 85, 063514 [Google Scholar]
 EliasLópez, A., Del Sordo, F., & Viganò, D. 2023, A&A, 677, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Federrath, C., RomanDuval, J., Klessen, R. S., Schmidt, W., & Mac Low, M. M. 2010, A&A, 512, A81 [Google Scholar]
 Ferrière, K. 2020, Plasma Phys. Controlled Fusion, 62, 014014 [Google Scholar]
 Franchetti, N. A., Gruendl, R. A., Chu, Y.H., et al. 2012, AJ, 143, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Guerrero, G., Del Sordo, F., Bonanno, A., & Smolarkiewicz, P. K. 2019, MNRAS, 490, 4281 [Google Scholar]
 Jones, C. A. 2011, Ann. Rev. Fluid Mech., 43, 583 [CrossRef] [Google Scholar]
 Kahniashvili, T., Brandenburg, A., Campanelli, L., Ratra, B., & Tevzadze, A. G. 2012, Phys. Rev. D, 86, 103005 [Google Scholar]
 Käpylä, P. J., Mitra, D., & Brandenburg, A. 2009, Phys. Rev. E, 79, 016302 [Google Scholar]
 Käpylä, P. J., Korpi, M. J., & Brandenburg, A. 2010, MNRAS, 402, 1458 [Google Scholar]
 Krause, F., & Raedler, K. H. 1980, Meanfield Magnetohydrodynamics and Dynamo Theory (Oxford: Pergamon Press) [Google Scholar]
 Mee, A. J., & Brandenburg, A. 2006, MNRAS, 370, 415 [NASA ADS] [Google Scholar]
 Monteiro, G., Guerrero, G., Del Sordo, F., Bonanno, A., & Smolarkiewicz, P. K. 2023, MNRAS, 521, 1415 [Google Scholar]
 Pencil Code Collaboration (Brandenburg, A., et al.) 2021, J. Open Source Software, 6, 2807 [Google Scholar]
 Rädler, K.H., & Rheinhardt, M. 2007, Geophys. Astrophys. Fluid Dyn., 101, 117 [Google Scholar]
 Rädler, K.H., Brandenburg, A., Del Sordo, F., & Rheinhardt, M. 2011, Phys. Rev. E, 84, 046321 [Google Scholar]
 Reichert, M., Obergaulinger, M., Aloy, M. Á., et al. 2023, MNRAS, 518, 1557 [Google Scholar]
 Seta, A., & Federrath, C. 2022, MNRAS, 514, 957 [Google Scholar]
 Skoutnev, V., Squire, J., & Bhattacharjee, A. 2022, MNRAS, 517, 526 [Google Scholar]
 Stefani, F., Stepanov, R., & Weier, T. 2021, Sol. Phys., 296, 88 [Google Scholar]
 Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253 [Google Scholar]
 Warnecke, J., KorpiLagg, M. J., Gent, F. A., & Rheinhardt, M. 2023, Nat. Astron., 7, 662 [Google Scholar]
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 EliasLó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_{0} = 10^{−6} as the seed field, Δt = 0.02 as interval between two different explosions, c_{s0}^{2} = 1 in the case of nonisothermal runs and an amplitude of the shearing profile A= 0.2 in the shearing cases.
Barotropic (i.e., isothermal) runs without shear.
Baroclinic (i.e., nonisothermal) runs without shear.
Models with sinusoidal shearing velocity profile. ND stands for the growth rates not determined, even though dynamo was present.
All Tables
Models with sinusoidal shearing velocity profile. ND stands for the growth rates not determined, even though dynamo was present.
All Figures
Fig. 1. Time evolution for ω_{rms} (left) and timeaveraged 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 k_{f} = 2/R. 

In the text 
Fig. 2. Different diagnostic quantities, ω_{rms}, u_{rms}, k_{ω}/k_{f}, and Re (from top to bottom), as a function of explosion width, R, for the forcing f_{acc}. These runs have Pm set to 1, i.e., Rm = Re, with a grid size of 256^{3} 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., k_{f} = 4, different resolutions, and varying the type of forcing. 

In the text 
Fig. 4. Time evolution (in units of t_{turn}) 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 t_{shear} 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 E_{mag} 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 dashdotted lines correspond to growth rates in time units of t_{shear} and solid lines in forcing turnover times, t_{turn}. Both time units are tabulated in Appendix A. The points marked with crosses are the corresponding dynamo runs with f_{mom}. 

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 ω_{rms}^{2}/2 (left axis and gray line) and vorticity growth terms (right axis and various colored lines for each term) for a nonisothermal 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 t_{turn}. However, the saturation regime does not present any relevant changes in time for more than 2000 t_{turn}, 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., ⟨F_{visc} ⋅ q⟩, overlap with ν⟨q^{2}⟩. 

In the text 
Fig. 8. Time evolution for ω_{rms}^{2}/2 and vorticity growth terms for a nonisothermal 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 ⟨F_{visc} ⋅ q⟩ overlap with ν⟨q^{2}⟩. 

In the text 
Fig. 9. Time evolution for ω_{rms}^{2}/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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.