Oleg Zaitsev
email: oleg.zaitsev@uni-bonn.de

Gas lasers with wave-chaotic resonators

Keywords: numbers: 42, 55, Zz, 42, 55, Ah

Semiclassical multimode laser theory is extended to gas lasers with open twodimensional resonators of arbitrary shape. The Doppler frequency shift of the lineargain coefficient leads to an additional linear coupling between the modes, which, however, is shown to be negligible. The nonlinear laser equations simplify in the special case of wave-chaotic resonators. In the single-mode regime, the intensity of a chaotic laser, as a function of the mode frequency, displays a local minimum at the frequency of the atomic transition. The width of the minimum scales with the inhomogeneous linewidth, in contrast to the Lamb dip in uniaxial resonators whose width is given by the homogeneous linewidth.

Introduction

Theory of gas lasers was developed some decades ago and is well covered in the literature [START_REF] Sargent | Laser Physics[END_REF][START_REF] Haken | Laser Theory[END_REF]. The gas lasers are characterized by their capability to emit light at frequencies that differ from the atomic-transition frequency by more than a homogeneous linewidth. This property is a consequence of the coherent amplification by moving atoms that interact with the electric field at Doppler-shifted frequencies. If a resonance frequency of the cavity is continuously changed, the intensity of the single mode drops within the homogeneous linewidth of the atomic frequency. This narrow local minimum in the inhomogeneously broadened line, the Lamb dip [START_REF] Lamb | [END_REF], can be observed in uniaxial cavities that support standing waves. The Lamb dip has been utilized for stabilization of laser frequency [4] and precise detection of atomic-transition frequency and line shape (Lamb-dip spectroscopy) [5].

The theory of lasers (including gas lasers) was originally intended for low-loss resonators with a well-defined direction of propagation. The continuing interest in random [6,7] and chaotic [8] lasers has motivated extensions of the conventional theory to irregular and strongly open systems [8][9][10][11][12][13]. These improvements relate to a number of aspects, such as construction of a quasimode basis of the open cavity, description of linear and nonlinear coupling between the lasing modes in irregular lasers, and incorporating the statistics of eigenfrequencies and eigenfunctions in the laser equations (see Ref. [14] for a review).

