The Tayler instability at low magnetic Prandtl numbers: between chiral symmetry breaking and helicity oscillations

The Tayler instability is a kink-type, current driven instability that plays an important role in plasma physics but might also be relevant in liquid metal applications with high electrical currents. In the framework of the Tayler-Spruit dynamo model of stellar magnetic field generation, the question of spontaneous helical (chiral) symmetry breaking during the saturation of the Tayler instability has received considerable interest. Focusing on fluids with low magnetic Prandtl numbers, for which the quasistatic approximation can be applied, we utilize an integro-differential equation approach in order to investigate the saturation mechanism of the Tayler instability. Both the exponential growth phase and the saturated phase are analyzed in terms of the action of the alpha and beta effects of mean-field magnetohydrodynamics. In the exponential growth phase we always find a spontaneous chiral symmetry breaking which, however, disappears in the saturated phase. For higher degrees of supercriticality, we observe helicity oscillations in the saturated regime. For Lundquist numbers in the order of one we also obtain chiral symmetry breaking of the saturated magnetic field.


Introduction
Electrical currents through an incompressible, viscous and resistive liquid conductor produce azimuthal magnetic fields which, beyond a critical field strength, become unstable to a non-axisymmetric, i.e. kink-type instability that we will call Tayler instability (TI) as a tribute to the seminal contributions of R.J. Tayler [6,7].
For a constant current density in an infinitely long cylinder Rüdiger et al had shown [8,9] that the governing parameter is the Hartmann number, Ha = B ϕ (R)R(σ/ρν) 1/2 , which has to exceed a value in the order of 20 for the TI to set in (B ϕ (R) is the azimuthal field at the outer radius R of the cylinder, σ, ρ and ν are the conductivity, density and viscosity of the fluid, respectively). This critical value of Ha is actually consistent with previous results [10,11,12,13] concerning the effects of viscosity and resistivity on the stability of various plasma z-pinches, if one leaves aside the effects of complicated boundary conditions and non-homogeneous material parameters in the plasma case. Note that at these early days [12] it was far from obvious that the governing parameter for the onset of the instability is Ha, rather than the Lundquist number S = HaP m 1/2 (with P m = νµ 0 σ denoting the magnetic Prandtl number, where µ 0 is the magnetic permeability constant).
Whilst the focus of fusion related pinch experiments was prominently on the plasma destabilization when the ratio of axial to azimuthal magnetic field (the so-called safety parameter) falls below a certain critical value [14], a recent liquid metal experiment, with uniform conductivity and viscosity as well as well-defined insulating boundary condition, has indeed confirmed the TI-threshold of Ha ≃ 20 [15]. From the application point of view, current-driven instabilities in liquid metals are presently considered a possible limitation for the integrity of large-scale liquid metal batteries. Such batteries are self-assembling stratified systems made of a heavy liquid metal or metalloid (e.g., Bi, Sb) at the bottom, a suitable molten salt mixture as electrolyte in the middle, and a light alkaline or earth alkaline metal (e.g., Na, Mg) at the top. While small versions of liquid metal batteries have already been tested [16,17,18,19], the occurrence of the TI could possibly present a serious problem for the stratification in larger batteries with prospected charging/discharging currents of some thousand amps. In [20] we had advised a simple trick to suppress the TI in liquid metal batteries by just returning the battery current through a bore in the centre. By the resulting change of the radial profile B ϕ (r) it is possible to avoid the (ideal) condition ∂(rB 2 ϕ (r))/∂r > 0 [7] for the onset of the TI.
In a follow-up paper [5], a numerical code has been presented that is capable of treating TI-problems at small values of P m as they are typical for liquid metals. This was achieved by replacing the solution of the induction equation for the magnetic field by applying the so-called quasistatic approximation [21]. This approximation allows to avoid the explicit time stepping of the magnetic field by computing the electrostatic potential by a Poisson equation, and deriving from this the electric current density. The induced magnetic field is then computed from the induced current density via Biot-Savart's law. This way one arrives at an integro-differential equation approach, as it had already been used by Meir and Schmidt for different magnetohydrodynamic (MHD) problems [22]. Our numerical scheme utilizes the open source CFD library OpenFOAM R [23], supplemented by an MPI-parallelized implementation of Biot-Savart's law. This code was then applied to a number of TI related problems, in particular for determining the scaling properties of the growth rate and the saturated velocity field, the dependence of the critical current on the geometric aspect ratio, as well as for validating various methods of preventing TI in liquid metal batteries [24]. Recently, our results were confirmed by another code working completely in the framework of the differential equation approach, by analyzing the scaling properties of the solutions with P m [25]. The authors also discussed carefully the limitations of the quasistatic approach for higher values of P m.
An interesting by-product of the battery-oriented simulations [5] was the observation of the transient occurrence, but ultimate disappearance, of helical structures during the evolution of the TI. On the first glance, the appearance of helical structures is surprising, since the underlying equations have no preference for left or right handed solutions. Yet, it is exactly this helical (or chiral) symmetry breaking that has gained considerable interest in various astrophysical problems. This applies in particular to the concept of the Tayler-Spruit dynamo [1] in which an azimuthal magnetic field is thought to become strong enough to drive the TI against the stable stratification in the radiation zone of a star. Combined with the usual differential rotation this effect might lead to a working dynamo. Despite of the attractiveness of the TI, in particular for explaining angular momentum transport in various types of stars [9,26,27], the concept of the actual Tayler-Spruit dynamo is not without caveats. Zahn et al [2] have argued that the TI-produced non-axisymmetric (m = 1) poloidal magnetic field alone would not be suited to close the dynamo loop (since the toroidal field wound up from it would have the same m = 1 dependence), but that some sufficiently large mean-field α effect would be needed to produce the necessary axisymmetric poloidal field.
It is exactly here where the question whether TI saturates with a finite helicity, produced by a finite α effect, becomes highly relevant. Some recent papers have answered this question affirmatively: Gellert et al [3] have found spontaneous chiral symmetry breaking of the TI in simulations with P m of 0.1, 1, and 10. Bonanno et al [4] got a similar result for very large P m = 10 7 . In addition to the numerical simulation, the latter authors developed a simple model of energy and helicity evolution resulting in an instructive phase portrait. The equations describing this behaviour can also be linked to a similar chiral symmetry breaking in biochemistry where it refers to the selection of one of two possible forms of bio-molecules (mainly sugars and amino acids) that are mirror images of each other [28].
With this background, the main motivation for the present paper is the discrepancy between the simulations of [3,4,29] and the preliminary result of our low P m simulations [5] showing that helicity starts to grow but ultimately decays to zero. Given the different P m at which the respective simulations were done, it is worthwhile to understand in detail the saturation mechanism of TI in dependence on P m.
Actually, helical states have a long history in plasma physics, tracing back to the early work of Lundquist on "Magneto-hydrostatic fields" [30]. Specializing general pressure-balanced fields to force-free fields that satisfy (∇ × B) × B = 0, he found, first, that fields with ∇ × B = a(r)B fulfill this demand, and second, that a(r) = const must be requested for the field to remain force-free during its time-evolution. For cylindrical geometry, Lundquist found that the force-free condition, i.e. the demand that the current is parallel to the field, is guaranteed by Bessel function profiles B z = AJ 0 (ar), B ϕ = AJ 1 (ar) (interestingly, the very same profiles for the velocity field turned later out to provide the most efficient dynamo of the Ponomarenko or Riga type [31]).
Soon after Lundquist's work, Chandrasekhar and Woltjer [32] interpreted this Bessel functions solution in terms of achieving "maximum magnetic energy for a given mean-squared current density" or, alternatively, as a "state of minimum dissipation for a given magnetic energy". Since Bessel functions also maximize the magnetic helicity for given magnetic energy (and magnetic helicity is a better conserved quantity than energy) a surge of work was devoted to understand how solutions of this kind can be achieved dynamically. This goes mainly under the notion of Taylor relaxation [33], and has found great interest in connection with the reversed field pinch. Quite a number of workers have tried to understand Taylor relaxation from different thermodynamic principles, such as minimum entropy production [34] or minimum dissipation rates [35,36,37,38,39,40,41].
One of the first applications of a general thermodynamic principle to plasma relaxation goes back to a note of Max Steenbeck [42] (relying, in turn, on an earlier idea of Compton and Morse [43]). Steenbeck's principle states that in real gas discharges at fixed current the heat power, and thus the voltage drop between the electrodes, is minimized (somewhat surprisingly, this minimum-dissipation principle corresponds perfectly with the maximum entropy production rate principle [44] if one considers the total system including the current-stabilizing external resistor [45].) Interestingly, it was also Steenbeck who was later to create the theoretical framework that nowadays allows for a deeper dynamical understanding of those somewhat vague thermodynamic principles. Mean-field magnetohydrodynamics (MHD) was originally developed to explain self-excitation of cosmic magnetic fields [46]. Its main idea is that certain correlations of the small-scale parts of velocity and magnetic field contribute to the dynamics of the large scale magnetic field [47]. For helical turbulence the authors introduced the celebrated α effect which drives an electromotive force parallel to a prevailing large scale magnetic field B. Similarly, turbulence leads to an increase of the resistivity by the β effect, so that the mean electromotive force can be written in the form E = αB − β∇ × B.
Nowadays, mean-field concepts play not only a role in dynamo theory but also in the description of magnetically driven instabilities. Flow-driven helical dynamos and magnetically dominated helical "dynamos" are presently considered as two different aspects of the very same mean-field MHD [48]. The detailed saturation mechanism of the TI, in particular its termination in a helical or non-helical state, is but one interesting application of mean-field MHD.
In this paper, we are going to study the exponential growth and the final saturation of the TI in finite cylindrical geometry for varying values of Ha and P m. On the basis of an axisymmetric (m = 0) base state with an homogeneous axial current J 0 that produces an azimuthal magnetic field B 0 , we compute the m = 1 TI-eigenmode comprising the velocity u and the induced magnetic field b from which we infer the mean electromotive force u × b (the overbar denotes the average over the azimuthal angle), and from this the mean-field coefficients α and β. As a product of two m = 1 modes, u × b comprises certain m = 0 components that drive an azimuthal current (by virtue of the α effect) and reduce the impressed axial current (by virtue of the β effect). Although there is not a big scale separation between the m = 0 base state and the m = 1 perturbation, mean-field theory perfectly applies here. The electromotive force E in direction of the mean field B (i.e. the α effect), will be interpreted in terms of its relation to the smallscale current helicity, E · B = −j · b/σ + e · b, that had been derived and utilized by different authors [49,50,51,52,53].
In case of S ≫ 1, the modified currents and fields could be expected to resemble the typical Bessel function structure as typical for Taylor relaxation. In this sense, the mean axial field produced by the α effect would follow from the principle of minimum dissipation [36].
However, this type of saturation mechanism, which relies on changing -by meanfield induction effects -the electromagnetic base state in such a way that it becomes just marginally stable against TI, does not apply for S ≪ 1. In this case the magnetic Reynolds number Rm of the TI-produced flow is much too small to induce any significant changes of the original applied magnetic field. The saturation must instead rely on a modification of the hydrodynamic base state, which we will discuss in detail. We will also evidence the occurrence of helicity oscillations, whose amplitudes and frequencies in dependence on Ha and P m we will characterize.
The paper closes with a discussion of the results, and with an outlook towards an application to stellar dynamo theory.

