Pierre Degond
email: pdegond@imperial.ac.uk

Hui Yu

Self-organized Hydrodynamics in an Annular Domain: Modal Analysis and Nonlinear Effects

Keywords: Collective dynamics, Self-organization, emergence, fluid model, hydrodynamic limit, symmetry-breaking, alignment interaction, polarized motion, spectral analysis, relaxation model, splitting scheme, conservative form, nonlinear mode-coupling AMS Subject classification: 35L60, 35L65, 35P10, 35Q80, 82C22, 82C70, 82C80, 92D50

The Self-Organized Hydrodynamics model of collective behavior is studied on an annular domain. A modal analysis of the linearized model around a perfectly polarized steady-state is conducted. It shows that the model has only pure imaginary modes in countable number and is hence stable. Numerical computations of the low-order modes are provided. The fully non-linear model is numerically solved and nonlinear mode-coupling is then analyzed. Finally, the efficiency of the modal decomposition to analyze the complex features of the nonlinear model is demonstrated.

Introduction

Self-organized collective dynamics is ubiquitous in the living world and emerges at all possible scales, from cell assemblies [START_REF] Yamao | Multi-Cellular Logistics of Collective Cell Migration[END_REF] to animal groups [START_REF] Tunstrom | Collective States, Multistability and Transitional Behavior in Schooling Fish[END_REF]. Collective motion happens when thousands of moving individual entities coordinate with each other through local interactions such as attraction and alignment. As a result, large-scale structures of typical sizes exceeding the inter-individual distances by several orders of magnitude are formed. One of the key questions is to understand how these self-organized structures spontaneously emerge from local interactions without the intervention of any leader. With this aim, Individual-Based Models (IBM), i.e. models that describe the behavior of each individual agent have been investigated [START_REF] Aoki | A simulation study on the schooling mechanism in fish[END_REF][START_REF] Chaté | Collective motion of self-propelled particles interacting without cohesion[END_REF][START_REF] Chuang | State transitions and the continuum limit for a 2D interacting, self-propelled particle system[END_REF][START_REF] Couzin | Collective Memory and Spatial Sorting in Animal Groups[END_REF][START_REF] Cucker | Emergent behavior in fmocks[END_REF][START_REF] Lukeman | Inferring individual rules from collective behavior[END_REF][START_REF] Peruani | Nonequilibrium clustering of self-propelled rods[END_REF][START_REF] Rejniak | Hybrid models of tumor growth[END_REF]. They consist of large systems of ordinary or stochastic differential equations the numerical resolution of which is computationally intensive. To describe large-scale structures coarse-grained models such as Fluid Models (FM) are needed. FM describe the dynamics of average quantities such as the mean density or mean velocity of the individuals [START_REF] Bertozzi | Blow-up in multidimensional aggregation equations with mildly singular interaction kernels[END_REF][START_REF] Toner | Long-range order in a two-dimensional dynamical XY model: how birds fly together[END_REF][START_REF] Toner | Hydrodynamics and phases of flocks[END_REF][START_REF] Topaz | A nonlocal continuum model for biological aggregation[END_REF]. Attempts to derive FM from IBM of collective motion can be found in [START_REF] Ratushnaya | Collective behavior of self-propelling particles with kinematic constraints: the relation between the discrete and the continuous description[END_REF]. An intermediate step in the hierarchy of models consist of kinetic models (KM) [START_REF] Bellomo | On the mathematical theory of the dynamics of swarms viewed as complex systems[END_REF][START_REF] Bertin | Boltzmann and hydrodynamic description for self-propelled particles[END_REF][START_REF] Bertin | Hydrodynamic equations for self-propelled particles: Microscopic derivation and stability analysis[END_REF][START_REF] Bolley | Mean-field limit for the stochastic Vicsek model[END_REF] which are Partial Differential Equations (PDE) describing the evolution of the probability density of the particles in phase-space. FM can be obtained as singular limits of the KM under the hypothesis that the individual scales are much smaller than the system scales. This PDE-based derivation of FM is referred to as the 'Hydrodynamic Limit'.

In [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF], the hydrodynamic limit of the Vicsek IBM [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] has been performed using an intermediate kinetic description [START_REF] Bolley | Mean-field limit for the stochastic Vicsek model[END_REF][START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF]. The Vicsek IBM describes a noisy system of selfpropelled particles interacting through local alignment. In [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF], it has been shown that the absence of conservation laws (such as momentum conservation) resulting from selfpropulsion can be overcome by introducing the new "Generalized Collision Invariant" concept. The resulting model, referred to as the "Self-Organized Hydrodynamics (SOH)" is written:

∂ t ρ + c 1 ∇ • (ρΩ) = 0, (1.1) ρ [∂ t Ω + c 2 (Ω • ∇)Ω] + ΘP Ω ⊥ ∇ρ = 0, (1.2) |Ω| = 1, (1.3)
where ρ(x, t) ≥ 0 and Ω(x, t) ∈ R d are the density and the orientation of the mean velocity of the particles, c 1 > 0, c 2 ∈ R and Θ > 0 are given parameters, and d is the spatial dimension. We let P Ω ⊥ = Id -Ω ⊗ Ω be the projection matrix onto the plane orthogonal to Ω. This model resembles the usual isothermal gas dynamics equations. Eq. (1.1) is the continuity equation expressing the conservation of mass. Eq. (1.2) describes how the velocity orientation evolves under transport by the flow (the second term) and the pressure gradient (the third term, where Θ is related to the noise in the underlying IBM and has the interpretation of a temperature). However, there are important differences, which arise from the fact the Ω is not a true velocity but the velocity direction, i.e. it is a vector of unit norm (which is expressed by (1.3)). To preserve this geometrical constraint, the pressure gradient has to be projected onto the normal to Ω, which is the reason for the presence of P Ω ⊥ . Other differences stem from the allowed discrepancy between the two constants c 1 and c 2 . While c 1 fixes the material velocity to c 1 Ω, the constant c 2 describes how Ω is transported. This discrepancy originates from the lack of Galilean invariance of the underlying IBM, itself resulting from self-propulsion [START_REF] Tu | Sound waves and the absence of Galilean invariance in flocks[END_REF]. This model has been extended into several directions [START_REF] Degond | Hydrodynamics of self-alignment interactions with precession and derivation of the Landau-Lifschitz-Gilbert equation[END_REF][START_REF] Degond | Hydrodynamic models of selforganized dynamics: derivation and existence theory[END_REF][START_REF] Degond | Hydrodynamics of the Kuramoto-Vicsek model of rotating self-propelled particles[END_REF][START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF][START_REF] Degond | Diffusion in a continuum model of self-propelled particles with alignment interaction[END_REF] and a rigorous existence result is established in [START_REF] Degond | Hydrodynamic models of selforganized dynamics: derivation and existence theory[END_REF].

This paper is devoted to the study of the SOH model in an annular domain. Annular geometries allow for simple observations of symmetry-breaking transitions induced by collective motion. When a transition from disordered to collective motion occurs, the system is set into a collective rotation in either clockwise or counter-clockwise directions. Annular geometries are a traditional design for salmon cages in sea farms [START_REF] Føre | Modelling of atlantic salmon (Salmo salar L.) behaviour in sea-cages: a lagrangian approach[END_REF][START_REF] Johansson | The Interaction between Water Currents and Salmon Swimming Behaviour in Sea Cages[END_REF] and for experiments with locusts [START_REF] Buhl | From disorder to order in marching locusts[END_REF][START_REF] Erban | From individual to collective behaviour of coupled velocity jump processes: a locust example[END_REF], pedestrians [START_REF] Moussaïd | Traffic Instabilities in Self-organized Pedestrian Crowds[END_REF] or sperm-cell dynamics [START_REF] Creppy | Dispositif de traitement d'un échantillon d'un fluide biologique actif[END_REF]. In all these examples, a polarized motion in one direction is observed. In the sperm-cell experiments, the observation of turbulent structures that superimpose to collective rotation motivates the present work. In pure semen, sperm-cells are mostly interacting through volume exclusion. But volume exclusion interactions of rod-like self-propelled particles result in alignment [START_REF] Peruani | Nonequilibrium clustering of self-propelled rods[END_REF]. This legitimates the use of the Vicsek model [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] and of its fluid counterpart, the SOH Model [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF], as models of collective sperm-cell dynamics. The Vicsek model in annular geometry has been shown to exhibit polarized motion in [START_REF] Czirok | Collective behavior of interacting self-propelled particles[END_REF]. Here, we focus on the SOH model and study its normal modes in annular geometry in both the linear and nonlinear regimes.

We first study the linear modes of the SOH model around a perfectly polarized steadystate in Sec. 2. One of the main results of this paper is that these modes are pure imaginary (and thus, stable) and form a countable set. In Sec. 3, we compute the eigenmodes and eigenfunctions numerically and investigate how the eigenmodes depend on the geometry of the annulus and on the parameters of the model. We then turn towards the nonlinear model with the aims of (i) validating the linear analysis for small perturbations, (ii) investigating how the nonlinearity of the model affects the modal decomposition of the solution and (iii) demonstrating the capabilities of the modal decomposition to analyze the complex features of the nonlinear model. In future work, the modal decomposition will be used to calibrate the model coefficients against experimental data. We first develop the scheme in Sec. 4 and then compare the results for the linear and nonlinear models in Sec. 5. Finally we draw conclusions and perspecives in Sec. 6.

Linear Modes of the SOH Model in Polar Coordinates

2.1 The SOH model in polar coordinates and perfectly polarized steady-states

