Triadic resonances in non-linear simulations of a fluid flow in a precessing cylinder

We present results from three-dimensional non-linear hydrodynamic simulations of a precession driven flow in cylindrical geometry. The simulations are motivated by a dynamo experiment currently under development at Helmholtz-Zentrum Dresden-Rossendorf (HZDR) in which the possibility of generating a magnetohydrodynamic dynamo will be investigated in a cylinder filled with liquid sodium and simultaneously rotating around two axes. In this study, we focus on the emergence of non-axisymmetric time-dependent flow structures in terms of inertial waves which - in cylindrical geometry - form so-called Kelvin modes. For a precession ratio ${\rm{Po}}=\Omega_p/\Omega_c=0.014$ the amplitude of the forced Kelvin mode reaches up to one fourth of the rotation velocity of the cylindrical container confirming that precession provides a rather efficient flow driving mechanism even at moderate values of ${\rm{Po}}$. More relevant for dynamo action might be free Kelvin modes with higher azimuthal wave number. These free Kelvin modes are triggered by non-linear interactions and may constitute a triadic resonance with the fundamental forced mode when the height of the container matches their axial wave lengths. Our simulations reveal triadic resonances at aspect ratios close to those predicted by the linear theory except around the primary resonance of the forced mode. In that regime we still identify various free Kelvin modes, however, all of them exhibit a retrograde drift around the symmetry axis of the cylinder and none of them can be assigned to a triadic resonance. The amplitudes of the free Kelvin modes always remain below the forced mode but may reach up to 6% of the of the container's angular velocity. The properties of the free Kelvin modes will be used in future simulations of the magnetic induction equation to investigate their ability to provide for dynamo action.


