Numerical evidence for a small-scale dynamo approaching solar magnetic Prandtl numbers

Magnetic fields on small scales are ubiquitous in the Universe. Although they can often be observed in detail, their generation mechanisms are not fully understood. One possibility is the so-called small-scale dynamo (SSD). Prevailing numerical evidence, however, appears to indicate that an SSD is unlikely to exist at very low magnetic Prandtl numbers (PrM) such as those that are present in the Sun and other cool stars. Here we have performed high-resolution simulations of isothermal forced turbulence using the lowest PrM values achieved so far. Contrary to earlier findings, the SSD not only turns out to be possible for PrM down to 0.0031 but also becomes increasingly easier to excite for PrM below about 0.05. We relate this behaviour to the known hydrodynamic phenomenon referred to as the bottleneck effect. Extrapolating our results to solar values of PrM indicates that an SSD would be possible under such conditions. High-resolution simulations of the small-scale dynamo (SSD) mechanism, with close-to-realistic parameters of deep stellar convection zones, indicate that SSDs are possible in the Sun and other cool stars, in contrast to previous theoretical expectations.


Results
We include simulations with resolutions of 256 3 to 4608 3 grid points and Re = 46 to 33000. This allows us to explore the parameter space from Pr M = 1 to 0.0025, which is closer to the solar value than has been investigated in previous studies. For each run, we measure the growth rate λ of the magnetic field in its kinematic stage and determine whether or not an SSD is being excited.
To afford an in-depth exploration of the effect of Pr M , we omit large-scale effects such as stratification, rotation and shear. We avoid the excessive integration times, required to simulate convection, by driving the turbulent flow explicitly under isothermal conditions. Our simulation setup consists of a fully periodic box with a random volume force, see Online Methods for details; the flow exhibits a Mach number of around 0.08. In Fig. 1, we visualize the velocity and magnetic fields of one of the highest resolution and Reynolds number cases. As might be anticipated for low Pr M turbulence, the flow exhibits much finer, fractal-like structures than the magnetic field. Note, that all our results refer to the kinematic stage of the SSD, where the magnetic field strength is far too weak to influence the flow and otherwise arbitrary.

Growth rates and critical magnetic Reynolds numbers
In Fig. 2 we visualize the growth rate λ as function of Re and Re M . We find positive growth rates for all sets of runs with constant Pr M if Re M is large enough. λ increases always with increasing Re M as expected. Surprisingly, the growth rates are significantly lower within the interval from Re = 2000 to The growth rates for Pr M = 0.1 match very well the ones from [10], indicated by triangles. From Fig. 2, we clearly see that the critical magnetic Reynolds number Re crit M , defined by growth rate λ = 0, first rises as a function of Re and then falls for Re > 3 × 10 3 , see the thin black line. Looking at Re crit M as a function of magnetic Prandtl number Pr M , it first increases with decreasing Pr M and then decreases for Pr M < 0.05. Hence, an SSD is easier to excite here than for 0.05 < Pr M < 0.1. We could even find a nearly marginal, positive growth rate for Pr M = 0.003125. The decrease of λ at low Pr M is an important result as the SSD was believed to be even harder [4,9] or at least equally hard [7,8] to excite when Pr M was decreased further from previously investigated values. The growth rates agree qualitatively with the earlier work at low Pr M [6][7][8].
For a more accurate determination of Re crit M , we next plot the growth rates for fixed Pr M as a function of Re M , see Fig. 3(a). The data are consistent with λ ∝ ln (Re M /Re crit M ) as theoretically predicted [36,37]. Fitting accordingly, we are able to determine Re crit M as a function of Pr M , see Fig. 3 . Extrapolating this to the Sun and solar-like stars would lead to Re crit M ≈ 40 at Pr M = 10 −6 , which means that we could expect an SSD to be present. For increasing Re, by decreasing ν, it would be reasonable to assert that the statistical properties of the flow and hence Re crit M become independent of Pr M . However, episodes of non-monotonic behavior of Re crit M when approaching this limit cannot be ruled out. The well-determined Re crit M dependency on Pr M together with its error bars and the power-law fit have been added to Fig. 2, and agree very well with the thin black line for λ = 0 interpolated from the growth rates.  In all cases the kinetic energy clearly follows a Kolmogorov cascade with E kin ∝ k −5/3 in the inertial range. When compensating with k 5/3 , we find the well-known bottleneck effect [38,39]: a local increase in spectral energy, deviating from the power-law, as found both in fluid experiments [40][41][42] and numerical studies [43,44]. It has been postulated to be detrimental to SSD growth [4,10]. For the magnetic spectrum on the other hand, yet clearly visible only for Pr M ≤ 0.005, we find a power-law following E mag ∝ k −3 . A 3/2 slope at low wavenumbers as predicted by [45] is seen only in the runs with Pr M close to one, while for the intermediate and low Pr M runs, the positive-slope part of the spectrum shrinks to cover only the lowest k values, and the steep negative slopes at high k values become prominent. A steep negative slope in the magnetic power spectra was also seen by [7] for Pr M slightly below unity. However, the authors propose a tentative power of −1 given that the −3 slope is not yet clearly visible for their Pr M values.

