Charles Poli

Dmitry V Savin

Olivier Legrand

Fabrice Mortessagne

Statistics of eigenfunctions in open chaotic systems: a perturbative approach

Keywords: numbers: 05.45.Mt, 03.65.Nk, 05.60.Gg

I. INTRODUCTION

In the domain of wave or quantum chaos [START_REF] Stöckmann | Quantum Chaos: an Introduction[END_REF], open systems are currently actively investigated both from experimental and theoretical points of view (see Refs. [START_REF] Kuhl | [END_REF]3] for recent reviews). Openness may be due to various physical mechanisms such as bulk absorption, coupling to the environment through physical channels as well as dissipative or radiative boundary conditions. Whatever the mechanism, openness results in spectral broadening ranging from the perturbative regime of non-overlapping (isolated) resonances to the so-called Ericsson regime of strong overlap. These mechanisms and their related spectral effects have been experimentally studied in various context: in microwave cavities [4,5,6,7,8], in optical microcavities [9,10,11], and in elastodynamics [12,13].

The most salient feature of open systems is the set of resonances which are quasibound states embedded in the continuum. A natural way to address them analytically is via the energy-dependent scattering matrix, S(E). Following the Heidelberg approach [14], the poles (i.e., resonances) of the S matrix turn out to be the complex eigenvalues of an effective non-Hermitian Hamiltonian H eff , whereas the bi-orthogonal eigenvectors of the latter determine the corresponding resonance states (quasimodes). Universal properties of resonance scattering in the chaotic regime can then be analyzed by applying random matrix theory (RMT) that amounts to replacing the actual non-Hermitian Hamiltonian with an RMT ensemble of the appropriate symmetry class [15]. The main advantage of such an approach is that it treats on equal footing both the spectral and scattering characteristics of open chaotic systems as well as that it is flexible enough to incorporate other imperfections of the system, e.g., disorder and losses [3].

By now, complex eigenvalues of such non-Hermitian random matrices have been studied quite systematically [16,17,18]. However, the statistical properties of the corresponding (left and right) eigenvectors are less understood. Quite a substantial progress in this direction has been achieved by Schomerus et al. [19], who studied mainly the systems with broken time-reversal symmetry. Other analytical results for a few physically interesting particular cases have also been re-ported in the literature recently [20,21,22,23]. Components of eigenvectors appear as residues of the S matrix at resonance positions and the understanding of their properties is thus important for many applications. For example, nonorthogonality of resonance eigenstates yields the enhancement (the socalled Petermann factor) of the line width of a lasing mode in open resonators [19] and influences branching ratios of nuclear cross-sections [24,25]. It features also in the particle escape from the scattering region [26] as well as in dissipative quantum chaotic maps [27].

This paper focuses on spectral and eigenvector statistics of such non-Hermitian random matrices describing effective Hamiltonians of open chaotic wave systems whose closed limit displays time-reversal symmetry (TRS). In this case, the quasimodes correspond to the complex-valued eigenvectors of H eff . To characterize this complexness, it is convenient to introduce [12,28] the ratio of the variances of the imaginary and real parts of the eigenvector as a single statistical parameter, hereafter called the complexness parameter [23]. One should note that this parameter is characteristic of the degree of non-orthogonality of the complex modes and, therefore, is closely connected to the Petermann factor mentioned above [19]. Other studies have considered the phase rigidity, another related parameter, introduced to characterize the degree to which a general scattering wave function is complex [22,29]. Both parameters are straightforwardly deduced from one another when the phase rigidity is calculated for a single eigenvector. The main advantage of considering the complexness parameter is to reveal a physical connection between spatial and spectral statistics [5,23].

In what follows, we study the probability distribution of the complexness parameter for a generic weakly open chaotic system and its connection with the distribution of resonance widths within the RMT approach. At the first stage, we derive an expression for the complexness parameter in the weak coupling regime and establish a general relation between its average and width fluctuations. Then accounting for the essential statistical feature of spectra in chaotic systems, namely, spectral rigidity, we investigate the case of a system whose closed limit is described by a pure picket-fence spectrum. An exact analytic prediction for the probability distribution of the complexness parameter is derived, depending on only two parameters: the number of open scattering channels and the mean resonance width. Finally, we consider the more realistic case of systems modeled by the Gaussian orthogonal ensemble (GOE). We derive an analytic expression for the probability distribution of the complexness parameter in this case and discuss the effect of spectral fluctuations.

II. EFFECTIVE HAMILTONIAN FORMALISM

A. Scattering approach

