S Boury

I Sibgatullin

E Ermanyuk

N Shmakova

P Odier

S Joubaud

L R M Maas

T Dauxois

Vortex cluster arising from an axisymmetric inertial wave attractor

We present an experimental and numerical study of the non-linear dynamics of an inertial wave attractor in an axisymmetric geometrical setting. The rotating ring-shaped fluid domain is delimited by two vertical coaxial cylinders, a conical bottom, and a horizontal wave generator at the top: the vertical cross-section of the fluid volume is a trapezium, while the horizontal cross-section is a ring. Forcing is introduced via axisymmetric low-amplitude volume-conserving oscillatory motion of the upper lid. The experiment shows an important result: at sufficiently strong forcing and long time scale, a saturated fully non-linear regime develops as a consequence of an energy transfer draining energy towards a slow two-dimensional manifold represented by a regular polygonal system of axially-oriented cyclonic vortices undergoing a slow prograde motion around the inner cylinder. We explore the long-term non-linear behaviour of the system by performing a series of numerical simulations for a set of fixed forcing amplitudes. We observe a rich variety of dynamic regimes, including a linear behaviour, a Triadic Resonance Instability (TRI), a progressive frequency enrichment reminiscent of weak inertial Wave Turbulence (WT), and the generation of a slow manifold in the form of a Polygonal Vortex Cluster (PVC). This vortex cluster is discussed in detail, and we show that it stems from the summation and merging of wave-like components of the vorticity field. The nature of these wave components, the possibility of their detection under general conditions and the ultimate fate of the vortex clusters at even longer time scale remain to be explored.

Introduction

Energy transfer in rotating fluids has received significant attention due to its relevance to geo-and astrophysical fluid dynamics and due to the rich complexity of the non-linear multi-scale interplay between coherent vortical structures, inertial waves and background small-scale nearly isotropic turbulence [START_REF] Greenspan | The Theory of Rotating Fluids[END_REF][START_REF] Hopfinger | Vortices in rotating fluids[END_REF][START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF][START_REF] Godeferd | Structure and dynamics of rotating turbulence: A review of recent experimental and numerical results[END_REF]. Inertial waves supported by rotating fluids, with the Coriolis force acting as restoring force, represent an essential ingredient of the transfer. The crucial role of inertial waves is assured by [i] the specific form of the dispersion relation, which contains no length scale, and [ii] the possibility of a cascade of wavewave interactions due to non-linear terms in the Navier-Stokes equations governing the dynamics of rotating fluids.

The dispersion relation of inertial waves obtained by seeking plane-wave solutions of the linearized inviscid Navier-Stokes equations reads ω = f k z /k = f cos α, where ω is the wave frequency, f = 2Ω is the Coriolis parameter with Ω the rate of the background rigidbody rotation of the fluid, and k z (respectively k) is the vertical component (respectively magnitude) of the wave vector k inclined at angle α to the vertical z-axis, which is taken as the axis of rotation. A similar type of dispersion relation ω = N k h /k = N sin α holds for internal waves in a stratified fluid, with the buoyancy frequency N replacing f , and k h the horizontal wave number. The absence of any length scale in the dispersion relation for inertial and internal waves implies that the global large-scale wave pattern depends on the geometry of wave generators and on the geometry that delimits the fluid volume -in particular, for the ocean, its bathymetry. Therefore, a rich variety of wave motions is encountered in rotating and stratified fluids as identified in early pioneering studies, and explored in detail in the subsequent literature: normal modes in bounded domains of simple geometry (sphere, axial cylinder, rectangular box) [START_REF] Aldridge | Axisymmetric inertial oscillations of a fluid in a rotating spherical container[END_REF][START_REF] Mcewan | Inertial oscillations in a rotating fluid cylinder[END_REF][START_REF] Mcewan | Degeneration of resonantly excited standing internal gravity waves[END_REF][START_REF] Maas | On the amphidromic structure of inertial waves in a rectangular parallelepiped[END_REF][START_REF] Bewley | Inertial waves in rotating grid turbulence[END_REF][START_REF] Lamriben | Excitation of inertial modes in a closed grid turbulence experiment under rotation[END_REF][START_REF] Boisson | Inertial waves and modes excited by the libration of a rotating cube[END_REF], wave beams emanating from isolated oscillatory sources [START_REF] Görtler | Über eine schwingungserscheinung in flüssigkeiten mit stabiler dichteschichtung[END_REF][START_REF] Hendershott | Impulsively started oscillations in a rotating stratified fluid[END_REF][START_REF] Thomas | A similarity solution for viscous internal waves[END_REF][START_REF] Mowbray | A theoretical and experimental investigation of the phase configuration of internal waves of small amplitude in a density stratified liquid[END_REF], and webs of wave beams (wave attractors) in bounded domains with sloping walls [START_REF] Stern | Trapping of low frequency oscillations in an equatorial boundary layer[END_REF][START_REF] Bretherton | Low frequency oscillations trapped near the equator[END_REF][START_REF] Stewartson | On trapped oscillations of a rotating fluid in a thin spherical shell[END_REF][START_REF] Stewartson | On trapped oscillations of a rotating fluid in a thin spherical shell ii[END_REF][START_REF] Maas | Geometric focusing of internal waves[END_REF][START_REF] Maas | Observations of an internal wave attractor in a confined stably stratified fluid[END_REF][START_REF] Manders | Observations of inertial waves in a rectangular basin with one sloping boundary[END_REF][START_REF] Klein | Inertial wave excitation and focusing in a liquid bounded by a frustum and a cylinder[END_REF]Wu et al. 2020b). Of particular interest is the latter configuration in the context of the present paper.

Due to the form of the internal and inertial wave dispersion relations, wave reflection on a solid boundary follows a very specific law and is, in general, non-specular [START_REF] Phillips | Energy transfer in rotating fluids by reflection of inertial waves[END_REF][START_REF] Eriksen | Observations of internal wave re1ection on sloping bottoms[END_REF][START_REF] Manders | On the three-dimensional structure of the inertial wave field in a rectangular basin with one sloping boundary[END_REF][START_REF] Maas | Wave attractors: linear yet non linear[END_REF]. To be more specific, in twodimensional domains, this law leads to a focusing or a defocusing effect of wave beams upon reflection at sloping walls [START_REF] Dauxois | Near-critical reflection of internal waves[END_REF]. In bounded or quasi-bounded two-dimensional domains, focusing prevails: the iterative process of subsequent wave reflections leads to the formation of a limit cycle, called a wave attractor, where the wave energy is concentrated [START_REF] Maas | Geometric focusing of internal waves[END_REF][START_REF] Maas | Observations of an internal wave attractor in a confined stably stratified fluid[END_REF][START_REF] Sibgatullin | Internal and inertial wave attractors: A review[END_REF]. Relevant to the topic of the present study is the first experimental observation of an inertial wave attractor in an elongated trapezoidal tank that showed the generation of a persistent mean flow, right above the location where the attractor was being focused over the sloping bottom. This mean flow was speculated to be the result of the breaking of focused inertial waves, leading to the mixing of the background radial stratification in angular momentum with which the solidly-rotating, homogeneous-density fluid is endowed [START_REF] Maas | Wave focusing and ensuing mean flow due to symmetry breaking in rotating fluids[END_REF]. In a three-dimensional setting, the variety of possible configurations is significantly enriched, involving the possibility of wave-energy trapping on a limit cycle located at certain preferential planes of motion, provided that there is a billiard pathway connecting this plane and the initial direction of the wave-energy propagation [START_REF] Hazewinkel | Tomographic reconstruction of internal wave patterns in a paraboloid[END_REF][START_REF] Pillet | Internal wave attractors in three-dimensional geometries: trapping by oblique reflection[END_REF]. The inertial-wave-ray billiard corresponding to the geophysically important case of a rotating spherical shell favors the formation of an attractor in the meridional plane [START_REF] Bretherton | Low frequency oscillations trapped near the equator[END_REF][START_REF] Stewartson | On trapped oscillations of a rotating fluid in a thin spherical shell ii[END_REF][START_REF] Friedlander | Internal waves in a contained rotating stratified fluid[END_REF][START_REF] Maas | Equatorial wave attractors and inertial oscillations[END_REF][START_REF] Rabitti | Equatorial wave attractors and inertial oscillations[END_REF]. Accordingly, the rich literature on the linear dynamics of inertial wave attractors in rotating spherical layers considers the motions in ring-shaped meridional slices [START_REF] Friedlander | Internal waves in a contained rotating stratified fluid[END_REF][START_REF] Dintrans | Gravito-inertial waves in a rotating stratified sphere or spherical shell[END_REF][START_REF] Rieutord | Wave attractors in rotating fluids: a paradigm for ill-posed cauchy problems[END_REF][START_REF] Rieutord | Inertial waves in a rotating spherical shell: attractors and asymptotic spectrum[END_REF][START_REF] Rieutord | Viscous dissipation by tidally forced inertial modes in a rotating spherical shell[END_REF] and disregards the azimuthal coordinate.