Consider the SOH model (1.1)-(1.3) in a two-dimensional annular domain D = {x ∈ R 2 | |x| ∈ (R 1 , R 2 )}. We introduce polar coordinates (r, θ) ∈ (R 1 , R 2 ) × [0, 2π
] where r = |x| and θ is the angle between x and a reference direction. We denote by (e r , e θ ) the local basis associated to polar coordinates, i.e. e r = x/|x| = (cos θ, sin θ) and e θ = e ⊥ r = (-sin θ, cos θ) where the exponent ⊥ indicates a rotation by an angle +π/2. Then, we let ρ = ρ(r, θ, t) and Ω = Ω(r, θ, t) = cos φ(r, θ, t) e r + sin φ(r, θ, t) e θ , where φ(r, θ, t) represents the angle between e r and Ω. We recall that the constants c 1 , c 2 and Θ are such that c 1 > 0, c 2 ∈ R, Θ > 0. For notational convenience, we introduce

α = c 2 Θ , (2.4)
and we note that α is of the same sign as c 2 and that c 2 /α = 1/Θ > 0. After easy algebra, the SOH model (1.1)-(1.3) is equivalent to the following system for ρ(r, θ, t) and φ(r, θ, t) with (r, θ) ∈ (R 1 , R 2 ) × [0, 2π] and t > 0,

∂ t ρ + c 1 r ∂ ∂r (rρ cos φ) + ∂ ∂θ (ρ sin φ) = 0, (2.5)
ρ ∂ t φ + c 2 cos φ ∂φ ∂r + sin φ r ∂φ ∂θ + sin φ r + Θ cos φ r ∂ρ ∂θ -sin φ ∂ρ ∂r = 0, (2.6)
subject to the boundary conditions

φ(R 1 , θ, t) = φ(R 2 , θ, t) = ± π 2 ,
ρ and φ periodic in θ.

(2.7)

The first boundary condition (2.7) imposes a tangential flow to the boundary ∂D and consequently ensures that there is no mass flow across this boundary. Now, we look for perfectly polarized steady states of the above system, i.e. steady states of the form (ρ s , φ s ) where ρ s is independent of θ and φ s = -π/2 in the whole domain (We have arbitrarily chosen a rotation in the clockwise direction but of course, the results would be the same, mutatis mutandis, with the opposite choice). We have the Lemma 2.1. The perfectly polarized steady-states form a one-parameter family of solutions given by

ρ s (r) = ρ * s r α , φ s (r, θ) = - π 2 ,
where α is given by (2.4) and ρ * s > 0 is any positive constant. Proof. Inserting φ s = -π 2 into (2.6) gives -c 2 r ρ s +Θ ∂ρs ∂r = 0. Therefore, there exists

ρ * s > 0 such that ρ s (r) = ρ * s r c 2 Θ = ρ * s r α .

Linearization about perfectly polarized steady-states

Next we study the linearization of (2.5), (2.6) about a perfectly polarized steady-state (ρ s , φ s ). Given ε > 0, a linear perturbation (ρ, φ) is given by

ρ = ρ s + ερ + O(ε 2 ), φ = φ s + ε φ + O(ε 2 ).
Expanding System (2.5), (2.6) about (ρ s , φ s ) and dropping terms of order ε 2 or higher, we deduce that the system satisfied by (ρ, φ) is given by:

∂ ∂t ρ φ + L ρ φ = 0, (2.8)
with

L = 0 c 1 ρ s c 2 αρs 0 ∂ ∂r + -c 1 r 0 0 -c 2 r ∂ ∂θ + 0 c 1 r (1 + α)ρ s -c 2 rρs 0 , (2.9)
supplemented with the boundary conditions:

φ(R 1 , θ, t) = φ(R 2 , θ, t) = 0, R 2 R 1 2π 0 ρr drdθ = 0, (2.10)
and ρ, φ periodic in θ. These bounday conditions are inherited from (2.7). The second Eq. in (2.10) is a normalization condition whose physical significance is that we are perturbing the steady-state keeping the total particle mass in the system fixed.

Looking for solutions (ρ, φ) in separation of variables form:

ρ(r, θ, t) = e λt ρ λ (r, θ), φ(r, θ, t) = e λt φ λ (r, θ),
we deduce that (ρ λ , φ λ ) must satisfy the following spectral problem:

(L + λ I) ρ λ φ λ = 0, (2.11)
supplemented with the boundary conditions (2.10), where I is the identity matrix. We now consider the decomposition of (ρ λ , φ λ ) into Fourier series, i.e. Then, (ρ n , φn ) satisfies the following spectral problem:

(L n + λ I) ρn φn = 0, (2.12)
with

L n = 0 c 1 ρ s c 2 αρs 0 ∂ ∂r + -in c 1 r (1 + α)c 1 ρs r -c 2 rρs -in c 2 r , (2.13)
supplemented with the boundary conditions:

φn (R 1 ) = φn (R 2 ) = 0, R 2 R 1 ρn r dr = 0. (2.14)
We first study the existence of non-trivial solutions (λ, ρn , φn ) to this spectral problem.

Study of the spectral problem for L n

We first prove the following Lemma 2.2. All the eigenvalues λ of (2.12) are pure imaginary.

Proof. We recall that ρ s (r) = ρ * s r α . We introduce the transformation:

ρ n (r) = r -α ρn (r), φ n (r) = ρ * s r α+1 φn (r), (2.15)
and find

∂ρ n ∂r + α r α+1 λ c 2 - in r φ n = 0, (2.16
)

∂φ n ∂r + λ c 1 - in r r α+1 ρ n = 0, (2.17)
subject to the boundary conditions:

φ n (R 1 ) = φ n (R 2 ) = 0, R 1 R 1 ρ r α+1 dr = 0. (2.

18)

Let λ = µ + iν where µ denotes the real part of λ and ν its imaginary part. Assume that µ = 0. We divide (2.17

) by λ c 1 -in r r α+1 = µ c 1 + i ν c 1 -n r
r α+1 = 0 and use the first equation (2.16) to get (remembering (2.4)):

∂ ∂r   1 r α+1 µ c 1 + i ν c 1 -n r ∂φ n ∂r   - 1 r α+1 µ Θ + i ν Θ - nα r φ n = 0. (2.19)
Multiplying (2.19) by φn (the complex conjugate of φ n ), integrating with respect to r, using the boundary conditions (2.18) and taking the real part of the so-obtained expression, we get

R 2 R 1 µ c 1 r α+1 µ c 1 2 + ν c 1 -n r 2 ∂φ n ∂r (r) 2 dr + R 2 R 1 1 r α+1 µ Θ |φ n (r)| 2 dr = 0.
Since Θ > 0 and µ = 0, we have φ n = 0, which shows that there cannot exist a non-trivial solution of the spectral problem when µ = 0.

We now determine the eigenvalues λ = iν, ν ∈ R of (2.12). Dropping the index n for simplicity, we introduce the following transformation:

u = c 2 c 1 α r α+1 2 ρ, v = i r α+1 2
φ.

(2.20)

From (2.16), (2.17), (2.18), (u, v) satisfies the spectral problem:

(A n + νI) u v = 0, (2.21)
with the operator A n acting on (u, v) ∈ L 2 (R 1 , R 2 ) 2 defined by

A n u v =     c 1 n r u - c 1 c 2 α 1 r α+1 2 ∂ ∂r r α+1 2 v c 1 c 2 α r α+1 2 ∂ ∂r u r α+1 2 + c 2 n r v     ,
with domain

D(A n ) = (u, v) ∈ H 1 (R 1 , R 2 ) 2 , v(R 1 ) = v(R 2 ) = 0 .
Note that ν ∈ R and that A n has real-valued coefficients. Without loss of generality, we look for real-valued eigenfunctions (u, v). For this operator, we have the Theorem 2.3. The spectrum of A n consists of a countable set of eigenvalues (ν nm ) m∈N , ν nm ∈ R associated to a complete orthonormal system of

L 2 (R 1 , R 2 ) 2 of eigenfunctions (u nm , v nm ) m∈N . Furthermore, |ν nm | → ∞ as m → ∞.
Proof. We first assume that n = 0 and drop the subindex n of A n for simplicity. We will show that there exists η ∈ R such that the resolvant R η = (A + ηI) -1 exists and is compact in

L 2 (R 1 , R 2 ) 2 . For this purpose, we consider (h, g) ∈ L 2 (R 1 , R 2 ) 2 and look for a solution (u η , v η ) ∈ D(A) of (A + ηI)(u η , v η ) = (h, g), i.e., 1 + ηr c 1 n c 1 n r u η - c 1 c 2 α 1 r α+1 2 ∂ ∂r r α+1 2 v η = h, (2.22)
c 1 c 2 α r α+1 2 ∂ ∂r u η r α+1 2 + 1 + ηr c 2 n c 2 n r v = g. (2.23) We take |η| < η 0 := c 1 |n| R 1 in such a way that 1 + ηr c 1 n > 0, ∀r ∈ [R 1 , R 2 ]. Multiplying (2.22) by (1 + ηr c 1 n ) -1 r -α-1
2 , taking its derivative with respect to r and using (2.23), we deduce that v η satisfies:

- ∂ ∂r 1 r α 1 + ηr c 1 n ∂ ∂r r α+1 2 v η - αn 2 r α+3 2 1 + ηr c 2 n v η = hη , (2.24)
with the boundary conditions

v η (R 1 ) = v η (R 2 ) = 0, where hη = α c 1 c 2 ∂ ∂r h r α-1 2 1 + ηr c 1 n - αn c 2 r α+1 2 g. Note that hη ∈ H -1 (R 1 , R 2 )
. Now the problem consists of showing the existence of a unique weak solution v η ∈ H 1 0 (R 1 , R 2 ) to (2.24), i.e. to the variational formulation:

a η (v η , ṽ) = hη , r α+1 2 ṽ H -1 ,H 1 0 , ∀ṽ ∈ H 1 0 (R 1 , R 2 ), (2.25)
with

a η (v η , ṽ) = R 2 R 1 1 r α 1 + ηr c 1 n ∂ ∂r r α+1 2 v η ∂ ∂r r α+1 2 ṽ dr (2.26) - R 2 R 1 αn 2 r 1 + ηr c 2 n v η ṽ dr,
and where the brackets at the right-hand side of (2.25) denote duality between the distribution hη ∈ H -1 and the function r

α+1 2
ṽ ∈ H 1 0 . We introduce σ > 0. Choosing σ large enough and |η| < η 0 there exists C, C ′ > 0 such that

σ - αn 2 r α+3 2 1 + ηr c 2 n ≥ C > 0, 1 r α 1 + ηr c 1 n > C ′ > 0, ∀r ∈ [R 1 , R 2 ].
Therefore, the bilinear form

a η σ (v, ṽ) = a η (v, ṽ) + (v, ṽ r α+1 2 ) where (•, •) is the usual inner product in L 2 (R 1 , R 2 ) is coercive on H 1 0 (R 1 , R 2 )
. Consequently, by the Lax-Milgram theorem, the variational formulation

a η σ (v, ṽ) = ℓ, r α+1 2 ṽ H -1 ,H 1 0 , ∀ṽ ∈ H 1 0 (R 1 , R 2 ),
has a unique solution for any ℓ ∈ H -1 (R 1 , R 2 ) and the dependence of v upon ℓ is continuous. This defines a continuous linear mapping T η σ :

H -1 (R 1 , R 2 ) → H 1 0 (R 1 , R 2 ), ℓ → v. Then v η ∈ H 1 0 (R 1 , R 2 ) is a solution of (2.25) if and only if v η = T η σ ( h + σv η ), i.e. (I -σT η σ )v η = T η σ hη . (2.27)
We note that

T η σ hη ∈ H 1 0 (R 1 , R 2 ) ⊂ L 2 (R 1 , R 2 ) and that T η σ , restricted to L 2 (R 1 , R 2 ), is a bounded operator from L 2 (R 1 , R 2 ) to H 1 0 (R 1 , R 2 ). After composition with the canonical imbedding of H 1 0 (R 1 , R 2 ) into L 2 (R 1 , R 2 ) (still denoted by T η σ ), T η σ is a compact operator on L 2 (R 1 , R 2 ). Therefore I -σT η
σ is a Fredholm operator. In addition, T η σ is self-adjoint because the bilinear form a η σ is symmetric. Thanks to the Fredholm alternative, we have

Im(I -σT η σ ) = Ker(I -σT η σ ) ⊥ ,
where Im denotes the range and Ker denotes the null space of an operator. Suppose that there exists η such that |η| < η 0 and Ker(I -σT η σ ) = {0}. Then, (I -σT η σ ) is invertible and there exists a unique solution

v η ∈ H 1 0 (R 1 , R 2 ) to (2.27
), or equivalently, to (2.24). Defining u η by (2.22) (remember that we suppose n = 0), then

u η ∈ L 2 (R 1 , R 2 ) since v η ∈ H 1 (R 1 , R 2 ) and h ∈ L 2 (R 1 , R 2 ). But, since v η satisfies (2.24) in the distributional sense, u η satisfies (2.23) in the distributional sense. From the facts that v η and g both belong to L 2 (R 1 , R 2 ), we get that u η ∈ H 1 (R 1 , R 2 ). Therefore, (u η , v η ) ∈ D(A) and by (2.22), (2.23), it satisfies (A+ηI)(u η , v η ) = (h, g). This shows that A + ηI is invertible. Furthermore, since D(A) is compactly imbedded into L 2 (R 1 , R 2 ) 2 and that R η = (A + ηI) -1 is a continuous linear map from L 2 (R 1 , R 2 ) 2 into D(A), the map R η is compact as an operator of L 2 (R 1 , R 2 ) 2 .
To prove that there exists η such that |η| < η 0 and Ker(I -σT η σ ) = {0}, we proceed by contradiction. We suppose that for all such η there exists a non-trivial

v η ∈ H 0 1 (R 1 , R 2 )
such that (I -σT η σ )v η = 0. Equivalently, Eq. (2.25) with right-hand side hη = 0 has a non-trivial solution v η ∈ H 1 0 (R 1 , R 2 ), which means that v η is an eigenvector for the eigenvalue 0 of the variational spectral problem: " to find λ ∈ R and

v λ ∈ H 1 0 (R 1 , R 2 ), v λ = 0, such that a η (v λ , ṽ) = λ(v λ , ṽ r α+1 2 ), ∀ṽ ∈ H 1 0 (R 1 , R 2 ). ′′
From the classical spectral theory of elliptic operators [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF], we know that the eigenvalues of this problem are isolated. Furthermore, 0 is a simple eigenvalue. Indeed, Eq. (2.24) is a linear second order differential equation. For a given η, consider two solutions

v 1 , v 2 in H 0 1 (R 1 , R 2 ) of (2.24) associated to hη = 0. The Wronskian v 1 ∂ r v 2 -v 2 ∂ r v 1 is zero because both v 1
and v 2 vanish at the boundaries. Therefore, v 1 and v 2 are linearly dependent and consequently the dimension of the associated eigenvectors is 1. We realize that the coefficients of a η given by (2.26) are analytic functions of η ∈ [-η 0 , η 0 ]. Then, from classical spectral theory again [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF], one can define an analytic branch of non-zero solutions

η ∈ [-η 0 , η 0 ] → v η ∈ H 0 1 (R 1 , R 2 )
. Now, from (2.25) with right-hand side hη = 0, it follows that for such v η , we have a η (v η , v η ) = 0. Taking the derivative of this identity with respect to η at η = 0, and using the fact that a η is a symmetric bilinear form, we get:

da η dη η=0 (v 0 , v 0 ) + 2 a 0 v 0 , dv η dη η=0 = 0.
Now, since v 0 is a variational solution of (2.24) for η = 0 with zero right-hand side, the second term is identically zero. Computing the first term, we get

-n c 1 n 2 R 2 R 1 1 r α-1 ∂ ∂r r α+1 2 v 0 2 dr + 1 Θ R 2 R 1 |v 0 | 2 dr = 0
The quantity inside the parentheses is a nonegative quantity which can only be 0 if v 0 is identically zero, which contridicts the hypothesis that v 0 is a non-trivial solution. This shows the contradiction and proves that there exists η ∈ R small enough such that Ker(I -σT η σ ) = {0}. In the case n = 0, it is an easy matter to see that the above proof can be reproduced or alternately, one can invoke directly the spectral theory of elliptic operators. Details are left to the reader. Now, for all n ∈ Z, there exists η n ∈ R such that R n = (A n + η n I) -1 exists, and is a compact self-adjoint operator of L 2 (R 1 , R 2 ) 2 . By the spectral theorem for compact selfadjoint operators, there exists a Hilbert basis (u nm , v nm ) m≥0 of L 2 (R 1 , R 2 ) 2 and a sequence (τ nm ) m≥0 of real numbers such that τ nm → 0 as m → ∞ and such that (u nm , v nm ) is an eigenfunction of R n associated to the eigenvalue τ nm . Then, (u nm , v nm ) m≥0 is a Hilbert basis in L 2 (R 1 , R 2 ) 2 of eigenfunctions of A n associated to the sequence of eigenvalues (ν nm ) m≥0 with ν nm = 1 τnm -η n . We have |ν nm | → ∞ as m → ∞, which concludes the proof.

We now come back to the original spectral problem (2.11). We define

ρnm = c 1 α c 2 r α-1 2 u nm , ψnm = 1 ρ * s r -(α+1) 2 v nm , (2.28)
where (u nm , v nm ) m≥0 is the Hilbert basis of eigenfunctions of A n found at Theorem 2.3, and (ν mn ) m≥0 is the associated sequence of eigenvalues. Thanks to the change of functions (2.15), (2.20), (ρ nm , i ψnm ) is a Hilbert basis of eigenfunctions of L n and (iν mn ) m≥0 is the associated sequence of eigenvalues. The system (ρ nm , i ψnm ) is orthonormal for the inner product

ρ φ , ρ φ = R 2 R 1 Θ c 1 r -(α-1) ρ(r) ρ(r) + (ρ * s ) 2 r α+1 φ(r) φ(r) dr. (2.29)
Furthermore, we have the following easy Lemma (whose proof is left to the reader):

Lemma 2.4. Let iν be an eigenvalue of L n associated to the eigenvector (ρ, i ψ), then -iν is an eigenvalue of L -n associated to the eigenvector (ρ, -i ψ).

As a consequence of this Lemma the eigenvalues for n = 0 come in opposite pairs and we number them such that ν 0 2m = -ν 0 2m-1 . Therefore, the sequence of eigenvalues of L 0 is (-iν 0 2m , iν 0 2m ) m≥1 . We note that 0 is not an eigenvalue of L 0 .

The spectral problem for L and resolution of the initial value problem

We now turn to the operator L defined on L 2 ((R 1 , R 2 ) × (0, 2π)) 2 by (2.9) with domain

