Md Delwar

Hossain Bhuyan
email: md-delwar-hossain.bhuyan@inria.fr

Michael Döhler
email: michael.doehler@inria.fr

Laurent Mevel
email: laurent.mevel@inria.fr

Statistical Damage Localization with Stochastic Load Vectors Using Multiple Mode Sets

Keywords: Damage localization, Load vectors, Covariance analysis, Statistical analysis, Subspace method

The Stochastic Dynamic Damage Locating Vector (SDDLV) method is an output-only damage localization method based on both a Finite Element (FE) model of the structure and modal parameters estimated from output-only measurements in the damage and reference states of the system. A vector is obtained in the null space of the changes in the transfer matrix computed in both states and then applied as a load vector to the model. The damage localization is related to this stress where it is close to zero. In previous works an important theoretical limitation was that the number of modes used in the computation of the transfer function could not be higher than the number of sensors located on the structure. It would be nonetheless desirable not to discard information from the identification procedure. In this paper, the SDDLV method has been extended with a joint statistical approach for multiple mode sets, overcoming this restriction on the number of modes. The new approach is validated in a numerical application, where the outcomes for multiple mode sets are compared with a single mode set. From these results, it can be seen that the success rate of finding the correct damage localization is increased when using multiple mode sets instead of a single mode set.

INTRODUCTION

Vibration based damage localization has become an important issue for Structural Health Monitoring (SHM) such as bridges, buildings and offshore structures. Sensors installed on the structure collect data and then the modal parameters (damping ratios, natural frequencies and mode shapes) can be estimated. Those parameters are meaningful for the monitoring of the structure and damage localization is possible when updating these changes in a Finite Element (FE) model of the structure.

The SDDLV approach [START_REF] Bernal | Load vectors for damage location in systems identified from operational loads[END_REF] is a vibration based damage localization technique using both finite element information and modal parameters estimated from output data. The estimates of the modal parameters are subject to variance errors [START_REF] Pintelon | Uncertainty calculation in (operational) modal analysis[END_REF][START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Carden | Challenges in developing confidence intervals on modal parameters estimated for large civil infrastructure with stochastic subspace identification[END_REF]. Based on that uncertainty information, a statistical extension of the SDDLV method was developed in [START_REF] Döhler | Statistical decision making for damage localization with stochastic load vectors[END_REF][START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF] for deciding if an element is damaged.

In [START_REF] Döhler | Statistical decision making for damage localization with stochastic load vectors[END_REF][START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF], the number of modes used in the computation could not be bigger than the number of sensors located on the structure. This is a restriction when there are more modes describing the structure than the available sensors. Here, the SDDLV method is developed with a joint statistical evaluation using multiple mode sets. It overcomes this limitation. It is demonstrated that the computation of stress for multiple mode sets increases the information content about the damaged or non-damaged elements of the structure. Finally, all stress values corresponding to each element are being tested for damage in an hypothesis test where the computed stresses are evaluated with their joint covariance. To derive such a test, the computation of the covariance of the resulting stress is necessary. Following [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF], the necessary covariance scheme is developed and extended for a joint statistical evaluation using multiple mode sets for the same or different Laplace variables.

Assessing the performance of the method is a requirement for showing the benefits of this approach. A proper criterion for the evaluation of the success rate is proposed based on Monte Carlo simulations.

The SDDLV approach is strongly dependent on the choice of the Laplace variable s where the transfer function is evaluated. Performance can be highly different in the classical SDDLV approach depending on the choice of the Laplace variable. Accommodating multiple s variables has been treated in [START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF]. Still, the choice of the Laplace variable is a complicated part of the procedure, even if past guidelines push for choosing this variable around the identified modes in the complex plane. That is why, in this paper, Monte Carlo simulations have been done for many s variables in the complex plane and results are evaluated on 2D grids. Performance of the proposed approach is tested against the requirement that such choice should not be critical.

This paper is formulated as follows. In Section 2, the SDDLV method is presented as a vibrationbased damage localization approach. In Section 3, the effect of removing the limiting restriction on the number of modes will be discussed. In Section 4, the statistical damage localization approach is derived using multiple mode sets. In Section 5, the new approach is applied on a numerical application to evaluate the success rate of the localization, and conclusion of the work is presented in Section 6.

DAMAGE LOCALIZATION APPROACH (SDDLV)

The Stochastic Dynamic Damage Locating Vector (SDDLV) approach is an output-only damage localization method based on interrogating changes δ G(s) in the transfer matrix G(s) of a system in both reference and damaged states [START_REF] Bernal | Load vectors for damage location in systems identified from operational loads[END_REF]. A vector is obtained in the null space of δ G(s) from system identification results using output-only measurements corresponding to both states. Then this load vector is applied to the Finite Element (FE) model of the structure for the computation of a stress field over the structure. Damage localization is related to this stress field where the computed stress is zero or close to zero in practice [START_REF] Bernal | Load vectors for damage location in systems identified from operational loads[END_REF][START_REF] Bernal | Flexibility-based damage localization from stochastic realization results[END_REF][START_REF] Bernal | Load vectors for damage localization[END_REF][START_REF] Bernal | Damage localization from the null space of changes in the transfer matrix[END_REF].

In this section, the deterministic computation of the stress field and the aggregation results is summarized, before deriving its statistical evaluation in Section 4.

Modeling of a mechanical structure

The behavior of a mechanical structure can be described by a linear time-invariant (LTI) dynamic system

M Ẍ (t) +C d Ẋ (t) + KX (t) = f (t) (1)
with M, C d , K ∈ R d×d are the mass, damping and stiffness matrices respectively, where t indicates continuous time and X ∈ R d defines displacements of the d Degrees Of Freedom (DOF) of the structure. The external force f (t) is not measurable and modeled as white noise. Let the dynamic system (1) be observed at r coordinates. Since f (t) is unmeasured, it can be substituted with a fictive force e(t) ∈ R r which is acting only in the measured coordinates and that regenerates the measured output. After replacing x = [X Ẋ ] T with (1) leads to the corresponding continuous-time state-space model

ẋ = A c x + B c e y = C c x + D c e (2)
with state vector x ∈ R n , output vector y ∈ R r , the state transition matrix A c ∈ R n×n and output matrix C c ∈ R r×n , where n = 2d is the system order and r is the number of outputs. Since input of the system is replaced by the fictive force e ∈ R r , the input influence matrix and direct transmission matrix are B c ∈ R n×r and D c ∈ R r×r respectively. However, only the system matrices A c and C c are relevant from output-only system identification and the non-identified matrices B c and D c will only be needed to obtain estimates of the transfer matrix. From Stochastic Subspace Identification (SSI) [START_REF] Van Overschee | Subspace Identification for Linear Systems: Theory, Implementation, Applications[END_REF], estimates A c and C c can be obtained from output only measurements, details are given in [START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF].

Computation of damage indicator

The transfer matrix G(s) ∈ C r×r of system (2) can be derived as

G(s) = R(s)D c , where R(s) = C c (sI -A c ) -1 C c A c C c † I 0 (3)
for the restriction of the system order to n ≤ 2r which is described in details in [START_REF] Bernal | Load vectors for damage location in systems identified from operational loads[END_REF][START_REF] Bernal | Flexibility-based damage localization from stochastic realization results[END_REF]. s is a Laplace variable in the complex plane, I is the identity matrix of size r, 0 is zero matrix, and † denotes the Moore-Penrose pseudoinverse. The difference of the transfer matrices in both damaged (variables with tilde) and healthy states is δ G(s) = G(s) -G(s). The matrices δ G(s) and δ R T (s) = RT (s) -R T (s) have the same null space [START_REF] Bernal | Load vectors for damage location in systems identified from operational loads[END_REF]. The desired load vector v(s) is obtained from the null space of the δ R T (s) from Singular Value Decomposition (SVD)

δ R T (s) = UΣV H = U 1 U 2 Σ 1 0 0 Σ 2 V 1 V 2 H , (4)
where U, Σ,V ∈ C r×r , Σ 2 ≈ 0 and H indicates conjugate transpose. Let n u be the dimension of image U 1 and (rn u ) be the dimension of the null space V 2 , where n u depends on kind and number of damaged elements [START_REF] Bernal | Load vectors for damage localization[END_REF]. Note that only output data is necessary for the computation of an estimate of v(s). To compute the stress field, load vector v(s) is applied to the FE model of the structure. This stress is obtained through a linear relation to v(s) by the matrix L model (s) ∈ C l×r from FE model of the structure [START_REF] Bernal | Load vectors for damage location in systems identified from operational loads[END_REF][START_REF] Döhler | Statistical decision making for damage localization with stochastic load vectors[END_REF] and satisfies

S(s) = L model (s)v(s). (5)
Theoretically, the stress vectors S(s) indicate damage [START_REF] Bernal | Load vectors for damage location in systems identified from operational loads[END_REF][START_REF] Bernal | Load vectors for damage localization[END_REF] where entries in S(s) close to zero indicate potentially damaged elements. However, these stresses are not exactly zero but small in practice because of modal truncation, model errors and uncertainties from measurements.

Stress aggregation for robustness

Let the s-values s x , x = 1, . . . , w be given. To minimize error, they should be chosen in the area of the identified poles of the system but not too close to them [START_REF] Bernal | Load vectors for damage location in systems identified from operational loads[END_REF][START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF]. After identification of the system matrices in both states, the computations of ( 4)-( 5) are repeated for each s x to get the respective stress vectors S(s x ). For multiple s-values, the stress aggregation is obtained for each entry j as follows [START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF] S

j = w ∑ x=1 |S j (s x )| (6)
In Section 3, the theoretical limitation on number of modes is discussed and then in Section 4, a new statistical scheme is proposed using multiple mode sets while considering uncertainties of the results at the same time.

REMOVING THE RESTRICTION ON NUMBER OF MODES

In practice, there may be more modes available from identification than number of sensors on the structure. It will be meaningful to be able to utilize this information completely from the identification procedure. In [START_REF] Bernal | Load vectors for damage location in systems identified from operational loads[END_REF], it was not possible to use all modes in that context due to theoretical restrictions, n ≤ 2r where n = 2m is the system order and r is the number of sensors. Note that m is the number of conjugated complex mode pairs identified from datasets where m has to satisfy the constraint m ≤ r.

In order to remove the restriction on the number of modes, the computation of stress from different mode sets is investigated. Note that for each mode set the current restriction needs to be satisfied. Now the computation of stress is developed at multiple mode sets by joint statistical approach which should increase the information content about damage localization on the structure.

For simplicity of the notation, consider only 2 sets of modes, while the following development can be easily generalized for arbitrary mode sets. Let M = (i, j) be two sets of modes of the system that contains m i and m j modes respectively with n i = 2m i and n j = 2m j being the respective model orders. In order to compute stress, the conditions m i ≤ r and m j ≤ r need to be satisfied. For each mode set, the system matrices are obtained using SSI performed on the same measurement in the healthy and damaged states [START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF]. The "transfer matrix" resultant R(s) in ( 4) is computed from both states and the stress vector is obtained separately for the respective load vector. For each mode set i and j, the system matrices A i c , C i c and A j c , C j c are computed and subsequently "transfer matrix" results R i (s), R j (s) are also computed. Then different load vectors are obtained as v i (s) and v j (s) from the same measurement data for the respective mode sets i and j. They are respectively applied to model for stress computation and it yields S i (s) = L model (s)v i (s) and S j (s) = L model (s)v j (s). Then both stresses S i (s) and S j (s) are used for joint stress evaluation.

There is a possibility to use same or different s-values for each of the stress. For simplicity, one s-value s i is used for mode set i and s j for mode set j and then S i (s i ) and S j (s j ) are computed. However, it can also be easily generalized for several s-values. Now, several stress values are calculated because of multiple mode sets and hence, robustness of the damage localization can be obtained by aggregating the stress results.

UNCERTAINTIES PROPAGATION AND COVARIANCE COMPUTATION

For the damage localization algorithm, estimates of the system matrices A c and C c are obtained in the damaged and undamaged states from limited data lengths using SSI [START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF]. The identification of system matrices are subject to variance errors because of unknown excitation, measurement noise and limited data length. The uncertainties in the estimates are penalizing the quality and precision of the damage localization results. For making decisions about damaged elements of the structure, these uncertainties need to be used to decide whether stress of an element is significantly zero or not. In [START_REF] Döhler | Statistical decision making for damage localization with stochastic load vectors[END_REF][START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF], uncertainty quantification of the stress vector S(s) in ( 5) was computed for a single mode set at one or several svalues. In this section, uncertainty computation of the stress vector S M (s) is derived for two mode sets M , though it can also be easily generalized to any number of mode sets.

All parameters derived in this method are estimated from measurement data through the computation of the H matrix, which is the Hankel matrix of the cov-SSI method. Details are given in [START_REF] Döhler | Statistical decision making for damage localization with stochastic load vectors[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF].

Let ΣH be the covariance of vec(H) and f be a function of H that is estimated from data. Its covariance can be approximated by

cov( f (H)) ≈ J f ,H ΣH J T f ,H . (7)
Note that sensitivity J f ,H = ∂ f (H)/∂ vec(H) in the above derivation. From perturbation theory, a first-order perturbation f is defined by operator of the function f , which gives

f = J f ,H vec( H). (8)
With this relationship, the desired sensitivity can be computed from [START_REF] Bernal | Flexibility-based damage localization from stochastic realization results[END_REF] and then it can be used for covariance computation in [START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF]. The vec(.) defines the column stacking vectorization operator. Since the Hankel matrix is dependent on datasets, the estimated covariance ΣH of the Hankel matrix H estimate can be propagated to any parameters, particularly to the modal parameters and stress estimate Ŝ(s). Note that some of the matrices are complex-valued variable in Section 4.1. To deal with uncertainties of the complex-valued matrices, define an equivalent real-valued notation for any matrix Q as follows

Q Re def = Real(Q) -Imag(Q) Imag(Q) Real(Q) , Q re def = Real(Q) Imag(Q)
Now, for multiple mode sets i and j, the covariance of the system matrices can be derived as follows

Σ A i c ,C i c def = J A i c ,H J C i c ,H Σ H J A i c ,H J C i c ,H T , Σ A j c ,C j c def = J A j c ,H J C j c ,H Σ H J A j c ,H J C j c ,H T (9)
where J A i c ,H , J C i c ,H and J A j c ,H , J C j c ,H are the sensitivity of the system matrices for mode sets i and j respectively and details are given in [START_REF] Reynders | Uncertainty bounds on modal parameters obtained from stochastic subspace identification[END_REF][START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF][START_REF] Döhler | Efficient multi-order uncertainty computation for stochastic subspace identification[END_REF] for subspace identification.

Covariance computation of the stress vector S(s)

For stress variance computation, the sensitivity J R(s),(A c ,C c ) of the matrix R(s) in ( 3) with respect to the system matrices A c and C c was derived in [START_REF] Döhler | Statistical decision making for damage localization with stochastic load vectors[END_REF][START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF]. The covariance of the matrices R(s) and R(s) are computed in the damaged and healthy states of the structure for a chosen s-value. These covariances are propagated to the load vector v(s) that belongs to the null space of δ R T (s) = RT (s) -R T (s) and subsequently propagated to the stress vector S(s).

For multiple mode sets M = (i, j), the sensitivity Proposition 4] and then it follows from ( 5)

