M Zribi

Member IEEE F Kotti

Z Lili-Chabaane

N Baghdadi

N Ben Issa

R Amri

ANALYSIS OF SOIL TEXTURE USING TERRASAR X-BAND SAR

Keywords: Radar, TERRASAR-X, soil texture, soil clay, modelling, mapping

In this paper, it is proposed to use TERRASAR-X data for analysis and estimation of soil surface texture. Our study is based on experimental campaigns carried out over a semi-arid area in North Africa. Simultaneously to TERRASAR-X radar acquisitions, ground measurements (texture, soil moisture and roughness) were made on different test fields. A strong correlation is observed between soil texture and a processed signal from two radar images, the first acquired just after a rain event and the second corresponding to dry soil conditions, acquired three weeks later. An empirical relationship is proposed for the retrieval from radar signals of clay content percent. Soil texture mapping is proposed over the study site, which includes bare soils and olive groves.

I. INTRODUCTION

The implementation of sustainable agricultural and environmental management requires an improved understanding of the soil, at increasingly finer scales. Conventional soil sampling and laboratory analyses cannot efficiently provide this information, because they are slow and expensive. Different methodologies have been proposed over the last decades, for the estimation of soil textures using specific multi-and hyper-spectral sensors, allowing the rapid and relatively inexpensive measurement of soil properties (eg. [START_REF] Selige | High resolution topsoil mapping using hyperspectral image and field data in multivariate regression modeling procedures[END_REF][START_REF] Lagacherie | Estimation of soil clay and calcium carbonate using laboratory, field and airborne hyperspectral measurements[END_REF][START_REF] Gomez | Continuum removal versus PLSR method for clay and calcium carbonate content estimation from laboratory and airborne hyperspectral measurements[END_REF]). Remote sensing can provide soil information over large geographical areas, but has the drawback of detecting, in most cases, only the top few millimetres at the surface of the soil profile.

Active remote sensing has shown considerable potential for the characterisation of different soil parameters, moisture and roughness in particular. In recent years, some studies have illustrated the potential of SAR measurements and very high resolution data for the retrieval of soil texture. [START_REF] Baghdadi | Analysis of TerraSARX data and their sensitivity to soil surface parameters[END_REF][START_REF] Paris Anguela | Analysis of local variation of soil surface parameters with TerraSAR-X radar data over bare agricultural fields[END_REF] have already observed TerraSAR-X SAR data variations due to soil texture within agricultural plots. This work was supported by two programs RESAMED (French MISTRAL Program) and the WASSERMED European project They demonstrated that the SAR signal is 3 dB weaker for a soil composition with more clay (32% clay, 64.5% silt and 3.5% sand) than crust variations. Soils with slaking crust have a greater hydric inertia than soil without crust. Consequently, following rainfall or dry events, soil moisture in the upper for a soil with less clay (17% clay, 79% silt and 4% sand). In a recent study, Aubert et al., [START_REF] Aubert | Analysis of TerraSAR-X data sensitivity to bare soil moisture, roughness, composition and soil crust[END_REF] show a high correlation between the radar signal and the hydric evolution of soil centimeters may differ between the two soil structures, resulting in variations in the TerraSAR-X signal within the field. Thus, it is sometimes possible to track surface degradation due to the slaking process using the TerraSAR-X sensor. Our objective in the present study is to estimate the soil clay percent using TERRASAR-X radar data over a semi-arid region situated in North Africa.

II. DATA BASE

A. Study site

The study area is located on the Kairouan plain [START_REF] Zribi | Multi-scale estimation of surface moisture over pluvial and irrigated fields in a semi-arid region using ASAR/ENVISAT radar data[END_REF], situated in central Tunisia (9°23' -10°17' E, 35° 1'-35°55' N). The climate in this region is semi-arid, with an average annual rainfall of approximately 300 mm per year, characterised by a rainy season lasting from October to May, with the two rainiest months being October and March. As is generally the case in semi-arid areas, the rainfall patterns in this area are highly variable in time and space. The landscape is mainly flat, and the vegetation in this area is dominated by agriculture (cereals, olive groves, fruit trees and market gardens).

B. Satellite database

Three TerraSAR-X images (~9.65 GHz) were acquired at high incidence angles (35°) in the horizontal polarisation. The images were taken with a horizontally polarized beam, with a ground pixel spacing of 1m.

Radiometric calibration of the TerraSAR-X signals was carried out [START_REF] Fritz | TerraSAR-X Ground Segment Level 1b Product Format Specification[END_REF]. All of the radar images were then geo-referenced using optically geo-referenced images.

C. Ground measurements over the test fields

The field campaigns presented in this study were carried out on the 8 th , 12 th , 19 th and 30 th of March, 2010. The first of these dates (March 8 th ) was one day after a period of heavy rain. The soil was wet on the first day, following which its surface progressively dried up due to the influence of strong solar radiation and wind. There was no further precipitation during the 21 day measurement period.

Simultaneously with the radar measurements, ground truth measurements of soil roughness, moisture content, bulk density and texture were carried out on several test fields. A larger number of field plots were sampled during the full set of campaigns.

-Soil texture Soil texture measurements were carried out over thirty test fields. For each test field, several soil samples were made and the measurements were processed in the laboratory. These indicated clay percents ranging between 12 and 53%. The sand percents lay in the range between 21% and 71%, with the highest values being observed in non-irrigated olive groves, which are particularly well adapted to sandy soil.

-Soil roughness

Soil roughness measurements were made using a pin profiler (with a total length of 1 m and a resolution of 2 cm). In order to improve the accuracy of the roughness computations, approximately 10 profiles were recorded for each field. As the surface height profile is considered to be ergodic and stationary, we can compute the correlation function for each profile [START_REF] Verhoest | On the soil roughness parameterization problem in soil moisture retrieval of bare surfaces from Synthetic Aperture Radar[END_REF], and derive two statistical parameters: the rms surface height (vertical scale of roughness), and the correlation length (l), which represents the horizontal scale over which similar roughness conditions are detected. The rms height values varied between 0.6cm and 1.5cm, and the correlation length between 3 and 8cm.

-Soil moisture Gravimetric soil moisture samples were collected as a function of the agricultural plot size: about 10 to 20 localized samples were taken per test field, at depth in the range 0-5 cm. The gravimetric soil moisture content was calculated by drying the samples at 105°C for 24 h. Table1 illustrates the mean value of 0-5cm soil surface volumetric moisture over the test fields, for each campaign day. We observe a mean decreasing of 16.65% between the 8 th and 30 th of March.

-Land use mapping

Land use classification is based on a decision tree as detailed in [START_REF] Zribi | Multi-scale estimation of surface moisture over pluvial and irrigated fields in a semi-arid region using ASAR/ENVISAT radar data[END_REF], using two types of satellite data: four SPOT images, and SRTM data. We established eight classes of land use: non-irrigated olive groves, irrigated olive groves, irrigated winter vegetables, irrigated summer vegetables, bare soils, urban areas, mountainous areas, water cover and areas of coastal salt flats ("sebkhas"). Validation of these remotely sensed classifications, based on ground verification over more than 100 fields with different types of land uses, reveals an accuracy of around 94%.

III. ANALYSIS OF THE RELATIONSHIP BETWEEN SOIL TEXTURE AND RADAR SIGNAL

A. Analysis of TERRASAR-X data behaviour

The Integral Equation Model "IEM" backscattering model [START_REF] Fung | Backscattering from a randomly rough dielectric surface[END_REF] is used for radar signal simulations. It has as input roughness, moisture and texture parameters. For roughness description, we consider the rms height, the correlation length and an exponential height correlation function. Fig. 1 provides a comparison between radar data and backscattering simulations. In general, we observe an over-estimation in the backscattering simulations, with respect to the real radar data. This over-estimation could be explained by various factors, including: the roughness description [START_REF] Verhoest | On the soil roughness parameterization problem in soil moisture retrieval of bare surfaces from Synthetic Aperture Radar[END_REF], instrumental errors [START_REF] Lievens | Error in Radar-Derived Soil Moisture due to Roughness Parameterization: An Analysis Based on Synthetical Surface Profiles[END_REF][START_REF] Verhoest | On the soil roughness parameterization problem in soil moisture retrieval of bare surfaces from Synthetic Aperture Radar[END_REF], and finally the soil moisture estimation. This backscattering behaviour is confirmed in Fig. 2, showing the temporal variation of the radar signals for four different test fields. For all test fields, we observed nearly stable radar signal between the 19 th and 30 th of March. This result is with the decreasing of soil moisture measured at depth of 0-5 cm, which ranged from approximately 15% down to 9% between these two dates. This apparent discrepancy could be explained mainly by the fact that, in the first millimetres of soil, which have the greatest influence on the X-band radar signals, there was a priori no variation in soil moisture between these two dates, as a consequence of the strong level of evaporation on the studied site. This result illustrates some limitations of X-band data for the estimation of soil moisture under such semi-arid climate conditions. In this context, we analyse the potential of TERRASAR-X data for the estimation of soil texture, as described in the following section. The main idea behind this analysis is to exploit differences in the rate at which the first millimetres of soil dry out, as a function of the percent of clay in the soil composition.

B. Developed approach

