Valentin Ageorges

Jorge Peixinho

Gaële Perret

Flow and air-entrainment around partially submerged vertical cylinders

In this study, a partially submerged vertical cylinder is moved at constant velocity through water, which is initially at rest. During the motion, the wake behind the cylinder induces free-surface deformation. Eleven cylinders, with diameters from D = 1.4 to 16 cm, were tested under two different conditions: (i) constant immersed height h and (ii) constant h/D. The range of translation velocities and diameters are in the regime of turbulent wake with experiments carried out for 4500 < Re < 240 000 and 0.2 < Fr < 2.4, where Re and Fr are the Reynolds and Froude numbers based on D. The focus here is on dragforce measurements and relatively strong free-surface deformation up to air-entrainment. Specifically, two modes of air-entraiment have been uncovered: (i) in the cavity along the cylinder wall and (ii) in the wake of the cylinder. A scaling for the critical velocity for air-entrainment in the cavity has been observed in agreement with a simple model. Furthermore, for Fr > 1.2, the drag force varies linearly with Fr.

I. INTRODUCTION

The flow past ships [START_REF] Rabaud | Ship Wakes: Kelvin or Mach Angle?[END_REF][START_REF] Moisy | Mach-like capillary-gravity wakes[END_REF], or an emerged body such as a bridge pillar, is a fundamental, familiar, and fascinating sight. Measurements and modeling of this simple flow can have relevance for the floating offshore structures and renewable energy systems [START_REF] Drew | A review of wave energy converter technology[END_REF][START_REF] López | Review of wave energy technologies and the necessary power-equipment[END_REF]. The disordered motion of the free-surface in the wake of the cylinder is related to the drag and lift effect experienced by the cylinder. The wake pattern is responsible for the surface waves and probably the air-entrainment. Many experiments and numerical works have considered a single vertical cylinder geometry in a fully submerged situation [START_REF] Schewe | On the force fluctuations acting on a circular cylinder in crossflow from subcritical up to transcritical Reynolds numbers[END_REF][START_REF] Braza | Turbulence properties in the cylinder wake at high Reynolds numbers[END_REF][START_REF] Achenbach | On vortex shedding from and rough cylinders in the range of Reynolds numbers 6 × 10 3 to 5 × 10 6[END_REF][START_REF] Breuer | A challenging test case for large eddy simulation: High Reynolds number circular cylinder flow[END_REF][START_REF] Tritton | Experiments on the flow past a circular cylinder at low Reynolds numbers[END_REF]. These works reported the drag-force measurements on the circular cylinders as a function of the Reynolds number Re defined as Re = U D/ν, with U being the velocity, D the cylinder diameter, and ν the kinematic viscosity of the fluid. In addition, Schewe [START_REF] Schewe | On the force fluctuations acting on a circular cylinder in crossflow from subcritical up to transcritical Reynolds numbers[END_REF] showed the occurrence of a drag crisis at Re ≈ 350 000, which separates two distinct regions of a nearly constant dimensionless drag coefficient, C D , defined later. In the subcritical regime, for Re < 350 000, the drag coefficient and the vortex-shedding frequency f , characterized by the Strouhal number St = f D/U , are almost constant at C D ≈ 1.2 and St ≈ 0.2. In this regime, the flow is characterized by a turbulent wake with a laminar boundary layer separation [START_REF] Sumer | Hydrodynamics around Cylindrical Strucures[END_REF].

In contrast to fully submerged situations, few studies have been done taking into account the effect of a free-surface, which modifies the pressure distribution around the cylinder. From a numerical point of view, the authors of Refs. [START_REF] Kawamura | Large eddy simulation of a flow past a free surface piercing circular cylinder[END_REF][START_REF] Yu | Large eddy simulation of flow past free surface piercing circular cylinders[END_REF][START_REF] Koo | Reynolds and Froude number effect on the flow past an interface-piercing circular cylinder[END_REF] have computed the flow around a wallmounted cylinder, the evolution of C D , and the free-surface elevation for Re varying from 27 000 up to 458 000. They have shown that the vorticity and C D vary with depth along the cylinder and diminish close to the free-surface. Additionally, C D is strongly influenced by free-end effects for finite-height cylinders [START_REF] Sumner | Wake structure of a finite circular cylinder of small aspect ratio[END_REF]. Specifically, Gonçalves et al. [START_REF] Gonçalves | Flow around circular cylinders with very low aspect ratio[END_REF] reported experimental results in a recirculating water flume for cylinders piercing the free-surface with the aspect ratio 0.1 < h/D < 2, where h is the immersed cylinder height. They focused on forces and particle image velocimetry measurements for 10 000 < Re < 50 000 and 0.07 < Fr < 0.36; the Froude number Fr is defined as In recent years, numerical studies have reproduced the flow and drag around free-end cylinders for different h/D [START_REF] Yu | Large eddy simulation of flow past free surface piercing circular cylinders[END_REF][START_REF] Koo | Reynolds and Froude number effect on the flow past an interface-piercing circular cylinder[END_REF][START_REF] Ducrocq | Flow and drag force around a free surface piercing cylinder for environmental applications[END_REF][START_REF] Benitz | CFD simulations and experimental measurements of flow past free-surface piercing, finite length cylinders with varying aspect ratios[END_REF], but there is a lack of data for Re around the drag crisis (Re ≈ 350 000) and almost no data for strong free-surface deformations and rupture. An important numerical output is the mean interface elevation [START_REF] Yu | Large eddy simulation of flow past free surface piercing circular cylinders[END_REF][START_REF] Koo | Reynolds and Froude number effect on the flow past an interface-piercing circular cylinder[END_REF]. Specifically, Keough et al. [START_REF] Keough | Time resolved measurements of wake characteristics from vertical surface-piercing circular cylinders[END_REF] analyzed the evolution of the fountain height ahead of the cylinder for Re varying between 45 000 and 380 000. They concluded a dependency of the fountain height with Fr 2 and revealed that this evolution follows the Bernoulli equation for small velocity. In addition, at the downstream side of the cylinder, the maximum cavity depth L is also believed to scale with Fr 2 [START_REF] Chaplin | Steady flow past a vertical surface-piercing circular cylinder[END_REF]. The cavity depth L is schematized in Fig. 1(c). Depending on the cylinder diameter, the free-surface downstream is locally subjected to strong deformations and air-entrainment in the liquid phase. Benusiglio [START_REF] Benusiglio | Indiscrétions aux interfaces[END_REF] recently demonstrated a velocity threshold for air-entrainment in the cavity of a translating vertical cylinder. This phenomenon can induce drag reduction effects and it is independent of the drag crisis introduced before. At larger Fr, Hay [START_REF] Hay | Flow about semi-submerged cylinders of finite length[END_REF] reported an extensive study of the ventilated regime, where air reaches the bottom of the cylinder and creates a ventilated cavity with large waves and sprays above the nominal free-surface. Modern studies of the artificially or atmospherically ventilated structures (plates or hydrofoils) close to a free-surface are concerned with the stability of ventilated cavities (e.g., Ref. [START_REF] Harwood | Ventilated cavities on a surface-piercing hydrofoil at moderate Froude numbers: Cavity formation, elimination and stability[END_REF]) or properties such as the void fraction measurements [START_REF] Mäkiharju | Time-resolved twodimensional X-ray densitometry of a two-phase flow downstream of a ventilated cavity[END_REF].

The goal of the present investigation is to quantify the onset of air-entrainment effects. In addition, the present study adds new experimental results and discusses new drag measurements carried out for 4500 < Re < 240 000 and 0.2 < Fr < 2.4 with a specific focus on the onset of air-entrainment [START_REF] Harwood | Ventilated cavities on a surface-piercing hydrofoil at moderate Froude numbers: Cavity formation, elimination and stability[END_REF][START_REF] Kiger | Air-entrainment mechanisms in plunging jets and breaking waves[END_REF][START_REF] Chirichella | Incipient air entrainment in a translating axisymmetric plunging laminar jet[END_REF][START_REF] Peixinho | Free surface flow between two horizontal concentric cylinders[END_REF]. This paper is composed of two parts. In the first part, the experimental setup is described together with an optical method for free-surface height reconstruction. In the second part, results are presented with the drag-force measurements and the determination of the velocity threshold for air-entrainment in the cavity. Finally, the air-entrainment results are compared with a simple model based on the balance between acceleration and gravity.

II. EXPERIMENTAL SETUP

The cylinders are hollow and made of polymethyl methacrylate (PMMA) with external diameters of 30, 40, 50, 60, 70, 80, 90, 110, and 160 mm. Their total height is 65 cm and their wall thickness is 5 mm. The Young modulus of PMMA is 3.1 GPa. There are additional hollow cylinders of 14 and 25 mm made respectively of copper and aluminum. All the cylinders are clogged at their ends. These cylinders are attached vertically, perpendicular to the water free-surface, and translated on a carriage riding horizontally along the flume of 34 m in length, 90 cm in width, and 120 cm in height. The flume is shown schematically in Fig. 1(a). The carriage (162 kg) is driven by a 600-V ac motor (34 kg), producing a maximum torque of 102 N/m. A protocol is applied where first the velocity increases linearly along a distance of 3.5 m. Then, the speed is maintained constant over a distance of 9.5 m, where force measurements are acquired in steady state. Finally, the speed decreases linearly over a distance of 2 m.

Two windows are located on the sidewalls of the flume, allowing local flow visualization around the cylinder. The determination of the critical velocity for air-entrainment behind the cylinder is based on images from a 2000 × 2000 pixel CCD camera (Basler aca2040) with focal length of 50 mm placed at 6 m from the start of the cylinder run. At this position along the flume, the force measurement (described later) has reached a steady-state value. The spatial resolution of the CCD camera is 300 μm per pixel and its acquisition frequency is 100 frames per second. Hence, bubbles smaller than 300 μm cannot be detected. A typical image of the cylinder passing through the field of view of the camera is shown in Fig. 1(c).

The working fluid is tap water and its height in the flume is constant: H = 40 ± 0.1 cm. The temperature is monitored; hence the viscosity, the density ρ, and the surface tension σ of water are corrected to take into account these variations. Indeed between all experimentations, the temperature of water changes from 16.1 to 21.3 • C. In this range of temperature, density variations represent 0.1% [START_REF] Tanaka | Recommended table for the density of water between 0 • C and 40 • C based on recent experimental reports[END_REF], while viscosity variations represent 13.3% [START_REF] Huber | New international formulation for the viscosity of H 2 O[END_REF]. The correction allows us to keep the error for Re between 0.8 and 9.8%. The Bond number [START_REF] Aristoff | Water entry of small hydrophobic spheres[END_REF], Bo = ρgD 2 /σ , varies between 26 and 3450, which implies that the capillary effects are weak.

Two series of experiments were performed: (i) at constant immersion depth h = 23 ± 0.1 cm leading to h/D ratios varying between 1.44 and 16.5 and (ii) at constant ratio h/D = 2.55 ± 0.05. In total, 372 runs were performed, as summarized in Table I, covering a wide range of Re and Fr. The maximum blockage width of the experiment, D/W , where W is the flume width, is 0.177, a value below the critical blockage width of 0.2 for the confinement effect on three-dimensional dynamics proposed by Griffith et al. [START_REF] Griffith | Vortex shedding and three-dimensional behavior of flow past a cylinder confined in a channel[END_REF]. In addition, the International Towing Tank Conference (ITTC) 2002 [START_REF] Müller-Graf | Six contributions to the hydrodynamics of high speed marine vehicles[END_REF] requires a vertical confinement, (Hh)/h > 0.8. For the current experiments, at constant h, the vertical confinement is 0.74, and for constant h/D, this criteria is satisfied for D < 9 cm.

A. Force measurements

