Chaotic transients and hysteresis in an $\alpha^{2}$ dynamo model

The presence of chaotic transients in a nonlinear dynamo is investigated through numerical simulations of the 3D magnetohydrodynamic equations. By using the kinetic helicity of the flow as a control parameter, a hysteretic blowout bifurcation is conjectured to be responsible for the transition to dynamo, leading to a sudden increase in the magnetic energy of the attractor. This high-energy hydromagnetic attractor is suddenly destroyed in a boundary crisis when the helicity is decreased. Both the blowout bifurcation and the boundary crisis generate long chaotic transients that are due, respectively, to a chaotic saddle and a relative chaotic attractor.


Introduction
Chaotic transients are a common phenomenon in fluids and plasmas, being usually associated with decaying turbulence, where an initially erratic fluid converges to a laminar state. A typical example is a pipe flow, where turbulent puffs can last for a long time, but eventually disappear if the pipe is sufficiently long and the Reynolds number energy remains high (the strong-field branch) when started with a strong field. Thus, the dynamo displays a hysteresis near the dynamo onset.
The motivation of the present study is to explore the origin of chaotic transients nearby the onset of an α 2 dynamo in the presence of hysteresis. For this, we consider a simple 3D MHD model of isothermal compressible fluid which is driven by a helical forcing function. The advantage of including an external driver is that we can vary the net helicity in our model to explore the dynamo transition. In section 2, the α 2 dynamo model is described; section 3 shows the main results of the paper and the conclusions are given in section 4.

THE MODEL
We adopt the model of α 2 dynamo employed in Refs. [20,21]. The fluid is assumed to be isothermal and compressible, with constant sound speed c s , constant dynamical viscosity µ, constant magnetic diffusivity η and constant magnetic permeability µ 0 . The governing equations are: where ρ is the density, u is the fluid velocity, A is the magnetic vector potential, J = ∇ × B/µ 0 is the current density, p is the pressure, f is an external forcing function and ∇p/ρ = c 2 s ∇ lnρ where c 2 s = γp/ρ is assumed to be constant. The magnetic induction equation (3) is written for the vector potential A to ensure a solenoidal magnetic field, since ∇·B = ∇·(∇×A) = 0. The logarithmic density is also adopted in Eqs. (1) and (2) for numerical reasons, since it varies spatially much less than density. The domain is a box with dimensions L x = L y = L z = 2π and periodic boundary conditions in all three directions for all variables. We adopt non-dimensional units with k 1 = c s = ρ 0 = µ 0 = 1, where ρ = ρ is the spatial average of ρ and k 1 is the smallest wave number in the simulation box. Time is measured in units of (c s k 1 ) −1 , space is measured in units of k −1 1 , u in units of c s , B in units of √ µ 0 ρ 0 c s , ρ in units of ρ 0 and the magnetic diffusivity η is in units of c s /k 1 . Equations (1)-(3) are solved with the PENCIL CODE ‡, a thoroughly tested MHD solver frequently employed in Astrophysical works (see, e.g., [22] and the journal issue dedicated to the physics and algorithms of the PENCIL CODE). The code adopts sixthorder finite differences in space and third-order variable step Runge-Kutta in time. The initial conditions are ln ρ = 0 and u = 0. The initial magnetic vector potential is modeled by noise with Gaussian distribution with zero mean and standard deviation ‡ http://pencil-code.nordita.org/ equal to 10 −3 . Kinetic energy is injected in the system by a forcing function f, which is defined by [20,23] f( where k(t) = (k x , k y , k z ) is a time dependent wave vector, x = (x, y, z) is position, and φ(t), with |φ| < π, is a random phase. Here, N = f 0 c s (kc s /δt) 1/2 , where f 0 is a nondimensional factor, k = |k|, and δt is the integration time step. We choose the forcing wave number k around k f = 5 and every time step, a vector k(t) with 4.5 < k < 5.5 is randomly selected from a set of 350 previously generated vectors with the given wave number. The operator f k is given by whereê is an arbitrary unit vector needed in order to generate a vector k ×ê which is perpendicular to k. Note that for σ = 1, |f k | 2 = 1 and the kinetic helicity density of this forcing function satisfies at each point in space. For σ = 0, the forcing function is nonhelical, so σ is a measure of the kinetic helicity of the forcing. This rather complex forcing function has been employed in several previous works for being able to generate turbulent statistics even for moderate Reynolds numbers [24] and for us it is interesting because it allows us to choose the level of helicity added to the flow, an important element for the generation of large-scale dynamos. The Reynolds number and magnetic Reynolds number are defined, respectively, as where ν = µ/ρ is the kinematic viscosity, u rms ≡ u 2 1/2 is the root-mean-square (r.m.s) velocity and spatial average is denoted by · .