J M v(s),δ R T of the load vector v(s) with respect to δ R T is computed based on v M (s) re = J M v(s),δ R T (vec( δ R(s) T )) M re with [5,
S M (s) re = J M S(s),δ R T (vec( δ R T (s))) M re ,
where

J M S(s),δ R T def = (L model(s) ) Re J M v(s),δ R T .
The above expression can be derived as in [START_REF] Döhler | Statistical decision making for damage localization with stochastic load vectors[END_REF]Theorem 5]. For multiple mode sets M , the above expression can be written for covariance computation of S M (s) in the selection of different mode sets i and j respectively

S M (s) re = J M S(s),δ R T J M R(s),( Ãc , Cc ) vec( ÃM c ) vec( CM c ) -J M S(s),δ R T J M R(s),(A c ,C c ) vec( A M c ) vec( C M c ) = J M S(s),( Ãc , Cc ) J M ( Ãc , Cc ), H vec( H) -J M S(s),(A c ,C c ) J M (A c ,C c ),H vec( H) (10)
with J M S(s),( Ãc , Cc )

def = J M S(s),δ R T J M R(s),( Ãc , Cc ) , J M ( Ãc , Cc ), H def = J ÃM c , H J CM c , H , J M S(s),(A c ,C c ) def = J M S(s),δ R T J M R(s),(A c ,C c ) andJ M (A c ,C c ),H def = J A M c ,H J C M c ,H .
After stacking the real and imaginary parts of the stress vector S M (s) for multiple mode sets M = (i, j), the total stress vector is derived as follows

S(s) def = S i (s) re S j (s) re . (11)
For multiple mode sets M = (i, j), S i (s) and S j (s) can be derived separately with [START_REF] Bernal | Damage localization from the null space of changes in the transfer matrix[END_REF]. After stacking all mode sets accordingly, it follows

S(s) = J i S(s),( Ãc , Cc ) J i ( Ãc , Cc ), H J j S(s),( Ãc , Cc ) J j ( Ãc , Cc ), H vec( H) - J i S(s),(A c ,C c ) J i (A c ,C c ),H J j S(s),(A c ,C c ) J j (A c ,C c ),H vec( H) . ( 12
)
Then covariance of the stress can be defined as cov(vec(S(s))) def = Σ S(s) and analytically, it yields