A piezoelectric sensor (Kistler 9327C) placed above the cylinder [see drawing in Fig. 1(b)] measures the axial force F x acting on the whole cylinder with a sensitivity of -7.8 pC/N. The sensor is connected to a charge amplifier equipped with a low-pass filter of 30 Hz. The contribution of air drag is neglected. Indeed, when h is minimum, the maximum error on C D is about 0.8%. Force signals are filtered using a moving average on 20 points. In order to measure the effect of air-entrainment on drag forces, a time average is performed on F x during 5 s at constant translating velocity. It is denoted Fx . The drag coefficient C D can thus be defined as

C D (t ) = 2 Fx ρhDU 2 .
(

) 1
The Strouhal number, defined earlier, is found to be close to 0.2 ± 0.025 for D < 6 cm, using the peak of f from the spectrum of lift-force signal. However, for larger diameters, the signal of lift force is noisy, presumably due to the turbulent flow or vibrations of the carriage. Moreover, the range of Re studied here and the materials used here suggest natural frequencies which are well separated from those associated with vortex-induced vibration, which is the motion induced on the elastic cylinder interacting with the turbulent fluid flow. Such vortex-induced vibration phenomena of the elastic cylinder deformation would act like an oscillator inducing several modes corresponding to different vortex sheddings observed by several authors [START_REF] Fujarra | Vortex-induced vibration of a flexible cantilever[END_REF][START_REF] Sarpkaya | A critical review of the intrinsic nature of vortex-induced vibrations[END_REF][START_REF] Williamson | Vortex-induced vibrations[END_REF][START_REF] Franzini | Comparison between force measurements of one and two degrees-of-freedom VIV on cylinder with small and large mass ratio[END_REF].

B. Optical method

The parameter that needs to be measured is the elevation of the free-surface. The height at rest H is well controlled. However, when the cylinder passes through the observation area, the light from a vertical laser sheet located above the flume and pointing to the free-surface, represented in red in Fig. 1(b), is shifted vertically from its initial position. The shift is the difference between the interface at rest and the deformed free-surface. It is proportional to the elevation of the free-surface. The images are captured by a camera placed on the top of the flume with a low incident angle to the free-surface. The proportionality is calibrated by recording images of the laser sheet at five known positions over 25 mm. This technique is able to detect free-surface heights with a precision estimated to ±0.1 mm.

Figure 2 presents examples of free-surface reconstructions using the optical method for Re = 8500. Two Fr are shown: (a) Fr = 0.5 and (b) Fr = 1.6. When the Fr increases, the free-surface experiences elevations of larger amplitude. The diagrams allow one to distinguish the fountain height upstream of the cylinders. It is also possible to see the angle of the V-shaped wakes. For Fr = 1.6 the V-shaped wake induces a planar wave downstream of the cylinder. Note that this method and other optical methods (e.g., Ref. [START_REF] Gomit | Développement de techniques de mesure de surfaces libres par moyens optiques: Application à l'analyse de l'écoulement généré par un modèle de bateau en bassin des carènes[END_REF]) are restricted by the dispersion of the laser due to the cylinder's turbulent wake. Hence, the present technique is limited for small Fr and Re.

III. RESULTS AND DISCUSSION

The results consist of drag measurements, surface elevation measurements, and image analysis. These have been obtained for different cylinder diameters translated at various speeds with constant immersion h and constant ratio h/D. The different modes of air-entrainment are depicted and used to detect the critical velocity for air-entrainment in the cavity. The scaling of the velocity can be explained using a simple model described below.

A. Drag measurements

In this section, force measurements are reported in a range of Re and Fr, such that the transition to air-entrainment is observed. Figure 3 presents drag coefficient measurements for different cylinder diameters as a function of Re, at constant h = 23 ± 0.1 cm and at constant ratio h/D = 2.55 ± 0.05. Solid (empty) symbols correspond to cases without (with) air-entrainment in the cavity L behind the cylinder. The error bars correspond to the standard deviation of C D . The results obtained are compared with the results of Schewe [START_REF] Schewe | On the force fluctuations acting on a circular cylinder in crossflow from subcritical up to transcritical Reynolds numbers[END_REF] for single-phase flow. Additional results of Refs. [START_REF] Yu | Large eddy simulation of flow past free surface piercing circular cylinders[END_REF][START_REF] Gonçalves | Flow around circular cylinders with very low aspect ratio[END_REF] on surface-piercing cylinders are represented for the same range of Re and h/D. In accordance with previous studies, the drag coefficient is found to be about 30% smaller for surface-piercing cylinders than for single-phase flows. The numerical results of Yu et al. [START_REF] Yu | Large eddy simulation of flow past free surface piercing circular cylinders[END_REF] at three different Re numbers and a constant ratio h/D = 4 are in good agreement with the present measurements. The results of Gonçalves et al. [START_REF] Gonçalves | Flow around circular cylinders with very low aspect ratio[END_REF] for a ratio of h/D = 2 show a nearly constant drag coefficient at higher values than the present study. In their experiments the cylinder is attached to the walls of the flume and the aspect ratio h/D is varied by varying the water depth. Moreover, their maximum Fr based on the cylinder diameter is Fr = 0.36, lower than the present values and below the threshold of air-entrainment. This may explain their highest C D values.

In the present study, the drag coefficient first increases with the Re and then decreases after air-entrainment has occurred in the cavity downstream from the cylinder. Air-entrainment is defined when the first bubble is observed or monitored from the camera images (i.e., bubble larger than 300 μm). In Fig. 3, for small h/D or large diameters D, C D values first increase with Re after the appearance of the first air bubbles. Then, C D values decrease from 10% for the 4-cm-diameter cylinder to 26% at maximum for the 11-cm-diameter cylinder. The existence of the drag reduction effect could be explained by the modification of pressure distribution along the cylinder.

In order to study the influence of the gravity effect, the drag coefficient C D is also plotted as a function of Fr for a constant immersion depth, h = 23 ± 0.1 cm, in Fig. 4(a) and for a constant ratio, h/D = 2.55 ± 0.05, in Fig. 4(b). As in Fig. 3, solid (empty) symbols correspond to measurements without (with) air-entrainment in the cavity. In both cases, the evolution of the drag coefficient clearly shows a transition occurring at Fr ≈ 1.2. For Fr < 1. with Fr. A similar linear behavior has also been observed in granular flows around an immersed cylinder [START_REF] Chehata | Dense granular flow around an immersed cylinder[END_REF].

B. Air-entrainment

Two mechanisms of air-entrainment are revealed: (i) in the wake thanks to vortices and (ii) in a cavity behind the cylinder. For 1.4 < D < 4 cm (5200 < Re/Fr < 25 000), air-entrainment occurs in the cavity along the cylinder exclusively. This cavity [START_REF] Yu | Large eddy simulation of flow past free surface piercing circular cylinders[END_REF][START_REF] Koo | Reynolds and Froude number effect on the flow past an interface-piercing circular cylinder[END_REF] forms viscous cusps where airentrainment is systematically injected. Once air bubbles form below the cavity, they are carried downstream into the wake. The bubbles created by air-entrainment in the cavity have an equivalent diameter ranging from 6 mm for D = 3 cm to 8.5 mm for D = 16 cm.

For 4 < D < 16 cm, air-entrainment first appears in the wake when Re is increased. Then, as the velocity increases, air-entrainment also appears in the cavity. For a cylinder of 5 cm in diameter, Fig. 5(a), at Re = 34 670, shows free-surface deformation without air-entrainment. In Fig. 5(b) at Re = 39 700, air-entrainment occurs in the cylinder wake. Finally, in Fig. 5(c), at Re = 45 700, the cavity induces air-entrainment more intensely. For D = 16 cm and higher Re, the mechanisms of air-entrainment presented in Fig. 6 are similar. Air-entrainment appears in the wake [see Fig. 6(a)] and then in the cavity [see Fig. 6(b)]. The corresponding movies are available as Supplemental Material [38]. A summary of all our observations is presented in Fig. 7, for both constant h and h/D. In both cases, several regimes are observed where the different modes of air-entrainment are (i) in the wake, (ii) in the cavity, and (iii) in the wake and the cavity. Boundaries can be drawn between the different regimes and it is the purpose of the model below to explain the nature of some of these boundaries. The main difference between the two diagrams is the area of the wake region, which is wider for large diameters and constant h. For small h/D, the onset of air-entrainment occurs at smaller Fr presumably because of the size of the recirculation below the free-surface [START_REF] Gonçalves | Flow around circular cylinders with very low aspect ratio[END_REF][START_REF] Benitz | CFD simulations and experimental measurements of flow past free-surface piercing, finite length cylinders with varying aspect ratios[END_REF]. Additionally, the fact that the two diagrams are almost indistinguishable suggests that the immersed height h FIG. 6. View from the side below-surface of the cylinder of 16 cm in diameter. (a) Air-entrainment in the wake at Re = 122 500 or Fr = 0.6; bubbles encircled in blue, see Supplemental Material, movie 4 [38]. (b) Air-entrainment in the cavity behind the cylinder at Re = 199 300 or Fr = 1; bubbles encircled in red, see Supplemental Material, movie 5 [38]. The dashed lines represent observed bubble trajectories. has little effect on the air-entrainment when h > 23 ± 0.1 cm and h/D > 2.55 ± 0.05. Hence, airentrainment is mainly due to the free-surface dynamics and deformation.

Indeed, the main feature associated with air-entrainment is the deformation of the cavity behind the cylinder. The optical method described earlier and the direct observations of graduations along the cylinder allowed us to measure the cavity depth L for various conditions. Specifically, the observations of graduated cylinders are reported in Fig. 8: L/D as a function of Fr 2 Re. L represents the mean value from ten snapshots [see Figs. 8(c) and 8(d)]. The error bars are the maximum observed deviation from the mean value. Clearly, the behavior is nonlinear, typically after air-entrainment takes place. However, in the range of Fr 2 Re before air-entrainment, L/D varies linearly with Fr 2 Re. The proportionality coefficients seem to decrease with the diameter. These measurements are insufficient to provide a definite scaling of the cavity height. In fact, the cavity is unsteady and the liquid motion is more disordered as the velocity increases.

Air-entrainment phenomenon in the cavity could be explained using the model developed by Benusiglio [START_REF] Benusiglio | Indiscrétions aux interfaces[END_REF] where the acceleration of the cavity γ is balanced by the gravitational acceleration g. γ is controlled by the vortex shedding frequency f , calculated from the Strouhal number. In the range of Re of the present study, St is assumed constant (∼0.2) [START_REF] Schewe | On the force fluctuations acting on a circular cylinder in crossflow from subcritical up to transcritical Reynolds numbers[END_REF], although some studies have shown its variability upon the cylinder aspect ratio [START_REF] Benitz | CFD simulations and experimental measurements of flow past free-surface piercing, finite length cylinders with varying aspect ratios[END_REF][START_REF] Farivar | Turbulent uniform flow around cylinders of finite length[END_REF][START_REF] Fox | Fluid-induced loading of cantilevered circular cylinders in a low-turbulence uniform flow. Part 1: Mean loading with aspect ratios in the range 4 to 30[END_REF]. The acceleration γ of the cavity of depth L is defined as follows:

γ = L f 2 .
(

) 2
In Eq. ( 2), L is determined using Bernoulli's equation between two positions upstream and downstream of the cylinder. At the upstream position the mean flow velocity is U , whereas at the downstream position the flow velocity is taken as the maximum velocity U of a Rankine vortex /2πa [START_REF] Li | Experimental and numerical investigation of free surface vortex[END_REF], where is the circulation around half of the cylinder DU /2 and a represents the core radius of the vortex. The parameter a scales with D/ √ Re, the boundary layer thickness [START_REF] Schlichting | Boundary-Layer Theory[END_REF]. With these assumptions, Benusiglio [START_REF] Benusiglio | Indiscrétions aux interfaces[END_REF] proposed a scaling for L:

L ∝ DU 3 gν , (3)
which has not been validated experimentally or numerically. However, our current measurements (see Fig. 8), using the direct depth observations, seem to be consistent with L/D ∝ Fr 2 Re, which is equivalent to Eq. (3). Finally, air-entrainment occurs when γ g; therefore Eq. ( 2) leads to the critical velocity U c for air-entrainment in the cavity L:

U c ∝ (νg 2 ) 1 5 D 1 5 . ( 4
)
Figure 9(a) shows the critical velocity U c for air-entrainment in the cavity L as a function of D at constant immersion h = 23 ± 0.1 cm and constant h/D = 2.55 ± 0.05. The error bars of U c correspond to the velocity step between the maximum velocity when air-entrainment is not observed and the velocity where the first bubble appears. The evolution of U c is fitted with power laws of D 0.2±0.02 for h/D 2.55 and D 0.22±0.02 for h 23 cm. Experimental results show good agreements with the theoretical model introduced with Eq. ( 4). The simple model based on energy conservation with Bernoulli's equation predicts the appearance of air-entrainment in the cavity. Additional data from Benusiglio [START_REF] Benusiglio | Indiscrétions aux interfaces[END_REF] are also plotted and fit well the present data for D 6 cm. For larger diameters, the critical velocity given by Benusiglio [START_REF] Benusiglio | Indiscrétions aux interfaces[END_REF] refers to air-entrainment in the cylinder wake like in Fig. 5(b), while our present data refer only to air-entrainment in the cavity. Figure 9(b) presents the evolution of Re c as a function of Fr c , based on the critical velocity U c . The dimensionless form of Eq. ( 4) is Re c ∝ Fr -4 c . The evolution of Re c is fitted with a power law of Fr -3.3±0.62 c for h = 23 ± 0.1 cm and Fr -4.7±0.51 c for h/D = 2.55 ± 0.05, in good agreement with the prediction of the model. The red shaded area refers to the range of dimensionless parameters where no air-entrainment in the cavity occurs. The black dashed-dotted line represents the threshold Re/Fr 25 000 above which air-entrainment occurs in the wake and in the cavity. Hence, air-entrainment requires large inertial effects or strong gravity effects. It is interesting to note that the stability diagram in 9(b) is consistent with the results of Kumagai et al. [START_REF] Kumagai | Bubble generation by a cylinder moving beneath a free surface[END_REF] for a horizontal cylinder translating beneath an air-water interface.

IV. CONCLUSIONS

Experimental results on the dynamics of the flow around a vertical translated cylinder have been presented. These deal with the turbulent regime before the so-called fully ventilated regime when air reaches the bottom of the cylinder. There is a critical velocity for the onset of air-entrainement. Then, three modes of air-entrainment are observed depending on the diameter and the translating velocity. For D 4 cm, air-entrainment only occurs in the cavity downstream of the cylinder. For D 5 cm, air-entrainment first appears in the wake and then also in the cavity when Re or Fr is increased. This study focuses on air-entrainment in the cavity, where forces, flow visualization, and critical velocities for air-entrainment are reported at constant h or h/D. C D show an increase of the drag coefficient with Re followed by a decrease after air-entrainment occurred in the cavity. A similar evolution is observed for drag-force measurements versus Fr. Moreover, for Fr > 1.2, the drags for all diameters overlap and vary linearly with Fr. The threshold velocity for air-entrainment in the cavity scales as U c ∝ D 0.2 , i.e., Re c ∝ Fr -4 c , which is in good agreement with the scaling of the model of Benusiglio [START_REF] Benusiglio | Indiscrétions aux interfaces[END_REF]. In the future, the shape of the cavity will be studied in more detail in order to quantify the contribution of wetting properties [START_REF] Aristoff | Water entry of small hydrophobic spheres[END_REF][START_REF] Truscott | Water entry of projectiles[END_REF] and the capillary-gravity waves.

FIG. 1 .

 1 FIG. 1. Experimental setup. (a) Schematic tilted view of the flume (drawn up to scale) with the following dimensions: 34 × 1.2 × 0.9 m. (b) Schematic of the arrangement around the cylinder, view from the side abovesurface. The shaded red area represents the laser sheet used for the optical method. (c) Typical picture side view. The cylinder is 2.5 cm in diameter and has been translated from left to right with a velocity of 45.8 cm/s, Re = 11 500, and Fr = 0.9. The depth of the cavity L is represented in the near cylinder wake.

FIG. 2 .

 2 FIG. 2. Free-surface height reconstruction for Re = 8500. (a) D = 3 cm and Fr = 0.5. (b) D = 1.4 cm and Fr = 1.6. The cylinders are indicated by black circles and the white bands correspond to shadow areas of the cylinders that cannot be analyzed.

4 FIG. 3 .

 43 FIG. 3. Drag coefficients, C D , as a function of the Reynolds number for different cylinder diameters D. (a) Constant immersion depth h = 23 ± 0.1 cm and (b) constant ratio h/D = 2.55 ± 0.05. Solid symbols represent the situation without air-entrainment and empty symbols represent flow with air-entrainment behind the cylinder. Additional results from the literature are added and the results without free-surface are highlighted using a gray band.

4 FIG. 4 .

 44 FIG. 4. Drag coefficients, C D , as a function of the Froude number, Fr, for different cylinder diameters D. Panel (a) corresponds to a constant immersion depth of h = 23 ± 0.1 cm and panel (b) corresponds to a constant ratio of h/D = 2.55 ± 0.05. Solid symbols represent the situation without air-entrainment and empty symbols represent flow with air-entrainment in the cavity. For Fr > 1.2, C D is linear and the black line is a linear fit.

FIG. 5 .

 5 FIG. 5. View from the side below-surface of the cylinder of 5 cm in diameter. (a) No air-entrainment at Re = 34 670 or Fr = 0.97, see supplementary movie 1 [38]. (b) Air-entrainment in the wake at Re = 39 700 or Fr = 1.11; bubbles encircled in blue, see Supplemental Material, movie 2 [38]. (c) Air-entrainment in the cavity behind the cylinder at Re = 45 680 or Fr = 1.3; bubbles encircled in red, see Supplemental Material, movie 3 [38]. The dashed lines represent bubble trajectories.

FIG. 7 .

 7 FIG. 7. Flow diagrams Re versus Fr depicting the different air-entrainment regimes: no air-entrainment (solid symbols), air-entrainment in the cavity (partially solid symbols), and air-entrainment in both the cavity and the wake (empty symbols) for different cylinder diameters or Re/Fr. (a) For constant immersion h = 23 ± 0.1 cm and (b) for constant h/D = 2.55 ± 0.05. The lines are guides to the eyes delimiting the different air-entrainment regimes.

FIG. 8 .

 8 FIG. 8. (a) Nondimensional cavity depth L/D as a function of Fr 2 Re, for different diameters. The inset (b) is a zoom restricted to the regime of no or weak air-entrainment, where the continuous lines are linear fits. Solid symbols represent the situation without air-entrainment and empty symbols represent flow with air-entrainment in the cavity. Panels (c) and (d) are snapshots of the cavity, view from the back, typical of no-air entrainment and air-entrainment, respectively. The continuous vertical red lines represent the emerged cylinder walls and the dashed blue line represents the free-surface at rest.

FIG. 9 .

 9 FIG. 9. (a) Critical velocity U c for air-entrainment in the cavity L as a function of the cylinder diameter D for constant immersion h = 23 ± 0.1 cm and constant h/D 2.55 ± 0.05. Results of Benusiglio [20] are plotted with blue crosses. The green dotted line represents the power law D 0.22±0.02 for the constant h/D case and the red solid line represents the power law D 0.20±0.02 for h = 23 ± 0.1 cm. (b) Critical Reynolds number, Re c , as a function of the critical Froude number, Fr c . The black dashed-dotted line represents Re/Fr 25 000. The shaded region corresponds to no air-entrainment. The inset (c) is the same data in log-log scale.

TABLE I .

 I Range of Re and Fr numbers for constant h and constant h/D runs.

			h = 23 ± 0.1 cm				h/D = 2.55 ± 0.05
	D (cm)	h/D	Re min	Re max	Fr min	Fr max	h (cm)	Re min	Re max	Fr min	Fr max
	1.4	16.50	7 470	12 460	1.45	2.41	3.6	4 340	12 160	0.87	2.41
	2.5	9.20	12 460	24 910	1.00	2.10	6.4	9 730	21 880	0.81	1.82
	3	7.77	17 440	29 900	1.08	1.84	7.6	15 950	27 900	0.98	1.72
	4	5.75	24 900	39 860	1.00	1.60	10.2	19 930	39 860	0.80	1.60
	5	4.60	30 630	59 800	0.86	1.71	12.8	28 460	54 550	0.86	1.64
	6	3.83	40 840	71 460	0.87	1.52	15.3	29 180	68 080	0.65	1.52
	7	3.28	29 900	97 260	0.52	1.74	17.9	32 820	97 260	0.58	1.72
	8	2.87	34 040	122 510	0.50	1.70	20.4	34 040	116 710	0.49	1.69
	9	2.55	34 880	119 600	0.41	1.42	23.0	34 880	119 600	0.41	1.42
	11	2.09	39 860	188 870	0.35	1.62	28.1	39 860	184 370	0.35	1.62
	16	1.44	40 860	239 190	0.20	1.20	39.9	39 860	239 190	0.20	1.20

ACKNOWLEDGMENTS

The authors acknowledge the financial support of the Normandy Council and the European Union ERDF funding's through the NEPTUNE project. We thank Antoine Bonnesoeur and Claude Houssin for his engineering support and Paul François for his help with the experiments. Our work has also benefited from exchanges with Yuji Tasaka (Hokkaido University) and Jimmy Philip (University of Melbourne).


Ángel M Arévalo-López

Elena Solana-Madruga

Eugenia P Arévalo-López

Dmitry Khalyavin

Michal Kepa

Antonio J Dos Santos-García

Regino Sáez-Puche

J Paul Attfield

Evolving spin periodicity and lock-in transition in the frustrated ordered ilmenite-type β-Mn 2 InSbO 6

Polar magnets are promising materials for new applications as multiferroics or in spintronics. In double corundum-related oxides, the cation ordering imposes a polar structure and the use of high pressure facilitates the insertion of magnetic cations into the compounds. Here we present the high-pressure synthesis of a new polar and ferrimagnetic corundum derivative Mn 2 InSbO 6 , which adopts the ordered-ilmenite-type structure. Neutron powder diffraction reveals that the high-temperature nearly collinear ferrimagnetic phase evolves to an incommensurate helical structure with k δ = [0 0 k z ] propagation vector, which then locks to the commensurate value of k z = 1/8. This complex magnetic behavior is likely to be related to magnetic frustration and the polar nature of the ordered double corundum structure.

I. INTRODUCTION

High pressure and high temperature treatments of ABO 3 oxides promote unusual cationic arrangements and interesting physical properties. [1] For instance, at ambient pressure, both MnXO 3 with X = V and Ti crystallise in the ilmenite-type structure. They consist in hexagonal closed-packed oxygens with two-thirds of the octahedral holes filled with Mn and V or Ti in an ordered-layer fashion along the zdirection with a non-polar R-3 space group. Under high pressure conditions, MnVO 3 -II shows a perovskite structure with itinerant 3d 1 V 4+ and localised 3d 5 Mn 2+ electrons, where the spins order in a helical magnetic structure. [2] MnTiO 3 -II presents a LiNbO 3 -type, where the cations order along the [111] direction of the cubic cell with the polar R3 space group. It presents multiferroic properties due to Dzyaloshinskii-Moriya exchange interactions (DM). [ 3 ] Following explorations of the double perovskites family, with more than 600 members, [4] double corundum derivatives A 2 BB'O 6 have demonstrated their potential for new multifunctional materials with magnetic ions into both A-and Bsites of Ni 3 TeO 6 -and ordered-ilmenite structural types with R3 polar space group. [5] However, despite the possible combinations, only 15 polar A 2 BB'O 6 double corundum-related compounds have been reported to date. Of these, three are theoretical predictions [ 6 , 7 ] , four synthesised at ambient pressures [8,9,10] and the rest can only be stabilised under high pressure and temperature conditions [11,12,13,14,15,16,17,18,19] (Supplementary Table S2).