Regions of dynamo excitation
Analyzing our simulations, we adopt the following strategy: For each spectrum, we determine the wavenumber of the bottleneck, k b , as the location of its maximum in the (smoothed) compensated spectrum, along with its starting point k bs < k b at the location with 75% of the maximum, see the middle-row panels of Fig. 4. We additionally calculate a characteristic magnetic wavenumber, defined as k M = k E mag (k)k dk/ k E mag (k) dk, which is often connected with the energy-carrying scale. Furthermore, we calculate the viscous dissipation wavenumber k ν = (ϵ K /ν 3 ) 1/4 following Kolmogorov theory, where ϵ K is the viscous dissipation rate 2νS 2 with the traceless rate-of-strain tensor of the flow, S. From the relations between these four wavenumbers (listed in Supplementary Table 2), we will draw insights about the observed behavior of Re crit M with respect to Pr M .
We plot k b /k ν and k bs /k ν as functions of Pr M in Fig. 5. As is expected, k b /k ν , or the ratio of the viscous scale to the scale of the bottleneck, does not depend on Pr M , as the bottleneck is a purely hydrodynamic phenomenon. The start of the bottleneck k bs should likewise not depend on Pr M , but the low Re values for Pr M = 1 to 0.1 lead to apparent thinner bottlenecks, hence an unsystematic weak dependency. The red shaded area between k b and k bs is the low-wavenumber part of bottleneck where the slope of the spectrum is larger (less negative) than −5/3 see Supplementary Table 2 for values of the modified slope α b and Supplementary, Section 1 for a discussion. We note that α b ≈ −1.3 . . . − 1.5 and can thus deviate significantly from −5/3. Overplotting the k M /k ν curve reveals that it intersects with the red shaded area exactly where the dynamo is hardest to excite (region II). This lets us conclude that the shallower slope of the low-wavenumber part of the bottleneck may indeed be responsible for enhancing Re crit M in the interval 0.04 ≤ Pr M ≤ 0.1. Using this plot, we can now clearly explain the three regions of dynamo excitation. For 0.1 ≤ Pr M ≤ 1 the low-wavenumber part of bottleneck and the characteristic magnetic scale are completely decoupled. This makes the SSD easy to excite (region I). For 0.04 ≤ Pr M ≤ 0.1, (grey, region II), the dynamo is hardest to excite because of the shallower slope of the kinetic spectra. In region III, where Pr M ≤ 0.04 the low-wavenumber part of bottleneck and the characteristic magnetic scale are again completely decoupled making the dynamo easier to excite.
Further, we find that the dependence of k M /k ν on Pr M also differs between the regions. In region I k M /k ν depends on Pr M via k M /k ν ∝ Pr 0.54 M and in region II and III via k M /k ν ∝ Pr 0.71 M . This becomes particularly interesting when comparing the characteristic magnetic wavenumber k M with the ohmic dissipation wavenumber which is defined as k η = k ν Pr 3/4 M . In region I, we find a significant difference of k M and k η in value and scaling. However, in region III the scaling of k M comes very close to the 3/4 scaling of k η . This behaviour can be even better seen in the inset of Fig. 5, where the ratio k M /k η is 0.3 for Pr M = 1 and tends towards unity for decreasing Pr M , but is likely to saturate below 0.75.