It is noteworthy that the purely geometrical mechanism of iterative focusing, which is linear, is at the origin of a spectacular forward energy cascade in wave attractors: the energy injected into the system at global scale (i.e. at the scale of the system itself) is transferred to the scale corresponding to the width of the attractor branches, which, even in laboratory experiments, can be an order of magnitude smaller than the global scale (Brouzet et al. 2017b). This small scale, or width of the inertial and internal wave beams in the linear regime, is set by the balance between geometric focusing and viscous dissipation and can be theoretically predicted with good agreement to experimental observations [START_REF] Rieutord | Inertial waves in a rotating spherical shell: attractors and asymptotic spectrum[END_REF][START_REF] Hazewinkel | Observations on the wavenumber spectrum and evolution of an internal wave attractor[END_REF][START_REF] Grisouard | Numerical simulation of a two-dimensional internal wave attractor[END_REF]. Other wave-damping mechanisms such as interaction of waves with convective motions, ohmic damping in presence of magnetic field in conducting fluids, and non-linear parametric decay into secondary waves of shorter wavelength have also been proposed by [START_REF] Ogilvie | Wave attractors and the asymptotic dissipation rate of tidal disturbances[END_REF], where a generic case with a weak inviscid "frictional" damping force has been considered. Further, it has been shown experimentally that at sufficiently high level of injected energy, internal wave attractors are prone to Triadic Resonance Instability (TRI) [START_REF] Scolan | Nonlinear fate of internal waves attractors[END_REF]. The replacement of purely viscous damping by the flux of energy carried by small-scale secondary waves (generated via TRI) away from the primary waves (i.e. from the beams of attractor) introduces a new non-linear scaling for the beam width (Brouzet et al. 2017b). Similar effects in inertial waves have been observed in numerical simulations [START_REF] Jouve | Direct numerical simulations of an inertial wave attractor in linear and nonlinear regime[END_REF]. Let us note in passing that [START_REF] Jouve | Direct numerical simulations of an inertial wave attractor in linear and nonlinear regime[END_REF] considered a two-dimensional setting, physically corresponding to a torus of infinite radius having a tilted-square "meridional" cross-section, so that any three-dimensional effects occurring in "equatorial" planes were completely excluded. The development of the energy cascade in wave attractors with the increase of injected energy leads to wave turbulence, with a significant occurrence of overturning events generating irreversible mixing (Brouzet et al. 2016a(Brouzet et al. , 2017a;;[START_REF] Davis | Succession of resonances to achieve internal wave turbulence[END_REF], and such a cascade reaches a statistically steady state when a balance is established between the injected and dissipated energy [START_REF] Jouve | Direct numerical simulations of an inertial wave attractor in linear and nonlinear regime[END_REF][START_REF] Davis | Energy budget in internal wave attractor experiments[END_REF]. This is to be put in context with some recent development on the understanding of inertial wave turbulence in rotating flows (see, for example, the recent works of Le [START_REF] Le Reun | Experimental study of the non-linear saturation of the elliptical instability: inertial wave turbulence versus geostrophic turbulence[END_REF] and [START_REF] Brunet | Shortcut to geostrophy in wave-driven rotating turbulence: the quartetic instability[END_REF] on the competition between the saturation of rotating turbulence in 3D wave turbulence and in 2D geostrophic turbulence).

It should be stressed that the overturning events and subsequent mixing are important constituents of the full energy cascade in internal wave attractors, and that they clearly fall apart from the wave turbulence formalism. Similarly, inertial wave turbulence plays an important but not exclusive role in the non-linear dynamics of inertial wave attractors described in the present paper. The full scope of dynamic events constituting the energy transfer in inertial wave attractors extends well beyond the wave turbulence framework and should be discussed in the rich context of the literature on turbulence in rotating fluids. The wealth of this literature is such that, in this paper, we restrict ourselves to a cursory discussion of effects directly relevant to the present study. The reader interested in the current state of the art is relegated to [START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF]; [START_REF] Godeferd | Structure and dynamics of rotating turbulence: A review of recent experimental and numerical results[END_REF].

The focus of the interest in rotating wave turbulence [START_REF] Galtier | Weak inertial-wave turbulence theory[END_REF][START_REF] Bellet | Wave turbulence in rapidly rotating flows[END_REF][START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF][START_REF] Godeferd | Structure and dynamics of rotating turbulence: A review of recent experimental and numerical results[END_REF] lays at the anisotropy of scales along the directions parallel and perpendicular to the axis of rotation, the presence of direct and inverse cascades of the key dynamically important quantities (energy, enstrophy, etc.), non-linear wave-wave interactions among inertial waves, and the development of coherent vortex structures aligned with the axis of rotation. The importance of these issues has been identified in early experimental studies with grid-generated turbulence in rotating tanks and their numerical counterparts (see for example [START_REF] Hopfinger | Turbulence and waves in a rotating tank[END_REF]; [START_REF] Godeferd | Direct numerical simulations of turbulence with confinement and rotation[END_REF]).

In order to study in isolation the effect of rotation on (initially isotropic) turbulence, considerable attention has been focused on theoretical investigations in domains of infinite extent and numerical simulations in triply periodic boxes (e.g. [START_REF] Waleffe | Inertial transfers in the helical decomposition[END_REF]; [START_REF] Cambon | Energy transfer in rotating turbulence[END_REF]). It has been shown that the anisotropy develops due to non-linear wavewave interactions modified by rotation and concentrates energy in the plane normal to the rotation axis at a slow two-dimensional manifold [START_REF] Cambon | Energy transfer in rotating turbulence[END_REF]. The relevance of the wave turbulence formalism and results of numerical simulations in triply periodic boxes to the experimental reality involving secondary currents, wall-induced vorticity, and formation of Ekman and Stewartson boundary layers remains an open issue. Indeed, various saturated turbulent regimes ranging between quasi-two-dimensional and wave turbulence can be obtained in numerical simulations in triply periodic domains depending on the damping mechanism imposed onto the geostrophic component to mimic the interaction with rigid boundaries (see, e.g. Le Reun et al. ( 2017)). Therefore, the experimental and numerical investigation of saturated turbulence regimes in a confined volume of rotating fluid attracts significant interest. For such studies, the choice of the range of parameters and of the geometric setup remains a non-trivial issue (see e.g. [START_REF] Godeferd | Structure and dynamics of rotating turbulence: A review of recent experimental and numerical results[END_REF]). Typically, to ensure the development of a fully non-linear energy cascade one needs to ensure a low value of the Ekman number E = ν/(2ΩL 2 ), where ν is kinematic viscosity and L is the global length scale, which for the bounded fluid has the meaning of the container size. Further, the effect of rotation must be sufficiently strong and therefore the global Rossby number Ro L = U/(2ΩL) (here U is a velocity scale) must be sufficiently low. However, Ro L cannot be vanishingly small since it is responsible for triggering non-linear effects. The micro-dynamics of the emerging vortex structures can be conveniently quantified by the micro-Rossby number Ro ξz = ξ z /(2Ω), where ξ z is the perturbative vertical vorticity measured in the rotating frame. The skewness of the Probability Density Functions (PDFs) of micro-Rossby numbers is known to reflect the symmetry breaking of cyclonic/anticyclonic motions, which is a well-known property of rotating systems [START_REF] Bradshaw | The analogy between streamline curvature and buoyancy in turbulent shear flow[END_REF][START_REF] Pedley | On the stability of viscous flow in a rapidly rotating pipe[END_REF][START_REF] Hopfinger | Vortices in rotating fluids[END_REF].

The goal of our joint experimental and numerical investigation is to study the nonlinear fate of an inertial wave attractor in an axisymmetric setting under experimental conditions which are compatible with the formation of a slow two-dimensional manifold coupled to the genuinely three-dimensional inertial wave field. This is achieved by designing a setup in the form of a rotating annulus having a trapezoidal (vertical) cross-section, thereby admitting a wave attractor structure in meridional planes and coherent vortex structures in equatorial plane. Note that this geometry is similar to the one considered in [START_REF] Klein | Inertial wave excitation and focusing in a liquid bounded by a frustum and a cylinder[END_REF] which describes in great detail the formation of wave attractors and the role of the Stokes and Ekman boundary layers in an annular rotating domain with central frustum under forcing due to libration. In the present paper we chose a configuration with a sloping (conical) bottom since, potentially, it may give rise to a topographic β-effect, consisting in the support of topographic Rossby waves [START_REF] Rossby | Relation between variations in the intensity of the zonal circulation of the atmosphere and the displacements of the semi-permanent centers of action[END_REF].

This work has been preceded by a preliminary numerical simulation of [START_REF] Sibgatullin | Direct numerical simulation of three-dimensional inertial wave attractors[END_REF] performed with the help of the spectral element method. The preliminary simulation has been run at a relatively weak forcing and qualitatively demonstrated the loss of axial symmetry and gradual build-up of inertial wave turbulence, but the saturated state, corresponding to the formation of a slow two-dimensional manifold, has not been reached. In the present paper, we address precisely this intriguing and previously overlooked issue. This manuscript is organised as follows. The experimental apparatus is described in section 2. In section 3 we describe the reference experiment demonstrating the emergence of a slow manifold in the form of a polygonal pattern of cyclonic vortices in the equatorial plane, which co-exists with a wave attractor (superimposed with wave turbulence) in the meridional plane. In section 4 we describe the results of direct numerical simulations performed with the help of the spectral element method for a set of forcing amplitudes. By using post-processing diagnostic tools we identify the key dynamic regimes, ranging from linear behaviour to Triadic Resonance Instability (TRI), and the generation of a slow manifold in the form of a Polygonal Vortex Cluster (PVC). The latter regime reproduces the essential features of the reference experiment described in section 3. In section 5 we summarise the key findings of the paper.

Experimental Apparatus

Figure 1 presents a schematic of the experimental apparatus in the vertical and horizontal cross-sections. The region of interest is bounded by two vertical acrylic cylinders, by the wave generator at the top, and by an acrylic conical surface at the bottom. The outer and inner radii of the domain are R 1 = 20.2 cm and R 0 = 5.0 cm, respectively. The generatrix of the conical bottom surface has an inclination of 45 • , and the apex of the cone points upwards. Note that this cone can also be reversed upon needs: the choice taken in the present study is explained in Section 3.

In a vertical (meridional) cross-section of the setup, two trapezoidal domains are facing each other as shown in figure 1: note that the amplitude of the wave generator is greatly exagerated and the upper bound of the fluid domain is actually nearly flat. The depth of fluid measured along the generatrix of the outer cylinder is H = 40 cm. In a horizontal (equatorial) cross-section, the experimental domain is a ring of width L = R 1 -R 0 . The whole setup is inserted into a square acrylic tank of 100 cm × 100 cm horizontal section and 65 cm height, as used in [START_REF] Boury | Excitation and resonant enhancement of axisymmetric internal wave modes[END_REF]. Each part of the facility is rigidly fixed to prevent any parasitic vibration when the whole setup, mounted on the rotating table, is brought to rotation at angular velocity Ω = 2π/T where T is the rotation period. The axis of rotation of the table coincides with the symmetry axis of the setup.

The axisymmetric wave generator [START_REF] Maurer | An axisymmetric inertia-gravity wave generator[END_REF], adapted from a previous planar version [START_REF] Gostiaux | A novel internal waves generator[END_REF], is used to produce inertial waves via a prescribed motion of ring-shaped elements discretizing the annular upper bound of the fluid domain. This device has been slightly modified to fit our needs, by lowering down the cylinders with 20 cm long aluminium rods. In the configuration presented in figure 2, the five inner cylinders have been removed. The motion amplitudes of the remaining eleven outer cylinders (grey boxes in the cross-section in figure 1(a)) have been adjusted to preserve the volume of fluid displaced during its motion, such that the profile z(r) of the generator (2.1)

In order to preserve the boundary condition of non pumping fluid at the cylinder edges, the radial velocity v r has to be zero at R 0 and R 1 . This condition writes, in terms of the profile z(r), as dz dr (r = R 0 ) = dz dr (r = R 1 ) = 0.

(2.2)

As shown in [START_REF] Boury | Excitation and resonant enhancement of axisymmetric internal wave modes[END_REF], this facility is efficient in producing modes 1 to 3 Bessel shaped profiles, though the discretization of the wave generator leads to lower resolved modes at high order. We therefore looked for the closest approximation of a radial mode 1 profile in such a confined geometry. The selected profile is a cubic shaped profile, as shown by the dashed line in figure 2, that sets the cam motion amplitudes. The highest amplitude for a cam is a = 2.5 mm next to the inner cylinder, low enough to ensure a gradual growth of non-linear effects.

The commonly used Particle Image Velocimetry (PIV) technique was implemented to visualise the velocity fields. Horizontal and vertical laser planes were created using a 2 W Ti:Sapphire laser (wavelength 532 nm) and a cylindrical lens. While filling the tank, hollow glass spheres and/or silver coated spheres of 10 µm diameter were added for the purpose of visualisation. Particle displacements were recorded at 40 Hz using a camera located either on the side of the tank (vertical-plane visualisation) or down below facing a 45 • mirror placed under the tank (horizontal-plane visualisation). The CIVx algorithm was subsequently used to process the PIV raw images and extract the velocity fields [START_REF] Fincham | Advanced optimization of correlation imaging velocimetry algorithms[END_REF]. An additional filtering of the wave fields is often performed, in order to isolate the component at a given discrete frequency peak ω = ω i , e.g. at the forcing frequency ω = ω 0 . This post-processing is performed as follows: first, we compute the Fourier transform of the temporal signal associated to the wave field in each space point; the requested frequency is then isolated using a narrow Hamming window; finally, we compute the inverse Fourier transform, giving us the filtered wave field around the frequency of interest. The slow motions are isolated using a low-pass filter, with a cut-off frequency ω c defining the width of the filter. It is important to note that the parameters of the frequency filters applied to the experimental and numerical data are identical. Below, we explicitly discuss the role of the width of the filters when appropriate. This issue is particularly relevant to the analysis of the numerical data.

Experimental Results

Before getting to the description of the experimental results, let us make a few notes on the geometry of the setup and the choice of the parameter range. As discussed in the introduction, the setup is designed to allow for the non-linear coupling between [i] inertial wave attractors in vertical (meridional) cross-sections and [ii] slow two-dimensional manifold in horizontal (equatorial) cross-section. We chose the geometric configuration of the conical bottom with the apex pointing upwards because, in such a geometry, the inertial waves undergo an additional focusing due to geometrical convergence of waves propagating from larger to smaller radial coordinate besides a primary focusing due to the reflection on the cone. This additional focusing favours the onset of instability close to the inner cylinder, as can be seen in the preliminary numerical study of [START_REF] Sibgatullin | Direct numerical simulation of three-dimensional inertial wave attractors[END_REF]. Furthermore, in the saturated regime we expect to localise the slow two-dimensional manifold in the vicinity of the inner cylinder, thereby facilitating the observation of a "vortex condensate" on top of an inertial wave background.

To reduce the effect of viscosity on the non-linear energy transfer we chose a rather high rotation rate of the setup, Ω = 2.093 rad • s -1 , so that the relevant value of the Ekman number

E = ν/(2ΩL 2 ) in the present experiments is reasonably low -considering ν = 10 -6 m 2 • s -1 , Ω = 2.1 rad • s -1 , L = (R 0 -R 1 ) = 0.15 m or 0.4 m (max depth H), leads to E = 1.1 • 10 -5 or E = 1.5 • 10 -6
. The Ekman dynamics typically relates to axial vertical processes, yielding a preference for the latter value. The global a priori Rossby number based on the horizontal scale of the fluid domain, Ro L = U/(2ΩL), can be defined using the maximum vertical speed of the generator rings as the velocity scale so that U = aω 0 , where ω 0 is the forcing frequency. For our experimental conditions, the amplitude a = 2.5 mm and frequency ω 0 = 1.7 rad • s -1 yield Ro L = 7 • 10 -3 so that the system is expected to be strongly affected by the Coriolis force. This low value of the global Rossby number, however, corresponds to a developed non-linear energy transfer so that the observed saturated regime is more appropriately characterized by a relevant micro-Rossby number based on vertical component of vorticity as discussed below.

Moreover, in this experimental section, we discuss only a single experiment that has been run at these parameters (a = 2.5 mm and ω 0 = 1.7 rad • s -1 ). The reasons for that are multiple. First, the range of accessible parameters is very narrow given our experimental apparatus: a ray tracing theory indicates that the forcing frequency should be between 1.47 and 2.14 rad • s -1 in order to observe a wave attractor, and the amplitude set for the generator profile is the optimum we can do. Indeed, on the one hand, since the cam amplitudes are set with a given error, the exact volume-conserving form of the forcing becomes difficult to prescribe as amplitude decreases: the local amplitudes of several cams may fall within the error margin. On the other hand, at high forcing amplitudes the non-linearity sets in too fast to allow the observation of the transient behaviour. At the instrumental level, additional limitations are imposed, especially in the horizontal plane, due to deposition of visualisation particles at the bottom of the test tank which decreases the imaging quality for visualisations in the horizontal plane. We therefore describe in detail only the experiment for which we have a complete data set. Nonetheless, we checked that the phenomenon is reproducible by performing experiments at different frequencies (ω 0 = 1.81 rad • s -1 and ω 0 = 1.92 rad • s -1 ), at a higher rotation rate (Ω = 2.723 rad • s -1 with ω 0 = 2.21 rad • s -1 , to keep the same ratio ω 0 /f ), in addition to many experiments at the same frequency and rotation rate as the one discussed and analysed below. In all these experiments, the same behaviour is qualitatively observed. To propose a more detailed study of the phenomenon, we then proceed with Direct Numerical Simulations (DNS), as presented in section 4, notably by investigating lower forcing amplitudes.

Linear Regime

In the linear regime we recover the classical dynamics: at the time-scale of order 10T 0 after the start of the forcing, where T 0 = 2π/ω 0 is the forcing period, iterative focusing downscales the wave motion from the global scale L to the scale associated with the width of the wave beams [START_REF] Rieutord | Inertial waves in a rotating spherical shell: attractors and asymptotic spectrum[END_REF][START_REF] Hazewinkel | Observations on the wavenumber spectrum and evolution of an internal wave attractor[END_REF][START_REF] Grisouard | Numerical simulation of a two-dimensional internal wave attractor[END_REF]. Typical wave patterns observed in the quasi-linear regime at t = 17T 0 in horizontal and vertical planes are presented in the upper rows of figures 3 and 4 in terms of the quantities filtered at ω = ω 0 and ω = 0, respectively.

Note that the width of the low-pass filter is set by the upper cut-off frequency, i.e. ω c = ω 0 /3. As discussed later, the essential low-frequency content of the experimental signal (see figure 8) lies well below the chosen cut-off frequency. The choice of the width of the low-pass filter is important for a unified description of numerical and experimental results (see sections 4.4. and 4.5 for details). For the experimental results as such the choice of the cut-off frequency is less important (i.e. a lower value of ω c could be prescribed). For the sake of brevity, the result of the low-pass filtering is below referred to as the "signal filtered at ω = 0".

For clarity, we visualise the fields of radial v r and azimuthal v θ velocity, and vertical vorticity ξ z in the horizontal (equatorial) plane, and the corresponding field of vertical velocity in the vertical (meridional) plane. It can be seen that the wave pattern observed at the forcing frequency (figure 3) in the horizontal plane is to a good approximation axisymmetric, while in the vertical trapezoidal cross-section we recover a classic pattern of the (1, 1) wave attractor [START_REF] Maas | Observations of an internal wave attractor in a confined stably stratified fluid[END_REF][START_REF] Maas | Wave focusing and ensuing mean flow due to symmetry breaking in rotating fluids[END_REF] in agreement with the ray tracing, whose branch width is due to an equilibrium between wave focusing and viscous dissipation [START_REF] Rieutord | Inertial waves in a rotating spherical shell: attractors and asymptotic spectrum[END_REF][START_REF] Hazewinkel | Observations on the wavenumber spectrum and evolution of an internal wave attractor[END_REF][START_REF] Grisouard | Numerical simulation of a two-dimensional internal wave attractor[END_REF]). The experimental signal filtered around ω = 0 remains weak at t = 17T 0 (see figure 4).

Non-linear Regime

The development of the fully saturated non-linear regime is illustrated in figures 3 and 4 by snapshots corresponding to t = 50T 0 and t = 150T 0 . The full vortex pattern representing the slow two-dimensional manifold is formed at the time scale of 100T 0 . It can be clearly seen that the initial axisymmetry observed at t = 17T 0 is lost while the slow manifold is gradually formed. The latter is represented by a regular polygonal system of eight cyclonic vortices. The vortices are nearly invariant in the vertical direction as attested by the right column of images in figure 4, representing the vertical velocity component. These vortex structures are reminiscent of the Taylor columns usually found in rotating systems but there is, however, a crucially important distinction: while the Taylor columns are normally formed as a consequence of a slow motion of a perturbation imposed on the rotating fluid, the coherent structures seen in figure 4 arise due to a non-linear process which drains energy from the wave field toward the slow twodimensional manifold. We note that the vertical velocity in the cyclonic vortices is directed downwards, corresponding to Ekman pumping, in agreement with existing experimental and numerical data (e.g. [START_REF] Hopfinger | Turbulence and waves in a rotating tank[END_REF]; [START_REF] Godeferd | Direct numerical simulations of turbulence with confinement and rotation[END_REF]). Interestingly, [START_REF] Lopez | Instability and mode interactions in a differentially driven rotating cylinder[END_REF] observed a similar breaking of a flow, axisymmetrically forced by a counter-rotating surface lid introducing azimuthal shear, into a discrete pattern of vortices distributed with an azimuthal periodicity. In their experiment, they relate the formation of such a manifold to the linear instability of the shear layer producing an azimuthal wave-like structure and, although much of the dynamics following the development of this instability is similar to the one we describe -e.g. the formation of a polygonal pattern of cyclonic vortices, intensified due to vortex stretching-, the mechanism is different. In our setup, the forcing introduces vertical motion in the fluid, and very little azimuthal motion is present at the beginning. As we will show in the numerical section, the instability observed in the non-linear regime (figures 3 and 4) is then likely due to a triggered Triadic Resonant Instability that grows and becomes itself unstable. An internal boundary layer is formed around the attractor, which is prime candidate for the mixing of angular momentum, resulting in a sheared mean flow and it is worth mentioning that even in the fully saturated regime one can still identify the branches of the inertial wave attractor in the signal filtered around the forcing frequency ω = ω 0 (see figure 3). The relevance of such experimental regime (where a wave attractor in the meridional plane co-exists with a polygonal vortex system in the equatorial plane) to geo-and astrophysical systems admitting the existence of inertial-wave attractors [START_REF] Dintrans | Gravito-inertial waves in a rotating stratified sphere or spherical shell[END_REF][START_REF] Rieutord | Wave attractors in rotating fluids: a paradigm for ill-posed cauchy problems[END_REF][START_REF] Rieutord | Viscous dissipation by tidally forced inertial modes in a rotating spherical shell[END_REF] represents an interesting direction for future research.

The visual evidence of the vortices is seen in figure 5 which presents the temporal evolution of the radial structure measured by sampling the azimuthal distribution of radial velocity at the radius r = 8 cm corresponding to the position of centres of vortices seen in figure 4. Figure 5 shows that coherent structures start to appear after roughly 25T 0 , and that further the vortex pattern self-organises itself so that new structures gradually appear and join the ensemble. After roughly 100T 0 all structures move at the same rate in cyclonic direction. This rate corresponds to half a turn of the vortex cluster around the inner cylinder which lasts about (88 ± 2)T 0 .

The vertical vorticity field in the horizontal plane can be characterized statistically, by measuring the Probability Density Function (PDF) of the micro-Rossby number Ro ξz = ξ z /(2Ω). The typical PDFs of Ro ξz corresponding to different stages of development of the coherent vortex structures are shown in figure 6. We calculate the PDFs for the raw signal and for the signal filtered around the forcing ω = ω 0 and ω = 0 frequencies. The PDFs are calculated over the surface of the ring-shaped zone between the inner and outer cylinders, and over the time-span of ±2T 0 around the time instances indicated in figure 6. It can be seen that at the beginning of the process, when the motion is represented essentially by the axisymmetric waves, the PDFs of Ro ξz have a sharp symmetrical form. As the non-linear energy transfer towards coherent vortex structures develops, there is a progressive evolution of the vorticity PDFs toward the shape characterized by asymmetric "shoulders", which indicates a well-pronounced cyclonic/anticyclonic asymmetry. This asymmetry is clearly seen in the PDFs calculated over the raw signal and the signal filtered around zero frequency ω = 0 (upper and lower images in figure 6), suggesting that a few strong cyclonic vortices seen in figure 4 are responsible for the asymmetry of the PDFs. The probability density functions calculated for the signal filtered at the forcing frequency (middle image in figure 6) remain approximately symmetrical at any time. A slight asymmetry visible in the curve corresponding to time around t = 30T 0 and 40T 0 can be tentatively attributed to the process of genesis of the regular Polygonal Vortex Pattern: new cyclonic vortices are emerging in the plane of visualization and are joining the ensemble. Nonetheless, since the asymmetry is weak one cannot exclude also a minor contribution from experimental noise.

Numerical Results

Numerical Procedure

The experimental section of this paper provides an important reference result: we explicitly show that under sufficiently strong forcing the motion in the system represents a non-linear combination of a complex inertial wave field, with an identifiable wave attractor at the forcing frequency, and a slow quasi-two-dimensional manifold represented by a Polygonal Vortex Cluster. In the numerical simulations described below, we faithfully reproduce the geometry of the setup shown in figure 1, and we use the experimental value of the rate of background rotation Ω = 2.093 rad • s -1 and the forcing frequency ω 0 = 1.70 rad • s -1 . We explore the effect of the forcing amplitude a on the long-term non-linear evolution of the inertial wave field in the setup. Our purpose is to identify a sequence of observed regimes, ranging from the linear to the non-linear regimes with a Polygonal Vortex Cluster described in the previous section.

The mathematical formulation of the problem consists in the Navier-Stokes equations and the continuity equation

∂v ∂t + (v • ∇) v = -∇p + ν∆v + 2 Ω × v, (4.1) p = p ρ - 1 2 |Ω × r | 2 , (4.2) ∇ • v = 0. (4.3)
The governing equations are written in the Cartesian system (x, y, z), co-rotating with the setup, where the z-axis points upwards and coincides with the axis of rotation. The origin of the coordinate system is taken at the center of the upper lid. Note that the fixed lid has no inherent rotation in the rotating system (x, y, z). In the fixed inertial nonrotating reference frame we define the anti-clockwise background rotation Ω as positive.

In the rotating reference frame the cyclonic vorticity is positive, and the sense of the prograde motion is positive -cyclonic is by definition in the direction of rotation; here we take a right-handed Cartesian coordinate system such that in cylindrical coordinates, (r, θ, z), the azimuthal coordinate θ increases in cyclonic direction. With this convention the cyclonic vortices appear in red colour both in numerical and experimental vertical vorticity patterns. If the motion in a vertical radial plane is considered (e.g. plane y = 0) we use the notations (r, z). Throughout the paper we visualize the fluid motion in the horizontal plane located at mid-depth of the fluid volume.

To complete the mathematical statement of the problem we prescribe the boundary conditions as follows: the no-slip condition is imposed at all rigid boundaries except at the fixed flat horizontal upper lid where a specific harmonic forcing is applied to simulate the experimental one. We thus require that both horizontal components of the fluid velocity vector at the upper lid are equal to zero as in standard no-slip condition, while the vertical component is prescribed as explained below. Figure 2 shows the discrete experimental profile of the amplitude of motion of the rings in the generator. This discrete stepwise form approximates the smooth profile z(r) shown in figure 2 by the dashed blue line. The forcing imposed in the numerical experiments at the upper lid is axisymmetric, with the vertical component of the fluid velocity given by v z (r, t) = az(r)ω 0 exp(ω 0 t) where a is the forcing amplitude and z(r) is a non dimensional profile (of unit amplitude at r = R 0 ) such that z(r) = az(r) corresponds to the dashed blue curve of figure 2. Since we consider a small-amplitude input perturbation (more precisely, a/ (R 1 -R 0 ) 0.0167), such approach seems to be justified. In the case of internal wave attractors with similar implementation of the input forcing, we observed a good qualitative agreement between numerical and experimental results (Brouzet et al. 2016b). However, an extension to rotating fluids stratified in angular momentum, and not in density, and driven by Coriolis force (which does no work), and not by gravity (which does), is not fully evident. Both this issue and the difficulty of precise evaluation of the efficiency of the experimental wave generator leaves a considerable margin of uncertainty regarding the possible correspondance between the experimental amplitude and its numerical counterpart.

The numerical simulations have been performed for the following set of 7 forcing amplitudes: a = 0.2; 0.5; 1.0; 1.8; 2.0; 2.4; 2.5 mm, the latter value corresponding to the experimental case. This corresponds to Rossby numbers Ro L ranging from 5 • 10 -4 to 7 • 10 -3 . Below, we use a reduced representative set of amplitudes to describe the key regimes observed. The typical duration of the numerical experiments was about 200T 0 , and in some cases up to 350T 0 (e.g. at a = 2.4 mm).

As evident from the experimental part of the paper, we need to model a strongly non-linear dynamical problem, where both viscosity and non-linearity play a role. The numerical simulation of transient and turbulent regimes is a challenge as we have to follow the development of small-scale structures during long time intervals. In this context, spectral or Galerkin decomposition is known to be a robust approach to tackle the non-linear effects without parasitic effects due to numerical viscosity. Classically, such a decomposition is possible only for simple geometry and boundary conditions. In the present work the direct numerical simulations are performed with the help of the spectral element method, using the open source code Nek5000 (see [START_REF] Fischer | Spectral element methods for large scale parallel navierstokes calculations[END_REF]; [START_REF] Fischer | An overlapping schwarz method for spectral element solution of the incompressible navierstokes equations[END_REF]; [START_REF] Deville | High-Order Methods for Incompressible Fluid Flows[END_REF]). This method combines the advantages of highorder decomposition with geometric flexibility, and permits to run long-term simulations of strongly non-linear dynamics. In the present study, we have used meshes with up to 100 thousand elements, with eighth-order polynomial decomposition within each element (up to 50 millions DoF).

Energy Spectra: a Preliminary Classification of the Observed Regimes

The snapshots of the simulated inertial wave fields in horizontal and vertical planes are shown in figure 7 for a set of forcing amplitudes a = 0.2; 1.0; 2.0; and 2.4 mm at t = 100T 0 (all images correspond to the same phase of the forcing). As the amplitude increases, one observes the increasing complexity of the inertial wave fields, in broad agreement with the known literature on the onset of Triadic Resonance Instability (Bordes et al. 2012;[START_REF] Bourget | Experimental study of parametric subharmonic instability for internal plane waves[END_REF]) and wave turbulence in internal and inertial wave attractors (Brouzet et al. 2016a,b). In particular, one can see [i] visible broadening of the attractor branches at higher forcing (Brouzet et al. 2017b), [ii] signature of discrete azimuthal symmetry in the vertical vorticity pattern in horizontal plane at a = 1.0 mm in agreement with preliminary simulations described in Sibgatullin et al. (2017), and [iii] emerging signature of discrete patches of vertical vorticity in the horizontal plane at a = 2.4 mm which are reminiscent of the experimental observations described in section 3. Below we apply a set of post-processing tools to analyse the numerical data in some detail.

The non-linear regimes observed in the numerical simulations of an axisymmetrically forced inertial wave attractor can be roughly classified by considering the development of the signal spectra with time for a set of forcing amplitudes. The typical time-frequency diagrams obtained at different values of a are shown in figure 8. These diagrams are calculated similarly to [START_REF] Bourget | Experimental study of parametric subharmonic instability for internal plane waves[END_REF], as follows

S r (ω, t) = +∞ -∞ v r (r, θ, τ )e iωτ h(t -τ ) dτ 2 rθ , (4.4)
where h is a Hamming window and v r is the radial component of the velocity field, and subscript rθ denotes a ring-shaped domain around the inner cylinder (θ ∈ [0, 2π], r ∈ [5, 10] cm). The calculations are performed with the Matlab toolbox described in [START_REF] Flandrin | Time-Frequency/Time-Scale Analysis, Time-Frequency Toolbox for Matlab[END_REF]. The time-frequency diagrams were also calculated for other variables (e.g. azimuthal velocity, vertical vorticity), demonstrating similar qualitative behaviour.

The signatures of inertial wave turbulence can be detected by analysing the energy spectra [START_REF] Yarom | Experimental observation of steady inertial wave turbulence in deep rotating flows[END_REF]. Note that the original approach implying fourdimensional Fourier analysis (three dimensions for space and one for time) has been proposed by [START_REF] Yarom | Experimental observation of steady inertial wave turbulence in deep rotating flows[END_REF] for a rotating system where the energy is injected via decorrelated random forcing thereby creating well-developed fully three-dimensional wave turbulence possessing no information on initial orientation of the input wave vector. In the present paper the energy is injected via axisymmetric deterministic forcing (with prescribed amplitude, length and frequency) producing a variety of dynamic regimes ranging from regular to turbulent (see figure 7). In order to adapt the method of postprocessing to the observed wave patterns, we perform the analysis described in [START_REF] Yarom | Experimental observation of steady inertial wave turbulence in deep rotating flows[END_REF]) in a vertical radial plane similar to Brouzet et al. (2016a), and, in addition, calculate the spatial energy spectrum in a horizontal plane. The latter is done to clarify the role of the integer-number azimuthal modes at different values of the forcing amplitude. Note that owing to different design of the experimental setup such modes were not present in the case studied by [START_REF] Yarom | Experimental observation of steady inertial wave turbulence in deep rotating flows[END_REF]. Before computing the energy spectra, the numerical results are re-sampled to yield the spatial resolution of 0.2 cm × 0.2 cm and the temporal resolution of 0.2 s. This is done to reduce the amount of data and to match the typical experimental resolution.

