1. Introduction
Large-scale bioreactors used in industry suffer from severe limitations due to slow cell growth (Garcia-Ochoa & Gomez Reference Garcia-Ochoa and Gomez2009). This problem is due to a lack of good blending, transfer and mixing of the gases and nutrients necessary for cell life. This is why, generally, an additional convective stirring is performed by introducing bubbles or using rotating blades. These methods partially solve the problem but also create a strong shear that inhibits the growth of many cells and can even be lethal for fragile cells (Cherry & Papoutsakis Reference Cherry and Papoutsakis1986; Doran Reference Doran1999). Recently, a soft mixer consisting of a tilted rotating cylinder (Meunier Reference Meunier2020) was proposed as a suitable bioreactor for the culture of microalgae. This set-up creates a flow similar to precessing flows, which are known to create an efficient mixing (Goto, Shimizu & Kawahara Reference Goto, Shimizu and Kawahara2014). But the soft mixer is simpler since the forcing of the flow is due to the free surface, which remains flat by gravity. The free surface can also create an efficient mixing in a rotating sphere (Watanabe & Goto Reference Watanabe and Goto2022). However, the cylindrical and spherical geometries limit the amount of light (necessary for photosynthesis) that can be brought inside the photobioreactor. The use of a cylindrical annulus instead of a cylinder is thus an excellent solution to increase the illumination which can enter through both the inner and the outer cylinders.
From a physical point of view, the soft mixer is a cylinder rotating around its own axis, tilted with respect to gravity. Thompson (Reference Thompson1970) showed that the free surface generates a flow that can be resonant for specific heights of fluids. In fact, this excitation is a particular case of the forcing by a tilted cover rotating at a different velocity than the rotating cylinder, as studied by McEwan (Reference McEwan1970). This type of forcing is very similar to the case of precessing flows, which have been largely studied due to their application to the precession of the Earth (Le Bars, Cébron & Le Gal Reference Le Bars, Cébron and Le Gal2015) and for their abilities to drive a dynamo (Giesecke et al. Reference Giesecke, Vogt, Gundrum and Stefani2018).
In a precessing cylinder, the flow becomes highly energetic and turbulent due to the resonance of global inertial modes, also called Kelvin modes (Kelvin Reference Kelvin1880), which are the equivalent of inertial waves in the cylindrical geometry. These modes are resonant when the height of the cylinder is equal to an odd multiple of the half-wavelength of the inertial mode (Manasseh Reference Manasseh1992, Reference Manasseh1994). The inertial modes are saturated at low Ekman number by surface Ekman pumping (Gans Reference Gans1970; Meunier et al. Reference Meunier, Eloy, Lagrange and Nadal2008) leading to an amplitude scaling as the square root of the Ekman number. The nonlinear saturation is harder to describe analytically because the inviscid nonlinear interaction of the inertial mode with itself does not generate any geostrophic mode (i.e. axisymmetric and axially invariant), as stated by Greenspan (Reference Greenspan1969). However, the viscous and nonlinear interaction of the inertial mode with itself leads to a spin-down of the fluid which detunes the resonance and creates a nonlinear saturation (Albrecht et al. Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021; Gao et al. Reference Gao, Meunier, Le Dizès and Eloy2021).
The nonlinear interactions can also trigger triadic instabilities (Kerswell Reference Kerswell1999), in which the forced inertial mode can couple with two free inertial modes that grow simultaneously (Lagrange et al. Reference Lagrange, Eloy, Nadal and Meunier2008). The viscous growth rate can be predicted analytically, in excellent agreement with the experimental results (Lagrange et al. Reference Lagrange, Meunier, Nadal and Eloy2011) and the numerical results (Albrecht et al. Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2015; Giesecke et al. Reference Giesecke, Albrecht, Gundrum, Herault and Stefani2015; Albrecht et al. Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2018; Lopez & Marques Reference Lopez and Marques2018). The transition to a fully turbulent flow depends on the parameters of the precession and is strongly related to the mean axisymmetric flow (Pizzi et al. Reference Pizzi, Giesecke, Šimkanin and Stefani2021). This transition may even by subcritical with a relaminarization after a probability distribution of lifetime as in pipe flows (Heraulta et al. Reference Heraulta, Gundrum, Giesecke and Stefani2015).
Despite the numerous works on precessing cylinders, there are very few studies on precessing annuli, although Lin, Noir & Jackson (Reference Lin, Noir and Jackson2014) showed that the global inertial modes and the triadic resonance instabilities are very similar. Recently, Xu & Harlander (Reference Xu and Harlander2020) studied experimentally the case of a rotating annulus tilted with respect to gravity. As in the cylinder, they found a resonance of the forced Kelvin mode for specific heights of fluid. The flow then destabilizes at low Ekman number thanks to a triadic resonance at small tilt angles and thanks to a shear instability at larger tilt angles.
The goal of this paper is to characterize this flow and its mixing properties. The experimental apparatus to study the flow is described in § 2. Experimental and theoretical characterization of the forced flow and of the instabilities are done in §§ 3 and 4, respectively. Finally we address preliminary results on mixing (§ 5) inside this new soft mixer.
2. Materials and methods
2.1. Experimental set-up
We are interested in studying a mixer formed by an annular cavity (see figure 1), with outer radius 
$R=9$ cm and inner radius equal to 
$R_i=3$ and 
$5$ cm. The peripheral walls and the bottom are made in polymethylmethacrylate (PMMA) for visualization. The annulus rotates around its own axis, defined as the 
$z$-axis, with angular velocity 
$\varOmega$ measured with an accuracy of 
$0.1$ r.p.m. thanks to a direct current (DC) brushless motor with a reduction of 3.7. The axis of rotation is tilted at an angle 
$\alpha$ with respect to gravity 
$\boldsymbol {g}$, 
$\alpha$ varying from 
$0.5^{\circ }$ to 
$10^{\circ }$ with an accuracy of 
$0.1^{\circ }$. In addition, the annulus is filled with water (with kinematic viscosity 
$\nu$, known within 1 % thanks to temperature measurement) at a height 
$H$, which is measured before the cylinder is tilted. A blob of dye, characterized by the molecular diffusivity in water 
$\kappa$, is introduced to study the mixing properties of the flow.
Experimental set-up of a tilted rotating annulus.
2.2. Non-dimensional parameters
The system is characterized by six non-dimensional parameters.
(i) The aspect ratio
$r_i=R_i/R$ is the parameter characterizing the flow geometry. Two different values ${r_i = 0.33; 0.56}$ are used in the experiments to observe how the flow depends on the confinement.(ii) The non-dimensional height
$h=H/R$ is the parameter that can be tuned to select the different resonances of the flow (see § 3). In the experiment, we consider values $0.5\le h \le 2$ to excite four different resonances.(iii) The tilting angle
$\alpha$ is the parameter describing the intensity of the forcing applied by the free surface. In the experiment, we consider values $1^\circ \le \alpha \le 10^\circ$.(iv) The Ekman number
(2.1)quantifies the ratio of the viscous term over the Coriolis force. In the experiment, it is varied from\begin{equation} E=\frac{\nu}{\varOmega R^2}\end{equation}$E=10^{-5}$ to $E=10^{-3}$.(v) The Froude number
(2.2)is defined as the square of the ratio of inertial velocity to surface gravity wave speed. This parameter is chosen to be small (\begin{equation} Fr=\frac{\varOmega^2R}{g}\end{equation}$Fr\lesssim 10^{-1}$) to neglect the curvature of the free surface in the theory.(vi) The Schmidt number
(2.3)is defined as the ratio of the viscous diffusion rate and molecular diffusion rate. In the experiment, where we used rhodamine B, we have a value of\begin{equation} Sc= \frac{\nu}{\kappa} \end{equation}$Sc=2.3\times 10^{3}$ such that the dye streaks are much thinner than the coherent structures of the flow.
2.3. Measurement techniques
The flow was characterized using three different methods, particle visualization (Gauthier, Gondret & Rabaud Reference Gauthier, Gondret and Rabaud1998), laser induced fluorescence and particle image velocimetry (PIV).
Particle visualizations are made with a vertical laser sheet and introducing flat mica particles (average diameter of 
$50~\mathrm {\mu }$m) in the fluid to qualitatively underline the coherent structures of the flow. This allowed us to observe the presence and the onset time of instabilities.
The PIV measurements are obtained using a double pulsed yttrium-aluminium-garnet laser (YAG) of 
$150$ mJ per pulse at 
$532$ nm with a frequency of 
$10$ Hz and with a time between pulses varying from 
$0.1$ ms to 
$100$ ms. The laser sheet is normal to the axis of the annulus and has a thickness close to 
$5$ mm. It is reflected with two vertical mirrors in order to prevent the shade made by the inner cylinder. The images are acquired by a 
$4$ megapixel camera, positioned below the annulus to decrease the distortion caused by the water, synchronized with the pulsing of the laser. Since the camera is fixed, it does not rotate with the experiment; the images of the second laser pulse are rotated, to remove the solid body rotation and measure the velocity in the annulus frame. The time interval could thus be enlarged in order to increase the accuracy of the velocity field measurements. In the end, vortices of the order of 
$1\,\%$ of the solid body rotation can be detected. The images are processed by a cross-correlating algorithm optimized for large velocity gradients (Meunier & Leweke Reference Meunier and Leweke2003).
Furthermore, using the temporal Fourier transform of the vorticity field filtered by the acquisition frequency, it is possible to measure the frequencies of the observed motions. Finally, the amplitudes of motion are measured by taking the scalar product of the experimental velocity field with the theoretical modes (see (B2)).
Laser induced fluorescence was made by inserting a blob of rhodamine B at the free surface and illuminating with a continuous 
$532$ nm and 
$100$ mW diode to create a laser sheet. The laser sheet is made parallel to the bottom thanks to a large spherical lens with a diameter of 
$20$ cm, in order to have a uniform intensity in the annulus. The intensity 
$I$ recorded by the camera is calibrated for each experiment as a function of the concentration by placing inside the annulus test tubes of rhodamine at a known concentration. The data are fitted, with an accuracy of 
$10\,\%$, by the empirical law (see figure 2)
\begin{equation} \frac{c(I)}{c_0}= \varGamma \frac{I-{\rm I}_0}{1-I/{\rm I}_{{sat}}}, \end{equation}
where 
${\rm I}_0$ and 
${\rm I}_{{sat}}$ are the intensity of background and saturation of dye, respectively, with the empirical parameter 
$\varGamma$.
Calibration between the intensity 
$I$ measured by the camera and the normalized concentration of rhodamine 
$c/c_0$.
3. Forced mode
Using 
$R$ and 
$\varOmega ^{-1}$ as length and time scales, we define the non-dimensional variables 
$r=r^*/R$, 
$t=\varOmega t^*$ and 
$z=z^*/R$ from the dimensional variables 
$(r^*,z^*,t^*)$. We use cylindrical coordinates for the position 
$(r,\theta,z)$ and for the velocity 
${\boldsymbol {u}}=(u,v,w)$.
For small Froude number (defined in (2.2)), the free surface is horizontal and motionless in the laboratory frame. In the reference frame rotating with the annulus, it is thus located at a height
\begin{equation} \zeta(r,\theta)=h-\alpha r \cos(\theta+t). \end{equation}
This position oscillates at a non-dimensional frequency 
$\varpi =1$ because the free surface rotates at angular velocity 
$-\varOmega$ (in dimensional units) in the rotating frame. The normal velocity 
$w={\rm d}\zeta /{\rm d}t$ of the free surface generates a forcing excitation with azimuthal wavenumber 
$m=1$:
\begin{equation} w(z=\zeta)=\alpha r \sin(\theta+ t). \end{equation}The flow is made by a constructive interference of inertial waves which lead to global Kelvin modes (Kelvin Reference Kelvin1880), as shown by Thompson (Reference Thompson1970) for a cylinder and by Xu & Harlander (Reference Xu and Harlander2020) for an annulus. The annular modes can be derived by considering the linearized Navier–Stokes (N–S) equations in the frame of the annulus
\begin{equation}
\left.\begin{array}{@{}c@{}}
\partial_t\boldsymbol{u}+2\hat{\boldsymbol{z}}\wedge\boldsymbol{u}={-}\boldsymbol{\nabla}p+E\nabla^{2}\boldsymbol{u},\\
\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0,
\end{array}\right\}\end{equation}
where we have introduced the reduced pressure 
$p=P-\rho gz+\rho r^2\varOmega ^2/2$.
This equation is solved in the limit of 
$E\to 0$ by deriving the equation for pressure. The easiest way to derive it is to take the divergence and curl of the equation, with the incompressibility condition and assuming that 
$(\boldsymbol {u},p)\to (\boldsymbol {u},p){\rm e}^{{\rm i}\varpi t}$, which combined give the Poincaré equation (Greenspan Reference Greenspan1968)
\begin{equation} \nabla^2p(r,\theta,z)-\frac{4}{\varpi^2}\partial_z^2p(r,\theta,z)=0.\end{equation} By decomposing the pressure in Fourier mode in 
$z$ and 
$\theta$ is it possible to show that the previous equation becomes the Bessel equation
\begin{equation} \frac{{\rm d}^2}{{\rm d}r^2}p(r)+\frac{1}{r}\frac{{\rm d}}{{\rm d}r}p(r)+ \left[k^2\left(\frac{4}{\varpi^2}-1\right)-\frac{m^2}{r^2}\right]p(r)=0,\end{equation}
where 
$m$ is the azimuthal wavenumber and 
$k$ is the axial wavenumber.
The general solution with 
$m=1$ and frequency 
$\varpi =1$ is given by (see Lin et al. Reference Lin, Noir and Jackson2014)
\begin{equation} \boldsymbol{u}=\sum_{j=1}^{\infty}{\rm Re}\left[A_j\begin{pmatrix} u_j(r)\cos(k_jz) \\ v_j(r)\cos(k_jz) \\ w_j(r)\sin(k_jz) \end{pmatrix}{\rm e}^{{\rm i}(\theta+ t)}\right] , \quad p=\sum_{j=1}^{\infty}{\rm Re} \left[A_j p_j(r)\cos(k_jz){\rm e}^{{\rm i}(\theta+t)}\right], \end{equation}with
\begin{gather}
\left.\begin{gathered}
p_j(r)=c_j^J{\rm J}_1(\sqrt{3}k_jr)+c_j^Y{\rm Y}_1(\sqrt{3}k_jr),\\
\displaystyle u_j(r)={-}\dfrac{{\rm i}}{3}\left(\partial_rp_j+\dfrac{2}{r}p_j\right) \\
\displaystyle ={-}\dfrac{{\rm i}}{3r}\left[\!c^{J}_j\left({\rm J}_{1}(\sqrt{3}k_jr)\,{+}\,\sqrt{3}k_j{\rm
J}_{0}(\sqrt{3}k_jr)\!\right)\,{+}\,c^{Y}_j\left(\!{\rm Y}_{1}(\sqrt{3}k_jr)\,{+}\,\sqrt{3}k_j{\rm Y}_{0}(\sqrt{3}k_jr)\!\right)\!\right],\\
\displaystyle v_j(r)=\dfrac{1}{3}\left(2\partial_rp_j+\dfrac{1}{r}p_j\right)\\
\displaystyle ={-}\dfrac{1}{3r}\left[c^{J}_j\left({\rm J}_{1}(\sqrt{3}k_jr)+2\sqrt{3}k_jr{\rm J}_{0}(\sqrt{3}k_jr)\right)\right.\\
\displaystyle \left.+\,c^{Y}_j\left({\rm Y}_{1}(\sqrt{3}k_jr)+2\sqrt{3}k_jr{\rm Y}_{0}(\sqrt{3}k_jr)\right)\right],\\
\displaystyle w_j(r)={-}{\rm i}k_jp_j={-}{\rm i}k_j\left(c_j^J{\rm J}_1(\sqrt{3}k_jr)+c_j^Y{\rm Y}_1(\sqrt{3}k_jr)\right),\end{gathered}\right\}\end{gather}where J and Y are Bessel functions of the first and second kind, and Re means the real part.
The ratio 
$c_j^J/c_j^Y$ of coefficients of the real solution is found from the boundary condition at 
$r=1$:
\begin{equation} u(r=1)=0,\quad \forall\, \theta,z. \end{equation}We choose (although this is valid within a multiplicative constant)
\begin{equation} c_j^J={\rm Y}_1(\sqrt{3}k_j) + k_j\sqrt{3}{\rm Y}_0(\sqrt{3}k_j) , \quad c_j^Y={-}\left[{\rm J}_1(\sqrt{3}k_j) + k_j\sqrt{3}{\rm J}_0(\sqrt{3}k_j)\right].\end{equation} In general, the slip boundary condition at 
$r=r_i$,
\begin{equation} u(r=r_i)=0,\quad \forall\, \theta,z, \end{equation}
gives the dispersion relation between the axial wavenumber 
$k$ and the frequency of the modes 
$\varpi$. In our case, the forcing frequency is equal to 
$\varpi =1$, which fixes discrete values of the axial wavenumber 
$k_j$ thanks to the dispersion relation
\begin{equation} c_j^J\left[{\rm J}_1(\sqrt{3}k_jr_i) + \frac{k_j\sqrt{3}}{r_i}{\rm J}_0(\sqrt{3}k_jr_i)\right]+ c_j^Y\left[{\rm Y}_1(\sqrt{3}k_jr_i) + \frac{k_j\sqrt{3}}{r_i}{\rm Y}_0(\sqrt{3}k_jr_i)\right]=0. \end{equation}
Numerical values for 
$r_i=0.333$ and 
$r_i=0.556$ are given in table 1.
Numerical values of the first four resonance wavenumbers and height of the 
$j$th Kelvin mode given by the dispersion relation in (3.11), for 
$r_i=0.333$ (a) and 
$r_i=0.556$ (b).
The structure of the first Kelvin mode is shown in figure 3. It exhibits a global flow towards negative 
$x$ in the upper part of the cylinder (see figure 3a) and a global flow towards positive 
$x$ in the lower part of the cylinder (see figure 3a). This is due to the choice of the height 
$h=1.306$ corresponding to half a wavelength of this Kelvin mode. This creates a meridional recirculation observed in figure 3(b).
Structure of the first Kelvin mode (3.6a,b) with azimuthal wavenumber 
$m=1$. (b) Visualization of the velocity field in two sections parallel to the 
$z$-axis showing the meridional recirculation. Two sections normal to the 
$z$-axis are shown for 
$z=0.13\equiv z_{{laser}}$ in (a) and for 
$z=1.18\equiv h_1-z_{{laser}}$ in (c). The amplitude 
$A_1$ is given by the viscous theory (3.17) for 
$r_i=0.333$, 
$E=4.1 \times 10^{-4}$, 
$h=1.31$ and 
$\alpha =3 ^\circ$.
The vorticity of the Kelvin mode can be calculated using the Kelvin modes property 
$\boldsymbol {\omega }=-2{\rm i}\partial _z\boldsymbol {u}$,
\begin{equation} \omega_{z}=\sum_{j=1}^{\infty}{\rm Re}\left[{-}2k_j^2A_j p_j(r)\cos(k_jz){\rm e}^{{\rm i}(\theta+t)}\right]. \end{equation} Near the free surface, the global flow creates an outer lobe of negative vorticity for 
$y$ positive and an outer lobe of positive vorticity for 
$y$ negative. The axial antisymmetry is simply due to the azimuthal wavenumber 
$m=1$. Additional small lobes of vorticity with opposite sign are located close to the inner cylinder.
The amplitude of the modes must be set in order to fulfil the boundary condition at the free surface (3.2). This boundary condition is also valid at 
$z=h$ rather than at 
$z=\zeta$ at first order in 
$\alpha$. Multiplying by 
$r\bar {\boldsymbol {u}}$ and integrating on 
$r$ gives (to find the amplitude we use the orthogonality rules 
$\int _{r_i}^{1}w_i(r)w_j(r)r\,{\rm d}r\propto \delta _{ij}$)
\begin{equation} A_j=\frac{\alpha F}{T\sin(k_j h)},\end{equation}with
\begin{equation} F=\int_{r_i}^{1}r^2\left[c_j^J{\rm J}_1(\sqrt{3}k_jr)+c_j^Y{\rm Y}_1(\sqrt{3}k_jr)\right]{\rm d}r \end{equation}and
\begin{equation} T=k_j\int_{r_i}^{1}r\left[c_j^J{\rm J}_1(\sqrt{3}k_jr)+c_j^Y{\rm Y}_1(\sqrt{3}k_jr)\right]^2\,{\rm d}r. \end{equation}
It is interesting to note that the term 
$F$ can be integrated analytically as given in (A16). Numerical values of 
$T$ and 
$F$ are given in table 2. This theory is similar to the one given in Xu & Harlander (Reference Xu and Harlander2020) except that we give analytical predictions of the amplitudes whereas they calculated the amplitudes numerically by fitting the theoretical solution with five modes to the top boundary condition (3.2).
Numerical values of the viscous correction terms, for 
$r_i=0.33$ and 
$r_i=0.56$ as a function of 
$h$.
From this expression it is clear that the amplitude diverges when 
$\sin (k_j h)=0$, i.e. for 
$k_j=n{\rm \pi} /h$. It means that the 
$j$th Kelvin mode is resonant when the height of fluid 
$h$ is an integer multiple of the half-wavelength (see table 1), i.e. when 
$h=n h_j$ with
\begin{equation} h_j={\rm \pi}/k_j. \end{equation}
In that case, the axial velocity of the mode vanishes at 
$z=h$ and is thus unable to contribute to fulfil the boundary condition (3.2) unless its amplitude is infinite. This is similar to the resonance of a string excited by an oscillation of its end when the end is located at the node of a mode.
This inviscid prediction is plotted in figure 4(a) as a green dashed line. It clearly exhibits a divergence at 
$h=h_1$. Outside the resonance (for 
$h<0.8$), the theoretical value is close (within 20 %) but slightly larger than the experimental value. Here, the experimental value is measured by projecting the total velocity field at 
$z_{{laser}}=0.13$ onto the inviscid Kelvin mode (see Meunier et al. (Reference Meunier, Eloy, Lagrange and Nadal2008) for further details).
Amplitude of the first mode of the first resonance 
$A_1/\alpha$ for 
$r_i=0.333$. It is plotted as a function of 
$h$ for 
$E=4.1\times 10^{-4}$ (a) and as a function of 
$E$ for 
$h = h_1$ (b). Experimental measurements are obtained at 
$\alpha = 3^{\circ }$ (black symbols) and at 
$\alpha = 1^{\circ }$ (blue symbols). Experimental data are represented by 
$\circ$ for stable experiments and 
$+$ for unstable experiments. The dashed green line corresponds to the inviscid theory (3.13) and the black solid line corresponds to the viscous theory (3.17). In (b) the dashed orange line is the viscous amplitude neglecting the volumic term (i.e. with 
$V=0$).
Obviously, the experimental amplitude does not diverge at the resonance because of viscous and nonlinear effects. As a first step, we have calculated the saturation due to the Ekman pumping and volumetric damping. It can be accounted for analytically by doing an analysis similar to the case of a cylinder (Gans Reference Gans1970; Meunier et al. Reference Meunier, Eloy, Lagrange and Nadal2008; Meunier Reference Meunier2020). As shown in Appendix A, the Ekman pumping 
$\tilde {\boldsymbol {u}}$ has a simple expression as a function of the 
$z$-derivative and the normal derivative of the normal velocity (A1). The amplitude is then found by a solvability condition including the Ekman pumping at the bottom, the inner and the outer cylinders (see Appendix A). It leads to an analytic expression for the amplitude
\begin{equation} A_j=\frac{\alpha F}{T\sin(k_j h)+\sqrt{E}\left(B+C_e+C_i\right)+E V},\end{equation}
where 
$F$ is the forcing term, 
$(B+C_e+C_i)$ describes the surface Ekman pumping term (bottom, external and internal cylinder) and 
$V$ is the volumetric term (see table 2 for reference values).
This theoretical amplitude is shown in figure 4(a) as a black solid line. It is in excellent agreement with the experiment around the resonance, with no fitting parameter. It should be noted that viscous corrections lead to a small shift of the maximum amplitude to 
${h=1.1h_{1}}$ due to the imaginary part of 
$B+C_e+C_i$.
When varying the Ekman number (figure 4b), the scaling of the amplitudes change from 
$A_1\propto E^{-1/2}$, at small 
$E$ where the surface term is dominant, to 
$A_1\propto E^{-1}$ when the volumetric term becomes dominant. There is again a good agreement between the experimental and numerical values.
Finally, we would like to note the fact that the free surface in the experiment cannot be considered completely clean and that the presence of contaminants changes the boundary condition on the free surface. However, it should be noted that in the extreme case of a no-slip boundary (where a liquid surface covered by a rigid layer of non-mobile surfactants) the amplitude of the flow is altered by a maximum of 
$15\,\%$ in absolute value. We have, therefore, consciously not accounted for surface contamination in the present analysis of the flow, but do acknowledge that it may need to be in other particular cases.
Figure 5 shows a comparison of the velocity and vorticity fields obtained experimentally and numerically (using a finite element code). They exhibit the same structures and amplitudes, with a mean flow towards positive 
$x$. This is in fair agreement with the theoretical prediction of figure 3(c) despite a larger vorticity by a factor of two. The main difference with theory is due to the presence of the Ekman layer.
Axial vorticity 
$\omega _z$ and velocity field at 
$z=0.13$ for 
$\alpha =3^\circ$, 
$E=4.1\times 10^{-4}$ and 
$r_i=0.33$. The experimental PIV measurement (a) is compared with a direct numerical simulation (DNS) obtained with the finite element code Comsol (b).
To conclude, we have given an analytic solution for the forced mode with no fitting parameter, in good agreement with the observation in experimental and numerical results.
4. Instabilities of the flow
As in a tilted cylinder and in precessing containers, the flow becomes unstable when 
$E$ decreases. This can be addressed experimentally by using PIV and visualizations.
4.1. Onset time of the instability
The experimental protocol consisted in getting the solid body rotation without tilt, and then tilting the annulus smoothly within a few seconds. Particle visualizations exhibited periodic coherent structures appearing in the laser sheet after a time depending on the height of fluid and on the Ekman number. This onset time 
$t_{{onset}}$ is plotted as a function of 
$h$ in figure 6 for the two aspect ratios 
$r_i=0.33$ and 
$r_i=0.56$. The Ekman number 
$E$ was chosen small enough for an instability to appear.
Inverse of the onset time 
$t_{onset}$ as a function of the height 
$h$ for 
$E=3.6\times 10^{-5}$, 
$\alpha =3^\circ$ and 
$r_i=0.56$ (a) and 
$r_i=0.33$ (b). Red dashed lines represent the theoretical values of the resonances, while experimental data are indicated by 
$\boldsymbol {\circ }$ for a triadic instability and 
$\boldsymbol {\diamond }$ for a centrifugal instability.
Because the instabilities are triggered from the amplitudes 
$A_j$ of the forced mode, the instability occurs faster when the height of water is near the resonant heights 
$h_j$. Indeed, experiments at 
$r_i=0.56$ show that 
$t_{{onset}}$ reaches its minimum value at the first and second resonances of the first Kelvin mode. In between these two resonances, the onset time is 10 times larger. For 
$r_i=0.33$, the flow is strongly unstable, but the onset time saturates at 
$1/0.017$ which correspond to the time needed for the forced mode to grow.
In the following we will focus only on the resonant heights.
4.2. Triadic instability
Figure 7 shows the experimental vorticity fields of the instability at midheight and quarter height for 
$r_i=0.33$ and 
$h=h_{1}$. The vorticity fields have been extracted from the mean forced mode by calculating the temporal Fourier component at the frequency 
$-9.58$ in the laboratory frame of reference. The instability exhibits a ring of alternate vortices with an azimuthal wavenumber 
$m'=9$ at quarter height and an azimuthal wavenumber 
$m''=10$ at midheight. Although these PIV measurements have been done separately, the difference in azimuthal wavenumber is not due to randomness since it always gives the same result. Xu & Harlander (Reference Xu and Harlander2020) have actually shown that the two modes are present at the same time. However, the Kelvin modes have an axial vorticity proportional to 
$\cos (kz)$ (such that 
$w \sim \sin (kz)$ vanishes at 
$z=0$). If the first mode has an axial wavenumber 
$k'={\rm \pi} /h_1$, it will not be visible at midheight. Similarly, if the second mode has an axial wavenumber 
$k''=2{\rm \pi} /h_1$ it will not be visible at quarter height. This explains the simultaneous presence of the two modes and their separate measurements in two well-chosen sections.
Vorticity field at the resonance, filtered at the frequency 
$-9.58$ in the laboratory frame of reference and averaged over one rotation. For 
$E=10^{-4}$, 
$r_i=0.33$, 
$\alpha =1^\circ$ and two different positions of the laser sheet: 
$h_{{laser}}=0.25$ (a) and 
$h_{{laser}}=0.5$ (b). The position of the laser sheet is chosen in order to separate the two free modes (respectively, at the nodes of 
$m'$ and 
$m''$ mode) thus revealing the azimuthal wavenumbers 
$m'=9$ (a) and 
$m''=10$ (b).
Theoretically, the first Kelvin modes with 
$(m',k')=(9,{\rm \pi} /h_1)$ and 
$(m'',k'')=(10,2{\rm \pi} /h_1)$ have frequencies 
$\varpi '=-0.348$ and 
$\varpi ''=0.649$ which give almost the same frequency (equal to 
$-9.348$ and 
$-9.351$) in the laboratory frame of reference (obtained by adding 
$-m'\varOmega$ or 
$-m''\varOmega$ with 
$\varOmega = 1$ from dimensionalization), very close to the experimental frequency 
$-9.58$. To conclude, the modes seem to respect the following conditions:
\begin{equation} k''-k'=k_1, \end{equation}
where 
$k_1={\rm \pi} /h_{1}$ is the axial wavenumber of the forced resonant mode;
\begin{equation} m''-m'=m, \end{equation}
where 
$m=1$ is the azimuthal wavenumber of the forced resonant mode;
\begin{equation} \varpi''-\varpi'=\varpi, \end{equation}
where 
$\varpi =1$ is the frequency of the forced resonant mode. This is fulfilled theoretically only for modes 
$m'=9$ and 
$m''=10$. These conditions are exactly similar to the case of precessing cylinders (Lagrange et al. Reference Lagrange, Meunier, Nadal and Eloy2011; Albrecht et al. Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2015) and of tilted cylinders (Meunier Reference Meunier2020) where modes of azimuthal wavenumbers 
$m'=5$ and 
$m''=6$ were observed. The higher azimuthal wavenumbers in the annulus probably come from the radial confinement. Moreover, it should be noted that theoretically these conclusions can also be drawn with arguments of symmetry of the problem (Marques & Lopez Reference Marques and Lopez2015).
These conditions indicate that a triadic resonance between the three modes is possible: the forced mode may feed simultaneously the two free modes by nonlinear interactions, leading to a triadic resonance instability. Although the theory has been done for a precessing cylinder (Lagrange et al. Reference Lagrange, Meunier, Nadal and Eloy2011) and for a tilted cylinder (Meunier Reference Meunier2020), it has never been done for an annulus despite experimental evidences for a precessing annulus (Lin et al. Reference Lin, Noir and Jackson2014) and a tilted annulus (Xu & Harlander Reference Xu and Harlander2020). We derive here the linear stability analysis for the annulus. Assuming that the amplitudes 
$A'$ and 
$A''$ of the two modes have a temporal variation of the first order (proportional to 
$A_1 \propto \alpha /\sqrt {E}$) and considering only the first-order terms in these amplitudes, it is possible to derive the amplitude equations (see Appendix B)
\begin{align}
\left.\begin{array}{@{}c@{}}
\partial_tA'=N'\bar{A}_1A''-D' A',\\
\partial_tA''=N'' A_1 A'-D'' A'',\end{array}\right\} \end{align}
where 
$N'$ and 
$N''$, defined in Appendix B, are analytic coefficients describing the nonlinear interaction of the forced mode and the free ones, and 
$\bar {A}_1$ is the complex conjugate of the amplitude 
$A_1$ defined in (3.17). They are given by integrals in Appendix B which can be calculated analytically in the case of a small gap using an asymptotic expansion (see Appendix C).
The parameters
\begin{equation} D'=\sqrt{E}s'+E v'+{\rm i}\Delta k'q',\quad D''=\sqrt{E}s''+E v''+{\rm i}\Delta k''q'',\end{equation}
are also analytical (see Appendix B for definition of the coefficients and numerical values). The first two terms describe the viscous corrections, with 
$s$ for the Ekman damping on the edges and 
$v$ for the volumetric damping. The last term is the detuning effect, due to the fact that the free modes are not exactly resonant with the forced mode.
The set of equations is exactly similar to (4.13a,b) in Lagrange et al. (Reference Lagrange, Meunier, Nadal and Eloy2011) and can be easily integrated with an exponential solution of the type 
$A',A''\propto {\rm e}^{\sigma t}$ leading to a growth rate 
$\sigma$ satisfying
\begin{equation} \left(\sigma+D'\right)\left(\sigma+D''\right)=\left|A_1\right|^2N'N''. \end{equation}
The flow is stable if 
${\rm Re}(\sigma )<0$ which leads to the threshold tilt angle (in degrees) as
\begin{equation} \alpha(E)=\frac{180}{\rm \pi}\frac{1}{\left|A_1/\alpha\right|\sqrt{N'N''}} \sqrt{ D'_{{{\rm Re}}{}}D''_{{{\rm Re}}{}}\left(1+\frac{\left(D'_{{{\rm Im}}{}}-D''_{{{\rm Im}}{}} \right)^2}{\left(D'_{{{\rm Re}}{}}+D''_{{{\rm Re}}{}}\right)^2}\right)},\end{equation}
where 
$D_{{{\rm Re}}{}}$ and 
$D_{{{\rm Im}}{}}$ are real and imaginary parts of the viscous correction.
This stability curve is plotted in figures 8–9 for the two different aspect ratios. It is in excellent agreement with the experimental data for 
$r_i=0.333$ with an accuracy of approximately 10 %. It should be noted that there are no fitting parameters in this viscous theory since all coefficients can be computed analytically. The agreement is not as good for 
$r_i=0.56$ but still correct with an error smaller than 30 %. The effect of the detuning is very small for 
$r_i=0.56$ and invisible for 
$r_i=0.33$.
Stability diagram in the parameter space 
$\alpha$ and 
$E$, for 
$r_i = 0.33$ and 
$h = h_{1}$. Colour symbols represent unstable experiments (red for triadic instability and green for centrifugal instability) and black symbols represent stable experiments. The experiments are visualized using PIV 
$_{\boldsymbol{\bigtriangleup}}$ and particle visualizations 
$\boldsymbol {\circ }$. The black solid line corresponds to the theoretical threshold of (4.7) with coefficients given in table 3 for the most unstable modes 
$(m',m'')=(9,10)$.
Stability diagram in the parameter space 
$\alpha$ and 
$E$, for 
$r_i = 0.55$ and 
$h = h_{1}$ (a) and 
$h = 2h_{1}$ (b). Red symbols represent unstable experiments and black symbols represent stable experiments. The theoretical triadic resonance critical curves of (4.7) with the coefficients given in table 3 are shown in black with detuning and in blue without detuning for the most unstable modes 
$(m',m'')=(17,18)$.
The scaling of the threshold curves can be explained asymptotically. First, we consider only the surface damping of the forced Kelvin mode (which is valid for 
$E\lesssim 10^{-3}$, see figure 4b). It leads to a scaling of the forcing amplitude as 
$A_1\sim F/\sqrt {E}S$ (where 
$S=(B+C_e+C_i)$). Then, (4.7) has three different scalings depending whether the detuning, the surface damping of the free modes or the volumic damping of the free modes is dominant (see also Lagrange et al. Reference Lagrange, Meunier, Nadal and Eloy2011)
\begin{equation}
\alpha(E)\sim\begin{cases}
E^{1/2}\left|\dfrac{S}{F}\right|\sqrt{\dfrac{|\Delta
k'q'\Delta k''q''|}{N'N''}}, & \textrm{for}\ E\lesssim
E^{(1)}_c, \\[12pt] E
\left|\dfrac{S}{F}\right|\sqrt{\dfrac{s'_{{\rm
Re}{}}s''_{{\rm Re}{}}}{N'N''}
\left[1+\left(\dfrac{s'_{{\rm Im}{}}-s''_{{{\rm
Im}}{}}}{s'_{{\rm Re}{}}+s''_{{\rm
Re}{}}}\right)^2\right]}, & \textrm{for}\ E^{(1)}_c\lesssim
E\lesssim E^{(2)}_c \\[12pt]
E^{3/2}\left|\dfrac{S}{F}\right|\sqrt{\dfrac{v'v''}{N'N''}},
& \textrm{for}\ E\gtrsim E^{(2)}_c, \end{cases},
\end{equation}
where 
$E^{(i)}_c$ are the critical values given by the formulae
$$\begin{gather} E^{(1)}_c=\left|\frac{|\Delta k'q'\Delta k''q''|}{s'_{{\rm Re}{}}s''_{{\rm Re}{}}} \left[1+\left(\frac{s'_{{\rm Im}{}}-s''_{{\rm Im}{}}}{s'_{{\rm Re}{}}+s''_{{\rm Re}{}}}\right)^2\right]^{{-}1}\right|, \end{gather}$$
$$\begin{gather}E^{(2)}_c=\left|\frac{s'_{{\rm Re}{}}s''_{{\rm Re}{}}}{v'v''} \left[1+\left(\frac{s'_{{\rm Im}{}}-s''_{{\rm Im}{}}}{s'_{{\rm Re}{}}+s''_{{\rm Re}{}}}\right)^2\right]\right|. \end{gather}$$ In the experiments, the Ekman number is always close to 
$E_c^{(2)}$ such that the scaling exponent is between 1 and 3/2.
4.3. Centrifugal instability
For large tilt angles the flow exhibited an instability different to the triadic resonance instability studied above. For 
$r_i=0.33$, this instability was found at tilt angles larger than 
$3^\circ$, as shown by the green symbols of figure 8. In this case, the particle visualizations and the PIV measurements highlighted the presence of small toroidal vortices of alternate sign close to the inner cylinder (see figure 10a). This structure is reminiscent of a centrifugal instability (Rayleigh Reference Rayleigh1917; Taylor Reference Taylor1923) occurring in Taylor–Couette flows.
(a) Instantaneous PIV measurements of the meridional velocity and the toroidal vorticity for 
$E=3.6\times 10^{-4}$ during the onset of the instability. (b) Azimuthally averaged tangential velocity in the frame of reference of the laboratory at 
$\alpha =5^{\circ }$, 
$h=h_1=1.31$, 
$r_i=0.33$ and 
$z=0.33$ for 
$E=7.3\times 10^{-4}$ in green, 
$E=3.6\times 10^{-4}$ in red and 
$E=2.2\times 10^{-4}$ in blue. The black line corresponds to the solid body rotation of the annulus.
For an axisymmetric flow, it is well known that the centrifugal instability occurs when the Rayleigh discriminant becomes negative
\begin{equation} \frac{{\rm d}}{{\rm d}r}\left(r^2v_\theta^2\right)<0. \end{equation} In our case, the forced mode is not exactly axisymmetric, but we can use the azimuthal average of the velocity 
$\langle v\rangle _{\theta }$ to define a Rayleigh discriminant. In the limit of vanishing tilt angles, the solid body azimuthal profile is equal to 
$r$. This prediction is plotted as a black line in figure 10(b). For large Ekman number (i.e. when the resonant Kelvin mode is strongly damped), the experimental value is close to this prediction. However, when the Ekman number decreases, the experimental value decreases because the resonant Kelvin mode decelerates the solid body rotation by self-nonlinear interaction (Albrecht et al. Reference Albrecht, Blackburn, Lopez, Manasseh and Meunier2021; Gao et al. Reference Gao, Meunier, Le Dizès and Eloy2021). The value of 
$v_\theta$ is thus a decreasing function of 
$r$ close to the inner cylinder (between 
$r=0.33$ and 
$r=0.4$) leading to a negative Rayleigh discriminant. We can thus expect a centrifugal instability close to the inner cylinder, which is consistent with the alternate toroidal vortices observed in the PIV measurements. It should be noted that the value of 
$v_\theta$ was not measured by PIV on the inner cylinder but it must be equal to 
$r_i$ for a no-slip boundary condition.
4.4. Turbulent flow regime
When the Ekman number 
$E$ is decreased further below the threshold, the flow becomes more and more turbulent. As found for a cylinder (Lopez & Marques Reference Lopez and Marques2016), the amplitude of the free modes first become time dependent. Then other azimuthal Fourier modes appear. It leads to smaller structures that can be observed (see figure 11a), although superimposed to the vorticity of the forced mode (positive on one side and negative on the other side). The maximal vorticity is close to three.
Turbulent fluctuations characterized by (a,d) the axial vorticity, (b,e) the strain rate and (c, f) the horizontal divergence. (a–c) Examples of instantaneous fields are shown for 
$\alpha =3^{\circ }$ and 
$E=6.9\times 10^{-5}$. (d–f) The root mean square (r.m.s.). values are plotted for 
$\alpha =1.5^{\circ }$ in black, 
$\alpha =3^{\circ }$ in red and 
$\alpha =10^{\circ }$ in blue. The solid lines correspond to the stable forced mode with the forced Kelvin mode at an amplitude given by (3.17). Here 
$r_i=0.33$; 
$h=h_1$.
Figure 11(b) shows the strain rate 
$s$ defined as the largest eigenvalue of the horizontal symmetric velocity gradient tensor
\begin{equation} B_{ij}=\tfrac{1}{2}\big[\partial_iu_j+\partial_ju_i\big],\quad i,j=1,2. \end{equation}The small-scale structures of strain are smaller than their vortical counterparts with a maximum shear of the order of one.
Finally, figure 11(c) shows the two-dimensional (2-D) divergence 
$\partial _x u + \partial _yv$ of the flow which exhibits again small-scale structures with a maximum value close to two.
In order to measure quantitatively the turbulent fluctuations, we have calculated the r.m.s. value of the vorticity (see figure 11d). It increases when the Ekman number decreases and saturates at a value close to one. These experimental values are larger by a factor of two than the theoretical prediction in the stable case. However, it was already observed for the forced mode that the vorticity maximum was larger in the experiments and the DNS than in the theory (see figures 5 and 3).
The r.m.s. value of the strain rate is plotted in figure 11(e). It also increases when the Ekman number decreases and saturates at a value close to 0.5. These experimental values are again larger by a factor of two than the theoretical prediction in the stable case.
Finally, the r.m.s. value of the 2-D divergence is plotted in figure 11( f). Again, it increases when the Ekman number decreases and saturates at a value close to 0.5.
The values found here for an annulus are larger by a factor of almost three than for a tilted cylinder. This is possibly due to the radial confinement, which creates smaller structures and thus increases the velocity gradients. Fortunately, these values are still one order of magnitude smaller than for a Rushton turbine, consisting of six vertical blades and a rotating horizontal disk (Nagata Reference Nagata1975) at the same angular velocity. This tilted annulus is thus possibly a good soft mixer if the mixing is sufficiently fast.
5. Mixing
To study quantitatively the mixing properties of this flow, we injected a blob of dye inside the fluid while the annulus is tilted. The dye is stretched and folded by the flow, leading to thin streaks of dye (see figure 12a). When the distance between two streaks is smaller than the diffusive thickness, the concentration becomes uniform if the streaks are homogeneously spread over the whole volume.
(a) Example of observed pattern of the injected dye inside the fluid after one rotation. (b) Temporal evolution of the variance of concentration 
$c_{rms}/c_0$ for 
$h=h_{1}=1.306$ in red and 
$h=1.57$ in blue, both with 
$r_i = 0.33$, 
$\alpha =3^{\circ }$ and 
$E = 2.1\times 10^{-5}$. The two black lines are exponential fits 
$c_0{\rm e}^{-t/\tau }$ with 
$\tau =8$ for 
$h=h_{1}$ and 
$\tau =170$ for 
$h=1.57$.
5.1. Variance of concentration
The efficiency of the mixing can be measured quantitatively by calculating the variance of the concentration
\begin{equation} \frac{c_{rms}}{c_0} = \frac{\sqrt{\langle c^2\rangle-\langle c\rangle^2}}{c_0}. \end{equation} Indeed, the variance decreases with time and eventually vanishes when the dye concentration is uniform. The variance is plotted in figure 12(b) for two different heights of fluid. It is clear that the variance decays much faster at the resonant height 
$h=h_1$ (red symbols) than at a height only 15 % larger (blue symbols). This highlights the fact that this flow is very sensitive to the height of the fluid due to the resonance of the forced Kelvin mode.
Figure 12(b) also indicates that the concentration variance decreases exponentially with time. This is classical in turbulent flows and it can be explained by looking at a single lamella subject (Villermaux Reference Villermaux2019) to a straining field 
$\gamma$. The length in the direction of the strain is proportional to 
${\rm e}^{\gamma t}$, which by conservation of mass leads to an exponential contraction of the thickness as 
${\rm e}^{-\gamma t}$ at early stages. When the thickness reaches the Batchelor scale 
$\sqrt {\kappa /\gamma }$ (Batchelor Reference Batchelor1959), the lamella fades away. The thickness remains constant and the maximum concentration 
$c$ must decay as 
${\rm e}^{-\gamma t}$ for mass conservation. The concentration variance is thus proportional to the length times the maximal amplitude squared, which decreases as 
${\rm e}^{-\gamma t}$. The characteristic time of mixing can thus be characterized empirically by the decay rate of 
$c_{rms}/c_0$, which is fitted by 
${\rm e}^{-t/\tau }$ for each experiment. In order to get a systematic picture of the efficiency of mixing, the time 
$\tau$ has been measured for different heights of fluids 
$h$, tilt angles 
$\alpha$ and Ekman numbers.
5.2. Decaying time as a function of $h$
A set of experiments is conducted at an Ekman number 
$E = 2.1 \times 10^{-5}$ sufficiently low for the instability to occur at the resonance (see figure 13). As explained above, the mixing is much more efficient at the resonance of the forced Kelvin mode (
$h=h_1$) than outside of the resonance. For an aspect ratio 
$r_i=0.33$, the decaying time 
$\tau$ decreases by a factor of 
$20$ at 
$h=h_{1}$ compared with its average value outside of the resonance. This band of fast mixing has a width approximately equal to 
$20\,\%$ of the resonant height. The faster decaying time is approximately equal to 10, which means that the variance decays by a factor 
$e$ in less than two rotation periods.
Decaying time 
$\tau$ of the concentration variance as a function of 
$h$ for (a) 
$r_i=0.33$ and (b) 
$r_i=0.56$. The red dashed lines indicate the theoretical resonance heights of the first mode. Here 
$E=2.1\times 10^{-5}$; 
$\alpha =3^{\circ }$.
For a larger aspect ratio 
$r_i = 0.56$ (see figure 13b), the decaying time also has two minima at the resonances of the first Kelvin mode (
$h=h_1$ and 
$h=2 h_1$). It also seems to have another weaker minimum around 
$h=1.2$, which could correspond to the third resonance of the second Kelvin mode (
$h=3 h_2=1.16$). For this aspect ratio 
$r_i = 0.56$, the decaying time 
$\tau$ at the resonance is generally larger by a factor of two than for 
$r_i = 0.33$. It may come from the weaker unstable flow due to the larger Ekman damping (see § 4). It may also come from the narrow shape of the annulus, which prevents an efficient advection of the dye from the top to the bottom.
5.3. Decaying time as a function of $E$
The effect of the Ekman number is then studied for various tilt angles 
$\alpha$ and at the resonant heights (see figure 14). For all experiments, the decaying time 
$\tau$ increases with the Ekman number. For example for 
$r_i=0.33$, 
$h=h_1$ and 
$\alpha =3^{\circ }$ (see figure 14a), the decaying time increases rapidly from 
$\tau \sim 10$ to 
$\tau =100$ when the Ekman number increases by factor of five from 
$E = 2\times 10^{-5}$ to 
$E \sim 10^{-4}$. The decaying time then seems to saturate when the Ekman number increases further. It should be noted that the flow becomes unstable at a critical value 
$E=3 \times 10^{-4}$ for these parameters. It means that the Ekman number must be at least 5 to 10 times smaller than the critical Ekman number for the flow to be sufficiently turbulent to create an efficient mixing. This behaviour is quantitatively similar to the case of a tilted cylinder (Meunier Reference Meunier2020), where it was shown that it mixes as fast as using a Rushton turbine with a mixing time 
$\tau \sim 10$. For a larger tilt angle 
$\alpha =10^{\circ }$, the decaying time 
$\tau$ evolves similarly than for 
$\alpha = 3^{\circ }$ but it is 10 times faster. This result is expected since the amplitude of the turbulent fluctuations increases with 
$\alpha$. At Ekman numbers smaller than 
$10^{-4}$ the decaying time seems to reach a plateau although experimental limitations prevent the study of smaller Ekman numbers.
Decaying time 
$\tau$ of the concentration variance as a function of 
$E$ for (a) 
$r_i=0.33$ and 
$h=h_1$, (b) 
$r_i=0.56$ and 
$h=h_1$ and (c) 
$r_i=0.56$ and 
$h=2h_1$. The symbols display different angles 
$\alpha =1^{\circ }$ (
$+$, red), 
$\alpha =3^{\circ }$ (
$+$, blue) and 
$\alpha =10^{\circ }$ (
$+$, green). The vertical dashed lines are positioned at the critical 
$E$ value for the various experiments passing from stable (on the right-hand side) to unstable (on the left-hand side).
Changing the aspect ratio from 
$r_i=0.33$ to 
$r_i=0.56$ does not change the behaviour of the results. But it leads to a decrease in efficiency as already observed in figure 13, with an increase by a factor of at least two in the decaying time. Finally, changing the height of the water from 
$h=h_1$ to 
$h=2h_1$ further increases the decaying time by a factor of approximately two.
To conclude, the mixing is as efficient as in a tilted cylinder or using a Rushton turbine for a moderate aspect ratio 
$r_i=0.33$ but slightly weaker for a large aspect ratio 
$r_i=0.56$.
6. Conclusion
The flow inside a rotating annulus tilted with respect to gravity has been characterized theoretically and experimentally. At low Froude numbers the free surface remains horizontal, which creates an anticlockwise forcing of azimuthal wavenumber 
$m=1$ in the frame rotating with the annulus. As already observed by Xu & Harlander (Reference Xu and Harlander2020), it excites global inertial modes which can become resonant if the height of fluid is equal to a multiple of the half-wavelength of the mode. Indeed, the excitation is then located at the node of the mode, such that the mode must have an infinite amplitude to satisfy the boundary condition imposed by the free surface.
The linear and viscous theory used for the cylinder (Meunier Reference Meunier2020) can be extended to the case of the annular cavity. The forced mode has been described analytically with no fitting parameter and is in excellent agreement with the experiments. At low Ekman numbers, the flow becomes unstable thanks to a triadic resonance instability where the forced mode triggers two free global modes of azimuthal wavenumbers 
$m=9$ and 
${m=10}$, as already observed experimentally by Xu & Harlander (Reference Xu and Harlander2020). We have developed a linear stability analysis which predicts the viscous threshold of the instability. It is in excellent agreement with the experiments despite the lack of fitting parameters. As far as we know, it is the first time that such an analysis has been done in a rotating annulus. The annular geometry also generates a different type of instability which seems to be related to a centrifugal instability. Indeed, at large tilt angles, the solid body rotation is slowed down by the forced Kelvin mode. The inner cylinder tends to accelerate the fluid, leading to a centrifugally unstable forced mode in a narrow region around the inner cylinder.
The mixing efficiency of this flow is then characterized experimentally by injecting small blobs of dye. The characteristic mixing time (defined as the decaying time of the variance of concentration) is as fast as in a tilted cylinder or using a Rushton turbine at the same angular velocity for a moderate aspect ratio 
$r_i=0.33$ (between the inner and outer radius). The mixing time is slightly slower for a larger aspect ratio 
$r_i=0.56$. Since the shear is slightly larger (by a factor of two) in a tilted annulus than in a tilted cylinder, it seems that this annular mixer is not as smooth as a tilted cylinder. However, the shear is still an order of magnitude smaller than using a Rushton turbine for the same characteristic mixing time.
This annular mixer is thus an excellent candidate for large photobioreactors which are used for example for the growth of microalgae. Indeed, it allows for a simultaneous illumination from the outer and the inner cylinders. It thus increases the total amount of light delivered to the photobioreactor while remaining below the maximum intensity that creates photobleaching of the microalgae. This geometry is especially efficient for large concentrations of microalgae since the light is strongly attenuated after a few centimetres. For such concentrations, a cylindrical photobioreactor would be limited to diameters of the order of 10 to 20 cm, corresponding to volumes of the order of 5 to 50 l. With an annular photobioreactor, it is possible to reach volumes of the order of 600 litres with inner and outer radii equal to 50 and 70 cm, respectively (for 
$h=2h_1$).
Funding
This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska–Curie grant agreement no. 956457 and from SATT Sud-Est.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Amplitude viscous correction
We have seen that the amplitudes in the linear theory are divergent. This is obviously a non-physical behaviour. In fact the viscous effects saturate the resonance at a finite value of the amplitude (for details, see Kerswell & Barenghi Reference Kerswell and Barenghi1995). This effect due to the Ekman pumping (Greenspan Reference Greenspan1968) near the walls where there is the Ekman boundary layer (Ekman Reference Ekman1905). In the layer, of depth 
$\propto \sqrt {E}$, there is an additional motion due to the viscous damping of order 
$\sqrt {E}$:
\begin{equation}
\tilde{\boldsymbol{u}}\boldsymbol{\cdot}\hat{\boldsymbol{n}}=\sqrt{E}\left[-\frac{\delta_-+\delta_+}{2}(\hat{\boldsymbol{n}}\boldsymbol{\cdot}\boldsymbol{\nabla})(\boldsymbol{u_b}
\boldsymbol{\cdot}\hat{\boldsymbol{n}})+\frac{\delta_--\delta_+}{\varpi}
\partial_{z}(\hat{\boldsymbol{n}}\boldsymbol{\cdot}\boldsymbol{u_b})\right], \end{equation}
where we have defined
\begin{equation}
\delta_+=\begin{cases}
\dfrac{1-i}{\sqrt{2(\varpi-f)}}, & {\rm for}\ \varpi > f
\\[9pt] \dfrac{1+i}{\sqrt{2|\varpi-f|}}, & {\rm for}\
\varpi < f \end{cases},\quad \delta_- =\begin{cases}
\dfrac{1-i}{\sqrt{2(\varpi+f)}}, & {\rm for}\ \varpi + f >
0 \\[9pt] \dfrac{1+i}{\sqrt{2|\varpi+f|}}, & {\rm for}\
\varpi + f < 0 \end{cases} \end{equation}
and 
$f=2 \cos (\hat {n}{\cdot }\hat {\varOmega })$.
To see that we consider the solubility condition, obtained by multiplying the N–S equation with the complex conjugate of the 
$j$th mode 
$\bar {\boldsymbol {u}}_j=(\bar {u}_j\cos (k_jz), \bar {v}_j\cos (k_jz), \bar {w}_j\sin (k_jz)){\rm e}^{-{\rm i}(\theta +t)}$ and integrating over the cylinder excluding Ekman boundary layers (
$\tilde {V}$)
\begin{equation} \iiint\limits_{\tilde{V}}\bar{\boldsymbol{u}}_j\boldsymbol{\cdot} \left[{\rm i}\boldsymbol{u}+2\hat{\boldsymbol{z}}\wedge\boldsymbol{u}+\boldsymbol{\nabla}p-E {\nabla}^2\boldsymbol{u}\right] {\rm d}^3r=0. \end{equation} To find a better expression we consider the scalar product of 
$\boldsymbol {u}$ and the N–S equation for 
$\bar {\boldsymbol {u}}_j$
\begin{equation} -{\rm i}\boldsymbol{u}\boldsymbol{\cdot}\bar{\boldsymbol{u}}_j+2 \boldsymbol{u}\boldsymbol{\cdot}\hat{\boldsymbol{z}}\wedge\bar{\boldsymbol{u}}_j={-} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\bar{p}_j, \end{equation}
the symmetry 
$2\boldsymbol {u}\boldsymbol{\cdot}\hat {\boldsymbol {z}}\wedge \bar {\boldsymbol {u}}_j=-2 \bar {\boldsymbol {u}}_j\boldsymbol{\cdot}\hat {\boldsymbol {z}}\wedge \boldsymbol {u}$ and the incompressibility condition to have (the surface integrals are carried out by considering 
$\boldsymbol {n}$ oriented towards the inside of the annular cavity)
\begin{equation} \oint_{\partial\tilde{V}}\left(\bar{p}_j\boldsymbol{u}+p\bar{\boldsymbol{u}}_j\right) \boldsymbol{\cdot}\boldsymbol{n}\,{\rm d}\varSigma-E\iiint\limits_{\tilde{V}} \bar{\boldsymbol{u}}_j\boldsymbol{\cdot}{\nabla}^2\boldsymbol{u}\,{\rm d}^3x=0. \end{equation} First consider the top surface (free surface) we have from the first term (
$\propto \bar {p}_ju$)
\begin{equation} -\int_0^{2{\rm \pi}}\int_{r_i}^{1}{\rm d}\theta r\,{\rm d}r\left[\bar{p}_j(r,z=h){\rm e}^{-{\rm i}(\theta+t)}\left(\alpha r \sin(\theta+t)\right)\right]\equiv{-}\alpha F, \end{equation}
and the second term (
$\propto p\bar {u}_j$)
\begin{equation} -\int_0^{2{\rm \pi}}\int_{r_i}^{1}{\rm d}\theta r\,{\rm d}r\left[p(r,z=h)\bar{w}_j(r,z=h){\rm e}^{-{\rm i}(\theta+t)}\right]\equiv A_j\sin(k_j h) T. \end{equation}The bottom contribution comes only from the first term
\begin{equation} \int_0^{2{\rm \pi}}\int_{r_i}^{1}{\rm d}\theta r\,{\rm d}r \left[\bar{p}_j(r,z=0){\rm e}^{-{\rm i}(\theta+t)}\tilde{\boldsymbol{u}}\boldsymbol{\cdot}\hat{n}\right], \end{equation}and using the expression of the Ekman pumping (A1)
\begin{equation} \int_0^{2{\rm \pi}}\int_{r_i}^{1}{\rm d}\theta r\,{\rm d}r \left[\bar{p}_j(r,z=0){\rm e}^{-{\rm i}(\theta+t)}\left(\frac{\sqrt{E}}{2}\left.\left({-}3\delta_{+}+\delta_{-}\right) \partial_zw\right|_{z=0}\right)\right] \equiv A_j \sqrt{E}B, \end{equation}
where 
$\delta _{+}=(1+i)/2$ and 
$\delta _{-}=(1-i)/\sqrt {6}$.
The same happens for the inner and outer cylinders,
$$\begin{gather} \int_0^{2{\rm \pi}}\int_{0}^{h}{\rm d}\theta \,{\rm d}z\left[\bar{p}_j(r=1,z){\rm e}^{-{\rm i}(\theta+t)} \left(\left.-\delta\partial_ru\right|_{r=1}\right)\right]\equiv A_j \sqrt{E}C_e, \end{gather}$$
$$\begin{gather}\int_0^{2{\rm \pi}}\int_{0}^{h}{\rm d}\theta\, {\rm d}z\left[\bar{p}_j(r=r_i,z){\rm e}^{-{\rm i}(\theta+t)}\left(-\left.\delta\partial_ru\right|_{r=ri}\right)\right]\equiv A_j \sqrt{E}C_i, \end{gather}$$
where 
$\delta =(1-i)/\sqrt {2}$.
Finally, the last term (volume integral) has an easy analytical expression
\begin{equation} \iiint\limits_{\tilde{V}}\bar{\boldsymbol{u}}_j\boldsymbol{\cdot}{\nabla}^2 \boldsymbol{u}\,{\rm d}^3x={-}4k_j^2\iiint\limits_{\tilde{V}}\bar{\boldsymbol{u}}_j\boldsymbol{\cdot}\boldsymbol{u} \,{\rm d}^3x\equiv A_j E V. \end{equation}With these definitions we can write the equation for viscous amplitude at the resonance
\begin{equation} A_j=\frac{\alpha F}{T\sin(k_j h)+\sqrt{E}\left(B+C_e+C_i\right)+ E V}. \end{equation}All the corrections calculated are 
$\propto E$ and in the limit 
$E\to 0$ we have the inviscid theory (see (3.13)) and from the definition it is easy to prove it.
These coefficients can be given analytically as
\begin{equation} C_e=\frac{2 \sqrt{2} h(1+i)(1+k_j^2)}{\rm \pi} \left[ 1+ \mathrm{sinc}(2 h k_j) \right] \end{equation}
with 
$\mathrm {sinc}(x)=\sin (x)/x$,
$$\begin{gather} C_i=\frac{\sqrt{2} h(1+i)(1+k_j^2 r_i^2) c_j^J \left[ c_j^J {\rm J}_1(\sqrt{3} k_j r_i) + c_j^Y {\rm Y}_1(\sqrt{3} k_j r_i) \right]} {r_i^2 \left[ {\rm Y}_1(\sqrt{3} k_j r_i) + \sqrt{3} k_j r_i {\rm Y}_0(\sqrt{3} k_j r_i) \right]} \left[ 1+ \mathrm{sinc}(2 h k_j) \right], \end{gather}$$
$$\begin{gather}F=2 {\rm i} {\rm \pi}\frac{ c_j^J \left[{\rm J}_2(\sqrt{3}k_j)-r_i^2 {\rm J}_2(\sqrt{3}k_j r_i) \right] + c_j^Y \left[{\rm Y}_2(\sqrt{3}k_j)-r_i^2 {\rm Y}_2(\sqrt{3}k_j r_i) \right] }{\sqrt{3}k_j} \cos(k_j h), \end{gather}$$
$$\begin{gather}T=2 {\rm i} {\rm \pi}k_j \cos(k_j h) I, \end{gather}$$
$$\begin{gather}B= {\rm i} {\rm \pi}k_j^2 \left( \frac{1-i}{\sqrt{6}} - 3 \frac{1+i}{\sqrt{2}}\right) I, \end{gather}$$with
\begin{equation} I= \int_{r=r_i}^1 \left[ c_j^J {\rm J}_1(\sqrt{3} k_j r) + c_j^Y {\rm Y}_1(\sqrt{3} k_j r) \right]^2 \,{\rm d}r. \end{equation}Appendix B. Amplitude equations for triadic instability
This result is achieved by doing a linear stability analysis, thanks to the fact that we have the small parameters 
$\alpha$ and 
$E$, of the complete N–S equation. To perform the linear stability analysis a perturbation is added to the forced Kelvin mode
\begin{equation} \boldsymbol{u}=\boldsymbol{u}_f+A'\boldsymbol{u}'+A''\boldsymbol{u}'', \end{equation}
where 
$\boldsymbol {u}_f$ is the forced Kelvin mode given in (3.6a,b) and 
$\boldsymbol {u}'$, 
$\boldsymbol {u}''$ are the two free Kelvin modes triggered by the forced mode.
In general, the free Kelvin modes are defined by (Xu & Harlander Reference Xu and Harlander2020)
\begin{equation} \boldsymbol{u}'= {\rm Re}\left[A'\begin{pmatrix} u'(r)\cos(k'z) \\ v'(r)\cos(k'z) \\ w'(r)\sin(k'z) \end{pmatrix}\exp({{\rm i}(m'\theta+ \varpi't)})\right],\end{equation}where
\begin{equation}
\left.\begin{array}{@{}c@{}}
p'(r)=c'^J{\rm J}_{m'}\left(\sqrt{\dfrac{4}{\varpi'^2}-1}k'r\right)+c'^Y{\rm Y}_{m'}\left(\sqrt{\dfrac{4}{\varpi'^2}-1}k'r\right),\\
u'(r)={-}\dfrac{{\rm i}}{4-\varpi'^2}\left(\varpi'\partial_rp'+\dfrac{2m'}{r}p'\right),\\
v'(r)=\dfrac{1}{4-\varpi'^2}\left(2\partial_rp'+\dfrac{\varpi'm'}{r}p'\right),\\
w'(r)={-}{\rm i}\dfrac{k'}{\varpi'}p',\end{array}\right\}\end{equation}and
\begin{align}
\left.\begin{aligned}
c'^J&=\dfrac{1}{2}\sqrt{\dfrac{4}{\omega'^2}-1}k'\varpi'\left[{\rm
Y}_{m'-1} \left(\sqrt{\dfrac{4}{\omega'^2}-1}k'\right)-{\rm
Y}_{m'+1}
\left(\sqrt{\dfrac{4}{\omega'^2}-1}k'\right)\right]\\
&\quad +\,2m'{\rm Y}_{m'} \left(\sqrt{\dfrac{4}{\omega'^2}-1}k'\right),\\
c'^Y&=\dfrac{1}{2}\sqrt{\dfrac{4}{\omega'^2}-1}k'\varpi' \left[{\rm J}_{m'+1}\left(\sqrt{\dfrac{4}{\omega'^2}-1}k'\right)-{\rm
J}_{m'-1} \left(\sqrt{\dfrac{4}{\omega'^2}-1}k'\right)\right]\\
&\quad -\,2m'J{_m}'\left(\sqrt{\dfrac{4}{\omega'^2}-1}k'\right).\end{aligned}\right\}\end{align} At the resonance of the 
$j$th forced Kelvin mode the amplitude scales as 
$A_j\sim \alpha E^{-1/2}$ and is taken as the parameter for the analysis. Assuming a small dependence in time on 
$A'$ and 
$A''$ (
$\partial _tA'\sim A_j$) and that the two free modes are not exactly resonant (
$k'=k'_{{res}}+\Delta k'$ with 
$\Delta k\sim O({A_j})$), we find the evolution equation (Lagrange et al. Reference Lagrange, Meunier, Nadal and Eloy2011)
\begin{equation}
\left.\begin{array}{@{}c@{}}
\partial_tA'=N'\bar{A}_jA''-D' A',\\
\partial_tA''=N'' A_j A'-D'' A'', \end{array}\right\}
\end{equation}
where the nonlinear terms are (we want to stress the fact that this result is different from the one found in Lagrange et al. (Reference Lagrange, Meunier, Nadal and Eloy2011) because we use a different basis)
$$\begin{gather} N'= \frac{-{\rm i}\left(\dfrac{k''}{\varpi''}-k\right) {\rm I}_{{coupling}} }{\displaystyle\int\nolimits_{r_i}^{1}\left( |u'|^2 +|v'|^2 +|w'|^2 \right) r\,{\rm d}r } , \end{gather}$$
$$\begin{gather}N''=\frac{{\rm i}\left(k-\dfrac{\pm k'}{\varpi'}\right) {\rm I}_{{coupling}} }{ \displaystyle\int\nolimits_{r_i}^{1}\left( |u''|^2 +|v''|^2 +|w''|^2 \right) r\,{\rm d}r } , \end{gather}$$with
\begin{equation} {\rm I}_{{coupling}} = \int_{r_i}^{1} \begin{vmatrix} u'' & u' & u \\ -v'' & v' & v \\ -w'' & \pm w' & w \\ \end{vmatrix} r\,{\rm d}r, \end{equation}
where 
$|.|$ is the determinant and the plus sign corresponds to the case 
$k''=|k+k'|$ and the minus sign corresponds to the case 
$k''=|k-k'|$.
These formulae can be linked to the values found in the literature in the case of the cylinder taking the limit 
$r_i\to 0$. In this limit the coefficient 
$c^{Y}$ vanishes such that the solution (3.7) is equivalent to the classical solution in a cylinder (Meunier Reference Meunier2020) except that it is 
$c^J/2$ times larger. It should be noted that the boundary condition 
$u=0$ at 
$r=r_i$ virtually disappears when 
$r_i$ tends to 0 since this boundary condition creates a perturbation of the flow only in a small region of size 
$r_i$, that disappears when 
$r_i$ tends to 0. The Kelvin modes of an annulus with very small 
$r_i$ are thus exactly the same as the Kelvin modes of a cylinder.
The nonlinear coefficients given by Meunier (Reference Meunier2020) are related to our coefficient by the relation
\begin{equation} \left.\begin{array}{@{}c@{}} N'_{{Meu}} = \lim_{r_i\to0} \dfrac{i2c^J_5}{c^J_6c^J_1} N',\\ N''_{{Meu}} = \lim_{r_i\to0} \dfrac{-i2c^J_6}{c^J_5c^J_1} N'',\end{array}\right\} \end{equation}
due to the observation 
$c^Y\ll c^J$ for 
$r_i\to 0$.
The dispersion relation gives the resonance for 
$m'=5$ and 
$m''=6$ with resonant height 
$h=1.99$, same value of the cylinder, and nonlinear coefficient
\begin{equation} \left.\begin{array}{@{}c@{}} N'_{{Meu}} ={-}1.26,\\ N''_{{Meu}} ={-}1.95,\end{array}\right\} \end{equation}
that gives a value of the growth rate with an error of 
$2\,\%$.
The other term
\begin{equation} D'=s'\sqrt{E}+v' E+{\rm i} \Delta k' q', \end{equation}describes the viscous corrections
\begin{gather} s'=\frac{\begin{array}{l}k'\left[-\left(\dfrac{1}{2}\left(\delta'_{+}+\delta'_{-}\right)+ \dfrac{1}{\varpi'}\left(\delta'_{+}-\delta'_{-}\right)\right)\right] \\ \quad \times \displaystyle\int\nolimits_{r_i}^{1}p'w'r\,{\rm d}r-\dfrac{\delta'}{2h} \left[\left.\left(p'\partial_{r}u'\right)\right|_{(r=r_i)} +\left(p'\partial_{r}u'\right)\big|_{(r=r_i)}\right]\end{array}}{\displaystyle\int\nolimits_{r_i}^{1}\left( |u'|^2 +|v'|^2 +|w'|^2 \right) r\,{\rm d}r}, \end{gather}
\begin{gather} v'=\frac{4k'^2}{\varpi'^2} \end{gather}and the detuning effect
\begin{equation} q'={-}{\rm i}k'\frac{\displaystyle\int\nolimits_{r_i}^{1}w'p'r\,{\rm d}r}{\displaystyle\int\nolimits_{r_i}^{1} \left( |u'|^2 +|v'|^2 +|w'|^2 \right) r\,{\rm d}r}. \end{equation}The numerical values of the correction terms used for comparison with the experiments are given in table 3.
Numerical values of the viscous terms for the triadic resonance, for 
$r_i=0.33$ (
$m'=9, m''=10$) and 
$r_i=0.56$ (
$m'=17, m''=18$) at the resonance 
$h=h_{1}$ and 
$h=2h_{1}$.
Appendix C. Small gap limit
We want to study the solution in the small gap limit (
$r_i\to 1$ or 
$\epsilon \to 0$, where 
$\epsilon =1-r_i$ is the gap). In this limit we can rewrite the pressure equation (3.5) with the rescaling 
$\tilde {r}=(r-1)/\epsilon$. In this variable the term 
$r^{-1}\partial _r$ is subdominant with respect to 
$\partial ^2_r$, so under the assumption of large wavenumber 
$k$ and 
$m$ (which are of order 
$\epsilon ^{-1}$) the equation becomes (when coming back to the variable 
$r$)
\begin{equation} (\partial^2_{r}+\lambda^2)p(r)=0,\quad \lambda=\sqrt{k^2\left(\frac{4}{\varpi^2}-1\right)-m^2}.\end{equation}
Adding the boundary conditions 
$u(1)=u(1-\epsilon )$ gives the solution
\begin{equation} p(r)=c\left(\cos[\lambda(r-1)]-\frac{2m}{\varpi\lambda}\sin[\lambda(r-1)]\right),\quad \sin(\lambda\epsilon)=0.\end{equation}We want to connect this solution with the forced mode in (3.7) by having the same amplitude in the limit of small 
$\epsilon$. Taking the asymptotic formula for 
${\rm J}_m$ and 
${\rm Y}_m$, in (B4) and (B3), we can show that the two relations are equal if 
$c=2\varpi /{\rm \pi}$.
Finally, using the relationships between flow and pressure, we find the complete solution
\begin{equation}
\left.\begin{array}{@{}c@{}}
p=\dfrac{2}{{\rm \pi}\lambda}\left(\varpi\lambda\cos[\lambda(r-1)]-2m\sin[\lambda(r-1)]\right),\\
u=\dfrac{2 {\rm i} (\lambda^2 \varpi^2+4 m^2)}{{\rm \pi} \lambda(4-\varpi^2)}\sin[\lambda(r-1)], \\
v={-}\dfrac{2}{{\rm \pi}\lambda(4-\varpi^2)} \left[m\lambda(4-\varpi^2) \cos[\lambda(r-1)]+2\varpi
(\lambda^2+m^2) \sin[\lambda(r-1)]\right],\\
w={-}{\rm i}\dfrac{2k}{{\rm \pi}\lambda\varpi}\left(\varpi\lambda\cos[\lambda(r-1)]-2m\sin[\lambda(r-1)]\right).
\end{array}\right\} \end{equation}
For the forced mode, 
$m=1$ and 
$\varpi =1$, the dispersion relation (C1a,b) gives the asymptotic axial wavenumber 
$k_j=j{\rm \pi} /\sqrt {3}\epsilon$, so the forced mode is resonant when 
$h$ is a multiple of half-wavelength 
$h_j={\rm \pi} /k_j$.
Assuming
\begin{equation}
\left.\begin{array}{@{}c@{}}
\varpi''-\varpi'=1,\\
m''-m'=1, \end{array}\right\}\end{equation}
we have the triadic resonant instability, with nonlinear coefficient 
$N'$ and 
$N''$ given by (B6)–(B7). Assuming the case 
$k''=|k+k'|$ for simplicity of description (at the end of the subsection we give the change in all the equations in the case 
$k''=|k-k'|$), these identities can be calculated for small 
$\epsilon$, by changing 
$r\,{\rm d}r\to dr$ in all integrals, and are
$$\begin{gather} N'=\frac{{\rm i}\left(\dfrac{k''}{\varpi''}-k\right){\rm \pi}^2\lambda'^2\varpi'^2}{16k'^2(k'^2+m'^2)\epsilon} {\rm I}_{{coupling}}, \end{gather}$$
$$\begin{gather}N''=\frac{-{\rm i}\left(k-\dfrac{k'}{\varpi'}\right){\rm \pi}^2\lambda''^2\varpi''^2}{16k''^2(k''^2+m''^2) \epsilon}{\rm I}_{{coupling}}, \end{gather}$$
where 
${\rm I}_{{coupling}}$ can be also calculated analytically.
If 
$\lambda ''+\lambda '+\lambda$ is an odd multiple of the fundamental wavelength 
${\rm \pi} /\epsilon$:
\begin{align} {\rm
I}_{{coupling}}&=\frac{-16}{{\rm \pi} ^3(\lambda
-\lambda'-\lambda'') (\lambda +\lambda'-\lambda'') (\lambda
-\lambda'+\lambda'')(\lambda
+\lambda'+\lambda'')}\nonumber\\
&\quad \times \left[(k'm''-k''m')
\left((-\lambda^2+\lambda'^2+\lambda''^2)+\frac{8k'k''}{\varpi'\varpi''}\right)
\frac{(k^2+m^2)}{\varpi}\right. \nonumber\\
&\quad +(k''m-km'')\left((\lambda^2-\lambda'^2+\lambda''^2)+\frac{8
kk''}{\varpi\varpi''}\right)\frac{(k'^2+m'^2)}{\varpi'}\nonumber\\
&\left.\quad +\, (k'm-km')
\left((\lambda^2+\lambda'^2-\lambda''^2)+\frac{8
kk'}{\varpi\varpi'}\right)\frac{(k''^2+m''^2)}{\varpi''}\right].
\end{align} Then if 
$\lambda ''+\lambda '+\lambda$ is an even multiple of the fundamental wavelength 
${\rm \pi} /\epsilon$, 
${\rm I}_{{coupling}}=0$ unless 
$\lambda ''=|\lambda \pm \lambda '|$; in that case
\begin{align} {\rm
I}_{{coupling}}&=\frac{-4\epsilon}{{\rm \pi}^3} \left[
(k'k''+m'm'')
\left(\frac{k''}{\varpi''\lambda''}+\frac{k'}{\varpi'({\pm}\lambda')}\right)\frac{(k^2+m^2)}{\lambda}
\right.\nonumber\\ &\quad - (k k''+m m'')
\left(\frac{k''}{\varpi''\lambda''}
+\frac{k}{\varpi\lambda}\right)\frac{(k'^2+m'^2)}{({\pm}\lambda')}
\nonumber\\ &\left.\quad - \,(k k'+m m')
\left(\frac{k'}{\varpi'({\pm}\lambda')}-\frac{k}{\varpi\lambda}
\right)\frac{(k''^2+m''^2)}{\lambda''}\right],
\end{align}
where the plus sign corresponds to the case 
$\lambda ''=|\lambda +\lambda '|$ and the minus sign corresponds to the case 
$\lambda ''=|\lambda -\lambda '|$. In the case of 
$k''=|k-k'|$ all the equations written before are valid under the transformation 
$k'->-k'$.
Some examples from the above formulae are given in table 4.
Examples of growth rate in the asymptotic theory for 
$h={\rm \pi} /k$. All of them shows that for 
$\epsilon \lesssim 0.01$ the value of 
$\epsilon ^2\sigma$ becomes constant. The table shows some different cases in the triad of 
$\lambda$ satisfying the condition (C4) and 
$k''=k+k'$.
The interesting case for comparing asymptotic and exact theory is 
$\lambda ''=\lambda '=\lambda ={\rm \pi} /\epsilon$; we can plot the numerical values of the growth rate as a function of 
$\epsilon$ of the exact and asymptotic theory, as shown in figure 15.
Comparison of asymptotic (blue) and exact theory (red) in growth rate values. It can be seen that the exact theory is asymptotically equivalent to the case where 
$\lambda ''=\lambda '=\lambda ={\rm \pi} /\epsilon$.
