D C Vaz

F Trier

A Dyrda
email: adyrdal@amu.edu.pl

A Johansson

K Garcia

A Barthélémy

I Mertig

J Barnaś

A Fert

M Bibes

Determining the Rashba parameter from the bilinear magnetoresistance response in a two-dimensional electron gas

Two-dimensional (2D) Rashba systems have been intensively studied in the last decade due to their unconventional physics, tunability capabilities, and potential for spin-charge interconversion when compared to conventional heavy metals. With the advent of a new generation of spin-based logic and memory devices, the search for Rashba systems with more robust and larger conversion efficiencies is expanding. Conventionally, demanding techniques such as angle-and spin-resolved photoemission spectroscopy are required to determine the Rashba parameter αR that characterizes these systems.

Here, we introduce a simple method that allows a quantitative extraction of αR, through the analysis of the bilinear response of angle-dependent magnetotransport experiments. This method is based on the modulation of the Rashba-split bands under a rotating in-plane magnetic field. We show that our method is able to correctly yield the value of αR for a wide range of Fermi energies in the 2D electron gas at the LaAlO3/SrTiO3 interface. By applying a gate voltage, we observe a maximum αR in the region of the band structure where interband effects maximize the Rashba effect, consistently with theoretical predictions.

I. INTRODUCTION

Magnetotransport phenomena in low-dimensional systems are extremely useful for applications in (nano)electronic devices. After a first generation of spintronic devices based on metallic magnetic multilayers [START_REF] Baibich | Giant Magnetoresistance of (001)Fe/(001)Cr Magnetic Snperlattices[END_REF][START_REF] Binasch | Enhanced magnetoresistance in layered magnetic structures with antiferromagnetic interlayer exchange[END_REF][START_REF] Camley | Theory of giant magnetoresistance effects in magnetic layered structures with antiferromagnetic coupling[END_REF][START_REF] Barnaś | Novel magnetoresistance effect in layered magnetic structures: Theory and experiment[END_REF][START_REF] Julliere | Tunneling between ferromagnetic films[END_REF][START_REF] Braduc | Spintronic Phenomena: Giant Magnetoresistance, Tunnel Magnetoresistance and Spin Transfer Torque in Giant Magnetoresistance (GMR) Sensors. From Basis to State-of-the-Art Applications[END_REF][START_REF] Barnaś | Giant Magnetoresistance and Applications in Magnetism of Surfaces, Interfaces, and Nanoscale Materials[END_REF], the next generation exploits the idea of a full electrical control of the spin degree of freedom via electricallyinduced spin torques of various origins. In this context, a key objective is to exploit spin-orbit coupling to achieve efficient charge-to-spin interconversion, owing to the spin Hall effect [START_REF] Dyakonov | Current-induced spin orientation of electrons in semiconductors[END_REF][START_REF] Hirsch | Spin Hall Effect[END_REF][START_REF] Engel | Theory of spin Hall effect in semiconductors[END_REF][START_REF] Sinova | Spin Hall effects[END_REF] or the current-induced spin density also known as the inverse spin-galvanic effect, or alternatively as the Edelstein effect [START_REF] Edelstein | Spin polarization of conduction electrons induced by electric current in two-dimensional asymmetric electron systems[END_REF][START_REF] Aronov | Nuclear electric resonance and orientation of carrier spins by an electric field[END_REF][START_REF] Ganichev | Conversion of Spin into Directed Electric Current in Quantum Wells[END_REF][START_REF] Ganichev | Spin-galvanic effect[END_REF].

While in the beginning of the study of the Edelstein effect great emphasis was put on semiconducting heterostructures, owing to their strong bulk and structural inversion asymmetries [START_REF] Awschalom | Challenges for semiconductor spintronics[END_REF][START_REF] Winkler | Spin-orbit Coupling Effects in Two-Dimensional Electron and Hole Systems[END_REF][START_REF] Fabian | Semiconductor Spintronics[END_REF], nowadays two dimensional (2D) systems show the most promise, with encouraging results found at the surfaces of topological insulators (TIs), two-dimensional graphene-like heterostructures and oxide interfaces [START_REF] Soumyanarayanan | Emergent phenomena induced by spin-orbit coupling at surfaces and interfaces[END_REF][START_REF] Varignon | A new spin for oxide interfaces[END_REF][START_REF] Song | Spin Proximity Effects in Graphene/Topological Insulator Heterostructures[END_REF][START_REF] Cummings | Graphene Spintronics in 2D Materials. Properties and Devices[END_REF][START_REF] Han | Graphene spintronics[END_REF]. In these systems, the Rashba spin-orbit coupling acts as an internal momentum-dependent magnetic field that ensures a fixed relative orientation between the electron spin and momentum, an effect known as spin-momentum locking.

Recently, a new kind of magnetoresistance, the bilinear magnetoresistance (BMR), was reported in materials with strong spin-orbit coupling [START_REF] Pan He | Bilinear magnetoelectric resistance as a probe of three-dimensional spin texture in topological surface states[END_REF][START_REF] Dyrda | Spin-Momentum-Locking Inhomogeneities as a Source of Bilinear Magnetoresistance in Topological Insulators[END_REF]. Here, two resistive states are observed depending on the relative orientation between the applied (bias) electric field and the magnetic field, as shown in Fig. 1(a). In contrast with other recently reported unidirectional magnetoresistance effects (e.g., spin Hall magnetoresistance in multilayered systems), the BMR appears in a single material. The effect was observed in topological insulators [START_REF] Pan He | Bilinear magnetoelectric resistance as a probe of three-dimensional spin texture in topological surface states[END_REF] and in materials having surface or interface two-dimensional electron gases (2DEG) with Rashba interactions. For the latter, clear experimental evidence has been shown at the surface states of Ge (111), where the BMR in a magnetic field of 1T represented 0.5% of the zero field resistance, and was larger than the magnetoresistance of standard symmetry [START_REF] Guillet | Observation of Large Unidirectional Rashba Magnetoresistance in Ge(111)[END_REF].

The BMR of the metallic surface states of Ge was interpreted by taking into account the locking α R (k × z) • σ between spin σ and momentum k in a Rashba 2DEG (Fig. 1b). As shown in Fig. 1c, the shift ∆k of the Fermi contour in the presence of a current j and the resulting non-equilibrium energy ∆ε ∼ α R (∆k×z)•σ, with ∆k ∼ j, can be described by the introduction of a magnetic field B j ∼ α R (j × z) acting on the spin σ. The field B j adds to the external field B ext , as shown in Fig. 1d. Depending the direction of the current, B j is added to (or subtracted from) the component of the applied field a b c [START_REF] Edelstein | Spin polarization of conduction electrons induced by electric current in two-dimensional asymmetric electron systems[END_REF][START_REF] Aronov | Nuclear electric resonance and orientation of carrier spins by an electric field[END_REF][START_REF] Ganichev | Conversion of Spin into Directed Electric Current in Quantum Wells[END_REF][START_REF] Ganichev | Spin-galvanic effect[END_REF]). The blue (white) arrow represents the orientation of accumulated (depleted) spins. (d) Effective magnetic field B ef f as a result of the current-induced magnetic field Bj (with its sign depending on the polarity of the applied current) and a rotating external field Bext (adapted from [START_REF] Guillet | Observation of Large Unidirectional Rashba Magnetoresistance in Ge(111)[END_REF]). perpendicular to the current, and the BMR appears in cross-terms between B ext and B j . Guillet et al. [START_REF] Guillet | Observation of Large Unidirectional Rashba Magnetoresistance in Ge(111)[END_REF] derived the Rashba parameter from the ratio between the BMR and the quadratic MR in B ext .

V g B ϕ V I d + = B ext B eff B j k x k y ε k y Δk x k x k F -k F j ϕ = 90°l ogic 0 ϕ = 270°I
In this paper, we propose a microscopic theory of the BMR for Rashba 2DEGs. The mechanism is related to the existence of a current-induced magnetic field acting on the electron spins in systems with spin-momentum locking, and it dominates when the Rashba coupling is much larger than the Zeeman field. We apply the model to the interpretation of experimental results in the LaAlO 3 /SrTiO 3 (LAO/STO) system. As will be discussed, we develop the theory for the simple case of a single pair of circular Fermi contours characterized by a Rashba parameter. While in LAO/STO the BMR can be measured at different Fermi levels with different Fermi contours and Rashba splittings, the analysis shows a good agreement with results from tight-binding calculations even in the multiband regime. The fit of the model with the experimental BMR in a given voltage range thus provides us with an effective value of the Rashba parameter in the corresponding Fermi energy range. We demonstrate that the BMR signal may serve as a useful tool to characterize spin-to-charge interconversion efficiencies and to probe the strength of spin-orbit coupling in 2D systems. According to our proposition, an effective Rashba parameter can be determined solely based on the relative amplitudes of the quadratic and bilinear magnetoresistance terms, carrier density, and electron effective mass.

II. THEORETICAL MODEL

The electronic states of a two-dimensional Rashba gas are described by the following Hamiltonian in the planewave basis:

ĤR = 2 k 2 2m * σ 0 + α R (k y σ x -k x σ y ) , (1)
where m * is the effective electron mass, α R is the Rashba parameter, k 2 = k 2 x + k 2 y , and σ j (j = {0, x, y, z}) are the unit and Pauli matrices acting in spin space. The Hamiltonian (1) above has two eigenvalues:

ε ± = 2 k 2 2m * ± α R k. A. Semiclassical picture
For an electric field in the x-direction (E x ) and a relaxation time τ , the momentum k acquires an extra componenet ∆k x = eE x τ / (e < 0) which, from Eq. 1, leads to a non-equilibrium term α R ∆k x σ y = -α R eτ E x σ y / . This term is equivalent to an interaction with a currentinduced Rashba field B j along the ŷ-axis, as shown in Fig. 1d, with gµ B B j y = b j y = -α R eτ E x / (more generally, b j = -α R eτ (ẑ × E)/ ). Here, b = gµ B B is given in energy units, where B is the magnetic field, g the electron g-factor, and µ B the Bohr magneton. It is also convenient to relate the field b j y to the non-equilibrium Edelstein spin polarization S (see Fig. 1c and Refs. [START_REF] Edelstein | Spin polarization of conduction electrons induced by electric current in two-dimensional asymmetric electron systems[END_REF][START_REF] Aronov | Nuclear electric resonance and orientation of carrier spins by an electric field[END_REF][START_REF] Ganichev | Conversion of Spin into Directed Electric Current in Quantum Wells[END_REF][START_REF] Ganichev | Spin-galvanic effect[END_REF]) induced by the current, S y = α R m * 2π ∆k x for E along x. A straightforward calculation leads to the relation b j y = J S y , where J = -2π m * and its sign changes for opposite signs of the effective mass m * . (more details in Ref.

[27]).

B. Microscopic description of BMR in a Rashba 2DEG

To calculate BMR we use the approach proposed recently by Dyrda l et al. [START_REF] Dyrda | Spin-Momentum-Locking Inhomogeneities as a Source of Bilinear Magnetoresistance in Topological Insulators[END_REF]. In the presence of external electric and magnetic fields, the total Hamiltonian can be written in the following form:

Ĥkk = Ĥ0 k + ĤA k δ kk + V imp kk (2)
Ĥ0 k = ĤR + Ĥb eff (3) Ĥb eff = b • σ + b j • σ ≡ b eff • σ (4)
where ĤA k = -ev • A describes the coupling of the charge carriers to an external in-plane dynamical electric field. Here, v is the velocity operator and A denotes electromagnetic vector potential, A(t) = A ω exp iωt/ with A ω = iω E ω . The term b • σ describes the coupling of the electron spin to an external in-plane magnetic field b = (b x , b y ). The term b j •σ expresses the coupling of the spin with the current-induced field described in Section II A. Without losing generality, one can assume that the external electric field is applied in the x-direction. Thus, ĤA k = -ev x A x , and b j • σ = J S y σ y . The relaxation occurs due to scattering from randomly-distributed point-like impurities. The potential V (r) is assumed to be short-range with zero average V (r) = 0, and the second statistical moment V (r)V (r ) = wδ(r -r ).

For impurities distributed randomly at points r i , the potential

V (r) = i v(r -r i ) = i v 0 δ(r -r i )
, where v 0 is the potential of a single impurity and w = n i v 2 0 . Thus, V imp. kk = V kk σ 0 .

C. Relaxation time and conductivity

The relaxation time and conductivity have been calculated within the Green's function formalism. The selfconsistent Born approximation is used to find the relaxation rate Γ (or relaxation time τ ). An important finding is that the relaxation rate/time is dependent on the external magnetic field and non-equilibrium spin density:

Γ(b) = 2τ b = Γ 0 1 + 3 J b y S y 4Γ 2 0 + b 2 8Γ 2 0 , (5)
where

Γ 0 = n i V 2 0 m * 2 2 ≡ 2τ0
is the relaxation rate in the absence of a magnetic field, Γ 0 = Γ(b = 0).

In turn, the longitudinal conductivity can be calculated from the expression [START_REF] Abrikosov | Methods of Quantum Field Theory in Statistical Physics[END_REF][START_REF] Mahan | Many-Particle Physics[END_REF]:

σ xx = e 2 2π Tr d 2 k (2π) 2 vx G R k vx G A k , (6)
where ... denotes the disorder average. In the simplest case, v x G R v x G A can be replaced by v x ḠR V x ḠA , where ḠR,A is the impurity averaged Green's function and V x is the renormalized velocity related to the vertex correction. Since the effective magnetic field is small and treated perturbatively, one can neglect its influence on the impurity vertex correction and consider the Hamiltonian of the 2DEG with only the Rashba term. In such a case, it is known that the impurity corrections to the velocity vertex function lead to cancellation of the socalled anomalous velocity, and V x = kx m * σ 0 (see, e.g., Refs. [START_REF] Inoue | Suppression of the persistent spin Hall current by defect scattering[END_REF][START_REF] Ol | Spin-Hall conductivity in a twodimensional Rashba electron gas[END_REF]).

Taking into account the impurity vertex correction and expanding Green's functions with respect to the effective magnetic field b eff , we arrive at the final expression for the diagonal resistivity

ρ xx = ρ 0 xx + 3π 4 h e 2 α R τ 0 |e| j x b sin φ ε 2 R + ε 2 F + ε F τ 0 h b 2 cos(2φ) ε 2 R + ε 2 F , (7)
where

ρ 0 xx = ρ xx (B = 0) = h e 2 Γ0ε F ε 2 F +ε 2 R
, ε F is the Fermi energy, and ε R = α R k 0 defines the Rashba energy (Rashba field) in the system (here k 0 = √ 2m * ε F / ). This expression clearly shows that even for a simple scalar scattering potential of point-like impurities, a bilinear magnetoresistance proportional do j and B appears, which reveals a sin(φ) angular dependence with an oscillation period of 2π.

D. Magnetoresistance

In the following we briefly discuss the behavior of the magnetoresistance M R = [ρ xx (B) -ρ xx (B = 0)] obtained within the model under consideration. First, we extract the bilinear and B-quadratic (symmetric) components of the magnetoresistance following the definitions:

BM R = [M R(B, j x = +j) -M R(B, j x = -j)] /2 and QM R = [M R(B, j x = +j) + M R(B, j x = -j)] /2.
Thus, the bilinear and quadratic magnetoresistance can be expressed as follows:

BM R = A BM R j x j sin φ, ( 8
) QM R = A QM R cos(2φ), (9)
with the amplitudes

A BM R = 3π 4 h e 2 gµ B |e| α R τ ε 2 R + ε 2 F jB, (10)
A QM R = 3π 4 (gµ B ) 2 e 2 ε F τ ε 2 R + ε 2 F B 2 . ( 11
)
When we use the normalized magnetoresistance, M R = [ρ xx (B) -ρ xx (B = 0)] /ρ xx (B = 0), then Eqs. ( 8) and (9) still hold, but the amplitudes of BMR and QMR take the forms:

A BM R = 3 2 π gµ B |e| α R τ 2 ε F jB, (12)
A QM R = 3 4 gµ B 2 τ 2 B 2 . ( 13
)
From Eqs. ( 8) and ( 9) follows that the QMR oscillates with the periodicity of π (as it is observed in usual anisotropic magnetoresistance experiments), whereas the BMR oscillates with the periodicity of 2π. This behaviour results in a well pronounced asymmetry between φ = π/2 and φ = 3π/2 in the angular dependence of the total magnetoresistance, as well as in an asymmetry of the total magnetoresistance for currents flowing in opposite directions. From Eq. ( 13), the amplitude of the QMR is expected to scale quadratically with an external magnetic field, a trend governed by the relaxation time τ in the system. From Eq. ( 12), the amplitude of the BMR is expected to scale linearly with both applied current and external magnetic field. Here, besides the dependence on external stimuli, A BM R also depends on material-dependent parameters, such as the relaxation time τ , the Rashba parameter α R and the Fermi energy ε F . Lastly, the ratio between both amplitudes:

A BM R A QM R ≡ Λ CS = 2π |e|gµ B α R ε F j B ( 14
)
gives a τ -independent relation from which α R can be found experimentally, provided that ε F and g are known.

Here, this ratio is expressed by universal constants, the ratio of externally-controlled parameters, j and B, and α R . Thus, this ratio, directly proportional to α R , is a quantity characterizing the Rashba coupling. Moreover, even in the case of a more complex multiorbital band structure strongly modified by spin-orbit interactions, the experimentally derived Λ CS can still provide a reasonable estimation of the magnitude of spin-orbit related effects, where an effective Rashba parameter may be considered [α R → α ef f in Eqs. ( 12) and ( 14)]. However, in this regime the model has to be used with caution, since a simple Rashba model cannot properly describe highly anisotropic Fermi contours. In the following section we describe why this approximation is reasonable for the system studied.

III. EXPERIMENTAL RESULTS

We have performed angle-dependent transport experiments in the prototypical Rashba 2DEG found at the interface of LaAlO 3 /SrTiO 3 (LAO/STO) oxide heterostructures. Hall bar devices [see Fig. 2(a)] were fabricated on LAO(1 unit cell)//STO through a combination of electron-beam lithography and a room-temperature deposition of a-LAO (30 nm) using pulsed laser deposition (more details in Ref. [START_REF] Trier | Electric-Field Control of Spin Current Generation and Detection in Ferromagnet-Free SrTiO3-Based Nanodevices[END_REF]). As shown in Fig. 2(b), at low temperatures the device exhibits a positive magnetoresistance in the doped-regime (applying a back-gate voltage of V g = 150 V) and a negative magnetoresistance in the depleted-regime (V g = -30 V), signalling the wellknown transition between weak antilocalization and weak localization regimes [START_REF] Caviglia | Tunable Rashba Spin-Orbit Interaction at Oxide Interfaces[END_REF]. Within this range of gate voltages, the sheet resistance R s increases from 200 Ω/ to 8 kΩ/ [inset of Fig. 2(b)]. Hall effect measurements as a function of V g , shown in Fig. 2(c), confirm the n-type metallic behaviour of the device, similarly to unpatterned samples. By fitting these curves with a model assuming two types of charge carriers (with distinct effective mass) contributing to the conduction, we are able to extract the carrier densities n 1 and n 2 , as shown in Fig. 2(d). For V g < 60 V, linear Hall curves are obtained, indicating that only one band contributes to transport. In this range, the carrier density can be tuned between about 0.5 and 1.5 × 10 13 cm -2 . Above this V g , heavy electron subbands start to be populated, resulting in slightly nonlinear Hall curves and a maximum total carrier density of 3 × 10 13 cm -2 (at V g = 150 V).

We now move to the angle-dependent magnetoresistance experiments. We start by measuring the longitudinal resistance under a rotating in-plane external mag- netic field B, as shown in Fig. 3(a). At φ = 0 • , B lies parallel to the applied DC current I DC . Under B = 9 T and V g = 80 V, the normalized magnetoresistance oscillates with a periodicity of 2π, in agreement with Eq. ( 9). For opposite polarities of I DC , a pronounced asymmetry between φ = 90 • and φ = 270 • is observed, a signature of the unidirectional and bilinear response predicted by our model and previously reported in the Rashba 2DEGs at LAO/STO interfaces [START_REF] Narayanapillai | Current-driven spin orbit field in LaAlO3/SrTiO3 heterostructures[END_REF][START_REF] Choe | Gate-tunable giant nonreciprocal charge transport in noncentrosymmetric oxide interfaces[END_REF], InAs quantum wells [START_REF] Choi | Spin-polarization-induced anisotropic magnetoresistance in a two-dimensional Rashba system[END_REF] and Ge(111) [START_REF] Guillet | Observation of Large Unidirectional Rashba Magnetoresistance in Ge(111)[END_REF]. We performed similar angle-dependent measurements at different magnetic fields and gate voltages values and extracted the BMR and QMR as described earlier. We plot in Fig. 3(b) A QM R as a function of the external magnetic field at V g = 150 V (angle-dependence in the inset). The data can be fitted with Eq. ( 13) (solid black line), confirming the quadratic dependence of A QM R with B. From the fitting, we extract a relaxation time of τ ∼ 1 ps, consistent with previous reports in this system within the same carrier density range [START_REF] Caviglia | Tunable Rashba Spin-Orbit Interaction at Oxide Interfaces[END_REF]. Next, we evaluate A QM R as a function of V g at different external magnetic fields [Fig. 3(c)]. By gating our device, the mobilities increase non-monotonically by one order of magnitude, leading to a similar modulation of the relaxation time, given that τ = (m * µ)/|e| (for carriers with a fixed effective mass m * ). From Eq. ( 13), A QM R is expected to scale with τ 2 , which explains its abrupt increase for larger mobilities at higher gate voltages. We use Eq. ( 14) to extract α R from the ratio A BM R /A QM R , ε F and the g-factor.

R BMR /R s (-) T = 2 K V g = 150 V I DC =
R xx /R s (-) T = 2 K B = 9 T V g = 80 V I DC =
R BMR /R s (-) V g = 150 V -20 V Global back-gate V g LAO//STO 2DEG B V g = 150 V d e f 0.5 T
The following results at different V g correspond to different Fermi energies, different bands, different Fermi contours and different Rashba parameters, as will be discussed further ahead. Thus, the interpretation of the BMR by our model in a given range of gate voltages characterizes an average value of the Rashba parameter for the corresponding Fermi contours. We start by calculating A BM R /A QM R . Since the A QM R fitting was only possible down to V g = 60 V, below this gate voltage α R was derived using Eq. ( 12) with τ so from magnetoresistance fittings to the Maekawa-Fukuyama formula in Ref. [START_REF] Caviglia | Tunable Rashba Spin-Orbit Interaction at Oxide Interfaces[END_REF], which show a good agreement when compared with τ QM R calculated using Eq. ( 13) up to V g ≈ 100 V. ε F can be estimated with ε F = ( 2 πn)/m * for single parabolic bands. We note that, together with g, the estimation of ε F is the biggest source of error when implementing this model, since it requires previous knowledge of m * . Nevertheless, using a fixed m * = 1.3 m (m is the electron mass), g = 0.5 (from Ref. [START_REF] Caviglia | Tunable Rashba Spin-Orbit Interaction at Oxide Interfaces[END_REF]), and the carrier densities showed in Fig. 2(d), we plot in Fig. 4(a) the calculated α R as a function of V g . We observe that α R increases from 2