D(L) = {(ρ, φ) ∈ H 1 ((R 1 , R 2 ) × (0, 2π)) 2 | φ(R 1 , θ) = φ(R 2 , θ) = 0, a.e. θ ∈ (0, 2π), (ρ, φ)(r, 0) = (ρ, φ)(r, 2π), a.e. r ∈ (R 1 , R 2 )}.
Using Theorem 2.3 and Lemma 2.4, we can state the following theorem (the proof of which is immediate and left to the reader):

Theorem 2.5. The spectrum of L, Spec L is discrete, and consists of

Spec L = n≥1,m≥0
{iν nm , -iν nm } m≥1 {iν 0 2m , -iν 0 2m } , associated to the following basis of eigenvectors n≥1,m≥0

{(ρ nm , i ψnm )e inθ , (ρ nm , -i ψnm )e -inθ } m≥1 {(ρ 0 2m , i ψ0 2m ), (ρ 0 2m , -i ψ0 2m )} , which is a Hilbert basis in L 2 ((R 1 , R 2 ) × (0, 2π)) 2 for the inner product ρ φ , ρ φ (2.30) = 1 2π 2π 0 R 2 R 1 Θ c 1 r -(α-1) ρ(r, θ) ρ(r, θ) + (ρ * s ) 2 r α+1 φ(r, θ) φ(r, θ) dr dθ.
From this theorem, we have the immediate Theorem 2.6. Let (ρ, φ)(r, θ, t) be the solution of the linearized model (2.8) with initial condition (ρ I , φI ) ∈ L 2 ((R 1 , R 2 ) × (0, 2π)) 2 . By standard semigroup theory, this solution belongs to

C 0 ([0, T ], L 2 ((R 1 , R 2 ) × (0, 2π)) 2 ∩ L 2 ([0, T ], D(L))
, for all time horizon T ∈ R.

Additionally, we assume that (ρ I , φI ) is real-valued. Then, (ρ, φ)(r, θ, t) can be expressed as:

ρ φ (r, θ, t) = n≥1, m≥0
k nm ρnm (r) cos(nθ + ν nm t + ϕ nm ) -ψnm (r) sin(nθ + ν nm t + ϕ nm )

+ m≥1 k 0 2m ρ0 2m (r) cos(ν 0 2m t + ϕ 0 2m ) -ψ0 2m (r) sin(ν 0 2m t + ϕ 0 2m ) ,
where the series converges in L 2 ((R 1 , R 2 ) × (0, 2π)) 2 and where k nm and ϕ nm are given, for all (n, m) with (n ≥ 1 and m ≥ 0) or (n = 0 and m ≥ 1), by:

ρ I φ I , ρnm e inθ i ψnm e inθ = 1 2
k nm e i ϕnm .

(2.31)

Remark 2.1. The mode indices n and m are related to the number of oscillations in the azimuthal and radial directions respectively. Below, we will refer to n as the azimuthal mode index and m the radial mode index.

3 Numerical Computation of the Eigenvalues and Eigenfunctions

Numerical method

First, we discuss the special case n = 0. From (2.12), (2.13) (with λ = iν), the function ψ0 = i φ0 is a solution to:

r ∂ 2 ψ0 ∂r 2 -(α + 1) ∂ ψ0 ∂r + αν 2 c 1 c 2 r ψ0 = 0, (3.32)
with homogeneous boundary conditions ψ0 (R 1 ) = ψ0 (R 2 ) = 0. This a classical Bessel equation. Its solution is found e.g. in [START_REF] Bowman | Introduction to Bessel Functions[END_REF], p. 117 and is given by:

ψ0 (r) = r α+2 2 [AJ ñ(βr) + BY ñ(βr)] ,
for integer ñ; r α+2 2 [AJ ñ(βr) + BJ -ñ (βr)] , for noninteger ñ.

(3.33)

Here, J ñ and Y ñ are the Bessel functions of the first and second kinds respectively,

β = αν 2 c 1 c 2 , ñ2 = (α + 2) 2 4 ,
and (A, B) are determined from the boundary conditions. For the sake of simplicity, we focus on the case where ñ is not an integer, but the extension of the considerations below to integer ñ would be straightforward. The boundary conditions lead to a homogeneous linear system of two equations for (A, B). The existence of a non-trivial solution ψ0 requires that the determinant of this system vanishes. This leads to the following relation:

B(ν) = J ñ(βR 1 )J -ñ (βR 2 ) -J ñ(βR 2 )J -ñ (βR 1 ) = 0. (3.34)
By finding the zeros of B, we obtain the eigenvalues ν and then the corresponding eigenfunctions ψ0 .

For n ≥ 1, introduce the following numerical scheme. Given an integer N, we define a uniform meshsize h

= R 2 -R 1 N on the interval [R 1 , R 2 ] and discretization points R 1 = r 0 < r 1 2 < r 1 < • • • < r j < r j+ 1 2 < • • • < r N = R 2 , where r j = R 1 + jh and r j+ 1 2 = R 1 + j + 1 2 h. For each n, {ρ j+ 1 2
} and {ψ j } denote the numerical approximation of ρn and ψn = -i φn on grid points r j+ 1 2 's and r j 's respectively. The numerical scheme is

         c 1 n r j+ 1 2 ρ j+ 1 2 - c 1 r α+1 j+ 1 2 ψ j+1 -ψ j h = νρ j+ 1 2 , 0 ≤ j ≤ N -1, (3.35a) c 2 r α+1 j α ρ j+ 1 2 -ρ j-1 2 h + c 2 n r j ψ j = νψ j , 1 ≤ j ≤ N -1. (3.35b)
No boundary condition is imposed on ρ j . Concerning ψ j , we have ψ 0 = ψ N = 0.

Remark 3.1. We can modify the scheme (3.35) and use it to compute the solution in the case n = 0 by adding ρ j+ 1 2 and ψ j on both sides of the equations respectively.

Eigenvalules

In the numerical tests, we choose a set of parameter values given by

c 1 = 0.89307, c 2 = 0.69757, Θ = 0.2, R 1 = 1.9, R 2 = 2.1. (3.36)
Accuracy tests (not reported here) have demonstrated that the numerical scheme is of order 2 for any value of n. In the case n = 0, we can illustrate the good accuracy of the scheme by comparing the computed value of the eigenvalue to its analytic expression (3.34). The comparision is given in Table 1. With N = 1280 mesh points, the scheme mentioned in Remark 3.1 gives almost the exact eigenvalues. 3.2 lists the eigenvalues corresponding to the first seven radial modes m = 1, . . . , 6, for the first four azimuthal modes n = 1, . . . , 4 computed by the numerical scheme (3.35). It shows that, starting from m = 1, the radial eigenmodes come by conjugate pairs of almost the same absolute value but opposite signs (see Columns (1,2), [START_REF] Bertin | Boltzmann and hydrodynamic description for self-propelled particles[END_REF][START_REF] Bertin | Hydrodynamic equations for self-propelled particles: Microscopic derivation and stability analysis[END_REF] and (5,6) in 3.2). The fact that they are not exactly opposite can be attributed to the breaking of the clockwise-anticlockwise symmetry due to the linearization about a steady-state with definite orientation (here the clockwise rotating steady-state has been chosen). Also, the difference between c 1 and c 2 plays a role in this discrepancy (see Sec. 3.4) where it is shown that varying c 2 may increase it).

Eigenfunctions

For illustration purposes, we plot some of the eigenmodes ρ nm (r, θ), -ψ nm (r, θ) := ρnm (r) cos nθ, -ψnm (r) sin nθ .

(3.37)

For a better interpretation of the results, we plot the perturbation density ρ nm and the orientation vector

Ω nm (r, θ) = cos - π 2 -εψ nm (r, θ) , sin - π 2 -εψ nm (r, θ) . (3.38)
While the former corresponds to the perturbation only, the latter corresponds to the total solution (steady-state plus perturbation). We use a fairly large value of ε in order to magnify the influence of the perturbation. Since the chosen annular domain is rather thin, we rescale the plot onto an artificially wider annulus. Again, the chosen set of parameters is given by (3.36). Fig. 1 displays the modes (n, m) = (0, 4) (Figs. 1 (a,b)) and (n, m) = (4, 1) (Figs. 1 (c,d)). The left figures (Figs. 1 (a,c)) show the color-coded values of the density perturbations ρ nm (3.37) as functions of the two-dimensional coordinates (x, y) in the annulus. The right figures (Figs. 1 (b,d)) provide a representation of the orientation vector field Ω nm (3.38). In the case of mode (n, m) = (0, 4) (Figs. 1 (a,b)), since n = 0, the solution does not vary in the θ direction and the density perturbation ρ nm has two zeros in the r direction. In the case of mode (n, m) = (4, 1) (Figs. 1 (c,d)), the solution displays four periods in the θ direction and has only one zero of the density perturbation ρ nm in the r direction. a,b), since n = 0, there is no variation in θ and the solution is plotted at t = 1 to make the perturbation visible.

3.4 Variation of the parameters R 1 , R 2 , c 1 , c 2 and Θ

We numerically investigate the influence of the parameters R 1 , R 2 , c 1 , c 2 and Θ on the eigenvalues. We take the parameter values (3.36) as references. We vary one of the five parameters (c 1 , c 2 , Θ, R 1 , R 2 ) at a time, fixing the other values to those of (3.36). Fig. 2 (a) shows the eigenvalues ν as functions of the parameters R 1 . The insertedinside Fig. 2 (a) display how the eigenvalues depend on c 1 . Fig. 2 (b) shows how the eigenvalues depend on c 2 . Four eigenvalues corresponding to the modes n = 2, m = 0, 1, 2, 3 are displayed We observe that when the annular domain becomes narrower, i.e. R 1 is larger and closer to R 2 , the absolute value of ν is getting larger. The influence of R 2 (not displayed) is similar. As a result, the phase velocities of the modes become faster in a thinner domain, except for m = 0, which corresponds to no oscillation in the radial direction. As a function of c 1 and Θ, |ν| is monotonically increasing for all values of n (see insert inside Fig. 2 (a) for c 1 . The behavior as a function of Θ is similar and not displayed). The effect of a variation of c 2 is different: ν itself (instead of |ν|) is increasing with respect to c 2 .