Open wave systems are commonly described using the socalled projection formalism [25,[START_REF] Mahaux | Shell-Model Approach to Nuclear Reactions[END_REF]. The exterior coupling is modeled by M scattering channels connected to N levels of a closed system. The coupling to the environment turns modes, with a infinite life time, into resonances, with a finite life time. Being initially introduced in nuclear physics, this formalism has been later applied successfully to wave billiards [START_REF] Dittes | [END_REF] for which antennas and absorption can be described by scattering channels [32]. In this approach, the resonance part of the Smatrix is given by:

S(E) = 1 -iV † 1 E -H eff V, (1)
where V is the coupling matrix of size N × M , the elements V c n of this matrix couple the n-th level to the c-th scattering channel. The poles of S are given by the eigenvalues of H eff . Assuming an independence of the coupling elements from the energy and neglecting direct processes [14], the effective Hamiltonian of the open systems is represented as follows:

H eff = H - i 2 V V † , ( 2
)
where H is the Hamiltonian of the closed system and the anti-Hermitian part i 2 V V † describes coupling to the channels. In the case of the systems with preserved TRS considered below, H is a real symmetric matrix of size N × N and V is also real. As usual, the limit N → ∞ is to be finally taken.

Since H eff is a non-Hermitian operator, the eigenvalue problems H eff |ψ n = E n |ψ n and ψn |H eff = E n ψn | define two sets of a priori independent eigenvectors, called right {|ψ n } and left { ψn |} eigenvectors associated to the same set of eigenvalues {E n }. These eigenvectors form a bi-orthogonal set which satisfies conditions of orthogonality, ψn |ψ m = δ nm , and completeness, n |ψ n ψn | = 1. Making use of the right eigenvectors, the diagonalization of H eff then reads:

H eff |ψ n = (E n - i 2 Γ n )|ψ n (3)
where E n and Γ n are, respectively, the energy and the width of the n-th resonance. Due to TRS present, H eff is a complex symmetric matrix; hence, the left and right eigenvectors are related by the transpose, ψn | = (|ψ n ) T [33].

The coupling to continuum, as described by the imaginary part of H eff , turns real eigenfunctions of the closed system into complex quasimodes of its open counterpart. In order to measure their complexness, we define the complexness parameter q 2 n as follows:

q 2 n = i (Im ψ i n ) 2 i (Re ψ i n ) 2 (4)
where ψ i n is the i-th component of the eigenvector (we note that the complexness parameter can be equivalently defined by means of the left eigenvectors). It is worth noting here that in contrast to the related Petermann factor [19], which is defined for a fixed value of the given resonance width, no additional constraints are imposed on (4). In chaotic systems, q 2 n reveals strong mode-to-mode fluctuations, which we describe through its probability distribution function to be derived below.

B. Statistical assumptions

Within the RMT approach, the universal statistical properties of closed chaotic systems with preserved TRS are described by GOE [START_REF] Stöckmann | Quantum Chaos: an Introduction[END_REF]. In this ensemble the joint probability distribution, P ({E i }), of the levels (the eigenvalues of H) is induced by a Gaussian distribution of the random real symmetric H with zero mean. The exact expression for P ({E i }) is well known to have the following form:

P ({E i }) ∝ n>m |E n -E m | exp - N π 2 8 n E 2 n . (5)
Here, we have chosen the variance of H such that it yields the mean level spacing ∆ = 1/N at the spectrum center, E = 0. The energy levels, as defined by Eq. ( 5), exhibit a linear level repulsion. As a result, the energy spectrum displays spectral rigidity which restrains the spectral fluctuations around the mean. This important feature can approximately be taken into account within the so-called picket-fence model of equidistantly spaced levels [34]. The usefulness of this model is in its simplicity that allows one to treat various resonance phenomena analytically, see, e.g., Refs. [35,36,37]. Here, we employ this model to single out a contribution to q 2 n due to fluctuations of the resonance widths.

As concerns the coupling amplitudes, the results are known to be model independent on statistical assumptions on V c n as long as the number of open channels is small compared to that of the levels [38,39]. The coupling amplitudes may be equivalently chosen as fixed [14] or random [33]. In order to preserve orthogonal invariance of H eff under (complex) orthogonal transformations [33], we consider the V c n 's as real Gaussian random variables with zero mean and

V c n V c ′ n ′ = (2κ∆/π)δ nn ′ δ cc ′ ≡ σ 2 δ nn ′ δ cc ′ . (6)
Henceforth, • • • stands for the statistical averaging over the ensembles. The coupling constant κ determines a transmission coefficient T = 1-| S | 2 = 4κ/(1+κ) 2 of the channels (assumed to be statistically equivalent). The cases of T ≪ 1 or T = 1 correspond, respectively, to weak or perfect coupling. In the weak coupling regime considered below, κ ≪ 1, all the resonances are almost isolated and Γ ≪ ∆.

III. PERTURBATIVE APPROACH A. Complexness parameter in the weak coupling regime

We now derive an expression for the complexness parameter of the eigenvectors for weakly overlapping resonances. The matrix representation of H eff in an arbitrary basis {|n } of the Hilbert space spanned by eigenvectors of H reads:

H eff = N n,p=1 |n H np p| - i 2 N n,p=1 M c=1 |n V c n V c p p| (7)
As we focus on the weak coupling regime, the imaginary part may be viewed as a perturbation of the Hamiltonian of the closed system. The repulsion of the energy levels exhibited by the systems under investigation allows us to consider the eigenenergies of H as nondegenerate. One can therefore apply first-order perturbation theory to obtain from (3) the eigenvalues and the eigenvectors of H eff straightforwardly. The eigenvalues read E n -i 2 Γ n , where the E n 's are the eigenvalues of H and the widths Γ n are given by:

Γ n = M c=1 (V c n ) 2 . ( 8
)
The perturbed eigenvectors of H eff written in the eigenbasis {|φ n } of H are easily found as follows:

|ψ n = |φ n -i p =n φ p |V V T |φ n 2(E n -E p ) |φ p . (9)
Splitting then the real and imaginary parts of |ψ n , the complexness parameter (4) of a given eigenvector reads

q 2 n = p =n Γ 2 np 4(E n -E p ) 2 , (10)
where we have introduced

Γ np = M c=1 V c n V c p .
These quantities are responsible for the coupling and interference of the resonance states due to the common decay channels [33].

In what follows, we study the statistical properties of the complexness parameter (10) for H being described by a picket-fence or belonging to GOE. It is worth noting here that expression (10) is a sum of correlated random variables which, therefore, does not obey the standard central limit theorem. Statistics of a similar kind of objects appears, e.g., in the study of the parametric level dynamics ("curvature") [40] and in the context of interference effects in neutron scattering from compound nucleus [41].

B. Rescaled parameters and their statistics

The complexness factor (10) contains two contributions of distinct types, one is due to the internal levels and the other is due to the coupling matrix elements Γ np . From a statistical point of view, these two are statistically independent of one another. We note, however, that the levels E n are mutually correlated. The quantities Γ np 's, unlike the original amplitudes V c n , are also not statistically independent. Although their joint distribution can be found from ( 6), the resulting expression is quite complicated [42], being of little practical use for actual calculations in the present context.

To overcome the difficulty of averaging over the coupling amplitudes, we follow Sokolov and Zelevinsky [33] and treat an arbitrary matrix element Γ np as a scalar products between M -dimensional vectors V n and V p of the coupling amplitudes {V c l } associated with the levels l = n and l = p. This suggests a natural parametrization for Γ np in terms of the angles θ np between the pairs of these N vectors,

Γ np = (V n • V p ) = Γ n Γ p cos θ np . (11)
The main advantage of this representation is that the angles θ np are mutually independent and also independent of Γ n .

The probability distribution of any angle (for M ≥ 2) can be easily found to be given by the expression for a solid angle in an M -dimensional space [33]:

p M (θ) = Γ(M/2) √ πΓ((M -1)/2) sin M-2 θ. ( 12
)
Note that Γ np = Γ n Γ p at M = 1. As concerns the lengths of these vectors, i.e., the widths (8), these are well-known to be independent and χ 2 distributed according to

P M (γ) = 1 2 M/2 Γ(M/2) γ M/2-1 e -γ/2 . ( 13
)
Henceforth γ n = Γ n /σ 2 stands for the dimensionless widths. This distribution function has the mean value γ = M and the variance

var(γ) = 2M = 2 M γ 2 . ( 14
)
Thus the widths cease to fluctuate as the number of open channels grows, with the average width being kept fixed.

It is now convenient to express all the quantities in their natural units and to consider a rescaled complexness parameter X n defined as follows:

X n ≡ ∆ 2 σ 4 q 2 n = γ n p =n ∆ 2 Z p 4(E n -E p ) 2 , (15)
where we have introduced the following quantities

Z p = γ p cos 2 θ np . (16)
Z p may be given a geometrical interpretation as (a square of) the projection of the vector σ -1 V p along the direction given by the vector V n . These projections are statistically independent, as is obvious from the above discussion. The probability distribution of any projection follows readily from Eqs. ( 12) and (13). Performing an integration first over γ and then over θ in the definition P (Z) = δ(Z -γ cos 2 θ) , one finds

P (Z) = 1 √ 2πZ e -Z/2 . ( 17
)
Thus, surprisingly, the distribution of Z p is independent of M , being given by the Porter-Thomas law at any M ≥ 1.

C. Average of X and width fluctuations

A general expression of the average value of the complexness parameter X can be readily found from Eq. ( 15) by making use of the mutual statistical independence between the widths {γ n }, the projections {Z n } and the levels {E n }. Noting that γ = M and Z = 1 = γ cos 2 θ , one obtains

X = M f , (18)
where the factor f depends on the statistical properties of the energies of the closed system only,

f = p =n ∆ 2 4(E n -E p ) 2 . ( 19
)
It is important to note that, generally, the nonzero values of the complexness parameter are solely due to fluctuations of the resonance widths. Indeed, in the extreme case of all the widths being equal, the anti-Hermitian part of H eff gets proportional to the unit matrix and, as a result, the complex (biorthogonal) eigenvectors become essentially real [23]. It is, therefore, instructive to take this explicitly into account and, in view of relation ( 14), bring Eq. ( 18) to the form:

X = f 2 var(γ) . ( 20
)
This expression relates the average complexness parameter to the natural measure of the width fluctuations, its variance.

Strong correlations between the complexness parameter and the spectral widths are already known. The proportionality between X and the average value of the fluctuating part of damping was recently found experimentally in a chaotic microwave billiard at room temperature, where this was also explained heuristically using a ray picture based on the ergodic character of the wave system [5]. Then this proportionality was established in [23] using a two-level RMT model and considering M ≫ 1 that was relevant for this experiment. Expression (20) readily provides this feature, in view of X = γ f /M , at any N and M . On the other side, it captures fluctuation properties of the widths properly, e.g., yielding the vanishing X in the absorptive limit of many weakly coupled channels with the average total width kept fixed, due to the vanishing variance (14). Therefore, we believe that relation ( 20) is a general feature of weakly open chaotic systems with non-degenerate spectrum in the perturbative regime. Figure 1 supports this suggestion through numerical simulations of the picket-fence and GOE models (with the details being given later in the next section).

A remark on the proportionality factor f is appropriate here. In the RMT limit N → ∞, this factor may be represented as follows f = 1 2 ∞ 0 ds s -2 R 2 (s), where R 2 (s) is the twopoint correlation function of the RMT. The main problem of the GOE case, already mentioned in [19,23], is an 'infrared' logarithmical divergency of f due to R 2 (s) ∼ s at s → 0. Practically, this divergence can be regularized by introducing a cut-off at small s, s ≥ ǫ, see Fig. 1. Without this cut-off the expression of the complexness parameter obtained using first 20) is represented by the solid line. In the GOE case (a), the proportionality factor is given by the regularized expression fǫ = 1 2 R ∞ ǫ ds s -2 R2(s). The results obtained with three different values of the cut-off ǫ are shown. In the picket-fence case (b), f = π 2 /12 as exactly given by Eq. (23).