Results
First we must choose the grid resolution for the numerical simulation of Eqs. (1)-(3). Figure 1 shows a comparison of the time series of the root mean square of the magnetic (B rms , black line) and velocity (u rms , red line) fields in log-linear axes for simulations using 64 3 (a) and 128 3 (b) grid points. The parameter values are f 0 = 0.07, ν = η = 0.005 and σ = 1. Initially, the seed magnetic field is too weak to impact the velocity field, but the velocity field has a strong impact on the induction equation (3), causing an exponential growth of the magnetic energy during the initial (kinematic) phase of the dynamo. The growth rate γ in this phase can be found as the slope of a fitted line in the log-linear plot, and it is approximately 0.05 for both resolutions.
As the magnetic flux grows, it starts to strongly affect the velocity field through the Lorentz force in the momentum equation (2), causing a rapid decrease in the kinematic energy around t = 200 for both resolutions. Eventually, the mean magnetic and kinetic energies reach a saturated value which is approximately the same for both resolutions, for which we have R e = R m ≈ 20. Based on that, and considering the large amount of long time series that must be computed in this work, we adopted the lower resolution of 64 3 points in the following sections. Next, we set f 0 = 0.07, ν = η = 0.005 and vary σ as a control parameter. Using σ as the control parameter is a natural choice, since it is known that the presence of kinetic helicity in the flow can be favorable for magnetic field generation [25], although it is not strictly needed for either small-or large-scale dynamo to operate (see [26,27] and references therein).
Since our forcing function has k around k f = 5, kinetic energy is injected at this scale in the flow, inducing the production of a series of eddies with that wave number in the physical space. When σ = 1, helicity is maximum in the flow, causing an inverse energy transfer from k = 5 to k = 1 that has been related to the α-effect in [20]. This can be observed by plotting the one-dimensional power spectra, as in Fig. 2, where the kinetic (red dashed line) and magnetic (black solid line) spectra are computed as the integrated energy along spherical shells in the k = (k x , k y , k z ) space for t = 1000, well inside the nonlinearly saturated dynamo regime. It is clear that the kinetic spectrum is peaked at k = 5 and the magnetic at k = 1. This causes the appearance of large scale magnetic structures that are illustrated in the upper panel of Fig. 3. Note that B y and B z display a large scale oscillation in addition to the small scale fluctuations. The lower panel shows the velocity field components, whose largest structures are at the k f scale.