The polar lattice of double corundum-type oxides can also stimulate the formation of modulated spin structures via a mechanism described by Lifshitz invariants in a Landau free-energy expansion. [20,21] This gradient energy term favour inhomogeneous magnetic states and is known to result in incommensurate magnetic ground states in MnSi [22] and BiFeO 3 [ 23 , 24 ] . Both of these compounds demonstrate rich phase diagrams including magnetic field induced skyrmion lattice [ 25 ] and hightemperature multiferroic properties respectively [26] .

In this work, we present the high pressure synthesis of a new polar and magnetic double corundum oxide Mn 2 InSbO 6 . The compound has been found to exhibit a remarkable evolution between complex magnetic structures towards a long-period modulated ground state compatible with the mechanism promoted by the gradient Lifshitz terms.

II. EXPERIMENTAL

Mn 2 InSbO 6 sample was synthesised under high pressure-high temperature conditions. The stoichiometric proportions of Mn 2 O 3 , In 2 O 3 and Sb 2 O 3 oxides were ground under acetone to prepare the precursor mixture and treated at 5 GPa and 1373 K during 20 minutes using a Walker-type multianvil apparatus. After the reaction time, the temperature was quenched and the pressure was slowly released. High resolution powder XRD data were collected for Mn 2 InSbO 6 at 80 K, 40 K, 30 K, 20 K, 17 K, 14 K, range with a 0.03º step size using λ = 0.467173 Å synchrotron radiation at the BL04-MSPD beamline at ALBA light source (Barcelona, Spain). For the refinement of the SXRD data, the background was modelled with a list of experimental equidistant points afterwards refined through the linear interpolation of intermediate values and the peak shape was fitted using a Thompson-Cox-Hastings type function. Neutron powder diffraction (NPD) profiles were collected for Mn 2 InSbO 6 at 2 K, 4 K, 6 K, 8 K, 10 K, 14 K, 17 K, 20 K, 30 K, 40 K and 80 K using the WISH instrument (ISIS facility, Oxford, U.K). The crystal and magnetic structures were refined using the FullProf suite [ 27 ] . The nuclear scattering was modeled in the R3 space group and the magnetic intensities were introduced in the refinement procedure using the propagation vectors formalism with k 1 = [0 0 0] (20 K), k 2 = [0 0 1/8] (2 K) and k δ = [0 0 k z ] (intermediate temperatures). Magnetic symmetry analysis was carried out using the program BasIreps [27] for R3 space group. Magnetic susceptibility measurements were performed on a Quantum Design XL-MPMS SQUID magnetometer under DC Zero-Field-Cooling (ZFC) and Field-Cooling (FC) conditions in the temperature range 2 K < T < 100 K under magnetic fields of 0.5 T and 0.01 T. Magnetization hysteresis loops were measured at 2 K, 10 K, 17 K, 20 K, 40 K and 60 K up to 3 T. Specific heat was measured using a Quantum Design PPMS under magnetic fields of 0 T, 3.5 T and 7 T. The dielectric response was measured on cooling with a Andeen-Hagerling Capacitance bridge at 1 kHz in a closed cycle refrigerator.

III. RESULTS AND DISCUSSION

A. Crystal Structure Mn 2 InSbO 6 prepared at ambient pressure crystallises in the Mg 3 TeO 6 -type structure, another corundum type derivative with R-3 centrosymmetric group (a = 8.93 Å and c = 10.71 Å) and all cations in octahedral coordination but with manganese and indium fully disordered. [ 28 ] The polycrystalline β-Mn 2 InSbO 6 polymorph reported here requires high-pressure synthesis at 5 GPa and 1373 K (see Experimental Section). Initial laboratory X-ray diffraction indicated a rhombohedral cell (a ≈ 5.38 Å, c ≈ 14.21 Å) and a small amount of unreacted In 2 O 3 (<5 % wt). Rietveld analysis of combined high-resolution synchrotron X-ray (SXRD) and neutron powder (NPD) diffraction data collected at 80 K (Supplementary Fig. S1) revealed that Mn 2 InSbO 6 crystallises in a doublecorundum-type derivative. It presents alternating xy honeycomb layers formed by ordered Sb / Mn and disordered In / Mn octahedra sharing edges, stacked along the z axis. This structure is known as Ordered-Ilmenite-type (OIL-type). [ 29 ] The contrast between Mn, In and Sb arising from their neutron scattering lengths (Mn = -3.73 fm, In = 4.06 fm, Sb = 5.57 fm) and X-ray form factors (Z = 25, 49, 51) allowed their antisite disorder to be quantified from combined NPD and SXRD refinement. One cation layer is fully ordered with Sb and Mn (site M2) and one is fully disordered between Mn and In (M1 and M3). The crystal structure is shown in Figure 1a and the final refinement results at 80 K are listed in Supplementary Table S3. The antisite exchange between M1 and M3 within the layer results in a larger Mn occupancy at both M1 (51.2 %) and M3 (55.5 %) sites in good agreement with the In deficiency matching the minor In 2 O 3 secondary phase. For M2 and Sb sites, refinements suggested no disorder (antisite allowance refined to <1% with no fit improvement), so hereafter we will refer to the compound as (Mn 0.5 In 0.5 ) 2 MnSbO 6 . This selective site exchange reflects the small size and charge mismatch between In 3+ (r VI = 0.8 Å) and Mn 2+ (r VI = 0.83 Å) compared to that between Sb 5+ (r VI = 0.6 Å) [ 30] and Mn 2+ . Bond valence sums (BVS) calculations with an interpolation method [31] result in +2.34, +2.39, +2.09, +5.10, for M1, M3, M2 and Sb in agreement with formal cationic oxidation states and anti-site disorder of (Mn 0.5 2+ In 0.5 3+ ) 2 Mn 2+ Sb 5+ O 6 . The proposed formal cation oxidation states are supported by X-ray photoelectron spectroscopy (XPS) measurements (Supplementary Table S1 and Fig. S4) [32,33,34,35,36,[START_REF] Hesse | Unifit Line Positions and Data Formats -Version[END_REF] . Bond distances and angles are registered in Supplementary Table S8.

B. Magnetic Properties and

Magnetic Structure The combination of magnetic cations in M1, M2 and M3 sites and anti-site disorder for M1 and M3 gives rise to a complex magnetic behavior. Temperature dependent magnetisation measurements for (Mn 0.5 In 0.5 ) 2 MnSbO 6 in 0.01 T and 0.5 T fields (Fig. 1b) show a ferrimagnetic transition at T N1 = 38(1) K and Curie-Weiss paramagnetism at higher temperatures with a paramagnetic moment of 6.1(1) µ B / Mn, close to the theoretical value of 5.92 µ B for localised S = 5/2 Mn 2+ , and a Weiss temperature of θ = -160(1) K showing that antiferromagnetic spinspin interactions are dominant but partially frustrated (|θ| / T N ≈ 4). A further transition is evidenced in the 0.01 T magnetic susceptibility curves at T N2 = 8(1) K, where zero field and field cooled data diverge from each other. Magnetic hysteresis loops and heat capacity measurements show also transitions at T N1 and T N2 , see Supporting Information.

NPD data were analysed to understand the lowtemperature magnetic behavior in (Mn 0.5 In 0.5 ) 2 MnSbO 6 . On cooling below T N1 = 38 K, magnetic diffraction intensity is observed, for instance on (003) and ( 101 commensurate to commensurate (k 1 to k 2 ) magnetic transition through a continuum of helical magnetic structures. Rietveld fits of both commensurate magnetic structures are shown on difference patterns of 20 K -40 K (k 1 ) and 2 K -40 K (k 1 + k 2 ) in Figure 2b left and the crystallographic details are summarised in Tables S4 andS5. Symmetry analysis allows a consistent description of the different magnetic structures. For stable refinements, the magnetic moment per Mn ion in the three different Mn sites were constrained to be the same. In the k 1 = [0 0 0] phase, the strongest magnetic interactions are the direct (AFM) M-M interactions between cations located into face-sharing octahedra, i.e. M2-M3 (d ~ 2.9 Å). Therefore, ferromagnetic layers M1/M3 and M2 stack antiferromagnetically with the spins ordered along the [110] direction with a refined magnetic moment of 1.9(1) μ B / Mn 2+ at 20 K (Figure 3a). Antisymmetric DM exchange interactions are allowed in the non-centrosymmetric space group R3. The spontaneous ferromagnetic component in this phase, observed in magnetisation measurements, probably comes from the uncompensated magnetization of the sublattices (ferrimagnetic) and canting of the spins (weak-ferromagnetism). Single crystal analysis will be needed to refine the ordered moments of the three magnetic sublattices independently. However, the d 5 -d 5 interactions within the (00l) layers via superexchange through oxygen at ~ 90 º should be AFM according to Goodenough-Kanamori-Anderson rules. Below T δ = 17 K, the interactions between the magnetically diluted M1/M3 sites strengthens and their honeycomb arrangement originates a strong magnetic frustration. This effect induces another transition, mainly in the M3 (blue) site, with a partial relief of frustration leading to spin reorientation as shown in Figure 3b.

All previously reported Mn 2 BSbO 6 corundum derivatives have only two independent magnetic sites, originating the rigid rotation of the FM layers. The OIL structure of Mn 2 InSbO 6 has three independent magnetic sites, what increases the degrees of freedom and allows for the independent rotation of the spins within the honeycomb layers. Below T N2 = 8 K the magnetic structure locks-in into a commensurate state with k 2 = [0 0 1/8], in accordance with the magnetisation and heat capacity measurements (Supporting Information). The structure describes M1 -M2 -M3 layers of spins with a relative orientation of 0º ~ 170º ~ 120º between their spin directions and with a 15º rotation between slabs (Figure 3b). The refined magnetic moment of 3.2(1) μ B per Mn at 2 K is significantly lower than the ideal value of 5 μ B for S= 5/2 Mn 2+ showing that frustration is still present. Magnetic diffuse scattering is present at 2 K, see Figure S3 in SI. The two magnetic structures coexist below 17 K and 11% of the k 1 phase persists at 2 K. Combined Rietveld refinements of SXRD and NPD patterns collected between 2 K and 40 K (Table S7) reveal magnetostriction in the cell parameters when the spin reorientation takes place, Figure 2b right.

C. Dielectric properties

The OIL-type structure of (Mn 0.5 In 0.5 ) 2 MnSbO 6 with polar R3 space group is expected to present a net spontaneous polarisation. However, no switching is possible since the loss of centre of symmetry comes from cation ordering and not from cooperative ion motion (i.e. the material is an improper ferroelectric). The variation of the permittivity and the dielectric loss (tan δ) thermal responses show an increase below T δ = 17 K, Fig. 4 (top). These variations resemble the magnetostriction observed in the cell parameters as shown in Figure 2b. The observed magnetostriction is coupled to the additional displacements of the cations along the c axis below T N . This ionic motion is expected to originate an increase of the polarization (P s ). The expected P s has been estimated using a point charge model from the experimental atomic displacements. When compared to the predicted value at 80 K (~ 12 µC / cm 2 ), P s shows a negligible change at 40 K, but it gradually increases below T N1 up to ~25% of ΔP s (~ 15 µC / cm 2 ) when it adopts an incommensurable magnetic structure at 14 K. This polarisation is directly related to the expansion (contraction) of M2 and Sb (M1 and M3) octahedra, as shown at the bottom of Figure 4, induced by the magnetic order. The octahedral distortion has been calculated in terms of their mean ellipsoid radii (<R>) and shape (S) parameters from Pieface software [START_REF] Cumby | [END_REF] . The results, shown in Fig. S6, confirm the described cooperative evolution of the octahedral distortions resulting from the magnetic transition at T N1 .

D. Relevant Free-Energy terms

The non-collinear orientation of the three magnetic sublattices in the helical phase (Figure 3b) strongly points to the presence of competing exchange interactions in the system. The inter-sublattice interactions are antiferromagnetic and comparable in magnitude, giving rise to magnetic frustration. Another important mechanism promoting incommensurate spin ordering is related to the polar nature of the crystal structure of Mn 2 InSbO 6 . This implies that the physically irreducible representation GM2GM3(η,η*) which transforms the magnetic order parameter in the high-temperature k=0 phase, allows the presence of gradient Lifshitz invariant in the Landau free-energy decomposition, . The component of the gradient operator is transformed by the totally symmetric representation and the transformational properties of the magnetic order parameter are specified in Table S9 and S10 of the Supplemental information. This term favours an inhomogeneous magnetic state [ 39 , 40 ] propagating along the Λ-line of symmetry (k = 0,0,k z ) and is expected to have a strong impact on the magnetic ground state in systems without strong single ion anisotropy. The mechanism involves relatively weak anti-symmetric DM exchange interactions at a microscopic level. The latter usually implies a longperiod modulation. A nearly quenched orbital angular momentum expected for Mn 2+ cations and the small value of k z found experimentally in (Mn 0.5 In 0.5 ) 2 MnSbO 6 (see Figure 2a) are both consistent with this mechanism.

It is very likely that frustration of the Heisenberg symmetric exchange and the antisymmetric interactions imposed by the polar structural distortions are both engaged into formation of the low-temperature modulated ground state. The former mechanism defines the relative magnetic phases of Mn spins between the three symmetry inequivalent sites M1, M2 and M3 (three sublattices) and the latter promotes the slow global spin rotation along the caxis. A possible reason why the modulated state does not develop below the first transition at T N1 = 38 K is a partial lifting of the frustration via adopting different magnitudes of the ordered moments for the distinct sublattices. In other words, the system can control the level of the frustration by varying the ordered moments of the sublattices, and fixing them at certain ratios. This degree of freedom, affecting both symmetric and anti-symmetric exchange, is available only at temperatures close to the transition, where the spins are far from saturation, and it is expected to vanish at lower temperatures. Note, this assumption can explain the magnetisation observed in M vs H hysteresis loops and the plateau in the temperature dependence of the ordered moment in the 15 < T < 30 K range (Fig. 3c). However, single crystal neutron diffraction will be needed in order to refine the magnetic moments of the three sublattices independently.

The lock-in transition with the commensurate value k z = 1/8 implies an activation of the free-energy invariant that selects the commensurate value for the propagation vector [START_REF] Saxena | [END_REF] .

In the case of (Mn 0.5 In 0.5 ) 2 MnSbO 6 , analysis of the allowed freeenergy terms revealed that the lock-in invariant indeed exists and it involves the order parameter components in an unprecedented 24 th power, ρ 24 +ρ * 24 +ξ 24 +ξ * 24 . Here, ρ,ρ* and ξ,ξ* are components of the complex magnetic order parameters of the modulated ground state, describing helices with opposite chirality (see Tables S9 andS10 in the Supplemental Information). This is a remarkably high-power term and to the best of our knowledge, (Mn 0.5 In 0.5 ) 2 MnSbO 6 is a unique example where lock-in takes place at such a small value of the propagation vector and implies a periodicity of ~114 Å (note that k 2 = [0 0 1/8] in the hexagonal setting of the rhombohedral cell corresponds to k 2 = [1/24 1/24 1/24], in the primitive setting).

IV. CONCLUSIONS (Mn 0.5 In 0.5 ) 2 MnSbO 6 OIL polymorph was synthesised under high pressure and high temperature conditions, stabilising a partial Mn/In cation order. This polymorph exhibits a sequence of magnetic phase transitions at 38 K, 17 K and 8 K. The first one gives rise to a commensurate magnetic structure with propagation vector k 1 = [0 0 0], where the Mn 2+ spins are arranged into honeycomb (M1 and M3 sites) and triangular (M2) sublattices with a net ferromagnetic component. The second transition results in an incommensurately modulated state with k = [0 0 k z ], which locks into the commensurate value k z = 1/8, below 8 K. The complex magnetic behavior is likely to be associated with magnetic frustration of the symmetric Heisenberg exchange and the polar nature of the crystal structure promoting DM interactions. The role of the anti-symmetric interactions in formation of the modulated ground state of (Mn 0.5 In 0.5 ) 2 MnSbO 6 is likely to be similar to that in other notable non-centrosymmetric magnets like skyrmion-forming MnSi and multiferroic BiFeO 3 .

Figure 1 .

 1 Figure 1. a) Crystal structure of (Mn 0.5 In 0.5 ) 2 MnSbO 6 along [110] direction. Grey and blue spheres represent the ordered Sb and Mn sites respectively. Blue-yellow spheres show the mixed Mn-In sites. b) Thermal evolution of the direct and inverse magnetic susceptibilities measured at 0.5 and 0.01 T in FC and ZFC modes. Dashed line represents the Curie-Weiss fit.

  ) peaks at ~ 4.7 Å and ~ 4.4 Å as shown in Figure 2a. All the magnetic diffraction peaks observed down to 17 K can be indexed with propagation vector k 1 = [0 0 0]. Below T δ = 17 K the intensity of the magnetic diffraction peaks decreases and new satellite peaks appear. Their position gradually changes, revealing a second transition with k δ = [0 0 k z ] that locks-in at a third transition with k z = 1/8 at T N2 = 8 K. The thermal evolution of k z is shown on the right panel of Figure 2a. k z increases from 0 to 1/8 on cooling and it evidences the gradual