In the present article, the semiclassical theory of gas lasers is generalized to lasers with open two-dimensional resonators of arbitrary shape. Equations for the amplitudes and frequencies of lasing modes are derived by expanding the atomic polarization in the powers of electric field and keeping the linear and third-order contributions. Since the dependence of the linear gain on the frequency is modified by the Doppler shifts in the frames of moving atoms, there exists an additional linear coupling between the modes. This effect is similar to the gain-induced coupling in systems with spatially inhomogeneous refractive index [11] or pumping [14]. We can argue, however, that the the magnitude of the coupling in gas lasers is negligible under reasonable assumptions. The nonlinear terms in the laser equations simplify for the weakly open cavities of chaotic geometry, where the eigenmodes can be modeled by a superposition of plane waves propagating in random directions [15,[START_REF] Berry | Chaotic Behaviour of Deterministic Systems[END_REF]. In this system, the intensity of a single lasing mode, as a function of its frequency, has a local minimum at the frequency of the atomic transition. In contrast to the Lamb dip in uniaxial resonators, the width of the minimum is comparable to the inhomogeneous linewidth.

Semiclassical theory of gas lasers

Semiclassical laser equations

We consider an open resonator of arbitrary shape, which, in general, can be filled with a homogeneous inactive medium characterized by a real dielectric constant .

The lasing action is provided by active atoms located inside the resonator. In the simplest model, only the the two atomic levels responsible for the lasing transition are described dynamically, whereas the information about the remaining degrees of freedom is contained in the empirical parameters: the homogeneous atomic linewidth γ ⊥ , the population-inversion decay rate γ , and the unsaturated population difference (per unit volume) ∆n 0 . Then the system of atoms can be characterized by the polarization P(r, t) and the population-inversion density (the difference between the populations of the upper and lower levels per volume) ∆n(r, t). The interaction of these quantities and the classical electric field E(r, t) in the resonator is specified by the semiclassical laser equations [START_REF] Haken | Laser Theory[END_REF]

c 2 ∂ 2 E ∂t 2 + ∇ × (∇ × E) = - 4π c 2 ∂ 2 ∂t 2 P(r, t), (1)
∂ 2 ∂t 2 + 2γ ⊥ ∂ ∂t + ν 2 P = -2ν d 2 E(r, t) ∆n(r, t), (2)
∂∆n ∂t -γ (∆n 0 -∆n) = 2 ν E(r, t) • ∂P(r, t) ∂t , (3)
where ν is the atomic-transition frequency and d is the atomic dipole matrix element.

The first equation follows from the Maxwell equations (written in the Gaussian units) and has the form of a wave equation with the polarization as a source. The second equation describes the dipole interaction with the field and the decay of polarization due to stochastic degrees of freedom. The third equation describes the energy exchange between the atoms and the field, as well as the relaxation of the population inversion.

If the field were absent, ∆n(r, t) would reach the stationary value of ∆n 0 . We will use this parameter as a measure of the pump strength. In the gas laser, the active atoms are moving with constant velocities v j , which are assumed to be Maxwell distributed. The polarization and the population-inversion density can be expressed in terms of the dipole moments p j (t) and the population differences ∆N j (t) of the individual atoms as

P(r, t) = N j=1 p j (t) δ r -r j (t) , (4) ∆n
(r, t) = N j=1 ∆N j (t) δ r -r j (t) , ( 5
)
where N is the number of atoms and r j (t) = r 0 j + v j t is the trajectory of the jth atom. The initial positions r 0 j are uniformly distributed within the resonator. The change of direction due to collisions or reflections at the boundary can be neglected if the time between these events is longer than the longest time scale in Eqs. ( 1)-(3), i.e., l, L v/γ , where l is the mean free path, L is the size of the system, and v is the root mean square velocity.

We will consider the case of a quasi-two-dimensional resonator, and assume that the field vector E(r, t) has a fixed polarization direction perpendicular to the plane of the resonator. This assumption is consistent with the Coulomb gauge condition ∇ • [ E(r, t)] = 0. Then the electric field and the dipole moments can be described in terms of their out-of-plane components, E(r, t) and p j (t).

Clearly, Eqs. ( 2) and (3) are local in the sense that the position r enters only as a parameter. Hence, the contributions of individual atoms can be separated by substituting the sums (4) and ( 5) and collecting the terms proportional to the respective delta functions:

∂ 2 ∂t 2 + 2γ ⊥ ∂ ∂t + ν 2 p j = -2ν d 2 E r j (t), t ∆N j (t), (6)
d∆N j dt -γ (∆N 0 -∆N j ) = 2 ν E r j (t), t dp j (t) dt , (7)
where ∆N 0 = ∆n 0 V /N is the unsaturated inversion per atom and V is the resonator volume.

Since the field in Eqs. ( 6) and ( 7) is evaluated at the positions of moving atoms, it is convenient to perform the Fourier transform

E(r, t) = 1 (2π) 2 dq E q (t) e iq•r , (8)
E q (t) = I dr E(r, t) e -iq•r , (9)
where the vectors r and q are two-dimensional. The transformation is applied only to the field inside the cavity, the integration being performed over the cavity interior I. In the moving frame of the atom, the Fourier component E q (t) is multiplied by an additional oscillating factor exp[iq • v j t], which expresses the Doppler shift for the frequency. In order to facilitate subsequent handling of lasing field oscillations, we carry out the temporal Fourier transform as well:

     E q (t) p j (t) ∆N j (t)      = 1 2π ∞ -∞ dω      E qω p jω ∆N jω      e -iωt , (10)
     E qω p jω ∆N jω      = ∞ -∞ dt      E q (t) p j (t) ∆N j (t)      e iωt . (11)
In the (q, ω) representation, the atom equations ( 6) and ( 7) assume the form:

-ω 2 -2iγ ⊥ ω + ν 2 p jω = -2ν d 2 1 (2π) 3 dq dω e iq•r 0 j E qω ∆N j,ω-ω +q•v j , (12)
(-iω + γ )∆N jω = 2πγ ∆N 0 δ(ω) - 2i ν 1 (2π) 3 dq dω ω e iq•r 0 j p jω E q,ω-ω +q•v j . (13)
Together with the wave equation (1), the above equations form a closed system, which, in principle, allows one to determine the electric field in the cavity.

Third-order perturbation expansion

In order to derive an equation for the electric field alone, the dipole moments p jω can be expressed in terms of the field using Eqs. ( 12) and ( 13) and then substituted in the polarization term of the wave equation [START_REF] Sargent | Laser Physics[END_REF]. To eliminate the inversion variable ∆N jω from the atom equations, we expand p jω and ∆N jω in powers of the field amplitudes E qω .

Here, the perturbation expansion is carried out up to the third-order terms, which provides the minimal description of the nonlinear effects and is a standard procedure in the laser theory [START_REF] Sargent | Laser Physics[END_REF][START_REF] Haken | Laser Theory[END_REF]. The third-order theory becomes invalid if the mode intensity is far above the lasing threshold, and a fully nonlinear approach along the lines of Refs. [12,13] has to be developed. The nonperturbative theories rely on the possibility to group terms according to their order of magnitude in the small parameter γ /δω, where δω is the frequency spacing between the lasing modes. These contributions can be traced back to oscillations of the population inversion, the zeroth order corresponding to the stationary-inversion approximation of Ref. [12]. The oscillations of the population at the beat frequencies in the multimode regime yield a correction to the polarization (linear in γ /δω and third order in the field), which oscillates at the same lasing frequency as the main contribution (zeroth order in γ /δω) to the polarization [START_REF] Sargent | Laser Physics[END_REF][START_REF] Zaitsev | [END_REF]. This is the level of approximation used in the present section. The leading-order corrections in γ /δω to the fully nonlinear theory of Ref. [12] were derived and numerically studied in Ref. [18]. A diagrammatic expansion in γ /δω, nonperturbative in the field, was proposed in Ref. [13]. In the frequency representation, the inversion pulsations are described by the function D (ω) defined below. In the laser with gaseous active medium, the Doppler frequency shifts mix the orders of magnitude of the pulsation terms, and additional procedures have to be introduced in order to approximately reclassify them. In Sec. 3.2, the stationary-inversion approximation is obtained for the nonlinear contribution of the third order in the electric field.

We start with the zeroth-order approximation to the population inversion by setting it equal to the unsaturated value: ∆N (0) jω = 2π∆N 0 δ(ω). Substitution in Eq. ( 12) yields the linear approximation to the dipole moment,

p (1) jω = -i d 2 γ ⊥ ∆N 0 D(ω) 1 (2π) 2 dq e iq •r 0 j E q ,ω+q •v j , (14)
where we introduced the notation

D(ω) = - 2iγ ⊥ ν ν 2 -ω 2 -2iγ ⊥ ω . ( 15
)
The real part of D(ω) provides the frequency dependence of the linear gain for atoms at rest. By substituting p

jω in Eq. ( 13), we obtain the quadratic correction to the inversion, ∆N

jω , which is substituted back into Eq. ( 12) to produce the nonlinear contribution,

p (3) jω = 2i νγ d 2 γ ⊥ 2 ∆N 0 D(ω) 1 (2π) 8 dq dq dq dω dω ω D(ω ) D (ω -ω + q • v j ) × e i(q +q +q )•r 0 j E q ,ω-ω -ω +(q +q )•v j E q ,ω +q •v j E q ω , (16)
where

D (ω) = 1 -i ω γ -1 . (17)
The resulting expansion for the dipole moment, p jω = p

(1)

jω + p (3)
jω , serves as a source in Eq. ( 1), which, after taking into account Eq. ( 4) and transformation to the ω representation, takes the form

ω 2 + c 2 ∇ 2 E ω (r) = -4π N j=1 1 (2π) 2 dq e iq•(r-r 0 j ) (ω -q • v j ) 2 p j,ω-q•v j . (18)
A possible means to solve this differential equation is to search for the solution E ω (r) in the form of an expansion in the eigenfunctions of the differential operator. In the case of an open cavity, one can either employ the continuous basis of eigenmodes of the infinite system or formulate a non-Hermitian eigenvalue problem by separating the functional spaces of the interior and exterior of the cavity. In the latter case, the complex eigenfrequencies provide the oscillation frequencies and lifetimes of the decaying quasimodes, whereas the right eigenfunctions yield their field distributions.

Here we make use of the quasimode basis, as the quasimodes supply, in most cases, the correct description of lasing modes near the lasing threshold. Specifically, we employ the constant-flux (CF) quasimodes [12], which are particularly suitable for laser problems (see Appendix A for more information). The CF modes are characterized by smooth wavefunctions ψ kω (r) depending on the real continuous parameter ω. In the asymptotic region outside the cavity, ψ kω (r) ∝ exp(iωr/c) describes a wave with frequency |ω| carrying constant energy flux. Within the cavity, the function satisfies the Helmholtz equation with the complex wavenumber Ω k (ω)/c. For a given ω, there exists an infinite set of eigenmodes labeled by the discrete (multi-)index k. The functions ψ kω (r) are the right eigenfunctions of a certain non-Hermitian operator. The conjugate left eigenfunctions, φ kω (r), behave as exp(-iωr/c) in the asymptotic region and have eigenfrequencies [Ω k (ω)] * inside the cavity. The CF eigenfunctions ψ kω (r) and φ kω (r) form complete biorthogonal bases of the open cavity.

Let us represent the electric field inside the cavity as a superposition of eigenfunctions ψ kω (r). After substitution of the expansion (A.6) in the left-hand side of Eq. ( 18), application of differential equation (A.2) and biorthogonality relation (A.4), we obtain an equation for the amplitude a kω of the mode k,

Ω 2 k (ω) -ω 2 a kω = 4π √ N j=1 1 (2π) 2 dq e -iq•r 0 j (ω -q • v j ) 2 φ * kqω p j,ω-q•v j , (19)
where φ kqω is the Fourier transform of φ kω (r), defined as in Eq. ( 9). The sum N j=1 can be substituted by the average N • • • over the positions r 0 j and velocities v j of active atoms. The distributions of the positions and velocities are uncorrelated: the former are distributed uniformly within the cavity and the latter obey the Maxwell distribution.

We approximate

e iq•r 0 j r 0 j 1 A I dr e iq•r = (2π) 2 A δ I (q), (20)
where (2π) 2 [δ I (q)] * is the Fourier transform of the characteristic function of the cavity interior and A is the cavity area. δ I (q) has a δ-function property dq δ I (q) f q = f q=0 , if f q is the Fourier transform of a function f (r) that vanishes outside of the cavity. By substituting the linear and third-order contributions to the dipole moment, Eqs. ( 14) and ( 16), we transform Eq. ( 19) to

Ω 2 k (ω) -ω 2 a kω = P (1) kω + P (3) kω , (21)
where the linear and third-order terms are, respectively, P

kω = -ic (1) 1 (2π) 2 dq φ * kqω E qω × (ω -q • v) 2 D(ω -q • v) v , (1)
kω = ic (3) 1 (2π) 8 dq dq dq dω dω ω D(ω ) φ * kqω E q-q -q ,ω × (ω -q • v) 2 D(ω -q • v) D [ω -ω -(q + q ) • v] E q ,ω-ω -ω -q •v E q ,ω +q •v v . (22) P (3)
In the velocity averages we dropped the index j labeling individual atoms. The coefficients are given by

c (1) = 4π ∆N 0 N √ A d 2 γ ⊥ , (24)
c (3) = 8π ∆N 0 N √ A 1 νγ d 2 γ ⊥ 2 . ( 25
)
The field amplitudes E qω in Eqs. ( 22) and ( 23) can be expanded in the CF quasimodes as in Eq. (A.6). Therefore, Eq. ( 21) represents a system of algebraic equations for the quasimode amplitudes a kω .

Lasing modes

The frequency integrals in Eq. ( 23) can be eliminated by approximating the electric field in the time domain with a superposition of constant-frequency oscillations, the lasing modes. It is well known that when the pumping increases above a certain threshold value, the first lasing mode appears. For higher pumping levels, additional modes start oscillating. In the multimode regime, the procedure of representing the field spectrum as a finite number of frequencies is an approximation: to preserve self-consistency of the laser equations, the oscillations at beat frequencies have to be neglected.

We start with an ansatz

E qω = 2π Nm l=1 E lq δ(ω -ω l ) + E * l,-q δ(ω + ω l ) , (26)
which assumes the existence of N m lasing modes with frequencies ω l > 0 and guarantees the realness of E(r, t). By expanding E lq in the CF quasimodes as

E lq = -1/2 k a lk ψ kqω l , (27)
we obtain the lasing-mode representation for the coefficients a kω ,

a kω = 2π Nm l=1 [a lk δ(ω -ω l ) + a * lk δ(ω + ω l )] . (28)
Usually, the quasimode with the longest lifetime starts oscillating at the lasing threshold. At higher pumping levels, the nonlinear terms in the laser equations mix the quasimodes, and the lasing mode no longer coincides with a single quasimode. Thus, the expansion ( 27) is nontrivial, as it contains, in general, more than one term.

After substituting the mode representations ( 26) and (28) in the laser equations ( 21) and collecting the terms proportional to δ(ω -ω l ) [or, equivalently, δ(ω + ω l )], we would obtain a rather cumbersome equation for the coefficients a lk . We write this equation in a simplified form by applying the standard approximation δω ν, where δω is the frequency spacing between the lasing modes. This approximation corresponds to neglecting the second time derivatives in Eqs. ( 1) and [START_REF] Haken | Laser Theory[END_REF]. Additionally, we neglect the terms containing D (2ν) ∼ γ /ν δω/ν. The resulting equation reads

(Ω k -ω l ) a lk = P (1) lk + P (3) lk , (29) with P
(1)

lk = - i 2 c (1) ν 1 (2π) 2 dq φ * kq D(ω l -q • v) v E lq , (30)
P (3) lk = i 2 c (3) ν 2 Nm l =1 1 (2π) 6 dq dq dq φ * kq E * l q × D(ω l -q • v) [D(ω l -q • v) + D * (ω l -q • v)] D [(q -q ) • v] v E l,q+q -q E l q + (1 -δ ll ) D(ω l -q • v) [D(ω l -q • v) + D * (ω l -q • v)] × D [ω l -ω l + (q -q ) • v] v E l ,q+q -q E lq . (31)
In the equations above we introduced the notation Ω k ≡ Ω k (ν), φ kq ≡ φ kqν , and

D(ω) = 1 -i ω -ν γ ⊥ -1 , (32)
which approximates D(ω) for positive frequencies. The field amplitudes E lq are assumed to be expanded in the CF wavefunctions ψ kq ≡ ψ kqν , as in Eq. ( 27). The nonlinear laser equations (29) can be solved for the frequencies and amplitudes of the lasing modes, e.g., by using the iteration procedure described in Ref. [19].

The term proportional to (1 -δ ll ) in Eq. ( 31) is usually omitted in the case of fixed atoms, v = 0. In such a system, the factor D (ω l -ω l ) ∼ γ /δω ∼ γ /γ ⊥ 1, arising due to oscillations of the population inversion, makes this term small compared to the first term in the braces. In Sec. 3.2 we argue that for a gas laser the second term is negligible as well.

Although the algebraic equation (29) has been derived above for the steady-state lasing oscillations, a lk exp(-iω l t), it can be easily generalized [14] to the case of slowlyvarying amplitudes a lk (t) by adding the term -i ȧlk to the left-hand side. The resulting first-order nonlinear differential equation can be used to study stability of the lasing modes.

Properties of the laser equations

In this section, we simplify the right-hand side of the laser equations ( 29) by taking into account the resonant dependence of the quasimode wavefunctions on q. In particular, we show that the linear gain-induced coupling of quasimodes (Sec. 3.1) and the nonlinear contribution proportional to (1 -δ ll ) (Sec. 3.2) are negligible.

Linear coupling of quasimodes

After the field amplitude E lq is expanded in the basis functions ψ kq , the linear contribution (30) can be written as

P (1) lk = - iν 2 √ c (1) k V kk (ω l ) a lk . (33)
The matrix elements

V kk (ω l ) = 1 (2π) 2 dq φ * kq D(ω l -q • v) v ψ k q (34)
describe possible linear coupling of the quasimodes via the Doppler shift in the frequency dependence of the linear gain. A similar coupling in the coordinate space arises in the systems with inhomogeneous refractive index or pump power distribution [11,14]. Below we argue, however, that the off-diagonal elements V kk (ω l ) are negligible, i.e., the quasimode wavefunctions correctly describe the lasing modes in the linear approximation. First, we evaluate the velocity average in Eq. (34). In the Doppler limit

|ω l -ν| ∼ νv/c γ ⊥ , (35)
typical for the gas lasers, the function (32) can be approximated as

D(ω) γ ⊥ πδ(∆ω) + iP 1 ∆ω , ∆ω ≡ ω -ν, (36)
where P denotes the principal value. Then the angular average is given by

D(ω -qv cos ϕ) ϕ = γ ⊥ θ (qv -|∆ω|) q 2 v 2 -∆ω 2 + i θ (|∆ω| -qv) ∆ω 2 -q 2 v 2 , ( 37
)
where θ(x) is the Heaviside step function. After the averaging over the two-dimensional Maxwell distribution,

W (v) dv = 1 πv 2 e -v 2 /v 2 v dv, (38)
we find

D(ω -q • v) v = γ ⊥ qv 1 2 √ π e -(∆ω/qv) 2 + i π F ∆ω qv , (39)
where F (x) ≡ exp(-x 2 )

x 0 exp(y 2 ) dy is the Dawson integral [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. As a function of q, D(ωq • v) v changes on the scale ∆q ∼ |∆ω|/v ∼ ν/c, for typical ∆ω.

Next, we note that ψ kq and φ kq are significant only when the magnitude of q lies within the band |q -q k | δq k , where q k = √ Re Ω k /c and δq k = -√ Im Ω k /c. Therefore, the matrix elements V kk (ω l ) vanish if |q k -q k | δq k . In the reverse case, |q k -q k | δq k , we can pull D(ω lq • v) v out of the integral (34), since this function is constant on the scale δq k ∆q. Then Eq. (A.5) yields the diagonality property,

V kk (ω l ) = D(ω l -q • v) v,q=q 0 δ kk , (40)
where the substitution of q = q 0 ≡ √ ν/c is consistent with the level of approximation in Eq. (29). We note that the linear-gain frequency dependence, given by the real part of Eq. (39), has a Gaussian shape, similar to one-dimensional gas lasers [see Eq. (B.3)].

Nonlinear terms

The resonant dependence of the wavefunctions ψ kq and φ kq on q, mentioned previously, makes it possible to simplify the nonlinear terms of Eq. (29) as well. The condition |q -q k | δq k implies one of the following restrictions on the integration domain in Eq. ( 31)

: (i) |q -q | δq k , (ii) |q -q | δq k , or (iii) |q + q | δq k .
Only in the first case we can estimate D [(qq ) • v] ∼ min(1, γ /δq k v). In the cases (ii) and (iii), we have D

[(q -q ) • v] ∼ γ c/νv min(1, γ /δq k v). Similarly, D [ω l -ω l + (q -q ) • v]
min(1, γ /δq k v) in all three cases. Thus, we can neglect the terms in Eq. ( 31) that are proportional to (1 -δ ll ) and rewrite this equation with the new integration variables Q = (q + q )/2 and δq = qq corresponding to the choice (i):

P (3) lk =iν 2 c (3) Nm l =1 1 (2π) 2 2π 0 dϕ 2π 0 dϕ Y (ϕ, ϕ , v) × D(ω l -q 0 v cos ϕ) L(ω l -q 0 v cos ϕ ) v , (41)
where L(ω) ≡ ReD(ω), and q = Q = q 0 is set in all nonresonant functions. The value of the correlator

Y (ϕ, ϕ , v) = q 2 0 (2π) 4 ∞ 0 dq ∞ 0 dQ d(δq) D (δq • v) × φ * kq E l,q+δq E * l ,Q+δq/2 E l ,Q-δq/2 (42)
is system specific. Here ϕ (ϕ ) is the azimuthal angle of the vector q (Q) counted from the direction of v. We can approximate D (δq • v) ≈ 1 under the assumption of long living quasimodes, such that δq k γ /v. In this case, the correlator depends only on the direction of v, but not on its magnitude. The absence of functions D from the laser equations is an indication that the approximation of constant population inversion is employed.

Chaotic gas laser

Nonlinear terms in the laser equations

We apply the general theory to a gas laser with wave-chaotic resonator. The notion of a chaotic system implies that the shape of the resonator does not permit a separation of variables in the wave equation and the relevant modes have a typical wavelength λ = 2π/q 0 L. Famous examples of chaotic resonators are the Sinai and stadium billiards and their desymmetrized versions [START_REF] Reichl | The Transition to Chaos: Conservative Classical Systems and Quantum Manifestations[END_REF]. The classical ray dynamics, defined in the shortwavelength limit, is assumed to be ergodic, i.e., the generic ray trajectories fill densely the available phase space. This assumption sets a limitation on the degree of openness of the resonator: the lifetimes of the CF modes contributing to the lasing modes are to be longer than the time needed for the ray trajectories to establish ergodicity. This time is of the order of inverse Lyapunov exponent. We will consider a weakly open resonator, such that its eigenmodes are well approximated by real eigenfunctions ψ k (r) ≈ φ k (r) independent of the parameter ω.

Let us examine a product of two wavefunctions evaluated at nearby points in the wavevector space:

ψ * k,Q+δq/2 ψ k ,Q-δq/2 = d(δr) e -iQ•δr I dR e iδq•R × ψ k R - δr 2 ψ k R + δr 2 . ( 43
)
Since the relevant magnitude of δq is much smaller than the optical wavenumber q 0 , the function e iδq•R changes on the length scale much longer than λ. On the other hand, the wavefunction product self-averages in the chaotic cavity over a correlation area of order λ 2 [15]. Consequently, the product can be substituted with its average value, expressed in terms of the correlation function

C kk (δr) = A ψ k R - δr 2 ψ k R + δr 2 R , (44)
where the wavefunctions are assumed to be normalized. Different eigenmodes, k = k , are uncorrelated in the chaotic system, whereas the autocorrelation function in two dimensions is approximated by a Bessel function [15],

C kk (δr) = J 0 (q k δr). (45)
This expression shows that the statistical characteristics of the eigenstates are isotropic (this approximation fails within a wavelength of the cavity boundary). Scars [START_REF] Kaplan | [END_REF], i.e., the excess probability density around unstable short periodic ray trajectories, should not significantly influence the correlation function. Indeed, the excess probability density depends on the inverse Lyapunov exponent, which scales with the system size. On the other hand, the fraction of the phase space occupied by a scar reduces with the wavelength, which, in gas lasers, is several orders of magnitude smaller than the system size. The product (43) is proportional to the Fourier transform of the autocorrelation function, (2π/q k ) δ(Q -q k ), which makes

ψ * k,Q+δq/2 ψ k ,Q-δq/2 = (2π) 3 q k A δ kk δ(Q -q k ) δ I (δq). (46)
We use this result to calculate the third-order contribution (41) to the laser equations. An important simplification that comes with the stochasticity of the wavefunctions is the isotropy of the correlator (42). Hence, this function can be pulled out of the angular integrals, and we find

P (3) lk = iν 2 3/2 A c (3) a lk Nm l =1 k |a l k | 2 × D(ω l -q 0 v cos ϕ) ϕ L(ω l -q 0 v cos ϕ) ϕ v . ( 47
)

Single-mode regime

We apply the laser equations (29), with the right-hand side given by Eqs. (33), (40), and (47), to determine the dependence of the mode intensity on the frequency in the single-mode regime. If the pumping level is not far above the threshold, only one quasimode will contribute to the lasing mode. Consequently, a single term with l = l and k = k remains in the summations in Eq. (47). Taking the imaginary part of Eq. ( 29), we obtain an expression for the intensity,

I k = |a lk | 2 / , I k (ω) I 0 = L(ν -q 0 v cos ϕ) 2 ϕ v L(ω -q 0 v cos ϕ) 2 ϕ v e - " ∆ω q 0 v " 2 - ∆N 0,k ∆N 0 , (48)
where the intensity unit

I 0 = A 8 √ π 2 γ ⊥ γ d 2 q 0 v L(ν -q 0 v cos ϕ) 2 ϕ -1 v (49)
is defined by the condition I k (ν) → I 0 , as ∆N 0 → ∞. The threshold value of the pumping strength ∆N 0 at ω = ν is given by the parameter

∆N 0,k = κ k A √ πN q 0 v d 2 ν , (50)
where

κ k = -Im Ω k is the decay rate of the quasimode. When evaluating L 2 ϕ v
, it is not possible to use the approximation (37) for the angular average, because the velocity integral would diverge. The angular integral can be calculated analytically for a finite value of γ ⊥ /q 0 v, however, the resulting expression is rather cumbersome. The remaining integral over the distribution W (v) has to be evaluated numerically. In the special case ω = ν the result can be expressed in terms of a special function:

L(ν -q 0 v cos ϕ) 2 ϕ v = x 2π e x E 1 (x), (51)
x ≡ γ ⊥ q 0 v 2 , (52)
where E 1 (x) =

∞

x dy e -y /y [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. The relative intensity I k (ω)/I 0 is plotted in Fig. 1 for several values of γ ⊥ /q 0 v in the Doppler limit (35). The frequency dependence is characterized by a local minimum of intensity at the center of the gain curve, ω = ν. The frequency interval where the intensity is reduced has a width of the order of q 0 v. Formally, this follows from the analysis of the denominator of Eq. (48). In fact, the integral

L(ω -q 0 v cos ϕ) ϕ = 1 π π 0 dϕ 1 + (∆ω -q 0 v cos ϕ) 2 /γ 2 ⊥ (53)
receives its major contribution from the angles in the interval δϕ ∼ γ ⊥ /q 0 v around ϕ = cos -1 (∆ω/q 0 v). Therefore, the integral changes its order of magnitude from γ ⊥ /q 0 v, for |∆ω| < q 0 v -γ ⊥ , to (γ ⊥ /∆ω) 2 , for |∆ω| > q 0 v + γ ⊥ , i.e., q 0 v is the characteristic frequency scale.

The results for the chaotic laser are compared in Fig. 1 with the relative mode intensity in a one-dimensional gas laser (Appendix B). In the former (latter) case, ∆ω is 2 (see the text). Parameters: relative pumping strength ∆N 0 /∆N 0,k = 1.3; relative homogeneous linewidth γ ⊥ /q 0 v = 0.01 (solid curves), 0.04 (dashed curves), and 0.07 (dash-dotted curves).

normalized by q 0 v ( √ 2q 0 v). Apart from the visual convenience, the difference in scaling is motivated by the relation, v 2D = √ 2 v 1D , between the root mean square velocities in two and one dimensions at a given temperature. Lasers with uniaxial resonators exhibit a much narrower intensity minimum, known as Lamb dip [START_REF] Sargent | Laser Physics[END_REF]. Its width is of the order of the homogeneous linewidth γ ⊥ q 0 v. In order to understand the origin of the dip, it is convenient to treat the lasing mode as a superposition of left-and rightmoving waves. The atoms moving with the velocity v x ≈ ∆ω/q x effectively interact with only one of the two waves, depending on the sign of ∆ω and v x . However, when |∆ω| γ ⊥ , the atoms at rest interact with both running waves. This leads to a stronger gain saturation near ∆ω = 0 than away from the line center. In the two-dimensional chaotic laser this effect is washed out, since a mode with the frequency ω draws its gain from a larger ensemble of atoms with velocities v > |∆ω|/q 0 . In particular, let us model the chaotic wavefunction locally by a random superposition of partial plane waves with wavevectors of fixed magnitude √ ω/c and isotropically distributed directions [START_REF] Berry | Chaotic Behaviour of Deterministic Systems[END_REF]. Then, all partial waves of the central mode ω = ν shall interact with the atoms at rest, but each partial wave interacts, in addition, with atoms moving perpendicular to its direction of propagation. Thus, the minimum of intensity results from a more subtle interplay of linear amplification and nonlinear saturation on the same frequency scale.

The results obtained above for a wave-chaotic resonator cannot be directly applied to resonators with mixed or regular classical ray dynamics, where the hypothesis of random, isotropically distributed, and uncorrelated eigenstates is no longer valid. Nevertheless, one can conjecture that increasing dimensionality of the system will make the intensity minimum shallower and wider (its depth and width being defined relatively to the intensity maximum and the inhomogeneous linewidth, respectively), but will not lead to its complete disappearance. As a practical consequence, the frequency stabilization [4] and spectroscopy [5] might still be possible in imperfect quasi-onedimensional resonators, where the lasing mode deviates from a superposition of two counter-propagating plane waves. It is worth mentioning that in both uniaxial and chaotic resonators the relative depth of the minimum decreases with increasing γ ⊥ .

Conclusions

In the present article, the semiclassical laser theory was extended to gas lasers with two-dimensional resonators of arbitrary shape. It was shown that the linear coupling between the modes, resulting from the Doppler shift in the gain frequency dependence, is negligible. The nonlinear coupling was considered in the third order of the perturbation theory in the electric field. The nonlinear terms in the laser equations that arise due to pulsations of the population inversion were identified and neglected. The criterion (ν/c) v γ , allowing one to approximate the population inversion as constant, is usually well fulfilled (ν is the atomic-transition frequency and v is the average velocity).

The general theory was applied to lasers with two-dimensional weakly open resonators of wave-chaotic geometry. In the short-wavelength limit, the eigenfunctions of the resonator can be locally approximated by a superposition of plane waves with a fixed wavelength propagating in random directions. The applicability condition for this description can be estimated as c/ν w L, where L is the size of the resonator and w is the width of the output window. The isotropic statistical characteristics of the eigenfunctions lead to decoupling of the field distribution from the Doppler-broadened gain curve in the nonlinear terms of the laser equations. The single-mode intensity, as a function of the laser frequency, has a local minimum at the frequency of the atomic transition. The width of the minimum scales approximately with the inhomogeneous broadening (ν/c) v and only weakly depends on the homogeneous linewidth γ ⊥ . This property distinguishes the intensity minimum in chaotic resonators from the Lamb dip in uniaxial resonators, which has the width of γ ⊥ (ν/c) v.

distribution with the outgoing flux and no incoming flux, the sign is chosen according to the requirements:

Re[Ω k (ω)] > 0, Im[Ω k (ω)] < 0 (ω > 0), (A.9) Ω k (-ω) = -[Ω k (ω)] * . (A.10)
In connection with the definition (A.8), this choice provides for negative Im[Ω k (ω)] for any ω, as is expected from the decaying quasimodes.

Appendix B. One-dimensional gas laser

We apply the general theory to a one-dimensional gas laser that was studied in detail in Ref. [START_REF] Sargent | Laser Physics[END_REF]. Specifically, we consider a resonator defined by = 1 in the region 0 < x < L, with a perfect mirror at x = 0 and a weakly transmitting mirror at x = L. The quasimode eigenfunctions have the form

ψ k (x) = e iΩ k x/c i √ 2L - e -iΩ k x/c i √ 2L ≡ ψ (+) k (x) + ψ (-) k (x), (B.1)
where we approximate q k ≡ Re Ω k /c ≈ πk/L (k = 1, 2, . . .) and assume δq k ≡ -Im Ω k /c q k . In principle, δq k can be determined by modeling the semitransparent mirror as a dielectric layer with large dielectric constant and solving Eqs. (A.1) and (A.2), but the specific value of δq k is not important for our purposes. In the q representation, the functions

ψ (±) kq = 1 √ 2L
1 -e -i(q∓q k ±i δq k )L q k -i δq k ∓ q (B.2) have a resonant dependence on q, peaked at ±q k and having a width of δq k . As in the two-dimensional case, the velocity average of the homogeneous gain curve can be computed within the approximation (36). Averaging over the one-dimensional Maxwell distribution yields

L(ω -qv) v = π 2 γ ⊥ qv exp - 1 2 ∆ω qv 2 . (B.
3)

The third-order contribution is obtained from Eq. (31), adapted to the one-dimensional case. The electric field can be represented as a superposition of the quasimode wavefunctions (B.2). The velocity averages depend on q on a scale larger than δq k , and can be pulled out of the integrals and evaluated at the resonant values of q. The integrals over the wavefunctions have a quite simple form in the coordinate representation: 1 (2π) 3 dq dq dq ψ In the single-mode near-threshold regime, when only one quasimode is excited, six of the 16 possible combinations of indices are resonant. However, only four combinations, (α, α , α , α ) = (+ -++), (-+ --), (+ ---), (-+ ++), yield q = q at resonance, which results in D = 1. In the limit δq k L 1, the resonant value of the integral (B.4) is equal to 1/4L. The velocity averages in Eq. ( 31) with the resonant values of q = ±q k and q = q ± q k sum up to 4 D(ω -q k v) [L(ω -q k v) + L(ω + q k v)] v . This quantity was calculated in Ref. [START_REF] Sargent | Laser Physics[END_REF]. The real part in the Doppler limit is given by

L(ω -qv) [L(ω -qv) + L(ω + qv)] v = π 8 γ ⊥ qv exp - 1 2 ∆ω qv 2 [1 + L(ω)] . (B.6)
The imaginary part of Eq. ( 29) yields an expression for the mode intensity I k = |a lk | 2 :

I k (ω) I 0 = 2 1 - ∆N 0,k ∆N 0 e 1 2 " ∆ω q 0 v " 2 1 + L(ω) , (B.7)
where we approximated q k ≈ q 0 . The parameters

I 0 = 2 γ ⊥ γ L

Figure 1 .

 1 Figure 1. Relative intensity of the single lasing mode for the two-dimensional (2D) chaotic laser [Eq. (48), thick curves] and a one-dimensional (1D) laser [Eq. (B.7), thin curves] as a function of the frequency separation ∆ω from the atomic-line center. The scaling of ∆ω for these systems differs by a factor of √2 (see the text). Parameters: relative pumping strength ∆N 0 /∆N 0,k = 1.3; relative homogeneous linewidth γ ⊥ /q 0 v = 0.01 (solid curves), 0.04 (dashed curves), and 0.07 (dash-dotted curves).

  α , α , α = +,and we took into account that φ * k (x) = ψ k (x). The phase of the integrand vanishes when the indices obey the resonance condition αk + α k -α k + α k = 0.(B.5)

  Sec. 4.2.

Acknowledgments

Financial support was provided by the Deutsche Forschungsgemeinschaft via the grant FOR557.

Appendix A. Constant-flux quasimodes

We consider an open cavity with a real dielectric constant (r) = const inside and = 1 outside of the system's boundary. The constant-flux (CF) modes [12] ψ kω (r), depending on the parameter ω > 0, are defined as solutions of a non-Hermitian eigenvalue problem.

Explicitly, the functions must satisfy the differential equation

in the exterior of the cavity with the outgoing-wave boundary conditions at infinity, while inside the system the same mode satisfies a different equation:

For any given value of ω, the complex eigenfrequency Ω k (ω) is quantized [k is a discrete (multi-)index labeling the modes], because the solutions are required to match smoothly at the interface.

The conjugate wavefunctions φ kω (r) obey Eq. (A.1) with the incoming-wave boundary conditions outside the cavity and the equation

inside the system. The CF quasimodes and their conjugates are biorthogonal, and can be chosen to satisfy the condition

where the integration is over the interior I. A similar relation,

is valid in the q representation, defined as in Eq. ( 9). Additionally, the wavefunctions possess the properties: φ kω (r) = [ψ kω (r)] * and φ kqω = ψ * k,-q,ω . A Fourier component E ω (r) of the lasing field can be expanded in the CF modes as

When continued to the exterior, this expansion yields a wave at the frequency ω propagating in the free space away from the system. In the stationary lasing regime ω becomes restricted to a finite number of values corresponding to the frequencies of lasing modes. To guarantee the realness of the field E(r, t), we need to extend the definition of the CF quasimodes to negative frequencies by

Since the eigenfunctions φ kω (r) depend only on the squared complex frequency Ω 2 k (ω), the frequencies Ω k (ω) are determined up to an overall sign. To describe the field


C Fiorio
email: fiorio@lirmm.fr

D Jamet
email: jamet@lirmm.fr

J.-L Toutant
email: toutant@lirmm.fr.

Discrete circles: an arithmetical approach with non-constant thickness

Keywords: arithmetical discrete geometry, discrete circle, discrete curve, Bresenham, norms

In the present paper, we introduce an arithmetical definition of discrete circles with a non-constant thickness and we exhibit different classes of them depending on the arithmetical discrete lines. On the one hand, it results in the characterization of regular discrete circles with integer parameters as well as J. Bresenham's circles. As far as we know, it is the first arithmetical definition of the latter one. On the other hand, we introduce new discrete circles, actually the thinnest ones for the usual discrete connectedness relations.

INTRODUCTION

Discrete geometry attempts to provide an analogue of euclidean geometry for the discrete space Z n . Such an investigation is not only theoretical, but has also practical applications since digital images can be seen as arrays of pixels, in other words, subsets of Z 2 .

Since the seventies, the discrete lines, namely the analogue of euclidean lines in the discrete space Z n , have been widely studied. On the one hand, J. Bresenham, [START_REF] Bresenham | Algorithm for computer control of a digital plotter[END_REF] A. Rosenfeld 15 and H. Freeman 7 followed an algorithmical approach and defined discrete lines as digitizations of euclidean ones. They provide tools for drawing and recognition. [START_REF] Freeman | Computer processing of line-drawing images[END_REF][START_REF] Hung | On the straightness of digital arcs[END_REF]15 On the other hand, J.-P. Reveillès initiated the discrete arithmetical geometry 14 and introduced the arithmetical discrete lines as subsets of Z 2 satisfying a double diophantine inequality (see Definition 3). Such an approach enhances the knowledge of the discrete lines. In addition to give new drawing 14 and recognition [START_REF] Debled-Rennesson | A linear algorithm for segmentation of digital curves[END_REF] algorithms, it directly links topological and geometrical properties of an arithmetical discrete line with its definition. For instance, the connectedness of an arithmetical discrete line is entirely characterized (see Theorem 1).

The other well-studied family of discrete objects is the one of discrete circles. Similarly as in the case of discrete lines, the first investigations were algorithmical. [START_REF] Bresenham | A linear algorithm for incremental digital display of circular arcs[END_REF]10,11,13 Discrete circles were considered as digitizations of euclidean ones. It is thus natural to ask whether J.-P. Reveillès' arithmetical approach is extendable to discrete circles, in order to supply them with a definition independent of euclidean geometry. If so, one can expect discrete circles to share topological and geometrical properties with euclidean circles: (i) the connectedness: a discrete circle is k-connected (see Definition 1), (ii) the thinness: a discrete circle is (1 -k)-minimal (see Definition 2), (iii) the tiling of the space: concentric discrete circles tile the discrete space Z 2 . 1 Such an extension has been proposed by É. Andres. [START_REF] Andres | Discrete circles, rings and spheres[END_REF] His discrete analytical circles, defined as subsets of Z 2 , verify a double diophantine inequality (see Definition 5) and are provided with a constant thickness, similarly to J.-P. Reveillès' arithmetical discrete lines. 14 In other words, a discrete analytical circle is nothing but a discrete ring, that is, a set a pixels at a bounded distance from an euclidean circle. In that case, one shows that the condition (iii) holds while (i) and (ii) are partially satisfied. [START_REF] Andres | Discrete circles, rings and spheres[END_REF] In the present paper, in order to deal with the properties of circles, we introduce an arithmetical definition of discrete circles with a non-constant thickness. In particular, we investigate different approximations of circles by discrete lines. Depending on the nature (naive, standard, . . . ) of these lines, we characterize new discrete circles, such as the k-minimal ones (see Definition 2), or already known discrete circles, such as the regular ones with integer parameters 1 or the Bresenham's ones. [START_REF] Bresenham | A linear algorithm for incremental digital display of circular arcs[END_REF] As far as we know, this is the first arithmetical characterization of the latter. This paper is sketched as follows:

• We begin with some recalls. First, Section 2 gives basic notions on discrete connectedness useful to fully understand the present paper. Second, The arithmetical discrete lines which inspire our approach are presented in Section 3. Finally, in Section 4, we focus on already known discrete circles. • In Section 5, we introduce a general arithmetical definition of discrete curves extending the artihmetical discrete lines. • In Section 6, we apply our general definition to the particular case of discrete circles. we characterize the already known circles of Section 4 and define new classes fo discrete circles.

BASIC NOTIONS

The aim of this section is to introduce the basic notions of discrete geometry which are used all along the present paper. Let n be an integer greater than 2 and let {e 1 , . . . , e n } denote the canonical basis of the euclidean vector space R n . Let us call discrete set any subset of the discrete space

Z n . If the point M = n i=1 M i e i ∈ R n , with M i ∈ R for each i ∈ {1, . . . , n}, then M is represented by (M 1 , . . . , M n ). A point of Z n is called a voxel. If n = 2, a voxel is also called a pixel.
Let N = (N 1 , . . . , N n ) and N = (N 1 , . . . , N n ) be two voxels and let k ∈ {0, . . . , n -1}. The voxels N and N are said to be k-neigbors if and only if:

N -N ∞ = max{|N 1 -N 1 |, . . . , |N n -N n |} = 1 and N -N 1 = n i=1 |N i -N i | ≤ n -k.
Definition 1 (k-connectedness) Let k ∈ {0, . . . , n -1}. A discrete set E is said to be k-connected if for each pair (N, N ) of voxels in E, there exists a finite sequence (N 1 , . . . , N p ) of voxels such that N = N 1 , N = N p and the voxels N j and N j+1 are k-neighbors, for each j ∈ {1, . . . , p-1}. For clarity issues, a 0-connected discrete set is said to be connected.

Let E be a discrete set, N be a voxel of E and k ∈ {0, . . . , n

-1}. The k-connected component of N in E is the maximal k-connected subset of E, for inclusion, containing N. The discrete set E is said to be k-separable if its complement E = Z n \ E has two k-connected components. Definition 2 (k-simple points, k-minimality) Let k ∈ {0, . . . , n -1}. A voxel N of a discrete set E is said to be k-simple in E if E \ {N} is k-separable. A k-separable discrete set without k-simple points is said to be k-minimal.
In the present paper, we work in the two-dimensional discrete space Z 2 for clarity issues but all the ideas and definitions presented extend in a natural way to Z n .

ARITHMETICAL DISCRETE LINES

In the early nineties, J.-P. Reveillès 14 has laid the foundations of a new geometry on Z n , the arithmetical discrete geometry. Its link with the usual euclidean geometry takes place at infinity * : a discrete lattice seen from an infinitely distant point seems to be continuous. Unfortunately, due to the discrete character of Z n , a single continuous object has several discrete representations, each one inheriting some of its properties. For instance, arithmetical discrete lines are defined as strips depending on a parameter w ∈ N, the so-called arithmetical thickness (see Figure 1

D(v, µ, ω) = (i, j) ∈ Z 2 | - ω 2 ≤ ai + bj + µ < ω 2 . ( 1
)
The choice of the inequalities in ( 1) is arbitrary, one thus can choose to define the arithmetical discrete line D(v, µ, w) of vector normal v = (a, b) ∈ R 2 , translation parameter µ ∈ R and arithmetical thickness ω ∈ R + as follows:

D(v, µ, ω) = (i, j) ∈ Z 2 | - ω 2 < ai + bj + µ ≤ ω 2 . (2)
The arithmetical discrete line D(v, µ, ω) (resp. D(v, µ, ω)) is called the lower (resp. the upper ) arithmetical discrete line of normal vector v ∈ R 2 , translation parameter µ ∈ R and arithmetical thickness ω ∈ R +

In fact, an arithmetical discrete line is the set of solutions of diophantine equations, that is, the ones with integer solutions. Given an arithmetical discrete line D(v, µ, ω), if v 1 and v 2 are rationnally dependent, then one can suppose 14 that v ∈ Z 2 , µ ∈ Z and ω ∈ N. Consequently, D(v, µ, w) is the set of solutions of at most ω diophantine equations. The arithmetical thickness w is also strongly related to the connectedness of arithmetical discrete lines (see Figures 1(b) and 1(c)). Definition 4 (Naive 14 and standard 6 discrete lines) An arithmetical discrete line D(v, µ, ω) is said to be naive (resp. standard

Theorem 1 (k-minimality and arithmetical discrete lines

14 ) The arithmetical discrete line D(v, µ, ω) is 0-minimal (resp. 1-minimal) if and only if ω = v 1 (resp. ω = v ∞ ). (i,j) (i-1,j+1) (i+1,j) (i',j') (i'-1,j'+1) (i',j'-1)
) if ω = v ∞ (resp. ω = v 1 ).
In other words, an arithmetical discrete line D(v, µ, ω) is naive (resp. standard) if and only if it is 1-minimal (resp. 0-minimal).

ALREADY KNOWN DISCRETE CIRCLES

In the present section, we go back over two different visions of discrete circles, an algorithmical one due to J. Bresenham [START_REF] Bresenham | A linear algorithm for incremental digital display of circular arcs[END_REF] and an analytical one due to É. Andres. 1

We first have a look at the algorithmical approach, the oldest and more studied one, before taking an interest in the analytical discrete circles, a vision arising from the arithmetical discrete geometry.

Algorithmical Approach

Required by digital plotters, many circle drawing algorithms, 4, 10, 11, 13 have been developped. All of them intend to find a thinnest discrete circle, ideally a 1-minimal one, and incrementally return the computed pixels. Moreover, most of them only compute in the first quadrant, the one where the pixels have all positive coordinates. The entire circle is deduced by axial symmetries.

Let us now remind J. Bresenham's algorithm [START_REF] Bresenham | A linear algorithm for incremental digital display of circular arcs[END_REF] for drawing discrete circles with centers in Z 2 and integer radii. Our choice is motivated by the fact that this discrete circle drawing algorithm is certainly the best known one. In addition, the outputs of the other ones 10, 11, 13 are close together.

In the first quadrant, a circle is a decreasing function of the abscissa. Consequently, in order to cover this quadrant clockwise, J. Bresenham's algorithm involves only three unitary moves, namely down, down-right and right (see Figure 2(a)).

Let R ∈ N and let us draw the Bresenham's circle B(0, R) of center O = (0, 0) and radius R, namely the digitization of the euclidean circle C(O, R) of center O = (0, 0) and radius R. Let ∆ be the following map:

∆ : N 2 -→ Z (i, j) -→ (i + 1) 2 + (j -1) 2 -R 2 .
Then, for each (i, j) ∈ N 2 , ∆(i, j) ≥ 0 if and only if the pixel (i + 1, j -1) is outside C(O, R).

The initial selected pixel is (0, R). Since R ∈ N , one deduces that the pixel (0, R) belongs to C(O, R).

Let (i, j) be a computed pixel in the first quadrant. J. Bresenham shows that ∆(i, j) ≥ 0 if and only if (i, j) is in the second octant of Z 2 , that is, the set {(i, j) ∈ Z 2 | 0 ≤ i ≤ j} and it remains two candidates for being the next pixel. Let us select the closest one to C(O, R). In the first octant, it is determinated by the difference δ 1 (i, j) between the distances from (i + 1, j) and (i + 1, j -1) to C(O, R):

δ 1 (i, j) = (i + 1) 2 + j 2 -R 2 -(i + 1) 2 + (j -1) 2 -R 2 .
(

) 3
If δ 1 (i, j) ≤ 0 then the pixel (i + 1, j) is selected, and so are the pixels (-i -1, j), (-i -1, -j), and (i + 1, -j) by symmetries. Otherwise, the pixels (i + 1, j -1), (-i -1, j -1), (-i -1, -j + 1), (i + 1, -j + 1) are selected.

In the second octant, the same holds for the pixels (i + 1, j -1), (i, j -1) and their symmetrical images.

As shown in the present section, J. Bresenham's algorithm allows to draw discrete circles of centers in Z 2 and integer radii (see Figure 2(b)). Most of the time, it returns a 1-minimal discrete circle (see Definition 2). Otherwise, the returned discrete circle contains exactly four 1-simple points located on the diagonals (spikes 9 or sharp corners 10 ). Z. Kulpa 9 characterized the concerned radii R:

R 2 = 2 R √ 2 2 - R √ 2 + 1. ( 4
)
whose first solutions are 4, 11, 134, 373, 4552, 12671, 154634, 430441, 5253004, 14622323, . . . .

Discrete Analytical Circles 1

The algorithmical approach of discrete circle drawing 4, 10, 11, 13 let several questions unsolved. For example, they do not compute a set of discrete circles tiling the discrete space Z 2 and do not allow to check whether a point belongs to a given discrete circle. É. Andres fills partially this lack by introducing the discrete analytical circles. [START_REF] Andres | Discrete circles, rings and spheres[END_REF] Roughly speaking, similarly to the arithmetical discrete lines (see Definition 3), the discrete analytical circles are defined as rings depending on a parameter ω ∈ R + called the arithmetical thickness (see Figure 3). Definition 5 (Discrete analytical circles 1 ) The discrete analytical circle C(M, R, ω) of center M 0 = (x 0 , y 0 ) ∈ R 2 , of radius R ∈ R + and of arithmetical thickness ω ∈ R + , is the following subset of Z 2 defined by:

C(M 0 , R, ω) = (i, j) ∈ Z 2 | R - ω 2 2 ≤ (i -x 0 ) 2 + (j -y 0 ) 2 < R + ω 2 2 . (5)
One obtains limited results on the relation between the arithmetical thickness and the connectedness of a discrete analytical circle.

Proposition 1 (Discrete analytical circles and connectedness

2 ) Let C(M 0 , R, ω) be a discrete analytical circle. If ω ≥ 1, then C(M 0 , R, ω) is connected. Moreover, for each point M 0 ∈ Z 2 and for each R ∈]0, 1], the set {C(M 0 , R + k, 1), k ∈ N} tiles the discrete space Z 2 .
Definition 6 (Regular discrete circles 2 ) Let C(M 0 , R, ω) be a discrete analytical circle. If w = 1 then C(M 0 , R, ω) is said to be regular (see Figure 3(b)).

One can notice that, contrary to the thickness of an arithmetical discrete line, the one of a discrete analytical circle does not entirely characterize its connectedness. For instance, if ω < 1, then we cannot state in the general case whether C(M 0 , R, ω) is connected. Moreover, when ω ≥ 1, we cannot conclude whether C(M 0 , R, ω) is 0-separable, 0-minimal or 1-connected.

(x 0 , y 0 ) R R + 1 2 R -1

ARITHMETICAL DISCRETE CURVES: A TANGENTIAL POINT OF VIEW

In Section 4, we notice that the algorithmical (Section 4.1) and analytical (Section 4.2) approaches do not characterize the same discrete circles. On the one hand, the Bresenham's discrete circles are often 1-minimal but do not tile the discrete space Z 2 . On the other hand, discrete analytical circles tile Z 2 but their thinness is not entirely controled by their parameters. On can see that discrete circles defined by J. Bresenham 4 and É. Andres 1 are quite far to satisfy the conditions (i), (ii) and (iii) page 1.

From theorem 1, given an arithmetical discrete line D(v, µ, ω) with normal vector v = (a, b) ∈ R 2 , the kconnectedness and the k-minimality of D(v, µ, ω) only depend on its thickness, expressed as a function of a and b, namely the coordinates of the vector v. In fact, those coordinates are nothing but the partial derivatives of the linear form ax + by + µ. In other words, v is the normal vector of the tangent to the euclidean line with equation ax + by + µ = 0 at each of its points. Thus, in order to define relevant discrete analogues of non-linear forms, it becomes natural to consider a differential approach. For clarity issues, let us first introduce notation.

Notation. -Let f : R 2 -→ R be a differentiable function. We denote by ∂ x f : R 2 -→ R and ∂ y f : R 2 -→ R the following functions:

∂ x f : R 2 -→ R and ∂ y f : R 2 -→ R (x, y) -→ ∂f ∂x (x, y) (x, y) -→ ∂f ∂y (x, y) Definition 7 (Arithmetical discrete curves) Let f ∈ R[x,
y] be a polynomial in two indeterminates over R and let ω : R 2 -→ R. The arithmetical discrete curve C(f, ω) of equation f (x, y) = 0 and thickness map ω is defined by:

C(f, ω) = (i, j) ∈ Z 2 | - ω (∂ x f (i, j), ∂ y f (i, j)) 2 ≤ f (i, j) < ω (∂ x f (i, j), ∂ y f (i, j)) 2 . ( 6
)
Let us notice that the arithmetical discrete line of normal vector v = (a, b) is canonically an arithmetical discrete curve and, for all (i, j)

∈ Z 2 , ω (∂ x f (i, j), ∂ y f (i, j)) = ω(a, b)
and ω can be replaced by a constant function.

As noticed for arithmetical discrete lines, the choice of inequalities in ( 6) is arbitrary and one can choose the opposite convention. The arithmetical discrete curve C(f, ω) of equation f (x, y) = 0 and thickness map ω is defined by:

C(f, ω) = (i, j) ∈ Z 2 | - ω (∂ x f (i, j), ∂ y f (i, j)) 2 < f (i, j) ≤ ω (∂ x f (i, j), ∂ y f (i, j)) 2 . ( 7
)

A PARTICULAR CASE: THE ARITHMETICAL DISCRETE CIRCLES

In the present Section, we focus on a particular case of arithmetical discrete curves. In fact, from the quadratic form f , (x -x 0 ) 2 + (y -y 0 ) 2 -R 2 for all (x, y) ∈ R 2 and (x 0 , y 0 , R) ∈ R 2 × R + , we investigate the arithmetical discrete circles.

General Definition

Definition 8 (Inner and outer discrete circles) Let M 0 = (x 0 , y 0 ) be a point of the discret space Z 2 , R ∈ R + and ω : R 2 -→ R.

The inner discrete circle C(M 0 , R, ω) of center M 0 , radius R and thickness function ω is the set of point (i, j) with integer coordinates defined by:

C(M0, R, ω) =  (i, j) ∈ Z 2 | - ω (2i, 2j) 2 ≤ f (i, j) < ω (2i, 2j) 2 ff . (8)
The outer discrete circle C(M 0 , R, ω) of center M 0 , radius R and thickness function ω is the set of point (i, j) with integer coordinates defined by:

e C(M0, R, ω) =  (i, j) ∈ Z 2 | - ω(2i, 2j) 2 < f (i, j) ≤ ω(2i, 2j) 2 ff . ( 9
)
We select appropriate thickness function to highlight relevant classes of arithmetical discrete circles. In the one hand, we provide the regular discretes circle 1 of integer parameters and the Bresenham's circles 4 with aritmetical definitions. In the other hand, we characterize new discrete circles, the 1-minimal and the 0-minimal ones.

Let us bring up again the relevant classes of arithmetical discrete lines (see Definition 4). The thickness of a naive (resp. standard) discrete line of normal vector v is the v ∞ (resp. v 1 ). The norms seems to control the connectedness of the arithmetical discrete lines. Moreover, the definition of the discrete analytical circles imply the same idea with the

• 2 : R - 1 2 ≤ (i -x 0 , j -y 0 ) 2 < R + 1 2 ⇔ R - 1 2 2 ≤ (i -x 0 ) 2 + (j -y 0 ) 2 < R + 1 2 2 .
Notation. -We denote by ω k : Z 2 -→ R + , the thickness function, related to the norm • k , as following:

ω k : Z 2 -→ R + (i, j) -→ ∂f (i,j)
∂i , ∂f (i,j) ∂j k .

Characterizations of Already Known Discrete Circles

The thickness function ω 2 provides the regular discrete circle 1 of integer parameters with a new arithmetical definition (see Figure 4(a)).

Proposition 2 (Regular discrete circles) Let N 0 = (i 0 , j 0 ) ∈ Z 2 and R ∈ N. The regular discrete circle 1 C(N 0 , R, 1) is the arithmetical discrete circle C(N 0 , R, ω 2 ) but also the arithmetical discrete circle C(N 0 , R, ω 2 ):

C(N 0 , R, ω 2 ) = C(N 0 , R, 1) = C(N 0 , R, ω 2 ) ( 10
)
Proof. On the one hand, the regular discrete circle 1 C(N 0 , R, 1) is characterized as follow:

C(N 0 , R, ω) = (i, j) ∈ Z 2 | R - 1 2 2 ≤ (i -i 0 ) 2 + (j -i 0 ) 2 < R + 1 2 2 .
Since N 0 = (i 0 , j 0 ) ∈ Z 2 and R ∈ N, the double inequality defining the set above amounts to: -R < (i

-i 0 ) 2 + (j -i 0 ) 2 -R 2 < R + 1.
On the other hand, the arithmetical discrete circle C(N 0 , R, ω 2 ) is characterized as follow:

C(N 0 , R, ω) = (i, j) ∈ Z 2 | -(i -i 0 ) 2 + (j -j 0 ) 2 ≤ (i -i 0 ) 2 + (j -j 0 ) 2 -R 2 < (i -i 0 ) 2 + (j -j 0 ) 2 .
Let k = (i -i 0 ) 2 + (j -j 0 ) 2 and g 1 : N -→ R be the map:

g 1 : R -→ R k -→ k 2 + k -R. g 1 is an increasing function. Moreover, g 1 (R -1) < 0 and g 1 (R) > 0. Let k min ∈ R such that g 1 (k min ) = 0. We have g 1 (R -1) < g 1 (k min ) < g 1 (R) and then, R < k min < R -1. Since (i -i 0 ) 2 + (j -i 0 ) 2 -R 2
is an integer, we look for integer bounds. Such an assumption leads to the lower bound -R < (i

-i 0 ) 2 + (j -j 0 ) 2 -R 2 .
Similarly, the upper bound (i

-i 0 ) 2 + (j -j 0 ) 2 -R 2 < R + 1 is obtained and finally, C(N 0 , R, 1) = C(N 0 , R, ω 2 ).
We can in a similar way prove C(N 0 , R, 1) = C(N 0 , R, ω 2 ).

The thickness function ω ∞ provides the Bresenham's circles with a new characterization. As far as we know, this is the first arithmetical definition of them. The original definition of the Bresenham's circle compute only circles of integer radius. Here, we use its extension to integral square radius introduced by M. McIlroy. 10 Proposition 3 (Arithmetical characterization of the Bresenham's circle)

Let N 0 = (i 0 , j 0 ) ∈ Z 2 and R such that R 2 ∈ N. The arithmetical discrete circle C(N 0 , R, ω ∞ ) is the Bresenham's circle B(N 0 , R): C(N 0 , R, ω ∞ ) = B(N 0 , R) (11)
Proof. First of all, The arithmetical discrete circle C(N 0 , R, ω ∞ ) contains the pixel initially selected by the bresenham's algorithm, (0, R). Now we show that the both definitions respect the same criterion to determine which will be the next pixel of the discrete circle. Let N = (i, j) ∈ N 2 be a pixel of the arithmetical discrete circle C(N 0 , R, ω ∞ ) such that 0 ≤ j ≤ i, namely in the second octant. Let consider that we cover the second octant clockwise. Because of our definition, and with I = i -i 0 , J = j -j 0 , two pixels are reachable: Thus, the sign of the term I 2 + J 2 -R 2 + 2I -J + 1 determined the next pixel to be selected.

• the pixel (i + 1, j -1) belongs to C(N 0 , R, ω ∞ ) if and only if -2I + J -1 ≤ I 2 + J 2 -R 2 , • the pixel (i + 1, j) belongs to C(N 0 , R, ω ∞ ) if and only if I 2 + J 2 -R 2 < -2I + J -1.
Bresenham's algorithm uses as selection criterion the minimization of the distance between the euclidean circle and the next pixel. It compares the two distances:

δ 2 (i, j) = (I + 1) 2 + (J) 2 -R 2 + (I + 1) 2 + (J -1) 2 -R 2 δ 2 (i, j) 2 = I 2 + J 2 -R 2 + 2I -J + 3 2 (12)
I 2 + J 2 -R 2 + 2I -J is an integer, thus the sign of the term I 2 + J 2 -R 2 + 2I -J + 1 is equivalent to the one of I 2 + J 2 -R 2 + 2I -J + 3 2 and determines the next pixel to be selected. Similarly, in the first octant, the selection criterion is the same. Then, starting from the same initialisation and computing with the same criterion, C(N 0 , R, ω ∞ ) = B(N 0 , R).

New Classes of Discrete Circles

Now, we introduce two new definitions of discrete circles related to their k-minimality.

Definition 9 (Naive discrete circles) By analogy with discrete lines, 1-minimal discrete circles are said to be naive.

Definition 10 (Standard discrete circles) By analogy with discret lines, 0-minimal discrete circles are said to be standard.

The thickness function ω ∞ allows to define naive discrete circles (see Figure 4(b)).

Proposition 4 (Naive discrete circles)

Let N 0 = (i 0 , j 0 ) ∈ Z 2 , R ∈ N and ∆ = R 2 -2 R √ 2 2 + R √ 2 .
1. The inner arithmetical discrete circle C(N 0 , R, ω ∞ ) is naive if and only if ∆ = 0.

2. The outer arithmetical discrete circle C(N 0 , R, ω ∞ ) is naive if and only if ∆ = 1.

So, at least one of them is a naive discrete circle.

Proof. First, we investigate the inner arithmetical discrete circle C(N 0 , R, ω ∞ ). Since it is the Bresenham's circle B(N 0 , R), it is 0-minimal except for radii verifying the equality 4. Then, the equivalence 1 is true.

Second, we investigate the outer arithmetical discrete circle C(N 0 , R, ω ∞ ). we have

C(N 0 , R, ω ∞ ) = C(N 0 , √ R 2 + 1, ω ∞ ). Hence, C(N 0 , R, ω ∞ ) is the Bresenham's circle B(N 0 , √ R 2 + 1) of center N 0 and radius √ R 2 + 1 ( √ R 2 + 1 2 is
an integer). One shows this discrete circle is 1-minimal except for radii verifying:

R 2 = 2 R √ 2 2 -R √ 2 .
Thus, the equivalence 2 is true.

Similarly, the thickness function ω 1 allows to define standard discrete circles (see Figure 4(c)).

Proposition 5 (Standard discrete circles) Let M 0 = (x 0 , y 0 ) ∈ R 2 and R ∈ R + . The arithmetical discrete circles C(M 0 , R, ω 1 ) and C(M 0 , R, ω 1 ) are standard discrete circles.

Proof. The arithmetical discrete circle C(M 0 , R, ω 1 ) includes the arithmetical discrete circles C(M 0 , R, ω ∞ ), since for all pixel N = (i, j), ω 1 (2i, 2j) > ω ∞ (2i, 2j), and is at least 1-minimal. With no loss of generality, let consider we cover clockwise the first quadrant of the inner arithmetical discrete circle C (M 0 , R, ω 1 ). A pixel N = (i, j) ∈ C -(M 0 , R, ω 1 ) satisfies the following double diophantine inequality:

-X -Y ≤ X 2 + Y 2 -R 2 < X + Y,
where X = i -x 0 and Y = j -y 0 . Under our assumptions, only three pixels, (i + 1, j), (i + 1, j -1) and (i, j -1), are canditates for being the next clockwise one belonging to the discrete circle C(M 0 , R, ω 1 ). It is sufficient that the pixel (i + 1, j) or the pixel (i, j -1) belongs to the discrete circle to have 0-separability (1-connectedness) and it is necessary that the both do not belong simultaneously to the discrete circle to avoid 0-simple points. From equation 6.3, we obtain the membership for respectively the pixels (i + 1, j) and (i, j -1) of the circle:

X -Y + 1 ≤ (X + 1) 2 + Y 2 -R 2 < 3X + Y + 1 -X -3Y + 1 ≤ X 2 + (Y -1) 2 -R 2 < X -Y + 1
The two inequalities could not be simultaneously valid, but one of the both is always verified. So, standard discrete circles are 0-minimal.

Same study leads to same results for outer standard discrete circles.

CONCLUSION

In the present paper, we introduce a general arithmetical characterization of discrete curves including and extending the definition of the arithmetical discrete lines. 14 The key point of our approach is the thickness function. As far as we now, this parameter has always been considered as a constant. It leads to a wide knowledge of discrete line 14 and partially improves the knowledge of discrete circles. [START_REF] Andres | Discrete circles, rings and spheres[END_REF] In the present paper, we decide to consider the thickness as a function of the local derivatives of the curve. It doesn't throw results on discrete lines back into question and improves again the characterization of discrete circles.

In fact, in the particular case of discrete circles with integer parameters, our definition includes the Bresenham's circles 4 and the analytical discrete circles. [START_REF] Andres | Discrete circles, rings and spheres[END_REF] Moreover, we reach new classes of discrete circles, the standard and naive ones. Similarly to arithmetical discrete lines, 14 all these characterizations are obtained by selecting the usual norms as thickness functions. We achieve the aim we state at the beginning of the present paper, namely define discrete circles sharing properties i, ii and iii (see page 1) with the euclidean circle.

In the more general case of discrete circles with real parameters, our definition only characterize standard discrete circles. But the definition of regular discrete circles by E. Andres 1 is valid for any center with real coordinates and any real radius, just as Bresenham's circles 4 since S. Pham 12 extended them. Tangents can be not an enough accurate approximation of circles. In our next papers we will try to overcome such a limit. Finally, our approach is valid for polynomial curves and does not dependent on the two-dimensional space. It is thus natural to investigate these both directions. regular discrete representations of most of the curves do not make any sense.then, Figure 5 shows only naive and standard discrete ellipses derivated from our definition. In Figure 5(a), naive discrete ellipses are drawn. In the more general case of polynomial curves, our definition does no characterize naive discrete representation since the two candidates C and C, can both have 1-simple points. In figure 5(b), one of the drawing contains 0-simple points and is no more a standard discrete representation as expected. This problem comes from the poor resolution of the discrete lattice compared to the local curvature of the euclidean ellipse. Figure 6 presents some results on spheres. It is possible to recover the regular discrete spheres 2 of E. Andres, figure 6(b) or, to define the thinnest discrete spheres, mamely the 1-minimal ones, figure 6(a). In the same way, the different hyperspheres, naive, regular and standard, are reachable.

  ω ax + by + µ = ω ax + by + µ = 0 (a) General definition. (b) Naive discrete lines. (c) Standard discrete lines.

Figure 1 .

 1 Figure 1. Arithmetical discrete lines.

Definition 3 (

 3 (a)). Arithmetical discrete lines 14 ) The arithmetical discrete line D(v, µ, ω) of normal vector v = (a, b) ∈ R 2 ,translation parameter µ ∈ R and arithmetical thickness ω ∈ R + is the subset of Z 2 defined by:

  (a) Incremental process of building (b) Bresenham's circles of integer radii between 1 and 10.

Figure 2 .

 2 Figure 2. Bresensenham's circles.

2 ω

 2 (a) General definition. (b) Concentric regular discrete circles of real coordinate center.

Figure 3 .

 3 Figure 3. Discrete analytical circles.

  (a) Regular discrete circles of integer radii between 1 et 10. (b) Naive discrete circles of integer radii between 1 and 10. (c) Standard discrete circles of radii 1, 3, 5, 7 and 9.

Figure 4 .

 4 Figure 4. Some classes of discrete circles we obtain.

  (a) Naive discrete ellipses. (b) Standard discrete ellipses.

Figure 5 .

 5 Figure 5. Discrete Ellipses with integer parameters.

  (a) Eighth of thinnest discrete spheres. (b) Regular discrete spheres tiling the space.

Figure 6 .

 6 Figure 6. Spheres with integer parameters.

* ". . . le désaccord entre le Discret et le Continu se mesure par un éloignement "infini" pris en compte dans l'un des axiomes fondateurs de la Mathématique. Le concept de nombre entier infiniment grand nous permet d'envisager une Mathématique Discrète Idéale intermédiaire entre la classique Mathématique du Continu et la Mathématique Discrète ; comme la pile d'un pont entre deux rives. . . " 14


Eloy Romero

Christoph Kowitz

Jose Enrique Roman

Florian Merz

Frank Jenko

Frank Jenko Preconditioning

Preconditioning Eigensolvers for Computing Unstable Modes of Plasma

Preconditioning Eigensolvers for Computing Unstable Modes of Plasma

List of authors: E. Romero1 C. Kowitz2 J. E. Roman3 F. Merz4 F. Jenko 5Numerical simulations are essential tools in the study of the processes occurring in magnetically confined plasmas, such as the anomalous transport (large amounts of energy and particles are transported out of the hot core zone of the plasma) which is the result of the microinstabilities in the plasma produced by gradients in the temperature, magnetic field and density. These microinstabilities can be characterized by the analysis of the unstable modes in the linearized gyrokinetic problem,

∂g ∂t = L[g], (1)
where L is a time independent, complex, non-Hermitian integro-differential operator. Time integration of (1) can only obtain the dominant unstable mode (with the largest growth rate). However, for detailed analyses also the subdominant unstable modes are required, which can be obtained via the partial eigendecomposition of L.

GENE [START_REF] Dannert | Gyrokinetic simulation of collisionless trapped-electron mode turbulence[END_REF] is a massively parallel code for plasma simulations that solves the linear and nonlinear gyrokinetic equations following an Eulerian approach. In GENE, the operator matrix is never computed explicitly, but implemented in a highly parallelized and efficient matrix-free form, allowing to deal with problems of sizes ranging from several hundred thousand for linear simulations up to a few billion for nonlinear problems, for reasonably accurate models.

In GENE, the eigenvalue computations of the linear operator are carried out by means of SLEPc, the Scalable Library for Eigenvalue Problem Computations [START_REF] Hernandez | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF]. SLEPc provides a collection of parallel iterative eigensolvers on top of PETSc [START_REF] Balay | PETSc users manual[END_REF], some of them allowing the use of a preconditioner for accelerating the convergence.