Σ S(s) = J i S(s),( Ãc , Cc ) J i ( Ãc , Cc ), H J j S(s),( Ãc , Cc ) J j ( Ãc , Cc ), H Σ H J i S(s),( Ãc , Cc ) J i ( Ãc , Cc ), H J j S(s),( Ãc , Cc ) J j ( Ãc , Cc ), H T + J i S(s),(A c ,C c ) J i (A c ,C c ),H J j S(s),(A c ,C c ) J j (A c ,C c ),H Σ H J i S(s),(A c ,C c ) J i (A c ,C c ),H J j S(s),(A c ,C c ) J j (A c ,C c ),H T . (13)
The covariance expression (13) yields a new statistical approach for damage localization using multiple mode sets based on a statistical test for each element t of the structure. In this approach, all stress components in S t corresponding to an element t are being tested in a hypothesis test, where computed stress values are related to their joint covariance (13). Since an estimate of the stress vector S t is asymptotically Gaussian distributed, a statistical evaluation is derived for each element t in a χ 2 t -test as

χ 2 t = S T t Σ -1 t S t . (14)
5. APPLICATIONS

Test case

The damage localization method has been applied to a mass-spring chain system as shown in Figure 1. The total number of Degrees Of Freedom (DOF) of the structure is 6. The damaged element in the model is simulated by decreasing stiffness by 10 % of its original value. For damaged and undamaged states, the acceleration data length for each set is N = 50,000. Data were generated from collocated white noise excitation using three sensors at elements 2, 4 and 6 in Figure 1 with sampling frequency of 50 Hz, 2% damping ratio and 5% white noise were also added to the output data. In this example, 6 modes can be identified from model of the structure using SSI. A subset of 3 modes could be used in the previous works [START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF] as the number of modes could not be bigger than the number of output sensors. For the new method, the identified modes are splitted in two mode sets namely i and j of 3 modes each.

In the application, the outcome of the damage localization results using multiple mode sets are compared with using only a single mode set. Recall that stress values close to zero indicate potentially damaged elements. From SSI, all estimated modes were chosen using a stabilization diagram procedure [START_REF] Döhler | Fast multi-order computation of system matrices in subspace-based system identification[END_REF]. Then, the system matrices in both reference and damaged states are computed [START_REF] Marin | Robust statistical damage localization with stochastic load vectors[END_REF] and for multiple mode sets, their variance is derived from the same dataset as described in Section 4.

Notice that there is always a possibility to use the same or different s-values for the computation of stress with the proposed method. Let the Laplace variable s i and s j be chosen within the range of the identified poles for mode sets i and j respectively. Then, real and imaginary parts of the stress vector and their joint covariance in (13) are computed for these Laplace variables s i and s j . Finally, 4 stress values are computed for each element t including real and imaginary parts of S i t (s i ) and S j t (s j ) for mode set i and j respectively. For each element t, these stress values are aggregated statistically in the χ 2 t -test.

Performance evaluation of damage localization

To analyse the performance of the proposed damage localization method using multiple mode sets, the results are compared with a single mode set experiment. 500 datasets are generated by Monte-Carlo simulations to evaluate the success rate of the damage localization. In order to indicate if an element is potentially damaged or not, the χ 2 t value is computed for each element. Successful damage localization means that the lowest χ 2 t value is at the damaged element. The success rate corresponds to the probability of detection or power of the test. Here, the success rate of the successful damage localization is defined as the number of occurrences of the smallest χ 2 t values at damaged elements. In order to see the influence of all s-values s = (p, q) ∈ C, a grid map in the complex s-plane has been proposed. Let real and imaginary parts of the s-value be defined as p and q. Define p = a : b : c while the discrete range [a, c] = [-3, 1] is chosen by using a step b = 0.5 and then q = d : e : f , when the discrete range [d, f]= [START_REF] Bernal | Load vectors for damage location in systems identified from operational loads[END_REF]130] is selected by step e = 2. The range of s-values has been chosen in the vicinity of the identified poles, λ i c ∈ C in Table 1. The χ 2 t -tests for all s-values are evaluated for each of the 500 datasets in order to evaluate the influence of different s-values. Then performance evaluation of the method is illustrated in a 3D-bar diagram where x-y axes indicate real-imaginary parts of s-value and z-axis corresponds to success rate of damage localization. In Section 5.3, the damage localization results are demonstrated in all elements computed at one s-value for one data set. Then, success rate of the damage localization results has been computed using a single mode set in Monte-Carlo simulation for 500 data set in Section 5.4. Finally, the success rate of the damage localization results is illustrated with the new proposed method using multiple mode sets in Section 5.5.

Localization of results in all elements at one s-value

The localization tests at all elements are computed using one of the Monte-Carlo datasets only in both damaged and healthy states. Recall that the damage position is inferred by the stress value closest to zero. For a single mode set, the computation of all stresses is done at s 1 = (-2, 51) ∈ C. In Figure 2(a), all stress values corresponding to healthy and damaged elements are presented, while the smallest stress value is correctly located in the damaged element at bar 4. In Figure 2

Success rate of the damage localization using a single mode set

In this section, the damage localization results are illustrated in Figure 3 and Figure 4 for both single mode sets i and j comprising the first and last three modes respectively. All these computations are done at Laplace variables, s = (p, q) ∈ C (i.e. recall p = 1 : 0.5 : -3 and q = 1 : 2 : 130) which are located in the vicinity of the identified poles (see Table 1). For the first 3 modes, in Figure 3, the success rate of the damage localization for the statistical χ 2 t -test is satisfactory only in the interval of the Laplace variables, s = (p, q) ∈ C with p = 1 : 0.5 : -3 and q = 8 : 2 : 64. After that, it is almost flat and cannot indicate the damage position due to modal truncation error in the region s = (p, q) ∈ C where p = 1 : 0.5 : -3 and q = 66 : 2 : 130. Note that the considered modes in this example are in the interval [START_REF] Bernal | Load vectors for damage location in systems identified from operational loads[END_REF]62] on the imaginary line (see Table1).

Similarly, for the last 3 modes, it can be seen in Figure 4 that the success rate of damage localization for the statistical χ 2 t -test is high enough in the area of the last two identified poles. In this case, in Figure 4, the good s-values are located in the interval [102, 128] on the imaginary line as are the last two identified poles, while there is no success of damage localization around the first four-identified poles because of modal truncation errors. Still, there is a slight bump around the fourth mode. So, choosing the s-value in the vicinity of the identified poles is not always perfect. In this section, the new approach is presented using mode sets to increase the success rate of damage localization in the s-plane. Hence, choose the range of the s-values s i and s j as (p i , q i ) ∈ C and (p j , q j ) ∈ C where p i = 1 : 0.5 : -3, q i = 1 : 2 : 64, p j = 1 : 0.5 : -3 and q j = 66 : 2 : 128 for the respective mode sets i and j. In the previous section, the success rate of the damage localization with SDDLV was not successful everywhere in the s-plane because of modal truncation and estimation errors. It motivates the use of multiple mode sets instead of using a single mode set alone. Therefore, the computation of stress for damage localization is derived by taking into account the information using multiple mode sets instead of using a single mode set. The joint stress evaluation derived in Section 4 is performed and its joint covariance is computed for multiple mode sets M = (i, j) containing the first three and last 3 identified poles (see Table 1).

Recall that the optimal choice of the Laplace variable for each mode set depends on the selected modes within this set. Suitable s-values have been found to be in the vicinity of the chosen modes. In the previous experiments, it has been seen that highest success rates of the damage localization are located in the vicinity of the identified poles (see Figures 3 and4) for the respective mode sets i and j. The impact of choice of s-values s i and s j has been investigated differently for the respective mode sets. The following two cases are now considered. t -test using multiple mode sets i, j in dependence of s [s i = (1 : 0.5 : -3, 1 : 2 : 64);s j = (1 : 0.5 : -3, 106)] Case 1: the s-value s i = (p i , q i ) ∈ C has been chosen with p i = 1 : 0.5 : -3 and q i = 1 : 2 : 64 in the vicinity of mode set i, while s j = (p j , q j ) ∈ C is chosen with p j = 1 : 0.5 : -3 and kept fixed at q j = 83 where it leads to a poor performance for the damage localization in both Figure 3 and Figure 4. In this case, in Figure 5, the maximum success rate is achieved in the region [s i = (1 : 0.5 : -3, 22 : 2 : 62);s j = (1 : 0.5 : -3, 83)], which is much larger and achieving higher success rate than both satisfactory region s(1 : 0.5 : -3, 20 : 2 : 62) in Figure 3 and s(1 : 0.5 : -3, 102 : 2 : 128) in Figure 4. Notice that the relatively lower rate in region [s i = (1 : 0.5 : -3, 1 : 2 : 20);s j = (1 : 0.5 : -3, 83)] of Figure 5 is still higher than the very same region in both Figures 3 and4.