meV Å up to 50 meV Å at V g = 120 V, where it reaches its maximum value. Beyond this, α R rapidly decreases, reaching a value of 20 meV Å at V g = 150 V. This nonmonotonic behaviour is reminiscent of recent results in SrTiO 3 -based 2DEGs, where the spin-to-charge current conversion efficiency [START_REF] Vaz | Mapping spin-charge conversion to the band structure in a topological oxide two-dimensional electron gas[END_REF] and spin current generation and detection [START_REF] Trier | Electric-Field Control of Spin Current Generation and Detection in Ferromagnet-Free SrTiO3-Based Nanodevices[END_REF] were modulated up and down using gate voltages. The values extracted are in good agreement with the α R reported in other studies for the LAO/STO system in similar carrier density regimes [START_REF] Caviglia | Tunable Rashba Spin-Orbit Interaction at Oxide Interfaces[END_REF][START_REF] King | Quasiparticle dynamics and spin-orbital texture of the SrTiO3 two-dimensional electron gas[END_REF]. Here, we also highlight that, contrary to the other methods, our model is applicable to a broader range of carrier densities, especially where α R takes the highest values.

We further clarify the origin of this modulation by deriving values of α R from an eight-band model Hamiltonian (four Rashba-split bands with opposite spin orientation) composed of two light d xy bands (lower in energy) and one heavy d yz and d xz band (higher in energy), as depicted in Fig. 4(b). Details of the tight-binding modelling can be found in Ref. [START_REF] Vaz | Mapping spin-charge conversion to the band structure in a topological oxide two-dimensional electron gas[END_REF]. We start by calculating α R = ∆k( 2 /2m * ) directly from the energy spectrum, where ∆k = k outer -k inner gives the difference of two neighbouring subbands. In order to respect the individual contribution of each band pair to the QMR and BMR signals, the mean value α R is calculated as an average of all α R from each band pair weighted according to their contribution to the charge conductivity. We observe in Fig. 4(c) that α R increases from 5 meV Å to a maximum of 40 meV Å at around ε F = -40 meV, coinciding with the crossing between the heavy d xz,yz and light d xy subbands [START_REF] King | Quasiparticle dynamics and spin-orbital texture of the SrTiO3 two-dimensional electron gas[END_REF][START_REF] Joshua | A universal critical density underlying the physics of electrons at the LaAlO3/SrTiO3 interface[END_REF], followed by a sharp decrease similar to what was observed experimentally in Fig. 4(a). We can identify a critical energy ε c (or critical V g ) up to which the splitting of the subbands may be directly estimated by a simple Rashba spin-orbit coupling, characterized by circular Fermi contours [left panel in Fig. 4(d)]. For our system, ε c ≈ -60 meV, where heavy d xz,yz subbands start to be populated (also known as Lifshitz transition, represented with a yellow line in Fig. 4). As follows from the TB model results in Fig. 4(c), below this Lifshitz point α R is mostly determined by the electronic states originating from the first pair of subbands. In this range, α R is sufficient to properly describe the Rashba system, and our analytical formulas derived for a simple Rashba model fit well to the experimental data. Above the Lifshitz point the energy spectrum becomes more complex, as displayed in the right panel of Fig. 4(d), and α R can no longer be interpreted as a simple Rashba parameter, but rather as an approximation of the spin-orbit coupling strength in the system, given by α ef f (as described in Section II D). Although in this regime the complex nature of the Fermi contours should in principle prevent us from using our model, we note that the majority of the conductivity above the Lifshitz point is in fact carried by the first pair of subbands, which has the highest density of states. Consequently, in Fig. 4(c) we observe that the magenta line (bands 1&2) has a very large weight on the mean contribution (black curve). We conclude that the calculated α ef f in this range can be approximated to the behaviour of simple Rashba bands, which gives us a fair comparison with respect to the experimental data.

IV. CONCLUSION

In summary, we have developed a model of bilinear magnetoresistance in Rashba systems, based on the existence of a current-induced magnetic field that acts on the electron spins in systems with spin-momentum locking. We employed this model to study the angle-dependent magnetoresistance measured in the Rashba 2DEG at the LAO/STO interface and derived the full gate dependence of α R . A maximum α R is observed at V g ≈ 120 V, coinciding with the crossing between light and heavy subbands, and -importantly -where other magnetotransport probes of Rashba physics are not applicable. Lastly, given the increasing interest of large spin-orbit coupling 2D systems for spin logic and computational applications [START_REF] Manipatruni | Scalable energy-efficient magnetoelectric spin-orbit logic[END_REF], this model provides a useful tool to study new Rashba systems with potentially larger α R , which may exhibit giant spin-to-charge current conversion or spinorbit torque effects.

logic 1 I

 1 FIG. 1. (a) Schematic of the BMR measurement setup under a rotating in-plane magnetic field B. Two different resistive states are expected when B (black arrows) is perpendicular to the applied current I (red arrows). (b) Schematic of the electronic and spin structure of a 2D Rashba system. Arrows represent the electron spin. (c) Schematic of the shift of a Fermi contour by a bias current j, generating a non-equilibrium spin density S (Edelstein effect[START_REF] Edelstein | Spin polarization of conduction electrons induced by electric current in two-dimensional asymmetric electron systems[END_REF][START_REF] Aronov | Nuclear electric resonance and orientation of carrier spins by an electric field[END_REF][START_REF] Ganichev | Conversion of Spin into Directed Electric Current in Quantum Wells[END_REF][START_REF] Ganichev | Spin-galvanic effect[END_REF]). The blue (white) arrow represents the orientation of accumulated (depleted) spins. (d) Effective magnetic field B ef f as a result of the current-induced magnetic field Bj (with its sign depending on the polarity of the applied current) and a rotating external field Bext (adapted from[START_REF] Guillet | Observation of Large Unidirectional Rashba Magnetoresistance in Ge(111)[END_REF]).

FIG. 2 .

 2 FIG. 2. (a) Optical microscope image of a LAO/STO Hall bar device (b) Magnetoresistance (M R) and sheet resistance RS (inset) measurements as a function of the back-gate voltage Vg. (c) Hall measurements as a function of Vg. (d) Carrier densities n1 and n2, extracted from a two-band model fit of the Hall curves (considering different effective electron masses).

FIG. 3 .

 3 FIG. 3. (a) Normalized longitudinal resistance under an in-plane rotating magnetic field B for positive (continuous) and negative (dashed) DC bias current. (b) Amplitude of the first harmonic contribution (angle-dependence in the inset) as a function of B. Here, an AC current was used since QMR is independent of the current direction. This quadratic/symmetric part is given by AQMR = {[RQMR(φ = 0 • ) -RQMR(φ = 90 • )]/Rs} × 100% and RQMR = (RI+ + RI-)/2. (c) AQMR as a function of a back-gate voltage for different B (angle-dependence in the inset). (d) Angle-dependence of the asymmetric contribution for different B, where RBMR = (RI+ -RI-)/2. (e) Amplitude of the bilinear/asymmetric part as a function of B, given by ABMR = [(|RBMR(φ = 90 • )| + |RBMR(φ = 270 • )|)/2Rs] × 100%. Inset: ABMR as a function of the applied current density j. (f) ABMR as a function of a back-gate voltage for different B (angle-dependence in the inset).

FIG. 4 .

 4 FIG. 4. (a)Rashba parameter as a function of the backgate voltage, extracted using the ABMR/AQMR ratio from Eq. (14) (full circles) and the ABMR from Eq. (12) with τso from Ref.[START_REF] Caviglia | Tunable Rashba Spin-Orbit Interaction at Oxide Interfaces[END_REF] (open circles). (b) Top: Band dispersion along the [100] direction, calculated through an eight-band tightbinding model. Bottom: Zoomed area with pronounced enhancement of the Rashba splitting. (c) Energy dependence of the Rashba parameter, calculated for each pair of Rashbasplit bands (colors match the bands in panel b) and for its mean contribution weighted by the conductivity of each band pair (black line). Note that for the orange band pair αR has a negative sign, not shown here for simplicity. (d) Fermi contours below (εF = -80 meV) and above (εF = -50 meV) the Lifshitz transition.

Acknowledgement. -This work received support from the ERC Consolidator Grant 615759 "MINT", ERC Advanced Grant 833976 "FRESCO" and partial support by the National Science Center in Poland as a research project No. UMO-2018/31/D/ST3/02351. A. D. acknowledges support from the French Embassy through the French Government Scholarship (BGF 2019). A. J. acknowledges support from grant CRC/TRR 227 of Deutsche Forschungsgemeinschaft (DFG). We also thank J.-M. George for insightful discussions.


M Gutiérrez

Df Reyes

D Araujo

Jean-Pierre Landesman

E Pargon

Analysis by HR-STEM of the Strain Generation in InP after SiNx Deposition and ICP Etching

published or not. The documents may come   L'archive ouverte pluridisciplinaire

1

Analysis by HR-STEM of the strain generation in InP after SiN x deposition and ICP etching M. Gutiérrez, D. F. Reyes and D. Araujo, University of Cadiz, Department of Materials Science, 11510 Cádiz, Spain.

J.P. Landesman, Univ Rennes, CNRS, IPR (Institut de Physique de Rennes) -UMR 6251, F-35000 Rennes, France E. Pargon, Université Grenoble Alpes, CNRS, LTM (Laboratoire des Technologies de la Microélectronique), F-38000 Grenoble, France

Abstract. Microelectronic processing to fabricate electronic devices on Si has been extensively studied and well characterized forseveral decades. This type of processing has been extrapolated for the fabrication of devices on other semiconducting materials such as Group III-Vs, SiGe and SiC. However, the structural damage and its functional implication are still not fully understood for sensitive materials like InP. The present work investigates the lattice strain induced by the ICP (inductively coupled plasma) dry etching process on InP stripes fabricated after masking the stripes with PECVD (plasma enhanced chemical vapor deposition) grown SiNx. HR-STEM (high resolution scanning transmission microscopy) observations allowed to deduce the strain at different location in the fabricated stripes. Local strains in the range of 10 -3 are observed which should strongly modify the material carrier behaviour.

Introduction. Plasma etching is widely used today in the fabrication process of nanotechnology platforms, both in research, development and production, for the manufacturing of a large number of components and integrated circuits. It is used for most the known semiconducting materials, from Si to III-Vs and even for more "exotic" semiconductors such as diamond or Ga2O3. A good example of a family of such components are photonic devices that can be manufactured using InP. An noted disadvantage of plasma etching is the possible structural and chemical damage in the vicinity of the etched surfaces as well as the introduction of lattice stress and strain [START_REF] Fouchier | Polarized cathodoluminescence for strain measurement[END_REF]. It is known that the mechanical stress in semiconductors alters the electrical and optical properties of the materials as a result of the bandgap modifications [START_REF] Adams | Strained-layer quantum-well lasers[END_REF]. Indeed, small atomic-position variations and stress in the lattice bonds change the semiconductor carrier dynamics. It results in device parameters, such as the emission or absorption wavelength of a laser or detector, spectral output, current threshold, differential gain, polarization state of the output of a laser, waveguide and laser guide properties, being strongly modified.

Therefore, methods to measure the material's stress in semiconductors have been reported in the literature.