The numerical scheme
The usual numerical schemes for the simulation of TI, which solve the Navier-Stokes equation for the velocity and the induction equation for the magnetic field, are working typically only for values of P m down to 10 −3 , although in a recent work by Herreman et al this limit has been challenged [25].
Here, we circumvent the usual P m limitations of these codes by replacing the solution of the induction equation for the magnetic field by invoking the socalled quasistatic approximation [21]. We replace the explicit time stepping of the magnetic field by computing the electrostatic potential by a Poisson equation, and deriving the electric current density. However, in contrast to many other inductionless approximations in which this procedure is sufficient, in our case we cannot avoid to compute the induced magnetic field, too. The reason for that is the presence of an externally applied electrical current in the fluid. Computing the Lorentz force term it turns out that the product of the applied current with the induced field is of the same order as the product of the magnetic field (due to the applied current) with the induced current. Here, we compute the induced magnetic field from the induced current density by means of Biot-Savart's law. This way we arrive at an integro-differential equation approach, as it had already been used by Meir and Schmidt [22].
In detail, the numerical model as developed by Weber et al [5] works as follows: it solves the Navier-Stokes equations (NSE) for incompressible fluidṡ with u denoting the velocity, p the (modified) pressure, f L = J ×B the electromagnetic Lorentz force density, J the total current density and B the total magnetic field. The NSE is solved using the PISO algorithm and applying no slip boundary conditions at the walls. Ohm's law in moving conductors allows to compute the induced current j by previously solving a Poisson equation for the perturbed electric potential ϕ = φ − J 0 z/σ: In the following, we will concentrate on cylindrical geometries with an axially applied current. Then, after subtracting the (constant) potential part J 0 z/σ, with z as coordinate along the cylinder axis, we use the simple boundary condition ϕ = 0 on top and bottom and n · ∇ϕ = 0 on the mantle of the cylinder, with n as the surface normal vector. The induced magnetic field can then be calculated by Biot-Savart's law Since this is a costly procedure, we modify here the method of [5] slightly. Actually, equation (4) is applied only at the boundary of the cylinder, while the magnetic field in the bulk is computed by solving the vectorial Poisson equation ∆b = µ 0 σ∇ × (u × B) which results from the full time-dependent induction equation in the quasi-stationary approximation.
Knowing b and j we compute the Lorentz force f L for the next iteration. A flow chart of this numerical procedure is shown in figure 1. For more details about the numerical scheme, see section 2 and 3 of [5].