order perturbation theory (15) does not yield finite moments, thus demanding for the characterization of fluctuations of X by means of its probability distribution.

IV. DISTRIBUTION FUNCTION

The probability distribution function of the rescaled complexness parameter X n is defined as follows

P M (X) = δ(X -X n ) , (21)
where the statistical averaging over the levels, the widths and the projections is performed with the help of Eqs. ( 5), (13), and (17), respectively. In the weak coupling regime, function (21) depends only on the number M of open channels. It is instructive first to consider the case of the completely rigid spectrum, which may be viewed as an approximation of the GOE spectrum where the fluctuations are neglected.

A. The picket-fence model

In this model the eigenenergies of the closed system are equally spaced, i.e. E n -E n±k = ±k∆, and the eigenvector components are random Gaussian variables. The complexness parameter is then given by

X n = γ n k =0 Z k 4k 2 . ( 22
)
This expression does not have any divergence problems of the GOE case, thus statistics of ( 22) can be also characterized by its moments. In particular, the average value is easily found to be exactly given by Eq. ( 20), with the factor f being

f = k =0 1 4k 2 = π 2 12 . (23)
Figure 1(b) illustrates the dependence X = π 2 24 var(γ). We now derive an exact expression for the probability distribution P pf M (X) in the picket-fence case. First we substitute in the definition (21) the Fourier representation of the delta function, δ(X -X n ) = dω 2π e iω(X-Xn) , where X n is given by Eq. ( 22). Then the integration over the projections Z k with the help of Eq. ( 17) becomes trivial, yielding