Relaxation model in cylindrical coordinates

In this section, we discuss the numerical resolution of the nonlinear SOH model (1.1)-(1.3), subject to the boundary conditions (2.7). Its numerical solution will be compared with the solution of the linearized problem found in Sec. 3. We will further analyze how the nonlinear model departs from its linearization when the perturbation of the steady-state becomes large. One of the difficulties in solving the nonlinear model is the geometric constraint |Ω| = 1 (1.3) and the resulting non-conservativity of the model, arising from the presence of the projection operator P Ω ⊥ in (1.2). We rely on a method proposed in [START_REF] Motsch | Numerical simulations of a nonconvervative hyperbolic system with geometric constraints describing swarming behavior[END_REF] where the SOH model is approximated by a relaxation problem consisting of an unconstrained conservative hyperbolic system supplemented with a relaxation operator onto vector fields satisfying the constraint (1.3). In this section, we introduce this relaxation system in cylindrical coordinates in the annular domain.

The relaxation model is given by:

   ∂ t ρ η + c 1 ∇ • (ρ η Ω η ) = 0, (4.39a) ∂ t (ρ η Ω η ) + c 2 ∇ • (ρ η Ω η ⊗ Ω η ) + Θ∇ρ η = ρ η η 1 -|Ω η | 2 Ω η , (4.39b)
where η ≪ 1 and Ω η ∈ R 2 is not constrained to be of unit norm. The relaxation term at the right-hand side of (4.39b) contributes to making |Ω η | ≈ 1. In cylindrical coordinates, let Ω η = (q η cos φ η , q η sin φ η ), q η ≥ 0. Dropping the superindex η for simplicity, (4.39) can be written as

∂ t ρ + c 1 r ∂ ∂r (rρq cos φ) + ∂ ∂θ (ρq sin φ) = 0, (4.40) ∂ t (ρq cos φ) + c 2 r ∂ ∂r (rρq 2 cos 2 φ) + ∂ ∂θ (ρq 2 sin φ cos φ) -ρq 2 sin 2 φ (4.41) +Θ ∂ρ ∂r = ρ η (1 -q 2 )q cos φ, ∂ t (ρq sin φ) + c 2 r ∂ ∂r (rρq 2 sin φ cos φ) + ∂ ∂θ (ρq 2 sin 2 φ) + ρq 2 sin φ cos φ (4.42) +Θ 1 r ∂ρ ∂θ = ρ η (1 -q 2 )q sin φ.
Of course, we request that (ρ, q, φ) are 2π-periodic with respect to θ. We supplement the relaxation system with similar boundary conditions as (2.7). First, we request that the mass flux vanishes on ∂D, implying that

(ρq cos φ(r, θ, t))| r=R 1 ,R 2 = 0, ∀θ ∈ [0, 2π], ∀t ∈ R + .
When η ≪ 1, the relaxation term forces q ≈ 1. Therefore, we assume the same boundary condition (2.7) as for the SOH model, supplemented with the condition that q = 1, namely

φ(r, θ, t))| r=R 1 ,R 2 = ± π 2 , q(r, θ, t))| r=R 1 ,R 2 = 1, ∀θ ∈ [0, 2π], ∀t ∈ R + . (4.43)
We have the following theorem, whose proof is analogous to that of Proposition 3.1 in [START_REF] Motsch | Numerical simulations of a nonconvervative hyperbolic system with geometric constraints describing swarming behavior[END_REF] and is omitted. 4.2 Relaxation system in conservative form.

The scheme developed in [START_REF] Motsch | Numerical simulations of a nonconvervative hyperbolic system with geometric constraints describing swarming behavior[END_REF] relies on writing the hyperbolic part of the relaxation system in conservative form. Indeed, the use of a non-conservative form may lead to unphysical solutions, which are not valid approximations of the underlying particle system [START_REF] Motsch | Numerical simulations of a nonconvervative hyperbolic system with geometric constraints describing swarming behavior[END_REF]. Introducing (m, u, v) defined by m = rρ, u = rρq cos(φ + θ), v = rρq sin(φ + θ), Eqs. (4.40), (4.41) can be rewritten in terms of the vector function U = (m, u, v) as follows:

∂ t U + ∂ ∂r F (θ, U) + ∂ ∂θ G(r, θ, U) = 1 η H(U), (4.44)
where

H(U) =       0 u 1 - u 2 + v 2 m 2 v 1 - u 2 + v 2 m 2       , F (θ, U) =       c 1 (u cos θ + v sin θ) c 2 u m (u cos θ + v sin θ) + Θm cos θ c 2 v m (u cos θ + v sin θ) + Θm sin θ       , (4.45) and G
(r, θ, U) = 1 r       c 1 (v cos θ -u sin θ) c 2 u m (v cos θ -u sin θ) -Θm sin θ c 2 v m (v cos θ -u sin θ) + Θm cos θ       . (4.46)
Of course, we request that (m, u, v) is 2π periodic in θ. The boundary conditions (4.43) translate into:

(u cos θ + v sin θ)| r=R 1 ,R 2 = 0, (-u sin θ + v cos θ)| r=R 1 ,R 2 = ±1,

Numerical method

We apply the method proposed in [START_REF] Motsch | Numerical simulations of a nonconvervative hyperbolic system with geometric constraints describing swarming behavior[END_REF], which consists in splitting (4.44) into a conservative step and a relaxation step. In the conservative step, we solve (4.44) with H = 0. In the relaxation step, we solve (4.44) with F = G = 0. When η ≪ 1 this last step can be replaced by a mere normalization of Ω i.e. changing (u, v) into 5 Comparison between the Linear and Nonlinear Models

Small perturbation

We take a pure eigenmode as initial condition and compare the numerical solution of the nonlinear model to that of the linearized model. We take an initial condition given by

(ρ I , φ I ) = (ρ s , φ s ) + εk nm (ρ nm , -ψ nm ), (5.47)
with (ρ nm , -ψ nm ) given by (3.37). Let (ρ, φ) denote the exact solution of the nonlinear model (2.5), (2.6) with boundary conditions (2.7), (ρ ℓ , φ ℓ ) the solution of the linearized system given by Theorem 2.6, and (ρ h , φ h ) the numerical solution of the nonlinear model computed thanks to the method summarized at Sec. 4.3. Consider ρ for example. Formally, we have ρ -ρ ℓ = O(ε 2 ) (we neglect the errors due to the numerical computation of the functions ρmn which are small), while ρ -ρ h = O(h) (since the scheme is of order 1). Consequently, we have

ρ h -ρ ℓ = O(ε 2 ) + O(h).
(5.48) Fig. 3 (a) shows the L 1 -distance (below referred to as the "error") between the numerical solution of the nonlinear model and that of the linearized system at time t = 0.5, as a function of the meshsize h for an initial condition (5.47) corresponding to mode (n, m) = (3, 2) and k nm = 0.01. Different perturbation magnitudes ε = 0.001 (red squares), ε = 0.0005 (green triangles), ε = 0.0001 (blue crosses) are used. The parameter values are those of (3.36). We notice that for a given value of ε, the error decreases with decreasing values of h until h reaches the approximate values h = 0.01 (for ε = 0.001 and ε = 0.0005) and h = 0.005 (for ε = 0.0001). When h is decreased further, the error stays constant but this constant is smaller for smaller ε. This suggests that, consistently with (5.48), the error is dominated by the linearization error for small values of h. This interpretation is also consistent with the observation that the threshold value of h under which the error saturates decreases when ε becomes smaller. However, the decay of the error seems to be first order in ε instead of being second order as inferred from (5.48). This suggests that nonlinear effects are rapidly moving the solution away from the linear regime. However, other diagnostics discussed in the section below show that the linearized model actually provides a very good approximation of the nonlinear model in practical situations.

Large perturbations

In this section, we take larger values of ε and quantify the difference between the solutions of the nonlinear and linearized models. Due to nonlinear mode coupling, it is expected that, even with a pure mode initial condition, new modes will be gradually turned on by the nonlinearity. Let (ρ h,ε , φh,ε ) = ε -1 ((ρ h , φ h ) -(ρ s , φ s )) denote the difference between the numerical solution of the nonlinear model and the steady-state, rescaled by the factor ε -1 . We define the energy E(t) of the perturbation as

E(t) = ρh,ε φh,ε , ρh,ε φh,ε = 1 2 n≥1, m≥0 k 2 nm (t) + 1 2 m≥1 k 2 0 2m (t),
where the double bracket refers to the inner product (2.30) and k nm (t) is given by (2.31) with (ρ I , ϕ I ) replaced by (ρ h,ε , φh,ε ). The quantity k 2 nm (t)/2 (respectively k 2 0 2m (t)/2) represents the energy stored in the modes (±n, m) (respectively in the modes (0, 2m -1) and (0, 2m)) at time t. In the purely linear case, k nm (t) is independent of t. In the nonlinear case, its variation with t provides a measure of how the nonlinearity affects the amplitude of the corresponding modes.