Results
For the sake of simplicity, we consider a cylindrical electrically conducting fluid with a ratio of height L to diameter 2R of 1.25. A current is applied from top to bottom, just by setting the electrical potential to constant (but different) values at the two faces. Note that we refrain from taking into account any currents in the electrodes at top and bottom which has been shown to lead only to minor modifications of the results [54]. The side walls of the cylinder are considered as electrically insulating. No-slip boundary conditions apply to the velocity at all boundaries.
In the following, we will focus on three different cases. The magnetic Prandtl number for the first two runs is P m = 10 −6 . Differing in the Hartmann number (60 and 100), these two runs will show a quite different behaviour of the helicity. Whereas for Ha = 60 the initially growing helicity ultimately goes to zero, for Ha = 100 we observe helicity oscillation in the saturation regime. For the third case, with the much higher P m = 10 −3 and Ha = 100, we will find a finite and non-zero value of the final helicity. At the end of this section, we will summarize the different ways of saturation.

Saturation with zero helicity
Here, we choose Ha = 60 and P m = 10 −6 which results in a Lundquist number S = 0.06 that is definitely low enough for applying the quasistatic approximation.   Figure 2a, to start with, shows the evolution of the averaged Reynolds number of the flow arising from the initial state at rest: Re = R/ν( u 2 ) 1/2 (here, ... denotes an average over the total volume rather than only over the azimuthal angle, which will always be indicated by an overbar). We chose a logarithmic scale in order to evidence the exponential growth of the TI. Approximately at t n = 0.27, the instability starts to saturate. We have added the respective energy contents of the various azimuthal wavenumbers m = 0, 1, 2. Evidently, the m = 1 mode is the dominating one throughout the evolution. However, approximately at t n = 0.17, both the m = 0 and m = 2 modes start to increase with the double growth rate as the m = 1 mode. These even modes result from the non-linear term of the NSE. Saturation sets in when the m = 0 and m = 2 modes have acquired an amplitude comparably to that of the m = 1 mode which ultimately brings the growth rate of the TI to zero. The corresponding evolution of the averaged induced magnetic field is depicted in figure 2b. Note that the m = 0 component is here significantly weaker than the m = 2 component, quite in contrast to the rather parallel evolution for the kinetic energy.
The kinetic helicity H u = u · (∇ × u) is the next quantity to be discussed (see figure 2c). Actually, we show here the helicity as normalized by the mean square of the velocity, i.e. H u n = u · (∇ × u) R/ u 2 . After an initial increase, this normalized helicity stays nearly constant for a while (i.e., the non-normalized helicity H u grows with the same growth rate as the kinetic energy), until it decays to zero when saturation is reached.
In Figure 2d we give an interpretation of the mean electromotive force E in direction of the large scale magnetic field in terms of the small-scale current helicity j · b, according to the relation E · B = −j · b/σ + e · b [49,50,51,52,53]. We see that this relation is perfectly fulfilled, showing that α is essentially proportional to the current helicity, with a minor correction coming from an electric field term (note that all quantities are normalized here by B 0 J 0 /σ). Figure 2e depicts separately the α effect, defined by α = (u × b) · B 0 /B 2 0 . Since α has the dimension of a velocity, we give it here in form of a magnetic Reynolds number which includes again a complete spatial average: α n = µ 0 σR (u × b) · B 0 /B 2 0 . We observe initially an (exponential) increase, though to a very small value in the order of 10 −10 , and then a decay to zero.
The induction effects of the mean field coefficients α and β are illustrated in figure  2g and 2h. Figure 2g shows the mean axial magnetic field b z which is produced by the azimuthal current that is driven, in turn, by α. Normalized to B 0 , we see again that this induction is negligibly small. Note also that b z , in contrast to α (figure 2e), does not completely vanish in the saturation regime. We attribute this to numerical inaccuracies which seem to "stray" some energy from the much stronger m = 1 mode into the m = 0 and m = 2 modes, an effect that is already visible in the non-physical parallelism of all modes in the beginning of the exponential growth regime (figure 2a). In contrast to this slight discrepancy, the behaviour of j z /J 0 (figure 2h) is nearly identical to that of β  For three selected instants in time, figure 3 illustrates the spatial structure of various features of the TI (the normalized values of u z , b z , (u × b) · B 0 , (u × b) · J 0 ). The left column depicts these quantities amidst the exponential growth phase in which we observe a clear helical structure. The middle and right columns show then the respective structures shortly before and during the saturated regime. Evidently, the helicity has completely disappeared here. In terms of the mode structure, we notice that the left handed spiral (m = 1, say) and the right handed spiral (m = −1, say) have grown to the same strength.
Having seen that, due to the low values of P m, neither α nor β is able to induce any relevant change of the electromagnetic base state that would lead to saturation of the TI, we have to look for an alternative saturation mechanism. Evidently, this can only be related to a change of the hydrodynamic state, i.e. the flow field. Let us return to figure 2a: after a long exponential growth period of the (more or less) pure m = 1 mode, at t n = 0.17 the m = 0 and m = 2 modes start to grow due to the action of the nonlinear terms in the NSE. As shown in figure 4, the m = 0 part of the velocity in the saturation comprises two poloidal vortices pointing outward in the equatorial plane (in the dynamo community this flow topology would be denoted by s2 + [55]). This axisymmetric state, together with the m = 2 component, changes the hydrodynamic base state for the TI so that it becomes just marginally stable. With the simultaneous appearance of an m = 0 and m = 2 component, it is no surprise that the saturation is also connected with a restoration of the chiral symmetry. According to the sum rule for the nonlinear interaction the m = 2 will produce, from any dominant m = 1 mode, a corresponding m = −1 mode, so that chiral symmetry is ultimately restored.