The wave energy spectra in the vertical radial plane are computed using the three dimensional Fourier transforms of the horizontal and vertical velocity fields, vr (k r , k z , ω) and vz (k r , k z , ω). The corresponding energy spectrum is defined as

E (k r , k z , ω) = |v r (k r , k z , ω)| 2 + |v z (k r , k z , ω)| 2 2A T , (4.5)
where A = 40 × 15 cm 2 is the area considered in the vertical plane and T = 50T 0 is the duration of the time-history sample. The calculations for all the cases shown below in figure 9 are performed for the numerical data obtained for t ∈ [100T 0 , 150T 0 ]. The spatial resolution of re-sampled data and the size of the fluid domain provide respectively the upper and lower bounds in wave numbers, of order k max = 8 rad • cm -1 and k min = 0.2 rad • cm -1 . Interpolation is performed to express the energy spectrum E (k r , k z , ω) as a function of E (k, α, ω), with k, the norm of the wave vector. Then, we perform integration over the entire range of resolved wave vectors [k min , k max ] as follows

E (α, ω) = kmax kmin E (k, α, ω)kdk. (4.6)
The calculated energy density E (α, ω) is normalised by the frequency energy density E (ω), obtained by integrating E (α, ω) over all directions. It should be stressed that it is only for a purely axisymmetric wave perturbation that E (α, ω), defined by (4.6), is a function of the true angle α between the vector of the phase speed and the horizontal plane. In the axisymmetric case, the vector (k r , k z , k θ ) has k θ = 0, where the direction of the component k θ is defined by a vector orthonormal to the (r, z) plane. In the general case, the vector (k r , k z , k θ ) is inclined at an angle α = arctan k z / k 2 r + k 2 θ 1/2 to the vertical axis. A projection of this vector onto the vertical (r, z) plane is seen at the apparent angle α * = arctan [k z /k r ]. Restricting, for brevity, our attention to small values of α and α * , the difference between the two angles can be quantitatively characterized as 1 + (k 2 θ /k 2 r ) 1/2 -1. For example, the relative difference |(α *α)/α| < 0.1 translates to |k θ /k r | < 0.46, thereby admitting the vectors whose azimuthal directions differ from the (r, z) plane by less than ±25 • . An additional analysis is performed in the horizontal plane located at the mid-depth of the setup. The wave energy spectrum in the horizontal plane is computed using the three dimensional Fourier transform of the horizontal components of the velocity field, vr (k r , n θ , ω) and vθ (k r , n θ , ω). Here n θ = 2π/∆θ is the non-dimensional azimuthal wave number, where ∆θ is the azimuthal wavelength (in radians). Note that using n θ is more convenient than using k θ (measured in rad • m -1 ) in view of the discrete azimuthal symmetry of the wave patterns in the horizontal plane clearly seen in figure (7) at a = 1.0 mm. The corresponding energy spectrum is defined as