Case 2: the choice of Laplace variable s i is similar as in case 1 but s j = (p j , q j ) ∈ C has been chosen with p j = 1 : 0.5 : -3 and kept fixed at q j = 106 such that it corresponds to achieving a good performance of the damage localization in Figure 4. In this case, it has been seen that the success rate of the damage localization with the new method has significantly improved the situation everywhere in the s-plane in Figure 6 compared to all previous results. When considering the statistical uncertainties of the modal parameters using multiple mode sets, the subsequent stress evaluation can improve the situation significantly.

CONCLUSIONS

In this paper, the damage localization with SDDLV has been extended considering multiple mode sets. The robustness of damage localization has been obtained by a joint statistical evaluation taking into account the information from all modes on the structure. The stress computation using multiple mode sets increases the information content about the damaged or non-damaged elements of the structure. In the application, several experiments have been carried out to evaluate the success of the damage localization in dependence of the s-value, where all stress values corresponding to an element are being tested for damage in a hypothesis test. The proposed method has increased the success rate of the correct damage localization almost everywhere in the complex s-plane compared to using a single mode set.

Figure 1 :

 1 Figure 1: Mass-spring chain system (6 DOFs)

Figure 2 :

 2 Figure 2: Localization using a single mode set: stress computation and statistical evaluation at s 1 = (-2, 51) ∈ C -three sensors, 5% output noise, 10% stiffness reduction at bar 4.

Figure 3 :Figure 4 :

 34 Figure 3: Success rate of the damage localization with statistical evaluation of χ 2 t -test using single mode set i in dependence of s

Figure 5 :Figure 6 :

 56 Figure 5: Case 1-Success rate of the damage localization with statistical evaluation of χ 2 t -test using multiple mode sets i, j in dependence of s [s i = (1 : 0.5 : -3, 1 : 2 : 64);s j = (1 : 0.5 : -3, 83)]

Table 1 :

 1 Identified Poles

	Mode i	1	2	3	4	5	6
	λ i c ∈ C (-0.25, 12.8) (-0.74, 37.4) (-1.16, 58.3) (-1.56, 77.8) (-1.96, 98.2) (-2.13, 106.6)

ACKNOWLEDGMENTS

Financial support from the Bretagne region is gratefully acknowledged.


Michael Döhler
email: michael.doehler@inria.fr

Sebastian Thöns

Efficient Structural System Reliability Updating with Subspace-Based Damage Detection Information

Keywords: reliability updating, performance models, damage detection, subspace methods

Damage detection systems and algorithms (DDS and DDA) provide information of the structural system integrity in contrast to e.g. local information by inspections or nondestructive testing techniques. However, the potential of utilizing DDS information for the structural integrity assessment and prognosis is hardly exploited nor treated in scientific literature up to now. In order to utilize the information provided by DDS for the structural performance, usually high computational efforts for the pre-determination of DDS reliability are required. In this paper, an approach for the DDS performance modelling is introduced building upon the non-destructive testing reliability which applies to structural systems and DDS containing a strategy to overcome the high computational efforts for the predetermination of the DDS reliability. This approach takes basis in the subspace-based damage detection method and builds upon mathematical properties of the damage detection algorithm. Computational efficiency is gained by calculating the probability of damage indication directly without necessitating a pre-determination for all damage states. The developed approach is applied to a static, dynamic, deterioration and reliability structural system model, demonstrating the potentials for utilizing DDS for risk reduction.

INTRODUCTION

Research focus within structural health monitoring (SHM) has been, amongst others, on the development of damage indicators [START_REF] Farrar | An introduction to structural health monitoring[END_REF] as well as on the utilization of SHM information in structural reliability analysis (e.g. [START_REF] Enright | Application of Probabilistic Fracture Mechanics to Prognosis of Aircraft Engine Components[END_REF] and [START_REF] Thöns | Monitoring Based Condition Assessment of Offshore Wind Turbine Structures[END_REF]) during the last decade. Both foci have been pursued mostly independently and only recently the research community started to interact, e.g. within the COST Action TU1402. For example, it is essential to know if a detected damage actually has a significant impact on the structural integrity. Conversely, given the information that no damage has been detected, the current structural performance is of no lesser interest. The connection of damage detection systems (DDS) and algorithms (DDA) with the evaluation of the structural performance is the aim of this study.

In order to utilize the information provided by DDS and DDA for the structural performance, usually high computational efforts for the pre-determination of DDS reliability are required. In this paper, an approach for the DDS performance modelling is introduced building upon the non-destructive testing reliability [START_REF] Gandossi | Probability of Detection Curves: Statistical Best-Practices[END_REF] which applies to structural systems and DDS containing a strategy to overcome the high computational efforts for the pre-determination of the DDS reliability. This approach takes basis in the subspace-based damage detection approach, see e.g. [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]- [START_REF] Döhler | Subspace-based damage detection under changes in the ambient excitation statistics[END_REF]. Analysing its mathematical properties regarding the probability of damage indication in connection with a structural model, we show how the reliability of the structural system can be updated using the damage detection information from measurement data. In this way, the probability of damage indication can be calculated directly without the necessity of pre-determination of the DDS reliability.

We introduce a deteriorating structural system and the DDA modelling in Sections 2 and 3. Section 3 contains a detailed description of the DDA properties which can be exploited for the approach to directly calculate the probability of indication. In Section 4, the structural system and the DDA models are utilised together for updating the structural system reliability. Section 5 contains a study on how to apply the developed approach and demonstrates which DDA characteristics lead to a reliability increase of the exemplary structural system. The conclusions (Section 6) highlight the potentials of DDA for risk reduction in general and the efficiency of the developed approach.

STRUCTURAL SYSTEM MODELLING

The performance of a structural system both in regard to system failure and system damage can be described by methods such as logical systems, Daniels Systems and Bayesian networks.

The probability ( ) S P F of a structural system failure can be calculated by integrating the joint probability density over the space of system failure S F Ω in dependency of the vectors of the system performance random variables X and the system degradation random variables D , i.e.

( ) ( )

, S F S S F P F f d d Ω = ∫ X D X D , (1)
where the system failure space S F Ω is defined with limit state functions.

Based on logical system modelling, a deteriorating structural system can be described as a series system with k n parallel subsystems consisting of j n components:

( ) ( ) ( ) , , , , 1 0 j k n n S R i j i j SL S i j j i P F P M R t M S =   = - ≤       ∪∩ with ( ) ( ) ( ) , , j,0 , 1 i j SL i i j SL R t R D t = - (2)
The vector of the system performance random variables X comprises then the resistance model uncertainties ( ) ( )

,1,1 1,1 1,1 , , , i j i j i j T R R,n ,n SL n ,n SL S n ,n M M R t R t M S S   =   X … … … . ( 3
)
The vector of the system degradation random variables contains the time dependent component damages ( ) ( )

1,1 i j T SL n ,n SL D t D t   =   D … including zero damages.

SUBSPACE-BASED DAMAGE DETECTION

Theoretical background

Automatic vibration monitoring aiming for damage detection is one of the most known and developed techniques for long term structural health monitoring (SHM) and has been recognized as an addition or alternative to visual inspections or local non-destructive testing performed manually. The rationale is that damages have an effect on the structural stiffness, and thus on the modal parameters (modal frequencies, damping ratios and mode shapes) which characterize the dynamics of the structure. A network of vibration sensors (usually accelerometers) is attached to the structure, measuring continuously the structural vibration response to ambient excitation like wind, traffic, waves or other sources. Changes in the measured signals with respect to the dynamic characteristics of the structure then indicate damage.

The subspace-based damage detection evaluates such changes with a statistical test value: vibration measurements from the current system are compared to a reference state in a subspace-based residual vector. In a hypothesis test, the uncertainties of the residual are taken into account and the respective 2 χ test statistic is compared to a threshold in order to decide if the structure is damaged or not. Based on these properties, the 2 χ test statistic is considered as the damage indicator value for damage monitoring. It is defined as follows.

The vibration behaviour of the monitored structure is assumed to be described by a linear time-invariant dynamical system

( ) ( ) ( ) ( ) F z t z t z t v t + + = C K Mɺɺ ɺ (4)
where t denotes continuous time, M, C and K

m m × ∈ ℝ
are the mass, damping and stiffness matrices, vector m z ∈ ℝ collects the displacements of the m degrees of freedom of the structure, and ( ) F v t is the external force which is usually unmeasured for long-term monitoring, being modelled as white noise. Observing the structural system (4) with a set of r acceleration sensors yields the measurements