Due to the difficulty of building an efficient preconditioner based on a dense operator matrix, which is not available in explicit form, previous works addressed the problem without preconditioning. In particular, moderate success was obtained with the Krylov-Schur method combined with the shift-and-invert spectral transformation [START_REF] Roman | Fast eigenvalue calculations in a massively parallel plasma turbulence code[END_REF] as well as employing the Jacobi-Davidson method [START_REF] Romero | A parallel implementation of the Jacobi-Davidson eigensolver and its application in a plasma turbulence code[END_REF]. In this work, we analyze different preconditioning strategies in the context of GENE.

The discretized operator L is obtained with a combination of spectral and finite difference techniques. As a consequence, it can be split into a dense part L φ and a sparse part L g . The latter matrix has a block diagonal sparse multi-banded structure because the differentiation is done locally.

This work considers the utilization of preconditioners based on the finite-difference part of the operator, L g , because it contains a representative part of the full operator and its structure is adequate for building efficient preconditioners. We will present results for the inexact shiftand-invert Krylov-Schur and Jacobi-Davidson methods, with several preconditioners including domain decomposition techniques (block Jacobi and additive Schwarz) with incomplete factorizations in each subdomain. We analyze the robustness and performance of the solvers in realistic plasma containment devices such as tokamaks and stellarators.