Saturation with helicity oscillations
We now increase the Hartmann number from Ha = 60 to Ha = 100. As before, figure 5 illustrates the time evolution of various quantities. While the behaviour of the Reynolds number and the β effect are not significantly different from the previous case with Ha = 60, the helicity and α look quite differently. Evidently, the TI does not anymore saturate with zero helicity, but instead produces a helicity oscillation in the final state. Figure 6 shows again the spatial structure of various quantities in the exponential growth phase and at two different instants in the saturated state. In the plots for (u × b) · B 0 , the helicity oscillation appears as a slightly changing asymmetry between positive (orange) and negative (blue) "blobs". For the range between Ha = 40 and Ha = 140 we characterize this helicity oscillation by its amplitude and frequency (figure 7).

Saturation with finite helicity
We leave now the realm in which the magnetic induction is irrelevant for saturation and move towards the parameter region which had already been explored by other codes. Actually, we increase the magnetic Prandtl number to P m = 10 −3 , and consider the case Ha = 100. This leads to a Lundquist number of S = 3.16 for which we come close to the edge of applicability of the inductionless approximation [25]. Again, figures 8 and 9 show the time evolution, and some snapshots, of various quantities. The main difference to the former runs at low P m is that now the helicity and α acquire non-zero values in the saturated state. This is also shown in the snapshots of figure 9. Evidently, we are now in a regime where the induced magnetic fields contribute significantly to the saturation mechanism. This involves also that the m = 0 component of the magnetic field (8b) becomes now comparable to the m = 2 part, quite in contrast to the former cases with P m = 10 −6 . This is illustrated in figure 10 which shows the dependence of the induced mean current j z /J 0 and the induced mean axial field b z /B 0 on P m (at fixed Ha = 100). According to the criterion of Kruskal-Shafranov [14], we know that b z tends to inhibit the TI, so that a finite value of α is indeed likely to appear in the saturation regime.