E ⊥ (k r , n θ , ω) = |v r (k r , n θ , ω)| 2 + |v θ (k r , n θ , ω)| 2 2A r,θ T , (4.7)
where A r,θ is the area considered (r ∈ [R 0 , R 1 ], θ ∈ [0, 2π]) and T = 50T 0 is the duration of the time-history sample. The spectrum is integrated over the frequency range ω ∈ [-f, f ] in which inertial waves can propagate, as follows

E ⊥ (k r , n θ ) = f -f E ⊥ (k r , n θ , ω)dω. (4.8)
Using (4.4), (4.6) and (4.8), we obtain figures 8 and 9, from which we can infer the information on qualitative evolution of the wave regime with amplitude. It can be seen in figure 8 that for the lowest amplitude considered (a = 0.2 mm), the frequency spectrum is monochromatic, with a weak but detectable component at twice the forcing frequency. This is confirmed by the energy spectrum for this amplitude represented in figure 9: the energy is localized on the linear dispersion relation at discrete "spots" corresponding to ω 0 /f = 0.39 and 2ω 0 /f = 0.78. The perturbation is nearly axisymmetric, a weak contribution of low azimuthal modes can be detected in terms of E ⊥ (k r , n θ ).

At higher forcing amplitude (a = 1.0 mm), one can observe (see figure 8) a classic signature of triadic resonance instability (TRI). One can see signal components at the frequencies of two most energetic secondary waves which satisfy the condition ω 0 = ω 1 + ω 2 , with ω 1 /ω 0 = 0.32 and ω 2 /ω 0 = 0.68 (here the subscripts 0, 1 and 2 denote the parameters of the primary wave, and two secondary waves). Also there is a significant energy content at frequencies ω 1 + ω 0 , ω 2 + ω 0 , and the multiples of the The white dashed lines correspond to the dispersion relation ± cos α = ω/f . Note that the energy peaks are not always localized at the dispersion curve, owing partially to the departure of the wave vectors from (r, z) plane and to the contribution of slaved (evanescent) waves (see discussion in the text). The lower four images show the energy spectra in the horizontal plane E ⊥ (kr, n θ ) on a logarithmic scale at different values of the forcing amplitude. By combining information from the spectra calculated in the vertical and horizontal planes one can see the evolution of wave regimes from nearly linear at a = 0.2 mm to wave turbulence evolving from a "discrete" form at a = 1.0 to a more "continuous" form at 2.0 and a = 2.4 mm. The latter is particularly well seen in terms of E ⊥ (kr, n θ ). In all cases the axisymmetric component of the wave field (corresponding to n θ = 0) is most significant.