Measurements of uniaxial or biaxial strain in material by using photoluminescence (PL) spectroscopy [START_REF] Landesman | Micro-photoluminescence for the visualisation of defects, stress and temperature profiles in high-power III-V's devices[END_REF][START_REF] Bassignana | Photoluminescence and double-crystal x-ray study of InGaAs/InP: Effect of mismatch strain on band gap[END_REF].

Raman spectroscopy [START_REF] Ureña | Raman measurements of uniaxial strain in silicon nanostructures[END_REF][START_REF] Neumann | Raman spectroscopy as probe of nanometre-scale strain variations in graphene[END_REF][START_REF] Gennari | A Raman study of the strain in InP/GaAs heterostructures grown by MOVPE[END_REF], transmission electron microscopy through dislocation characterisation [START_REF] Araújo | Dislocation behavior in InGaAs step-and alternating step-graded structures: Design rules for buffer fabrication[END_REF][START_REF] González | Work-hardening effects in the lattice relaxation of single layer heterostructures[END_REF][START_REF] Perumal | Investigations of the characteristics of strain-free oxidation on InAlAs epilayer lattice matched to indium phosphide[END_REF], X-Ray diffraction [START_REF] Van Der Sluis | Determination of strain in epitaxial semiconductor structures by high-resolution X-ray diffraction[END_REF][START_REF] Lankinen | Crystal defects and strain of epitaxial InP layers laterally overgrown on Si[END_REF], cathodoluminescence (CL) [START_REF] Fouchier | Polarized cathodoluminescence for strain measurement[END_REF][START_REF] Landesman | Mapping of mechanical strain induced by thin and narrow dielectric stripes on InP surfaces[END_REF][START_REF] Rich | Linearly polarized and time-resolved cathodoluminescence study of strain-induced laterally ordered (InP)2/(GaP)2 quantum wires[END_REF], using even polarized CL have been reported. However, these techniques employ indirect methods to determine strain and, in general, have low spatial resolution compared to the requirements of recent device technology processing. Micro-and nano-technology processing of components such as photonic waveguides on InP involve several steps, which can introduce significant crystal strain in, or close to, the active areas. In particular, dielectric stripes at the surface, obtained after PECVD of SiN x or SiO x thin films at moderate temperatures followed by lithography and etching generate mechanical stress which have been extensively described through the well-known edge force model [START_REF] Hu | Film-edge-induced stress in silicon substrates[END_REF]. Recently, detailed measurements have been published for the SiN x /InP system, using PL and CL [START_REF] Fouchier | Polarized cathodoluminescence for strain measurement[END_REF][START_REF] Landesman | Mapping of mechanical strain induced by thin and narrow dielectric stripes on InP surfaces[END_REF]. In these studies, information on the local crystal strain around SiN x stripes is deduced from both the spectral shift of the luminescence lines and from measurement of the degree of polarization of the luminescence emission.

The present contribution brings additional data on this important technological question from high-resolution scanning transmission microscopy (HR-STEM) using procedures, which have been proposed for the analysis of strain at the nanoscale level [START_REF] Hÿtch | Quantitative measurement of displacement and strain fields from HREM micrographs[END_REF][START_REF] Bárcena-González | Strain mapping accuracy improvement using super-resolution techniques[END_REF][START_REF] Chung | Practical and reproducible mapping of strains in Si devices using geometric phase analysis of annular dark-field images from scanning transmission electron microscopy[END_REF][START_REF] Cooper | Strain mapping of semiconductor specimens with nm-scale resolution in a transmission electron microscope[END_REF][START_REF] Grieb | Determination of the chemical composition of GaNAs using STEM HAADF imaging and STEM state analysis[END_REF][START_REF] Zhu | Interface lattice displacement measurement to 1 pm by geometric phase analysis on aberration-corrected HAADF STEM images[END_REF]. Very local evaluation of the strain is carried out in a region of nanometerextension. Such strain may affect surface carrier behaviour as for example the channel of MOSFETs.

The studied sample consists in silicon nitride stripes, deposited on an InP substrate where ICP etching is carried out after SiN x stripe deposition.

Experimental technique. The samples for this study are etched in bulk InP (100) wafers (n-type, S-doped, n=1.8x10 18 cm-3) on which a 500 nm thick SiNx hard mask layer is realized by PECVD. The etched features are rectangular-shaped lines, having lateral dimension L = 0.5, 1, 6, 10 and 20 µm, and aligned along the (110) direction of InP. These lines can be isolated, or grouped in series where the spacing between lines is S = 0.5, 1, 2, 3, 6 or 10 µm. After resist deposition, e-beam lithography, opening of this mask layer forming stripes (using a CF4/CH2F2/Ar plasma) and stripping of resist (using an O2 plasma), the SiNx patterns are transferred into the InP substrate. This is done with a CH4/Cl2/Ar plasma operated in an ICP reactor equipped with a hot cathode. During this step the samples are heated to a temperature of 200 °C. The etch rate under the conditions of this study is approximately 900 nm/min. After etching, the SiNx hard mask and the passivation layers formed on the InP sidewalls are removed by a dual step cleaning procedure using O2 plasma followed by 49% HF dip. Thus, the lattice strain revealed in the present contribution results from both technological steps: SiN mask deposition/removal and the effect of the ICP dry etching. Cross-sectional TEM lamellas were prepared by focus ion beam (FIB) and for the strain analysis a FEI TITAN3 Themis 60-300 microscope, in HR-STEM mode, was used.

The quantitative strain analysis in HR-STEM images is acquired at the [110] zone axis by using geometrical phase analysis (GPA) [START_REF] Hÿtch | Quantitative measurement of displacement and strain fields from HREM micrographs[END_REF][START_REF] Bárcena-González | Strain mapping accuracy improvement using super-resolution techniques[END_REF][START_REF] Chung | Practical and reproducible mapping of strains in Si devices using geometric phase analysis of annular dark-field images from scanning transmission electron microscopy[END_REF][START_REF] Cooper | Strain mapping of semiconductor specimens with nm-scale resolution in a transmission electron microscope[END_REF][START_REF] Grieb | Determination of the chemical composition of GaNAs using STEM HAADF imaging and STEM state analysis[END_REF][START_REF] Zhu | Interface lattice displacement measurement to 1 pm by geometric phase analysis on aberration-corrected HAADF STEM images[END_REF][START_REF] González | Quantitative analysis of the interplay between InAs quantum dots and wetting layer during the GaAs capping process[END_REF][START_REF] Reyes | Impact of N on the atomic-scale Sb distribution in quaternary GaAsSbN-capped InAs quantum dots[END_REF]. The lasts steps of the FIB procedure lamella thinning were carried out at 1kV ion beam accelerated voltage to avoid a possible generation of Ga implantation strain artefacts [START_REF] Balboni | Improved TEM Sample Preparation by Low Energy FIB for Strain Analysis by Convergent Beam Electron Diffraction[END_REF], that can occur above 5kV beam voltage.