Figure 2 .

 2 Figure 2. a) 2D contour plot showing the temperature evolution of the (003) and (101) magnetic diffraction maxima. The green lines at 38 K, 17 K and 8 K mark the transitions between the different magnetic structures. The panel on the right shows the change on the z component of the propagation vector k δ = [0 0 k z ]. b) Magnetic structure fits to the difference between 20 K and 40 K (top) and 2 K and 40 K (bottom) WISH NPD profiles of (Mn 0.5 In 0.5 ) 2 MnSbO 6 . The vertical ticks indicate the positions of the Bragg reflections for k 1 = [0 0 0] (20 K) and k 2 = [0 0 1/8] (2 K) (lower row) models. The panel on the right shows the thermal evolution of the cell parameters, revealing magnetostriction induced by the evolving magnetic structure.

Figure 3 .

 3 Figure 3. [001] and [110] projections of the magnetic structures of (Mn 0.5 In 0.5 ) 2 MnSbO 6 . Yellow, red and blue arrows represent the Mn 2+ spins at M1, M2 and M3 sites respectively (Y, R, B). a) k 1 = [0 0 0] at 20 K. b) Half-cell of the k 2 = [0 0 1/8] structure showing the frustration relieved by B rotation; there is a 15º rotation between Y-R-B spins slabs. c) Thermal evolution of the magnetic

Figure 4 .

 4 Figure 4. Thermal evolution of the dielectric response of (Mn 0.5 In 0.5 ) 2 MnSbO 6 (top) showing a transition in the dielectric permittivity below 8 K. The variation in the polarization (middle), determined from the experimental

ACKNOWLEDGEMENTS

The authors thank EPSRC, STFC and the Royal Society for support. MINECO and Comunidad de Madrid are also acknowledged for funding through projects MAT2017-84385-R, and S-2013/MIT-1275 respectively. A.M.A.L. thanks a Marie-Słodowska Individual Fellowship in the Horizon 2020 research and innovation programme. We thank ISIS and ALBA-Synchrotron for awarding beam-time in WISH and BL04-MSPD beam-lines respectively. Thanks to Dr. O. Vallcorba for help on SXRD experiments and to L. Huerta Arcos for help on the XPS experiments.


Francois Bondu

Hen Goulc'

Ludovic Loas

Cyril Frein

Anthony Hamel

Guillaume Carré

Jean-Francois Ducournau

Virginie Lampin

Francoise Daburon

Binet

An optical instrument for the submillimeter spectroscopy of the volatile metabolome

Abstract-We are designing and realizing an instrument for the spectroscopy of the volatile metabolome of the soil microorganisms. The use of optical telecommunication elements makes it possible to realize a compact source with a wide range of accessible frequencies (100 GHz -1 THz range; 3-33 cm⁻ ¹) with an accuracy of 10⁵ ⁻⁶ ⁻ Hz (3.10 cm ¹), as well as a low bandwidth detection system. The instrument will discriminate among complex polar molecules. We report on the preliminary performances of parts of the source system.

I. INTRODUCTION

ACKLING the complex ecological situations of living bodies adequate tools. This complexity has multiple aspects. No vegetal or animal living organism is on its own: it always needs the contribution of diverse communities of microorganisms. Most species of microorganisms can not be grown alone in Petri dishes, making it hard to pinpoint them. Scale laws in ecology suggest that the ratio of known species of microorganisms to unknown ones might be a ratio of the order of 10⁻ ⁵ [START_REF] Locey | Scaling laws predict global microbial diversity[END_REF]. There is a difficulty to match the results on microorganism population presence given with the genomics techniques compared to the techniques that observe living material. The soil itself is a large reservoir of microorganisms, with a number of 10 to 10 per gram.

⁶ ⁶ How the bacterial diversity and function in soil are linked to biotic control is an open issue [START_REF] Monard | Relationship between bacterial diversity and function under biotic control: the soil pesticide degraders as a case study[END_REF]. The use of new instruments targeting the volatile metabolomes of microorganisms may contribute to the global comprehension of the possible equilibrium or dynamical evolution of complex communities, a goal that looks way ahead today.

The submillimeter spectroscopy has been historically developed for astronomy purposes [START_REF] De | The submillimeter: A spectroscopist's view[END_REF], leading to the identification of molecules as complex as alanine in the inter stellar medium. This technique has not been used, up to now, for the analysis of the metabolism of microorganisms. One possible reason might be the complexity of operation and space requirements of the existing apparatus. Even in the lab, the change from one frequency band to another one, with various multipliers, amplifiers, connectors, waveguides, is tedious. However, the submillimeter spectroscopy would present advantageous characters. Contrary to the chromatography -mass spectrometry instruments and similar apparatus, it does not require heat or high pressures, thus promises a better integrity of the complex biological molecules. Even more, at low pressure and at room temperature, the absorption linewidths reduce to the range of some hundreds of kHz, thus simultaneously making a clear discrimination between molecules types and conformations and improving the signal to noise ratio.

The analysis of the metabolism of a variety of microorganisms requires sensitivity to a wide range of organic molecules. The access to a wide range of attainable frequencies seems thus to be a desirable target: for one given organic molecule, public tables (HITRAN, JPL, CDMS) indicate that the absorption peak that gives the maximum sensitivity can not be constrained to a restricted frequency band. For that prospect, the techniques based on optical photomixing [START_REF] Hindle | Widely tunable THz synthesizer[END_REF] seem to be more promising than the techniques based on up-conversion of microwave signals [START_REF] Neese | Compact Submillimeter/Terahertz Gas Sensor With Efficient Gas Collection, Preconcentration, and ppt Sensitivity[END_REF].