forcing frequency 2ω 0 and 3ω 0 . Thus, both TRI and the Two-Wave Interactions (TWI) [START_REF] Beckebanze | Experimental evidence of internal wave attractor signatures hidden in large-amplitude multi-frequency wave fields[END_REF][START_REF] Boury | Experimental generation of axisymmetric internal wave super-harmonics[END_REF])) are present in the system and play a role. The corresponding energy spectrum shown in figure 9 demonstrates that the linear dispersion relation attracts the energy maxima. However, some energy content falls apart from the dispersion relation, producing a pattern of weak horizontal stripes as consequence of the difference between angles α and α * arising from quasi-two-dimensional analysis (4.6). In terms of E ⊥ (k r , n θ ) we observe a contribution of low azimuthal modes, where significant peaks can be distinguished at n θ = ±12. As it will be discussed in section 4.3 the latter is related with spiral waves propagating in retrograde and prograde sense which are likely to be a consequence of triadic resonance. The data analysis presented in figures ( 8) and ( 9) for a = 1 mm supports the idea that this regime corresponds to the onset of weak (discrete) inertial wave turbulence.

As the forcing amplitude increases further (a = 2.0, 2.4, and 2.5 mm), one can observe a significant increase of the energy content in the continuous part of the frequency spectrum, accompanied with an enrichment of the discrete part (figure 8). The peaks corresponding to TRI-generated secondary waves ω 1 and ω 2 (detected at a = 1.0 mm) remain persistent at higher forcing amplitudes. In general, the overall trends seen in figure 8 at higher forcing amplitudes correspond well to the effects described in the literature on internal wave attractors in stratified fluids in the non-linear regime (Brouzet et al. 2016a(Brouzet et al. , 2017a)). There is, however, an important distinction: in the case of rotating fluid we observe a highly complicated frequency content of the signal at the low-frequency end of the spectrum which does not have a straightforward interpretation in terms of TRI. In the context of this paper this issue is of central interest (see section 4.5 for discussion). Note that some discrete frequency components seen in the numerical calculations performed at a = 2.4 and 2.5 mm appear also in the experimental spectrum. In the numerical simulations, the transients seem to have a significantly longer duration as compared to the experiment. This issue is discussed in section 4.5.

The energy spectra for higher forcing amplitudes shown in figure 9 seem to be consistent with the concept of inertial wave turbulence. As a increases, the character of wave turbulence gradually evolves from a "discrete" form at a = 1.0 to a more "continuous" form at a = 2.0 and 2.4 mm. The latter is particularly well seen in terms of E ⊥ (k r , n θ ). The presence of a significant continuous component in E ⊥ (k r , n θ ) is a direct consequence of the loss of the discrete azimuthal symmetry observed in figure 7 at a = 2.0 mm and 2.4 mm.

As the energy spectra E (α, ω) are calculated in the vertical plane we can see that the structure of the low-frequency zones (where the dispersion relation crosses the line ω/f = 0) evolves considerably as a increases. In the quasi-linear case (a = 0.2 mm) there is no detectable energy component at ω/f 0. As a increases, a complex discrete structure emerges in the vicinity of ω/f 0. At higher a this structure evolves toward a smoothed energy distribution. The latter can be interpreted as a trend toward merging of energy peaks corresponding to discrete frequencies.

It has been already noted that some energy content falls apart from the dispersion relation due to the difference between α and α * as a consequence of the two-dimensional analysis in (4.6). Let us note that there is an additional mechanism for the concentration of such peaks, namely along the vertical lines at fixed α corresponding to the forcing frequency ω 0 . This means that forced higher harmonics are generated which propagate at an angle belonging to that of the fundamental frequency, as opposed to free higher harmonics that follow the dispersion curves. These "slaved" higher harmonics have been encountered in internal wave attractors before [START_REF] Lam | Internal wave focusing revisited; a reanalysis and new theoretical links[END_REF]). An interested reader is relegated to [START_REF] Davis | Attracteurs d'ondes internes de gravité: des résonances en cascade. Une approche expérimentale des régimes linéaire et non-linéaire[END_REF], in which a similar phenomenon is observed and discussed in more detail.

Triadic Resonant Instability in a Rotating Annulus

Relevant to our study, let us consider in more detail the Triadic Resonance Instability (TRI) in the rotating annulus. TRI has been observed in various configurations of rotating flows involving precession, and often led to symmetry breaking and to the formation of vortical structures distributed periodically [START_REF] Albrecht | Triadic resonances in precessing rapidly rotating cylinder flows[END_REF][START_REF] Marques | Precession of a rapidly rotating cylinder flow: traverse through resonance[END_REF][START_REF] Albrecht | On triadic resonance as an instability mechanism in precessing cylinder flow[END_REF][START_REF] Lopez | Rapidly rotating precessing cylinder flows: forced triadic resonances[END_REF]Wu et al. 2020a). For example, in the case of a rotating cylinder (e.g. [START_REF] Marques | Precession of a rapidly rotating cylinder flow: traverse through resonance[END_REF]), this instability gives rise to structures aligned with the rotation axis and distributed along θ.

The data presented in figure 8 for a = 1.0 mm show well-localised discrete frequency components. Using the technique of Hilbert transform filtering introduced in [START_REF] Mercier | Refection and diffraction of internal waves analysed with the hilbert transform[END_REF], one can separate the key components of the wave patterns observed in vertical and horizontal planes as shown in figures 10 and 11. The filtered wave field components seen in the vertical plane remind of the patterns already described in literature for internal wave attractors (see e.g. [START_REF] Scolan | Nonlinear fate of internal waves attractors[END_REF]; Brouzet et al. (2017a)). The length of the wave vector components can be evaluated with reasonable accuracy as location of the maximum of corresponding probability density functions calculated over the domain of interest. The corresponding vector triad in vertical plane is shown in figure 12 (left panel) together with the curves defining the admissible wave triads. It can be seen that the observed triad is consistent with the theoretical curve, with reasonable experimental accuracy.