Bifurcation diagram
The onset of dynamo action is shown in Figure 4, which represents the bifurcation diagram as a function of σ for the time-averaged root mean square magnetic fieldB rms . For each value of σ, a weak seed magnetic field is used as initial condition, the initial transient is discarded and then, time averages of B rms are plotted. The magnetic energy of the attractor is zero when σ is below 0.21, implying that there is no dynamo action and the attracting state is purely hydrodynamic. For σ > 0.21 dynamo action takes place and there is a sudden jump in the B rms of the attractor. The average magnetic energy of this hydromagnetic attractor keeps growing from then on, up toB rms ≈ 0.275 for σ = 1. Note that at the critical transition, the hydrodynamic attractor loses stability and small amplitude magnetic perturbations are sufficient to drive the system toward the hydromagnetic attractor. In previous works [28], transition to dynamo in MHD simulations in a periodic box with the helical ABC-forcing was shown to be due to a nonhysteretic blowout bifurcation. In this type of bifurcation, the dynamical system has a smooth invariant manifold within which lies a chaotic attractor for parameter values less than a critical value. As the parameter value is increased, a blowout bifurcation takes place, in which the manifold loses its attracting property and the chaotic set on it ceases to be an attractor. Right after the transition, solutions display on-off intermittency, i.e., they spend long times very near the manifold, then are "blown out" of it in fast bursts where they move far from the manifold. After each burst, trajectories go back to the vicinity of the manifold and the process repeats intermittently [29]. We tried in vain to find intermittent bursts in the transition to dynamo using the forcing function in the form (4). Figures 5(a)-(b) show that just before the transition, small magnetic perturbations decay and the solutions converge to the purely hydrodynamic (B = 0) manifold. Thus, this manifold is attracting and there is a chaotic hydrodynamic attractor in it (velocity field fluctuations are always chaotic in our work). Very near the transition, some magnetic bursts may occur as in Fig. 5(b) for σ = 0.21374, but they have a tiny amplitude and, eventually, the hydrodynamic manifold attracts the solution and no more bursts are observed. For σ = 0.21379, right after the transition, the solution stays near the manifold for a long time before suddenly jumping toward a chaotic attractor with high magnetic energy, the hydromagnetic attractor. This shows that the hydrodynamic manifold has lost transversal stability, since initial conditions with B = 0 stay on the manifold for all parameter values, but even small nonzero values of B (i.e., perturbations that are transversal to the hydrodynamic manifold) are able to expel trajectories away from the manifold and toward the hydromagnetic attractor. This means that the previous hydrodynamical attractor has also lost its transversal stability and what is left is a transient chaotic hydrodynamic phase. It could be said that the hydrodynamic chaotic attractor has become a chaotic saddle, but since this chaotic set still attracts all initial conditions on the hydrodynamic manifold, which itself has become unstable, we refer to this hydrodynamic chaotic set as a relative chaotic attractor, adopting a nomenclature introduced by Skufca et al. [30]. Note that the stable manifold of the relative attractor in question is not a fractal structure, but the whole hydrodynamic subspace defined by B = 0. The reason for the absence of high-amplitude intermittent bursts in our dynamo transition is explained in the next section.