The proposed approach is based on change detection analysis. We analysed the difference between radar signals measured on the 8 th and 30 th of March, the first of these dates being just one day after a strong rainfall event (P=40mm), and the second date being three weeks later, with no rainfall events occurring between these two dates. We considered two types of test field: bare soils and non-irrigated olive groves. In both cases, the difference is considered as a function of soil moisture variations only. For all test fields, we verified that there was no temporal variation in roughness.

In the case of the olive groves, the trees are sparsely distributed, with the vegetation cover corresponding to approximately 5% of the soil surface area. In order to remove the influence of the vegetation, and to consider bare soil measurements only, we used a high-resolution (2.5m) optical image to identify, and then eliminate the contribution of each of the olive trees. A classification system was then applied to the images, to eliminate the tree-covered zones, thus conserving the bare soil only for further analysis, as illustrated by the example shown in Fig. 3, corresponding to two test fields. We included approximately two additional pixels around the classified tree areas, in order to take the influence of the radar incidence angle into account.

  clay f   0  (2)
From ground measurements, the clay percent ranged between 12% and 53%. It can be seen that the processed backscattered radar signals have a clear dependence on the clay content (%) of the soil, with a high correlation coefficient, R 2 , of approximately 0.6. This relationship arises from the dependence of the radar signals on soil moisture content variations. An increase in the clay percent leads to higher moisture values one day after the first rainfall, and therefore to a stronger change in moisture in the first millimetres than in the case of a low clay percent. This strong difference in soil moisture, due to large differences in clay content, certainly led to strong variations in backscattered radar signal levels, between the two observation dates.

In order to validate the theoretical robustness of the behaviour observed in Fig. 4, we ran backscattering coefficient simulations with the IEM model, using as inputs moisture and texture measurements, together with various roughness values taken within the domain of validity of the IEM (s=0.3 cm, s=0.5 cm, s=0.7cm and l=6cm).

C. Validation of our methodology

Validation of our soil texture retrieval algorithm was based on comparisons between ground-truth estimations of clay percent and measurements derived from TERRASAR-X observations using the inversion model:

D. Soil texture mapping

Following the process described above, it was possible to apply our validated model to TERRASAR-X images. On the basis of a land-use map, we took only bare soil fields and olive groves into account. In the case of the olive groves, we removed the trees from the radar images, using SPOT high resolution (2.5 m) images. In order to eliminate the effects of local terrain heterogeneities and speckle in the radar signals, the soil texture was estimated over cells defined by 50 x 50 pixel areas, corresponding to 50 x 50 m cells. Fig. 6 shows the resulting soil texture map over the studied area. As confirmed by ground measurements, a generally lower percent of clay is found in the soil of olive groves. This paper proposes a methodology which can be successfully applied to the mapping of soil texture, using TERRASAR-X radar data. Our approach is based on the fact that, after a rainfall event, the variation moisture content of the soil surface depends on the soil"s texture.

Our methodology is based on the detection and analysis of change, under different conditions: -No change in roughness or vegetation cover: in this case, the period between two successive radar images should be short.

-Soil texture estimations are derived from processed radar signals, and then averaged over cells with 50 x 50 sample spacing.

-The soil texture is retrieved by inverting the relationship between the processed signals and the soil clay percent.

Fig. 1 .

 1 Fig. 1. Inter-comparison between TERRASAR-X measurements and backscattered signal simulations with the IEM model, over various test fields

Fig. 2 .

 2 Fig. 2. Measured radar signals as a function of time, over four different test fields. The points correspond to the mean radar signal observed over the test fields, on three different dates (8 th , 19 th and 30 th of March).

Fig. 3 .

 3 Fig. 3. Classification of one olive grove: inside each olive field, two classes of coverage can be identified: olive trees in blue and bare soil in yellow.

Fig. 4

 4 Fig. 4 illustrates the relationship between processed radar signals and the clay percent in the soil, for different test fields:

Fig. 4 .

 4 Fig. 4. Observed relationship between clay percent and processed radar signals, over different test fields.

Fig. 5 .

 5 Fig. 5. Inter-comparison between ground clay measurements and estimated values, based on the proposed empirical model

Fig. 6 .

 6 Fig. 6. Illustration of soil clay percent mapping over the studied site IV. CONCLUSIONS


Keywords: Causality, Causal inference, Coronary heart disease, Degradation model, Epidemiology, Evidence synthesis, Inverse Gaussian distribution, Myocardial infarction, Stochastic processes

and their risk factors. The fit of the model using the evidence synthesis approach was satisfactory and the effects of the four risk factors were highly significant.

Evidence synthesis through a degradation model applied to myocardial infarction

Received: date / Accepted: date Abstract We propose an evidence synthesis approach through a degradation model to estimate causal influences of physiological factors on myocardial infarction (MI) and coronary heart disease (CHD). For instance several studies give incidences of MI and CHD for different age strata, other studies give relative or absolute risks for strata of main risk factors of MI or CHD. Evidence synthesis of several studies allows incorporating these disparate pieces of information into a single model. For doing this we need to develop a sufficiently general dynamical model; we also need to estimate the distribution of explanatory factors in the population. We develop a degradation model for both MI and CHD using a Brownian motion with drift, and the drift is modeled as a function of indicators of obesity, lipid profile, inflammation and blood pressure. Conditionally on these factors the times to MI or CHD have inverse Gaussian (IG) distributions. The results we want to fit are generally not conditional on all the factors and thus we need marginal distributions of the time of occurrence of MI and CHD; this leads us to manipulate the inverse Gaussian normal distribution (IGN ) (an IG whose drift parameter has a normal distribution). Another possible model arises if a factor modifies the threshold. This led us to define an extension of IGN obtained when both drift and threshold parameters have normal distributions. We applied the model to results published in five important studies of MI and CHD

Introduction

Most often in epidemiology, risks of events are modeled using a proportional hazard assumption. Degradation models may be closer to physiological mechanisms in many cases; indeed it is often the case that an event occurs when a degradation process reaches a certain threshold [START_REF] Doksum | Gaussian models for degradation processes-part i: Methods for the analysis of biomarker data[END_REF][START_REF] Aalen | Understanding the shape of the hazard rate: a process point of view[END_REF][START_REF] Hashemi | A latent process model for joint modeling of events and marker[END_REF][START_REF] Lee | Threshold Regression for Survival Analysis: Modeling Event Times by a Stochastic Process Reaching a Boundary[END_REF][START_REF] Aalen | Survival and event history analysis: a process point of view[END_REF]. A good example is myocardial infarction (MI). It has been shown that MI is most commonly due to occlusion (blockage) of a coronary artery following the rupture of an atherosclerotic plaque [START_REF] Hansson | Inflammation, atherosclerosis, and coronary artery disease[END_REF][START_REF] Nicholls | Relationship between LDL, HDL, blood pressure and atheroma progression in the coronaries[END_REF]. A Brownian motion with positive drift seems well adapted to describe the progressive growth of atheroma, with MI occurring when this process reaches a certain threshold. One advantage of this approach is that we can link MI with broader coronary heart disease (CHD) events which happen before MI when occlusion is not complete but the heart already suffers from hypoxia. These CHD events may occur when the atheromatous process reaches a threshold below that required for MI. One of the reason to use such a model is to express the effect of risk factors. Several risk factors are already known for MI. Most analyses focus on one particular risk factor rather than presenting a global model. Few works have attempted to develop more global dynamic analysis: [START_REF] Wilson | Prediction of coronary heart disease using risk factor categories[END_REF], using data from the Framingham study, developed prediction scores using indicators of lipid profile, diabetes, obesity, blood pressure and tobacco consumption as explanatory variables in a conventional Cox model; [START_REF] Gamborg | Dynamic path analysis in life-course epidemiology[END_REF] used the dynamic path analysis of [START_REF] Fosen | Dynamic path analysis -a new approach to analyzing time-dependent covariates[END_REF] to take into account the possible evolution of obesity and blood pressure, but they did not take into account other factors.

However, for calibrating complex dynamic models, cohorts having recorded all the relevant risk factors may be lacking; in particular, a factor like C-reactive pro-tein (CRP), a marker of inflammation, has attracted interest only recently [START_REF] Ridker | Inflammation, aspirin, and the risk of cardiovascular disease in apparently healthy men[END_REF] so that this marker is not available in historical studies. So, there is an interest in developing a method for synthesizing evidence from different studies with the following potential advantages: (i) incorporating recently discovered factors, (ii) investigating to which extent the results are consistent across studies, (iii) gaining power. Synthesis studies are different from meta-analysis [START_REF] Van Houwelingen | Advanced methods in metaanalysis: multivariate approach and meta-regression[END_REF]: while meta-analyses aim at giving a global conclusion based on several studies which have estimated the same parameter (the effect of treatment for instance), synthesis analyses aim at incorporating information of different types for estimating the parameters of a global model (parameters which may not appear in any single study). Evidence synthesis (other than meta-analyses) is chiefly represented by "Bayesian synthesis" [START_REF] Raftery | Inference from a deterministic population dynamics model for bowhead whales[END_REF], although a likelihood approach was suggested by [START_REF] Schweder | Bayesian synthesis or likelihood synthesis-what does borel's paradox say?[END_REF]; this has also a link with so-called indirect methods [START_REF] Jiang | The indirect method: inference based on intermediate statistics -a synthesis and examples[END_REF]. The literature on evidence synthesis based on dynamical models is rather scarce; an example in the field of HIV epidemiology is [START_REF] Presanis | Bayesian evidence synthesis for a transmission dynamic model for HIV among men who have sex with men[END_REF] who used a Bayesian approach.