Universidad Politécnica de Valencia, Instituto I3M, C. de Vera s/n, 46022 Valencia, Spain, elroal@upv.es

Max-Planck-Institut für Plasmaphysik

Universidad Politécnica de Valencia

Max-Planck-Institut für Plasmaphysik

Max-Planck-Institut für Plasmaphysik

Acknowledgement We acknowledge computer resources provided by the Barcelona Supercomputing Center (BSC). This work was supported by the Spanish Ministerio de Ciencia e Innovación under project TIN2009-07519.


Gwenol

Scalable Preconditioners for Transient Navier-Stokes Equations in Hemodynamics Simulations

teaching and research institutions in France or abroad, or from public or private research centers.

Scalable Preconditioners for Transient Navier-Stokes Equations in Hemodynamics Simulations

List of authors:

G. Grandperrin 1 A. Quarteroni 2

S. Deparis 3

We are interested in the resolution of transient Navier-Stokes equations on large scale parallel architectures. For this reason, we are developing preconditioners for the linear system associated to the linearized problem. Weak and strong scalabilities are not the only measures for the parallel performances. In fact, the number of iterations to solve the linear system shall not increase with increasing number of parallel tasks. This can be achieved by devising specific scalable preconditioners.

We present a family of Pressure-Convection-Diffusion (PCD) preconditioners and their parallel performances. We also reformulate the preconditioned linear system in order to lower the condition number of the system. We then compare the scalability of the presented methods.