The pattern of triadic resonance in a ring-shaped domain seen in the horizontal plane posesses a discrete symmetry as shown in figure 11. This is in agreement with earlier observations described in [START_REF] Sibgatullin | Direct numerical simulation of three-dimensional inertial wave attractors[END_REF], and with raw snapshots shown in figure 7 for a = 1.0 mm in the present paper. Thus, the azimuthal pattern in the horizontal plane is reminiscent of 'modal' triadic resonance in a rectangular box described in [START_REF] Mcewan | Degeneration of resonantly excited standing internal gravity waves[END_REF]. However, the modal pattern of internal waves in the vertical plane [START_REF] Mcewan | Degeneration of resonantly excited standing internal gravity waves[END_REF]) is clearly compatible with the dispersion relation, while we are not aware of any theoretical work predicting the number of the expected azimuthal mode of the secondary waves in a ring-shaped domain. Assuming the triadic resonance as the underlying key mechanism and a discrete symmetry of the azimuthal wave pattern, we can expect that the purely axisymmetric wave (i.e. zero azimuthal mode) at the forcing frequency ω 0 should give rise to two secondary waves propagating in the opposite azimuthal directions (cyclonic and anti-cyclonic) and corresponding to the same azimuthal mode. This is confirmed by the data presented in figure 11 and by the construction of the projection of the wave vector triad on the horizontal plane shown in figure 12 (right panel). Both secondary waves seen in the horizontal plane correspond to 12-th azimuthal mode. The secondary waves corresponding to ω 1 and ω 2 propagate in cyclonic and anti-cyclonic directions, respectively. Summing up, we can characterise each wave component by three numbers (k r , k z , n θ ), where k r , k z are conventional wave vector components measured in r, z-plane (in

-5 -4 -3 -2 -1 0 1 2 3 4 5 6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 (a) k 0 k 2 k 1 k r (rad • cm -1 ) k z (rad • cm -1 ) -2 0 2 4 6 8 10 12 14 -2 -1 0 1 2 3 4 5 (b) k 0 k 2 k 1 n θ k r (rad • cm -1 )
Figure 12. Verification of the triadic resonance: (a) curve of admissible wave-vector triads in the vertical plane (kr, kz) with superimposed measured wave vectors; (b) wave-vector triad in the horizontal plane (kr, n θ ). The wave-vector triads shown in the left and right panels represent their projections on the vertical and horizontal planes respectively, and correspond to the patterns depicted in figures 10 and 11. The solid, dashed, and dash-dotted lines in panel (a) represent the locus of the tips of the resonant wave vectors according to the dispersion relation.

rad • cm -1 ), while n θ = 2π/∆θ is non-dimensional integer azimuthal mode number, where ∆θ is the azimuthal wave length (in radians). For the vector triads depicted in figure 12, we have (1.38, 0.94, 0), (-1.88, -0.43, -12) and (3.02, 1.26, 12) for the primary and two secondary waves, respectively. With reasonable accuracy we have

k 0 r ≈ k 1 r + k 2 r , k 0 z ≈ k 1 z + k 2 z , and n 0 θ = 0 = n 1 θ + n 2 θ .
It should be stressed that the wave vector components can be strictly defined only for spatially monochromatic fields (or approximately monochromatic as in [START_REF] Bourget | Experimental study of parametric subharmonic instability for internal plane waves[END_REF]) while for narrow wave beams the objective measurement of wave-vector components raises some problems as discussed in [START_REF] Fan | Finite-amplitude instabilities of thin internal wave beams: experiments and theory[END_REF]. However, experimentally it is often possible to construct the Probability Density Functions (PDFs) for the wave vector components measured in a zone of interest and estimate the length of the wave vectors from the positions of PDF's maxima. This approach is taken in the present paper. In light of the results presented in figure 12, we note that the resonance conditions for TRI are satisfied here. Moreover, it is interesting to point out that the primary wave is three-dimensional but axisymmetric, and therefore mostly lives in a vertical plane, contrary to the two secondary waves, which are non-axisymmetric and fully three-dimensional. At this stage, however, the discussion of this 3D TRI is purely exploratory. The reader interested in a more thorough development is referred to Boury (2020) (section 7.5.2).

Transition to a Polygonal Vortex Cluster

Let us now consider how the regime in which a Polygonal Vortex Cluster emerges in the numerical simulations when we systematically increase the forcing amplitude. By considering the results obtained at a = 1.0 mm we have already made an important observation that the low-frequency behaviour of the system is represented by cyclonic propagation of a secondary wave generated by triadic resonance. In other words, we can say that, after the onset of TRI, this low-frequency component of the wave field represents a "germ" of the slow manifold. Therefore, to compare the flow patterns representing the slow manifold in the horizontal plane it is reasonable to use the low-pass filtering with the cut-off frequency set around ω c = ω 0 /3. The width of this filter corresponds to the width of the filter used in the experimental part of the paper to separate the lowfrequency signal. Moreover, the time-frequency diagrams shown in figure 8 allow us to conclude that this width of the low-pass filter captures the essential features of the lowfrequency behaviour of the system. By applying this filter we can identify the variation of the structure of the vertical vorticity field at different forcing amplitudes, as shown in figure 13. The snapshot of the pattern is complemented by a video in supplementary material. The snapshots and the video are taken in vicinity of the time instant t = 100T 0 . It can be clearly seen that the cyclonic wave motion at a = 1.0 mm is replaced by a highly complicated pattern which exhibits cyclonic/anti-cyclonic rotation close to the inner/outer cylinder, respectively. We illustrate such a regime by the data obtained at a = 2.0 mm. This regime is observed in a certain range of the forcing amplitudes (around a ∈ [1.8 mm; 2.0 mm]). Thus, the transition from the regime with wave field possessing discrete azimuthal symmetry (a = 1.0 mm) to the regime with polygonal vortex cluster (a = 2.4 mm) is highly non-trivial. The intermediate regime observed at a = 2 mm deserves a special study which falls out of scope of the present paper.

As the forcing amplitude increases further, one can clearly identify the patches of cyclonic vorticity arranged in polygonal fashion, with a slow drift of the vortex cluster in prograde direction. This regime is shown in figure 13 for a = 2.4 mm. A very similar regime is observed at a = 2.5 mm, which suggests that the regime is sufficiently robust and can be reproduced in a certain range of the forcing amplitudes. The numerically obtained vortex pattern with 7 vortices arranged at the vertices of a regular polygon can be compared to the experimental vortex pattern with 8 vortices. We see that qualitatively the patterns are similar, and, since the colour scale for both patterns is the same, at quantitative level there is a reasonable agreement between the magnitudes of the vertical vorticity. However, the mechanism of the formation of the vortex patches and the longterm evolution of the pattern in the numerical calculations remain to be identified. To do this we performed a long series of calculations specifically for this regime at a = 2.4 mm. The results of these simulations are described below.

Long-term Fate of the Polygonal Vortex Cluster in Numerical Simulations

The long-term time-frequency diagram corresponding to the case a = 2.4 mm is presented in figure 14. The diagram is calculated over 350 forcing periods. Window 6 shows schematically the domain where the low-pass filter with the cut-off frequency ω c = ω 0 /3 has been applied to obtain the vortex pattern depicted in figure 13 for a = 2.4 mm. This pattern remains virtually the same over a long time span. In particular, the 7-vortex cluster can be observed when window 6 is shifted in time to cover the narrow peaks 3 and 4. These peaks 3 and 4 are well-localized, making possible narrowband filtering of the corresponding wave components of the vertical vorticity field. The result of such filtering is represented in figure 14 by the patterns corresponding to ω 3 = 0.36 rad • s -1 (ω 3 /ω 0 = 0.21) and ω 4 = 0.16 rad • s -1 (ω 4 /ω 0 = 0.09). The former represents a wave of 7-th azimuthal mode propagating in cyclonic (prograde) direction in the close vicinity of the inner cylinder, the positive and negative patches of vorticity in this wave have the same magnitude. The latter represents a wave of 5th azimuthal mode propagating in anti-cyclonic (retrograde) direction. Note that the wave corresponding to ω 4 /ω 0 = 0.09 is significantly weaker in magnitude than the wave corresponding to ω 3 /ω 0 = 0.21. Apart from this well-defined discrete components there is also a contribution from the continuous part of the spectrum. Thus, the 7-vortex pattern seen in figure 13 can be interpreted as azimuthal wave of 7-th mode propagating in prograde sense over weak background vorticity pattern. In this background pattern we can identify 5-th azimuthal retrograde wave mode while other discrete wave contributions are difficult to identify. The observable result of such superposition is the polygonal 7vortex pattern.

It should be noted that at long time scale the observable vortex pattern evolves. It can be seen in figure 14 that after approximately 225T 0 the discrete frequency components at the low-frequency end of the spectrum exhibit a strong trend towards merging. Towards the time around 275T 0 it becomes difficult to distinguish individual low-frequency wave components. The low-pass filtering applied to the domain schematically shown by rectangle 5 returns the vorticity pattern which develops as result of merging of the lowfrequency discrete wave components. The resulting slow manifold is represented by a Polygonal Vortex Cluster with 8 cyclonic vortices drifting in cyclonic direction. This pattern is very similar, in terms of strength and spatial arrangment of vortices, to the experimental pattern shown in figure 13. The evolution of the azimuthal drift of the vortex cluster with time is illustrated in figure 15 representing the numerical counterpart of figure 5. It can be seen that numerical simulations yield 6 to 3 times faster drift than the experimentally measured one.

An important remark should be made regarding the slow background azimuthal current observed in the horizontal plane at mid-depth of the setup. The evolution of this current is shown in figure 16 which represents the time-history of radial distribution of mean azimuthal velocity component V θ θ,t (measured in cm/s). Here • θ,t denotes the azimuthal averaging performed over θ ∈ [0, 2π] and the temporal averaging performed in a moving window of width 20T 0 . The corresponding non-dimensional quantity is introduced as V θ = V θ θ,t • T 0 /(2πr), which physically corresponds to the portion of the full circle passed by the mean azimuthal current during one forcing period. The time-history covers the same span of time as figures 14 and 15. Figure 16 the radial distribution of mean non-dimensional azimuthal velocity corresponding to t = 100T 0 in the numerical simulations performed at different values of the forcing amplitude a = 0.2; 1.0; 2.0; 2.4 mm. This time t = 100T 0 corresponds to the patterns shown in figure 13. The profiles displayed in figure 16(b) show that the mean azimuthal "wind" in the system considerably increases with the forcing amplitude. This effect is quantitatively investigated in figure 16(c), showing that the maximum of |V θ | increases proportionally to the square of the forcing amplitude a, a result consistent with other studies (see, e.g., the recent work of [START_REF] Cebron | Mean zonal flows induced by weak mechanical forcings in rotating spheroids[END_REF]). No relevant scaling could be found, however, for the radial location of these maxima. There exists a significant literature on the effect of the mean current on TRI in two-dimensional problems, where the wave vectors and the mean velocity vector belong to the same vertical plane (see e.g. [START_REF] Richet | Impact of a mean current on the internal tide energy dissipation at the critical latitude[END_REF]; [START_REF] Fan | Effect of background mean flow on psi of internal wave beams[END_REF]). In our case the mean azimuthal flow is perpendicular to the primary wave motion which occurs in a vertical (meridional) plane. We conjecture that a sufficiently strong mean azimuthal flow has an important effect on the radial extent of the observed azimuthal wave modes depending on the prograde/retrograde sense of motion of the flow and the waves. In the left panel of figure 16, we superimpose the evolution of the mean position of the vortex centers (marked by black dots) on the spatio-temporal diagram of the mean azimuthal flow. It can be seen that vortices are initially close to the surface of the inner cylinder as shown in figure 13 for the case a = 2.4 mm. At t around 150T 0 a narrow zone of prograde current starts to develop near the inner cylinder (visually it corresponds to the yellow "tongue" seen in the left panel of figure 16, see also the dashed profile in the right panel). Owing to interactions with this current, the vortex centers gradually move away from the inner cylinder while sitting at the radial coordinate roughly corresponding to the border between the prograde and retrograde currents. When the mean radial position of the vortex centers increases, the 7-vortex cluster evolves toward the 8-vortex configuration.