P pf M (X) = ∞ -∞ dω 2π e iωX ∞ 0 dγ P M (γ) ∞ k=1 1 1 + i ωγ 2k 2 . (24)
The infinite product here can be evaluated explicitly [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]. Making use of the explicit expression (13) for P M (γ) and applying the change of variables γ = 2|z| 2 , Eq. ( 24) can then be cast in the following form:

P pf M (X) = 1 Γ(M/2) +∞ -∞ dω 2π × +∞ -∞ dz |z| M-1 e iωX-z 2 √ iωzπ sinh( √ iωzπ) . (25)
As one can easily check, this expression is properly normalized to unity. It is also worth noting that the integrand of Eq. ( 25) is an analytic function in ω except for the poles located on the upper part of the imaginary axis at

ω k = i(k/z) 2 , k = 1, 2, . . . , ∞. This readily implies that P pf M (X) = 0 at X < 0 identically.
The details of the subsequent calculations of P pf M (X) are given in Appendix A. The final expression reads: with J ν (x) being the Bessel function of order ν. In the case of an odd number of channels, M = 2n + 1, n = 0, 1, . . . , this expression can be integrated further to yield an attractively simple formula

P pf M (X) = 2π( √ X) M/2-1 Γ(M/2) × ∞ 0 dz z M/2+1 sinh(zπ) J M/2-1 (2 √ Xz) , (26
P pf 2n+1 (X) = √ πX n-1/2 2Γ(n + 1/2) - ∂ ∂X n 1 cosh 2 ( √ X) . (27)
In particular, the single-channel distribution P pf 1 (X) reads

P pf 1 (X) = 1 2 √ X 1 cosh 2 ( √ X) . (28)
It is interesting now to study in details the case of the large number of weakly open channels M ≫ 1. In view of the scaling (18), we consider the limiting probability distribution of x = X/M defined as

p(x) = lim M→∞ M P M (M x) . (29)
Expression ( 26) is actually not very convenient for evaluating this function. However, one can note that in the limit considered, the distribution P M (γ), Eq. ( 13), tends to the Dirac distribution, δ(γ -M ). Then, starting from Eq. ( 24), the integration over γ is trivial and the probability distribution of x reads:

p pf (x) = 1 2π ∞ -∞ dω e iωx ∞ k=1 1 1 + iω/(2k 2 ) (30)
Using the residue theorem, one readily gets:

p pf (x) = 4 ∞ k=1 (-1) k+1 k 2 e -2k 2 x (31)
and finally

p pf (x) = -2e -2x d dx ϑ 4 (0, e -2x ) (32)
where ϑ 4 is a Jacobi theta function [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF].

The above analytical predictions concerning the average value of the complexness factor and its probability distribution have been checked through numerical simulations of random matrices, see Figs. 1 and2. Numerical simulations are based on the diagonalization of the effective Hamiltonian (2) viewed as a random non-Hermitian matrix. We have considered resonances in the bulk only, i.e. resonances with a large number of neighbors on the left and on the right of the spectrum. This restriction is introduced to neglect the edge effects whose contribution tends to vanish as N → ∞.

The picket-fence Hamiltonian is built such that the eigenenergies are equally spaced and the eigenvectors are random Gaussian variables. This is readily done by following a procedure adapted from [START_REF] Gorin | [END_REF] where the authors used it to generate the POE ensemble. Thus, in a basis deduced from its eigenbasis through an arbitrary orthogonal transformation O with random Gaussian variables, the Hamiltonian H is given by:

H = O diag{E n }O t (33)
where E n = n/N , such that ∆ = 1/N , and

O ij = 0, O 2 ij = 1/N (34)
Statistics were performed with 100 matrices of size 1000×1000. In order to make the calculated distributions insensitive to edge effects, 100 levels at each end of the spectrum were discarded. In all the simulations the mean spectral width is kept fixed and equal to Γ /∆ = 10

B. The GOE model

The probability distribution in the GOE case can be found by making use of group integration methods and results obtained in [19]. Outlining the details of the computation in Appendix B, we state the final result here:

P goe M (X) = π 2 M 24X 2 1 + π 2 (3 + M )/(4X) [1 + π 2 /(4X)] M/2+2 . (35)
To check our findings, the same kind of numerical simulations as in the picket-fence model have been performed. The closed Hamiltonian H now belongs to GOE, its elements being defined by their first two moments:

H ij = 0, H 2 ij = 4/(N π 2 ), i = j 2/(N π 2 ), i = j, (36)
where N is the size of the matrix. Like in the picket-fence case, the normalization is chosen such that ∆ = 1/N . Statistics were obtained with 150 matrices of size 1000×1000.

Only levels near E = 0 for which spacings deviate less than 5% from ∆ were kept. The agreement between numerical and analytical results is flawless, as shown in Fig. 3. The comparison between the probability distribution of X in the picket-fence model and for GOE illustrates the effects of the fluctuations of the spectrum on the complexness parameter. The maximum of both distributions are close to each other. This is mainly due to the spectrum rigidity in both ensembles. But at large X the statistical weight is larger for GOE than for the picket-fence model. This difference is introduced by the behavior of the levels at small distance: the spacing of two eigenenergies can be very small, the corresponding contribution to the complexness parameter is large, then the tail of P M (X) is larger for GOE than for the picket-fence model. This feature is most explicitly seen by comparing the corresponding limiting distributions at M ≫ 1. The distribution (29) of x is easily obtained from Eq. ( 35) and reads

p goe (x) = π 2 24x 2 1 + π 2 4x exp - π 2 8x (37)
In contrast to the asymptotic exponential behavior in the picket-fence case, p pf (x) ∝ e -2x , see Eq. ( 31), the tail of the distribution (37) follows a power-law decay: p goe (x) ∝ x -2 .

V. CONCLUSION

In this paper, we have studied the statistics of complex wavefunctions associated to the resonances of weakly opened wave chaotic systems with the preserved time-reversal symmetry. More specifically, in the perturbative regime, we have considered the case of the completely rigid spectra defined through the picket-fence model and that of the GOE displaying spectral fluctuations. One of the key features of this study relies on the proportionality between the average of the complexness parameter and the variance of the resonance widths, which we believe is valid for generic nondegenerate spectra. We have also derived the exact probability distribution of the complexness parameter in these two cases.

To check the validity of the present results, recent experiments in elastodynamics are available. In particular, in the case of vibrating plates, a complete knowledge of the eigenfunctions can be obtained through noninvasive measurements [45] even for moderate overlap of resonances. Indeed, the understanding of the statistics of eigenfunctions beyond the perturbative regime still remains an open problem. (We note that some relevant interesting numerical results for microwave billiards with large openings were recently reported [46].) Finally, one should also note that the complexness parameter may be considered as a sensitive probe of the crossover from localized to extended states in open disordered systems [47].

25 FIG. 1 :

 251 FIG.1:The average rescaled complexness parameter versus the width variance for the GOE and picket-fence models. The symbols correspond to the results of numerical simulations performed at M = 1, 2, 3, 5, 10 (see the text for details). The linear dependence predicted by Eq. (20) is represented by the solid line. In the GOE case (a), the proportionality factor is given by the regularized expression fǫ = 1

3 FIG. 2 :

 32 FIG. 2: (Color online)The distribution of the rescaled complexness parameter for the picket-fence model at M = 1, 5 and 10 (top, middle and bottom, respectively). The analytical result(25) is plotted in the solid line while the histograms correspond to numerics. Insets show the tail of the distribution in a semi-log scale.

FIG. 3 :

 3 FIG. 3: (Color online)The distribution of the rescaled complexness parameter for the GOE model at M = 1, 5 and 10 (top, middle and bottom, respectively). The analytical result(35) is shown in the solid line and compared to that (25) of the picket-fence case (dashed line) while the histograms correspond to numerical simulations. Insets show the tail of the distribution in a log-log scale.

Acknowledgments

We are grateful to H.-J. Sommers for his instructive advice and help with evaluating Eq. ( 26). One of us (D.V.S.) acknowledges gratefully the generous hospitality of LPMC in Nice and the financial support of University of Nice during his stay there. The partial financial support by BRIEF grant (D.V.S.) is also acknowledged.

APPENDIX A: DERIVATION OF EQS. (26) AND (27)

We first note that the integrand of Eq. ( 25) is a symmetric function in z that allows us to restrict the z-integration to the positive axis. Then we deform the contour of integration over ω from the real to imaginary axis by putting Ω = iω. Performing after that the scaling transformations of the integration variables, first z → z/ √ Ω and then Ω → Ω/X, and interchanging the order of integrations over z and Ω, we may cast Eq. ( 25) in the following form

To calculate here the last integral over Ω, we expand e -Xz 2 /Ω into a series and evaluate the result termwise

where we have used

. Making now use of the well-known series representation for the Bessel function [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF], one can immediately recognize the r.h.s. of (A2) to be equal to

Collecting all the factors together, we finally arrive at Eq. (26).

Further progress is possible in the case of odd M . It is instructive first to start with the case of M = 1, which turns out to play the central rôle in this calculation. We may use the known relation J -1/2 (z) = 2z/π cos(z)/z in this case [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF], thus

allows us to perform the integration in Eq. (A1) analytically:

Taking now into account the (omitted) factor 2 π/X, we obtain P pf 1 (X), Eq. ( 28). The general case of odd M = 2n + 1 may be reduced to that of M = 1 considered above, if one notices that the term z M /Ω M/2 e -Xz 2 /Ω in the integrand of Eq. (A1) can be generated by a differentiation with respect to X as follows:

(A4) Substituting this representation into Eq. (A1) and changing the order of the integrations and differentiation there, we see that the resulting integral is already given by Eq. (A3) that readily yields the expression (27) of Sec. IV A.

APPENDIX B: DERIVATION OF EQ. (35)

We use the recent result by Schomerus et al. [19], who calculated the joint probability distribution P (A, B) of

where {α p } are the statistically independent real Gaussian variables distributed according to

and {E n } are taken from the GOE. They found the following expression for P (A, B):

We note that the above expression (B3) was obtained in [19] for the particular case of one open channel. The key fact which allows us to apply this result to our M -channel case is the representation (15) in terms of projections with the distribution (17). The later corresponds to the Gaussian distribution (B2) with κ = 1 and π 2 κ∆ α 2 p = Z p , thus giving a connection X = π 2 4 γB. Correspondingly, the distribution function of X in the GOE case can be found from

by averaging over A, B and γ. Substituting the explicit form (B3), it is convenient first to integrate out B that yields

dA(1 + π 2 A 2 )e -(aγ/2)π 2 A 2 , (B5) with a = π 2 /4X. The Gaussian integration over A is now straightforward and gives

where we have substituted expression (13) for P M (γ). The remaining integration yields Eq. ( 35).


Daniel Garcia-Cortes

Catherine Xuereb

Patricia Taillandier

Ulises J Jauregui-Haza

Joel Bertrand

Flow Induced by Dual-Turbine of Different Diameters in a Gas-Liquid Agitation System: the Agitation and Turbulence Indices

Keywords: dual-impeller, disk style turbine, agitation index, turbulence index, LDV, gas-liquid

Flow induced by a dual turbine stirred tank was characterized measuring local velocities with a LDV and drawing the main velocity fields and the maps of turbulence intensities. The hydrodynamic regime studied in all the experiments was the so-called merging flow regime. Two impeller configurations were studied. In the first one, two disk style turbine of the same dimensions (configuration A) were used, while in the second one, the dimensions of the upper turbine were 20 % proportionally smaller than those of the lower turbine (configuration B). The agitation and turbulence indices were used to evaluate, as a first order approximation, the power consumption distribution between convective and turbulent flows. The comparison of the two-phase agitation systems studied showed that configuration B seems to be more efficient than configuration A, since both induce a similar global convective flow, but the first one assures a significant reduction of power consumption. The distribution of power consumption between convective and turbulent flows was evaluated using the agitation index and a new global parameter: turbulence index.

1.

INTRODUCTION

The gas -liquid dispersion process in stirred tanks with standard geometry, equipped with a ring sparger and dual-turbine agitation system, where the liquid height is larger than the tank diameter (H > T) has been widely reported [START_REF] Abrardi | Sparged vessels agitated by multiple turbines[END_REF][START_REF] Hudcova | Gas liquid dispersion with dual Rushton turbine impellers[END_REF][START_REF] Roman | Hydrodynamics of a gas-liquid vessel stirred with Rushton turbines in newtonian and non-newtonian systems[END_REF][START_REF] Vasconcelos | Scale-up of mixing in gassed multi-turbine agitated vessels[END_REF]. However, in the literature is scare the experimental data concerning aeration-agitation system with all the following characteristics: a) H = T [START_REF] Kuboi | The power drawn by dual impeller systems under gassed and ungassed conditions[END_REF], b) baffles fixed away from the wall and the bottom, c) tanks with curved bottom [START_REF] Bouaifi | Bubble size and mass transfer coefficients in dual-impeller agitated reactor[END_REF], and d) dual-turbine impellers. This geometrical configuration is common in laboratory and pilot scale bioreactors [START_REF] Garcia-Cortes | Effect of dual impeller-sparger geometry on the hydrodynamics and mass transfer in stirred vessels[END_REF]. In our knowledge, configuration system with two disk turbines with the upper impeller smaller than the lower one has not been before studied.