Between chiral symmetry breaking and helicity oscillations
In the following, we will summarize the three different saturation mechanisms. We will be guided by the simple and instructive model of chiral symmetry breaking model that had been worked out by Bonanno et al [4].
The authors started from left and right handed helical modes for the velocity and the magnetic field, fulfilling the Beltrami relation which can be realized by appropriate linear combinations of Chandrasekhar-Kendall functions J m (r λ 2 + n 2 π 2 /h 2 ) cos(mφ) cos(nzπ/H). Invoking some symmetry arguments, the authors "guessed" the simplest form of a Lagrangian which then leads to the following evolution equations for the energy of the left and right handed helical modes The last terms on the r.h.s. of equations (6,7) describe the so-called mutual antagonism between the two chiralities that has been used extensively in the theory of homochirality of bio-molecules [56,57,28].
From here, one arrives at the following evolution equations for the total energy E = E R + E L and the helicity H = E R − E L : In figure 11 we illustrate the phase portrait of this equation system showing a clear chiral symmetry breaking both in the exponential growth phase as well as in the  Figure 11. Phase portrait of the coupled equation system (8,9), for the parameters γ = 2.71, µ = 3.0, µ * = 5.7. In this typical situation, S 1 is a repeller, S 2 and S 3 are attractive points, while S 4 is a saddle point.
the phase portrait of the current helicity (right column). To make contact with figure  11, we have now chosen a different normalization so that both helicities start at zero.
We start, in the first row of figure 12, with P m = 10 −6 and Ha = 60. Evidently the exponential growth phase looks very similar as in figure 11, but ultimately the system runs into a state with zero helicity.
The second case with P m = 10 −6 , Ha = 100 is similar but terminates with a helicity oscillation around zero. It is remarkable that this helicity oscillation proceeds without any significant oscillation of the energy. The third case, P m = 10 −3 , Ha = 100 has indeed some resemblance with the above model of Bonanno et al and terminates with a finite, non-zero helicity.
In the fourth row, we add here also the plot for P m = 10 −2 and Ha = 100 which amounts to S = 10 which is beyond the simple applicability of the quasistatic approximation and should be, therefore, considered with caution. Evidently, the magnetic helicity shows now a phase portrait that is similar to figure 11, despite the fact that we now observe an "overshoot" to the helicity with the other sign, and also a remaining oscillation in the saturated state. While we do not claim that this is a reliable result, the two last rows of figure 11 at least suggest that we are now approaching the usual saturation scheme as already discussed by [4].