Introduction
Instabilities and waves in rotating fluids are important in numerous technical applications. In most cases it is important to avoid these phenomena in order to ensure the stability of fast spinning liquid filled bodies like gyroscopes (Pfeiffer 1974), spacecrafts (Bao & Pascal 1997) or projectiles (Miller 1982). Waves in rotating fluids are also of general interest because of their fundamental character in oceanographic and atmospheric flows (Rossby 1939) like, e.g., zonal flows in the atmospheres of Saturn or Jupiter that may develop via a modulation instability of meridional Rossby waves (Connaughton et al. 2010). Moreover, the identification of inertial waves in the Earth's fluid outer core (Aldridge & Lumb 1987) allows conclusions on the core dynamics and the associated dynamo process that is responsible for the generation of the Earth's magnetic field. Nowadays, it is believed that in general planetary dynamos are driven by rotating thermal or compositional convection (King et al. 2010, M. Breuer et al. 2010, D. Breuer et al. 2010). On the other hand, other types of flows cannot be ruled out, and, for example, a superposition of random inertial waves in a rotating conducting fluid is capable of transferring energy in a magnetic field as well (Moffatt 1970). However, in the model of Moffatt (1970), field and flow must decay in the long term because of the lack of energy sources required for a steady driving. This problem can be overcome by more complex models in which, for example, the equatorial heat flux in the Earth's outer core provides a persistent excitation mechanism for inertial waves with sufficient helicity for the generation of planetary magnetic fields (Davidson 2014). Inertial waves can be excited by precessional forcing as well, and it has long been discussed whether precession of the Earth can provide the necessary power to drive the geodynamo (Stewartson & Roberts 1963, Busse 1968, Malkus 1968. In a next-generation dynamo experiment currently under development at Helmholtz-Zentrum Dresden-Rossendorf (HZDR) a precession driven flow of liquid sodium will be used to validate its suitability to drive a dynamo (Léorat et al. 2003, Léorat 2006, Stefani et al. 2012. Since the precessional forcing provides a natural driving mechanism, it represents a qualitatively new approach compared to previous dynamo experiments, where an optimized flow was driven either by impellers (Gailitis et al. 2000, Monchaux et al. 2007) or by a system of electromagnetic pumps (Stieglitz & Müller 2001).
A small-scale water experiment is currently running at HZDR in order to examine the hydrodynamic properties of a precession driven flow in a cylindrical container with strong forcing and a large precession angle. Measurements of pressure fluctuations and of axial velocity profiles show three distinct flow regimes in dependence of the precession ratio: a laminar state (with only the forced mode) followed by a non-linear regime with various modes superimposed on the fundamental mode, eventually leading to a permanent chaotic state when a critical precession ratio is exceeded (Herault et al. 2015). Presumably, the ability of the flow to drive a dynamo will be different in the three regimes and although the precessional dynamo experiment at HZDR will allow magnetic Reynolds numbers of the order of Rm ≈ 700 (based on the rotation velocity and the radius of the cylindrical container) it is a priori not obvious whether a magnetic field will be self-generated at all. There are promising indications for dynamo action driven by precession both from liquid sodium experiments by Gans (1971) which -despite of the rather small precession ratio -achieved an amplification of an applied field by a factor of three, as well as from simulations in spheres (Tilgner 2005), spheroids (Wu & Roberts 2009), ellipsoids (Ernst-Hullermann et al. 2011), cubes (Krauze 2010) and cylinders (Nore et al. 2011). However, kinematic dynamo simulations demonstrate that the primary forced mode alone, which consists of a single, non-axisymmetric mode with azimuthal wave number m = 1 and axial wave number k = 1, cannot drive a dynamo at a reasonable Rm . This coincides with Herreman & Lesaffre (2011) who showed that a linear wave packet of inertial waves cannot drive a mean-field dynamo. However, this does not generally preclude the possibility of quasi-laminar (in the sense of Dudley & James 1989) or small scale dynamo action driven by inertial waves or another precession induced instability with sufficient intricate topology.
Experimental studies in a weakly precessing cylinder filled with a non-conducting fluid showed a large diversity of flow phenomena with complex structures frequently resulting in a vigorous turbulent state (Gans 1970, Manasseh 1992, Manasseh 1994, Kobine 1995, Manasseh 1996. Other experiments revealed large scale structures like a system of intermittent, cyclonic vortices that may provide a strong source of helicity and hence be beneficial for a dynamo (Mouhali et al. 2012), or free Kelvin modes with an azimuthal wave number m = 5 and m = 6 propagating around the symmetry axis of a weakly precessing cylinder , Lagrange et al. 2011.
Free Kelvin modes are the natural eigenmodes in a rotating cylinder (Thomson 1880, Greenspan 1968, McEwan 1970 and are promising candidates for driving quasilaminar dynamo action in case of precessional forcing because at least a subclass of these modes has a structure similar to the columnar convection cells that are responsible for dynamo action in convection driven models of the geodynamo (Lorenzani & Tilgner 2003). Free Kelvin modes may emerge from non-linear interactions of a forced Kelvin mode with itself or as a parametric instability that involves the forced mode and two free Kelvin modes (Kerswell 1999, Lagrange et al. 2011). The latter case requires appropriate combinations of wave numbers, frequencies and a geometry such that all three modes become resonant simultaneously which usually is named a triadic resonance. Triadic resonances have been observed in simulations of precession driven flow in a spheroid (Lorenzani & Tilgner 2003) and, experimentally, in weakly forced precession in a cylinder ) and a cylindrical annulus (Lin et al. 2014).
In the present study we conduct numerical simulations of a precessing flow far below the transition to the chaotic state observed in water experiments, but with sufficient forcing to excite free Kelvin modes with azimuthal wave numbers m > 1. We focus on the impact of the aspect ratio which -at least in the linear approximation -determines whether an inertial wave becomes resonant. We identify different free Kelvin modes by Figure 1. Sketch of the problem setup. The cylinder rotates with Ω c and the rotation axis precesses with Ω p . In the simulations we varied the aspect ratio Γ = H/R and fixed the angle between rotation axis and precession axis at α = π/2. means of their spatial structure, and we compute their frequencies from the azimuthal drift motion in order to conclude whether the free Kelvin modes constitute a triadic resonance. The results are used to further constrain the modes that may be observed in the water experiment at HZDR ‡ and will be applied in future simulations of the magnetic induction equation in order to proof whether free Kelvin modes are suitable to drive a dynamo.

Theoretical background
The flow of a fluid with kinematic viscosity ν in the frame of an enclosed precessing cylinder is described by the Navier-Stokes equation including terms for the Coriolis force and the Poincaré force (see e.g. Tilgner 1998): Here, u denotes the velocity field, which additionally obeys the incompressibility condition ∇ · u = 0, r is the position vector, P is the reduced pressure (including the centrifugal terms), Ω c is the angular frequency of the cylindrical container and Ω p denotes the precession (see figure 1). The precession axis is time dependent in the cylinder frame, and the unit vector Ω p that denotes the orientation of the precession ‡ Pressure measurements in the water experiment at HZDR show a periodic signal with two frequencies close to the values expected from the dispersion relation for free Kelvin modes (see below). However, the measurements do not yet allow a unique identification of azimuthal, axial or radial wave numbers. axis is given by with α the angle between rotation axis and precession axis. We consider a cylinder with aspect ratio Γ = H/R, where H and R are height and radius so that r ∈ [0; R] and z ∈ [0; ΓR]. In the following we non-dimensionalize all quantities using R as length scale and Ω −1 c as time scale. In the inviscid limit a linear solution that fulfills the boundary conditionsẑ · u = 0 at the top and the bottom is given by a Kelvin mode (Thomson 1880, Greenspan 1968): U k,m,n (r, z, ϕ, t) = u k,m,n (r, z)e i(Ω k,m,n t+mϕ) + c.c. ( with a frequency Ω k,m,n that depends on the azimuthal wave number m, the axial wave number k and a third number n that counts the roots of the dispersion relation in which J m denotes the Bessel function of order m and the triple index (k, m, n) is replaced by j. Positive (negative) frequencies correspond to retrograde (prograde) propagation. Both signs correspond to different radial wave numbers (different solutions of the dispersion relation) and hence to a different spatial structure of the corresponding Kelvin mode. The dispersion relation ensures the fulfillment of the radial boundary conditionŝ r · u = 0 at the sidewalls and although λ j is not an integer it plays a role similar to a radial wave number with its position in the sequence of zeros corresponding to the number of half-cycles in the radial direction. The spatial structure of a Kelvin mode in a meridional plane u j (r, z) is given by with λ j and Ω j taken from the solutions of the dispersion relation (4).
In the inviscid linear approximation the total flow of the forced mode with m = 1 is given by where we replaced the index j by the individual contributions and A 1kn represents the amplitude of the forced mode which can be computed by a projection of the applied forcing onto a particular mode (Manasseh 1994). For a precession ratio (or Poincaré number) Po = Ω p /Ω c and a precession angle α the corresponding calculation yields (Liao & Zhang 2012) If the wave length that corresponds to an axial wave number of a mode matches exactly the height of the cylinder, this mode becomes resonant with an eigenfrequency Ω j = 1 and the expression for the inviscid amplitude (7) diverges. For each Γ there exist, in principle, an infinite number of resonant modes. From the dispersion relation (4) we find that the primary forced mode, i.e. the mode with m = 1, the axial wave number k = 1 and a radial wave number corresponding to the first root, becomes resonant at Γ = 1.9898174 which is rather close to the so called spherical cylinder with height equals diameter (Γ = 2). Further resonances for increasing radial wave number occur at Γ = 0.9560735 (for n = 2) and at Γ = 0.6206981 (for n = 3). The computation of the amplitude at resonance requires a consideration of viscous effects in the bulk and in terms of boundary layers with associated Ekman layer suction. Then the amplitude can be calculated by matching the Ekman layer suction to the precessional force (Gans 1970). A more general approach for the computation of the amplitudes of the response created by the precessional forcing is given by Liao & Zhang (2012) who derived an asymptotic solution including viscous effects that is valid at and away from resonance.
The precessional forcing only excites modes with m = 1 and k odd. Higher azimuthal modes or modes with even axial wave number must be triggered by nonlinear interactions, e.g., in terms of triadic resonances that involve this forced Kelvin mode and two free Kelvin modes (m a , k a ) and (m b , k b ) or by the non-linear interaction of the forced Kelvin mode with itself . In the following we assume that the resonant case is most promising in view of the dynamo problem due to the expected larger amplitudes of the velocity field.
The free Kelvin modes are solutions of the linearized Navier-Stokes equation with a Coriolis term as a restoring force but without a precessional driving on the right hand side: The principle of the interaction of free and forced modes becomes evident when we consider an exact triad with a forced mode U f , with m f = 1 and k f = 1 and two free Kelvin modes U a with m a and k a and U b with m b and k b . We ignore viscosity and further assume a constant amplitude A f of the forced mode so that we can write the total flow as where a(t) and b(t) denote the amplitude modulations of the free mode (without the harmonic oscillations Ω a and Ω b ) which are given by Kerswell (1999) The square brackets ·, · denote a scalar product defined as U a , U b = U * a · U b dV , and a normalization of the free modes is assumed such that U i , U j = δ ij . In order to achieve non-trivial solutions of (10), and to constitute a triadic resonance, the wave numbers and frequencies of the free modes a and b must fulfill the conditions where m f = 1, k f = 1 and Ω c = 1 denote the properties of the forced mode. In contrast to sine and cosine-functions, there are no corresponding addition theorems for Bessel functions so that no further restrictions are imposed on the interaction of different radial wave numbers.
In the following, we solve the dispersion relation (4) for k = 1, 2 and m = 1, . . . , 7. We assume that higher radial modes are damped and only consider the first root for each combination of m and k. Figure 2 shows the frequency difference δΩ of two free Kelvin modes with δm = 1 and δk = 1 versus the aspect ratio in the vicinity of the critical value δΩ = 1. The crossings of the curves with the horizontal line at δΩ = 1 mark aspect ratios for which an exact triadic resonance may occur. These aspect ratios at which δΩ = 1 are listed in table 1 together with the corresponding frequencies of the resonant free Kelvin modes.   Table 1. Possible triadic resonances with δm = 1, δk = 1 and δΩ = 1 and corresponding aspect ratios. Negative frequencies present pro-grade modes and positive frequencies denote retrograde modes. Only the first radial wave number has been considered for each value of m and k.
A triad is always formed by free Kelvin modes consisting of one prograde and one retrograde mode (Waleffe 1993). Note that the frequencies of free Kelvin modes that constitute different triads often are close (if we restrict to triads with k a = 1 and k b = 2, see table 1) making it difficult to reliably identify individual Kelvin modes in realistic setups with potential frequency shifts due to viscous and/or non-linear effects.

Numerical model
We conduct numerical simulations in a cylindrical geometry with the code SEMTEX that applies a spectral element Fourier approach for the numerical solution of the Navier-Stokes equation (Blackburn & Sherwin 2004). In order to simplify the implementation of the precessional forcing we switch to the precessional frame in which the precession is stationary and the cylinder rotates with a frequency Ω c . Thus only a term for the Coriolis force appears in the Navier-Stokes equation which now reads Due to the rotation of the cylinder in the precessing frame the boundary conditions for the azimuthal flow change to u ϕ = rΩ c (at endcaps and sidewall) whereas no-slip conditions are still applied for u z and u r . The problem is described by four parameters, the Reynolds number defined with the angular velocity of the container Re = Ω c R 2 /ν, the precession ratio (or Poincaré number) Po = Ω p /Ω c , the aspect ratio Γ = H/R and the precession angle α with cos α = Ω p · Ω c /|Ω p ||Ω c |. In the present study, we keep the precession ratio fixed at Po = 0.0141 and the precession angle is set to α = π/2. All simulations are performed at Re = 6500 except for one run at Re = 10000 to briefly examine the impact of increasing Re. The aspect ratio is varied in the range Γ ∈ [1.7; 2.4] in order to find the container geometry for which individual modes or triads become resonant.

Pattern of the total flow
The simulations are started from an initial state with pure solid-body rotation u = r e ϕ . After switching on the Coriolis force at t = 0 the fluid shows a direct response in form of a large scale mode with m = 1. In dependence of the aspect ratio this is only a transient state until after a finite time (which may take up to 1000 rotation periods or more) a quasi-steady state is reached. The typical flow pattern for Γ = 2 and Re = 6500 is presented in figure 3 where the nested isosurfaces show a snapshot of the axial velocity u z at 30% (60%, 90%) of its maximum value. The flow is clearly dominated by a Volume rendering of the axial velocity u z at Re = 6500 and Γ = 2.
The nested isosurfaces show u z at 30% (60%, 90%) of its maximum value. An animation of the flow's behavior in the statistically stationary state is available at https://www.hzdr.de/db/VideoDl?pOid=45098. velocity mode with m = 1 and k = 1. This is the primary Kelvin mode driven by the precessional forcing. Superimposed on the m = 1 mode we see some non-axisymmetric (time-dependent) disturbances that will be analyzed below.

Kinetic energy
4.2.1. Forced mode In the following we discuss the behavior of the flow by means of azimuthal Fourier modes u m (r, z) with the total flow given by u(r, z, ϕ) = 2π/(N ϕ + 1) Nϕ m=0 u m (r, z)e imϕ . Figure 4 shows the history of the kinetic energy of the m = 1 mode for different aspect ratios. A final state arises either after a short transitional period (Γ = 2.000) or after one or more bifurcations (see e.g. the blue curve for Γ = 2.150). The final state can be chaotic (Γ = 2.000, black curve), oscillatory (Γ = 2.15, blue), quasi-periodic with collapses (Γ = 2.2, red), quasi-periodic with bursts (Γ = 2.225, green) or stationary (Γ = 1.8, yellow). We also find transient periods with constant energy (e.g. Γ = 2.15 for t = 300...600), and it may take up to 1200 revolutions of the cylinder (Γ = 2.15, blue curve) till a final quasi-stationary state is reached. This saturation time also depends on the Reynolds number and decreases with increasing Re.
The bifurcations shown by the m = 1 mode go along with the emergence of nonaxisymmetric modes with m ≥ 2 (figure 5). These modes become unstable only after quite a long time, e.g., for Γ = 1.825 it takes up to 600 rotation periods of the cylinder till modes with m ≥ 2 start to grow, and it may take up to twice that time till a quasisteady state is reached. The non-axisymmetric modes with m ≥ 2 essentially proceed similar to the behavior of the m = 1 mode, i.e., we see a chaotic behavior when the m = 1 mode behaves chaotic (Γ = 2.0, figure 5, central panel) and we see periodic growth and decay when the m = 1 mode is periodic (Γ = 1.825, top panel in figure 5 and Γ = 2.200, bottom panel in figure 5). In the periodic regime the amplitude of the energy oscillations of the modes with m ≥ 2 is rather large so that the kinetic energy may vary in time by more than one order of magnitude.
In the following we compare analytic expressions for the energy of the forced mode with the time-averaged energy taken from the quasi-stationary period of the final state in our simulations. This comparison is only partly justified because the analytical expressions are based on a number of serious simplifications, such as the linearization of  the Navier-Stokes equation, the limit Po → 0, and the neglect of the time dependent part of the Coriolis force which is justified only for small precession angles α. Furthermore, the final state in our simulations is instable and additional caution is advised when comparing simulations with a linear time independent solution. Nevertheless, such an analysis is helpful to identify the regimes in which, e.g., a maximum response of the flow can be expected, or to localize the aspect ratios at which triadic resonances are possible.
Analytically, the energy of the forced Kelvin mode is given by (Liao & Zhang 2012) where the index j now denotes a combination of the axial wave number k and the radial wave number n, and A j is the inviscid linear amplitude given by (7) with the azimuthal wave number fixed at m = 1. A j can also be calculated by the method of Liao & Zhang (2012) which includes viscosity and is valid at and off the resonance. However, in that case equation (13) is accurate only up to an error of O(Re − 1 2 ) since equation (13) does not consider the energy from the flow in the boundary layers. Figure 6 shows the time-averaged energy taken from the simulations (black curve) compared with the analytic solutions (green curve, inviscid solution and blue curve, viscous solution). For the calculation of the analytical solutions equation (13) is truncated at n = 6 and k = 6 which is sufficient to reach convergence (off-resonance in the inviscid case) and avoids the occurrence of too many resonances in the inviscid case (green curve) which anyway  Figure 6. Time-averaged kinetic energy of the fundamental mode (m = 1) versus aspect ratio. Re = 6500, Po = 0.014, α = π/2. The dashed red lines denote the resonances with larger k and/or n. For the computation of the energy, equation (13) has been truncated at k = 6 and n = 6 (with prograde and retrograde modes respectively). Note that the inviscid solution (green curve) diverges at the resonances whereas the viscous solution has a finite maximum of E max kin ≈ 0.021 at Γ ≈ 1.9814.
are suppressed in the more realistic viscous or non-linear computations. Comparing the numerical solutions (black curve) with the linear solutions we see significant deviations around the primary resonance (k = 1, n = 1, Γ res ≈ 1.9898) while outside of this regime the agreement of all curves remains good (say for Γ < ∼ 1.85 and for Γ > ∼ 2.05). Two key features characterize the behavior of the mean energy in our simulations. On the one hand we see a clear shift of the maximum to smaller aspect ratios and the maximum energy arises rather close to the aspect ratio expected for the resonance with k = 3 and n = 3. However, this resonance is suppressed by viscous effects and we do not see strong contributions with k = 3 and/or n = 3 in our simulations (see also section 4.5.1). Therefore the correlation is either coincidental or we have to assume an indirect impact of the (k = 3, n = 3) resonance. The second scenario is supported by the abrupt transition from the absolute maximum to a rather low energy state right below Γ ≈ 1.871. This jump is also connected to a modification of the character of the m = 1 mode which changes its behavior below Γ ≈ 1.86 from the chaotic regime to the periodic regime. Abrupt transitions of the energy of higher azimuthal modes can also be seen at the other (theoretical) resonances of the forced mode (k = 5, n = 5 and, a little further away k = 5, n = 4, see figure 6) which will be revisited in the next paragraph. Beside the shift of the maximum to smaller aspect ratios the energy obtained in the simulations remains considerably smaller compared to the linear viscous solution. The maximum of the energy based on the amplitude of Liao & Zhang (2012) is E max kin ≈ 0.021 (at Γ ≈ 1.9814) which is roughly two times larger than the maximum kinetic energy obtained in our simulations (E max kin ≈ 0.0093 at Γ = 1.871). The differences are not surprising because the prerequisites for the computation of the amplitudes A j are not well met in our simulations. In particular, we observe a significant impact of nonlinear interactions in the simulations in terms of higher non-axisymmetric modes and the forming of an azimuthal shear flow that both draw energy from the forced m = 1 mode.

Higher azimuthal wave numbers
We only find higher non-axisymmetric modes with a noteworthy amount of energy between Γ ≈ 1.81 and Γ ≈ 2.24. This corresponds approximately to the regime in which the energy of the m = 1 mode deviates from the linear prediction (see previous paragraph). For the sake of brevity we will call this regime the non-linear regime, assuming that the higher modes are triggered by non-linear effects. The kinetic energy of the modes with m ≥ 2 is roughly 2 orders of magnitude smaller than the energy of the forced mode but essential characteristics such as the maximum at Γ = 1.871 and the abrupt drop of kinetic energy below this threshold are also reflected in the behavior of the higher non-axisymmetric modes (figure 7). The  sudden drop of the energy for Γ ≤ 1.871 pertains to all modes which indicates that this breakdown might be a global flow property that is not connected to an individual mode. A similar steep decrease of energy can also be found on the two outer edges of the non-linear regime. It is striking that all the abrupt transitions are close to theoretical resonances of the forced mode (k = 3, n = 3 ⇒ Γ res ≈ 1.8621, k = 5, n = 5 ⇒ Γ res ≈ 1.8142, and with some greater distance k = 5, n = 4 ⇒ Γ res ≈ 2.2911 , see figure 6). These resonances are, however, suppressed and cannot be directly observed in the simulations. So far, we cannot conclude whether the correlation of the energy jumps and the occurrence of (theoretical) resonances is a coincidence, or whether they may be some systematic relationship.
Away from the primary resonance of the forced mode we find two further regimes with local maxima of the higher modes. For 1.80 < ∼ Γ < ∼ 1.85 we see a narrow window with the maximum around Γ ≈ 1.825 in which the m = 6 and m = 7 mode dominate (blue and yellow curve in figure 7). Likewise, but more explicit, we find the m = 4 and m = 5 mode becoming dominant for 2.05 < ∼ Γ < ∼ 2.25 with a maximum around Γ ≈ 2.2 (green and red curves in figure 7).
The local maxima around Γ ≈ 1.825 and Γ ≈ 2.2 for Fourier modes with δm = 1 indicate the presence of corresponding triadic resonances and roughly agree with predictions from the dispersion relation (see table 1 and figure 2). From the linear predictions we would also expect triads with m = 5 and m = 6 in the intermediate regime (for 1.85 < Γ < 2.05) but obviously there are other flow contributions with at least comparable energy which may disguise the energetic signature of the involved free Kelvin modes (e.g., we see relatively high energies of the m = 4 and the m = 5 modes, whereas in comparison the m = 6 mode is clearly suppressed).

Spatial Structure of the Fourier modes
A qualitative impression of the flow structure is provided by isosurfaces of the individual Fourier modes of the axial velocity shown in figure 8 for Γ = 2 (for an impression of the temporal behavior see the snapshots of the m = 5 mode in figure 9 and the movie at https://www.hzdr.de/db/VideoDl?pOid=45104). The m = 1 mode is obviously dominated by a k = 1 contribution and shows comparatively small variations in time. Regarding the higher non-axisymmetric modes we see a complex pattern for the m = 2 and m = 3 mode with irregular temporal and spatial fluctuations and a more regular behavior of the modes m = 4, 5 and 6 with typical signatures of k = 1 and k = 2. It is striking that the m = 6 mode, which should become resonant around Γ ≈ 2, remains weak compared to the m = 4 or the m = 5 mode.
The amplitude of the individual contributions with m = 4, 5, 6 varies in time and the resulting flow perpetually changes its axial structure. The azimuthal mode resulting from the superposition of k = 1 and k = 2 contributions has a distinct and time dependent asymmetry with respect to the equatorial plane of the cylinder, so that the associated flow is concentrated alternately in both halves of the cylindrical container (see time series of the m = 5 mode in figure 9). The typical timescale for this process is of the order of the rotation period 2π/Ω c but we can not derive a reliable periodicity from our simulations.
Away from the primary resonance of the forced mode the structure of higher azimuthal modes is much more regular and -except for periodic variations of amplitude -essentially remains constant over time. Figure 10 shows the dominant azimuthal modes (beyond m = 1) at Γ = 1.825 (m = 6 and m = 7) and at Γ = 2.200 (m = 4 and Figure 8. Snapshot of the axial velocity field u z filtered to m = 1 . . . 6 (upper left to lower right) for Γ = 2 (isosurface at 50% of the respective maximum value). An animation of the flow pattern is available at https://www.hzdr.de/db/VideoDl?pOid=45104. The temporal behavior of the modes m = 2 and m = 3 is complex and chaotic. The modes m = 4, 5 and 6 show a more regular behavior essentially resulting from the superimposition of free Kelvin modes with k = 1 and k = 2. Re = 6500, Po = 0.014, α = π/2. Figure 9. Time series of the m = 5 mode that visualizes the axial symmetry breaking. The isosurfaces show u m=5,z (r, ϕ, z) at 50% of its maximum value at various characteristic time steps (form upper left to lower right: t=2617.5418, 2617.8489, 2618.0845, 2618.1958, 2618.8738, 2618.9599 (time in units of the rotation time); note that the time gaps are not equidistant). m = 5). In both cases, the axial behavior of the modes shows clear indications of axial wave numbers k = 1 and/or k = 2 and we will show below that these modes indeed fulfill all conditions of equation (11) and thus constitute triadic resonances. In contrast to the behavior at Γ = 2, the geometric structure of the dominant modes exhibits only minor temporal variations except strong oscillations of the amplitude (see animation at https://www.hzdr.de/db/VideoDl?pOid=45105).

Radial dependence
In order to quantify the visual observations made in the previous section and to estimate amplitude and frequency of individual Kelvin modes we use a discrete sinetransformation (DST) applied to each Fourier mode with N z the number of points in axial direction. Figure 11 shows radial profiles for the dominant contributions to the axial velocity at Γ = 1.825 (a) and at Γ = 2.200 (b) where we only present the three leading modes for each case. The solid thick curves present the time average and the bright colors in the background represent the variation of the radial profiles in time. For comparison, the thick dashed curves show the analytical behavior ∝ J m (λ j r) predicted by the linear in-viscid approximation (5) with λ j the first root of the dispersion relation (4). The differences between the numerical data and the analytical curves are surprisingly small with deviations mainly in the boundary layers which can be explained by the different boundary conditions for the simulations (no slip) and the in-viscid solution (free slip). Except for (considerable) changes in the amplitude the pattern of the radial profiles exhibits only minor temporal changes and shows no significant contributions from higher radial modes (see color coded contour plots in figure 11 that show the variations of the radial profiles in time).  (5)). The colored contour plots below each profile show the corresponding variation of the radial profiles in time.
The behavior becomes more complex at Γ = 2 at which we see many modes with different m and k with nearly the same amplitude and clear changes in the radial structure. Figure 12 shows characteristic radial profiles for the mode (m = 5, k = 1) and their temporal behavior as a typical example. The amplitude exhibits considerable fluctuations and, unlike in the previous cases, we see additional changes in the radial structure. When the amplitude is around a (local) maximum, the radial velocity profile is in accordance with a wavenumber that corresponds to the first root of the dispersion relation and the radial profiles can nicely be described by a function ∝ J 5 (λ 1 r) (see red curves in figure 12a). However, when the amplitude is weak we see clear signatures of higher radial wave numbers (see figure 12 b and d). Even though the higher radial modes remain weak, they might represent a channel for the transfer of energy from the (periodically occurring) maxima of a mode into smaller scales. Figure 12. Snapshots of radial profiles u k m,z (r) with m = 5 and k = 1 at Γ = 2. (a) snapshots taken at local maxima of the (m = 5, k = 1) mode dominated by n = 1 (marked by the the vertical red lines in panel (c)). The blue curve represents the analytical behavior ∝ J 5 (λ j r) with λ j the first root of the dispersion relation. (b), (d): snapshots of the radial profiles that match to a higher radial wave number. The red curves denote the numerical data and the blue curves denote a fit ∝ J 5 (λ j r) with λ j corresponding to the second, third and fourth root of the dispersion relation (4). The corresponding time steps are marked by the green lines in panel c. (c) temporal variations of the radial profiles and of the amplitude | u 1 5,z | at r = 0.67.

Amplitudes and frequencies
4.5.1. Height equal diameter (Γ = 2) Figure 13 shows the temporal behavior of | u k m,z | at r = 0.67 for m = 1 to m = 7 (from top to bottom) and for k = 1 to k = 4 (blue, red, green, yellow curve). The peak values of the amplitudes are additionally listed in table 2 (Γ = 2) and in table 3 (Γ = 1.825 and Γ = 2.2) together with the frequencies obtained from the time derivative of the azimuthal phase of each mode. Note that we only list modes with a regular behavior of the frequencies, i.e., contributions which exhibit a steady and unique azimuthal drift that allows a conclusive computation of the time derivative of the azimuthal phase. The results quantitatively complement the observations made in the previous paragraph. For Γ = 2 (central column and table 2) we see a dominant forced mode with m = 1 and k = 1 and the contributions with k > 1 are negligible (smaller by a factor of 20). The higher azimuthal modes (m = 2 . . . 7) reach approximately 10% of the amplitude of the m = 1 mode with k = 1 and k = 2 slightly prevailing over k = 3 and k = 4. A striking property of the flow at Γ = 2 is the common orientation of the azimuthal drift motion of all modes (see  find modes with a retrograde azimuthal drift so no combination is possible that fulfills the requirements for a triadic resonance (in all cases δΩ = 1, see table 2). The presence of distinct frequency signals that fit to higher radial modes (which however are hardly evident in the radial profiles for the (m = 5, k = 1) mode, see figure 12) shows that there must be further contributions with chaotic behavior which are probably dominant and having no regular azimuthal drift.
We also performed simulations at a slightly larger Reynolds number Re = 10000 (at Γ = 2) and found only minor changes with respect to the run at Re = 6500 except for the modes (m = 4, k = 1) and (m = 6, k = 1) which have lower frequencies at Re = 6500. However, so far our data is not sufficient to explain the impact of Re on the drift frequencies or to establish a conclusive scaling towards more realistic Re that will be reached in the experiment. Hence, we refrain from any further discussion of the properties of the flow at larger Re.  table 3). The emergence of a triadic resonance is less explicit at Γ = 1.825. However, at this aspect ratio we find even two triadic resonances (see figure 13) with (m, k) = (6, 1) and (7, 2) and a second resonance with (m, k) = (5, 1) and (6, 2) (with much weaker amplitude, see  Table 3. Amplitudes and frequencies from the simulations for the dominant modes at Γ = 1.825 and Γ = 2.200. The columns denote (from left to right): aspect ratio, azimuthal wave number, axial wave number, maximum of the amplitude, azimuthal drift frequency from numerical data, and theoretical frequency from the dispersion relation (4) at the aspect ratio at which the linear inviscid approximation predicts a triadic resonance (always supposing a radial wave number corresponding to the first root of the dispersion relation which is confirmed in figure 11).
triadic resonance strongly oscillates with a maximum of up to 50% of the forced mode (Γ = 2.2, m = 4, k = 2) which corresponds to roughly 6% of the angular velocity of the cylindrical container.

Azimuthal shear flow
The growth of the triads goes along with a growth of an axisymmetric azimuthal shear flow that is mostly geostrophic (see left panel in figure 14). The axisymmetric mode emerges with a slight delay with respect to the free Kelvin modes which indicates a saturation process by a detuning of the resonance frequencies (see right panel in figure 14). The induced axialsymmetric flow component is negative on average and, hence, causes a breaking of the solid body rotation which increases with increasing precession ratio. In the extreme case the breaking of the induced axisymmetric flow entirely cancels the solid body rotation, giving the impression that in the laboratory frame the cylindrical container is rotating around a standing fluid (Herault et al. 2015). This effect is closely connected to the abrupt transition to a chaotic flow which takes place at a critical precession ratio (Herault et al. 2015, Kong et al. 2014.

Conclusions
We performed numerical simulations of a precession driven flow in cylindrical geometry. The simulations confirm that the energy that can be injected via precessional forcing is very sensitive to the aspect ratio of the cylindrical container. Significant contributions of higher non-axisymmetric components appear only in a limited range of aspect ratios (for Γ ∈ [1.81; 2.24]) and it turns out that in that regime the behavior of the forced mode (with m = 1) cannot be described by a linear theory. Nevertheless, in all cases the forced Kelvin mode (m = 1) dominates and free Kelvin modes with higher azimuthal wave number appear as distortions of the forced mode. The maximum response of the flow takes place at Γ max ≈ 1.871 which is remarkable far away from the resonance predicted by the linear in-viscid approximation (at Γ = 1.98982). At the resonance maximum our results yield an amplitude for the z-component of the flow of |u max z | ≈ 0.27, whichassuming isotropy -roughly agrees with the estimations of Léorat (2006). The amplitude of the forced mode decreases with increasing distance from the resonance maximum, and in the regimes with triadic resonances, we still observe values of |u max z | ≈ 0.10 (in terms of the azimuthal velocity of the container). In real units (assuming Ω c = 2π · 10Hz and R = 1 m) this would correspond to fluid velocities of 15 m/s (at resonance) and 6 m/s (off resonance) which is of the same order as the typical flow speed in the Riga dynamo experiment. Regarding the higher Kelvin modes, the saturated amplitudes of the strongest free Kelvin modes are mostly independent of the forced mode and achieve values of about 0.05 to 0.06 in terms of the angular velocity of the container corresponding to ∼ 3.5 m/s in the experiment.
At the aspect ratio envisaged for the experiment (Γ = 2 with H = 2 m and R = 1 m) we may expect an energy of the fluid flow at roughly two thirds of the maximum value at the optimum aspect ratio (Γ = 1.871). However, this optimum aspect ratio probably depends on the Reynolds number and on the forcing, so that we believe that the chosen geometry for the dynamo experiment is a good compromise to ensure an efficiently driven fluid flow in a container with a ratio of diameter to height that at least approximately reflects the geometry of planetary bodies.
Although at present we do not know whether the flow fields that emerge in our simulations will be able to drive a dynamo, we may try to estimate typical magnetic field strengths that can be expected in the planned dynamo experiment. We assume a saturation of the magnetic energy at roughly 5% of the kinetic energy of the hydrodynamic flow which is substantiated by the saturation behavior of the Riga dynamo (Gailitis et al. 2008). We further presuppose that the internal velocity caused by the precessional forcing linearly scales with the rotation of the container (and hence with Re) and assume an angular velocity of Ω c = 2π · 10 Hz, a radius R = 1 m, a height H = 2 m and a density of liquid sodium of ρ ≈ 930 kg/m 3 (at ∼ 400 K). After the change over to real units the kinetic energy obtained in our simulations at Γ = 2 (E kin ≈ 0.006) corresponds approximately to a magnetic field strength B ≈ 0.05µ 0 ρE kin Ω 2 c ≈ 40 mT. This is in the range of the Riga dynamo which is not particularly surprising since the typical velocities obtained in our simulations indeed match the velocity arising in the Riga dynamo. However, given that the corresponding m = 1 mode has produced no dynamo in the kinematic simulations it might be more realistic to refer to the free Kelvin modes with higher m which have less kinetic energy. In that case the typical magnetic field strength amounts to only B ≈ 5 mT which, however, is still in the range of the values obtained at the Von-Kármán-Sodium dynamo (Monchaux et al. 2007). We expect a stronger response and hence a larger saturation field strength for increasing Po, at least as long as we remain below the critical Poincare number for the sudden transition to a chaotic state (Po crit ≈ 0.0725 at Re = 5.65 × 10 5 , see Herault et al. 2015). However, unless we know whether or not the precession driven flow will have the right structure and sufficient magnitude to excite dynamo action at all, these estimations are only useful to demonstrate that the reference values assumed for the planned dynamo experiment will allow magnetic fields with a reasonable field strength.
Regarding the structure of the flow, we find various manifestations of nonaxisymmetric contributions in terms of free Kelvin modes with higher azimuthal wave numbers m. We find free Kelvin modes in resonance with the forced mode only when the forced mode is sufficiently far away from its primary resonance, i.e., when the amplitude of the forced mode is not too strong. Triadic resonances can be clearly identified using a combined Fourier-Discrete Sine transformation in the azimuth and along the axis with structure and frequencies close to predictions from linear theory. Actually, we were also expecting triads around Γ ≈ 2 with m = 5 and m = 6. Indeed, we do find patterns of free Kelvin modes in that regime, but no combination of them satisfies the triadic resonance conditions. Instead, we see free Kelvin modes solely with retrograde drift motion and for each azimuthal wavenumber m the contributions with k = 1 and k = 2 are approximately of equal strength so that in sum we see a clear breaking of the equatorial symmetry of the flow. This symmetry breaking has been essential for the functioning of a dynamo in simulations of Tilgner (2005) in a precessing sphere whereas it seemed less important in the direct numerical simulations of Nore et al. (2011). However, preliminary kinematic dynamo simulations we conducted recently with an analytic flow of Kelvin modes suggest that the joint appearance of contributions with k = 1 and k = 2 considerably facilitates the occurrence of dynamo action.
The triadic instability requires a rather long time to emerge and saturates to a final state with slow periodic growth and decay. This period is likely related to the precession time scale but a conclusive correlation requires further simulations at different Po. The periodic growth followed by a fast decay of the free Kelvin modes in the regime with triadic resonances reminds of the so called resonant collapse (McEwan 1970, Manasseh 1992. However, in our simulations the energy is essentially distributed among very few, large scale modes. These modes are the fundamental forced mode and two free Kelvin modes that establish the triad and an axisymmetric, essentially geostrophic mode. We do not spot any turbulence-like behavior during the decay of the free Kelvin modes, i.e. there is no transition from the dominance of the large scale mode into small scale turbulent flow which would be characteristic for the collapse phenomenon. This behavior might change when the Reynolds number approaches larger (i.e. more realistic) values. This is substantiated by the simulations at Γ = 2. Here we see no triadic resonance, but rather a quasi-periodic occurrence of free Kelvin modes that exhibit a trend towards larger radial wave numbers during their decay which may represent (part of) a cascade to small scale structures.
In principle the occurrence of very few dominant modes calls for a low-dimensional model. However, the ansatz (9) used in this study to explain the occurrence of triadic resonances is too simplistic to allow for a reasonable description of the non-linear behavior and the saturation obtained in the simulations presented above. A more sophisticated approach is presented by Lagrange et al. (2011) who use a low-dimensional model with four coupled ordinary differential equations that describe the amplitude of the forced mode, the free Kelvin modes and an axisymmetric geostrophic mode. This model allows for saturation via the axisymmetric mode and for a reinforcement of the forced mode by the free Kelvin modes and qualitatively reproduces the results obtained in our simulations. However, in order to yield a quantitative agreement, advanced models are still required.