We aim to build a compact instrument based on optical elements designed for the optical telecommunications. A challenge with optical techniques is that the accessible range for the frequency of a laser is often a compromise between the linewidth and a robust and repeatable operation. We demonstrate here that these key elements are simultaneously under control, still with simple and compact elements.

II. RESULTS

The submillimeter wave synthesizer source is based on three commercial fibered distributed feedback laser diodes (DFB) at a wavelength of 1.5 μm, cf. figure 1. The frequency difference of lasers 1 and 2 is adjustable between DC and 753 GHz with constant laser output powers, with the change of the temperature of the free-running DFB chips. The laser lines are locked to two resonances using the Pound-Drever-Hall technique [START_REF] Hallal | Synthesis of a 30-Hz Linewidth Wave Tunable Over 500 GHz[END_REF], with 10's of MHz bandwidth home-made electronics. The frequency difference is then a multiple of 1 GHz, the free spectral range of the commercial optical resonator. The frequency difference of lasers 2 and 3 can be swept continuously in a 5-10 GHz range, when this difference is controlled with a phase locked loop. Ultimately, the submillimeter synthesizer frequency will be the difference between lasers 1 and 3, continuously tunable by contiguous segments between DC and 763 GHz. A "UTC photodiode" and associated antenna, designed and realized for the project, have a bandwidth of several hundreds of GHz.

François Bondu 1 , Goulc'hen Loas 1 , Ludovic Frein 1 , Cyril Hamel 1 , Anthony Carré 1 , Guillaume Ducournau 2 , Jean-François Lampin 2 , Virginie Daburon 3 , Françoise Binet Up to now, we have demonstrated experimentally the locking loops between laser 1 and 2 and between laser 2 and 3 separately. Figure 2 shows the performance of the phase noise of the beat notes that are measured by photomixers. The phase noise will be therefore negligible in the final system, as the linewidth of the synthesizer line is of the order of a few 10's of Hz.

The required accuracy, for spectroscopy, is of the order of 10 kHz for the radiated submillimeter wave. It seems to be at reach. The commercial optical resonator has a 1.6 GHz/K sensitivity for one optical carrier. Temperature stabilization at 1/10 th of a degree will be enough to disambiguate the interrogated resonance. The beat note accuracy at the 10 kHz range should be attained with a temperature control of the cavity at a few mK level.

III. NEXT STEPS

We are now aiming at showing the performances of the full optical submillimeter synthesizer, and then at the detection system.

TFig. 1 .

 1 Fig. 1. Principle of the submillimeter spectroscopy instrument. On the left side, the submillimeter wave synthesizer, based on commercial telecommunication optical elements, at a wavelength of 1.5 µm. A dedicated UTC photodiode and antenna convert the optical radiation into a submillimeter wave in free space.

Fig. 2 .

 2 Fig. 2. Performance of the phase noise of two subsets of the frequency synthesizer source.

  3 1 Institut FOTON -UMR CNRS 6082, Univ Rennes, CNRS, , Rennes, France 2 IEMN -UMR CNRS 8520, Université de Lille 1, Villeneuve d'Ascq, France3 Ecobio -UMR 6553, Univ Rennes, CNRS, FranceAn optical instrument for the submillimeter spectroscopy of the volatile metabolome


Tanguy Choupin

Bruno Fayolle
email: bruno.fayolle@ensam.eu

Gilles Régnier

Christophe Paris

Jacques Cinquin

Benoît Brulé

A more reliable DSC-based methodology to study crystallization kinetics: Application to poly(ether ketone ketone) (PEKK) copolymers

Keywords:

Method is based on derivative expression for crystallization Hillier model.

• Model fitting directly on the crystallization heat flux measured by DSC.

• The method avoids the extrapolation of the crystallization peak beginning.

• PEKK Crystallization kinetic modeling including primary and secondary crystallization.

Introduction

High-performance poly(ether ketone ketone) (PEKK) thermoplastics are currently studied with great interest by the aeronautic industry as a matrix for carbon-fibre reinforced structural parts. In fact, PEKK composites can be consolidated out of autoclave, while they have lower processing temperatures than poly(ether ether ketone) (PEEK) composites. The final mechanical performances of composite parts are closely linked to their processing. Depending on the thermal processing cycle, the matrix properties as well as the interaction between the matrix and carbon fibres are modified. During the manufacturing of PEKK composite parts, the composite is first heated above the matrix melting temperature to decrease the viscosity of the matrix in order to evacuate all porosity under pressure. This step is known as consolidation. The composite part is then cooled to room temperature. During this step, the matrix may crystallize with a different crystallinity ratio depending on the cooling conditions that directly impact the final mechanical performances of the composite part [START_REF] Talbott | The effects of crystallinity on the mechanical properties of PEEK polymer and graphite fiber reinforced PEEK[END_REF][START_REF] Cebe | Effect of thermal history on mechanical properties of polyetheretherketone below the glass transition temperature[END_REF][START_REF] Peacock | The mechanical properties of random copolymers of ethylene: force-elongation relations[END_REF][START_REF] Gao | Cooling rate influences in carbon fibre/PEEK composites. Part 1. Crystallinity and interface adhesion[END_REF][START_REF] Diani | On the relevance of the micromechanics approach for predicting the linear viscoelastic behavior of semi-crystalline poly(ethylene) terephtalates (PET)[END_REF][START_REF] Pawlak | Plastic deformation of crystalline polymers: the role of cavitation and crystal plasticity[END_REF]. For this reason, PEKK crystallization kinetics should be controlled in order to optimize the processing parameters in terms of the mechanical performances and processing times of the composite parts.

The crystallization mechanisms of PEKK copolymers have been studied by different authors [START_REF] Gardner | Structure, crystallization and morphology of poly (aryl ether ketone ketone)[END_REF][START_REF] Hsiao | Isothermal crystallization kinetics of poly(ether ketone ketone) and its carbon-fibre-reinforced composites[END_REF][START_REF] Quiroga Cortés | Morphology and dynamical mechanical properties of poly ether ketone ketone (PEKK) with meta phenyl links[END_REF][START_REF] Choupin | Isothermal crystallization kinetic modeling of poly(etherketoneketone) (PEKK) copolymer[END_REF] who have shown that similarly to PEEK, PEKK copolymers crystallize with two different crystallization processes, namely a primary crystallization that appears at the early stage of crystallization and a secondary crystallization corresponding to the enhancement of crystallization with the growth of a crystalline structure within lamellae [START_REF] Gardner | Structure, crystallization and morphology of poly (aryl ether ketone ketone)[END_REF]. This subsidiary crystallization should be taken into account in the crystallization kinetic model. Among the different crystallization kinetic models that consider this two-step crystallization, it was previously argued [START_REF] Choupin | Isothermal crystallization kinetic modeling of poly(etherketoneketone) (PEKK) copolymer[END_REF] that Hillier [START_REF] Hillier | Modified avrami equation for the bulk crystallization kinetics of spherulitic polymers[END_REF] proposed a good model with secondary crystallization kinetics directly dependant on the primary crystallized volume. Hsiao et al. [START_REF] Hsiao | Isothermal crystallization kinetics of poly(ether ketone ketone) and its carbon-fibre-reinforced composites[END_REF] proposed an extension of this model by considering that the Avrami exponent of secondary crystallization n 2 could have values other than 1; the total relative crystallinity α(t) thus evolves as:

= + t w t w t ( ) ( ) ( ) 1 1 2 2 (1)
with the respective primary and secondary crystallization expressions:

= t K t ( ) 1 exp ( ) n 1 1 1 (2) = [ ( ) ] [ ] t K n K t K t d ( ) 1 exp ( ) exp ( ) t n n n 2 2 2 0 1 1 2 1 2 2 (3)
where n 1 and n 2 correspond to the Avrami exponents, K 1 and K 2 to the crystallization rate constants, and w 1 and w 2 to the weight factor for the primary and the secondary crystallization, respectively. These weight factor values provide the relative importance between both crystallizations. At any time, we have w 1 +w 2 = 1. For an infinite time, the total relative crystallinity is the unity that corresponds to the Avrami model hypotheses [START_REF] Avrami | Kinetics of phase change. I general theory[END_REF].

To identify the crystallization kinetic parameters, the Hillier model has to be fitted to the relative crystallinity calculated from the integration of crystallization heat flows measured using differential scanning calorimetry (DSC):

= t Q t Q t dt ( ) ( ) ( ) t t 0 0 (4)
where Q(t) is the heat flow measured at time t, and t ∞ is the time when the polymer is fully crystallized.

As shown in the literature for other polymers [START_REF] Hsiao | Isothermal crystallization kinetics of poly(ether ketone ketone) and its carbon-fibre-reinforced composites[END_REF][START_REF] Jonas | Thermal stability and crystallization of poly(aryl ether ether ketone)[END_REF][START_REF] Lorenzo | DSC isothermal polymer crystallization kinetics measurements and the use of the Avrami equation to fit the data: guidelines to avoid common problems[END_REF][START_REF] Tardif | Experimental study of crystallization of polyetheretherketone (PEEK) over a large temperature range using a nano-calorimeter[END_REF], the integration of the PEKK crystallization peak is rendered difficult due to the truncation of the crystallization peak beginning, which has been attributed to instrumental thermal inertia when switching from the cooling to the isothermal step (Fig. 1). This phenomenon is more significant when the isothermal crystallization temperature decreases or the polymer has faster crystallization kinetics.

In a previous article [START_REF] Choupin | Isothermal crystallization kinetic modeling of poly(etherketoneketone) (PEKK) copolymer[END_REF], we applied the usual method by extrapolating the beginning of the crystallization peak with a straight line (dashed line in Fig. 1) until the base line. However, this method implied the introduction of a crystallization induction time, which is not consistent with overall crystallization kinetic theories. Secondly, the Avrami exponent in this method was found to be close to 2, which is associated with two-dimensional crystallization growth according to the Avrami theory if we consider instantaneous nucleation. However, this value of 2 is not consistent with the spherulite entities detected by microscopic observations as reported in the literature [START_REF] Gardner | Structure, crystallization and morphology of poly (aryl ether ketone ketone)[END_REF][START_REF] Choupin | Isothermal crystallization kinetic modeling of poly(etherketoneketone) (PEKK) copolymer[END_REF].

This paper aimed to propose a new methodology to quantify the isothermal crystallization kinetics of PEKK without extrapolating the crystallization peak beginning. We show that this extrapolation could lead to physically inconsistent values for the Avrami exponent in terms of crystalline morphologies. Indeed, microscopic observations of quiescent crystallization put in evidence an instantaneous nucleation and a spherulitic growth corresponding to an Avrami exponent equal to 3. This new method will be applied to two PEKK copolymers with different crystallization kinetics, namely PEKK 6002 and 7002, to identify their crystallization kinetic parameters. Let us recall here that the proposed kinetic modeling can be applied only for some processing conditions in which strain (elongation or shear) rate during the crystallization can be neglected such as additive manufacturing or welding processes. Indeed, high strain rate conditions mainly promote the nucleation mechanisms by increasing nuclei number [START_REF] Bustos | Effect of molecular architecture on quiescent and shear-induced crystallization of polyethylene[END_REF]. Finally, TTT diagrams of the relative volume crystallinity will be built from the overall crystallization kinetic modeling.

Experimental

Materials

The PEKK homopolymer is prepared using diphenyl ether (DPE) and terephthalic acid (T) (Fig. 2a). Terephthalic acid with para phenyl links can be substituted by isophthalic acid (I) with meta phenyl links (Fig. 2b) and form copolymers comprising two different isomers, namely terephthaloyl (Fig. 2a) and isophthaloyl isomers (Fig. 2b). These isomers can form two different diads: TT diads compounded with two terephthaloyl isomers and TI diads compounded with one terephthaloyl isomer and one isophthaloyl isomer.