The aim of the paper is to present an approach based on a degradation model for synthesis of information of different types coming from different studies and to apply it to the epidemiology of MI and CHD. The paper is organized as follows.

In section 2 we present a frequentist approach to evidence synthesis. Section 3 proposes a degradation model for MI and (CHD); on our way we come upon a new distribution, the inverse Gaussian normal-normal (IGN N ) distribution. Section 4 presents the application of the approach to several studies of MI and CHD. Section 5 concludes.

2 An approach to evidence synthesis through a dynamical model Consider a dynamical statistical model describing the causal influences of different processes. It is expressed in terms of a possibly multivariate stochastic process, the possible laws of which are indexed by a parameter vector θ, θ ∈ Θ. The parameter θ could be estimated from observations of one study. It may happen that no study contains information on all the processes involved in the model. In that case it may be necessary to synthesize the information of several studies. This has the advantage of giving a more robust result (not relying on a single study) and also gives the opportunity of examining whether the different studies have consistent results. In a synthesis analysis we assume that study k gives an estimate Qk of a quantity Q k , accompanied by a standard error σ k that is assumed known for sake of simplicity. For instance Q k may be the incidence of a certain condition in one study, the incidence of another condition in another study, a relative risk in yet another study. So let us admit that we can define Q k and express it as a function of the parameters θ: we have Q k (θ).

We propose to consider that Qk is like an observation of Q k (θ * ), where θ * would be the true parameter value if the model was well-specified. Typical Qk are proportions or maximum likelihood estimates of relative risks. In both cases they are asymptotically normal. Thus we can write a contribution to a pseudologlikelihood:

L k (θ) = - ( Qk -Q k (θ)) 2 2σ 2 k . ( 1
)
If a study gives several contributions, Qk is a vector which is assumed to have a multinormal distribution with variance matrix Ω k . The contribution to the pseudolikelihood is then:

L k (θ) = - 1 2 ( Qk -Q k (θ)) T Ω -1 k ( Qk -Q k (θ)). (2)
Summing over k we get the total pseudo-loglikelihood: L(θ) = K k=1 L k (θ). The estimate θ maximizes L(θ). If L(θ) has a maximum the problem is identifiable.

If it does not, we can either add new studies, or add a priori knowledge on the parameters or other quantities Q k in a Bayesian spirit.

It is good here to give an example. We have developed a degradation model (described in section 4) which gives the joint distribution of the time of occurrence of MI and of risk factors; these distributions are indexed by a parameter vector θ which includes in particular the effects of the risk factors on the drift, the threshold parameter(s), and the correlations between risk factors (see Table 9). The Atherosclerosis Risk In Communities (ARIC) surveillance study (National Heart Lung and Blood Institute, 2006) provided an estimate of the incidence of MI among American men. For instance, the estimated incidence in the 55-64 age stratum based on the observation of 32,572 person-years was Qk = 6.26‰. Here we have:

Q k (θ) = 1 10 F θ (65) -F θ (55) 1 -F θ (55) , (3)
where F θ is the marginal cumulative distribution function of time of occurrence of MI in our model given the value θ of the parameters. The variance of Qk in the person-years method is classically estimated by

σ 2 k = Q k 32 572
; this is justified by assuming that the observed number of cases has a Poisson distribution [START_REF] Clayton | Models for temporal variation in cancer rates. i: age-period and age-cohort models[END_REF]. In fact the ARIC study yielded estimates of both MI and CHD, so we had to use formula (2).