Results. During the acquisition of a STEM image there are different scanning distortions. The most important source of error in the measurement of the distortion, comes from the displacement of the beam during scanning. That is, when the beam reaches the end of the i th -line, (𝑥end, 𝑦i) it should then start at the beginning of the next i th +1-line (𝑥 , 𝑦 ). This error in the return of the beam is essentially 1-D, and thus can be eliminated by choosing a scanning direction parallel to the direction in which the deformation is to be studied [START_REF] Bárcena-González | Strain mapping accuracy improvement using super-resolution techniques[END_REF][START_REF] Chung | Practical and reproducible mapping of strains in Si devices using geometric phase analysis of annular dark-field images from scanning transmission electron microscopy[END_REF][START_REF] Cooper | Strain mapping of semiconductor specimens with nm-scale resolution in a transmission electron microscope[END_REF][START_REF] Zhu | Interface lattice displacement measurement to 1 pm by geometric phase analysis on aberration-corrected HAADF STEM images[END_REF], making it possible to obtain highly accurate GPA deformation maps. Therefore, in order to reduce such scanning artefacts, we followed the same experimental approach, i.e. for each studied region, two micrographs perpendicular to each other are recorded. To reduce the effects of drift, sufficient time was allowed to pass between finding a suitable position for the analysis and acquisition of the micrograph. Additionally, the acquisition time was limited to 2s, thus limiting contributions due to the drift. After minimizing the drift in the sample, we acquired pairs of micrographs, perpendicular to each other, in each area of study, aligning the scan with the direction of growth, which we have named parallel (//), and perpendicular () to the growth plane. Figure 2 shows an example of two STEM images perpendicular to each other, acquired in the area taken as reference and the deformation maps extracted from each of them, using GPA. For the choice of the reference vectors 𝑔 𝑔 (inset figure 2a), which will be used for the comparison of the displacements, a rectangular area has been chosen at the bottom of the image. In figure 2c and 2e artefacts due to the scanning are revealed (see arrows), which indicates that those observations must be discarded and keeping only those collinear to the scanning direction, i.e. the 𝜀 and 𝜀 for the // and , respectively. Normally, this type of analysis is carried out near interfaces or areas where there are changes in composition or crystallinity, and the reference area for choosing the vectors 𝑔 and 𝑔 , are within the image one want to analyse. However, because of the significant distance between one region and another in the present study, a different procedure should be here applied. The reference vectors 𝑔 and 𝑔 used in the GPA are carried out from the area labelled here "reference" area (see figure 1), located far from the stripe (>2µm far in the direction of the substrate), where the effects of the lithography processes are below the strain noise. This allows us to estimate the strain variations in the different study areas (numbers 1 to 9) that are far from each other. The calculation of the strain at each point in the micrographs is given by:

𝜀 , , = , , , (1)
where i refers to the  or // direction with respect to the plane surface, n to the region where the strain is estimated (number 1 to 9), 𝑎 , to the  or // lattice parameter in the region n, 𝑎 , to the one of the reference region and 𝜀 , , to the estimated strain  or // to the surface plane and in the region n. show the FFT pattern of the image and the localization of vector 𝑔 and 𝑔 for strain analysis.

In Fig. 3, an example of the strain mapping, recorded on the region 6 is shown. From the atomic displacement recorded in two HAADF micrographs of the same region, acquired at 0 and 90 degrees to the direction of growth, in order to scan the areas parallel to the strain that is to be measured, xx and yy are recorded in Figs.3b 3c respectively. Fig. 3b has more regions where the average strain is positive and is measured to be 222 x 10 - 5 while Fig. 3c has an average strain that is negative and was found to be -164 x 10 -5 . Focusing the attention to the strain observed in the reference zone, one can notice the presence of a non-zero strain (Table I).. To measure these fluctuations, which will indicate the accuracy of the measurement, a small region of image has been taken where no noticeable fluctuations were observed. This sub-region was used to determine the average strain in the reference region, and the value shown in the Table I (reference) is the average strain of all image. At this point it should be noted that 𝑎 , in equation one refers to 𝑎 for the sub region in the reference area.

Therefore, the value reported in table I for the reference micrograph give the method accuracy (± 5 × 10 ).

The latter depends strongly on the stability of the sample below the e-beam. If no drift, no charging effects and resistance to the high energetic incident electron bombardment of the e-beam occur, then larger integration times could be used and lower noise levels can be obtained. Here, the scanning acquisition used an estimated convergence angle of 16 mrad using a camera length of 46 mm. However, these are sufficient to detect important strain variations induced by the SiNx deposition and the ICP process. Fig. 4 shows the general tendency of the strain at different locations. The black arrows indicate if the strain is compressive () or tensile () in both x and y directions. In the lateral locations as in regions 4, 6, 7, 9, yy is compressive while xx is tensile. The results are consistent between the four locations. At the edges as in locations 1 and 3, yy is tensile while xx is compressive. This is also consistent in both locations. In the middle of the stripe (location 5), between locations 4 and 6, the vertical trend is similar to what is observed on the sides but with lower values. Lowering the location towards the "bulk" material as location 8, the strain appears lower below the accuracy of the method, i.e. as low as the reference one. Note that at the location 1, xx is the inverse to that of both location 2 and 3 (see arrow in Fig. 4). However, the values are very low (see table I) and in the expected experimental noise [-60 x 10 -5 , +60 x 10 -5 , i.e. 10 -4 range) even though the accuracy of the technique is much better (10 -5 range). Indeed, the location where the micrograph is recorded is not rigorously symmetric between position 1&3 as the latter is closer to the centre (position 2) of the stripe than the position 1. Such little variation can also affect the deduced strains values. It is impressive that strain value as high as several 10 -3 is observed. It is worth noting that the critical thickness determined by Matthews and Blakeslee (mechanical equilibrium) where plastic strain begins is below 100nm for such strain values. Obviously, this surface elastic strain is a surface effect and its value diminishes deeper into the material probably as a result of the surface bend. However, the carrier dynamic could be strongly affected by such values of strain. The possible origins for the local crystal deformation evidenced in this work are diverse. A first possibility is the fact that the SiNx stripe which is used as the hard mask layer for ICP etching of the InP is itself highly strained, due to the difference in the thermal expansion coefficient of SiNx compared to that of InP (the temperature at which the PECVD process is carried for this layer is typically 250 °C) and to the intrinsic stress generated by ion bombardment during PECVD. However, as mentioned, this SiNx hard mask stripe has been removed prior to the HR-STEM measurements. In addition, separate micro-PL and CL experiments with InP stripes obtained with SiNx hard mask layer grown under different built-in stress conditions (both compressive and tensile) indicate that the state of crystal deformation observed in the final InP stripe -after SiNx removal -does not change with the nature of the built-in stress in the dielectric layer used as the hard mask for etching. In other words, there does not seem to be a "history effect" related to built-in stress in the hard mask dielectric layer.

A second possibility that could explain the residual crystal strain observed here is related to some penetration of atoms / ions / molecules active during the ICP plasma etching process into the InP crystal. Note that the sample temperature during the ICP etching is too low for classical atomic diffusion inside the InP material. However, as previously reported indicate, for the etching of GaAs or InP, channelling or high energetic ion scattering penetration can take place [START_REF] Germann | Energy dependence and depth distribution of dry etching-induced damage in III/V semiconductor heterostructures[END_REF][START_REF] Rahman | Channeling and diffusion in dry-etch damage[END_REF]. In the present case, the dominant species in the ICP plasma, probably atomic Cl or Cl -, penetrate the InP material. This second possibility is more likely to explain our results. As a general observation, the deeper the ICP etching is, the larger is the strain difference detected either by HRSTEM or by CL and micro-PL experiments. These observations tend to confirm that ICP-related ion implantations is responsible of the induced strain observed. The details of the mechanism leading to the spatial variation of the strain (as evidenced in the present study) still remains to be proposed and established quantitatively. This will be the topic for future work, which will include finite element simulations of local crystal deformation under the presence of inhomogeneous concentration of such foreign species in InP after etching features like stripes.

Conclusion.

Subsurface local elastic strain has been evaluated by HR-STEM on ICP fabricated stripes using geometrical phase analysis (GPA) method. Consistent results are obtained all along the analysed surface. A compressive lattice strain is revealed parallel to the surface with an extensive strain perpendicular to the surface, the strain perpendicular to the surface was much greater than the strain parallel to the surface (factor 3-6). In the present case, the temperature of the sample during ICP etching is too low for classical diffusion in the InP. Therefore, the observed behaviour is tentatively explained by the implantation of ICP-related ions of the dominant species in the ICP used plasma (Cl or Cl-). Strain values in the range of some percent are deduced which highlights the importance of evaluating such effects on each lithographic approach when designing device architectures.

Figure 1 Fig. 1 :

 11 Fig.1: STEM-HAADF and EDS micrographs of the studied stripe. The locations where the strain is examined are indicated in the HAADF micrograph. The EDS mapping shows the In, P and Pt chemical distributions. The latter has been deposited during the FIB-dual beam sample preparation to protect the stripe.

Fig. 2 :

 2 Fig.2: HRSTEM images with scanning // and  to growth plane (a) and (d), respectively. The maps (b) and (f) corresponds strain maps 𝜀 and 𝜀 , while the maps (c) and (e) are those that should be discarded because the scanning direction is parallel to that of the strain determination, scanning artefacts are clearly observed (see black arrows). The inset in (c) is a scheme of cross-sectional view of FIB lamella axis election. The inset in (a)

Fig. 3 :

 3 Fig.3: (a) HR-STEM micrograph in location 6, where strain evaluation has been carried out, (b) and (c) Strain distribution obtained by the atomic displacement recorded in HR-STEM micrographs scanned parallel to the strain evaluation direction.

Fig. 4 :

 4 Fig.4: Strain tendency inside the ICP fabricated InP stripe. Surface strain is evidenced (see table I data).

TABLE I :

 I Strain values obtained in the different locations indicated in Fig.4

	Region	xx/10 -5	yy/10 -5
	1	60	32
	2	-60	361
	3	-22	258
	4	171	-71
	5	48	-66
	6	222	-164
	7	179	-83
	8	8	3
	9	157	-60
	reference	-5	-2

Acknowledgments

We thank the Ministerio de Economía y Competitividad (MINECO) of the Spanish Government for funding under Grants No. TEC2017-86347-C2-1-R (DiamMOS project) and No. ESP2017-91820-EXP (Diam-Air, EXPLORA project) from FEDER funds. We thank to EU Frame-work Programme for research ¬ Innovation -H2020 for funding under Grant No. SEP-210184415 (GreenDiamond project) of H2020-LCE-2014-1.640947 program.


Yonghui Xu
email: yoghui.xu@sorbonne-universite.fr

Ivan Delbende
email: ivan.delbende@sorbonne-universite.fr

Daniel Fuster
email: daniel.fuster@sorbonne-universite.fr

Maurice Rossi
email: maurice.rossi@sorbonne-universite.fr

Instability of a swirling bubble ring

Keywords: Instability, Vortex ring, Two-phase flow, Direct numerical simulation

A toroidal bubble or a cylindrical gas jet are known to be subjected to the Rayleigh-Plateau instability. Air bubble rings produced by beluga whales and dolphins however are observed that remain stable for long times. In the present work, we analyse the generation of such toroidal bubbles via numerical simulations, in particular how the process depends on surface tension. Their stability properties are then briefly analysed. For the estimated Reynolds and Weber numbers relative to the bubbles produced by these animals, the presence of a vortex inside and around the bubble is found to strongly stabilize the Rayleigh-Plateau instability.

Introduction

This paper is written in memoriam Y. Couder. Yves was a colleague of deep thoughts and a humanist. In everyday life, he was prone to enthusiasm and passion for knowledge. He was an extraordinary observer and many of his breakthroughs originated from careful observations. As a modest homage, this paper presents a topic that takes its origin from an observation on Youtube! In his turbulence experiments [START_REF] Couder | Observation expérimentale de la turbulence bidimensionnelle dans un film liquide mince[END_REF][START_REF] Cadot | Characterization of the low pressure filaments in three-dimensional turbulent shear flow[END_REF], Yves used bubbles to visualize vortices. The present work concerns a topic somewhat related: the generation and dynamics of an air bubble in the form of toroidal rings. It was some kind of surprise to see dolphins [3] or beluga whales [4] generate bubble rings from their mouth or blowholes and thereafter play for a long period of time with rings in which any sign of instability was absent. Indeed an axisymmetric columnar jet of fluid 1 embedded in a fluid 2, when slightly disturbed, is known to break into small droplets or bubbles. This subject was first studied by Lord Rayleigh [START_REF] Rayleigh | On the instability of a cylinder of viscous liquid under capillary force[END_REF][START_REF] Rayleigh | Xix. On the instability of cylindrical fluid surfaces[END_REF]. He computed how infinitesimal perturbations (more specifically normal modes of axial k and azimuthal m wavenumbers) evolve in time and obtained analytically the complex eigenvalue s(k, m) = s r + is i , and precisely s r their growth rate as a function of wavenumbers k and m. For an inviscid two-phase column of radius a 0 , axisymmetric perturbations m = 0 are the most dangerous ones and the dimensionless eigenvalue s is a function of the dimensionless wavenumber k: s2 1 + ρ (2) ρ (1) K 0 ( k)I 0 ( k)

K 1 ( k)I 0 ( k) = k 1 -k2 I 0 ( k) I 0 ( k)
with s ≡ a 3 0 ρ (1) σ s, k ≡ ka 0 .

(

) 1
For 0 k 1, s is real with the two possible signs, meaning that the jet is unstable. For k 1, s = is i is purely imaginary indicating a stable jet. This result can be explained by the theory of minimal energy: a multiphase system tends to possess a minimum capillary energy at rest. This instability known as Rayleigh-Plateau instability is due to a capillary mechanism. It breaks the jet into several drops or bubbles when k 1 so that the interface area is minimized. In order to stabilize Rayleigh-Plateau instability, an added rotation is an interesting possibility for applications, e.g. liquid atomization, spray generation, combustion processes. Previous works on the effect of rotation are those of [START_REF] Hocking | The stability of a column of rotating liquid[END_REF][START_REF] Kubitschek | The effect of viscosity on the stability of a uniformly rotating liquid column in zero gravity[END_REF][START_REF] Kubitschek | Helical instability of a rotating viscous liquid jet[END_REF][START_REF] Gillis | The stability of a rotating viscous jet[END_REF][START_REF] Weidman | On the instability of inviscid, rigidly rotating immiscible fluids in zero gravity[END_REF]. In our case, we introduce this aspect to understand the video of dolphins playing with ring bubbles!

The present paper covers the three stages observed during the dolphin's video inside three separate sections. First, Section 2 introduces the generation of bubble rings through blowhole or mouth. To do so, one considers an idealized axisymmetric nozzle containing initially a gas bubble which is rapidly expelled. This part uses direct numerical simulation. A second stage in which the toroidal bubble is stretched is analyzed in Section 3 by analytical means. As an end result, this provides the vorticity profile reached after the stretching period. The third stage is discussed in Section 4 in which we study the stability of this new profile. The ring is replaced by a straight columnar bubble embedded in a vortex and a standard linear stability method is employed to investigate how superimposed infinitesimal perturbations evolve in time.

Two-phase ring generation

A bubble ring may appear in various ways: a first method consists in generating tiny gas bubbles in advance in a region which is thereafter crossed by a vortex ring produced away from this location [12]. The tiny bubbles migrate because of centrifugal force inside the vorticity ring core and then coalesce to form a toroidal gas bubble transported by the pre-formed vortex ring. In a second method, the vortex is generated at the same spot where tiny bubbles are located [13]. An other method consists in initializing a unique gas bubble inside a nozzle. Ejecting the fluid then causes the bubble to deform and break, yielding a toroidal bubble at the same time and location in which vorticity rolls up, e.g. near the trailing edge of the nozzle. This latter case is probably pertinent for vortex ring bubbles produced through the mouth or blowholes by belugas or dolphins. The presence of the tongue and its motion probably facilitate the toroidal bubble formation as mentioned and demonstrated by a scuba diver [14] but it is shown here not to be an essential feature.

We study below the toroidal bubble production through the mouth or blowholes in an idealized and simplest fashion without introducing the effect of the tongue. The geometry is an axisymmetric nozzle of radius R 0 and of thickness 2b 0 = R 0 /4 with an edge at 26.5 • (see Figure 1a). The nozzle edge enables to fix the region of boundary layer separation and ends at x = 2.8R 0 . Initially the two fluids are at rest and inside the tube, an air bubble of axial length L B is located C. R. Mécanique, 0000, 1, n 0, 000-000 within the axial interval x ∈ [2.6R 0 -L B , 2.6R 0 ] (see Figure 1c). The numerical domain (see Figure 1b) has been checked to be large enough to avoid boundary effects during the time of simulation: it is a square box (x, r ) ∈ [0, 16R 0 ] × [0, 16R 0 ]. Outflow boundary condition (stress free) are imposed at the downstream domain boundary x = 16R 0 , impenetrability and slip condition on the lateral wall r = 16R 0 and no-slip everywhere else except for the region x = 0 and r R 0 where the velocity is imposed. The air bubble is suddenly expelled from the nozzle and injected inside the ambient water. This occurs together with the phenomenon of roll-up of the boundary layer vorticity. The velocity field far inside the nozzle (at x = 0) mimics this rapid ejection

u x (x = 0, r, t ) = U 0 erf t τ erf(η), u r = 0 with η = R 0 -r ν (w) t ≥ 0, (2)
where erf(x) ≡ 2 π x 0 exp(-s 2 ) ds stands for the error function. The center velocity is changing from rest to a constant value U 0 within a time scale τ. For this problem we assume that τ R 0 /U 0 : the center velocity reaches its constant value rapidly compared to the time of roll-up itself. In our simulations, this parameter is fixed at τ ≡ τU 0 /R 0 = 0.1. Moreover the flow is almost spatially uniform inside the tube except near the tube border r = R 0 where a boundary layer is present to ensure no-slip at the wall. This solution is close to what would be observed in a long tube subjected to a sudden pressure gradient. Kinematical viscosity of air and water are taken to be respectively (dimensional values are given in SI units) ν (a) ∼ 1.8 × 10 -5 and ν (w) ∼ 9 × 10 -7 and air and water density respectively ρ (a) = 1, ρ (w) = 10 3 . Finally, air-water surface tension is assumed to be σ = 73.4 × 10 -3 . Based on the nozzle radius R 0 as characteristic length scale and U 0 as the velocity scale, as well as water density and viscosity, this flow is defined by several dimensionless numbers: a Reynolds number, a Weber number and the relative length r B of the bubble:

Re = U 0 R 0 ν (w) , W e = ρ (w) U 2 0 R 0 σ , r B ≡ L B R 0 . ( 3
)
Note that Re/W e 1/2 = (ρ (w) σR 0 ) 1/2 /µ (w) is independent on velocity and is equal to Re/W e 1/2 ≈ 10 4 R 0 1/2 when R 0 is expressed in SI units. The boundary layer thickness δ 0 verifies δ 0 /R 0 = t /Re, where a nondimensional time variable t ≡ U 0 t /R 0 is introduced. This thickness remains small for advective times much less than the Reynolds number t

Re. The largest vorticity occurs at the tube wall. In the bubble, for times such that t is of order one, the dimensionless vorticity maximum is of order

SI units R 0 = 0.02,U 0 = 0.1 R 0 = 0.1,U 0 = 0.1 R 0 = 0.02,U 0 = 0.3 R 0 = 0.
R 0 U 0 max(ω) ∼ - 2 π Re ν (a) /ν (w) ∼ - Re 15.7 . ( 4
)
For Re = 5000, the maximum is close to the value 17.8. Using dimensionless quantities, the entrance condition reads

u x = erf t τ erf(η), u r = 0, with η = 1 -r t /Re . ( 5
)
The axisymmetric flow is described by the velocity fields u i (x, r, t ) with i = 1, 2 related to the pseudo-scalar vorticity component ω θ (x, r, t ), the plane (x, r ) being oriented by the normal unit vector e z . In the present work, we focus on the transient period once the ejection has started: we hence study how the roll-up is modified by surface tension. For the formation process of the ring, axisymmetry is an adequate hypothesis but later on, one expects three-dimensional instabilities to appear on the starting bubble ring for a monophasic [START_REF] Saffman | Vortex Dynamics[END_REF][START_REF] Lim | Vortex rings[END_REF] as for a two-phase vortex [START_REF] Saffman | Vortex Dynamics[END_REF]. The numerical code we use in the present work is the Basilisk code [START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations[END_REF]. It has been largely tested and used to simulate two-phase flow problems. One problematic aspect is the motion of the triple line between bubble, solid and liquid, present during the initial ejection stage. We will not go into numerical details here but let us simply mention en passant that Basilisk numerically mimics the presence near the triple line of a Navier condition, i.e. the presence of a slip length condition which is used in the phenomenology of triple line motion. Here the slip length is half the mesh size near the solid boundary, namely 1 2 ∆ c = 0.0039 R 0 . The present work is not intended to be an extensive investigation of the parameter space of this problem, which in itself would cover several articles. We simply provide here some evidence about the generation of a bubble ring. Table 1 gives estimated Reynolds and Weber numbers for the toroidal bubble produced by scuba divers, beluga whales or dolphins. As a consequence, we use a realistic Reynolds number which is still affordable by numerical simulations running on a standard workstation. For most simulations, the Reynolds number is set at Re = 5000. In the monophasic case, the standard vortex is formed by vorticity sheet roll-up. In the two-phase situation, this mechanism is still active but the bubble expelled from the tube may break up in the same time interval. More importantly, the bubble contains and is surrounded by positive vorticity transported from the boundary layer at the nozzle wall but also from positive and negative vorticity produced at the interface. This interfacial vorticity (related in particular to capillary waves and bubble retraction) is shed in the bulk in the form of vortices. At later times, such vortices may induce a dynamics on the bubble as do the vortices generated by the nozzle wall vorticity roll-up.

The first set of simulations corresponds to a small Weber number W e = 3. For a length r B = 0.6 (Figure 2), the toroidal vortex is not formed but a bubble remains with vortices inside. This means that when surface tension is too strong, inertial effects are too weak to deform the bubble and thus to generate a toroidal bubble. The second simulation (Figure 3) corresponds to a larger width r B = 2.4. In this instance, the bubble ring is again not formed but the starting period is quite complicated with phase of disconnection followed by reconnection. In this process, multiple C. R. Mécanique, 0000, 1, n 0, 000-000 vortices shed at the interface tend to favor this reconnection process. For larger bubbles, the starting vortex is shed inside the air bubble. As a conclusion, it can be checked that for small Weber numbers (e.g. W e = 3), the bubble that is formed is not toroidal!

The second set displayed on Figures 4-5 corresponds to a large Weber number W e = 100 for r B = 0.6 and r B = 2.4 respectively. In the first case (r B = 0.6, Figure 4), the bubble becomes a thin layer near the centre and a rim forms near r = R 0 with the vorticity roll-up ( t = 2.8). The layer breaks thereafter into small fragments since surface tension is too weak and a toroidal bubble is formed near the edge attached to the nozzle edge ( t = 4.8). The flow then detaches the toroidal bubble ( t = 7.2) and the system evolves through capillary waves dynamics and interactions with vortices ( t = 15.2). For a larger r B (here r B = 2.4, Figure 5), the boundary layer is not thin but the roll-up of the vortex sheet is still located near r = R 0 ( t = 4.8). At the center, the two interfaces get nearby ( t = 7.2) and then reconnect forming a unique toroidal bubble ( t = 8). Later the bubble is stretched ( t = 11.2) before being detached ( t = 13.6). The result is a larger toroidal bubble ( t = 19.6). The process is identical for r B = 8 (Figure 6). For even larger values of r B , the bubble does not break and the vorticity roll-up occurs inside the gas bubble. The largest toroidal bubble appears for a value of r B which is clearly dependent on the Weber number.

The intermediate case W e = 20 is an in-between situation (Figures 789) in which a fine balance between the classical roll-up of vorticity and the bubble break-up occurs. For r B = 2.4 (Figure 7), the same process occurs as at larger W e: a toroidal bubble is formed, first attached to the edge ( t = 9.6) then detached ( t = 14.4). Thereafter its cross-section becomes almost circular ( t = 20, t = 24). For r B = 1.2 (Figure 8), the bubble is detached ( t = 5.6) before being broken ( t = 9.6). This is due to the presence of a vortex that was previously shed by the boundary layer. More generally, vorticity production at the interface makes the process quite complex since the toroidal bubble interacts with the vorticity generated at its interface beforehand, or generated by the wall boundary layer (see Figure 8 at times t = 9.6, 20 and 30). In many simulations for W e = 20, the chaotic interaction between such vortices and the toroidal bubble brings this bubble towards the symmetry axis, merging it into a simple bubble. This process is seen here on a simulation for r B = 0.6 (Figure 9): the toroidal bubble is formed at t = 4.8 but its reconnects towards the center at t = 8. This was also observed for r B = 2.4 by slightly changing the Navier condition for the triple line. It is worth mentioning that modifying the Reynolds number, as illustrated in Figure 10, and/or the position of the bubble may change some dynamical features but does not modify the overall picture described above. When the toroidal bubble is formed in the ejection problem, the structure of the vorticity field is complex inside the bubble and around it (see Figure 11). This is also seen on vorticity profiles (see Figure 12). One may define a bubble-ring total length L B which is of order 2πR 0 , as well as the size a B of the toroidal bubble section S and the mean circulation Γ B inside the bubble: Based on such quantities, one computes the mean vorticity inside the bubble Γ B /(πa 2 B ). For a uniform vorticity inside a bubble of radius a B , the azimuthal velocity at the interface is equal to Γ B /(2πa B ). A dimensionless parameter called the swirl number q S can thus be defined as the ratio of this azimuthal velocity and the characteristic capillary velocity

a B = I 1 π R 0 and Γ B = I 2 U 0 R 0 , with I 1 = S dx dr, I 2 = S ω dx dr. ( 6
)
U cap ≡ [σ/(ρ (w) a B )] 1/2 : q S ≡ Γ B 2πa B U cap . (7)
After a transient period, this swirl seems to fluctuate around an average value, as illustrated in Figure 13.

The action of stretching on the bubble ring

In the videos [3,4], dolphins or beluga whales play with the vortex rings and are able by a motion of their head to increase the bubble length. This is typically known as stretching mechanism of vortex tubes. Starting from a vorticity structure similar to the field obtained numerically at the end of the simulations in the previous section (that is those corresponding to vorticity fields in Figure 11), one should model the effect of the dolphin head. In order to remain simple, we analyse the effect of this stretching mechanism by a simplified model which enables us to use analytical methods. An air column of radius a B replaces the two-phase bubble ring with a vorticity field vanishing away from the bubble and the effect of stretching for the ring is introduced in the associated columnar bubble via an axial axisymmetric unsteady stretching

u = γ(t )x, u y (y, z, t ) -1 2 γ(t )y, u z (y, z, t ) -1 2 γ(t )z . ( 8
)
The vorticity structure at t = 0 is characterized by a unique vorticity component ω x (y, z) along x which is identical to ω θ (r, z) obtained at the end of previous simulations (see for instance Figure 11). The stretching occurs during a finite period of time between time t = 0 and t = T S . The stretched solution is again characterized by the unique vorticity component ω x (y, z, t ). This field satisfies the governing equation

Dω x Dt = γ(t )ω x + ν∆ 2D ω x , (9)
where ∆ 2D stands for the two-dimensional Laplacian and D • /Dt for the material derivative:

∆ 2D ≡ ∂ 2 ∂y 2 + ∂ 2 ∂z 2 , D Dt ≡ ∂ ∂t - 1 2 γ(t ) y ∂ ∂y + z ∂ ∂z + u y ∂ ∂y + u z ∂ ∂z . ( 10
)
Across any point at the interface, the velocity field is continuous i.e. u (1) y = u (2) y and u (1) z = u (2) z , and the tangential stress is continuous as well. Let us now introduce the following change of variables [START_REF] Gibbon | Dynamically stretched vortices as solution of the 3D Navier-Stokes equations[END_REF] τ ≡ t 0 S(t ) dt , with S(t ) ≡ exp t 0 γ(t ) dt , χ = S(t )y, η = S(t )z, [START_REF] Weidman | On the instability of inviscid, rigidly rotating immiscible fluids in zero gravity[END_REF] as well as a rescaling of vorticity and velocity fields

ω x (χ, η, τ) = ω x /S(t ), u y (χ, η, τ) = u y / S(t ), u z (χ, η, τ) = u z / S(t ). ( 12
)
After such changes, the field [ u y (χ, η, τ), u z (χ, η, τ)] has the dynamics of an unstretched solution with the same initial condition:

D ω x Dτ = ν ∆2D ω x , with ∆2D ≡ ∂ 2 ∂χ 2 + ∂ 2 ∂η 2 and D Dτ ≡ ∂ ∂τ + u y ∂ ∂χ + u z ∂ ∂η , (13)
associated to equivalent boundary conditions. This trick allows one to obtain the effect of stretching on a vortex profile in two steps: first, one solves the pure two-dimensional advection-diffusion equation ( 13) starting with the velocity field at τ = t = 0 and ending the simulation at τ = τ S ≡ t =T S 0 S(t ) dt ; second, one applies the scalings ( 11)-( 12) back to the initial variables pertinent to describe the physical state. In the videos [3], the bubble ring length L B increases in a substantial manner: assume this increase to be by a factor of F and γ to be constant, this means that F = exp(γT S ) and thus τ S ≡ t =T S 0 S(t ) dt yields (F -1)/γ. Hence, the ratio τ S /T S = (F -1)/ ln F is large: for a factor F = 10, its value is around 4. The first stage of the simulation is thus a pure two-dimensional advection diffusion equation applied during a time τ S larger than T S . In the monophasic case, differential rotation in such an equation is known to accelerate diffusion [START_REF] Bernoff | Rapid relaxation of an axisymmetric vortex[END_REF], so that the vorticity profile rapidly tends towards a smoother axisymmetric one. More generally, two-dimensional vortical flows localized in a finite region are known to decay towards an axisymmetric monopole or, in case the total circulation is zero, towards a dipole [START_REF] Mcwilliams | The emergence of isolated coherent vortices in turbulent flow[END_REF][START_REF] Robert | Statistical equilibrium states for two-dimensional flows[END_REF]. The present case is slightly different because of the presence of two phases. Since ν (a) /ν (w) = 20 is large, vorticity actually diffuses much more rapidly inside the bubble than outside. This implies that vorticity tends to become uniform inside the bubble. Vorticity outside becomes smoother than initially and the interface becomes axisymmetric by damping of capillary waves. This evolution is partially seen on Figure 14. In order to go back to the vortex solution at t = T S (second step), one applies the transformation inverse of ( 11)- (12). The stretching transformation leaves the vorticity uniform in the bubble, the section of which is now circular with radius a new a(t = 0). Indeed, by conservation of mass, a new /a(t = 0) is of the order of 1/ F . If we assume that the circulation Γ B inside the bubble is not changing much, this stretched vortex also possesses a new swirl number

q new S ∼ a(t = 0) a new q S (t = 0) ∼ F q S (t = 0). ( 14
)
The swirl number is thus expected to increase when the bubble ring length increases.

Instability of a columnar interface in the presence of surface tension and rotation

After the stretching phase, a toroidal gas bubble inside a vortex with an almost circular section is generated. In order to understand its stability, we proceed again by studying the same simplified problem: an air column of radius a 0 = a new replaces a two-phase bubble ring with a vorticity field vanishing away from the bubble. This simplified problem is much more tractable since the base flow is homogeneous in the axial direction x. Whenever necessary, the notation Q (p) is explicitly used to represent the quantity Q in the inner phase p = 1 or the outer phase p = 2. The jump of a field Q across a point x of the interface is denoted by The theoretical relation [START_REF] Couder | Observation expérimentale de la turbulence bidimensionnelle dans un film liquide mince[END_REF] indicates that the presence of an external liquid does not change the stability criterion for various density ratios ρ (2) /ρ (1) (Figure 15a). Here, the density ratio is ρ (w) /ρ (a) = 10 3 , which simply makes the system of a columnar fluid 1 more stable: while the Laplace pressure pushes fluid 1 outward, fluid 2 will oppose a force, which slows down the instability and shift the most unstable wavenumber k to a smaller value [START_REF] Eggers | Physics of liquid jets[END_REF]. The significance of dynamical viscosity µ (1) (fixing the ratio µ (2) /µ (1) to the value µ (w) /µ (a) ≈ 50) compared to surface tension σ is described by the dimensionless Ohnesorge number Oh = µ (1) /(ρ (1) σa 0 ) 1/2 , a larger Ohnesorge number indicating an increased influence of viscosity. Since µ (1) /(ρ (1) σ) 1/2 ≈ 6.64 × 10 -5 , Oh ≈ 6.64 × 10 -4 for a radius a 0 = 1 cm. Figure 15b presents the classical dispersion relation for an air column surrounded by water when viscosity is taken into account. Again it is found that viscosity makes the system less unstable and shifts the most amplified wavelength to a longer on from Oh = 10 -2 on. The added inertia and viscosity are hence incapable to bring the system from an unstable to a stable configuration.

[[Q]] ≡ Q (1) (x, t ) -Q (2) (x, t ).
We assume here that an additional rotational motion stabilizes the toroidal bubble generated by dolphins: with an appropriate azimuthal speed, the centrifugal force stabilizes the Rayleigh-Plateau instability. The stability analysis for a swirling jet has been already studied by Hocking [START_REF] Hocking | The stability of a column of rotating liquid[END_REF]. In the case of an inviscid fluid column with no external fluid and subjected to rotation of constant angular rotation rate Ω 0 (or equivalently of constant vorticity 2Ω 0 ), a stability criterion for azimuthal wavenumber m at k = 0 was obtained using the Weber number based on the velocity Ω 0 a 0 , i.e. W e = ρ (1) a 3 0 Ω 2 0 /σ (not to be confused with the Weber number used in (3) which was based on U 0 and R 0 ). This criterion implies that a mode with k = 0 and m = 0 is stable when W e m(m + 1). The planar mode (k = 0, m = 1) being a displacement mode, it is neutral and stability for (k = 0, m = 0) modes is ensured when W e W e c = 6. In follow-on studies, Gillis [START_REF] Gillis | The stability of a rotating viscous jet[END_REF] provided a general stability criterion for three-dimensional disturbances on a viscous fluid column with no external fluid

W e (ka 0 ) 2 + m 2 -1. (15)
Weidman [START_REF] Weidman | On the instability of inviscid, rigidly rotating immiscible fluids in zero gravity[END_REF] showed how the dominant azimuthal mode for a two-phase axisymmetric rotating system depends on the Reynolds number Re. Finally the linear stability of a uniformly rotating viscous liquid column has been investigated by Kubitschek and Weidman [START_REF] Kubitschek | The effect of viscosity on the stability of a uniformly rotating liquid column in zero gravity[END_REF] and validated experimentally [START_REF] Kubitschek | Helical instability of a rotating viscous liquid jet[END_REF]: the dominating mode depends on Re and on the rotation speed. In the present work, we study how perturbations evolve when superimposed on an axisymmetric air column in an infinite ambient fluid, namely water. We consider the case where surface tension, viscosity, density contrast and centrifugal force act on the instability growth rate, and in which the angular velocity profile U θ (r )/r is not uniform outside the bubble. To the best of our knowledge, this stability problem has not yet been performed.

The governing equation for a columnar flow with rotation

We adopt the model of an infinitely long two-phase capillary jet of radius a 0 . The inner and outer fluid are both incompressible, immiscible and viscous. Using now radius a 0 , [ρ (1) a 3 0 /σ] 1/2 and ρ (1) as characteristic scales for length, time and density, we write the Navier-Stokes equations inside each fluid p = 1 or 2:

∇ • ū(p) = 0 (16) ρ(p) D D t ū(p) = -∇ p(p) + μ(p) Re B ∇ 2 ū(p) , 1
Re B = Oh = µ (1)
ρ (1) σa 0 [START_REF] Popinet | A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations[END_REF] where the Reynolds number Re B of the vortex is the inverse of the Ohnesorge number. Ratios μ(p) = µ (p) /µ (1) , and ρ(p) = ρ (p) /ρ (1) are fixed : μ(2) = 50 and ρ(2) = 10 3 and by definition μ(1) = ρ(1) = 1. In the following, we use only dimensionless quantities except when specified and the bar notation is assumed. Since the inner fluid p = 1 is air and the outer fluid p = 2 is water, it is clear that

Re B = 1 ρ (2) a 0 R 0 Re W e ≈ 300 a 0 ( 18
)
when a 0 is the bubble radius in SI units. When a 0 ∼ 1 cm, this implies Re B ∼ 30, and when a 0 ∼ 10 cm, this implies Re B ∼ 100. Alternatively, if one assumes that a 0 /R 0 ∼ 0.1, for simulations Re = 10000 and W e = 100, one gets Re B = 10; for simulations Re = 10000 and W e = 1, one gets Re B = 100.

C. R. Mécanique, 0000, 1, n 0, 000-000

The Navier-Stokes equations are written in cylindrical coordinates. The radial position of the interface between the two fluids can be described by a function ζ in cylindrical coordinates: r = ζ(θ, x, t ). This interface is characterized by a outward normal n and two tangential vectors t x and t θ :

t x = ∂ζ ∂x , 0, 1 1 + ∂ζ ∂x 2 , n = 1, -1 ζ ∂ζ ∂θ , - ∂ζ ∂x 1 + ∂ζ ∂x 2 + 1 ζ ∂ζ ∂θ 2 , t θ = t x × n, ( 19
)
forming an orthonormal basis (t x , n, t θ ). The velocity is continuous across the interface [[u]] = 0.

The normal velocity must be compatible with the interface motion

Dζ Dt = v r (ζ, θ, x, t ). ( 20
)
The dynamic conditions express that normal stress and tangential shear must be balanced by the surface tension effect:

-p + 2 µ Re B n • e • n = 1 R 1 + 1 R 2 , [[2µt θ • e • n]] = 0 with e (p) i j = 1 2   ∂u (p) i ∂x j + ∂u (p) j ∂x i   . (21)

The rotating base state

Here, we mainly study a basic state u = W θ (r, t )e θ which possesses the feature described in Section 3: it is characterized by a non-uniform angular velocity W θ /r and a non-uniform axial vorticity ω = [W θ /r + ∂W θ /∂r ]e x in phase 2 only, and the interface is circular ζ(x, θ, t ) = 1.

If one neglects the unsteadiness due to bulk viscosity, such profile is a possible solution. In physical terms, this approximation is pertinent since instability occurs a priori much faster than the diffusion of the base state. To be mathematically correct, one introduces a body force to counterbalance the bulk diffusion for the basic state. However viscous effects are still present in the boundary conditions at r = 1:

W (1) θ (1) = W (2) θ (1), (22)
µ (1) ∂W (1) θ ∂r (1) -W (1) θ (1) = µ (2) ∂W (2) θ ∂r (1) -W (2) θ [START_REF] Couder | Observation expérimentale de la turbulence bidimensionnelle dans un film liquide mince[END_REF] . (23) This latter condition can be rewritten as well as ω (2) x (1) = µ (1) µ (2) ω (1) x (1) -2

µ (1) µ (2) -1 W (1) θ (1). ( 24
)
In each phase (p), the basic pressure Π(r ) should satisfy equation ρW 2 θ (r )/r = ∂Π/∂r and a jump condition at r = 1. This leads to a pressure of the following form:

Π (1) (r ) = ρ (1) r 0 [W (1) θ (r )] 2 r dr , Π (2) (r ) = Π (1) (1) + ρ (2) r 1 [W (2) θ (r )] 2 r dr -1. ( 25
)
The vorticity profile is assumed to be uniform in fluid 1 that is W (1) θ (r ) = qr, ω (1) x = 2q for r < 1, where the dimensionless parameter q is equal to the ratio of the dimensional azimuthal velocity at the interface with respect to the characteristic capillary velocity [σ/(ρ (1) a 0 )] 1/2 -this quantity is equal to 0.27/a 0 1/2 when a 0 is expressed in SI units. Note that

q ≡ ρ (1)
ρ (2) q S = ρ (a) ρ (w) q S . (26) Vorticity in phase 1 being uniform, Equation (24) imposes that ω (2) x (1) = 2q. Apart from this constraint, let us assume a vorticity profile ω (2) (r ) = ω (2) x (1) f (r ) in phase 2, where function f is chosen as f (r ) = exp[-(r -1) 2 /δ 2 ] and decays in the radial direction, thus generating a vortex ribbon of radial size δ in phase 2 surrounding the interface. It is easily seen that W (2) θ (r ) = q r + 2 q r r 1 r f (r ) dr , ω (2) x (r ) = 2q f (r ) for r > 1.

This expression, which can be written as well as

W (2) θ (r ) = q r 1 + δ 2 1 -e -(r -1) 2 /δ 2 + δ πerf r -1 δ for r > 1, (27)
is displayed in Figure 16.

The linear instability results

We study the linear stability of the above base flow using standard but cumbersome methods. Indeed, the problem is an eigenvalue problem s( k, m) for given axial k and azimuthal m wavenumbers, that contains viscous diffusion in the bulk as well as on the interface. A code has been specifically developed, that solves coupled Orr-Sommerfeld equations in cylindrical coordinates for the two phases. For obvious reasons, details on the implementation will be published elsewhere, and only results for the air-water columnar bubble are presented here. Figure 17 displays the growthrate s as a function of wavenumber k for the axisymmetric mode m = 0, which is known to be the most unstable mode for q = 0. The presence of rotation is found to stabilize axisymmetric modes, and for q larger than a critical value q c ≈ 0.03, waves become stable for all wavelengths. This observation holds irrespective of the Reynolds number and δ within the studied interval δ ∈ [0.3, 2]. Actually the structure in the corona at r > 1 seems to play a minor role for the axisymmetric mode. for m > 0 and m < 0. The equality exactly holds when q = 0 for symmetry reasons, and the fact that it holds also roughly for q = 0 can be explained by the uniform vorticity inside the bubble. Interestingly, rotation does not influence mode |m| = 1, that remains stable. For sufficiently low Reynolds number and small δ, other azimuthal modes m = 0 for small k can be destabilized. For instance, the mode m = 2 is found unstable for Re B = 10 and δ = 0.3 (Figure 18b). The growth rate decreases as the Reynolds number increases (Figure 20b). This is reminiscent of planar modes k = 0 in Ref. [START_REF] Kubitschek | The effect of viscosity on the stability of a uniformly rotating liquid column in zero gravity[END_REF]. Note that for larger values of δ, this mode is restabilized: for δ = 1, all modes tested are stabilized by rotation, this conclusion holding irrespective of the Reynolds number (Figures 19 and21).

Finally, going back to the bubble ring, the axisymmetric Rayleigh-Plateau modes are stabilized by rotation. Indeed, from the simulations in Section 2, the pertinent values of q (Figure 13 and (26)) are found to be q ∼ 10 -3/2 × 30 ∼ 0.94. When a stretching factor is applied, q generally increases (Equations ( 14) and (26)) and is thus larger than the stability threshold of q c ≈ 0.03.

Conclusions

This paper covers three different aspects of the dynamics of toroidal bubbles such as those produced by scuba divers, beluga whales or dolphins. During the initial production of such bubble, vorticity trailing from the nozzle or from the opening rolls up into a toroidal vortex ring able to trap gas in its core; to do so, it has to overcome surface tension which would favor the formation of a single spherical bubble. Simulations show that vortices shed from the nozzle or generated at the bubble interface itself can influence the formation dynamics in a non trivial way. Once formed, the bubble ring may be stretched by an appropriate flow so as to increase its length. A theoretical model where curvature effects are overlooked shows that, during the stretching process, the vorticity in the gas tends to become uniform, the cross-section tends to circular while the core radius is highly reduced and the swirl is enhanced. This latter finding is important since a linear stability study of a gas columnar vortex surrounded by water eventually shows that the dominant axisymmetric Plateau instability is stabilized as swirl increases above a C. R. Mécanique, 0000, 1, n 0, 000-000 well-defined critical value. Put all together, these pieces of investigation give some serious clues to physically understand the surprising stability of swirling bubble rings.

Figure 1 .

 1 Figure 1. (a) Scheme of the numerical nozzle. (b) The whole numerical domain. In the yellow region, the mesh size is ∆ c = 0.125R 0 and in the red region, ∆ c = 0.0078R 0 . (c) The three initial cases: the flow is at rest and a gas bubble (in red) of length L B is contained in the nozzle. From top to bottom: L B = 0.6R 0 , 1.2R 0 and 2.4R 0 .

Figure 2 .

 2 Figure 2. Simulations at Re = 5000, W e = 3 and r B = 0.6. From left to right, the snapshots correspond to times t = 2, 3.2 and 6.4. The upper half displays phases and interfaces (red is gas and blue water) and the lower half the vorticity field. The vorticity value is colored with a maximum scale of ω = 3.75, which is one fifth of the absolute vorticity maximum at the wall boundary layer (see (4)).

Figure 3 .

 3 Figure 3. Simulations at Re = 5000, W e = 3 and r B = 2.4. From left to right, first row: t = 4, 6, 8; second row: t = 8.4, 10 and 16. Visualizations are similar to Figure 2.

Figure 4 .

 4 Figure 4. Simulations at Re = 5000, W e = 100 and r B = 0.6. From left to right: t = 2.8, 4.8, 7.2 and 15.2. Visualizations are similar to Figure 2.

Figure 5 .

 5 Figure 5. Simulations at Re = 5000, W e = 100 and r B = 2.4. From left to right, first row: t = 4.8, 7.2, 8, 11.2; second row: t = 13.6 and 19.6. Visualizations are similar to Figure 2.

Figure 6 .

 6 Figure 6. Simulations at Re = 5000, W e = 100 and r B = 8. From left to right: t = 14, 20, 24 and 40. Visualizations are similar to Figure 2.

Figure 7 .

 7 Figure 7. Simulations at Re = 5000, W e = 20 with r B = 2.4. From left to right, first row: t = 4.4, 7.2, 9.6; second row: t = 14.4, 20 and 24. Visualizations are similar to Figure 2.

Figure 8 .

 8 Figure 8. Simulations at Re = 5000, W e = 20 with r B = 1.2. From left to right, first row: t = 4.4, 5.6, 8, 9.6; second row: t = 20, 30. Visualizations are similar to Figure 2.

Figure 9 .

 9 Figure 9. Simulations at Re = 5000, W e = 20 with r B = 0.6. From left to right: t = 2, 3.6, 4.8 and 8. Visualizations are similar to Figure 2.

Figure 10 .

 10 Figure 10. Simulations with Re = 10000 and r B = 2.4. From left to right, top row: W e = 3 at t = 8, 10, 14.6; bottom row: W e = 100 at t = 8, 10 and 14.6. Visualizations are similar to Figure 2.

Figure 11 .

 11 Figure 11. Simulations with Re = 5000, r B = 4. Vorticity field near the bubble for (a) W e = 20 at t = 20; (b) W e = 100 at t = 26.4.

Figure 12 .

 12 Figure 12. Profiles of vorticity ω θ (x, r ) at Re = 5000, W e = 100, r B = 4 and t = 17.6 (corresponding to the second case of Figure 11) on two orthogonal lines passing through the center of section S: (a) along the x direction, (b) along the r direction. The red dots are located inside the gas bubble, blue dots are outside.

Figure 13 .

 13 Figure 13. Swirl number q S as a function of nondimensional time t . Simulations at Re = 5000, W e = 100 with (a) r B = 4, (b) r B = 5.

Figure 14 .

 14 Figure14. Evolution of the vorticity field ω x (χ, η, τ) near the bubble, starting from the initial condition corresponding to the second case in Figure11and evolving through the two-dimensional advection-diffusion equation (13), represented (left to right and top to bottom) at times τ = 0, 8, 16, 24, 32 and 40.

Figure 15 .

 15 Figure 15. rate s as a function of wavenumber k for axisymmetric modes m = 0. (a) Inviscid column of liquid 1 surrounded by a liquid 2 at various density ratios ρ (2) /ρ (1) = 10 3 , 10 2 , 10, 1, 10 -1 and 10 -3 . (b) Viscous air column surrounded by water at various Ohnesorge numbers Oh between 10 -3 and 10 -1 with µ (2) /µ (1) = 50. The inviscid curve almost concides with that of case Oh = 10 -3 .

Figure 16 .

 16 Figure 16. (a) Normalized vorticity ω x (r )/(2q) and (b) normalized velocity W θ (r )/q profiles for the base state with ρ(2) /ρ(1) = 10 3 , µ(2) /µ(1) = 50 for δ = 0.3 (left), δ = 1 (center) or δ = 2 (right).

Figure 17 .

 17 Figure 17.Growth rate s as a function of wavenumber k for axisymmetric mode m = 0 at various swirl numbers q = 0, 0.01, 0.02, 0.03. Reynolds number is (a) Re B = 10, (b) Re B = 100. Viscous air column surrounded by water with δ = 0.3 (left) or δ = 2 (right), with ρ (2) /ρ(1) = 10 3 and µ(2) /µ(1) = 50.

Figure 18 .

 18 Figure 18. Growth rate s as a function of wavenumber k for a viscous column of air surrounded by water at δ = 0.3, Re B = 10 and various swirl q = 0.2, 0.4, 0.6. (a) m = ±1; (b) m = ±2; (c) m = ±3. Positive m: symbols; negative m: bold line.

Figure 19 .

 19 Figure 19. Same as Figure 18 but for δ = 1 (Re B = 10).

Figures 18 -

 18 Figures 18-21 provide linear instability results for non axisymmetric modes m ± 1, ±2 and ±3. The first observation is that, at given k and |m| values, the growth rate is found almost the same

Figure 20 .

 20 Figure 20. Same as Figure 18 but for δ = 0.3 and Re B = 100.

Figure 21 .

 21 Figure 21. Same as Figure 18 but for δ = 1 and Re B = 100.

Table 1 .

 1 Estimated typical Reynolds and Weber numbers for scuba divers, beluga whales or dolphins

C. R. Mécanique, 0000, 1, n 0, 000-000