( ( ) ) ) ( z e y t t L t + = ɺɺ (5)
where r y ∈ℝ is the measurement vector, matrix

r m L × ∈ ℝ
indicates the sensor locations and e is the measurement noise. Measurements are taken at discrete time instants ( )

k y y kτ =
, where τ is the sampling rate.

From a set of outputs

T N S ζ =
H , which deviates from zero if the system is damaged. This residual vector is asymptotically Gaussian (for large N) with zero mean in the reference state and non-zero mean in the damaged state, see also details in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]- [START_REF] Döhler | Subspace-based damage detection under changes in the ambient excitation statistics[END_REF]. Let θ be a vector of system parameters that shall be monitored, and let θ0 be its value in the reference state. Then, the respective test statistic writes as

1 1 1 1 ( ) T T T s ζ ζ - - - - = Σ Σ Σ J J J J ( 6
)
where J is the residual sensitivity with respect to θ and ( )

T ζζ Σ = E
is the residual covariance matrix. Both can be computed from measurement data in the reference state. The test statistic (6) is asymptotically 2 χ distributed with a non-centrality parameter in the damaged state.

To decide if the monitored structure is damaged or not, the test statistic is compared to a threshold λ. The threshold is typically chosen so that the probability of false alarms (type I error) is below some chosen level. In the following, the event of indication is defined as a detection of damage, i.e. when the test statistic exceeds the threshold. Hence, the test statistic s is our damage indicator value (DIV).

Performance properties

The probability distribution of the test statistic s is well-known (see [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]- [START_REF] Döhler | Subspace-based damage detection under changes in the ambient excitation statistics[END_REF]). Its 2 χ distri- bution has the following parameters: its degrees of freedom are the dimension of θ, and its noncentrality parameter in the damaged state is given by

1 T T δ δ - Σ J J
, where vector δ is linked to the change in the system parameter due to damage and to the number of data samples by 0 (

)

N δ θ θ = - . (7)
In Fig. 1 an example of the respective distributions in the reference and damaged state is given. With these known parameters, the probability of indication (PoI) can be calculated for a given damage vector δ. It holds 2 χ distribution, which is given by the gener- alized Marcum Q-function [START_REF] Shnidman | The calculation of the probability of detection and the generalized Marcum Q-function[END_REF] and can be easily evaluated. Note that an alternative to the direct computation of the PoI in [START_REF] Thöns | Structural reliability updating with stochastic subspace damage detection information[END_REF] are Monte Carlo simulations of the DIV s as in [START_REF] Thöns | Structural reliability updating with stochastic subspace damage detection information[END_REF]. In this way, the empirical distribution of the DIV is obtained, from which the PoI can be evaluated. While this approach is applicable to other DIV's than the one presented in this paper, generalizing the approach, it is computationally much more complex. In particular, it becomes unfeasible for an increasing dimension of the system parameters when the evaluation of PoI's for many damages is necessary.

T PoI s p x P x χ λ δ λ θ δ δ ∞ - = > = Σ ∫ J J 1 ( ( ;di ) ( ) ), m( )d T

UPDATING THE STRUCTRAL SYSTEM RELIABILITY WITH DAMAGE DE-TECTION INFORMATION

In order to update the structural system reliability with the damage detection information, the structural damages have to be related to the performance of the damage detection system. The performance of the damage detection is dependent on the static and dynamic properties of the system (Equation ( 4)) and the damage detection algorithm.

Each structural damage D defined in Section 2 affects the structural components. This means that if one or several system components are damaged, then the respective damage D corresponds to a particular set of matrices describing the monitored system in (4) with modifications of the matrices M, C and K. Vice versa, the structural reliability of the system is also affected by the damage D.

In order to model this relation, the system parameter vector θ in Section 3 is linked to the possible damages D, i.e. a damage D corresponds directly to a change in θ and is thus related to the damage vector δ in [START_REF] Shnidman | The calculation of the probability of detection and the generalized Marcum Q-function[END_REF]. In general this relation may be described as a functional relationship, i.e. δ = f(D). Such a "compatible" modelling can be assured when the structural components defined in Section 2.1 correspond to the elements of finite element model in (4), and the parameter vector θ is the collection of the selected damage-sensitive parameter for each structural component.

In the remainder of this paper, the damage D is modelled as a stiffness loss in the components of a structure. These components are the elements of the finite element model described in (4), and system parameter θ contains the element stiffnesses. Then, it follows from (7) that a damage D is linked to the damage vector δ by