The initial data is a perturbation of the steady-state by a pure eigenmode, i.e. ρh,ε φh,ε t=0 = k n 0 ,m 0 ρn 0 ,m 0 cos(n 0 θ) -ψn 0 ,m 0 sin(n 0 θ) , (5.49) with (n 0 , m 0 ) = (3, 2), k n 0 ,m 0 = 0.01. We test different values of ε. For this initial condition, Figs. 3 (b,c) show k nm (t) as a function of t in log-log scale for ε = 1 and 1.5 respectively. The initial mode (n 0 , m 0 ) = (3, 2) is represented with blue X's. In Fig. 3 (b) corresponding to a moderate perturbation ε = 1, only modes (n, m) = (0, 4) (red squares) and (n, m) = (6, 4) (purple triangles) appear. Mode (0, 4) appears first but saturates while mode (6, 4) appears later but reaches higher intensities. Both modes eventually saturate. Likewise, the initial mode decays as higher order modes (not represented in the figure) are turned on by the nonlinearity. The initial growth of modes (0, 4) and [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF][START_REF] Bertin | Hydrodynamic equations for self-propelled particles: Microscopic derivation and stability analysis[END_REF] is linear in log-log scale, which corresponds to a power law growth in time. The two modes have comparable growth rates (the two increasing parts of the curves are parallel straight lines). In the case of a larger perturbation ε = 1.5 displayed in Fig. 3 (c) the situation is strikingly more complex, with a wealth of other modes appearing. In addition to modes (n, m) = (0, 4) (red squares) and (6, 4) (purple triangle), we notice mode (6, 3) (cyan diamonds) and (3, 1) (green circles). Mode (6, 3) which was absent from Fig. 3 (b) now overtakes mode [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF][START_REF] Bertin | Hydrodynamic equations for self-propelled particles: Microscopic derivation and stability analysis[END_REF] at the beginning, but the latter reaches a higher intensity after some time. The decay of the initial mode [START_REF] Bertin | Boltzmann and hydrodynamic description for self-propelled particles[END_REF][START_REF] Bellomo | On the mathematical theory of the dynamics of swarms viewed as complex systems[END_REF] is also more pronounced. It should be noted that some modes stay extinct all the time. This shows that some pairs of modes are only weakly coupled by the nonlinearity.