The flows generated in stirred tanks with H = T provided with dual -impeller agitation systems of D = T/3 were studied by [START_REF] Rutherford | Hydrodynamics characteristics of dual Rushton impeller stirred vessels[END_REF]. They identified three stable flow patterns: parallel, merging and diverging. Studying these flow patterns, [START_REF] Mahmoudi | Velocity and mixing characteristics of stirred vessels with two impellers[END_REF] determined that the merging flow pattern led to a mixing time lower by around 20 % compared to the two others patterns. This make the merging flow pattern more attractive for mixing operations. Taking into account this fact, the impeller configuration selected in this work guarantees the presence of the merging flow pattern in the stirred tank.

In the aeration -agitation systems with dual impellers the power drawn by the lower impeller is smaller than the power drawn by the upper one because of the direct influence of the gas flow from the sparger over the lower impeller [START_REF] Kuboi | The power drawn by dual impeller systems under gassed and ungassed conditions[END_REF]. One way to achieve similar power consumption in both impellers is reducing the dimensions of the upper impeller. The aims of this work were to establish the liquid flow patterns in the two -phase system with the upper impeller smaller than the lower one and to analyse the overall energy distribution between convective and turbulence flows by means of calculation of the agitation index and a new global parameter called turbulence index.

EXPERIMENTAL