Hysteresis and chaotic transients
Motivated by Karak et al. [19], we search for hysteresis in this dynamo system. Although their work has an imposed uniform large-scale shear flow that is absent in our simulations, we still managed to find a hysteresis in our α 2 -dynamo model. Recall from section 3.1 that the transition to the hydromagnetic attractor takes place at σ = σ c ≈ 0.21379, where random seed magnetic fields are amplified; for σ < σ c , the seed fields decay to zero, as the solutions approach a hydrodynamic attractor. However, if one takes as initial condition a state with a high energy magnetic field, it may not decay to the hydrodynamic attractor, as shown in Fig. 6. Here, the initial condition is a state taken from the hydromagnetic attractor at σ = 0.3; when the control parameter is reduced to σ = 0.2 ( Fig. 6(a)) and σ = 0.199 ( Fig. 6(b)), the solutions remain in the hydromagnetic attractor, with the magnetic energy still following the upper branch in Fig. 4. That is a signature of a hysteresis, as the saturation of the magnetic field amplitude depends on the previous history of the control parameter. For lower values of σ, the hydromagnetic attractor loses stability and becomes a hydromagnetic chaotic saddle, leaving a chaotic transient in the region of phase space previously occupied by the attractor. Two of such long chaotic transients are exhibited in Fig. 7 in the form of time series of B rms , and in Fig. 8 in the form of the spatiotemporal evolution of the B y component averaged in the horizontal plane (B x , B y ) as a function of z and time.
The hydromagnetic chaotic attractor is shown in Fig. 9(a) for σ = 0.199, where the chaotic trajectories represent the temporal variation of the components of the magnetic field vector B = (B x (x 0 , y 0 , z 0 , t), B y (x 0 , y 0 , z 0 , t), B z (x 0 , y 0 , z 0 , t)) computed at the origin of the spatial domain (x 0 , y 0 , z 0 ) = (0, 0, 0). A hydromagnetic chaotic saddle is shown in Fig. 9(b) for σ = 0.198. It hints that both chaotic sets occupy approximately the same region in the phase space, but a small variation in the parameter σ is sufficient to considerably reduce the size of the chaotic set, following the rapid decay of the magnetic energy in the upper branch of the bifurcation diagram in Fig. 4.  As illustrated by Fig.  7, the average duration of the chaotic transients decreases with decreasing σ. This observation suggests that the destabilization of the Figure 9.
Trajectory of initial conditions on invariant sets in the (B x (x 0 , y 0 , z 0 , t), B y (x 0 , y 0 , z 0 , t), B z (x 0 , y 0 , z 0 )) space for the hydromagnetic attractor at σ = 0.199 (a) and for the hydromagnetic chaotic saddle at σ = 0.198 (b). The point (x 0 , y 0 , z 0 ) is at the origin of the spatial domain.
hydromagnetic chaotic attractor at σ = σ bc ≈ 0.198 is due to a boundary crisis, where a chaotic attractor collides with the boundary of its own basin of attraction, leading to the destruction of both the attractor and its boundary. A chaotic saddle is then left in place of the chaotic attractor and, as shown by Grebogi et al. [31], the average duration τ of the chaotic transients near a boundary crisis decays with the distance from the crisis parameter value σ bc following the law for a negative γ. We computed τ for a set of values of σ close to σ bc and obtained the results shown in Fig. 10, where the fitted line has slope γ = −0.04. The following procedure was adopted to produce this figure. First, a set of 100 initial conditions are selected from the hydromagnetic chaotic attractor at σ = 0.3 > σ bc ; then, those initial conditions are used to generate transient chaotic time series for a given σ < σ bc ; the transient time for each initial condition is recorded when B rms is below a certain threshold (B rms < 0.01) and the average transient time τ is computed from the 100 time series. This process is repeated for the 10 values of σ < σ bc shown in the figure.
The fitted line was obtained by using linear regression.
Having established the existence of hysteresis in the current α 2 dynamo model, we go back to the previously raised question of why strong intermittent bursts are not observed in the transition to dynamo near σ = 0.21. As mentioned before, the transition in the ABC dynamo is due to a nonhysteretic blowout bifurcation, where intermittency is present [28]. Differently, the hysteresis in our model mean that there is a bistability region in the parameter space for 0.199 ≤ σ ≤ 0.213 where two chaotic attractors coexist, the hydrodynamic and the hydromagnetic attractors. Our results reveal that a hysteretic blowout bifurcation is responsible for the dynamo transition. According to Ott and Sommerer [29], in this type of bifurcation there is an attracting chaotic set in the (hydrodynamic) invariant manifold and this chaotic set loses stability as the system parameter σ increases through a critical value σ c (as in Fig. 4). For σ < σ c orbits not in the basin of the attractor on the invariant manifold go to some other (hydromagnetic) attractor off the manifold (as in Fig. 6). As σ increases through σ c the attractor loses stability. For σ slightly greater than σ c , points started near the manifold can experience a chaotic transient in which their orbits initially closely resemble those on the σ < σ c manifold attractor, but almost all points started near the manifold eventually move off to the other attractor off the manifold (as in Fig. 5(c)). Thus, on-off intermittency is not present in a hysteretic blowout bifurcation.
Our findings are summarized in Fig. 11. The horizontal line atB rms = 0 represents the hydrodynamic manifold, which is attracting for σ < σ c (solid red line) and contains a chaotic attractor. For σ > σ c , the hydrodynamic chaotic attractor loses stability in a hysteretic blowout bifurcation and becomes a hydrodynamic "relative attractor" (dashed red line), where hydrodynamic chaotic transients are observed for initial perturbations near the hydrodynamic manifold. In the upper branch, the solid black line represents the hydromagnetic chaotic attractor, which loses stability in a boundary crisis at σ bc , giving rise to a hydromagnetic chaotic saddle (dashed black line), responsible for chaotic transients involving the magnetic field. There is a bistability window between σ bc and σ c where the hydromagnetic chaotic attractor coexists with the hydrodynamic chaotic attractor.

Conclusion
We have found a transition to dynamo in a hysteretic blowout bifurcation in direct numerical simulations of an MHD α 2 dynamo model. Our transition differs from [28] in that theirs is due to a non-hysteretic blowout bifurcation and it differs from [19] in that they adopted an α − ω dynamo model. To our knowledge, this is the first time that a hysteretic blowout bifurcation is observed in an α 2 dynamo model. Previously, a nonhysteretic blowout bifurcation had been observed in a truncated mean field dynamo model [32]. The hydromagnetic chaotic attractor in the upper hysteresis branch is destroyed in a boundary crisis. Although we have not formally characterized the blowout bifurcation and the boundary crisis, the behavior of the chaotic transients generated by our transitions strongly suggests that this is, indeed, the case. Thus, our work illustrates the importance of chaotic transients in revealing the complex dynamical transitions present in turbulent flows.