Conclusions
In this paper, we have utilized an integro-differential equation solver for addressing the problem of chiral symmetry breaking in the exponential growth phase and in the saturation phase of the Tayler instability. The advantage of this code is its easy applicability for small magnetic Prandtl numbers, while its suitability for problems at S > 1 is at least questionable (at least it has to be checked case by case whether a final steady state with finite and non-oscillatory, i.e. static, helicity could still be treated with our quasistatic scheme). Our simulation have allowed to identify three different saturation regimes.
To start with the last regime, for a comparably large value of S = 3.16 we have confirmed a similar type of chiral symmetry breaking as it was previously evidenced by Gellert et al [3] and Bonanno et al [4]. Depending on the random initial conditions, the TI grows with one of the two possible helicities which does not disappear in the saturated regime. The helicity is intrinsically connected with a non-zero α effect that generates a current parallel to the applied azimuthal magnetic field. At the same time the mean-field e.m.f. contains also a significant β effect that changes the axial current. Both effects together work against the TI. In the ultimate case of high S (which is, probably, not accessible by our code) one could expect a sort of Taylor-relaxation into a helicity maximizing state [33]. Whether for those large values of S one reaches a regime of helicity oscillation around a finite value (as suggested by the fourth row of figure 12) is still to be validated by complementary codes.
The saturation mechanism, which relies on changing (by mean-field induction effects) the electromagnetic base state in such a way that it becomes just marginally stable against TI, does not apply for S ≪ 1. This can already be seen from the general scaling Re ∝ Ha 2 , or Rm ∝ S 2 which applies both to S < 1 and S > 1 [5]. For S ≪ 1 the final Rm becomes much too small in order to induce any significant changes of the original applied magnetic field.
In this parameter region, the saturation relies instead on the non-linear appearance of an m = 0 and an m = 2 velocity component, which now changes the hydrodynamic base state of the TI in a way that the growth rate of the TI vanishes.
Perhaps the most interesting result of our study is the observation of a robust and systematic helicity oscillation whose amplitude and frequency dependence on Ha has been worked out. Interestingly this helicity oscillation is not connected with any significant energy oscillations.
Based on the latter observation, we would like to conclude with some, admittedly, very speculative ideas. There is a long tradition in trying to link the various frequencies of the solar dynamo action to corresponding periodicities of planetary motion. Tracing back to a paper by Jose [58], who had related the 11.86 years Jupiter orbit with the 22.08 year solar cycle, some refinements of this connection in terms of a combined torque and gravity action of Earth, Venus and Jupiter have been discussed recently [59]. Other papers have tried to link periodicities of the Jupiter-Saturn orbits to various longer-time cycles of the solar dynamo [60,61], with possible connections to the climate on Earth. However, in all those cases, it was noticed that the planetary forces are much to weak to compete with the typical acceleration forces in the convection zone. The only viable way of influencing the solar dynamo was speculated to rely on the action of gravity on the shape, or local rotation rate, of the tachocline. Yet, this would imply that the solar dynamo works indeed as some sort of Tayler-Spruit dynamo [1], in which the transformation from poloidal to toroidal field is traditionally realised by differential rotation, the reverse mechanism, however, by some α effect due to the TI.
It is exactly here where helicity oscillations, and their possible synchronization with planetary forces and torques, might come into play. In particular, since the oscillations of α are not connected to any significant changes of the energetic content, very minor changes of the state of the tachocline might just open the "α-bottleneck" for the Tayler-Spruit dynamo (which is, in any case, still a sort of α − Ω like dynamo). Even if the α oscillations are around some non-zero mean values prevailing in the two solar hemispheres, it still might give rise to dynamo oscillations.
Note that a parametric resonance and synchronization of m = 1 dynamo eigenmodes with m = 2 velocity perturbations has been observed both for galactic dynamos (swing excitation, [62]), and for a VKS-like dynamo [63]. Whether a similar effect may actually be at work for synchronizing the solar dynamo with periodic planetary forces via their action on the tachoclinic state, will remain a topic for future investigations.