Figure 1. Aeration -agitation system.

The experiments were performed in a curved bottomed cylindrical tank with internal diameter of T = 0.19 m and a radii of bottom curvature of 0.19 m (Figure 1). The mean diameter of the ring sparger was 0.042 m. The sparger was located with a bottom clearance of 0.13 T. The reactor was filled with water until H = T and sparged with air at 8.7 10 -5 m 3 /s (1 vvm). The diameters of the disk style turbines used in these experiments were 0.065 and 0.052 m. In the studied agitation systems the turbines were positioned with a similar bottom and impeller clearance equal to 1/3 T. The vessel was placed inside a transparent square tank, filled with water, in order to minimize refraction at the cylindrical surface during the measurements of local velocities in a grid of 5 x 11. The measurements were performed in the vertical plane at 45° between the baffles. The radial and axial velocities were measured with a two-beam Dantec Fiberflow system operating in a backscattering mode. Seeding was found to be necessary; for this purpose, small amount of Iriodin 111 (Merck) particles was added to the contents of the vessel.

The torque on the agitator shaft was also measured by a torque-measuring coupler (Staiger and Mohilo; precision 0.2 %) mounted on the agitator shaft [START_REF] Baudou | Agitation par de systemes axiaux simples ou multi-etages. Obtention de l'hydrodynamique par velocimetrie laser à effet Doppler[END_REF]. This experimental value of torque was corrected for friction losses using measurements performing in an empty vessel and was then used to calculate the power consumption. The experiments were conducted at a Reynolds number (Re) of 52 500, corresponding to an impeller rotational speed of 500 ± 1 rpm. More details about the equipment and the methodology of measurements can be found elsewhere [START_REF] Mavros | Determination of 3-D flow fields in agitated vessels by laser-Doppler velocimetry: Effect of impeller type and liquid viscosity on liquid flow patterns[END_REF].

RESULTS AND DISCUSSION

Mean Velocity Vector Fields and Turbulence Maps

The map of velocity vectors obtained from the velocity measurements in the single liquid phase system with two disk style turbines of the same dimensions (configuration A) is shown in Figure 2a. The two characteristic ring vortices of this kind of systems are observed here. However, the results obtained differ from that reported in [START_REF] Kuboi | The power drawn by dual impeller systems under gassed and ungassed conditions[END_REF] and [START_REF] Hudcova | Gas liquid dispersion with dual Rushton turbine impellers[END_REF], where the impeller streams follow an almost straight-line orientation toward one another until they merge at an elevation approximately midway between the impellers. In this agitation system the stream from the lower impeller is not observed and the upper ring vortices of both impellers seem to be fused, forming one vortex with the centre in the midway between the impellers. The sense of fluid rotation of this vortex corresponds with the usual sense of the upper vortex in agitation systems provided with a turbine impeller. This situation presupposes a predominance of the lower impeller on the upper one during the interaction between both impellers. The lower vortex of the upper impeller is reduced and located near the shaft, between the impellers. This structure variation of the merging flow pattern may arise due to the specific combination of the baffles position and the curved bottom form used in this work, and it remains to be evaluated whether other configuration reveals a similar behaviour. Four focus of maximum turbulence are observed in the turbulence map (Figure 2b). This result differs from that reported by [START_REF] Rutherford | Hydrodynamics characteristics of dual Rushton impeller stirred vessels[END_REF], who found two focuses of turbulences near the impeller tips in the region between the impellers. In general there is a more organized flow in the lower half of the tank and more turbulent near the upper impeller tip and in the vicinity of the impeller shaft.