N δ = D . (9)
The probabilities of indication for such damages are required to update the structural system reliability given DDS information as follows. The probability of structural system failure given the DDS information of no-indication, (

∩ = = . ( 10
)
The marginal probability of no-indication ( ) S P I can be calculated based on the developed approach for the DDS performance calculation and following [START_REF] Hong | Reliability analysis with nondestructive inspection[END_REF] with [START_REF] Wall | Modelling of NDT reliability (POD) and applying corrections for human factors[END_REF]. For each realization of D, the PoI's are readily computed using Equ. [START_REF] Thöns | Structural reliability updating with stochastic subspace damage detection information[END_REF], where δ depends on D e.g. as described in [START_REF] Jcss | Probabilistic Model Code[END_REF]. In contrast, the pre-computation of the PoI's for all possible D results in non-polynomial computational demands (see e.g. [START_REF] Thöns | Structural reliability updating with stochastic subspace damage detection information[END_REF]). The computation "on the fly" during the computation of the updated probability in [START_REF] Wall | Modelling of NDT reliability (POD) and applying corrections for human factors[END_REF] as introduced here removes the necessity of the PoI computation in advance. In particular, this has the advantage of only computing PoI's for the actually needed realizations of D, instead of computing all possible combinations of D. The latter becomes computationally almost infeasible when the dimension of D grows due to the exponential growth of the possible damage combinations.

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

APPLICATION

To illustrate the developed approach, a simplistic structural system of two components subjected to deterioration with a DDS comprising two sensors and a subspace based DDA (Figure 2) is considered. The structural system is described with its static, dynamic, deterioration, reliability characteristics.

( ) ( ) ( ) 2 , 1 0 c n S R i i SL S i i P F P M R t M S = =   = - ≤     ∪ (12)
The formulation contains the number of components 2 c n = with the random variables com- ponent resistance i R (dependent on the time SL t ), system loading S and their associated model

( ) i SL D t , i.e. ( ) ( ) ( )
,0 1

i SL i i SL R t R D t = - (13)
For clarity, the temporal dependence of the damage is neglected in the further.

The structural reliability model is summarised in Table 1. The system loading is represented with a Weibull distributed random variable S which results by equilibrium in the component loading i S . The loading and resistance model uncertainties and the resistance model are deter- mined according to [START_REF] Jcss | Probabilistic Model Code[END_REF] as Lognormal distributed with a standard deviation of 10%. The component probability of failure is calibrated to 3 1 10 - ⋅ by adjusting the mean of the component resistance in the undamaged state. The correlation of the resistances and the deterioration are modelled with a coefficient of correlation of 0.5 or varied.

The static and dynamic system properties are modelled with distributed component stiffness and mass subjected to a structural damping of 2% for each mode, see Table 1. The static and dynamic system behaviour is calculated with the Finite Element method. The DDS is modelled with the acceleration sensors 1 s in X-direction and 2 s in Y-direction recording the responses using the subspace-based DDA. Based on the dynamic structural system model, a reference dataset of length N = 10000 at a sampling frequency of 50 Hz is simulated in the undamaged state for both sensors from white noise excitation in order to compute the reference parameters of the damage detection method (null space S, residual sensitivity with respect to both component stiffnesses J and covariance Σ , threshold λ for type I error of 0.01). For the computation of the probabilities of indication for the damages in [START_REF] Wall | Modelling of NDT reliability (POD) and applying corrections for human factors[END_REF], a data length of N = 10000 is assumed. The distribution of the DIV takes into account the uncertainties related to the measurement data of finite length N, the unknown ambient excitation and the measurement uncertainty. Human errors in the application and operation are accounted by the multiplication of the PoI with a factor of 0.95, see [START_REF] Wall | Modelling of NDT reliability (POD) and applying corrections for human factors[END_REF].

DDA performance: probability of indication

From Equ. ( 8), the probabilities of indication of the DDA (without accounting for human errors) are determined for the damages D that are required to evaluate the probabilities in Equ. [START_REF] Wall | Modelling of NDT reliability (POD) and applying corrections for human factors[END_REF] for updating the structural reliability. In Figure 3, the respective damages are shown (black dots) at which the probability of indication is evaluated during the computation of [START_REF] Wall | Modelling of NDT reliability (POD) and applying corrections for human factors[END_REF]. A correlation of 0.5 between the damages in both components was assumed in this example. Furthermore, Figure 3 shows the theoretical PoI curve for all possible damages. It is observed that the probability of indication is very high for the majority of the damage samples whereas the size of the damage samples is on a low level. This is caused by the high sensitivity of the utilised subspace based DDA.

Updated structural reliability with DDA information

The DDS information are utilised to update the structural reliability of the deteriorated structural system. The deteriorated series system without updating by DDS information shows a slight decrease of the structural system failure probability for both increasing the initial resistance and the deterioration correlation 0 R ρ and D ρ (Figure 4). The decrease of the system reliability is slightly higher for the resistance correlation considering that the system is not highly deteriorated.

When utilizing the DDS information, the system failure probability drops and decreases with increasing correlation for both sensors 1 s and 2 s . With increasing resistance correlation, the decrease rate is comparable to decrease rate without monitoring. The decrease rate for increasing the damage correlation is significantly higher. This can be explained by the higher probability of indication for correlated damages as the dynamic system properties change is more pronounced and can thus be more reliably detected. Sensor 2 s has a better performance as the system the absolute stiffness and thus the stiffness change is in the vertical direction higher and can thus be more reliably detected.

CONCLUSIONS

An approach for the DDS and DDA performance modelling has been introduced which applies to structural systems and contains a strategy to overcome the high computational efforts for the pre-determination of the DDS reliability.

The introduced approach facilitates (1) that the probability of damage indication can be calculated directly without necessitating a pre-determination and (2) the updating of structural reliability with DDS information on structural system level. The developed approach requires a consistent structural performance modelling in terms of the static, dynamic and deterioration characteristics in order to derive the structural reliability and the DDS performance models.

With an example it is demonstrated that the DDS and DDA performance is evaluated on a system level for both the DDS and the structural system. As such the DDS performance modelling accounts for the dependencies of the structural system damage states and encompasses the measurement system (number of sensors, sensor positions, precision of the system) and the employed damage detection algorithms.

In the perspective of further research, the approach taken here can be used to design efficient DDS and DDA, quantifying the value of DDS and DDA and optimising their performance.

where 2 (

 2 ;dof , nc) p x χ is the probability density function of the non-central 2 χ distribution with dof degrees of freedom and non-centrality parameter nc. The integral is linked to the cumulative distribution function of the non-central 2

Figure 1 :

 1 Figure 1: Scheme of probability density functions (pdf) of damage detection test statistic χ 2 in the reference and in a damaged state.

Figure 2 :

 2 Figure 2: Structural system with sensor locations and sensor measurement directions

M

  , respectively. The time dependent resistance ( ) i SL R t of the component i is modelled with the initial resistance ,0 i R and the time dependent damage

Figure 3 :

 3 Figure 3: Probability of indication for damages in both structural components.

Figure 4 :

 4 Figure 4: Prior and posterior system probability of failure ( ( ) S P F and ( ) | S P F I ) in dependency of the resistance correlation, the damage correlation and the probability of component failure for different sensor positions.

Table 1 :

 1 Static and dynamic and structural reliability model properties

	Random variable		Distribu-	Mean Standard devia-
				tion		tion
	Mass per component		-	0.5	-
	Stiffness of component 1: 1 EA	-	1000	-
	Stiffness of component 2: 2 EA	-	2000	-
	Damping ratio			-	2 %	-
	Loading S			WBL	3.50	0.1
	Model uncertainty S M	LN	1.00	0.1
	Component resistance in undamaged state 0 ,i R	LN	Cali-	0.1
					brated
	Model uncertainty	M	, R i	LN	1.00	0.1
	Damage i D			N	0.07	0.03

ACKNOWLEDGMENTS

The support of COST Action TU1402 on Quantifying the Value of Structural Health Monitoring is gratefully acknowledged.


L Safari

P Amaro

J P Santos

F Fratini

Angular and polarization analysis for two-photon decay of 2s hyperfine states of hydrogenlike Uranium

Keywords: numbers: 31.10.+z, 32.10.-f, 32.30.-r, 32.80.Wr, 32.10.Fn Keywords: two-photon decay, hyperfine interaction, angular correlation, polarization correlation

The amplitude of two-photon transitions between hyperfine states in hydrogenlike ions is derived based on relativistic Dirac equation and second order perturbation theory. We study angular and linear polarization properties of the photon pair emitted in the decay of 2s states, where spin-flip and non-spin-flip transitions are highlighted. We pay particular attention to hydrogenlike uranium, since it is an ideal candidate for investigating relativistic and high-multipole effects, such as spin-flip transitions. Two types of emission patterns are identified: i) non-spin-flip transitions are found to be characterized by an angular distribution of the type W (θ) ∼ 1 + cos 2 θ while the polarizations of the emitted photons are parallel; ii) spin-flip transitions have somewhat smaller decay rates and are found to be characterized by an angular distribution of the type W (θ) ∼ 1 -1/3 cos 2 θ while the polarizations of the emitted photons are orthogonal, where θ is the angle between photons directions. Deviations due to non-dipole and relativistic contributions are evaluated for both types of transitions. This work is the first step toward exploring the effect of nucleus over the the angular and polarization properties of the photon pairs emitted by two-photon transitions.

I. INTRODUCTION

Two-photon decay in atoms and ions was introduced by Max Born's PhD student Goeppert-Mayer in 1931 [1]. Since then, many aspects of such a process, like the total decay rate and the spectral distribution, have been extensively investigated in the context of few-electron atoms and ions, both in theory and experiments [2][3][4][5][6][7][8][9][10][11][12]. Recently, some interest has been also devoted to the relativistic effects on angular and polarization properties of the two emitted photons [13][14][15][16][17][18][19], and to electronelectron interaction effects on the total decay rate [20]. Apart from fundamental interest, two-photon transitions revealed themselves as a useful tool for investigation of different physical areas and applied science. Already in 1940, for instance, Breit and Teller derived that the double photon emission was the principal cause of the decay of interstellar hydrogen atoms from their metastable 2sstate [21], while, more recently, polarization properties of the emitted photons have been employed to successfully explore quantum entanglement [22][23][24]. Furthermore, two-photon transitions have been proposed as a tool to measure weak interaction properties [25,26].

In this article, the angular and polarization properties of the photon pair emitted by the two-photon decay of 2s hyperfine states in hydrogenlike ions are presented. Par- * laleh.safari@oulu.fi ticular attention is paid to hydrogenlike uranium (U 91+ ). In fact, due to its strong electromagnetic field, high multipoles contributions that lead to spin-flip transitions are enhanced in this system. This makes hydrogenlike Uranium an ideal candidate for our studies. The angular and polarization analysis of the emitted light is carried out within the Independent Particle Approximation (IPA), i.e., by coupling the spin-angular momenta of electron and nucleus and by neglecting any hyperfine interaction between electron and nucleus. We derive the analytical expression for the transitions amplitude within IPA. This work is a first step toward exploring the effect of the nuclear angular momentum (spin) on the angular and polarization properties of the emitted photons. Our analysis may pave the way for a new route to get information on the direction and the magnitude of the spin distribution inside the nucleus (which is still quite an unraveled problem) by using two-photon angular and polarization correlations.

II. THEORETICAL BACKGROUND

A. Construction of the overall set of states

The presence of the nuclear spin has a twofold effect on the states of hydrogenlike systems. First, the energies of the atomic metastable states are slightly shifted, mainly due to the magnetic dipole interaction that nucleus and electron experience. This energy correction can be described by using first order perturbation theory with additional contributions, such as the relativistic, Bohr-Weisskopf, Breit-Rosenthal and QED contributions [27,28]. Since this energy correction does not influence the angular and polarization properties of the emitted radiation, it will be totally neglected in the following. Second, the atomic states acquire a new quantum number, usually denoted by F , that represents the total angular momentum of the overall -nucleus plus electron-system.

The overall atomic state can be described by coupling the nucleus and electron angular momenta (referred to as IPA), i. where n, κ and j are the (Bohr) principal, the Dirac and the angular momentum quantum numbers of the electron respectively, while I represents the nuclear spin. On the other hand, m I , m j and m F are the projections of the nuclear, electronic and total (nucleus plus electron) angular momenta onto the quantization axis, respectively. Finally, β is a collective label that denotes any other quantum number needed to specify the nuclear state apart from I and m I . Using standard notation, j, m j , I, m I |F, m F are Clebsch-Gordan coefficients.

To further proceed, we suppose that the nucleus does not interact with the radiation field. In the language of quantum mechanics, this equates to considering that the interaction Hamiltonian couples only electron fields through the photon emission, while it does not act on the quantum space of nuclear states. This hypothesis holds for decays which involve bound states of neutral atoms, since the energy released in such decays is far lower than the nuclear excitation energies (there are few exceptions to this, like the nucleus of Th 229 , where the first metastable excited state is 10 eV above the ground state). For highly charged ions, nuclear excitations are of the order of ∼ MeV while photon energies can take values up to a hundred of keV. This supposition represents, therefore, a first approximation, since the inclusion of nuclear-field interactions in atomic decays is anyhow of higher order in perturbation theory [29]. As a result of this assumption, we shall find in the next subsection that the radial part of the decay amplitude is characterized by only electron state components. On the other hand, we shall see that the angular part of the decay amplitude is characterized by both electron and nucleus states components, due to the coupling of their angular momenta. We shall see that the value for total spin quantum number will directly determine the shape of the angular and polarization distributions in the atomic transitions.