Let us note that the data presented in figure 15 and figure 16 permits to quantify that the prograde propagation speed of the vortex cluster, varying from 1/30 to 1/56 cycle per period T 0 , is approximately equal to or smaller than the prograde average flow speed observed in the narrow yellow "tongue" (see left panel of figure 16) toward the end of the simulation. This observation carries the suggestion that the components of the vortex cluster can be identified as Topographic Rossby Waves (TRWs) propagating in retrograde direction with respect to the current (while in prograde direction with respect to the rigid boundaries of the rotating setup).

In the experimental part of the paper (section 3) we have demonstrated that the cyclonic/anti-cyclonic asymmetry can be conveniently illustrated by the PDFs of the vertical vorticity components. Indeed, the vertical vorticity PDF provides a diagnostic from which the importance of redistributing processes can be obtained; its width (variance) tells how much vorticity is present in vortices and in inertial and topographic Rossby waves relative to the background ("planetary") vorticity, its skewness testifies about the dominance of concentrated strong cyclonic vorticity (in cyclones) relative to a more widespread background of weak anticyclonic vorticity. As seen in figure 6, in the experiment, the key contribution to the asymmetry is due to the low-frequency component of the vorticity field filtered with the help of the low-pass filter. In figure 17 we show the PDFs calculated for the numerically simulated vertical vorticity field. Note that the PDFs of the experimentally measured vertical vorticity component at a = 2.5 mm (black line) are fully consistent with the PDFs obtained for the numerically simulated vertical vorticity component obtained for a = 2.5 mm (green line) and a = 2.4 mm (purple line). Interestingly, the widest PDFs are obtained in numerical simulations performed at a = 2.0 mm, which is consistent with the snapshots shown in figures 7 and 13 (note that a larger range of the colorbar is used in these figures precisely for this amplitude).

The result of the direct numerical simulations is satisfying in the sense that the experimental pattern is successfully reproduced. It means that formation of a Polygonal Vortex Cluster in the geometric setup under consideration is a robust phenomenon. Moreover, the numerically observed merging of low-frequency wave components into a regular pattern of cyclonic vortices seems to be a plausible physical scenario for the formation of a slow manifold.

Many important issues remain unclear. In particular, it is not clear why the timescale of transient evolution of the flow toward a regular vortex cluster in the numerical simulations appears to be significantly longer than in the experiments. This is evident from the comparison of figure 15, where we show the temporal evolution of the instantaneous radial velocity as function of the azimuthal angle, with its experimental counterpart, figure 5. The formation of the 7-vortex cluster in numerical calculations takes about 100T 0 , roughly twice the time needed for formation of the experimental 8-vortex cluster. The 8-vortex cluster in simulations appears after roughly 300T 0 . This also raises a question: do we observe a fully saturated regime in the physical experiment and in the numerical simulations? Also note that the rate of slow cyclonic drift of the vortex cluster around the inner cylinder in numerical calculations is roughly 6 to 3 times higher than in experiment. Considering numerical calculations, we also notice that the rate of cyclonic drift systematically decreases with time, and the precession rate of the 8vortex cluster is significantly lower (by roughly 20%) as compared to that of the 7-vortex cluster. This decrease may be attributed to slow evolution of the mean radial position of vortices sitting at the border between the prograde and retrograde currents as illustrated in figure 16.

In the description of the numerical setup (see section 4.1) we have mentioned that the boundary condition (imposing a prescribed vertical velocity component at the fixed horizontal upper lid of the fluid domain) is not fully identical to the experimental situation with the deformable upper lid. This issue requires further investigation.

Conclusions

The investigation of inertial wave attractors in rotating fluids offers a number of possibilities concerning the geometric setup of the problem. A number of experimental studies [START_REF] Maas | Wave focusing and ensuing mean flow due to symmetry breaking in rotating fluids[END_REF][START_REF] Manders | Observations of inertial waves in a rectangular basin with one sloping boundary[END_REF][START_REF] Manders | On the three-dimensional structure of the inertial wave field in a rectangular basin with one sloping boundary[END_REF][START_REF] Brunet | Linear and nonlinear regimes of an inertial wave attractor[END_REF]) consider a rectangular box with one sloping wall placed at a rotating table which is similar to the geometry used to reveal the linear [START_REF] Maas | Observations of an internal wave attractor in a confined stably stratified fluid[END_REF]) and non-linear [START_REF] Scolan | Nonlinear fate of internal waves attractors[END_REF]Brouzet et al. 2016aBrouzet et al. , 2017a;;[START_REF] Davis | Succession of resonances to achieve internal wave turbulence[END_REF]) dynamics of internal wave attractors in stratified fluids. It has been realized that the inherent three-dimensionality of inertial waves is responsible for considerable secondary currents [START_REF] Maas | Wave focusing and ensuing mean flow due to symmetry breaking in rotating fluids[END_REF][START_REF] Manders | On the three-dimensional structure of the inertial wave field in a rectangular basin with one sloping boundary[END_REF]) and for a number of notable changes in the scenario of Triadic Resonance Instability [START_REF] Maurer | Generation and stability of inertia-gravity waves[END_REF][START_REF] Brunet | Linear and nonlinear regimes of an inertial wave attractor[END_REF]. In contrast, there is a rich theoretical literature which considers linear viscous regimes of inertial wave attractors in spherical liquid shells, where the flow is studied in the meridional cross-section only while the structure of the flow in the equatorial cross-section is supposed to be trivial [START_REF] Rieutord | Inertial waves in a rotating spherical shell[END_REF][START_REF] Rieutord | Wave attractors in rotating fluids: a paradigm for ill-posed cauchy problems[END_REF][START_REF] Rieutord | Viscous dissipation by tidally forced inertial modes in a rotating spherical shell[END_REF].

In the present paper, we consider numerically and experimentally the geometric setup which builds a bridge between the two above statements: the experiments are carried out in a (horizontal) annular and (vertical) trapezoidal domain which admits the existence of inertial wave attractor structures in meridional planes while leaving the freedom for the formation of a slow two-dimensional manifold which drains energy from the genuinely three-dimensional inertial wave field as result of an energy transfer. The experimental system is subject to axisymmetric forcing, and experiments are performed at low values of the global Rossby (Ro is of order 10 -3 ) and Ekman numbers (E is of order 10 -5 ). The main finding of the experimental part of the present study is the formation of a slow two-dimensional manifold in the saturated regime in the equatorial plane, which co-exists with an inertial wave attractor in the meridional plane. The two-dimensional manifold is represented by an ensemble of eight cyclonic vortices in a regular polygonal arrangement. The vortex cluster undergoes a slow cyclonic motion around the axis of rotation of the experimental system.

In the numerical part of the paper we reproduce the experimental setup and perform a series of simulations for a set of forcing amplitudes in order to obtain the experimental reference result. We observe a sequence of regimes: [i] linear regime with axisymmetric inertial wave attractor, [ii] onset of Triadic Resonance Instability, [iii] wave turbulence, [iv] formation of the slow manifold in the form of a regular polygonal vortex cluster, and [v] the slow cyclonic drift of the cluster. The Triadic Resonance Instability is observed both in the vertical plane (where it appears to be compatible with the exisiting knowledge) and in the horizontal plane where the secondary waves form the patterns possessing discrete azimuthal symmetry, with the secondary waves propagating cyclonically/anti-cyclonically. The concept of wave turbulence appears to provide a useful general framework for interpretation of the observed phenomena. However, there is an interesting specific feature: the motion observed in the meridional vertical planes seems to correspond to classic wave turbulence, evolving from "discrete" to "continuous" form, while the motion seen in the horizontal (equatorial) planes shows a strong trend toward development of discrete clearly defined azimuthal modes. An important issue for further investigation is the possible cross-interaction between prograde/retrograde azimuthal modes and prograde/retrograde mean azimuthal currents in the case of strong forcing when the currents are appreciable. It is noteworthy that the probability density functions of vertical vorticity calculated for experimental and numerical data are fully consistent, showing significant cyclonic/anti-cyclonic asymmetry.

Special attention has been paid to the development and long-term behaviour of the vortex cluster, which was analysed with the help of low-pass and, where appropriate, narrow-band filters. The results of the analysis support the idea that the vortex cluster emerges as the result of summation and/or merging of wave-like vorticity components in the presence of a mean current. We can speculate that the observed phenomena may be related to topographic Rossby waves. Indeed, the presence of propagating vortices in a typical topographic β set-up warrants a discussion in terms of topographic Rossby waves. The issue is a bit delicate for two reasons. On the one hand, for the conical shape of the bottom that we used (decreasing in depth inwards), the topographic Rossby waves are expected to propagate in retrograde direction (against the rotation sense of the tank), counter to what was actually observed in the laboratory experiment as well as in matching numerical experiments. As in the classical case studied originally by Rossby in the Earth's atmosphere, these waves may still propagate in prograde direction when advected by a strong background flow [START_REF] Rossby | Relation between variations in the intensity of the zonal circulation of the atmosphere and the displacements of the semi-permanent centers of action[END_REF]. The numerical computations that support the laboratory experiments indeed show that at the radius where the vortices sit, a strong prograde mean flow has formed. On the other hand, the mean flow is not necessarily formed by the mixing of angular momentum due to breaking of focused and amplified inertial waves, but could also be generated by non-linear interactions of periodic flows within a viscous boundary layer (see, e.g., the works of [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF]; [START_REF] Sauret | Fluid flows in a librating cylinder[END_REF]). In our study, it is important to highlight that this mean flow itself has a complex radial structure that contains strong retrograde as well as eventually strong prograde parts. We believe this aspect, the precise generation and evolution of the radial mean flow structure by focused inertial waves, to be sufficiently complex to require further investigation. In numerical simulations (not reported in this paper) we have also made a limited preliminary study on robustness of the observed effects by considering a conical bottom with a smaller slope (30 • instead of 45 • ). The calculations have demonstrated that a Polygonal Vortex Cluster also arises at comparable input forcing. Additionally, the form of forcing has been changed from the cubic profile shown in figure 2 to a linear profile (under the condition that the forcing should be volume-conserving). Again, the computations have demonstrated the emergence of a vortex cluster, and we can therefore conclude that the observed phenomenon is indeed robust.