We do not expect that the model for the Qk is well specified. We assume that the Qk are independent. Let θ 0 the value which maximizes the expectation of the pseudo-likelihood. Under mild regularity assumptions, an extended theory of Mestimators ( [START_REF] Van Der Vaart | Asymptotic statistics[END_REF][START_REF] Freedman | On the so-called huber sandwich estimator and robust standard errors[END_REF] ensures the consistency of θ for θ 0 . The variance of θ cannot be estimated directly by the inverse Hessian of the pseudo-loglikelihood (H -1 ), but by the sandwich estimators:

V ar θ = H -1 K k=1 U k U T k H -1 , (4)
where

U k = ∂L k ∂θ |
θ is the score for study k.

When a study gives several contributions we should evaluate the covariance matrix Ω k and use formula (2). In many cases however the correlations between contributions are very small so that they can be treated as independent: this is the case of incidence estimates for different age strata. In other cases the correlations can be difficult to estimate. We can then treat them as independent when writing the pseudo-likelihood but not in the sandwich estimator, an approach similar to the GEE [START_REF] Liang | Longitudinal data analysis using generalized linear models[END_REF]. In principle we could use this approach for all the studies but we would need a large number of studies in order to get reliable estimates of the variances. In this paper we will examine whether the contributions of the same study can be considered as independent and the L k , k = 1, . . . , K will be the independent contributions to the pseudo-loglikelihood.

A degradation model for myocardial infarction and CHD

The degradation model; modeling the drift

The atheromatous process A(t) can be modeled as a Brownian motion with drift defined by the stochastic differential equation: dA(t) = λdt + dB A (t), where B A is a Brownian motion and λ is the drift. The time parameter was taken as age (in years) minus 20 and we take A(0) = 0. A basic degradation model is that MI happens when the atheromatous process reaches a certain threshold η, so that the time T at which MI occurs is the first hitting time of A(t); for fixed η, T has an inverse Gaussian (IG) distribution with parameters (η/λ, η 2 ); its density is:

f (t) (λ,η) = η 2 2πt 3 1/2 exp    -λ 2 t - η λ 2 2t    1 {t≥0} .
What is interesting from an epidemiological point of view is to model the drift as a function of physiological conditions suspected to play a role in the atheromatous process. Here we will take into account four of them: obesity, represented by body mass index (BMI), lipid profile represented by low density lipid concentration (LDL), inflammation represented by C-reactive protein concentration (CRP) and blood pressure represented by systolic blood pressure (SBP); in this paper the risk factors are assumed constant. We assume a linear model for λ:

λ = λ 0 + β BMI BMI + β LDL LDL + β CRP CRP + β SBP SBP + ε λ , (5)
where λ 0 is a baseline drift and ε λ has a normal distribution with zero expectation and may represent unmeasured risk factors. We will in fact use transforms of BMI, LDL, CRP and SBP in order to get an approximately multinormal distribution (see section 4.2). Conditional on the values of these factors and the random error term ε λ , the time of occurrence of MI has an IG distribution.

However for our synthesis we wish to use results of studies which have not recorded all of these factors. Thus, we potentially need all the distributions of T obtained by conditioning or not on these factors, and for computing them, all the corresponding distributions of λ. We assume that the indicators (BMI, LDL, CRP, SBP) have a multinormal distribution; equation ( 5) gives us the conditional distribution of λ given (BMI, LDL, CRP, SBP). From this we can compute analytically the marginal distribution of λ (which is normal), as well as the distribution of λ conditional on BMI for instance (but marginal on LDL, CRP and SBP). Then the marginal distributions of T are inverse Gaussian Normal (IGN ) [START_REF] Whitmore | Normal-gamma mixtures of inverse Gaussian distributions[END_REF]. If the drift parameter λ has a N (m λ , s 2 λ ) distribution, then the hitting time T has the distribution IGN (m λ , s 2 λ , η), with density:

f (t) (m λ ,s 2 λ ,η) = η 2 2πt 3 (1 + s 2 λ t) 1/2 exp -(η -m λ t) 2 2t(1 + s 2 λ t) 1 {t≥0} (6)
This is an improper distribution in that P (T = ∞) > 0; this is not a problem in our model since not everybody develops a MI.

The model can be extended for defining two hitting times. It often happens that the progression of the atheromatous process first produces symptoms related to hypoxia (like angina pectoris) before the completion of MI; CHD includes these symptoms as well as MI. Thus two thresholds that we will denote η CHD and η MI can be defined and determine the distribution of the time of occurrence of CHD,

T CHD , and of MI, T MI . This is illustrated in Figure 1.

A graph of the causal influences between processes as suggested in Commenges and Gégout-Petit ( 2009) is represented in Figure 2.

Modeling drift and threshold

Another model arises if the threshold rather than the drift varies with the value of an explanatory variable. Several authors [START_REF] Aalen | Survival and event history analysis: a process point of view[END_REF][START_REF] Saebø | Disease resistance modelled as first-passage times of genetically dependent stochastic processes[END_REF][START_REF] Pennell | Bayesian random-effects threshold regression with application to survival data with nonproportional hazards[END_REF] have used a degradation model where both drift and starting point may depend on covariates; these authors have also introduced random effects for these parameters and Whitmore (1986) considered random drift and variance. All these authors parameterized the model in term of starting point with fixed threshold, while we assume that the starting point is 0 and our parameter is the threshold. This is just two different parameterizations but with different interpretations. Parameterizing in term of starting point is adapted if the subjects start at different levels of the pathological process. In our application it is natural to consider that all subjects start with a very low level of atherosclerosis, thus A(0) = 0, and to parameterize in term of threshold. For instance it may be asked whether SBP has a cumulative effect represented by an effect on the drift, or an t (age in years) instantaneous effect represented by an effect on the threshold, or both. SBP could favor the occurrence of MI by increasing the drift of the atheromatous process and also for a given state of the atheromatous process by lowering the threshold. The model for the threshold could be:

A(t)
η = η 0 + β ′ SBP SBP + ε η .
The graph of causal influences would then be as in Figure 3. Conditional on the explanatory variables and ε λ and ε η , the distribution of T is still IG. As before, for the synthesis analysis we need marginal distributions for which both λ and η are normal and may be correlated. We call the resulting distribution inverse Gaussian normal-normal (IGN N ), which to the best of our knowledge, has not been previously described.

It happens that the density of this distribution has an analytic form, given in Appendix A. This distribution may also be useful in Bayesian computations if we put normal priors on baseline drift and thresholds. The Physicians Health Study (PHS) is a randomized study of the prevention of cardiovascular diseases based on 22,071 subjects followed-up between 1982 and 1995. We use a case-control study nested in the PHS. Among these subjects, 543 men who developed a MI during the follow-up and 543 controls were chosen for studying the effect of CRP [START_REF] Ridker | Inflammation, aspirin, and the risk of cardiovascular disease in apparently healthy men[END_REF]. The study gives estimates of the relative risks for CRP strata with respect to a reference stratum (Table 3); it is not clear whether these are relative risks or odds-ratios but the two computations yield nearly the same result. The only information about age is the mean, which is 59. For simplifying the computations we attribute this age to all the subjects.

(θ) = -1 2 ( Qk -Q k (θ)) T Ω -1 k ( Qk -Q k (θ)), k = 1, 2, 3, 4. The
Q k (θ) = (Q k,1 (θ), Q k,2 (
The follow-up period was 13 years, thus we computed the risk for the age-period (59 -72). The relative risk of the CRP stratum ]0.055;0.115[ with respect to the stratum ]0;0.055] was computed as : [START_REF] Ridker | Inflammation, aspirin, and the risk of cardiovascular disease in apparently healthy men[END_REF]. Approximate computation of the correlations between the estimates of the three relative risks gives values lower than 0.01 and thus we consider that the PHS study gives three independent contributions L 5 (θ), L 6 (θ), L 7 (θ) of the form given by formula (1).

RR 1 CRP = F θ|CRP∈
The Health Professionals Follow-up Study (HPFS) is a prospective cohort study bearing on 42,351 subjects. [START_REF] Flint | Excess weight and the risk of incident coronary heart disease among men and women[END_REF] have estimated the risk of CHD on a 12-year period as a function of the BMI at inclusion; see Table 4. The absolute risk for the 12 year follow-up period for subjects presenting a BMI lower than 23 was computed assuming that all subjects were 53 years old at inclusion (the mean age at inclusion) as: The Framingham Heart Study (FHS) is a cohort study especially designed for cardiovascular epidemiology. [START_REF] Kannel | Risk stratification of obesity as a coronary risk factor[END_REF] have estimated the effect of BMI on CHD risk and give the number of CHD events over a 16 years period and the number of persons at risk. This allowed us to build Table 5. From the results of [START_REF] Wilson | Prediction of coronary heart disease using risk factor categories[END_REF] we can construct the same kind of tables for LDL and SBP:

Q 8 (θ) = F θ|BMI∈[0;23[ (65) -F θ|BMI∈[0;23[ (53) 1 -F θ|BMI∈[0;23[ (53) . ( 7
see Tables 6 and7. The estimators of absolute risks for the different strata of BMI, LDL and SBP were computed by the same type of formula as ( 7) presented for the PHS study. These estimators however can not be considered as independent because they rely on observations of the same subjects; thus we can treat them in the GEE spirit, considering that for the computation of the sandwich estimator there is only one independent contribution L 13 (θ) (the sum of the loglikelihood contributions for the 10 observations).

Construction of normal physiological indicators

The NHANES study is a repeated transversal surveys study recording a large number of physiological and behavioral factors on large representative samples of the American population. We used the 2007-2008 survey giving observations of BMI, CRP, LDL and SBP on n = 2871 men aged 20 or older. To make the computation of the marginal distributions of the event times feasible we need an underlying normal distribution for the physiological factors. Inspection of the empirical distribution of these variables reveals that they do not have normal distributions. In order to construct indicators with a distribution close to the normal we used Box-Cox transformations [START_REF] Box | An Analysis of Transformations[END_REF]) defined as X κ -1 κ κ = 0 and ln(X)

if κ = 0. The parameter κ can be estimated by maximum likelihood [START_REF] Velilla | A note on the multivariate Box-Cox transformation to normality[END_REF][START_REF] Yeo | A new family of power transformations to improve normality or symmetry[END_REF]. This is implemented in the R package CAR [START_REF] Fox | An R Companion to Applied Regression[END_REF]. We found the optimal values of κ to be -0.58, 0.41, -0.04, -1.27 for BMI, LDL, CRP and SPB respectively.

As shown in Figure 4 the optimal Box-Cox transformations yield new variables with a distribution much closer to the normal than the original ones. Of course marginal normal distributions do not mean that we have a multinormal distribution but we expect to get closer to it [START_REF] Kowalski | The performance of some rough tests for bivariate normality before and after coordinate transformations to normality[END_REF]. We applied a further transformation to standardize the indicators (zero mean and unit variance), so as to compare more easily the effects of these physiological factors.

From this study we could estimate the covariance (or correlation) matrix of our transformed variables:

corr(LDL, BMI, CRP, SBP) =        
1.00 0.07 0.09 -0.05 0.07 1.00 0.34 0.10 0.09 0.34 1.00 0.11 -0.05 0.10 0.11 1.00

       
We consider the empirical correlation as observations. To get approximate normal distribution we used Fisher transformation. For the correlation between LDL and BMI for instance, we considered Q14,1 = arctanh[corr(LDL, BMI)] as normal with mean arctanh(ρ LDL,BMI (θ)) (where ρ LDL,BMI (θ) is computed from the model) and variance 1 n-3 . However the Q14,j , j = 1, . . . , 6 may be correlated in a way which is difficult to analyze. Thus we consider that the NHANES study gives just one independent contribution to the pseudo-loglikelihood

L 14 = 6 j=1 n-3 2 ( Q14,j - Q 14,j (θ)) 2 .

Results

Finally, the global pseudo-loglikelihood was the sum of 14 independent contributions : L(θ) = 14 k=1 L k (θ). The pseudo-loglikelihood was computed using a Fortran program and was maximized using a Marquardt algorithm [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF] (also programmed in Fortran). The algorithm converged in 29 iterations (it was verified that different starting points led to the same convergence point). Table 8 gives the fit of the model for the information of the studies about MI and CHD. We do not show the fits for the correlations of the physiological factors in the NHANES study: they were very good and the contribution was L 14 = -0.003.

For the other quantities, we consider that the fit is rather satisfactory given the heterogeneity of the studies. For instance the risk observed in the HPFS study for the 27-30 BMI stratum is around 8%; it is to be compared with the risk computed with the model of developing CHD between age 53 and 65 for this BMI stratum, which is around 10%: we get the good order of magnitude. Of course this is not a good fit in view of the small standard deviation attributed to the observed risk; if the model was well specified -L k would have a chi-squared distribution with one degree of freedom: the value of 17.5 is much too large. We know from the beginnning that our model cannot be well specified and this is why we used the term "pseudo-loglikelihood" rather than "loglikelihood".

Table 9 gives the parameter estimates together with sandwich estimates of their standard deviations; the estimate of the standard deviation of ε λ was very close to zero and is not shown in the table. We note that the four regression parameters are positive, meaning that increasing the level of these factors increases the risk, which is consistent. Since the values are greater than 2 times their standard deviations, they are significant. There is a positive baseline drift and, as expected η CHD < η MI .

From the sandwich covariance matrix we can compute the standard deviation of ηCHD -ηMI . We found 0.051 so that ηCHD -ηMI is significantly different from zero. We see that the latter profile has very low risk to develop MI during lifespan while profile (i) has high risk, for instance a probability of 0.80 to develop MI before 80.

Finally we have fitted the model described in section 3.2, where SBP could modify both thresholds (we took the same regression coefficient for both). The pseudo-loglikelihood for this model was -212.56. This is not significantly better than that of the basic model (without influence on the thresholds) which was

Conclusion

We have shown that a complex degradation model could be used for evidence synthesis of disparate results coming from different studies. The results are consistent in that we found a positive baseline drift, positive effects of the risk factors, positive values of the thresholds with higher threshold for MI than for CHD. We could have restricted the parameter space to impose the constraint η MI > η CHD . However at the maximum of the likelihood, either we would find ηMI > ηCHD or ηMI = ηCHD .

In the former case the constraint is inactive so that the result is the same whether we impose the constraint or not; in the latter case the constraint is active and this would throw suspicion on the model adequacy because without the constraint we would have ηMI < ηCHD .

We are limited in the epidemiological and clinical interpretation of these results because of the small number of studies that we have included. We have made the assumption of a linear relationship between λ and the risk factors. This is necessary to keep the computations relatively simple. A non-linear relationship would entail a non-normal joint distribution and we would lose the analytical marginal or partially conditional distributions. This is feasible but at the price of heavier computations. Also we have not included tobacco consumption; this factor raises a methodological problem due to the difficulty of finding a normalizing distribution. As smokers and non-smokers form a dichotomy, the situation calls for a two-component mixture of IG distributions rather than a normal mixture (IGN)

over smoking status. Another potential problem is that we implicitly assume that death acts as an independent censoring. There may be selection effects due to deaths by diseases which share some common physiological process: for instance oxidative stress may play a role in both CHD and cancer.

The synthesis analysis can be viewed as an extension of meta-analysis. Indeed fixed-effect meta-analysis is a particular case when the Q k (θ) are the same parameter for all k. Often meta-analyses put a random effect per study. This would also be possible in the synthesis approach and would certainly improve the fit, at the price of more computations. This is a "proof-of-concept" paper, which shows that such an approach is very promising and could be developed in several directions. It could be applied to a large number of studies of MI or CHD for producing reliable results, allowing to understand and to link the results of the literature on the subject. The evidence synthesis approach would become even more important if we moved toward lifecourse epidemiology for which information about the evolution of physiological factors would be incorporated, as suggested in [START_REF] Commenges | The stochastic system approach to causality with a view toward lifecourse epidemiology[END_REF]. This approach could also be applied to other topics as well.

Appendix A: The IGN N distribution

The inverse Gaussian normal-normal (IGN N ) distribution is an IG where parameters λ and η have a bivariate normal distribution with marginals N (m λ , s 2 λ ) and N (m η , s 2 η ) and correlation coefficient ρ. The probability density function of this distribution, although more complicated than that of IGN , has an analytic form that we have derived with the help of Maple TM :

f (t) (m λ ,m η ,s 2 λ ,σ 2 η ,ρ) = - s η 1 + ts 2 λ (1 -ρ 2 ) π √ t

Fig. 1 Fig. 2 Fig. 3

 123 Fig.1Example of the trajectory of the atheromatous process A(t). η CHD and η M I are the two thresholds; here, the subject has CHD symptoms at 62 and develops an infarctus at 68.

4

  Application to five studies 4.1 Studies giving information on CHD and MI risks We have used five large studies. The National Health and Nutrition Examination Survey (NHANES) has been used to estimate the joint distribution of the physiological indicators in the population (see section 4.2). Four studies give information

  θ)) T were computed by equation (3), using the IGN (or IGN N ) distributions of T CHD and T MI . For instance in formula (3) we have to compute F θ (65), the value at 65 of the marginal cumulative distribution of the time of occurrence of MI.

  ]0.055;0.115[ (72) -F θ|CRP∈]0.055;0.115[ (59) 1 -F θ|CRP∈[0;0.055] (59) F θ|CRP∈[0;0.055] (72) -F θ|CRP∈[0;0.055] (59) 1 -F θ|CRP∈]0.055;0.115[ (59) , where F θ|CRP∈]a;b[ (.) is the cumulative distribution function of the T MI conditional on CRP ∈]a; b[. To avoid the numerical integration required for the computation of F θ|CRP∈[a;b] (t), we approximated it by F θ|CRP=c (t) where c is the median value of CRP in the sub-sample [a, b] of the NHANES study. The standard deviations σ k were computed from the confidence intervals given in the publication. In Ridker et al (1997) the confidence intervals have been derived from a normal approximation of the distribution of the regression coefficient. Thus the Qk were the estimated values of this coefficients and the standard deviations σ k were deduced from the

  ) Here, F θ|BMI∈[a;b[ (.) is the cumulative distribution function of T CHD given BMI ∈ [a; b[. The variance was estimated by Qk a k . Similar formulas were applied for the other BMI strata. The estimates for the five BMI strata are approximately independent, so the HPFS study yielded five independent contributions L k (θ), k = 8, 12.

Fig. 4

 4 Fig. 4 Quantile-quantile plot of raw and transformed physiological variables. Left: original indicators; right: after Box-Cox transformation.

Figure 5

 5 Figure 5 displays the cumulative distribution of the times of occurrence of MI (probability to have MI before age t) for three profiles: (i) mean level for all risk factors: λ = λ 0 = 0.99, (ii) above the mean level by one standard deviation for all factors: λ = λ 0 + β LDL + β BMI + β CRP + β SBP = 0.354, (iii) below the mean level by one standard deviation for all factors: λ = λ 0 -β LDL -β BMI -β CRP -β SBP = -0.155.

Fig. 5

 5 Fig. 5 Probability of developing MI as a function of age for three risk factors profiles: λ = 0.099: all risk factors at their mean value; λ = -0.155: all risk factors at one standard deviation below the mean; λ = 0.354: all risk factors at one standard deviation above the mean.

Table 1

 1 The studies used for the synthesis: NHANES: National Health and Nutrition Examination Survey; ARIC: Atherosclerosis Risk In Communities; PHS: Physicians Health Study; HPFS: Health Professionals Follow-up Study; FHS: Framingham Heart Study; # Contributions gives the number of independent contributions to the pseudo-loglikelihood, and in parenthesis the number of elementary contributions.

	Study	Nature of Q k	Event	Risk factors	# Contributions
	NHANES	Correlations	-	BMI, LDL, CRP, SBP	1 (6)
	ARIC	Incidences	MI and CHD	-	4 (8)
	PHS	Relative risks	MI	CRP	3 (3)
	HPFS	Absolute risks	CHD	BMI	5 (5)
	FHS	Absolute risks	CHD	BMI, LDL, SBP	1 (10)

on CHD and MI risks. Table

1

gives the list of the studies used together with the type of information they bring.

The ARIC (Atherosclerosis Risk In Communities) surveillance study gives estimates of the incidence of both MI and CHD for four age strata; see Table

2

, taken from Tables 4.9 and 4.10 of National Heart

Lung and Blood Institute (2006)

. The estimates of the MI incidences of different age strata can be considered as approximately independent; this is also the case for CHD incidences. We also make this approximation for incidence of MI and CHD for different age strata. This approximation however is not tenable for incidence estimates of MI and CHD for the same age stratum because the CHD cases observed include the MI cases. For stratum k we have Y k,CHD = Y k,MI + Y k, MI , where Y k,CHD , Y k,MI and Y k, MI are the observed numbers of CHD, MI and non-MI cases respectively. The incidences are estimated by respectively, where a k1 and a k2 are the respective numbers of person-years. Assuming independence between Y k,MI and Y k, MI we have cov(Y k,CHD , Y k,MI ) = var(Y k,MI ). Assuming Poisson distributions for these numbers,

Table 2

 2 Incidence of CHD and MI in men in USA, ARIC Surveillance, 1987-2001.

	Age	CHD sample	CHD Incidence	MI sample	MI Incidence
	35-44	56, 457	1.19‰	61, 554	1.05‰
	45-54	42, 257	3.55‰	45, 831	3.22‰
	55-64	29, 606	7.32‰	32, 572	6.26‰
	65-74	20, 796	12.20‰	23, 049	9.64‰

Table 3

 3 Relative risk of MI in men according to the quartile of CRP concentration,PHS, 1982PHS, -1995. . Table4CHD risk in men in USA according toBMI strata, HPFS, 1988-2004.

		Quartile of CRP concentration (mg • dL -1 )
	≤ 0.055	0.056 -0.115	0.116 -0.210	≥ 0.211
	Relative risk	1.0	1.7			2.6	2.9
	95% CI	-	1.1-2.9			1.6-4.3	1.8-4.6
						BMI
			< 23	23 -24.9	25 -26.9	27 -29.9	≥ 30
	Number of events		349	625		821	643	333
	Number of subjects at risk	7,669	12,104	11 466	7,712	3,400
	confidence intervals of the relative risks as σ k = (log( RR k CRP ) -log(RR k inf ))/1.96 where RR k inf is the lower bound of the confidence interval for RR k CRP given in

Table 5

 5 Risk of CHD in US men according toBMI strata, FHS, 1971[START_REF] Clayton | Models for temporal variation in cancer rates. i: age-period and age-cohort models[END_REF]

			BMI
		< 25	25 -29.9	≥ 30
	Number of events	99	188	72
	Number of subjects at risk	772	1,306	375
	Table 6 CHD risk in US men according to LDL strata, FHS, 1971-1983.
		LDL (mg • dL -1 )
		< 130	130 -159	≥ 160
	Number of events	104	124	155
	Number of subjects at risk	929	866	719

Table 7

 7 CHD risk in US men according toSBP strata, FHS, 1971-1983.

			SBP (mm Hg )
		< 130	130 -139	140 -159	≥ 160
	Number of events	110	77	115	81
	Number of subjects at risk	1127	526	556	303

Table 8

 8 Fit of the information from the four studies. RR 1 CRP , RR 2 CRP , RR 3 CRP are the relative risks for the CRP strata 0.56-1.14, 1.15-2.10, ≥ 2.11 respectively, versus the reference stratum CRP ≤ 0.55, and AR stands for absolute risk. Last column: -L k , where L k , k = 1, . . . 13 are the first 13 independent contributions to the pseudo-loglikelihood. For ARIC the contributions comes from the bivariate observations of incidences of CHD and MI in age strata. For FHS the global contribution is given. The 14th contribution brought by the correlations is not shown.

	Study	Event	Nature of Q	Q( θ)	Q	σ	-L k
	ARIC	CHD	Incidence 35-44	0.0013	0.0012	0.00015
	ARIC	MI	Incidence 35-44	0.0011	0.0010	0.00013	3.5
	ARIC	CHD	Incidence 45-54	0.0045	0.0035	0.00029
	ARIC	MI	Incidence 45-54	0.0041	0.0032	0.00026	6.2
	ARIC	CHD	Incidence 55-64	0.0070	0.0073	0.00050
	ARIC	MI	Incidence 55-64	0.0063	0.0063	0.00044	14.9
	ARIC	CHD	Incidence 65-74	0.0080	0.0122	0.00077
	ARIC	MI	Incidence 65-74	0.0077	0.0096	0.00065	32.6
	PHS	MI	log RR 1 CRP	0.47	0.53	0.27	0.0
	PHS	MI	log RR 2 CRP	0.77	0.96	0.26	0.3
	PHS	MI	log RR 3 CRP	1.17	1.06	0.24	0.1
	HPFS	CHD	ARBMI(< 23)	0.034	0.046	0.0024	10.6
	HPFS	CHD	ARBMI(23-25)	0.056	0.052	0.0021	2.6
	HPFS	CHD	ARBMI(25-27)	0.077	0.072	0.0025	2.2
	HPFS	CHD	ARBMI(27-30)	0.103	0.083	0.0033	17.5
	HPFS	CHD	ARBMI(≥ 30)	0.159	0.189	0.0074	7.9
	FHS	CHD	ARBMI(< 25)	0.042	0.128	0.013
	FHS	CHD	ARBMI(25-30)	0.080	0.144	0.010
	FHS	CHD	ARBMI(≥ 30)	0.136	0.192	0.023
	FHS	CHD	ARLDL(< 130)	0.037	0.112	0.011
	FHS	CHD	ARLDL(130-159)	0.100	0.143	0.013
	FHS	CHD	ARLDL(≥ 160)	0.174	0.216	0.017
	FHS	CHD	ARSBP(< 130)	0.035	0.098	0.0093
	FHS	CHD	ARSBP(130-139)	0.087	0.146	0.0167
	FHS	CHD	ARSBP(140-160)	0.128	0.207	0.0172
	FHS	CHD	ARSBP(≥ 160)	0.212	0.267	0.0297	114.4

Table 9

 9 Parameter estimates: ρ X,Y is the correlation between X and Y , the β's are the regression coefficients of the standardized transformed physiological factors, λ 0 is the baseline drift, the η's are the thresholds.

		Estimate	Sandwich SD
	ρ BMI,LDL	0.070	0.0015
	ρ BMI,CRP	0.343	0.0003
	ρ BMI,SBP	0.096	0.0014
	ρ LDL,CRP	0.093	0.0007
	ρ LDL,SBP	-0.045	0.0024
	ρ CRP,SBP	0.107	0.0006
	β BMI	0.053	0.014
	β LDL	0.092	0.015
	β CRP	0.024	0.006
	β SBP	0.087	0.015
	λ 0	0.099	0.008
	η MI	17.50	0.81
	η CHD	17.06	0.79

  t 2 s 2 λt(2s λ ρs η -1) + s 2 η -2s λ m η ρs η m λ t + m 2 η 2 (1 + ts 2 λ (1ρ 2 )) s 2 η + (s λ ρs η m λs 2 λ m η )t + s λ ρs η m ηs 2 η m λm η 2t 2 s 2 λ + (2 -4s λ ρs η )t + 2s 2 η 2 √ π t 2 s 2 λt(2s λ ρs η -1) + s 2 η λ ρs η m λs 2 λ m η )t + s λ ρs η m ηs 2 η m λm η s η 1 + ts 2 λ (1ρ 2 ) 2t 2 s 2 λ + (2 -4s λ ρs η )t + 2s 2 η

		e	-	s 2 λ m 2 η + s 2 η m 2 λ 2
				
	e			
	where erf is the error function (erf(x) =	2 √ π

-(tm λm η ) 2 2 t 2 s 2 λt(2s λ ρs η -1) + s 2 η erf   √ t (s x 0 t -2 dt).


N Fabbri
email: fabbri@lens.unifi.it

S D Huber

D Clément

‡ L Fallani

C Fort

M Inguscio

E Altman

H P Büchler

Phys. Rev. B G Blatter

] G D Mahan

Quasiparticle dynamics in a Bose insulator probed by inter-band Bragg spectroscopy

Keywords: numbers: 37.10.Jk, 67.85.Hj, 67.85.De

We investigate experimentally and theoretically the dynamical properties of a Mott insulator in decoupled one-dimensional chains. Using a theoretical analysis of the Bragg excitation scheme we show that the spectrum of inter-band transitions holds information on the single-particle Green's function of the insulator. In particular the existence of particle-hole coherence due to quantum fluctuations in the Mott state is clearly seen in the Bragg spectra and quantified. Finally we propose a scheme to directly measure the full, momentum resolved spectral function as obtained in angle-resolved photoemission spectroscopy of solids.

The observation of the superfluid (SF) to Mott insulator (MI) transition of bosons in optical lattices [1] has received considerable attention as a paradigmatic example of a quantum phase-transition driven by interactions. The properties of lattice bosons in this strongly correlated regime have been probed using several methods [2][3][4][5][6][7][8][9][10]. For example, time-of-flight experiments were used to study the development of spatial first-order coherence over increasing length-scales inside the Mott state upon approaching the transition to the SF state [3]. In a quantum system, the emergence of such spatial correlations must go hand in hand with increasing temporal correlations. Near the quantum critical point the precise relation between the two is determined by the dynamical critical exponent of the transition (see e.g. Ref. [START_REF] Sachdev | Quantum Phase Transitions[END_REF]). Away from the transition, where critical properties are not yet apparent, the temporal first-order coherence lends insight on the nature of the quasi-particle excitations of the strongly correlated state.

So far however, most of the dynamical experiments in the Mott regime using schemes of lattice modulation and Bragg spectroscopy have focused on excitation frequencies matching transition within the lowest-energy Bloch band [2,7]. In this case, the external perturbation is coupled to density fluctuations and in the linear response regime the absorption spectrum is directly related to the dynamical structure factor S(q, ω) of collective excitations, or particle-hole spectra [START_REF] Iucci | [END_REF][13][14].

In this Letter, we show how inter-band Bragg spectroscopy [7,15] supplemented by a theoretical model of the Mott insulator can be used to extract properties of single-particle excitations in the many-body state. We then suggest a refined approach for directly measuring the single-particle Green's function in a modelindependent way.

The MI state is realized in decoupled one-dimensional (1D) chains of interacting bosons, as represented in Fig. 1. We excite the system with two simultaneous laser pulses (Bragg beams) which induce an energy transfer ω = (ω 1 -ω 2 ) (ω 1,2 being the laser beams frequen- cies) and a momentum transfer q = q 1q 2 = qe x along the axis of the 1D chains (q 1 and q 2 being the wave-vectors of the Bragg photons). We measure the energy absorption spectrum D(ω) at a fixed momentum transfer q π/a, where a is the periodicity of the lattice along the chains. We show how, with the precise knowledge of the particle dispersion in the high band [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF][START_REF] Zwerger | [END_REF], it is possible to obtain information about quasiparticle structure and dynamics in the MI. Moreover, a refined scheme would give access to a momentum resolved absorption rate D(k, ω), which we show is directly related to the single-particle spectral function A(k, ω) in the lowest-energy band. From the latter, one obtains the Green's function (GF) of the MI, lending information on both spatial and temporal coherence of quasiparticles. Such a measurement is analogous to angle-resolved photoemission spectroscopy used in solid-state [18], recently extended to ultracold gases through Raman [19,20] and rf spectroscopy [21].

We start by briefly describing the experimental setup used to obtain the inter band spectra (details can be found in [7,22]). We load a Bose-Einstein condensate of 87 Rb in a 3D optical lattice at the wavelength λ L = 830 nm. The amplitudes V i of the lattices along each axis i = x, y, z are expressed in units of the recoil energy

E R = h 2 /(2mλ 2 L ), V i = s i E R
, where m denotes the mass of 87 Rb. The optical lattices are ramped up to their final values s i with an exponential ramp of duration 140 ms and time constant 30 ms. Two lattice amplitudes (s y = s z = 35) are fixed to create an array of 1D chains. The amplitude of the third lattice s x is varied to tune the ratio between the on-site interaction energy U and the tunneling amplitude J 1 between Wannier states of the lowest Bloch band in each 1D MI chain from U/2J 1 ≃ 7 to U/2J 1 ≃ 42.

The Bragg beams are derived from a laser at 780 nm, detuned by ∼ 200 GHz from the D 2 line of 87 Rb. In this work we fix q = 0.96π/a and we measure the amount of excitations induced by the Bragg beams as a function of their frequency difference ω. The measured quantity D(ω) is the mean square width of the zero-momentum peak in a phase-coherent lattice-gas obtained after lowering the 3D optical lattices. Timing and details of the experimental procedure can be found in [22], in particular the way D(ω) is re-scaled with the parameters of the Bragg beams to allow a relative comparison of the different spectra. In [22], we also verified that D(ω) is proportional to the energy transferred to the gas so that it can be written as D(ω) = C ωS(q, ω) where C is a constant independent of the lattice strength [23].

In Fig. 2(a) we show an example of the spectra obtained for lattice of amplitude s x = 9 at frequencies resonant with transitions to the second and third Bloch bands. The total spectral weight W = dωD(ω) of transitions to the third band as a function of s x is shown in Fig. 2(b). The suppression of spectral weight is due to reduction of the matrix element, or Frank-Condon overlap, between wavefunctions of the two bands with increasing lattice strength. The distribution of spectral weight within each band is primarily determined by two factors (i) the density of final states (DOS) in that band and (ii) the particle hole coherence in the Mott ground state driven by quantum fluctuations about the classical state with precisely n (integer) particles on each site. The peaks seen at the band edges, e.g., are the result of the divergent DOS there. Furthermore we quantify the asymmetry of the spectra about the band centers through the skewness (third moment) [START_REF]The skewness is defined as γ =[END_REF] of S(q, ω). This is shown in Fig. 2(c) and (d) for the second and third band as a function of s x . A positive skewness corresponds to an imbalance towards lower energies. Using a theoretical model we shall relate the skewness and its variation with band index and lattice amplitude, to the particle-hole coherence in the Mott state. The Bragg perturbation couples to the particle-density and, in the linear response regime, the excitation rate D(ω) is directly related to the dynamic structure factor

S(q, ω) = m | m|ρ q |0 | 2 δ(ω -ω m0 ) = Im[Π(q, ω)], (1)
where ρ q is the density operator at momentum q [START_REF] Stenger | [END_REF]. The sum runs over excited states and ω m0 is the associated excitation energy.

The response function Π(q, ω) is represented graphically by the bubble diagram in Fig 3(a): (i) the vertex describes the coupling of the Bragg beams to the density operator ρ q ; (ii) the full line is the GF of the hole produced in the lowest band by the Bragg excitation; (iii) the dashed line represents the GF of the particle excited to the nth band; (iv) the filled area is the T -matrix for scattering of the upper-band particle with the hole in the lowest band. We explain how we evaluate each part of the diagram, showing that for excitations to the higher bands the contribution of the final-state interaction can be neglected, hence the process can be well described by the bare bubble shown in Fig. 3(c). Since without this interaction the propagation of the upper-band particle can be calculated as the one of a free particle, the experiment effectively probes the remaining unknown which is related to the single-particle GF.

Γ ω, q ω, q (a) (b) (c) ω -, q -k , k ω
Let us first consider the edge vertices which correspond to the matrix element for the excitation of a particle from the lowest to the upper band by the density operator. To find this matrix element we express the field operators in terms of creation operators in Wannier states of the lattice sites, so that ρ(x) = ij,nm w * mj (x)w ni (x)b † ni b mj , where n, m are band indices, i, j site indices and w ni (x) the respective Wannier functions. From now on we assume that the lattice is sufficiently deep that the overlap between Wannier functions of neighboring sites can be neglected in the density operator. In addition, we focus on the component of the density operator that induces transitions from the lowest to the nth band:

ρ n (q) = ρF 1n (q) i e -iqRi b † ni b 1i , (2)
where R i is the position at site i and ρ the filling of the MI. F 1n (q) is the matrix element

F 1n (q) = dx w 1i (x)e -iqx w ni (x) ≈ (-iql 0 ) n √ 2 n n! e -(ql 0 ) 2 4 .
(3) To get a simple expression for the dependence of F 1n (q) on the lattice strength we approximated the Wannier states by those of a harmonic well with oscillator length l 0 = λ/(2πs

1/4
x ); however for the comparison with experimental results we use the exact Wannier functions [START_REF] Zwerger | [END_REF].

To describe the hole propagation in the MI we use the generalized Bogoliubov theory [14,26], which accounts for quadratic quantum fluctuations (i.e. virtual particle and holes) about the classical (Gutzwiller) groundstate. The particle operator in the lowest band can be represented in terms of the Bogoliubov quasi-particle and quasi-hole excitations as b k = f (k)(β † h,k + β p,k ), which are a combination of an added particle (doublon) and removed hole on the Mott background, e.g.,

β † p,k = u k p † k + v k h -k and f (k) = (u k -v k ) 2 .
The ground state is the vacuum of the operators β p/h,k . Therefore f (k) = b † k b k is simply the momentum distribution. Accordingly the single particle Green's function in the lowest Bloch band is given by:

G 1 (k, iω) = f (k) iω -ω p (k) + f (k) -iω -ω h (k) ≡ G p + G h . (4)
ω p (k) (resp. ω h (k)) denotes the dispersion relation of a particle (resp. a hole) in the lowest Bloch band. The momentum distribution f (k) stems from quasiparticle coherence factors and within the Bogoliubov theory it is

f (k) = 1 1 -J1 Jc cos k . ( 5
)
Here J c is the critical hopping strength at the Mott transition. We shall restrict ourselves to using the Bogoliubov theory deep in the Mott insulator where it provides a good approximation of the single particle Green's function [START_REF]On approaching the transition the Bogoliubov expression of the single particle Green's function suffers from being a non-conserving approximation[END_REF]. In this regime f (k

) ≈ 1 + 1 2 (J 1 /J c ) cos k to leading order in J 1 /J c .
The GF of a single particle in the nth upper band is taken to be that of a free particle with appropriate band dispersion ωn (k). We take into account a slightly renormalized dispersion due to interaction with the background of filled sites of the MI [28].

The interaction between the particle in the upper band and the hole in the lowest band is included in the full Tmatrix (filled box in Fig. 3(a)). In general this leads to a complicated sum of diagrams including all possible sequences of multiple collisions through the interaction term

U 1n b † ni b ni (p † i p i -h † i h i ).
Here we represented the interaction in terms of actual particles and holes in the classical ground state, where U 1n is the interaction matrix element between Wannier states of the lowest and the nth band. The interaction looks more complicated when expressed in terms of the Bogoliubov quasi-particles and quasi-holes. However, the sum simplifies in the strong lattice limit when v k ≪ u k . Then, to leading order in v k we can include only the ladder diagrams shown in Fig. 3(b), which are easily summed up as a geometric series [28]. The result of the interaction, treated by the ladder summation, is to induce a bound state between the upper band particle and the hole in the MI [28]. For the experimental parameters the weight carried by this bound state (< 1%) is too small to affect the measurements. We conclude that to an excellent approximation we can use the bare bubble diagram shown in Fig. 3(c) to compute the structure factor.

The structure factor computed from the bare bubble diagram is given by

S(q, ω) = ρ2 |F 1n | 2 dk 2π , A h (k -q, ω -ωn (k)) (6)
where

A h (k, ω) = -π -1 Im G h (k, ω + i0 + )
is the hole spectral function in the Mott insulator. In particular, the spectral function obtained from the generalized Bogoliubov theory of the Mott insulator is

A h (k, ω) = f (k)δ[ω -ω h (k)].
This is the formula we use to compare with the experiment. However to clearly reveal the two important factors in the spectra it is worth making another simplification. In the regime of interest, of strong optical lattice, the bandwidth of the hole in the lowest band can be neglected compared to the dispersion of a particle in the excited band and we can take ω h (k) = ω 0 . This leads to

S(q, ω) ≈ ρ2 |F 1n (q)| 2 ρ n (ω -ω 0 )f (k n (ω -ω 0 ) -q), (7)
ρ n (ω) being the single particle DOS in the nth band and k n (ω) the inverse function of the dispersion in that band ω n (k).

In the limit of infinitely deep lattice f (k) → 1 and the observed line-shape is determined solely by the single particle DOS. With reduction of the lattice amplitude (increased hopping) f (k) becomes more strongly peaked near k = 0 and therefore contributes more significantly to the line-shape. Specifically it gives increased weight to frequencies resonant with transitions that create a Mott hole near k = 0 and an excited nth band particle with quasi-momentum q. If we take q ≈ π/a, as in the experiment, this effect enhances the weight of transitions that create a higher band particle near the Brillouine zone edge. From the band structure shown in Fig. 1 it is clear that in this way the momentum distribution skews the spectra of the second band toward lower energies and those of the third band toward higher energies. It should be noted however that the single-particle DOS is itself not symmetric about the band centers. In particular for both the second and third band the peak of the DOS is skewed toward lower energies (positive skewness). Therefore, the effect of coherence in the Mott insulator is to increase the skewness of the second band and decrease it in the third band spectra.

Using Eq. ( 6) we calculate the experimental observable D(ω) = CωS(q, ω). The proportionality constant C is fixed by matching the integrated spectral weight of excitations to the third band W for a single value of the lattice strength (s x = 10). We use the same constant to compute the spectra for all other lattice amplitudes and for all the bands. In addition, we broaden the deltafunction in Eq. ( 6) to effectively account for the trap confining potential.

The spectrum obtained in this way is presented and compared to the experimental results in Fig. 2(a) for s x = 9. Note that there are no free parameters except the overall proportionality constant C which was calibrated once. We attribute the relative shift of the spectra to the systematic uncertainty in the actual lattice amplitude in the experiment [29]. The calculated total weight of absorption W , shown in Fig. 2(b), is in good agreement with the experimental data. The reduction of W with increasing lattice strength is due to suppression of the Frank-Condon factor F 1n .

As discussed above, the skewness of the structure factor relative to that of the pure single particle density of states is a direct measure of the quasi-particle coherence factor in the Mott insulator. Fig. 2 (c) and (d) compares the measured skewness [START_REF]The skewness is defined as γ =[END_REF] of S(q, ω) in the second and third bands to that calculated from the theoretical spectra (thick lines). Both are compared to the skewness of the single particle DOS in the corresponding bands. As anticipated, the actual skewness is consistently higher than the pure DOS effect in the second band and lower than the DOS effect in the third band and this effect is observed in a systematic way for different lattice amplitude. This is a clear indication of coherence effects inside the MI.

To conclude, we have shown that the inter-band Bragg absorption spectrum, in the linear response regime, gives information on particle Green's function. In particular we have quantified the single hole coherence through analysis of the asymmetry of the spectra.

It is important to note that the structure is related, through Eq. ( 6) to a rather complicated weighted sum over the hole spectral function and is not proportional to the spectral function itself. For this reason input from a theoretical model of the Green's function in the Mott insulator was needed to extract a measure of the hole coherence. For this we used a Bogoliubov-like theory that takes into account zero-point fluctuations around the mean field Mott wave function. The non-trivial hole coherence stems from these zero point fluctuations.

It would be interesting to directly measure the single hole spectral function and thereby obtain quasi-particle energies, coherence factors and decay times in a model independent way. We suggest that in principle this can be done using a band mapping technique [START_REF] Greiner | [END_REF]. By counting how many particles are excited to a particular k state in the upper band we would eliminate the k integral in Eq. ( 6). Then the response function corresponding to the excitation rate per final momentum k would be

S 1n (q, k, ω) = ρ2 |F 1n (q)| 2 A h (k -q, w -ωn (k)). (8)
This theoretical description of the measurement process is identical to the description of angle-resolved photoemission spectroscopy [31,[START_REF] Mahan | Many-Particle Physics[END_REF], a method that is currently used extensively to measure the electronic spectral function of interesting materials [18]. Our proposed scheme could also be implemented on more complex many-body ground-states in a lattice. is enhanced by a few percent. This comes about as the combined effective potential is less confining than the pure sinusoidal optical lattice, see Fig. 1(a).

uv -u 2 u 2 (g) (f) (d) (e) (c) (b) w 2 1 ( x ) V latt V e (
By including the Hartree shift we obtain a new effective Hamiltonian for the particles in the upper bands:

H ef f = n>1,k e n (k)b † nk b nk (6)
Note that the Bose operators in the upper band are defined using the Wannier states in the effective potential V ef f :

ψ † (x) = n>1 wni (x)b † ni . (7)
Quasi-particle scattering. We now turn to the residual interactions which include scattering of particles in the upper band on particles and holes in the lowest band, which is represented by

H res = n>1 U 1n (p † i p i -h † i h i )b † ni b ni , (8)
with

U 1n = 2π ρa s 2 m dx w 2 1,i (x) w2 n,i (x). (9)
Since we are interested in one particular Bloch band, we now suppress the index n from the Bose operators of the upper band and we write V = U 1n . We can now write the inter-band interaction using the Bogoliubov quasiparticles of the lowest band defined in 5. This gives

H res = V N 2
k1,k2,q

s k1+q/2,k1 b † k2-q/2 b k2 (10)
× (β † p,k1+q/2 β p,k1 -β † h,-k1-q/2 β h,-k1 ),

+ V N 2
k1,k2,q r k1+q/2,k1 b † k2-q/2 b k2 × (β † p,k1+q/2 β † h,-k1 + β p,k1+q/2 β h,-k1 ).

where the interaction coefficients are given by s p,q = u p u q -v p v q J1/U→0 = 1, (11) r p,q = u p v q -u q v p J1/U→0 = 0. [START_REF] Iucci | [END_REF]

FIG. 1 :

 1 FIG.1:(color online) An array of 1D gases created by a 2D optical lattice is driven into the MI by a third OL in the Ox direction. The energy band structure in the 1D lattice is depicted on the left with red solid lines. The lowest band corresponds to particles in the MI and above it are the higher single-particle energy bands. Two laser beams (Bragg beams, green) excite the initial MI by transferring a particle to a high-energy band and leaving a hole in the many-body state.

FIG. 2 :

 2 FIG. 2: (color online) Inter-band Bragg spectra.(a) Energy absorption rate D(ω) over the energy range of the second and third Bloch band for a Bragg momentum transfer q = 0.96π/a and for a lattice depth sx = 9. The (blue and green) dots are the experimental data with error bars indicating statistical uncertainties after averaging over 4 to 5 experimental acquisitions, the black lines are the theoretical predictions for D(ω). (b) Integrated spectral weight W in the third band. (c) and (d) Skewness (third moment) of S(q, ω) in the second (c) and third (d) band. The measured values are compared to the theoretical prediction including particle-hole coherence (thick lines) and to the contribution of the singleparticle density of states alone (thin lines). The shaded area indicates the systematic uncertainty in determination of sx in the experiment.

FIG. 3 :

 3 FIG. 3: Response functions. (a) Diagram describing the Bragg excitation in the linear response. The black squares describe the coupling of the Bragg beams to the density fluctuations. The full (dashed) lines denote the GF of a hole in the MI (of an upper band particle). The gray area Γ stands for the final-state interactions between the upper band particle and the MI hole. (b) In a T -matrix approximation the final state interaction gives rise to a ladder diagram where the wiggly lines denote the interaction between upper band particle and the hole. (c) In the absence of final state interactions the bare bubble describes the experimental response.

  uv

FIG. 1 .

 1 FIG. 1. (color online) Final state interactions. (a) The effective potential for higher Bloch bands due to the underlying Mott insulator. (b)-(e) Scattering processes due to the inter-band interactions.A full line denotes a particle in the excited Bloch band, a black (gray) dashed line a quasi-hole (quasi-particle) in the Mott insulator. Note that quasi-holes and excitedband particles attract each other. (f) Polarization bubble corresponding to Sn(q, ω). (g) Interaction corrections to Sn(q, ω) leading to a bound state, a Mahan-exciton.

The different scattering processes included in the above equations are illustrated in Fig 1

We note that processes (b) and (c) are proportional to u 2 whereas (d) and (e) scale as uv and are therefore strongly suppressed deep in the Mott phase. Furthermore process (c) vanishes identically because there is only a quasi-hole in the lower band. So, only process (b) contributes significantly to the final state interaction. The contribution of this scattering process to Π(q, w) can be computed exactly by resummation of the ladder diagrams shown in Fig. 1(g).

The resummation reveals the presence of a bound state of a quasi-hole with the high-band particle [2] at energy of order V below the bottom of the two-particle continuum. The spectral weight of the bound state scales as (V /J n ) 2 ≪ 1, where J n is the band-width of the n-the Bloch band. We can therefore safely neglect its contribution to S(q, ω).

This work was supported by ISF, MAE-MIUR (Joint Laboratory LENS-Weizmann), ECR Firenze, by ERC through the DISQUA grant, and by MIUR through PRIN2009. S.D.H acknowledges support by the Swiss Society of Friends of the Weizmann Institute of Science.

Supplementary material for the manuscript "Quasiparticle dynamics in a Bose insulator probed by inter-band Bragg spectroscopy"

N. Fabbri,

MEAN-FIELD THEORY FOR THE MOTT INSULATOR

The Mott insulator which is realized in the deep lattice limit U ≫ J 1 can be described using a strong-coupling mean-field approach [1]. In spirit it is similar to the Bogoliubov theory for weakly interacting particles

Here, β † p/h,k create a quasi-particle and a quasi-hole excitation, respectively. The parameters u k , v k behave as u k → 1 and v k → 0 for J 1 /U → 0 and the combination u k v k reflects the coherence in the Mott insulator. The operators p † i (h † i ) create a state with an additional (missing) particle at site i on top of the classical Gutzwiller state ∝ i (b † 1,i ) ρ|vac . Within a mean-field approach the effective Hamiltonian for the Mott insulator reads

where J c ≈ 8U ρ is the critical hopping for the formation of the Mott state with filling ρ and µ denotes the chemical potential. Note, that since the gap cuts off the infra-red divergencies characteristic of gapless 1d systems, the mean field theory is expected to work well inside the Mott phase. Details about the mean-field approach can be found in Ref. [1].

FINAL STATE INTERACTIONS

Let us now turn to the higher Bloch bands and the interaction effects therein. We treat the interactions by splitting the field operators into Wannier states for the lowest band and keep continuum states for the higher bands

. The part of the interaction which describes scattering between the final state and the remaining particles in the lower band is:

For later convenience let us also express the density operator in the lowest band in terms of the particle and hole fluctuation operators:

Hartree shift. The Hartree approximation of the final state interaction consists of replacing the density operator in the lower band by its average (i.e. neglecting the particle and hole fluctuations.

The effective potential obtained in this way is shown for 87 Rb at s x = 10 and ρ = 2 in Fig. 1(a). The addition of the Hartree term has two effects. First, it shifts the excited Bloch bands to slightly higher energies (by ∼ 0.2 E R for the third band) due to the repulsive interaction with the Mott insulating state. Second, the bandwidth of the high bands