B. Second order transition amplitude

The theory of two-photon decay is based on the secondorder transition amplitude and has been discussed in a number of recent papers [17,18,26]. One of the characteristic features of such amplitude is that it contains a summation over the intermediate atomic states which runs through the whole atomic spectrum, including a summation over the discrete part as well as an integration over the (positive and negative) continuum. For the problem under consideration, such a summation splits up into summations over: i) the principal quantum number n ν , ii) the Dirac quantum number κ ν , iii) the total angular momentum F ν , and iv) its projection onto the quantization axis m Fν .

By using (1) and by taking into account the orthonormality of the nuclear states, the amplitude for twophoton transitions between hyperfine states takes the form

M λ1λ2 (i → f ) = -(2π) T T ′ κν mI mj ν L1 L2 M1 M2 p1 p2 Λ1 Λ2 (λ 1 ) p1 (λ 2 ) p2 [L 1 , L 2 ] 1/2 i -L1-L2-p1-p2 ξ p1 L1 Λ1 ξ p2 L2 Λ2 P T P T ′ × D L2 * M2 λ2 (ϕ 2 , θ 2 , 0)D L1 * M1 λ1 (ϕ 1 , θ 1 , 0) U T T ′ Λ1 Λ2 χ f T ν T mI mj ν χ ν T ′ i T ′ mI mj ν + 1 ↔ 2 , (2)
where λ j and k j are the helicity and wavevector of the jth photon. The term D Lj Mj λj (θ j , ϕ j ) stands for the Wigner rotation matrices of order L with angle cordinates (θ j , ϕ j ). The notation [L] stands for 2L + 1 and Λ j runs from L j -1 to L j + 1. T, T ′ = L, S denote the large (L) and small (S) components of the electron Dirac spinor, for which the factor P T is defined as P L = 1 and P S = -1. Furthermore, p 1,2 = 0, 1 and the function ξ p L Λ is given by

ξ 0 L Λ = δ L, Λ , ξ 1 L Λ =          L+1 2L+1 for Λ = L -1 - L 2L+1 for Λ = L + 1 0 otherwise . (3)
The radial part of the amplitude in Eq. ( 2) is represented by the integral U T T ′ Λ1 Λ2 , which reads

U T T ′ Λ1 Λ2 = drdr ′ r 2 r ′2 j Λ1 (k 1 r ′ )j Λ2 (k 2 r)g T * f g T T ′ Ei+ω1 g T ′ i , (4
) where g T f,i are the small and large radial components of the final and initial electron state, while

g T T ′ Ei+ω1 = nν g T ν g T ′ * ν E ν -E i -ω 1 , (5)
is the radial Green function of the process. Here E i,ν are the energies of the initial and intermediate atomic states, while T refers to the reverse radial component of T , i.e. T = L for T = S and vice versa.

The integral in Eq. ( 4) involves only electron state components. However, its evaluation is not an easy task due to the (infinite) summation over the principal quantum number n ν contained in the radial Green function. In the present work, such integral has been computed by using the Greens library [30]. Other computational techniques could be also used, such as B spline finite basis set method [31,32] The angular part of the amplitude in Eq. ( 2) is represented by the elements χ f T ν T mI mj ν and χ ν T ′ i T ′ mI mj ν therein contained and can be computed analytically:

χ f T ν T mI mj ν = mj f j f , m j f , I, m I |F f , m F f × κ f , l T f , m j f | σ • T * L2 Λ2 M2 |κ ν , l T ν , m jν χ ν T ′ i T ′ mI mj ν = mj i j i , m ji , I, m I |F i , m Fi × κ ν , l T ′ ν , m jν | σ • T * L1 Λ1 M1 |κ i , l T ′ i , m ji , (6)
where σ are Pauli matrices. The elements

κ f , l T f , m j f | σ • T * L2 Λ2 M2 |κ ν , l T ν , m jν and κ ν , l T ′ ν , m jν | σ • T * L1 Λ1 M1 |κ i , l T ′
i , m ji have been already discussed elsewhere [14,33] and will not be here recalled.

The initial and final states involved in the two-photon transitions which we shall analyze below are unpolarized. It has been shown that, for this case, we may arbitrarily adopt the quantization axis (ẑ) along the momentum of the "first" photon: ẑ k1 [14][15][16][17][18]. We furthermore adopt x such that the xz-plane is the reaction plane (plane spanned by the photons directions). Figure 1 sketches the geometry we consider for the decay. Within this geometry, the Wigner matrices in Eq. ( 2) simplify to

D L1 * M1 λ1 (ϕ 1 , θ 1 , 0) = δ M1, λ1 and D L2 * M2 λ2 (ϕ 2 , θ 2 , 0) = d L2 M2 λ2 (θ), where d L M λ (θ)
is the reduced Wigner matrix and θ ≡ θ 2 is the polar angle of the second photon, which coincides, in the chosen geometry, with the angle between the photons directions (opening angle). Hence, the relative photons directions are uniquely determined by assigning the opening angle θ, which will be the independent variable for plotting the angular distributions. FIG. 1: (color online). Geometry considered for the twophoton emission. The propagation direction of the first photon is adopted as z direction. x is chosen such that xz is the reaction plane (plane spanned by the photons directions). θ is the angle between the photons directions, while angles χ1,2 define the linear polarizations of the first and second photon respectively with respect to their respective polarization planes. The polarization plane of the first (second) photon is denoted by A (B) and represents the plane orthogonal to the photon direction.

Since part of this work is devoted to analyze photons linear polarizations, further details concerning the detection geometry must be provided before proceeding with the analysis. In Fig. 1, we show how the photon polarizations may be defined in a case experiment. The polarization of each photon is measured in the "polarization plane", which is the plane orthogonal to the photon direction. In Fig. 1, the polarization planes of the first and second photon are denoted by A and B, respectively. Each detector is supposed to have a transmission axis, along which the linear polarization of the photon is measured. Such a transmission axis is rotated by an angle χ with respect to the reaction plane shown in Fig. 1 by the red dashed lines. Finally, each detector is supposed to work as a filter: Whenever a photon hits it, the detector either gives or does not give off a "click", which would respectively indicate that the photon has been measured as having its linear polarization along the direction χ or χ + 90 • .

C. Definition of angular and polarization correlations

Within IPA, equation (2) represents the relativistic transition amplitude for the two-photon decay between hyperfine states in hydrogenlike ions. It contains the complete information on the emitted radiation. Assuming that the ion is initially unpolarized and that the polarization of the final atomic state remains unobserved, taking into account the axes geometry chosen for the twophoton emission, and using the well-known relations be- tween linear and circular polarization bases [START_REF] Rose | Elementary Theory of Angular Momentum[END_REF], we can write the polarization-dependent differential decay rate as a function of the opening angle θ [3]:

W χ1 χ2 (θ) ≡ dw χ1χ2 d cos θ = 8π 2 2F i + 1 mF i mF f λ1λ2 λ ′ 1 λ ′ 2 × dω 1 ω 1 ω 2 4(2π) 3 c 2 e i(λ1-λ ′ 1 )χ1 e i(λ2-λ ′ 2 )χ2 × M λ1λ2 M λ ′ 1 λ ′ 2 * . (7)
In this article, the integration over the photon energies is numerically carried out by using the trapezoidal rule method. The number of points we used for the numerical integration has been checked to provide a precision of one percent. Hereafter, the function W χ1 χ2 (θ) shall be called "polarization correlation". It represents the probability density of detecting the emitted photons at the opening angle θ with defined linear polarizations χ 1 and χ 2 . Finally, by summing over the photons polarizations, we define the "angular correlation" as

W (θ) ≡ dw d cos θ = 8π 2 2F i + 1 dω 1 ω 1 ω 2 (2π) 3 c 2 × mF i mF f λ1λ2 M λ1λ2 (i → f ) 2 , (8)
which represents the probability density of detecting the emitted photons at the opening angle θ, irrespectively of their polarizations.

W χ 1 χ 2 (θ) s -1 W χ 1 χ 2 (θ) s -1 χ 2 = 0 (deg) χ 2 = 90 (deg)
FIG. 3: (color online). Polarization correlations in hydrogenlike 235 92 U ion. The function W χ 1 χ 2 (θ) is shown for the transitions 2s 1/2 (F = 3) → 1s 1/2 (F = 3) (solid-black curve) and 2s 1/2 (F = 3) → 1s 1/2 (F = 4) (dashed-green curve). The four polarization configurations (χ1, χ2 = 0 • , 90 • ) are displayed.