The implementation is done with LifeV (http://www.lifev.org), a finite element library for High Performance Computing (HPC) based on Trilinos (http://trilinos.sandia.gov). All the methods are illustrated using simulations relevant to hemodynamics.

CMCS, EPFL, Switzerland

CMCS, EPFL, Switzerland

CMCS, EPFL, Switzerland 1


Jérémie Schiettecatte

Settlement Process in Ancient Hadramawt

). This area was characterized by a common language and the dominant worship of the deity Sayîn

The evolution of local populations was conditioned by their very specific environmental contexts, that is by their irrigation capacities and the opportunities to interact with their neighbours. The first irrigation systems appeared towards 3000 BC, mostly made of simple structures such as diversion walls on the slopes or check-dams 2 . The consequences of an increase of the aridity during the 4 th and 3 rd millennium BC combined to the drying-up of some springs and the shrinking of pasture areas could have only been twofold: the gathering of populations down to the basins, the lower valleys or close by the rare perennial sources on the first hand; and on the other hand, the development of sheep and goat farming, to the detriment of cattle, allowing a greater mobility and an easier access to scattered resources.

During the 3 rd and 2 nd mill. BC, as a response to the environmental stress caused by increasing aridity, small cultivated areas were exploited. The irrigation structures were not widely used as long as the karstic network and the last active springs compensated the water requirements.

) shows a town together with one or several temples located at every confluence of the Masîla with its major tributaries. The steeply sided setting of the wâdî Masîla into the silt plain made it impossible for local population to use its water to supply the irrigation systems. Indeed the