Although the experimental reference regime is numerically reproduced, and the role of the forcing amplitude is clarified, the present study raises a number of interesting general issues, in particular, [i] how the vortex cluster regime evolves when the key parameters of the problem (Rossby and Ekman numbers, geometric aspect ratio and the particular type of forcing) are varied in a broad range, and [ii] whether or not the observed regime might be relevant to realistic geo-and astrophysical systems.

Cyclonic clusters arranged in the form of regular polygons have been reported for the polar regions of large planets, e.g. Jupiter [START_REF] Adriani | Clusters of cyclones encircling jupiters poles[END_REF], demonstrating remarkably persistent long-term behavior [START_REF] Adriani | Twoyear observations of the jupiter polar regions by jiram on board juno[END_REF]. This phenomenon has a different physical origin as compared to our experiments. It has been argued that the pattern is captured by shallow-water models, e.g. [START_REF] Cho | The emergence of jets and vortices in freely evolving, shallow-water turbulence on a sphere[END_REF]; [START_REF] Scott | Forced-dissipative shallow-water turbulence on the sphere and the atmospheric circulation of the giant planets[END_REF]. However, the key puzzle of vortex clusters remains [START_REF] Adriani | Clusters of cyclones encircling jupiters poles[END_REF]): "The manner in which the cyclones persist without merging and the process by which they evolve to their current configuration are unknown". Recently, [START_REF] Reinaud | Three-dimensional quasi-geostrophic vortex equilibria with m-fold symmetry[END_REF] has shown numerically that a system of m quasi-geostrophic vortices equally distributed over a ring whose center is already occupied by a vortex can be stable under certain conditions. In our experiment, the inner cylinder may play a similar role to this central vortex. Hence, the "toy system" proposed here may help to shed light on the stability of polygonal vortex clusters and their possible existence not only in a shallow "atmosphere", but also in liquid interiors of rotating natural systems geometrically compatible with the existence of inertial wave attractors.

Figure 1 .Figure 2 .

 12 Figure1. Schematic of the experimental apparatus in (a) a vertical cross-section, and (b) a horizontal cross-section. The hashed area is out of the experimental domain. The trapezoidal region of interest in the vertical plane is located between the inner and the outer cylinders, of radii R0 = 5 cm and R1 = 20.2 cm, respectively.

Figure 3 .

 3 Figure 3. Fields of radial vr and azimuthal v θ velocity and vertical vorticity ξz in the horizontal plane (at 20 cm depth) and of vertical velocity vz in the vertical plane. The presented quantities are filtered around ω = ω0. Positive vorticity corresponds to cyclonic motion. The parallelogram with arrows in the vertical plane shows the theoretical attractor in which the energy propagates clockwise.

Figure 4 .

 4 Figure 4. Fields of radial vr and azimuthal v θ velocity and vertical vorticity ξz in the horizontal plane (at 20 cm depth) and of vertical velocity vz in the vertical plane. The presented quantities are filtered around ω = 0 with the low-pass filter having the upper cut-off frequency ωc = ω0/3. Positive vorticity corresponds to cyclonic motion.

Figure 5 .

 5 Figure 5. Temporal evolution of the radial velocity as a function of azimuth, θ, obtained by taking profiles of radial velocity at radius r = 8 cm in the horizontal plane at 20 cm depth. To guide the eye, we added a solid line showing that the cluster rotates half a turn (π rad) in (88 ± 2)T0.

Figure 6 .

 6 Figure 6. PDFs of the vertical vorticity component at different times in the experiment. The curves (from top to bottom) correspond to processing of the raw signal, and the signal filtered around ω = ω0 and ω = 0.

Figure 7 .

 7 Figure 7. Snapshots of vertical vorticity field, obtained by numerical simulations, in horizontal plane at mid-depth z = -20 cm (upper row) and corresponding vertical velocity fields in vertical plane y = 0 (lower row) taken at t = 100T0 at different values of the forcing amplitude a. Note that both in experiments and in numerical computations we denote cyclonic vorticity and prograde azimuthal velocity as positive. Thus the experimental and numerical patterns have the same colour coding and can be directly compared regardless of the direction of the background rotation.

Figure 8 .

 8 Figure 8. Time-frequency diagram of the logarithm (color) of the normalised spectrum of the radial velocity, Sr(ω, t)/ S0 , where S0 is the time averaged spectrum at ω/ω0 = 1, calculated with the help of equation (4.1) for a ring-shaped region around the inner cylinder (r ∈ [5, 10] cm). The last image in the lower row corresponds to the time-frequency diagram calculated for the experimental data obtained at a = 2.5 mm. The white dashed line indicates the cut-off frequency f .

Figure 9 .

 9 Figure 9. The upper four images present the normalized energy spectra E (α, ω)/E (ω) at different values of the forcing amplitude. Colors indicate the levels of normalized energy spectra.The white dashed lines correspond to the dispersion relation ± cos α = ω/f . Note that the energy peaks are not always localized at the dispersion curve, owing partially to the departure of the wave vectors from (r, z) plane and to the contribution of slaved (evanescent) waves (see discussion in the text). The lower four images show the energy spectra in the horizontal plane E ⊥ (kr, n θ ) on a logarithmic scale at different values of the forcing amplitude. By combining information from the spectra calculated in the vertical and horizontal planes one can see the evolution of wave regimes from nearly linear at a = 0.2 mm to wave turbulence evolving from a "discrete" form at a = 1.0 to a more "continuous" form at 2.0 and a = 2.4 mm. The latter is particularly well seen in terms of E ⊥ (kr, n θ ). In all cases the axisymmetric component of the wave field (corresponding to n θ = 0) is most significant.

Figure 10 .

 10 Figure 10. The components of the inertial wave field (displayed in terms of the vertical velocity component) filtered at frequencies ω0, ω1 and ω2 (upper row), with the corresponding phase patterns (lower row). Only one branch of the wave attractor is shown. The corresponding vector triad is shown in figure 12 (a).

Figure 11 .

 11 Figure 11. The components of the vertical vorticity field ξz filtered at frequencies ω0, ω1 and ω2 (upper row), with corresponding filtered wave fields as function of radial coordinate r and azimuthal coordinate θ, with amplitude shown in the middle row and phase shown in the lower row. Note that the primary wave is axisymmetric and propagates radially, while the secondary waves corresponding to frequencies ω1 and ω2 propagate azimuthally in the prograde (cyclonic) and retrograde (anti-cyclonic) directions, respectively, and the azimuthal components of the wave vectors have the same length. The corresponding vector triad is shown in figure 12 (b).

Figure 13 .

 13 Figure13. The low-frequency pattern of the vertical vorticity field filtered by low-pass filter with the cut-off frequency ωc = ω0/3 at different values of the forcing amplitude in the DNS and for the experimental case. The images correspond to the time t = 100T0. In the supplementary material we provide a short video of this pattern demonstrating prograde motion of the wave pattern at a = 1.0 mm, prograde/retrograde motion close to inner/outer cylinder respectively at a = 2.0 mm and prograde motion of the vortex cluster at a = 2.4 mm (in numerical simulations) and at a = 2.5 mm (in experiments). Note that for experimental data we take a mirror image, so that the background rotation and prograde motion are anti-clockwise. The numerically calculated wave/vortex patterns are rotating around the central cylinder while interacting with the mean azimuthal currents shown in the right panel of figure16.

-Figure 14 .

 14 Figure 14. The time-frequency diagram corresponding to the case a = 2.4 mm (upper image). Lower row (from left to right) shows the snapshots of the vertical vorticity fields filtered in narrow frequency windows 3 and 4 around t = 170T0, and via low-pass filter schematically shown as 5 around t = 310T0. Patterns 3 and 4 propagate in prograde and retrograde directions, respectively. Note that the discrete low-frequency components seen in the time-frequency diagram exhibit a trend toward merging at t > 225T0, and filter 5 is applied to the result of the merging process. The result of filtering with 6 (low pass-filter, earlier times) is shown in figure 13 (c).

Figure 15 .

 15 Figure 15. Same as experimental figure 5 but computed from the numerical data. It can be seen that in the numerical simulations the prograde motion of the vortex cluster is significantly faster than in the experiment case. The three lines indicate, from left to right, an azimuthal drift of the vortex cluster at the rate of half a turn in (15 ± 2)T0, (22 ± 2)T0, and (28 ± 2)T0. The rate of the azimuthal drift thus decreases with time.

Figure 16 .

 16 Figure 16. (a) Spatio-temporal diagram representing the numerically computed time-history of radial distribution of mean non-dimensional azimuthal velocity component V θ = V θ θ,t • T0/(2πr). Positive/negative values correspond to prograde/retrograde current. Black dots show the mean radial position of the centres of vortices in the vortex cluster. (b) Radial distribution of mean azimuthal velocity corresponding to fixed moment of time t = 100T0 in numerical simulations performed at different values of the forcing amplitude a. The "slow manifold" patterns shown in figure 13 and in supplementary material correspond to t = 100T0. These patterns are slowly rotating while interacting with the mean azimuthal current. In (b), for a = 2.4 mm we present also the radial profile of the azimuthal current taken at t = 310T0 (purple dashed line): note the presence of the prograde current corresponding to the yellow "tongue" in the left panel (a). The measurement times t = 100T0 and t = 310T0 are indicated by white dashed lines in the left panel. In both figures, a grey rectangle is added to represent the fluid domain at rest contained in the inner cylinder. (c) Scaling of the maximum of V θ with the square of the forcing amplitude, a 2 .

Figure 17 .

 17 Figure17. PDFs of the vertical vorticity components calculated for the time instant around t = 100T0 at different values of the oscillation amplitude: (a) corresponds to the raw signal, and (b) corresponds to the signal filtered via the low-pass filter. The PDF for the experimental data obtained at a = 2.5 mm is shown in light blue. The PDFs for the numerical data obtained at a = 2.5 mm (green) and a = 2.4 mm (purple) are fully consistent with the experimental result (drawn in black).

Acknowledgments

Vortex cluster arising from an axisymmetric inertial wave attractor This work was supported by the grant ANR-17-CE30-0003 (DisET) and by the LABEX iMUST (ANR-10-LABX-0064) of Université de Lyon, within the program Investissements d'Avenir (ANR-11-IDEX-0007), operated by the French National Research Agency (ANR). This work was also supported by a grant from the Simons Foundation (651475, TD). It has been achieved thanks to the resources of PSMN from ENS de Lyon. EE gratefully acknowledges his appointment as a visiting scientist at ENS de Lyon during the experimental campaign. EE and NS acknowledge support from the Russian Science Foundation (Project 20-11-20189) during work on the manuscript and complementary data processing. The direct numerical simulations were performed at the Lomonosov Supercomputer Center. The authors deeply thank the anonymous reviewers whose constructive comments greatly helped to improve the overall quality of this manuscript and the scope of our study. Declaration of Interests. The authors report no conflict of interest.