III. RESULTS AND DISCUSSIONS

Here we analyze the angular and polarization correlations defined in Eqs. ( 8) and ( 7) respectively, for decays of hyperfine 2s states in hydrogenlike ions, with special attention to hydrogenlike 235 92 U ion, whose nuclear spin is I = 7/2 [START_REF] Stone | [END_REF]. The function W (θ) obtained for 2s 1/2 (F = 4, 3) → 1s 1/2 (F = 4, 3) transitions in hydrogenlike 235 92 U ion is displayed in Fig. 2. The full multipoles and the electric dipole (E1E1) contributions are separately displayed. Within the dipole approximation, the angular correlation for non-spin-flip transitions can be well described by the familiar shape W E1E1 (θ) ∼ 1 + cos 2 θ [10]. However, the full-multipole calculation shows some asymmetric deviations from such a shape. This effect is already known from the past literature, where it has been showed that high multipoles contribute with terms of the type ∼ cos θ to the angular correlation in 2s 1/2 → 1s 1/2 transitions in highly charged ions [13].

On the other hand, for spin-flip transitions, the angular correlation within the dipole approximation is well described by the function W E1E1 (θ) ∼ 1 -1/3 cos 2 θ. This emission pattern is typical for two-photon transitions of the type J T OT = 1 (0) → J T OT = 0 (1), where J T OT is the total angular momentum of the system which undergoes the decay. The two-photon decay (1s 2s) 3 S J=1 → (1s 1s) 1 S J=0 in heliumlike ions, where J is the total angular momentum of the two-electron system, shows approximately the same behavior [17]. As in the previous case, the full-multipole calculation shows remarkable asymmetric deviations from the symmetric shape. Quantitatively, the ratio W (π)/W (0) is ≃ 1.09 for non-spin-flip transitions and ≃ 1.35 for spin-flip transitions.

For low-charged ions, we find that the angular correlation is fully described by the functions ∼ 1 + cos 2 θ and ∼ 1-1/3 cos 2 θ for non-spin-flip and spin-flip transitions, respectively (i.e., the full-multipole calculations coincides with the calculations performed within the dipole approximation).

From the figures, we also notice that spin-flip transitions are overall suppressed with respect to non-spin-flip transition, which is to be expected in view of the fact that the electric dipole, which is the leading multipole here, conserves the electron spin if evaluated nonrelativistically. This entails that the curves we obtained for spin-flip transitions are fully determined by relativistic and high-multipole contributions.

We now turn to analyze the polarization correlations for some of the hyperfine transitions considered in Fig. 2. In Fig. 3, we plot the function W χ1χ2 (θ) as obtained for the transitions 2s 1/2 (F = 3) → 1s 1/2 (F = 3) and 2s 1/2 (F = 3) → 1s 1/2 (F = 4). We see that, in general, photons coming from spin-flip transitions and nonspin-flip transitions have mainly orthogonal and parallel linear polarizations respectively. We find that this polarization scheme holds perfectly (i.e., without deviations) for low-Z hydrogenlike ions. However, for hydrogenlike Uranium, as well as for any highly charged ions, sizeable deviations are evident, as can be seen in Fig. 3.

Hydrogenlike heavy ions can be nowadays efficiently produced in storage rings [36]. The energy of the emitted radiation in hydrogenlike Uranium is in the range of hard X-rays. An experimental polarization analysis of such energetic photons would be nowadays possible through the use of Compton polarimeters [37][38][39][40][41]. By analyzing the decay spectrum, a conventional photon-photon coincidence measurements enables one to distinguish twophoton decay events from the dominant single-photon M1 decay channel [42,43]. Therefore, information on the polarization state of two photons can be achieved by selecting events which have been recorded in coincidence by two polarimeters and which have the desired scattering angle [26,[44][45][46][47][48].

IV. SUMMARY AND PERSPECTIVES

In summary, the amplitude for two-photon transitions between hyperfine states in hydrogenlike atoms has been calculated. By using such amplitude, the angular and linear polarization properties of the photon pair emitted in two-photon decays of 2s hyperfine states have been investigated within second-order perturbation theory and the Dirac relativistic framework. Special attention has been paid to hydrogenlike 235 92 U ion. Results have been showed for the transitions 2s 1/2 (F = 4, 3) → 1s 1/2 (F = 4, 3). It has been possible to identify two emission patterns: i) two-photon non-spin-flip transitions are found to be characterized by an angular distribution approximately of the type W (θ) ∼ 1 + cos 2 θ and by photon polarizations apparallel one to another; ii) two-photon spinflip transitions have somewhat smaller decay rate and are found to be characterized by an angular distribution approximately of the type W (θ) ∼ 1 -1/3 cos 2 θ as well as by photon polarizations approximately orthogonal one to another. Deviations to this patterns come from highmultipoles and relativistic contributions, are negligible for low-charged ions, and are of size 1 to 25% in hydrogenlike 235 92 U ion. This article is the first step toward exploring the nuclear spin effect on the angular and polarization properties of the photon pair emitted in two-photon decays. This study might pave the way for a new route to get information on the direction and the magnitude of the spin distribution inside the nucleus, which is still quite an unraveled problem [49].

  e by |n, β; F, I, κ, m F = mI , mj j, m j , I, m I |F, m F |n; κ, m j |β; I, m I , (1)

FIG. 2 :

 2 FIG. 2: (color online). Angular correlations hydrogenlike235 92 U ions. The function W (θ) is shown for the transitions 2s 1/2 (F = 4, 3) → 1s 1/2 (F = 4, 3). The dashed-red curve refers to the electric dipole approximation while the solidblack curve refers to the full multipoles contribution.

V. ACKNOWLEDGMENTS F.F. acknowledges support by Fundação de Amparo à Pesquisa do estado de Minas Gerais (FAPEMIG) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). F.F. and L.L. acknowledge support by the Research Council for Natural Sciences and Engineering of the Academy of Finland. J. P. S. and P. A. acknowledge support by FCT -Fundação para a Ciência e a Tecnologia (Portugal), through the Projects No. PEstOE/FIS/UI0303/2011 and PTDC/FIS/117606/2010, financed by the European Community Fund FEDER through the COMPETE -Competitiveness Factors Operational Programme.


Pauline Collon
email: pauline.collon@univ-lorraine.fr

Guillaume Rongier

Pierre Anquez

Pauline Mourlanette

Karst simulation with Lindenmayer-systems and ODSIM

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Abstract

Karstic systems are geological structures that strongly impact underground flows. Despite intensive explorations by speleologists, they remain partially described as many conduits are not accessible to humans. Paleokarsts are buried karstic systems with a significant reservoir potential. But they are not easily identifiable on seismic images. In those contexts, a huge uncertainty subsists on the network location and the conduit geometry. Stochastic simulations help to better assess that uncertainty. The difficulty is to reproduce the system connectivity at different scales while integrating as much geological knowledge as possible without involving poorly constrained parameters (e.g. paleo-climate, boundary conditions...).

In this paper we propose to work on two aspects and scales of karstic systems. At large scale, we stochastically simulate karst network skeletons with a new method based on a formal grammar, the Lindenmayer-system. Based on an alphabet, an axiom and user-defined rules, the method puts together segments to build the network skeleton. The definition of proper rules and the introduction of karst-dedicated parameters generate curves reproducing the complex architectures encountered in those systems, mixing branchwork and anastomotic patterns. At the conduit scale, we propose to build a 3D envelope around these skeletons with an enhanced Object Distance based Simulation Method. It uses a custom distance field from the skeleton which takes into account geological features influencing karstogenesis (horizons, faults or fractures). This controls the first-order shape of the conduits. It is then combined to a custom random threshold controlling finer-scale features of the conduits. This threshold is generated with several parameter values depending on the involved geological structures. This workflow is demonstrated on a synthetic case, showing the potentialities of the approach at both scales. Results are encouraging and various improvements are in focus. Data conditioning, both to karst observations and local shape information has to be enhanced. The network simulation has the advantage to be grid-free, meaning that no background grid is needed to perform the simulation. Thus, it avoids the stair-step effect that can be observed in other techniques. On the opposite, the method used to simulate the conduit shapes relies on a grid, necessary to compute the distance fields and to perform the threshold geostatistical simulation. For detailed conduit geometry, the grid requires a high resolution, which impacts directly the computational efficiency. Finally, it would be interesting to test the approach on a real dataset and to develop a coupling with a flow simulator to evaluate the impact of the shape and of the network connections on the flow response.

Acknowledgements

This work was performed in the frame of the RING project at Université de Lorraine. We would like to thank the industrial and academic sponsors of the Gocad Research Consortium managed by ASGA for their support. We also thank Paradigm for providing the SKUA-GOCAD software and API.