In order to illustrate the successive turn on of the various modes, we have arbitrarily fixed a threshold value k t = 0.0005 (represented by the horizontal dashed blue lines on Figs. 3 (b,c)). In Fig. 3 (d), we have reported the first time t 1 at which k nm (t) reaches the values k t and plotted it as a function of ε in log-log scale, for modes (n, m) = (6, 4) (blue X's), (n, m) = (0, 4) (blue squares) and (n, m) = (6, 3) (red circles). The corresponding times t 1 are also indicated explicitly on Figs. 3 (b,c)). Fig. 3 (d) shows that for small ε, mode [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF][START_REF] Bertin | Hydrodynamic equations for self-propelled particles: Microscopic derivation and stability analysis[END_REF] is the earliest one to turn on. But as ε increases, this feature changes and mode (0, 4) (which was extinct for smaller value of ε) appears earlier. When ε is increased further, mode (6, 3) also appears, later than (0, 4) but earlier than [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF][START_REF] Bertin | Hydrodynamic equations for self-propelled particles: Microscopic derivation and stability analysis[END_REF]. This illustrates that the nonlinear mode coupling can exhibit rather complex features and nonmonotonic behavior as a function of the perturbation intensity ε. However, even for these large perturbation cases, the amplitude of the initial mode always remains one order of magnitude larger than those of the successively excited modes. This shows that the linear model still provides a fairly good approximation of the solution of the nonlinear model.

We now investigate the qualitative features of the solution in a large amplitude case. Figs. [START_REF] Bertin | Hydrodynamic equations for self-propelled particles: Microscopic derivation and stability analysis[END_REF] shows the numerical solution corresponding to a pure mode initial data (5.49) with (n 0 , m 0 ) = (4, 1), k n 0 ,m 0 = 1 and ε = 0.01. It displays the density ρ at times t = 0 Figure 3: (a) small perturbation case. L 1 -distance between the numerical solution of the nonlinear model and that of the linearized system at time T = 0.5, as a function of the meshsize h for an initial condition (5.47) corresponding to mode (n, m) = (3, 2): ε = 0.001 (red squares), ε = 0.0005 (green triangles), ε = 0.0001 (blue crosses). We use ∆t = 0.0001, N r = N θ = 10, 20, 40, 80, 160, 320 mesh points. The parameter values are those of (3.36). (b, c, d) large perturbation case. (b) k nm as a function of t in log-log scale for ε = 1 and for (n, m) = (3, 2) (the initial mode, blue X's), (0, 4) (red squares) and (6, 4) (purple triangles). (c) k nm as a function of t in log-log scale for ε = 1.5 and for (n, m) = (3, 2) (the initial mode, blue X's), (0, 4) (red squares), (6, 4) (purple triangles), (6, 3) (cyan diamonds) and (3, 1) (green circles). (d): Mode turn-on time t 1 as a function of ε for modes (n, m) = (6, 4) (blue X's), (0, 4) (blue squares) and (6, 3) (red circles). The mode turn-on time t 1 is the first time for which k nm (t) reaches the threshold value k t represented by the horizontal dashed blue line on Figs. (left) and t = 2 (right) as a function of the two-dimensional position coordinates (x, y) in the annulus, in color code (color bar to the right of the figure). We observe that the solution remains π/2-periodic in the θ-direction (as the linear mode would be) but the density contours have lost their sinusoidal shape. Instead, oblique shock waves have formed and are reflected by the boundary. These simulations suggest the existence of unsmooth periodic solutions of the nonlinear SOH model in this geometric configuration. We now investigate a large perturbation amplitude case with a random intitial data. More precisely, the initial data is given by a random combination of eigenmodes such that n ≤ 12 and m ≤ 12 as follows:

ρh,ε φh,ε t=0 = 0≤n,m≤12 (n,m) =(0,0) k nm ρnm cos(nθ + ϕ nm ) -ψnm sin(nθ + ϕ nm ) , (5.50)
where k nm and ϕ nm are randomly sampled in the intervals (0, 1] and [0, 2π] respectively, according to the uniform distribution. The numerical simulation is performed with N r = N θ = 640, ∆t = 0.0005 and ε = 0.0025. We observe a very good agreement between the linearized and nonlinear solutions, in spite of a fairly large perturbation amplitude. By looking carefully, one notices that the nonlinear solution has slightly lower maxima and larger minima, due to the action of numerical diffusion (which is absent from the linearized solution). The nonlinear solution also exhibits steeper gradients due to nonlinear shock formation. The use of the linearized solution results in considerable computational speed-up compared to that of the nonlinear one. Indeed, the computation of the eigenmodes and their summation to construct the solution is almost instantaneous on a standard laptop. By comparison, the computation of the nonlinear solutions takes of the order of an hour. Therefore, given the considerable computation speed-up, we consider that the performances of the linearized model are excellent. These performances make the linearized model a model of choice to perform parameter calibration on experimental data. Indeed, parameter calibration involves the iterative resolution of a minimization problem which consists of finding the set of parameters which minimize the distance between the solution and the data. With the linearized model, this calibration phase can be expected to require very little computational time. This is important, since this set of parameters is expected to change from one experiment to the next and consequently, the calibration phase must be performed for each experiment. A real-time analysis of an experiment therefore requires a very efficient algorithm. We observe a very good agreement between the linearized and nonlinear solutions, in spite of a fairly large perturbation amplitude.

Conclusion

In this paper, we have studied the SOH model on an annular domain. We have linearized the system about perfectly polarized steady-states. and shown that the resulting system has are only pure imaginary eigenvalues and that they form a countable set associated to an ortho-normal basis of eigenvectors. A numerical scheme for the fully nonlinear system has been proposed. Its results are consistent with the modal analysis for small perturbations of polarized steady-states. For large perturbations, nonlinear mode-coupling has been shown to result in the progressive turn-on of new modes in a complex fashion. Finally, we have assessed the efficiency of the modal decomposition to analyze the complex patterns of the solution. In future work, we will gradually include more physical effects in the model such as adding a repulsive force between the particles to prevent the formation of large concentrations, or immersing the particles in a surrounding fluid to give a better account of the dynamics of active particle suspensions like sperm. Finally, we plan to use the modal analysis to accurately calibrate the model against experimental observations of collective motion.

ρ

  λ (r, θ) = n∈Z ρn (r)e inθ , φ λ (r, θ) = n∈Z φn (r)e inθ , r, θ)e -inθ dθ.

  ρ nm : (n, m) = (0, 4), ν nm = 13.29. Ω nm : (n, m) = (0, 4), ν nm = 13.29. ρ nm : (n, m) = (4, 1), ν nm = -5.07. Ω nm : (n, m) = (4, 1), ν nm = -5.07.

Figure 1 :

 1 Figure 1: Density perturbation ρ nm (r, θ) (3.37) (Figs. (a, c)) and orientation vector Ω nm (3.38) (Figs. (b, d)) for mode (n, m) = (0, 4) (Figs. (a, b)) and mode (n, m) = (4, 1) (Figs. (b, d)), as functions of the two-dimensional cartesian coordinates (x, y) ∈ R 2 in the annulus (color online). The values of the density are color-coded according to the color-bar to the right of the figure. The orientation vector field is represented by blue arrows. The parameter set is given by (3.36) and N = 400 mesh points in the radial direction have been used. In Figs. (a, b), since n = 0, there is no variation in θ and the solution is plotted at t = 1 to make the perturbation visible.

  Varying R 1 (main figure) and c 1 (insert)

Figure 2 : 3 4

 23 Figure 2: Eigenvalues ν as functions of the parameters R 1 (Fig. (a)), c 1 (insert inside Fig. (a)) and c 2 (Fig. (b)). While varying one parameter, the other parameters are fixed to the values given by (3.36) and N = 400 mesh points in the radial direction have been used. We display four eigenvalues corresponding to the modes n = 2, m = 0, 1, 2, 3

Theorem 4 . 1 .

 41 The relaxation model (4.40)-(4.42) with boundary conditions (4.43) converges to the original model (2.5), (2.6) with boundary conditions (2.7) as η goes to 0.

  log(t 1 ) versus log(ε)

  Figure3: (a) small perturbation case. L 1 -distance between the numerical solution of the nonlinear model and that of the linearized system at time T = 0.5, as a function of the meshsize h for an initial condition (5.47) corresponding to mode (n, m) = (3, 2): ε = 0.001 (red squares), ε = 0.0005 (green triangles), ε = 0.0001 (blue crosses). We use ∆t = 0.0001, N r = N θ = 10, 20, 40, 80, 160, 320 mesh points. The parameter values are those of(3.36). (b, c, d) large perturbation case. (b) k nm as a function of t in log-log scale for ε = 1 and for (n, m) = (3, 2) (the initial mode, blue X's), (0, 4) (red squares) and (6, 4) (purple triangles). (c) k nm as a function of t in log-log scale for ε = 1.5 and for (n, m) = (3, 2) (the initial mode, blue X's), (0, 4) (red squares), (6, 4) (purple triangles), (6, 3) (cyan diamonds) and (3, 1) (green circles). (d): Mode turn-on time t 1 as a function of ε for modes (n, m) = (6, 4) (blue X's), (0, 4) (blue squares) and (6, 3) (red circles). The mode turn-on time t 1 is the first time for which k nm (t) reaches the threshold value k t represented by the horizontal dashed blue line on Figs. (b) and (c). The parameter values are those of (3.36) and ∆t = 0.0005, N r = N θ = 400.

= 2 Figure 4 :

 24 Figure4: Density ρ as a function of the two-dimensional position coordinates (x, y) in the annulus. The initial condition is given by (5.49) with (n 0 , m 0 ) = (4, 1), k n 0 ,m 0 = 1 and ε = 0.01. The solution is represented at time t = 0 (left) and t = 2 (right). The density is color coded according to the color bar to the right of the figure. The parameters are given by(3.36), ∆t = 0.0005, and N = 400 mesh points in both the radial and azimuthal directions have been used.

Fig. 5

 5 shows the numerical solution at time t = 2 (which approximately corresponds to the rotation of the fluid by a quater of a circle). It displays the density ρ as a function of the two-dimensional position coordinates (x, y) in the annulus, in color code (color bar to the right of the figure). Fig. 5 (a) shows the solution of the linearized model, obtained by summation of the corresponding eigenmodes, while Fig. 5 (b) displays the numerical solution of the nonlinear model with the same initial condition.

Figure 5 :

 5 Figure5: Density ρ as a function of the two-dimensional position coordinates (x, y) in the annulus. The initial condition is given by (5.50) with k nm and ϕ nm randomly chosen in the intervals (0, 1] and [0, 2π] respectively according to the uniform distribution, and ε = 0.0025. (a) Linearized solution. (b) Nonlinear solution. For the latter, the numerical simulation is performed with N r = N θ = 640 and ∆t = 0.0005. The solution is represented at time t = 2. The density is color coded according to the color bar to the right of the figure. We observe a very good agreement between the linearized and nonlinear solutions, in spite of a fairly large perturbation amplitude.

Table 1 :

 1 The eigenvalues ν for azimuthal mode n = 0 and various values of the radial mode index m.

	Comparison between the method using the Bessel functions (formula (3.34))
	and the scheme mentioned in Remark 3.1 with N = 1280 mesh points. The parameter
	values are given by (3.36).
	m	1	2	3	4	5	6
	ν(Bessel)	-6.6631 6.6631 -13.2895 13.2895 -19.9240 19.9240
	ν(Finite Difference) -6.6631 6.6631 -13.2895 13.2895 -19.9240 19.9240

Table 2 :

 2 The first seven radial modes m = 0, . . . , 6 for the first four azimuthal modes n = 1, 2, 3, 4, computed by the finite difference method (3.35) with N = 400 mesh points. The parameter values are given by(3.36). Starting from m = 1, the radial eigenmodes appear in conjugate pairs of almost (but not equal) absolute value and opposite signs (compare the pairs (m 1 , m 2 ) = (1, 2),[START_REF] Bertin | Boltzmann and hydrodynamic description for self-propelled particles[END_REF][START_REF] Bertin | Hydrodynamic equations for self-propelled particles: Microscopic derivation and stability analysis[END_REF],[START_REF] Bowman | Introduction to Bessel Functions[END_REF][START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]).

	m	0	1	2	3	4	5	6
	n
	1	0.4452 -6.2647 7.0618 -12.8913 13.6876 -19.5256 20.3217
	2	0.8905 -5.8668 7.4608 -12.4935 14.0860 -19.1277 20.7199
	3	1.3357 -5.4692 7.8603 -12.0958 14.4845 -18.7299 21.1182
	4	1.7810 -5.0720 8.2601 -11.6983 14.8833 -18.3323 21.5167

m √ u 2 +v 2 (u, v). The conservative step is solved by classical shock-capturing schemes (see[START_REF] Motsch | Numerical simulations of a nonconvervative hyperbolic system with geometric constraints describing swarming behavior[END_REF] for details). We take uniform meshes for r and θ. Careful accuracy tests (not reported here) have demonstrated that this method is of order 1.

This work was supported by the ANR contract 'MOTIMO' (ANR-11-MONU-009-01). The first author is on leave from CNRS, Institut de Mathématiques, Toulouse, France. He acknowledges support from the Royal Society and the Wolfson foundation through a Royal Society Wolfson Research Merit Award and by NSF Grant RNMS11-07444 (KI-Net). The second authors wishes to acknowledge the hospitality of the Department of Mathematics, Imperial College London, where this research was conducted. Both authors wish to thank F. Plouraboué (IMFT, Toulouse, France) for enlighting discussions.


R Artinyan
email: remy.artinyan@ec-lyon.fr

A Benamrouche

C Belacel

M Kozubova

A Berthelot

A M Jurdyc

G Beaudin

V Aimez

P Rojo-Romeo

J L Leclercq

R Peretti

X Letartre

S Callard

Optical characterization of active Photon Cages

Keywords: Near-infrared microcavity, photonic-crystal, air-mode, colloidal quantum dots

Recently, we developed a new family of 3D photonic hollow resonators which theoretically allow tight confinement of light in a fluid (gaz or liquid): the photon cages. These new resonators could be ideal for sensing applications since they not only localize the electromagnetic energy in a small mode volume but also enforce maximal overlap between this localized field and the environment (i.e. a potential volume of nano-particles). In this work, we will present numerical and experimental studies of the interaction of a photon cage optical mode with nano-emitters. For this, PbS quantum dot emitters in a PDMS host matrix have been introduced in photon cages designed to have optimal confinement properties when containing a PDMS-based active medium. Photoluminescence measurements have been performed and the presence of quantum dot emitters in the photon cages has been demonstrated.

INTRODUCTION

The general approach to achieve strong confinement of photons consists in high index contrast structuring of space at the wavelength scale. Numerous configurations have been used demonstrating the storage of light during a long time (high Q factor) into microcavities (small mode volume), leading to high Purcell effect. In particular, design and fabrication of multidimensional architectures for functional optical devices are of great interest for bio and health applications (such as sensing, fluorescence imaging) where it is aimed at strengthening the interactions with biological molecules or nanoparticles. For these sensing operations, hollow 3D resonators are ideal since they not only maximize the 3D light confinement but also enforce maximal overlap between the localized field and the environment (i.e. a potential volume of nano-particles). Recently, we proposed an original photonic concept to design new 3D photonic hollow resonators which theoretically allow tight confinement of light in vacuum or in a fluid (gas or liquid). The idea is to use a very thin and highly reflective membrane made of a 1D or 2D photonic crystal in semiconductor material as a non-absorbing broadband mirror. This effect, which has been demonstrated for 1D Fabry-Pérot cavities is the keystone of the resonator principle. In previous work, we focused on Si/air PC membranes. Simulations showed that reflectivity properties of these Photonic Crystal Membranes (PCM) are well preserved with curvature radius up to a few µm. We have therefore proposed to "bend" the photonic membrane in air, in order to achieve cylindrical optical cavities. Numerical simulations showed that for a 6µm diameter microcavity of 10µm-high silicon cylindrical pillars (diameter 0.27µm), strong confinement of light in air is obtained in the central hollow part of the cylinder (see Fig. 1) for relatively high Q modes (~1400 at 1560nm) [START_REF] Sieutat | 3D harnessing of light with photon cage[END_REF][START_REF] Sieutat | Strong confinement of light in low index materials[END_REF]. In this work, we propose to exploit these properties to enhance light-matter interaction by coupling an active material to the optical modes of the photon cage. For this we investigated the case of a photon cage filled with a layer containing PbS quantum dots. Such emitters can be inserted into a PDMS matrix deposited inside a photon cage. The photon cage properties, though, are strongly depending on the contrast index between the pillars and the background medium. We therefore designed new photon cages adjusted to the refractive index of the PDMS medium. In this paper, we will present the designing process of active photon cages, the activation (insertion of an active medium inside the cage) protocol and the first characterization results.

DESIGN OF ACTIVE PHOTON CAGES

To preserve the quality factor of the mode in the wavelength range around 1.5µm for a Si/low index medium cage, the geometrical parameters of the cage need to be optimized. In our case, the low index media is a PDMS resin in which variable concentration of PbS colloidal quantum dots have been inserted. The first step consists in designing the reflective membrane based on a 1D photonic crystal that allows the best reflecting properties; The period and filling factor of the PCM are then adjusted to optimize the properties of the 2D photon cage .The last step involves optimizing the photon cage diameter (or number of pillars) to obtain the best possible confinement with 10 m-high silicon pillars (light confinement in photon cages increases with the height of the pillars [START_REF] Sieutat | Strong confinement of light in low index materials[END_REF], and 10 m is the maximum height achievable with our reactive ion etching process while maintaining sufficient accuracy on the other parameters during fabrication). The photon cage designing process was conducted by 2D and 3D-FDTD simulations using the Lumerical and Harminv (by S.G.Johnson) [START_REF] Mandelshtam | Harmonic inversion of time signals[END_REF] softwares.

Design of a Si/low index medium 1D-PC mirror

In order to obtain mirrors with high reflectivity (>0.995) over a broad wavelength range (1500-1650 nm), the PCM parameters are optimized with a refractive index contrast between PDMS (n~1.5) and silicon (n=3.48). We simulated the PCM by defining one period of the photonic crystal, using periodic Bloch conditions at the boundaries. As the basic pattern of the photonic crystal, we chose ellipses instead of circles to be able to adjust its geometry more precisely and still be able to build the designed structures. We optimized the geometrical parameters of the PCM : the ellipse major and minor radiuses a and b, and the photonic crystal period p (see Fig. 2), and obtained a broadband mirror with a reflectivity close to 1 over a wavelength range of more than 150nm. Optimal reflectivity plateau was achieved for p=0.88 m, a=0.142 m, and b=0.161µm.

Two-dimensional photon cage design

The previously obtained PCM was then used as the basis for a 2D-cavity (see Fig. 3). Parameters a, b and p were first adjusted to obtain the best possible light confinement for a given 39-pillar photon cage. We obtained results quite close to the parameters of the optimal 1D-mirror : p=0.88 m, a=0.153 m, and b=0.159µm. It can be observed that the ellipses defining the pattern of the PCM present very small eccentricity. In the case of the 2D photon cage, the numbers of pillars (or the diameter of the cage) must also be optimized to obtain the best quality factor, leading to the best light confinement. We then proceeded to investigate the influence of the number of pillars (or photon cage diameter). We simulated the response of the photon cage to an excitation by a broadband dipole (center wavelength 1550nm, width 200nm) polarized along the axis of the cylindrical cavity and located at its center. The highest quality factor of all modes in a photon cage is overall increasing with the number of pillars of the cage, and it reaches ~21000 for a 76-pillar (or 21.3 m diameter) photon cage. The wavelength variation of the mode with the highest quality factor in a cage as a function of the number of pillars of the cage was also studied. It increases with the number of pillars up to 1613nm, and then drops to 1517nm and starts increasing again. These results are presented in Fig. 4. The constant increase of the mode wavelength and quality factor with the photon cage diameter is similar to the behavior of a Fabry-Pérot resonator. The sudden drop in the mode wavelength is explained by the optimal reflectivity plateau of the PCM. Indeed, as long as its wavelength stays on the PCM optimal reflectivity plateau, the mode with highest quality factor retains the same number of radial orders. Then, as the PCM is no longer efficient at the mode wavelength, the quality factor of another mode with more radial orders becomes higher than that of the previously best mode. The wavelength of this new mode increases again starting from the bottom end of the PCM optimal reflexion range (see Fig. 5).

Three-dimensional photon cage

The work realized for 2D photon cages only accounts for radial losses of the structures. We also simulated 10 m-high cages from the 2D-parameters obtained in the previous part (see Fig. 6). The 3D quality factor reaches a maximum for a 64-pillar cage (diameter 17.9 m ; Q>2600) To estimate the losses in the third dimension, we also plotted the difference between the losses in the 3D-cage (~1/Q 3D ) and the losses of the 2D-cage(~1/Q 2D ) . The difference should represent the losses in the vertical direction = (1/Q 3D )-(1/Q 2D ). For a mode with a given number of radial orders, the losses in the third dimension increase with the cavity diameter. Indeed, 2D losses become negligible compared to losses in the third dimension. A 64-pillar photon cage has therefore been chosen as a suitable candidate for fabrication and activation.

ACTIVATION OF A PHOTON CAGE

The challenge addressed in this part is to deposit a thin layer (10 m-high) of PDMS containing PbS quantum dots inside a photon cage. The PbS concentration in this layer needs to be homogeneous and controlled.

In order not to damage the photon cages and to obtain an homogeneous repartition of emitters inside and around the photon cages, we chose to deposit a drop of liquid solution containing the light emitters on photon cages. The PDMS is diluted in heptane, thus creating a fluid whose viscosity is small enough to be uniformly deposited on the sample containing photon cages. PbS quantum dots are then added to the solution so that a desired concentration of light emitters is reached. The result is then deposited onto the photon cages, and after drying and heating steps, the heptane is evaporated and we obtain a 10 m-high solid layer of homogeneous active medium on a big area of the sample (see Fig. 7). To prevent blue-shifting of the PbS dots caused by oxidation, we coat the active photon cage with a 10nm-thick layer of alumina (Al 2 O 3 ) deposited by pulsed laser deposition (PLD).

CHARACTERIZATION RESULTS

Photoluminescence experiments have been performed on samples like the one shown in Fig. 7. The photon cages used for the first characterization attempts, though, were composed of 39 cylindrical pillars. We used the experimental setup described in Fig. 8. We focus a laser beam at =980nm on the sample and measure the photoluminescence of the PbS quantum dots for laser spots inside and outside of photon cages. The results show an enhancement of measured photoluminescence inside photon cages (see Fig. 9). The measured emission spectrum, though, is not structured by coupling with photon cage modes.

Fig. 9 : Photoluminescence of PbS quantum dots inside photon cages (blue) and outside of photon cages (red). Despite photoluminescence enhancement, the measured spectrum of the light emitters seems to be unaffected by the cavities.

One possible explanation is that the observed photon cages are non-optimized 39-pillar photon cages, and the quality factors of our photon cage modes in the range delimited by the PbS emission spectrum may therefore be quite low. This underlines the importance of the photon cage designing process described above. Still, the PbS emission spectrum measured inside of photon cages proves that we successfully inserted an active medium in our microcavities.

CONCLUSION AND FUTURE PROSPECTS

We have designed new photonic micro-cavities adapted to the refractive index contrast between silicon and an homogeneous medium containing PbS quantum dots. A process has been established to reliably insert such a medium inside photon cages without damaging them. The first photoluminescence measurements performed on active photon cages show emission from the PbS quantum dots inside of photon cages. This suggests light-matter interaction between the photon cages and an active medium. This would be instrumental in using the photon cages as sensors, for example with the insertion of a micro-fluidic channel inside a photon cage. The photon cage modes, though, are yet to be pinpointed, probably due to their weak quality factors in the observed samples. Alternative techniques to characterize them include studying the losses from photon cage modes, producing active NSOM probes or using metal-coated multimode NSOM [START_REF] Grosjean | Polarization-dependant extraction properties of bare fiber probes[END_REF] probes on active photon cages. These ways are currently being explored in our team.

Fig. 1 :

 1 Fig. 1 : Left : SEM image of a 18-pillar, 6µm diameter photon cage. The pillars are 10µm high and their diameter is 0.27 µm; Right : Electrical field magnitude along a diameter of the cage (the dots are located on the boundaries of the photon cage) at =1560nm (mode quality factor Q=2500) .

Fig. 2 :

 2 Fig. 2 : Left : simulated PCM. Right : reflectivity obtained for p=0.88 m, a=0.142 m, and b=0.155 µm, b=0.161µm and b=0.175µm. The green curve shows good reflectivity on a broad wavelength band.

Fig. 3 :

 3 Fig. 3 : Cross-section of a 2D photon cage with 18 pillars, p=0.88 m, a=0.153 m, and b=0.159µm. The pillars are infinitely high.

Fig. 4 :

 4 Fig. 4 : Left : best mode quality factors of the optimized 2D-photon cage against the number of pillars. Right : wavelength of mode with the highest quality factor for cages of given numbers of pillars.

Fig. 5 :

 5 Fig.5: Top : Electrical field magnitude in the 18-pillar photon cage described in Fig.1. Insert : |E|-field mapping at =1560nm (Q~50000). Bottom : Electrical field magnitude in the 39-pillar photon cage with p=0.88 m, a=0.153 m, and b=0.159µm. Insert: |E|-field mapping at =1564nm (Q~3100).

Fig. 6 :

 6 Fig. 6 : Left : best mode quality factors of the optimized 3D-photon cage against the number of pillars. Right : 3D losses minus 2D losses for cages of given numbers of pillars.

Fig. 7 :

 7 Fig. 7 : Active photon cage. The final sample covers an area of several square millimeters containing multiple photon cages embedded in a PDMS layer containing PbS quantum dots.

Fig. 8 :

 8 Fig. 8 : Experimental photoluminescence characterization setup. The sample is similar to Fig. 7. As a light source, we used a Ti-Sa laser at =980nm. The laser spot on the sample has a diameter of 2 m.

ACKNOWLEDGEMENTS

We thank L.Lalouat and B.Gonzalez-Acevedo for their help and advice regarding FDTD computing and A.Pereira for performing the pulsed laser deposition of alumina on the active cage samples. IMUST is acknowledged for the funding of this research.