Discussion
In conclusion, we find that the SSD is progressively easier to excite for magnetic Prandtl numbers below 0.04, in contrast to earlier findings, and thus is very likely to exist in the Sun and other cool stars. Provided saturation at sufficiently high levels, the SSD has been proposed to strongly influence the dynamics of solar-like stars: previous numerical studies, albeit at Pr M ≈ 1, indicate that this influence concerns for example the angular momentum transport [19,20], and the LSD [21][22][23][24][25]. Our kinematic study, however, only shows that a positive growth rate is possible at very low Pr M , but not whether an SSD is able to generate dynamically important field strengths. As the Re M of the Sun and solar-like stars is several orders of magnitude higher than the extrapolated Re crit M value of 40, we yet expect dynamically important SSDs as indicated by Pr M = 1 simulations [15]. However, numerical simulations with Pr M down to 0.01 show a decrease of the saturation strength with decreasing Pr M [46].
The results of our study are well in agreement with previous numerical studies considering partly overlapping Pr M ranges [6][7][8]10]. Those found some discrepancies with the Kazantsev theory [45] for low Pr M , for example the narrowing down of the positive Kazantsev spectrum at low and intermediate wavenumbers, and the emergence of a negative slope instead at large wave numbers [7]. We could extend this regime to even lower Pr M and therefore study these discrepancies further. For Pr M ≤ 0.005 we find that the magnetic spectrum shows a power-law scaling k −3 , which is significantly steeper than the tentative k −1 one proposed in [7] for 0.03 ≲ Pr M ≲ 0.07 (but only for 8thorder hyperdiffusivity). This latest finding of such a steep power law in the magnetic spectrum challenges the current theoretical predictions and might indicate that the SSD operating at low Pr M is fundamentally different from that at Pr M ≈ 1.
Secondly, we find that the growth rates near the onset follow an ln(Re M ) dependence as predicted by [36,37], and not a Re 1/2 M one as would result from intertial-range-driven SSD [1,7]. We do not observe a tendency of the growth rate to become independent of Re M at the highest Pr M either, which could be an indication of an outer-scale driven SSD, as postulated by [7]. Furthermore, we find that the pre-factor of γ ∝ ln(Re M /Re crit M ) is nearly constant with its mean around 0.022, in agreement with 0.023 of [10]. A constant value means that the logarithmic scaling is independent of Pr M and seems to be of general validity.
Thirdly, we find that the measured characteristic magnetic wavenumber k M is always smaller than the estimated k η , and furthermore, k M is not always following the theory-predicted scaling of k η ∝ Pr 3/4 M with Pr M . For the region I, where Pr M is close to 1, this discrepancy is up to a factor of three and the deviation from the expected Pr M -scaling is most significant here. These discrepancies have been associated with the numerical setups injecting energy at a forcing scale far larger than the dissipation scale, i.e. k f ≪ k η [1]. Furthermore, our runs in region I also have relatively low Re and therefore numerical effects are not dismissible. In region III (low Pr M ), k M /k η is approaching the constant offset factor 0.75. Hence, the scaling of k M /k ν with Pr M gets close to the expected one. This result again indicates that the SSD at low Pr M is different from that at Pr M ≈ 1.
An increase of Re crit M with decreasing Pr M followed by an asymptotic levelling-off for Pr M ≪ 1 was expected in the light of theory and previous numerical studies. Instead, we found non-monotonic behavior as function of Pr M ; we could relate it to the hydrodynamical phenomenon of the bottleneck. If the characteristic magnetic wavenumber lies in the positive-gradient part of the compensated spectrum, where the spectral slope is significantly reduced from −5/3 to ∼ −1.4, the dynamo is hardest to excite (0.1 ≥ Pr M ≥ 0.04). For higher or lower Pr M , the dynamo becomes increasingly easier to excite. The local change in slope due to the bottleneck has often been related to an increase of the "roughness" of the flow [1,10,43], which is expected to harden dynamo excitation based on theoretical predictions [4,9] from kinematic Kazantsev theory [45]. In line with theory, the roughness-increasing part of the bottleneck appears decisive in our results, however, only when k M is used as a criterion. The usage of k η would in contrast suggest that the peak of the bottleneck is decisive [10]. Such interpretation appears incorrect, as the rough estimate of k η employed here does not represent the magnetic spectrum adequately and the peak of the bottleneck does not coincide with the maximum of "roughness".