water supply for cultivation did not came from the main wâdî but rather from the diverting of floods drained out of tributary valleys. Communities were settled where the tributaries open out so that they could control its access, handle the floods and develop an irrigation system that stretched over both the central valley and upstream from the tributary mouth. Each of these archaeological sites was the centre of a cultivated area. The growth of these various settlement depends partly -not to say mainly -on the size of the drainage basin watering the cultivated land.

Thus, the settlement pattern in the valley depended on environmental constraints that in turn affected the urban growth in this region: the growth of each site was quickly limited by the agricultural capacity that the diverted floods made possible. The genuine towns were therefore extremely rare and most of the settlements located at the mouth of tributaries resembled rather villages than towns. Moreover, according to the map of the region, each cultivated area allowed only a single major settlement. Due to their limited capacities, the runoffs ensured very seldom an important density of the settlement network. The rare intermediate facilities between the main sites were but farms or hamlets exploiting the land.

In the western part of the Hadramawt valley, the settlement pattern differs from the one above-mentioned. This is partly due to a different landscape and to the greater width of the Hadramawt valley. According to the available data (fig. 2), the seldom urban settlements did not expand at the mouth of tributaries of the wâdî Hadramawt but rather on the middle/lower reaches of the main tributaries. There, flood was easier to divert and control. Set back from the main valley, none of these sites display sign of an influence from a centralized power. These towns do not seem having properly integrated the political life of the kingdoms located around the Ramlat as-Sab'atayn however, at least until the end of the 1st mill. BC. Political conflict and power struggle certainly took place among the inhabitants of these sites or between neighbouring settlements. But that no supra-communal intervention is mentioned before the end of the 1 st mill. BC says a lot about the stability of political territories and land-holdings. It seems possible to assume that their extension mainly -not to say only -depended on technical capacities and on the time required to reach the fields.