The PEKK used in this study is provided by Arkema and referenced as PEKK 6002 and 7002. The ratios of terephthalic acid with para phenyl links and isophthalic acid with meta phenyl links are 60/40 and 70/30 for PEKK 6002 and 7002, respectively. The two PEKK grades have an average weight (M w ) and number average molecular mass (M n ) of about 70 kg.mol -1 and 30 kg mol -1 , respectively. Glass transition temperatures of both PEKKs are measured at around 155 °C, while the melting temperatures for PEKK 6002 and 7002 are 305 °C and 332 °C, respectively, using DSC with a heating stage of 10 °C/min.

Differential scanning calorimetry (DSC)

Isothermal crystallization analyses are carried out using TA Instruments Q2000 with granules of about 7-8 mg. All specimens are first heated at 10 °C.min -1 from room temperature to 360 °C for PEKK 6002 and 380 °C for PEKK 7002 for 5 min. These temperatures are above the PEKK equilibrium melting temperatures (T m0 = 340 °C [START_REF] Gardner | Structure, crystallization and morphology of poly (aryl ether ketone ketone)[END_REF]) and erase the thermal history. Specimens are then cooled at 40 °C.min -1 to annealing temperatures. Annealing temperatures are maintained for 240 min, and specimens are finally cooled to room temperature at 40 °C.min -1 . The chosen cooling and heating rates are high enough to ensure that the polymer does not crystallize before annealing, while the annealing time is long enough to allow the polymer to crystallize fully. For all specimens, a heat scan at 10 °C.min -1 to 400 °C is carried out after crystallization to measure the glass transition, melting temperature, and melting enthalpy induced by the crystallization cycle.

Results and discussion

Derivative Hillier method

This new identification method based on the Hillier model consists of directly fitting the crystallization heat flow measured by DSC without extrapolating the heat flux at low times and considering only reliable DSC measurements (Fig. 1). The fitting can be done using multi-nonlinear regression available in various numerical tools; we used the "fit" function of the Scilab ® development tool. Moreover, directly identifying the heat flux by multi-nonlinear regression allows the best baseline to be easily determined according to the model. By deriving Eq. ( 4), it can thus be written:

= = = d dt Q t Q t dt Q t H Q t H ( ) ( ) ( ) ( ) o t c m ( 5
)
where α is the relative crystallinity, Q(t) the DSC heat flux per mass unit, ΔH c the crystallization enthalpy, and ΔH m the melting enthalpy after crystallization, which was obtained by rapid heating to hinder any fusion or recrystallization process. The derivative of the primary crystallization kinetic equation (Eq. ( 2)) can be expressed as follows: = ( )

d dt K n t K t exp n n 1 1 1 1 1 1 1 (6)
It can be observed that at time close to 0, the derivative of the primary crystallization kinetic equation (Eq. ( 6)) tends toward:

d dt K n t n 1 1 1 1 1 (7)
Thus, the choice or imposition of curve type to extrapolate the DSC signal will have an influence on the value of the Avrami exponent. If n 1 is equal to equal to 2, for example, the expression of the derivative of the primary crystallization corresponds to the equation of a linear curve with a slope equal to K 1 n 1 . This observation is very interesting, because if we assume that the primary crystallization is the main physical phenomenon during the beginning of crystallization, then our extrapolation of the DSC signal by a linear curve primarily explains the Avrami exponent value close to 2 for the primary crystallization found in the first approach [START_REF] Choupin | Isothermal crystallization kinetic modeling of poly(etherketoneketone) (PEKK) copolymer[END_REF]. However, this value is inconsistent with microscopic observations previously reported in Refs. [START_REF] Gardner | Structure, crystallization and morphology of poly (aryl ether ketone ketone)[END_REF] and [START_REF] Choupin | Isothermal crystallization kinetic modeling of poly(etherketoneketone) (PEKK) copolymer[END_REF]: it can be observed that all crystallite entities have the same size, so instantaneous crystallization can be easily assumed. Moreover, it is difficult to imagine that crystalline entities are not spherulites, and thus according to the Avrami theory [START_REF] Avrami | Kinetics of phase change. I general theory[END_REF], the primary crystallization Avrami exponents n 1 could be logically fixed to a value of 3.

The secondary crystallization rate has the following form:

= t f t d ( ) ( , ) t 2 0 (8)
Its derivative can be calculated as follows:

= + d dt f t t f t t d ( , ) ( , ) t 2 0 (9)
Consequently, the derivative of the secondary crystallization kinetic equation (Eq. ( 3)) calculated using Eq. ( 9) can be expressed by:

+ K n K t K n K K t n t K n t d [ 1 exp( )] [ 1 exp( )]exp [ ( ) ] { ( 1) ( ) ( ) } d dt n t n n n n 2 2 1 2 2 0 1 2 2 2 2 2 2 2 2 1 1 2 2 2 ( 10
)
The fitting of the derivative Hillier model (Eq. ( 6) and ( 10)) to the heat flow measured by DSC and divided by the melting enthalpy ΔH m (Eq. ( 5)) for the isothermal crystallization of PEKK 6002 at 230 °C and 200 °C is presented in Fig. 3. Heat flux for the primary (Eq. ( 6)) and secondary crystallization mechanisms (Eq. ( 10)) based on the derivative Hillier model fitting are also plotted. Assuming that the crystallization process begins at the same time as the isothermal stage, it can be observed that no induction times are found for this approach compared to our previous work [START_REF] Choupin | Isothermal crystallization kinetic modeling of poly(etherketoneketone) (PEKK) copolymer[END_REF].

It is noteworthy that the best fits shown in Fig. 3 for the derivative Hillier model determined with non-linear regression fitting are obtained with a value of 1 for the secondary crystallization Avrami exponent (n 2 = 1). In this case, Eq. ( 10) can be simplified as follows:

[ ( )] [ ( )] d dt K K t K K K t d 1 exp 1 exp exp [ ( )] n t n 2 2 1 1 2 2 0 1 1 2 (11)
As a result, the derivative of the Hillier crystallization kinetic model can be expressed by:

= + w K n t K t w K K t w K K K t d exp( ) [ 1 exp( )] [ 1 exp ( )]exp [ ( )] d dt n n n t n 1 1 1 1 1 1 1 2 2 1 1 2 2 2 0 1 1 2 (12)
Finally, the simplified expression of the derivative Hillier model is, according to our knowledge, given below for the first time:

= + d dt w K n t w K (1 ) ( ) n 1 1 1 1 1 2 2 1 2 1 ( 13
)
Although this last relation is not useful for fitting, it allows us to understand well the physical meaning of the Hillier model by showing the dependence of the secondary crystallization on the primary crystallization. Nevertheless, this formulation is easily implementable in numerical models.

Finally, we can observe in Fig. 3 that the overall derivative Hillier model fits well with the heat flow and allows us to build the beginning of the crystallization peak truncated by DSC. Besides, we can witnessed the first part of the derivative Hillier model is first associated with the primary crystallization, which shows that the very beginning of the crystallization process corresponds only to the growth of the primary crystallization. This is then followed by the secondary crystallization growth in relation to the enhancement of the first crystallization stage. The crystallization kinetics of PEKK 7002 is also studied to validate the derivative Hillier model on a PEKK copolymer that crystallizes faster than PEKK 6002. Fig. 4 presents the result of the fitting of the derivative Hillier model (Eqs. ( 6) and ( 10)) with the heat flow measured by DSC for PEKK 7002 crystallized at 250 °C compared to PEKK 6002 crystallized at 230 °C. It is noteworthy the crystallization kinetic for PEKK 6002 is not complete after 13 min, especially the secondary one which is more visible in Fig. 3a. The same primary crystallization Avrami exponent (n 1 = 3) as PEKK 6002 is chosen for PEKK 7002 regarding the optical microscopic observations. The same secondary crystallization Avrami exponent (n 2 = 1) as PEKK 6002 is also used corresponding to the best fit of the derivative Hillier model. The chosen crystallization temperatures correspond to the fastest temperatures for the crystallization kinetics of PEKK 6002 and 7002. It can be observed that even if the crystallization peak for PEKK 7002 is very narrow compared to PEKK 6002 (Fig. 4), the derivative Hillier model fits well with the experiments and allows us to rebuild the crystallization peak beginning.

Hillier parameter modeling

Following the general approach reported in the literature for nucleation and crystallization growth, the crystallization rate constants (K 1 and K 2 ) can be expressed as a function of the initial number of potential nuclei (N 0 ) and crystal growth rate (G) depending on the crystallization dimensionality and nucleation mode. For primary crystallization, shown to have three-dimensional crystallization growth with an instantaneous nucleation, the crystallization rate constant can be written as:

= K N G 4 3 1 01 1 3 ( 14
)
where N 01 is the initial number of potential nuclei and G 1 the crystallization growth rate for primary crystallization.

For secondary crystallization, shown to have one-dimensional crystallization growth with instantaneous nucleation, the crystallization rate can be written as:

= K N G 2 02 2 (15)
where N 02 is the initial number of potential nuclei and G 2 the crystallization growth rate for secondary crystallization.

According to the Hoffmann and Lauritzen theory [START_REF] Turnbull | Rate of nucleation in condensed systems[END_REF][START_REF] Hoffman | Crystallization of bulk polymers with chain folding: theory of growth of lamellar spherulites[END_REF][START_REF] Hoffman | Kinetic of crystallization from the melt and chain folding in polyethylene fractions revisited: theory and experiment[END_REF], crystal growth can be expressed by:

= G T G U R T T K T T f ( ) exp ( ) exp i i gi 0 (16)
where G 0i is a pre-exponential factor independent of temperature. The first exponential term contains the contribution of the diffusion process to the growth rate, where U * is the activation energy of the molecular transferring through the melt crystal interface, T ∞ is the temperature below which diffusion stops (T ∞ = T g -30 K) with T g = 155 °C, and R is the gas constant.

The second exponential term is the contribution of the nucleation process, where Kgi is the activation energy of nucleation for a crystal of critical size, ΔT is the degree of supercooling (ΔT = T m0 -T) with T m0 as the equilibrium melting temperature, and f is a correction coefficient for the temperature dependence of the melting enthalpy (f = 2T/ (Tm0+T)). The index i for G i , G 0i , and K gi is equal to 1 for primary crystallization and 2 for secondary crystallization.

Finally, the primary crystallization rate constant can be written as: And the secondary crystallization rate constant can be written as: Regarding the evolution of the primary crystallization weight factor w 1 with the crystallization temperature (Fig. 6), it was chosen to model w 1 by a linear equation:

= K T K U R T T K T T f ( ) exp 3 ( ) exp 3 g 1 01 1 (17)
= K T K U R T T K T T f ( ) exp ( ) exp
= + w T a T b ( ) 1 1 1 (19)
where a 1 and b 1 are two constants independent of temperature.