Numerical setup
For our simulations, we use a cubic Cartesian box with edge length L and solve the isothermal magnetohydrodynamic equations without gravity, similar as in [5,47].
where u is the flow speed, c s is the sound speed, ρ is the mass density, and B = ∇ × A is the magnetic field with A being the vector potential. J = ∇ × B/µ 0 is the current density with magnetic vacuum permeability µ 0 , while ν and η are constant kinematic viscosity and magnetic diffusivity, respectively. The rate-of-strain tensor S ij = (u i,j + u j,i )/2 − δ ij ∇ · u/3 is traceless. The forcing function f provides random white-in-time non-helical transversal plane waves, which are added in each time step to the momentum equation, see [5] for details. The wavenumbers of the forcing lie in a narrow band around k f = 2k 1 with k 1 = 2π/L. Its amplitude is chosen such that the Mach number Ma = u rms /c s is always around 0.082, where u rms = ⟨u 2 ⟩ V is the volume and time-averaged root-mean-square value. The Ma values of all runs are listed in Supplementary Material Table 1. To normalize the growth rate λ, we use an estimated turnover time τ = 1/u rms k f ≈ 6/k 1 c s . The boundary conditions are periodic for all quantities and we initialise the magnetic field with weak Gaussian noise. Diffusion is controlled by the prescribed parameters ν and η. Accordingly, we define the fluid and magnetic Reynolds numbers with the forcing wavenumber k f as We performed numerical free decay experiments (see Supplementary Material, Section 7), from which we confirm that the numerical diffusivities are negligible.
The spectral kinetic and magnetic energy densities are defined via where B rms = ⟨B 2 ⟩ V is the volume-averaged root-mean-square value and ⟨ρ⟩ V the volume-averaged density. Our numerical setup employs a drastically simplified model of turbulence compared to the actual one in the Sun. There, turbulence is driven by stratified rotating convection being of course neither isothermal nor isotropic. However, these simplifications were to date necessary when performing a parameter study at such high resolutions as we do. Nevertheless, we can connect our study to solar parameters in terms of Pr M and Ma. Their chosen values best represent the weakly stratified layers within the bulk of the solar convection zone, where Pr M ≪ 1 and Ma ≪ 1. The anisotropy in the flow on small scales is much weaker there than near the surface and therefore close to our simplified setup.

Data availability
Data for reproducing Figs. 2, 3, and 5 are included in the article and its supplementary information files. The raw data (timeseries, spectra, slices, and snapshots) are provided through IDA/Fairdata service hosted at CSC, Finland, under https://doi.org/10.23729/ 206af669-07fd-4a30-9968-b4ded5003014. From the raw data, Figs. 1 and 4 can be reproduced.

Code availability
We use the Pencil Code [48]

Author Contributions Statement
JW was leading, but all the authors contributed, to the design and performing the numerical simulations. JW was leading the data analysis. MJKL was in charge of acquiring the computational resources from CSC. All the authors contributed to the interpretation of the results and writing up the manuscript.

Discussion on the roughness of the flow
We calculate the slope in the low-wavenumber part of the bottleneck by fitting E mag ∼ k α b in the interval k bs . . . k b . As shown in Supplementary Fig. 1 Seeking further insight into the dynamo operating in the three regimes of Pr M , we look at the spectral magnetic energy transfer functions. We follow the approach of [2], but using the convention by [3], where the contribution of compressibility is subsumed in the stretching/shearing and advection terms, T Str and T Adv , respectively, as follows (1) Here, the hats indicate the Fourier transform and c.c. refers to the complex conjugate expressions. In Supplementary Fig. 4a, we show these functions after shell integration, i.e. as functions of k = |k|, for six runs with 0.005 ≤ Pr M ≤ 0.8. As expected, the curves peak at higher wavenumbers for higher Pr M . We also look at the ratio of the wavenumbers at which T Str and T Adv has its maximum and minimum, respectively, k Str /k Adv , which is determined from the curves after smoothing. We find that this ratio is maximal where the dynamo is hardest to excite, shown as grey shade in the Supplementary Fig. 4b. However, we should note that the simulations with the two lowest Pr M have a higher resolution than the other four, which could also influence the result. It remains unclear whether, or how, this enhanced ratio would relate to the difficulty in exciting the dynamo at 0.

Numerical diffusivity
In our low Pr M simulations, η ν, hence numerical diffusion might only play a role for the velocity field. We estimate the numerical viscosity in our simulations in two ways. Firstly, we follow the approach of [2], to estimate the fluid Reynolds number from the power spectra. For this we determine the Taylor microscale and the integral scale for turbulent motions These scales can be used to estimate the effective Reynolds number as This is often used to estimate Re eff from observations or in simulations that use only numerical diffusivities, although the proportionality constant is unknown. As shown in Supplementary Fig. 5, the effective Re based on Equation (5) scales linearly with Re based on the prescribed value of ν. Assuming now plausibly numerical diffusion to enter mainly in the form of hyperdiffusion, we infer that from this linear relationship a notable influence of numerical diffusion can be ruled out.
Supplementary Figure 5: Effective Reynolds number proxy (5) as a function of Re from the prescribed ν. The red line indicates (λ int /λ TM ) 2 = 0.026 Re.
In addition, we performed a flow decay experiment to estimate the numerical diffusion from the flow decay rate in the absence of all other terms. The resulting effective viscosity is to high precision the same as the prescribed ν up to wavenumbers k 1000 k 1 , covering safely all scales relevant for the magnetic field. From these two results we conclude that our numerical simulations are not affected by numerical diffusion, and that the Pr M regimes are accurately identified.