Characteristic mean velocity vectors in the same agitation system, but in a two-phase gas-liquid system, are shown in Figure 3a The liquid circulation vortex situated between the impellers becomes more remarkable, presumably helped by the bubbles that ascend next to the shaft from one impeller to another. In Figure 3b The LDV measurements in the same previous two-phase aeration-agitation system, where the upper impeller dimensions were diminished in 20 % with regard to the lower impeller (configuration B), are shown in Figure 4. This modification caused a drastic change in the flow pattern, in such a way that it is not possible to identify any vortex in the vertical plane at 45° between baffles (Figure 4a). In the turbulence map of this system (Figure 4b), there are four maxima of turbulence intensity which are located by order of importance: under the sparger, between the impellers near the shaft, near the lower part of the upper impeller tip, and near the bottom of the vessel. In comparison with the aeration-agitation system with both impellers of the same dimensions (Figure 3b), the turbulence maximum in this agitation system (Figure 4b) is higher, although the other turbulence maxima are weaker.

Quality of Agitation

Turbulence flow has been defined by different authors [START_REF] Hinze | Turbulence[END_REF][START_REF] Davies | Turbulence phenomena[END_REF][START_REF] Tennekes | A first course in turbulence[END_REF][START_REF] Bradshaw | Introduction to turbulence and its measurements[END_REF][START_REF] Lesieur | Turbulence in fluids[END_REF]. According to Davies's definition: "as the flow rate is increased, the flow becomes unsteady, with chaotic movements of parts of the liquid in different directions superimposed on the main flow of the liquid. Such movement of any particular element of fluid is very complicated, and it can only be described in terms of averages". [START_REF] Baldyga | Turbulence mixing and chemical reactions[END_REF] pointed out that the main flow may be identified with the flow described by the average velocity and that the chaotic movements of fluid elements are equivalent to the turbulent fluctuations of velocity (turbulent flow), interpreted as the local instantaneous values minus the average one. In this way, the convection flow generated by the impeller may be identified with the main flow and the turbulence flow with the fluctuating velocities in order to study the performance of impellers [START_REF] Zhou | Distribution of energy between convective and turbulent flow for three frequently used impellers[END_REF].

To assess the quality of agitation of single phase and two-phase gas -liquid systems described previously, two global parameters were used, the agitation and turbulence indices. These parameters allow quantify objectively the extent of mixing (macromixing and micromixing) achieved in the stirred vessel based on LDV velocity measurements.

The agitation index (I g ), defined by [START_REF] Mavros | Quantification of the performance of agitators in stirred vessels: Definition and use of an agitation index[END_REF], represents the global mean velocity expressed as a percentage of the tip velocity:

[ ] % U U Î tip g 100 = ( 1 )
where U tip is the agitator tip velocity and Û is the volume weighted average velocity.

The principal movements of the fluid inside the vessel are studied using most often 2-D velocity measurements, although it is well known that this flow is 3-D. Nevertheless, [START_REF] Mavros | Quantification of the performance of agitators in stirred vessels: Definition and use of an agitation index[END_REF] demonstrated that it is possible to use this reduced sets of data, if comparisons are made consistently among agitation indices obtained with 2-D velocities. In this way, the volume weighted average velocity was determined following the postulate that each composite mean local 2-D (U ij ) velocity corresponds to a 3-D cell with a liquid volume, which is related to the vessel dimensions and the grid point coordinates (Figure 5):