Many sites located between Shabwa and Makaynûn area offered a defended fallback position for inhabitants of the area. Nevertheless, these defensive systems were seldom able to endure a siege. They were most probably designed to ensure the security of people against facing with raids or neighbourhood disputes. During the 1 st mill. BC, the kingdoms bordering Hadramawt do not seem to have felt drawn to the Inner valley; the Hadrami king himself did not leave a mark there. This region does not appear as a potential target before the end of the 1 st mill. BC. Its agricultural production capacities cannot feed more than small communities. As far as we know, no mineral resources have been exploited. These unassuming communities acted as a weak magnet on their neighbours and their geographical enclosing probably did not encourage these neighbours to raid them except maybe for nomadic tribes from the plateau. There would be no need for these settlements to get strong walls in view of improbable foes before the end of the 1 st mill. BC at least.

Until the end of the 1 st mill. BC, the Hadramawt entity is defined by a cultural unity based on common language, worships and material culture. In spite of this unity, inner Hadramawt did not have integrated the political life of the Hadramawt kingdom. This kingdom was centred on western Hadramawt (Shabwa, al-Binâ', Naqb al-Hagar), on the fringe of the Ramlat as-Sab'atayn (fig. 2). About the 2 nd and 1 st cent. BC, several conflicts set Hadramawt and Qatabân kingdoms against one another. Many settlements disappeared at that time. Some of them went obviously through dramatic events (the fire at Raybûn for example) ; others declined quickly without any historical explanation to be put forward such as Makaynûn.

On the opposite, the desert and coastal edges of Hadramawt merged quickly into a territorial entity governed from a single capital, Shabwa (fig. 2). Since the 7 th cent. BC, the settlements in the wâdî Jirdân constituted a possible resting place on the track from Shabwa to Qatabân kingdom. From the 4 th -3 rd cent. BC onwards, southern and coastal towns made their appearance. From this period at the latest, Hadramis settled on the Omani shore, at Khawr Rûrî; the Hadrami king made a kabîr -go-between of the king -responsible of the building of fortifications at Naqb al-Hajar, a stopping place between the shore and the Hadrami capital city. Moreover, a settlement is known at Shihr-East8 , an inhabited mooring place. These facts could be explained by king's will to control new maritime commercial routes that developed at that time with an increase of coastal shipping practice. Frankincense from Zufâr may have been transported by coastal shipping and then by land along the wâdî Mayfa'a, via the new grown city of Naqb al-Hajar, up to the capital Shabwa and later towards the north by the caravan road.

During the 1st cent. BC, the port of Bi'r 'Alî was founded and the pass of al-Binâ' fortified. These facts emphasize the primacy of the maritime route for incense handling from the growing areas of Zufâr and Mahra to Shabwa. By the middle of the 1st cent. AD, while Khawr Rûrî was rebuilt in its entirety and Bi'r 'Alî initiated a quick growth, the maritime trade definitely overtook the caravan trade.

Therefore, contrary to inner sites, most of the Hadrami settlements of the desert edge and coastal plain of the kingdom were soon integrated in a territory that the king controlled directly, because of their strategic function near the border of the Qatabânian kingdom and their importance on the trade routes. At the beginning of the Christian era, the power of the king significantly increased, thanks to the collapse of the Qatabânian kingdom. The monarch was addressed as mukarrib in the inscriptions, a title taken by the king when its kingdom had become the most powerful of South-Arabian kingdoms. He extended its control over the aromatics' trade. It is only at that very moment that we get clues of the integration of inner Hadramawt into this political entity.

Fig. 1 :Fig. 2

 12 Fig. 1 : Relief of Hadramawt

A. Rougeulle, 1999. « Coastal settlements in southern Yemen: the 1996-1997 survey expeditions on the Hadramawt and Mahra coasts », Proceedings of the Seminar for Arabian Studies 29, pp. 123-135.