The comparison of the logarithm of the Hillier parameters K 1 and K 2 for PEKK 6002 and 7002 is plotted in Fig. 5 and fitted with the Hoffman and Lauritzen model (Eqs. ( 16) and ( 17), respectively). Unlike PEKK 6002, it is not possible to carry out isothermal crystallization measurements over a large temperature range, but only until 260 °C. As PEKK 7002 crystallizes very quickly, below this temperature, the polymer crystallizes before the isothermal crystallization stage. A solution could be to use a flash DSC to reach the lower crystallization as performed by Tardif et al. [START_REF] Tardif | Experimental study of crystallization of polyetheretherketone (PEEK) over a large temperature range using a nano-calorimeter[END_REF] for PEEK (although this tool implies very small samples that may not be representative of a voluminous crystallization process). Results for K 0i and K gi are presented in Table 1. As expected, we observed that both plots have typical bell shapes in accordance with the Hoffman and Lauritzen growth theory [START_REF] Hoffman | Crystallization of bulk polymers with chain folding: theory of growth of lamellar spherulites[END_REF]. Crystallization kinetic parameters K 1 and K 2 are much higher for PEKK 7002 than PEKK 6002. Crystallization kinetics are the fastest around 240 °C for PEKK 6002 and 245 °C for PEKK 7002. Values of K gi are similar to those reported in the literature for other poly(aryl ether ketones) [START_REF] Liu | Isothermal melt and cold crystallization kinetics of poly(aryl ether ether ketone ketone) (PEEKK)[END_REF][START_REF] Jenkins | Structure evolution in melt crystallised PEEK[END_REF].

The primary crystallization weight factor (w 1 ) for PEKK 6002 and 7002 is plotted as a function of temperature in Fig. 6 with the fitting of the linear model (Eq. ( 18)). Results for a 1 and b 1 are presented in Table 1. We observed that for PEKK 6002, w 1 decreases with temperature, knowing that w 1 +w 2 = 1 (see eqs. ( 1)-( 3)), it shows that the formation of the secondary crystallization is promoted by the low

Table 1

Parameters of the kinetic models for the crystallization rate constants (K 1 and K 2 ) and the weight factor (w 1 ) for poly(ether ketone ketone) (PEKK) 6002 and 7002 crystallized from the melt. degree of supercooling. This observation is also reported by other authors [START_REF] Quiroga Cortés | Morphology and dynamical mechanical properties of poly ether ketone ketone (PEKK) with meta phenyl links[END_REF][START_REF] Dasriaux | Evolution of the amorphous fraction of PEEK during annealing at atmospheric and high pressure above the glass transition temperature[END_REF], with an increase in the first endothermic peak, which occurs 10 °C above the annealing temperature attributed to the melting of the secondary crystallization. In fact, for high annealing temperatures, the diffusion process is very important, with a high mobility of chains allowing the crystalline phase to be enhanced. For PEKK 7002, w 1 seems to be quasi-constant as a function of temperature with a value of around 0.78. PEKK 7002 wt factors are higher than for PEKK 6002, which could be attributed to a more organized crystallization process that leaves less room for the growth of secondary crystallization.

Time-temperature-transformation (TTT) diagrams

As the relative crystallinity of PEKK 6002 and 7002 depending on time can be determined for any temperatures with the modeling of the Hillier parameters, we built TTT diagrams of the relative volume fraction crystallinity for the overall crystallization for PEKK 6002 (Fig. 7a) and 7002 (Fig. 7b). This kind of representation allows the crystallinity of the material to be determined for any isothermal crystallization process as well as the corresponding crystallinity ratios of each crystallization mechanism. As expected, crystallization kinetics are fastest around 240 °C for PEKK 6002 and 245 °C for PEKK 7002. The TTT diagrams show the very fast crystallization kinetics of PEKK 7002 compared to PEKK 6002. For example, it takes 10 min for PEKK 7002 to fully crystallize at 240 °C, whereas only about 50% of PEKK 6002 crystallizes in the same conditions.

Conclusion

The modeling of PEKK crystallization kinetics proves difficult due to a significant secondary crystallization, which has to be taken into account in the process modeling. Furthermore, to typically assess the crystallization enthalpy based on the DSC thermogram, the integration of the crystallization heat flow implies an extrapolation of the beginning of isothermal crystallization stage. If the extrapolation is done using a straight line, this generally involves a truncation of the beginning of all isothermal crystallization peaks. As a result, this truncation leads to underestimating the crystallization enthalpy value, constraining the Avrami exponent to a value close to 2, and sometimes introducing an induction time to fit the data with overall crystallization kinetic modeling.

To solve this problem, we established a new methodology from the derivative of the Hillier isothermal crystallization kinetic model, which can be directly fitted on the crystallization heat flux measured by DSC. The Hillier model allows us to deconvolute the contribution of the primary and secondary crystallization processes. This approach has been successfully applied to PEKK 6002 and 7002 for isothermal crystallizations from the glassy state and from the melt between 200 °C and 300 °C. PEKK 7002 was shown to crystallize much faster than PEKK 6002, with rapid crystallization kinetics around 245 °C. According to the linear modeling of the primary crystallization weight factor w 1 , the secondary mechanism was found to be enhanced for a low degree of supercooling due to the high diffusion process for PEKK 6002, whereas it was observed to be constant for PEKK 7002. Finally, the isothermal TTT diagrams for crystallinization built for PEKK 6002 and 7002 provide a useful tool for PEKK processing.

Fig. 1 .

 1 Fig. 1. Differential scanning calorimetry thermogram of poly(ether ketone ketone) 6002 crystallized from the melt during isothermal crystallization at 230 °C with a linear extrapolation of the crystallization peak beginning (dashed line) inducing a crystallization induction time t i (dotted line).

Fig. 2 .

 2 Fig. 2. Poly(ether ketone ketone) a) terephthaloyl and b) isophthaloyl isomers.

Fig. 3 .

 3 Fig. 3. Differential scanning calorimetry thermograms of neat poly(ether ketone ketone) 6002 crystallized from the melt at (a) 230 °C and(b) 200 °C with the fitting of the derivative Hillier model (solid line) and primary (dash curve) and secondary crystallization (dash-dot curve) modeling.

Fig. 4 .

 4 Fig. 4. Comparison of the fitting of the derivative Hillier model for neat poly (ether ketone ketone) (PEKK) 6002 and 7002.

  independent of temperature.

Fig. 5 .

 5 Fig. 5. Logarithmic plot of the (a) primary and (b) secondary crystallization rate constants versus temperature for neat poly(ether ketone ketone) (PEKK) 6002 and 7002 isothermal crystallization from the melt using the Hoffman and Lauritzen model fitting.

Fig. 6 .

 6 Fig.6. Weight factor of the primary crystallization versus temperature isothermal crystallization for poly(ether ketone ketone) (PEKK) 6002 and 7002 from the melt using the linear model fitting.

Fig. 7 .

 7 Fig. 7. Time-temperature-transformation diagram of the relative crystallinity of poly(ether ketone ketone) (PEKK) (a) 6002 (b) and7002 crystallized from the melt.

Acknowledgments

The authors wish to thank the Arkema Company for providing the PEKK materials. ANRT (CIFRE N°2014/0797) is also gratefully acknowledged for financial support.


L Meyrand
email: meyrand@lma.cnrs-mrs.fr

E Sarrouy
email: sarrouy@lma.cnrs-mrs.fr

B Cochelin
email: cochelin@lma.cnrs-mrs.fr

Low-dimensional Nonlinear Modes computed with PGD/HBM and Reduced Nonlinear Modal Synthesis for Forced Responses

This work proposes an algorithm allowing to perform a fast and light computation of branches of damped Nonlinear Normal Modes (dNNMs). Based on a previous work about undamped NNMs (uNNMs), it couples Proper Generalized Decomposition (PGD) features, harmonic balance and prediction-correction continuation schemes. After recalling the main contributions of the method applied on an example with cubic nonlinearities, the issue of a reduced nonlinear modal synthesis is briefly addressed.

The differential equations governing the motion of a nonlinear dynamical system can usually take the following form after a spatial discretization, where the nonlinear efforts f nl (x, ẋ) are separated from the linear ones:

Mẍ(t) + C ẋ(t) + Kx(t) + f nl (x(t), x(t)) = f e (t) (dim.: N) (1)
Assuming the Rosenberg's framework [START_REF] Rosenberg | On nonlinear vibrations of systems with many degrees of freedom[END_REF], an undamped NNM is a set of limit cycles of Eq. ( 1) from which dissipative and external forcing terms are put to zero. Shaw and Pierre extended this first definition of NNMs to the case of dissipative systems: a (damped) NNM is a two-dimensional invariant manifold in the phase space [START_REF] Shaw | Non-linear normal modes and invariant manifolds[END_REF]. Hence, one wants to obtain a dNNM by computing a set of pseudo-periodic solutions of Eq. ( 1) with f e (t) = 0.

A frequential reduced algorithm coupling a PGD approach with an harmonic balance method (HBM) is implemented to make a quick and compact dNNM computation. The first step of the PGD process is separating the variables, here space and time:

x(t) ≈ m j=1 p j q j (t) ⇔ x(t) ≈ Pq(t) with P = [p 1 , . . . , p m ] (2)
m ≪ N is a positive integerand denotes the number of PGD modes (p j , q j (t)) used for x decomposition. P is the (N × m)-sized matrix of the m PGD mode shapes p j , and q(t) is the vector containing the time dependence of each PGD mode. Then a spatial subproblem S m and a temporal subproblem T m are defined from specific weak formulations [START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF]. The calculus, its notations and operators are detailed in [START_REF] Meyrand | Nonlinear normal mode continuation through a Proper Generalized Decomposition approach with modal enrichment[END_REF]. Given the spatial matrix P, T m is a set of m ordinary differential equations of order 2. A complex HBM [START_REF] Joannin | A reduced-order modeling technique for the study of non-linear vibrations in dissipative systems[END_REF] is then implemented to obtain another nonlinear algebraic system for T m from the following solution form:

q(t) = a 0 √ 2 + H k=1 e -kβt [a k cos (kωt) + b k sin (kωt)] (3)
β-dependant linear terms can easily be separated from the classic undamped HBM matrices. Given the temporal part q(t), the spatial subproblem S m is a N × m nonlinear algebraic system which takes into account the damping ratio β into the definition of its operators. The PGD/HBM solver is obtained by integrating S m and T m into an alternated directions solver. This PGD/HBM solver is eventually embedded into a continuation scheme as a corrector in order to build the dNNM branch. The choice of predictor is left to the user.

The first point of the branch is described with only one PGD mode which contains the shape and the frequency of the underlying linear damped mode with a null amplitude. When the error criterion can no longer be met through the continuation, the size of the description m is incremented and some new spatial and temporal information is added to the PGD description. Unlike Grolet and Thouverez [START_REF] Grolet | On the use of the proper generalised decomposition for solving nonlinear vibration problems[END_REF] who initialized q(t) with random values, we propose to process T m first based on the shape on the next damped linear modes. The full algorithm adds new modal data on the fly, only when it is necessary with respect to the error criterion.

It should be noted that two variants have been implemented: oPGD (optimized PGD) recomputes the whole P matrix at each point of the branch whereas pPGD (progressive PGD) only computes the last shape p m+1 when it is introduced. Although oPGD generally needs less PGD modes than pPGD, the computational cost is higher.

The method is here applied on a cantilever beam with a cubic spring at its free end investigated in [START_REF] Meyrand | Nonlinear normal mode continuation through a Proper Generalized Decomposition approach with modal enrichment[END_REF]. A 1% modal damping is added on each linear mode. The Frequency-Energy Plot given on Fig. 1 illustrates the hardening effect of the cubic nonlinearity while the conservative mechanical energy grows and the damping-energy dependence of the NNM.

where s is an index for the dNNM branch and φ is the phase of the response. Although this decomposition is slightly different of the one given in [START_REF] Joannin | A reduced-order modeling technique for the study of non-linear vibrations in dissipative systems[END_REF], the two equations used to solve for s and φ are similar.

Figure 1 :

 1 Figure 1: (a) Schematic diagram of beam+cubic spring (b) Main branch of dNNM1 (c) Damping ratio with respect to energy. Squares: Solution points where a PGD mode is added.Eventually, this compact description of NNMs can be embedded into a reduced modal synthesis solver in order to quickly build Frequency Response Functions by looking for x solution of Eq. (1) as follows:x(t) = P(s) a 0 (s) √ 2 + H