∑∑ ∑∑ = i j ij L i ij j L V U V U ˆij (2)
The denominator in Equation ( 2) corresponds approximately to a vessel volume minus the volume swept by the agitator (calculation methodology of and [START_REF] Baudou | Agitation par de systemes axiaux simples ou multi-etages. Obtention de l'hydrodynamique par velocimetrie laser à effet Doppler[END_REF]. The composite mean local 2-D velocity for each grid point (g ij ) was calculated using measurements of the components of mean local velocities (V ij,k ) carried out with the LDV:

V Lij published in Mavros r 1, z 2 r 1, z 1 r 2, z 1 V L ij V L ij
2 1 2 / k k , ij ij V U ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = ∑ , k = r, z ( 3 )
The turbulence index (I t ) is a new global parameter proposed by [START_REF] Garcia-Cortes | Contribution à l'étude de l'hydrodynamique et du transfert de masse en réacteurs de type cuve agitée par une approche globale et locale en système gas-liquide[END_REF] to quantify the magnitude of the turbulent flow in the agitation system. The procedure to calculate this index is similar to the one proposed for I g . First, the composite root mean square (r.m.s.) local 2-D velocities (U ij *) for each grid point (g ij ) is calculated using LDV measurements of the components of r.m.s. local velocities (u ij,k ):

2 1 2 / k k , ij ij u * U ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = ∑ , k = r, z ( 4 )
Further, a volume -weighted average r.m.s. velocity was determined following a similar postulate, that each such composite r.m.s.velocity (U ij *) corresponds to a volume of liquid (V Lij ), which is related to the vessel dimensions and the grid point coordinates:

∑∑ ∑∑ = i j ij L i j ij ij L V * U V * U ˆ(5)
where the denominator corresponds approximately to the vessel volume excluding the agitator-swept region. Finally, expressing this volume -average r.m.s. velocity as a percentage of the agitator tip velocity (U tip ) yield the turbulence index.

[ ]

% U * U Î tip t 100 = ( 6 )
The turbulence index represents the global r.m.s. velocity and gives a global measure of the turbulence intensity in the vessel bulk, contrary to the parameter turbulence intensity, which is a local parameter that represents a relative measure of the fluctuations of the local velocity.

The agitation index obtained in the single phase system with dual-turbines of the same diameters (D = T/3) (Table 1) was compared with the value of this parameter (Ig(rz) = 14.5 %) reported by [START_REF] Mavros | Quantification of the performance of agitators in stirred vessels: Definition and use of an agitation index[END_REF]. This last report was done for a single phase agitation system too, obtained in the same experimental setup, with the following geometrical configuration: four baffles positioned without clearance to the bottom and to the wall; and a single Rushton turbine with a diameter of D = T/2 located with a clearance to the bottom of T/3. The value of the agitation index reported by Mavros and Baudou for the one-impeller system is greater than the value obtained in this work for the dual-impeller system. This is presumably due to the important influence of the impeller diameter (D = T/2) on main flow that, in turn, determines that the average convective flow in an agitation system with a single impeller is greater than in the dual-impeller system.

The agitation index calculated for configuration A, but in the two-phase gas -liquid system showed a decrease of the convective flow in 21 %. This result was expected taking into account the decrease of the energy input caused by formation of the air cavities behind the turbine blades. However, the reduction of the upper impeller dimensions in 20 % with regard to the lower one, in the same two-phase system, did not cause a significant variation of the agitation index, and hence, of the convective flow. The distribution of liquid volumes related to specific mean velocity subranges, in the same way proposed by [START_REF] Mavros | Quantification of the performance of agitators in stirred vessels: Definition and use of an agitation index[END_REF], is shown in the Figure 6. In the single liquid phase system (Figure 6a), near 60 % of the liquid volume, excluding the agitator-swept region, flows with an average velocity lower than 0.1 U tip , while in this agitation system with an aeration of 1 vvm, this volume fraction increases until 70 %. However, in the same two -phase system when the dimensions of the upper impeller were reduced in 20 % with regard to the lower one, the liquid volume fraction that flows with an average velocity lower than 0.1 U tip is also of ~ 70 %, while the distribution profile resembles more to the distribution profile of the single phase system with dual impellers of the same dimensions.

The turbulence index was calculated in the studied agitation systems by means of Equations 4 -6 (Table 1). The inclusion of a second phase in the agitation system with configuration A caused a slight decrease of the turbulent flow. However, the reduction of the dimensions of the upper impeller in 20 % with regard to the lower one in the two-phase gas -liquid system produce a diminution of the turbulent flow in 18 %.

The distribution of liquid volumes related to specific r.m.s. velocity subranges is shown in the Figure 7. In agitation systems, while more narrow is this distribution, more homogeneous is the turbulence in the vessel bulk. In this way, comparing the graphics corresponding to the three studied systems, it can be concluded that the turbulence is more homogeneous in the two-phase systems with impellers of different dimensions.

3.3

Energy distribution between convective and turbulent flows [START_REF] Zhou | Distribution of energy between convective and turbulent flow for three frequently used impellers[END_REF] determined the way in which the impellers distribute energy between convective and turbulent flows using the macroscopic energy balance equation applied to the control volumes. Nevertheless, they remarked that efforts to extend this analysis to regions of the tank far removed from the impeller were not successful.

As a first order approximation of energy distribution in the bulk volume of the tank between convective and turbulent flows, we propose to use the agitation and turbulence indices that are global parameters. For example, as was analysed previously, the decrease of the energy input caused by formation of air cavities behind the turbine blades is reflected by the agitation index and more slightly by turbulence index (for the turbulence the reduction of power consumption caused by aeration is partly compensated by the turbulent kinetic energy that the gas expansion transfers to the liquid).

However, the decrease of power consumption due to reduction of the upper impeller dimensions, in the same two-phase system, no caused a significant diminution of the agitation index, and hence, of the global convective flow. In this case the occurrence of the flow pattern variation (Figure 3 and4) led to a variation of the power consumption distribution between convective and turbulent flow where evidently the decrease of power consumption has been at the expense of the decrease of the global turbulent flow reflected in the diminution of the turbulence index.

Finally, from the comparison of the two gas-liquid agitation systems studied (Table 1), it can be concluded that the system with two disk style turbine with six square blades where the upper impeller is smaller than the lower one is more appropriate for processes where it is important the fluid circulation in the vessel, since it induces an average convective flow in the bulk of the reactor similar to that induced by the agitation system with both impellers of the same dimensions, but with a smaller energy consumption (~36 %).

CONCLUSIONS

The flow patterns induced by dual -impeller agitation systems in a non-standard vessel have been established from laser Doppler velocimetry in a liquid and gas -liquid systems. In addition, a dual -turbine system with the upper impeller smaller than the lower one was characterized in a gas -liquid system. This hydrodynamic study allowed to draw the fields of velocity vectors and the maps of turbulence intensity in the vessel bulk of the three studied agitation systems, as well as to evaluate the distribution of the power consumption between convective and turbulent flows by means of two global parameters, the agitation and turbulence indices.

The comparison of the two-phase agitation systems studied showed that the system with two disk turbines with the upper impeller smaller than the lower one seems to be more efficient than the agitation system with the impellers of the same dimensions, since both induce a similar global convective flow, but the first one assures a significant reduction of power consumption.

NOTATION

Figure 2 .

 2 Figure 2. Measurements in the plane r -z, midway between the baffles in the single liquid phase agitation system with dual -impellers of the same dimensions: a) 2-D map of mean velocity vectors, b) Contours of three component fluctuating velocities.

Figure 3 .

 3 Figure 3. Measurements in the plane r -z, midway between the baffles in the two -phase gas -liquid agitation system with dual -impellers of the same dimensions: a) 2-D map of mean velocity vectors, b) Contours of three component fluctuating velocities.

Figure 4 .

 4 Figure 4. Measurements in the plane r -z, midway between the baffles in the two -phase gas -liquid agitation system with dual -impellers of different dimensions: a) 2-D map of mean velocity vectors, b) Contours of three component fluctuating velocities.

Figure 5 .

 5 Figure 5. Illustration of part of the volume grid of the stirred tank. An example of 3-D cell of liquid volume.

Figure 6 .

 6 Figure 6. Distribution of mean velocity-related vessel volumes a) single liquid phase system with configuration A; b)two -phase gas -liquid system with configuration A; c) two -phase gas -liquid system with configuration B.

Figure 7 .

 7 Figure7. Distribution of r.m.s. velocity-related vessel volumes a) single liquid phase system with configuration A; b) two -phase gas -liquid system with configuration A; c) two -phase gas -liquid system with configuration B.

Table 1 .

 1 Flow distribution between convective and turbulent flows.

	Systems	Dimensions of the two impellers	I g , %	I t, %	P, W
		(configuration)
	Single phase	Equal (configuration A)	10.6	15.1	4.2
	Two-phase	Equal (configuration A)	8.4	14.4	2.8
	Two-phase	Different (configuration B)	8.8	11.8	1.8

ACKNOWLEDGEMENTS

DGC would like to acknowledge ALFA program project for the award of the fellowship that made this work possible.