Yuyao Yang

Robert Chahine

Robert Rubinstein

Wouter J T Bos

Mixing in modulated turbulence. Numerical results

Keywords:

Direct numerical simulations are carried out to investigate scalar mixing in an isotropic turbulent flow with a time-periodic forcing. For high amplitudes of the modulation, it is shown that the average mixing rate is negatively affected at low frequencies. In this limit the mixing time scale increases, whereas the typical velocity timescale decreases. We further determine the frequency response of scalar statistics to a periodic scalar-forcing.

I. INTRODUCTION

The mixing rate of a scalar quantity advected by a fluid is a key quantity in a wide range of applications. Increasing the mixing rate by changing the flow properties can have far reaching consequences in process-optimization. Whereas the mixing in laminar flows can often be studied analytically, and the mixing rate can be greatly enhanced by changing the boundary conditions [START_REF] Gouillart | Walls inhibit chaotic mixing[END_REF] or the time-dependence of the flow [START_REF] Aref | Stirring by chaotic advection[END_REF], the turbulent case is in general far more complicated. If any understanding of the modification of turbulent mixing through time-dependent forcing is to be obtained, we think it is compulsory to look at the most simplified case. We consider therefore the academic case of periodically forced isotropic turbulence, advecting a passive scalar.

The response of a turbulent velocity field on a time-periodic isotropic forcing has received a considerable interest since the beginning of the 2000s. The initial studies aimed at identifying a possible resonance in the energy transfer process [START_REF] Lohse | Periodically kicked turbulence[END_REF][START_REF] Der Heydt | Response maxima in modulated turbulence[END_REF][START_REF] Der Heydt | Response maxima in modulated turbulence. II. Numerical simulations[END_REF]. Direct numerical simulations [START_REF] Kuczaj | Response maxima in time-modulated turbulence: Direct numerical simulations[END_REF][START_REF] Kuczaj | Turbulence modification by periodically modulated scale-dependent forcing[END_REF] and experiments [START_REF] Cadot | Observation of resonances in modulated turbulence[END_REF][START_REF] Cekli | Resonant enhancement of turbulent energy dissipation[END_REF] were carried out to systematically investigate the response of the kinetic energy and dissipation rate on the forcing frequency. An analytical study, using two-point closure techniques, allowed to explain the different scaling regimes of the time-dependent quantities [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF]. The ability of engineering models to reproduce the different features of periodically forced turbulence was investigated in [START_REF] Rubinstein | On the unsteady behavior of turbulence models[END_REF].

The application of particular forcing schemes to influence turbulent mixing was considered in [START_REF] Kuczaj | Mixing in manipulated turbulence[END_REF]. In that study the forcing was introduced in different wavenumber bands in Fourier space to mimic the complex nature of turbulent flows generated by realistic objects.

The influence of the so-generated flow on turbulent mixing was assessed by monitoring the wrinkling of level-sets of an advected scalar. Those results inspired several experimental investigations with application to turbulent combustion (see for instance [START_REF] Verbeek | Optimization of combustion in gas turbines by applying resonant turbulence[END_REF][START_REF] Verbeek | Fractal turbulence enhancing low-swirl combustion[END_REF]). In [START_REF] Kuczaj | Mixing in manipulated turbulence[END_REF] the spatial character of the forcing was modified, but no temporal modulation was applied.

In the experiments the modulation was both spatial and temporal, and it is not straightforward to disentangle the different effects, so that it is not clear whether the observed effects were caused by the time-periodic nature of the experimental inlet conditions or the spatial complexity of the latter.

In the present study we consider the influence of a temporal modulation of the forcing on the mixing of a passive scalar. Two different cases are investigated. The first case corresponds to a stationary injection of scalar into a periodically forced velocity field. The second one is the case where, in a statistically steady isotropic turbulence, the scalar is introduced by a modulated injection. We will investigate both cases here by means of direct numerical simulation. In particular, we will show that in the case of a strong velocity modulation not only the periodic part, but also the time-average of the scalar variance and kinetic energy depend on the modulation frequency.

II. DEFINITIONS AND EQUATIONS

We consider the Navier-Stokes equations for incompressible flow u, mixing a passive scalar θ:

∂ ∂t u + u • ∇u = -∇P + ν∇ 2 u + f (1)
∇ • u = 0 (2) ∂ ∂t θ + u • ∇θ = D∇ 2 θ + g (3)
where P is the pressure (normalized by a uniform density), ν and D are kinematic viscosity and diffusivity, respectively. The flow and the scalar field are kept in a statistically stationary state through an energy and scalar variance input f , g. The forcing is applied to the smallest wavenumbers of the Fourier-transformed velocity and scalar field, denoted by ûi and θ, respectively,

fi (κ, t) = 1 N f ûi (κ, t) | û(κ, t)| 2 (p + p cos(ωt)) (4) ĝ(κ, t) = 1 N f θ(κ, t) θ(κ, t) 2 (p θ + pθ cos(ωt)), |κ| ≤ κ F , (5)
with N F the total number of forced modes. The quantities p and p θ denote the kinetic energy and scalar variance injection rates. Overlined quantities denote time-averages and tilded quantities denote periodic fluctuations around the temporal average. In applying such forcing schemes, a statistically steady state is observed where the phase averaged energy and scalar injection fluctuate periodically around the time-averaged values,

f i u i = p + p cos(ωt), (6)
gθ = p θ + pθ cos(ωt). (7)
These forcing schemes will result in a statistically isotropic velocity and scalar field. In our case the flow domain is a spatially periodic box. In such a setting the evolution equations for the kinetic energy k = 1 2 |u| 2 and the variance of the scalar

k θ = 1 2 θ 2 reduce to dk dt = p -ǫ (8)
dk θ dt = p θ -ǫ θ . (9)
In these equations ǫ and ǫ θ are the phase averaged dissipation of kinetic energy and scalar variance, respectively. All the different statistical quantities of interest will in the following be decomposed into a time-averaged and a periodic component. The time-averaged balance equations for the kinetic energy and scalar variance are given by 0 = P -ǫ (10)

0 = P θ -ǫ θ , (11)
and the periodic quantities evolve as

-ω k sin(ωt + φ k ) = p cos(ωt) -ǫ cos(ωt + φ ǫ ), (12)
-ω kθ sin(ωt + φ k θ ) = pθ cos(ωt) -ǫθ cos(ωt + φ ǫ θ ).

In these expressions, we have assumed that all quantities will periodically oscillate around a mean value with a period ω. In previous works, [START_REF] Der Heydt | Response maxima in modulated turbulence[END_REF][START_REF] Der Heydt | Response maxima in modulated turbulence. II. Numerical simulations[END_REF][START_REF] Kuczaj | Response maxima in time-modulated turbulence: Direct numerical simulations[END_REF][START_REF] Kuczaj | Turbulence modification by periodically modulated scale-dependent forcing[END_REF][START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF], the frequency dependence of k and ǫ was investigated.

We will start by reproducing these results in section IV. After that we will investigate the case where the relative periodic forcing amplitude α p = p/p is large enough to not only influence the periodic quantities, but also the time-averaged kinetic energy and scalar variance, k and k θ . These results show that by adding a large-amplitude modulation to the forcing, one can directly influence the transfer-rate of the kinetic energy χ, and the mixing rate of the passive scalar χ θ . These quantities can be defined as the inverse of the integral timescale and scalar timescale, respectively,

χ = ǫ k and χ θ = ǫ θ k θ . (14)
In section V, we will determine the frequency dependence of the periodically fluctuating scalar quantities, kθ , ǫθ for the case where the scalar injection rate is modulated.

A standard pseudospectral method is used to compute the velocity and scalar field in a space-periodic cubic domain of size 2π. A conventional 2/3 wavenumber truncation is used to eliminate the aliasing error and a third order Runge-Kutta, Total Variation Diminishing scheme is used as time discretisation. The same code was used in ref. [START_REF] Bos | On the scaling of temperature fluctuations induced by frictional heating[END_REF] to study the mixing of temperature fluctuations in isotropic turbulence. A challenge in the study of the frequency response of turbulent flows is the convergence of the statistics. At low forcing frequencies the simulations become very long if a sufficient number of periods is to be resolved. At high frequencies the response to a periodic forcing will be shown to be small, so that also in this case long simulations are needed, not to resolve sufficient periods, but to be able to distinguish the frequency response from the turbulent fluctuations. Obtaining converged statistics is therefore challenging in both the small and large frequency limits.

Previous investigations [START_REF] Kuczaj | Response maxima in time-modulated turbulence: Direct numerical simulations[END_REF][START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] focused in particular on the linear response of turbulence on a periodic modulation. In this limit linearized equations around a given equilibrium allow to analytically derive certain results. The verification of such results is not straightforward in the nonlinear regime, where the perturbation is large. Ideally, to investigate the linear response of a turbulent flow, the amplitude of the forcing should be chosen small compared to the amplitude of the steady part of the forcing (α p ≡ p/p ≪ 1). However, since the turbulent fluctuations are in this case much larger than the periodic response, very long simulations should be carried out to obtain an estimate of the frequency response. In particular at large frequencies, where the frequency response will be shown to drop rapidly as a function of frequency this would impose prohibitively long computations. A compromise is to consider a larger modulation amplitude. In this study, as in [START_REF] Kuczaj | Response maxima in time-modulated turbulence: Direct numerical simulations[END_REF], we use α p ≡ p/p = 0.2. This allows to obtain converged statistics for a large range of frequencies at a reasonable computational cost for low Reynolds number (R λ = 32). For higher Reynolds number this already leads to prohibitively long simulations. Therefore we have carried out another set of simulations with a relative modulation amplitude α p = 1. Even though this certainly violates the linear perturbation assumption, we will show that the frequency-response of the modulated quantities is not quantitatively altered. We will further show that this has an interesting direct influence on the time-averaged quantities.

Before extracting the frequency response of the simulations, the flow was simulated for approximately 10 eddy turn-over times to obtain a statistically steady state. It proved convenient to determine the amplitude of the periodic response by using a Fourier-transform of the signal. Before Fourier-transforming the time-series of a given quantity, a hanning window function is applied to the signal to eliminate the aliasing error at high frequencies due to the finite length of the signal. In the frequency spectra, if the simulations are carried out for a sufficiently long time-interval, the amplitude of the periodic response is easily identified by a sharp peak. Comparing the value of this peak to the neighbouring values in the spectrum gives a direct estimate of the signal-to-(turbulent)-noise ratio. In all simulations the value of the peak at the considered frequency was at least ten times the value of the neighbouring values in the spectra. For the phase averaged amplitudes, error-bars are added to the datapoints in the figures, computed from the signal to noise ratio. In most cases, this error-bar is smaller than the size of the symbols used in the figures and is then omitted. The time-averaged value is conveniently estimated from the ω = 0 component of the spectrum.

The different simulations we have carried out are documented in table I.

IV. RESPONSE OF TURBULENCE AND MIXING ON A PERIODIC KINETIC ENERGY INPUT

In this section we consider the case where we only modulate the kinetic energy,

p = p + p cos(ωt), (15)
p θ = p θ . (16)
It will be shown that the modulation p of the velocity field does also affect the mixing of the passive scalar.

A. Frequency response of the modulated kinetic energy and dissipation

The frequency response of k and ǫ is shown in Fig. 1 for R λ = 32. We compare in this figure the frequency responses for two different relative forcing amplitudes, α p = p/p = 0.2 and α p = 1. In order to compare the frequency response for the different forcing amplitudes, we plot in these figures the quantities

k * = α -1 p k k and ǫ * = α -1 p ǫ ǫ ( 17
)
as a function of frequency. Several observations can be made. Firstly, the 20% and 100%

relative forcing amplitudes give results that superpose at almost all frequencies. In the following we will therefore focus on α p = 1 results, which allow to obtain results at a lower computational cost, and therefore, at higher Reynolds number. Secondly, the powerlaw dependence proportional to ω -1 observed in [START_REF] Kuczaj | Response maxima in time-modulated turbulence: Direct numerical simulations[END_REF] and [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] is clearly reproduced both for k and ǫ. At small frequencies both k and ǫ tend to constant values, as predicted in [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF], but unlike the DNS results in [START_REF] Kuczaj | Response maxima in time-modulated turbulence: Direct numerical simulations[END_REF] at these frequencies, perhaps due to unconverged statistics in their simulations. The quality of the collapse of the data at different values of α p is not so good for k * in the low frequency limit. This is due to the value of k in this limit, which is affected by the strong modulation, as will be illustrated below, in section IV B.

Figure 2 illustrates the influence of the Reynolds number on the modulated kinetic energy and dissipation. It is observed that this influence is small for the modulated kinetic energy.

However, for the dissipation this influence is larger, as was explained in [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] by the fact that the ω -1 asymptote is inversely proportional to the Reynolds number, since it corresponds to the direct influence of the viscous damping on the forced scales. The intermediate zone between the low frequency plateau and the high frequency asymptote is significantly steeper than the prediction that it should be proportional to ω -3 for large Reynolds numbers.

Whether this is a low Reynolds number effect, or caused by the strong forcing, or due to something else can not be concluded from the present observations.

In Figure 3 (a) the phase-shifts of the kinetic energy -φ k and dissipation -φ ǫ with respect to the forcing are plotted. It is observed that -φ k evolves almost continuously from 0 at small frequencies, where everything is in phase, to the expected high frequency asymptote of 90 • . A small overshoot is observed before reaching the high frequency asymptote, as also reproduced using closure simulations [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] and in DNS [START_REF] Kuczaj | Response maxima in time-modulated turbulence: Direct numerical simulations[END_REF]. The phase-shift -φ ǫ also behaves as expected, with a peak value around the integral frequency. The difference between the two phases is proportional to ω for small frequencies (Fig. 3 (b)), reflecting the finite cascade-time between the large and small scales of the flow [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF]. Interestingly a small local maximum is observed at ωT ≈ 16. We have currently no explanation for this observation and since is does not occur in the linear predictions, we suspect this maximum to be related to a violation of the linear-response assumptions.

Overall, these results show that the observations from previous studies are reproduced even though the relative forcing amplitude α p is significantly higher in the present study.

B. Frequency response of the time-averages

Within the framework of linear response theory, we consider that we measure the response of a system to an infinitesimal perturbation. In the case of a periodic perturbation the response is assumed to be at the same frequency as the perturbation, around the unaltered system. When the perturbation is large, the system itself can be affected importantly, and linear-response theory is no longer valid. In our case, a perturbation of the injection with an amplitude equal to the average injection cannot possibly be considered infinitesimal. It is therefore interesting to see if such a perturbation modifies the time-averaged properties of the flow. We recall that in the present section pθ = 0: only the forcing of the velocity contains a time-periodic contribution. Figure 4 shows the effect of modulated energy input on average components. It is observed that for low frequencies a quantitative deviation is observed from the average observed at high frequencies. This deviation is significantly larger than the error in the convergence of the statistics. It is in particular observed that k(ω)/k(∞) decreases approximately to 93% for small frequencies. On the contrary, the scalar variance increases more than 20%, i.e., k θ (ω)/k θ (∞) ≈ 1.2. This effect is only observed for a forcing amplitude α p = 1. For α p = 0.2 the average seems unaffected, in agreement with the assumption of linear response. Also, two peaks in the frequency response of the average scalar variance are observed, one around the integral frequency, another one around ωT ≈ 20. The timeaveraged dissipation ǫ and scalar dissipation ǫ θ (not shown) are not measurably affected. This is understandable since the forcing protocol is designed to maintain a well-defined value of these quantities.

Since the mixing rate is defined as χ θ = ǭθ / kθ (see equation ( 14)), the results in fig 4 show that the mixing-rate is negatively affected by the periodic forcing. But this effect disappears at high frequencies. The opposite is shown for the velocity field: the kinetic energy transfer rate is enhanced at low frequencies. These results show that a periodic forcing can modify the mixing rate and the transfer rate when the modulation amplitude is large and the frequency is low. The requirement that the relative amplitude should be large to affect the averages is an indication that the effect is a nonlinear correction to the linear response prediction, where the time-average is not affected.

C. Frequency response of the modulated scalar variance and dissipation

Since the scalar forcing does not contain a time-periodic component, the equation for the modulated scalar variance is,

-ω kθ sin(ωt + φ k θ ) = -ǫ θ cos(ωt + φ ǫ θ ). ( 18
)
From this equation it is not obvious that the scalar should contain a periodic component.

Indeed from (18), taking the modulus, we only find that

ω| kθ | = |ǫ θ |. (19)
No further information can be obtained from these equations about the frequency dependence of kθ . Obviously, kθ = 0 is a possible solution of these equations. But since the velocity field advecting the scalar contains a periodic contribution, it is not excluded that a periodic contribution is observed in the scalar dynamics. This is indeed the case, as observed in Figure 5. It is observed that the scalar variance contains a periodic component which is constant at low frequencies, and rapidly drops off at high frequencies, following a powerlaw proportional to ω -3 . The periodic part of the dissipation is then determined by relation (19), as seems to be confirmed in figure 5. The origin of the ω -3 powerlaw is not easily obtained from the single point equations.

In Figure 6 the phaseshifts φ kθ and φ ǫθ are shown. As expected from equation (18) the relative phaseshift between kθ and ǫθ is always equal to π/2. We also consider the phaseshift between the signals of the modulated scalar variance and the modulated kinetic energy. If we evaluate this relative phaseshift at the two frequencies were the time-averaged scalar variance showed a clear peak (Fig. 4) it is found that this corresponds to -(φ k θ -φ k ) ≈ 0 or ≈ π. It seems that a kind of resonance is observed in the scalar response, when the scalar variance and the kinetic energy are completely in phase, or completely out of phase.

V. MODULATION OF THE SCALAR INJECTION

In the foregoing, we have considered the case where the velocity was stirred at large scales by a periodic force term, but where the injection of the scalar was not modulated. In this section, we consider a modulated scalar injection in a non-time-periodic velocity field. In Figure 7 we show the results for the modulated scalar variance kθ / kθ and dissipation ǫθ /ǭ θ .

A clear resemblance with the frequency behaviour of the modulated kinetic energy and dissipation, in Figure 2, is observed. In particular the small and large frequency asymptotes are identical. Indeed, the reasonings leading to the prediction of the frequency behaviour of the kinetic energy and dissipation [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] can be extended to the case of the passive scalar, in particular in the limit of the linear response approximation. A further resemblance is observed in the phase-shift of the scalar variance and dissipation, shown in figure 8. It seems that the results of the periodically forced velocity can be transposed to the periodically forced scalar, without major modifications.

VI. CONCLUSION

Perhaps the most interesting outcome of the present investigation is that the periodical modulation of a turbulent velocity field can directly affect its mixing rate. This effect was only observed when the relative forcing frequency was large (α p = 1), suggesting that the effect is nonlinear. The average mixing rate was shown to decrease by approximately 20%, for two distinct Reynolds numbers, whereas the energy transfer rate increases by approximately 10%. Indeed, in the linear limit this effect is absent.

We further showed that in this nonlinear regime the frequency response of the amplitudes of the kinetic energy and dissipation were still in agreement with the outcome from linear perturbation analysis. Interestingly, the modulated velocity also directly induces a modulation of the scalar field. It is therefore, according to our investigation possible to modify the time-averaged and phase-averaged scalar by acting on the velocity field.

By adding a modulation on the scalar injection rate, it was shown that the frequency response of the periodic contributions of the scalar variance and the scalar dissipation behave

  Simulations are carried out at two different Reynolds numbers. First, low resolution simulations at a spatial resolution of 64 3 are performed with kinematic viscosity ν = 0.009, corresponding to a Taylor Reynolds number Re λ = 32, with eddy turn-over time T = 1.949. The resolution allows to resolve the smallest scales upto k max η = 0.93. Another set of simulations is carried out at a resolution of N 3 = 256 3 gridpoints, kinematic viscosity ν = 0.0009; Taylor Reynolds number Re λ = 105; eddy turn-over time T = 2.317 and k max η = 0.97. In all simulations the Schmidt number Sc ≡ ν/D = 1. The forced modes correspond to 0.5 < |κ| ≤ 2.5.

FIG. 1 :

 1 FIG. 1: The frequency dependence of (a) k * = k k p p -1 and (b) ǫ * = ǫ ǫ p p -1 as a function

FIG. 2 :

 2 FIG. 2: Frequency response of (a) the modulated kinetic energy k/k and (b) modulated dissipation ǫ/ǫ for R λ = 32( * ) and R λ = 105(×). The relative forcing amplitude is α p = 1.

FIG. 3 :

 3 FIG. 3: (a) Phase shifts -φ k ( * ) and -φ ǫ (•) as a function of ωT (in degrees) for R λ = 105 and α p = 1. (b) Relative Phase shift -(φ ǫ -φ k) as a function of ωT for R λ = 105 and α p = 1.

FIG. 4 :

 4 FIG. 4: Normalized time-averaged kinetic energy k/ k(ω = ∞) (a) and scalar variance kθ / kθ (ω = ∞) (b) as a function of ωT . Results for R λ = 32, α p = 0.2 (•); R λ = 32, α p = 1, ( * ); and R λ = 105, α p = 1 (×).

FIG. 5 :

 5 FIG. 5: Frequency response of (a) the modulated scalar variance kθ /k θ and (b) modulated scalar dissipation ǫθ /ǫ θ for R λ = 32 ( * ) and R λ = 105 (×). The relative forcing amplitude is α p = 1.

FIG. 6 :

 6 FIG. 6: (a) Phaseshifts -φ k θ ( * ) and -φ ǫ θ (•) as a function of ωT . (b) Phaseshift -(φ k 2 θ -φ k ). All results are given for R λ = 105 and α p = 1.

FIG. 7 :FIG. 8 :

 78 FIG. 7: Amplitudes of the modulated scalar variance kθ / kθ (a) and dissipation ǫθ /ǭ θ (b) as a function of ωT for the case of a modulated scalar injection. Results for R λ = 32( * ) and R λ = 105(×), both at α p = 1.


Keywords: Modèle linéaire généralisé, Effet aléatoire, Modèle de mélange, Algorithme EM, Algorithme de Metropolis-Hastings

Modèles de mélange fini pour des données exponentielles répétées Résumé : Nous nous intéressons à un modèle de mélange pour des données répétées de loi exponentielle. Les composants du mélange traduisent différents états possibles des individus. Pour chacun de ces composants, on modélise la dépendance et l'extra-variabilité dues à la répétition des données par l'introduction d'effets aléatoires. Dans ce modèle de mélange exponentiel mixte, la distribution marginale n'étant pas accessible, l'utilisation de l'algorithme EM n'est pas directement envisageable. Nous proposons alors une première méthode d'estimation des paramètres basée sur une linéarisation spécifique à la loi exponentielle. Nous proposons ensuite une méthode plus générale puisque s'appuyant sur une étape de Metropolis-Hastings pour construire un algorithme de type MCEM. Cet algorithme est applicable pour un mélange de modèles linéaires généralisés mixtes quelconques.

Introduction

In the past decades, finite mixture models have been extensively developed in the literature. Surveys on these issues can be found in Titterington, Smith and Makov [START_REF] Titterington | Statistical analysis of finite mixture distributions[END_REF], McLachlan and Basford [START_REF] Mclachlan | Mixture models: inference and application to clustering[END_REF] and McLachlan and Peel [START_REF] Mclachlan | Finite mixture models[END_REF]. In such finite mixture models, it is assumed that a sample of observations arises from a specified number of underlying groups or classes with unknown proportions and according to a specific form of distribution in each of them. A large number of distributions from the exponential family have been considered such as normal, Poisson, exponential (Hasselblad [6]). Wedel and DeSarbo [START_REF] Wedel | A mixture likelihood approach for generalized linear models[END_REF] have proposed a mixture of generalized linear models which contains the previously proposed mixtures as special cases, as well as a host of other parametric specifications theretofore not dealt with in the literature. The use of these mixture models can be, in particular, a way to consider unexpected variance in GLM's or also a way to take into account underlying unobserved latent variable forming groups or classes. More recently, Celeux, Martin and Lavergne [START_REF] Celeux | Mixture of linear mixed models for clustering gene expression profiles from repeated microarray experiments[END_REF] have proposed a mixture of linear mixed models (LMM) in a microarray data analysis context. The introduction of random effects allowed them to take into account the variability of gene expression profiles from repeated microarray experiments. In our work, we consider the analysis of finite mixture for exponential repeated data. The mixture components correspond to different possible states of the statistical units. Dependency and variability of exponential repeated data are taken into account through exponential mixed models defined for each mixture component. In the field of the Health Sciences, applications may concern the modelling of lengths of repeated hospital stays for patients belonging to unknown clusters. Another example is the analysis of elimination times after repeated absorptions of a drug by patients not being controlled a priori. Concerning parameter estimation in the proposed mixture of exponential mixed models, the use of the EM algorithm which allows to take into account the incomplete structure of the data is considered (McLachlan and Krishnan [START_REF] Mclachlan | The EM algorithm and extensions[END_REF]). But the algorithm presented by Celeux et al. [START_REF] Celeux | Mixture of linear mixed models for clustering gene expression profiles from repeated microarray experiments[END_REF] for a mixture of linear mixed models cannot here be transposed because the marginal distribution of each mixture component cannot be analytically derived. Thus, we propose two parameter estimation methods. The first one uses a linearisation specific to the exponential distribution hypothesis associated with each mixture component. The second approach is adapted from the algorithm presented by McCulloch [START_REF] Mcculloch | Maximum likelihood algorithms for generalized linear mixed models[END_REF] for generalized linear mixed models (GLMM) [START_REF] Mcculloch | Generalized, linear and mixed models[END_REF] and uses a Metropolis-Hastings step (Hastings [7]) to allow construction of an MCEM algorithm. This algorithm can be adapted to a mixture of any generalized linear mixed models. The paper is organized as follows. After a description of the model hypotheses in section 2, we outline the EM algorithm presented by Celeux et al. for a mixture of linear mixed models in section 3. In section 4, we describe the two proposed parameter estimation methods. Finally in section 5, we study the behaviour of these approaches on simulations.

INRIA 2 Mixture of exponential mixed models : model definition

Consider y = (y 1 , . . . , y I ) a vector of observations where y i is associated with the ith statistical unit. Each y i contains the n i repetitions y ij . Consider also different components C k , k = 1, . . . , K, corresponding to different possible states of the statistical units. We assume that all repeated measures of a statistical unit belong to the same component and we define the indicator vectors z i = (z i1 , . . . , z iK ), i = 1, . . . , I, with z ik = 1 if unit i ∈ C k and 0 otherwise.

To take into account dependency and variability of repeated data, we consider for each component an exponential mixed model with a random effect associated with each unit. This leads to a mixture of exponential mixed models and the density of Y i may be written as follows:

f (y i |θ, p) = K k=1 p k f k (y i |θ k )
where the p k 's are mixing weights with 0 < p k < 1 for k = 1, . . . , K and K k=1 p k = 1, and f k (.|θ k ) denotes the density function of the marginal distribution associated with the exponential mixed model with unknown parameters θ k = (β k , σ 2 k ). Note that this marginal distribution cannot be analytically derived.

More precisely, given the mixture component C k from which unit i arises and given the unobserved random effect ξ i , Y i is assumed to be exponentially distributed:

(Y i |ξ i , Z ik = 1) ∼ Exp(µ k ξ,i ) with µ k ξ,i = exp(X i β k + U i ξ i ) (ξ i |Z ik = 1) ∼ N (0, σ 2 k )
where

• ∀i, i ∈ {1, . . . , I} 2 i = i , ξ i and ξ i are assumed to be independent,

• β k is the q × 1 fixed effect parameter vector associated with component C k ,

• σ 2 k is the variance of the random effect associated with component C k ,

• X i =    x i1 . . .
x ini    and U i = (u i1 , . . . , u ini ) are the n i × q and n i × 1 known design matrices.

Thus, we focus here on a mixture model-based approach to the clustering of exponential repeated data (McLachlan and Peel [START_REF] Mclachlan | Finite mixture models[END_REF]).

EM algorithm for a mixture of linear mixed models

In this section, we outline the maximum likelihood estimation approach for a mixture of linear mixed models using the EM algorithm presented by Celeux et al. [START_REF] Celeux | Mixture of linear mixed models for clustering gene expression profiles from repeated microarray experiments[END_REF]. The EM methodology takes into account the incomplete structure of the data (Dempster, Laird and Rubin [START_REF] Dempster | Maximum likelihood for incomplete data via the EM algorithm[END_REF]). Missing data are here of two types : the indicator vectors z i , i = 1, . . . , I of unit memberships to the mixture components and the random effects ξ i , i = 1, . . . , I.

Given the mixture component C k from which unit i arises, Y i is here assumed to be normally distributed:

(

Y i |Z ik = 1) = X i β k + U i ξ i + ε i
where

• (ξ i |Z ik = 1) ∼ N (0, σ 2 k ),
• ε i is the n i × 1 error vector assumed to be normally distributed: ε i ∼ N (0, τ 2 I ni ) with I ni the identity matrix of order n i .

• ∀i ∈ {1, . . . , I}, ε i and ξ i are assumed to be independent and ∀i, i ∈ {1, . . . , I} 2 i = i , ξ i and ξ i , respectively ε i and ε i , are assumed to be independent,

• β k is the q × 1 fixed effect parameter vector associated with component C k ,

• σ 2 k is the random effect variance associated with component C k ,

• X i =    x i1 . . . x ini  
 and U i = (u i1 , . . . , u ini ) are the n i × q and n i × 1 design matrices.

Thus the distribution of Y i is a mixture of linear mixed models defined by

f (y i |θ, p) = K k=1 p k f k (y i |θ k )
where p = (p 1 , . . . , p K ) are the mixing weights, θ = (θ 1 , . . . , θ K ) with θ k = (β k , σ 2 k , τ 2 ) the linear mixed model parameters associated with component C k , and f k (y i |θ k ) denotes the density function of the distribution of Y i i.e. a Gaussian distribution with mean

X i β k and variance matrix Γ k,i = τ 2 I ni + σ 2 k U i U i .
In their paper, Celeux et al. [START_REF] Celeux | Mixture of linear mixed models for clustering gene expression profiles from repeated microarray experiments[END_REF] consider a mixture model where all parameters are dependent on component C k . We consider here a mixture model where the parameters β k and σ 2 k depend on component C k whereas the residual variance τ 2 is the same for all mixture components.

The log-likelihood associated with the complete data (y, z, ξ) is given by

L(θ, p|y, z, ξ) = I i=1 K k=1 z ik ln p k + ln f (y i , ξ i |z ik = 1, θ k ) INRIA where ln f (y i , ξ i |z ik = 1, θ k ) can be written as ln f (y i |ξ i , z ik = 1, θ k ) + ln f (ξ i |z ik = 1, θ k ) with • ln f (y i |ξ i , z ik = 1, θ k ) = -1 2 n i ln 2π + n i ln τ 2 + ε i ε i τ 2 , = -1 2 n i ln 2π + n i ln τ 2 + (yi-Xiβ k -Uiξi) (yi-Xiβ k -Uiξi) τ 2 . • ln f (ξ i |z ik = 1, θ k ) = -1 2 ln 2π + ln σ 2 k + ξ 2 i σ 2 k At iteration [t + 1]
, the E step consists of computing the expectation of the complete data log-likelihood given the observed data and a current value of the parameters (θ [t] , p [t] ):

Q(θ, p|θ [t] , p [t] ) = E L(θ, p|y, z, ξ)|y, θ [t] , p [t] = I i=1 K k=1 t [t] k (y i ) ln p k - 1 2 I i=1 K k=1 t [t] k (y i ) (n i + 1) ln 2π + n i ln τ 2 + ln σ 2 k + E [t] Ck (ε i ε i ) τ 2 + E [t] Ck (ξ 2 i ) σ 2 k where E [t] Ck (.) = E(.|y i , z ik = 1, θ [t] k ) and t [t] k (y i ) = P (Z ik = 1|y i , θ [t] , p [t] ) = p [t] k f (y i |z ik = 1, θ [t] k ) f (y i |θ [t] , p [t] ) = p [t] k f k (y i |θ [t] k ) K l=1 p [t] l f l (y i |θ [t]
l ) denotes the posterior probability that unit i arises from component C k .

The M step consists of maximizing Q(θ, p|θ [t] , p [t] ). It leads to the following explicit expressions for k = 1, . . . , K:

p [t+1] k = I i=1 t [t] k (y i ) I β [t+1] k = I i=1 t [t] k (y i )X i X i -1 I i=1 t [t] k (y i ) τ 2[t] X i Γ [t]-1 k,i (y i -X i β [t] k ) + X i X i β [t] k σ 2[t+1] k = 1 I i=1 t [t] k (y i ) I i=1 t [t] k (y i ) σ 4[t] k (y i -X i β [t] k ) Γ [t]-1 k,i U i U i Γ [t]-1 k,i (y i -X i β [t] k ) +σ 2[t] k -σ 4[t] k tr(Γ [t]-1 k,i U i U i ) τ 2[t+1] = 1 n I i=1 K k=1 t [t] k (y i ) τ 4[t] (y i -X i β [t] k ) Γ [t]-1 k,i Γ [t]-1 k,i (y i -X i β [t] k ) + n i τ 2[t] -τ 4[t] tr(Γ [t]-1 k,i ) .
Details can be found in Celeux et al. [START_REF] Celeux | Mixture of linear mixed models for clustering gene expression profiles from repeated microarray experiments[END_REF].

Two parameter estimation methods

We consider here the parameter estimation for the mixture of exponential mixed models presented in section 2. In this context, the use of the EM algorithm is not directly possible.

The complete data log-likelihood associated to this model is given by

L(θ, p|y, ξ, z) = I i=1 K k=1 z ik ln p k + I i=1 K k=1 z ik ln f (y i |ξ i , z ik = 1, θ k ) + I i=1 K k=1 z ik ln f (ξ i |z ik = 1, θ k ) with • ln f (y i |ξ i , z ik = 1, θ k ) = ni j=1 ln f (y ij |ξ i , z ik = 1), = - ni j=1 x ij β k + u ij ξ i + y ij exp(x ij β k + u ij ξ i )
because the y ij 's are independent conditionally on ξ i .

• ln f (ξ i |z ik = 1, θ k ) = - 1 2 ln 2π + ln σ 2 k + ξ 2 i σ 2 k
In this case, at iteration [t+1], the EM algorithm leads to formulae depending on conditional expectations

E [t] Ck (ξ 2 i ), E [t]
Ck exp(-u ij ξ i ) and posterior probabilities t

[t]
k (y i ), i = 1, . . . , I, k = 1, . . . , K. Because of the non-availability of the marginal distribution for each mixture component, probabilities t

[t]
k (y i ) cannot be derived in closed form. Furthermore, neither the conditional expectation

E [t] Ck (ξ 2 i ) nor E [t]
Ck exp(-u ij ξ i ) can be computed too since these calculations involve the unknown conditional distribution of ξ i given y i . We propose two parameter estimation methods which allow to get round these problems.

A method based on linearisation

This first approach is a conceptually simple method which involves two steps: a linearisation specific to the exponential mixed model (Gaudoin, Lavergne and Soler [START_REF] Gaudoin | A generalized geometric de-eutrophication software reliability model[END_REF]) associated with each mixture component and the use of the EM algorithm for parameter estimation in a mixture of linear mixed models.

Knowing the component C k , the distribution associated with statistical unit i is given by

(Y i |ξ i , Z ik = 1) ∼ Exp(µ k ξ,i ),
or equivalently :

Y i µ k ξ,i ∼ Exp(1), thus ln(Y i ) -ln(µ k ξ,i ) ∼ Gumbel
, where the Gumbel density function is defined by ∀t ∈ IR f (t) = exp(texp(t)) with mean γ = -0.57722 and variance π 2 6 . This enables us to write:

ln(Y i ) -ln(µ k ξ,i ) = γ + ε i where E(ε i ) = 0 ni and var(ε i ) = π 2 6 I ni .
Defining the variable D i = log(Y i )γ, we end up with the linearized model:

D i = X i β k + U i ξ i + ε i
with 0-mean error vector ε i and known variance matrix π 2 6 I ni , which is viewed and considered as a linear mixed model M k for the data d i = log(y i )γ given the component C k .

Finally, we use the EM algorithm to estimate the parameters of the mixture of linear mixed models defined by

h(d i |θ, p) = K k=1 p k h k (d i |θ k )
where d i = ln(y i )γ and h k (d i |θ k ) is the density function of the Gaussian distribution with mean vector

X i β k and variance matrix Γ k,i = π 2 6 I ni + σ 2 k U i U i .
In this approach, note that vector d i is derived from the data y i whatever the component C k from which unit i arises and without any use of the current value of the parameters.

The parameter estimation for this mixture of linear mixed models using the EM algorithm described in section 3 leads to the following expressions for k = 1, . . . , K:

p [t+1] k = I i=1 t [t] k (d i ) I β [t+1] k = I i=1 t [t] k (d i )X i X i -1 I i=1 t [t] k (d i ) π 2 6 X i Γ [t]-1 k,i (d i -X i β [t] k ) + X i X i β [t] k σ 2[t+1] k = 1 I i=1 t [t] k (d i ) I i=1 t [t] k (d i ) σ 4[t] k (d i -X i β [t] k ) Γ [t]-1 k,i U i U i Γ [t]-1 k,i (d i -X i β [t] k ) +σ 2[t] k -σ 4[t]
k tr(Γ

[t]-1 k,i U i U i )

An MCEM algorithm

The proposed algorithm is adapted from the MCEM algorithm presented by McCulloch [START_REF] Mcculloch | Maximum likelihood algorithms for generalized linear mixed models[END_REF] for generalized linear mixed models. Since expectations

E [t]
Ck (ξ 2 i ) and

E [t]
Ck exp(-u ij ξ i ) and posterior probabilities t k (y i ) cannot be derived in closed form, our goal is to form Monte Carlo approximations of these quantities. To this aim, we incorporate a Metropolis-Hastings step into the EM algorithm which does not require specification of the marginal distribution of Y i . This leads us to draw values from the unknown conditional distribution of ξ i given Y i , Z ik = 1 and the current value θ k . One can then calculate Monte Carlo approximations of the two required expectations. In the same way, we draw values from the known distribution of ξ i given Z ik = 1 and the current value θ k (y i ). Before presenting the proposed algorithm in section 4.2.2, we recall the Metropolis-Hastings step applied to our specific case in section 4.2.1.

The Metropolis-Hastings step

The Metropolis-Hastings algorithm [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] is certainly one of the most famous MCMC methods (Robert and Casella [START_REF] Robert | Monte Carlo statistical methods[END_REF]). The aim of the MCMC methods is to generate samples from a target distribution π unavailable in closed form. To this end, a candidate distribution h (called the instrumental or proposal distribution) must be specified from which potential new values are drawn. Among samples generated from h, Metropolis-Hastings selects representative samples of the target distribution π using an acception/rejection method.

To define the proposed Metropolis-Hastings step, we need to specify the candidate distribution h. We propose to take h equal to the marginal distribution in class C k of ξ i given the current value θ

ρ(ξ [m] i , ξ * i ) = min 1, f (ξ * i |y i , z ik = 1, θ [t] k )h(ξ [m] i ) f (ξ [m] i |y i , z ik = 1, θ [t] k )h(ξ * i )
where the second term simplifies to:

f (ξ * i |y i , z ik = 1, θ [t] k )h(ξ [m] i ) f (ξ [m] i |y i , z ik = 1, θ [t] k )h(ξ * i ) = f (ξ * i |y i , z ik = 1, θ [t] k )f (ξ [m] i |z ik = 1, θ [t] k ) f (ξ [m] i |y i , z ik = 1, θ [t] k )f (ξ * i |z ik = 1, θ [t] k ) = f (y i |ξ * i , z ik = 1, θ [t] k ) f (y i |ξ [m] i , z ik = 1, θ [t] k )
.

By choosing h equal to the random effects distribution, probability ρ is simplified since the obtained formula only involves the specification of the conditional distribution of Y i given ξ i and the component C k from which unit i arises.

The proposed MCEM algorithm

Incorporating this Metropolis-Hastings step into the EM algorithm gives the following Monte Carlo EM (MCEM) algorithm at iteration [t + 1]:

1. For i = 1, . . . , I and k = 1, . . . , K, draw:

-M values ξ [1] i , . . . , ξ [M ] i from the distribution of ξ i |Y i , Z ik = 1 given the current value θ [t]
k using the Metropolis-Hastings algorithm described above and use them to form Monte Carlo approximations of the two required expectations in the function Q(θ, p|θ [t] , p [t] ):

E [t] Ck (ξ 2 i ) 1 M M m=1 ξ [m]2 i E [t] Ck exp(-u ij ξ i ) 1 M M m=1 exp(-u ij ξ [m] i ) -N values ξ [1] i , . . . , ξ [N ] i
from the known distribution of ξ i given Z ik = 1 and the current value

θ [t]
k in order to approximate the marginal distribution:

f k (y i |θ [t] k ) = f (y i |z ik = 1, θ [t] k ) = ni j=1 f (y ij |ξ i , z ik = 1, θ [t] k ) f (ξ i |z ik = 1, θ [t] k ) dξ i ≈ 1 N N n=1 ni j=1 f (y ij |ξ [n] i , z ik = 1, θ [t] k ) RR n°6119
and to obtain an approximation of the posterior probability t

[t]

k (y i ).

2. Then maximise the function Q(θ, p|θ [t] , p [t] ) to obtain new parameter values θ [t+1] and p [t+1] .

5 Simulation results

Preliminary results

In order to study the behaviour of the MCEM algorithm developed in section 4.2, we first consider its use in the gaussian case. Indeed, in this case, the performances of the MCEM algorithm can easily be compared to those of the EM algorithm. We simulate 100 data sets. We give in Table 1 the mean and standard deviation of the 100 estimated values obtained with both EM and MCEM for each parameter. We consider here a two-component mixture model. We set the number of statistical units I equal to 100 and we consider the same number of repetitions for each unit: ∀i = 1, . . . , I n i = J = 6. The mixing parameters are p 1 = 0.6 and p 2 = 0.4. The random effect variances are σ 2 1 = 0.2 and σ 2 2 = 0.8 and the residual variance is τ 2 = 2. We consider a unique fixed effect parameter by component: 1 clearly shows that the MCEM algorithm performs close to the EM algorithm. However, it is important to note that the MCEM algorithm is numerically intensive. For instance, the EM algorithm implemented using R requires here only a few minutes whereas the MCEM algorithm implemented in C takes a few hours.

β 1 = -2 and β 2 = 2.

INRIA

Comparison of the two proposed methods

Simulations are performed to assess the ability of the proposed methods to estimate mixture parameters and to highlight the interest of taking into account repetitions. We consider a two-component mixture model. We set the number of statistical units I equal to 100. The mixing parameters are p 1 = 0.6 and p 2 = 0.4 and the random effect variances are σ 2 1 = 0.2 and σ 2 2 = 0.8. We also consider one fixed effect parameter by component and we generate samples from:

• model (A) defined by β 1 = -3 and β 2 = 3,

• model (B) defined by β 1 = -1 and β 2 = 1.

It allows us to assess the ability of the methods to correctly estimate parameters when the mixture components are more (model (A)) or less separated (model (B)). In order to study the impact of the number of repetitions on the quality of the estimations, we also consider different number of repetitions J: we take J = 4 and J = 8. Table 2 provides the mean and standard deviation of the estimations obtained from 100 samples generated from model (A). Table 3 gives the results obtained for model (B). Table 2 shows, in each situation (J = 4 and J = 8), that both methods provide accurate parameter estimations. As expected, the precision of the estimation depends on the random effect variances : the greater the variance, the greater the estimation's standard deviation. Note that the results obtained with the MCEM algorithm are slightly better than those obtained with the method based on linearisation. Table 2 also shows that the number of repetitions has an influence on the quality of the estimations. The results shown in Table 3 are obtained when the mixture components are less separated. They are not as adequate as the first case but they are still reasonable. Remarks similar to those made for the first case can be made. We just note that the impact of the number of repetitions is more important in this case.

Table 4 provides the correct classification rate for each model using the maximum a posteriori (MAP) decision rule from the estimate parameter values p, θ. The MAP decision rule consists of assigning all the measures of unit i to the mixture component C k such as k = argmax 1≤l≤K t l (y i ) with t l (y i ) = P (Z il = 1|y i , p, θ). The obtained results are globally satisfactory. Table 4 shows that the correct classification rate decreases when the random effect variance increases. It also clearly shows that the correct classification rate increases with the number of repetitions. Finally, we also note that the rates obtained with the MCEM algorithm are slightly better than those obtained with the method based on linearisation.

Discussion and conclusions

In this paper, we define a new class of models for repeated data: mixtures of generalized linear mixed models. These models allow us to introduce a notion of heterogeneity in the GLMM. They take dependency and variability of repeated data into account through random effects defined for each mixture component. We proposed two parameter estimation methods for these models: the MCEM algorithm which can be used for mixtures of any generalized linear mixed models and the method based on linearisation specific to a mixture of exponential mixed models. These two methods are adaptations of the EM algorithm getting round problems related to the direct use of it.

The simulations performed in the exponential case show that the two proposed parameter estimation methods globally perform well. They also show that the MCEM algorithm gives slightly better results than the method based on linearisation. This behaviour difference is even greater in difficult situations with less separated mixture components or low number of repetitions. However, it is important to note that, like all MCMC approaches, the MCEM algorithm is numerically intensive since a large number of simulations is required at each iteration. In practice, a compiled programming language had to be used to reduce computation times. Moreover, even if this algorithm seems to perform well in practice, it is necessary to note that we still have not established theoretical results for convergence. On the contrary, the implementation of the method based on linearisation is fast and can easily be done with R for instance. Nevertheless, the use of this method is restricted to the exponential case.

Coming back to the numerically intensive problem of the MCEM algorithm, it would be interesting to propose an intermediate version using simulation via a stochastic approximation in order to avoid calculations. A future work could adapt the method developed by Kuhn and Lavielle [START_REF] Kuhn | Coupling a stochastic approximation version of EM with an MCMC procedure[END_REF]. In their paper, Kuhn and Lavielle proposed an algorithm combining the stochastic approximation version of EM (SAEM) [START_REF] Celeux | A stochastic approximation type EM algorithm for the mixture problem[END_REF] [3] with a Markov chain Monte Carlo procedure.

Finally, in this paper, it is assumed that the number of components is known. However, in practical situations, this is mostly not the case. It thus becomes a part of the estimation process. To determine the appropriate number of components, it would be interesting to consider the model selection problem for mixtures of generalized linear mixed models. Model selection criteria proposed by Martinez [START_REF] Martinez | Modèles linéaires généralisés à effets aléatoires : contributions au choix de modèle et au modèle de mélange[END_REF] and Lavergne et al. [START_REF] Lavergne | Empirical model selection in generalized linear mixed effects models[END_REF] for generalized linear mixed models could be adapted to mixtures of generalized linear mixed models.

INRIA

INRIA

k

  in order to approximate marginal distribution f k (y i |θ [t] k ) by Monte Carlo methods and to calculate posterior probability t [t]

  draw from the conditional distribution of ξ i |Y i , Z ik = 1 given the current value θ [t] k . The probability of accepting the INRIA new value ξ * i generated using the candidate distribution h is given by

Table 1 :

 1 Parameter estimation results obtained with EM and MCEM in the gaussian case on 100 simulated data sets

		EM		MCEM
	Simulated values	mean	s.d.	mean	s.d.
	p 1 = 0.6	0.6006 0.0171	0.6006 0.0170
	C 1 β 1 = -2 σ 2 1 = 0.2	-1.9872 0.1141 -1.9873 0.1140 0.1994 0.1201 0.1991 0.1166
	p 2 = 0.4	0.3994 0.0171	0.3994 0.0170
	C 2 β 2 = 2 σ 2 2 = 0.8	2.0289 0.1737 0.7657 0.2897	2.0292 0.1733 0.7605 0.2855
	τ 2 = 2	2.0210 0.1340	2.0216 0.1339
	Table

Table 2 :

 2 Parameter estimation results from 100 simulated models defined by β 1 = -3 and β 2 = 3

		Model (A)	Model (A')
		J = 4	J = 8
	Simulated values Linear.	MCEM	Linear.	MCEM
	p 1 = 0.6	0.6017	0.5999	0.6002	0.6001
		(0.0039) (0.0023) (0.0012) (0.0010)
	C 1 β 1 = -3	-2.9990	-2.9945	-3.0154	-3.0057
		(0.0878) (0.0815) (0.0921) (0.0817)
	σ 2 1 = 0.2	0.2166	0.1960	0.1986	0.1977
		(0.1227) (0.0834) (0.0749) (0.0645)
	p 2 = 0.4	0.3983	0.4001	0.3998	0.3999
		(0.0039) (0.0023) (0.0012) (0.0010)
	C 2 β 2 = 3	3.0181	3.0042	3.0105	3.0121
		(0.1765) (0.1641) (0.1555) (0.1588)
	σ 2 2 = 0.8	0.7419	0.7953	0.7798	0.7893
		(0.2490) (0.2588) (0.2455) (0.2280)

Table 3 :

 3 Parameter estimation results from 100 simulated models defined by β 1 = -1 and β 2 = 1

		Model (B)	Model (B')
		J = 4	J = 8
	Simulated values Linear.	MCEM	Linear.	MCEM
	p 1 = 0.6	0.6634	0.6536	0.6481	0.6330
		(0.1340) (0.0851) (0.0856) (0.0797)
	C 1 β 1 = -1	-0.8908	-0.9334	-0.9485	-0.9601
		(0.1845) (0.1494) (0.1599) (0.1371)
	σ 2 1 = 0.2	0.2760	0.2376	0.2523	0.2330
		(0.2251) (0.1466) (0.1324) (0.1119)
	p 2 = 0.4	0.3366	0.3464	0.3519	0.3670
		(0.1340) (0.0851) (0.0856) (0.0797)
	C 2 β 2 = 1	1.2909	1.2202	1.2036	1.1504
		(0.5317) (0.3047) (0.3426) (0.2895)
	σ 2 2 = 0.8	0.5437	0.6173	0.6195	0.6795
		(0.4568) (0.3591) (0.3257) (0.3274)

INRIA

Table 4 :

 4 Correct classification rates (%) from 100 simulations

		Model (A) Model (A') Model (B) Model (B')
		β = (-3, 3) β = (-3, 3) β = (-1, 1) β = (-1, 1)
		J = 4	J = 8	J = 4	J = 8
	Linear. C 1	99.98	100.00	93.37	96.38
	C 2	99.55	99.95	67.87	76.95
	MCEM C 1	99.96	100.00	96.15	96.52
	C 2	99.92	99.97	73.17	79.95

RR n°6119


Myron Thomas Curtright

Thomas Curtright Five

Five Talks with Myron

Thomas Curtright, University of Miami Please go along with me on a short excursion through some of the marsh sanctuaries of my memory. Rather than five conversations with Myron Bander, as you might infer from my title, I will describe vignettes involving five seminars where he and I were both present, and I will try to put some aspects of these brief scenes into a broader historical context. The places and dates of the five seminars are 1) Caltech 1974 (with Feynman & Gell-Mann) 2) Irvine 1976 (with Pellam) 3) Irvine 1976 (by Schonfeld) 4) Irvine 1977 (by Freedman) 5) Fort Lauderdale 2012 (by Myron) I met Myron when I was a graduate student at Caltech, in the early 1970s. Subsequently he hired me for my first post-doctoral position here at UC Irvine, 1976-78.

Those were important formative years for me, of course, and I greatly appreciated the opportunities, guidance, and advice that Myron gave me during this period.

As most of you know, in addition to being a first-rate theoretical physicist with excellent intuition, Myron had a well-developed, keen sense of humor that ranged from the darkly sardonic to the gleefully juvenile, but which was always poignant. I liked this very much. I will always remember him for it.

1)

The first vivid memory of Myron I can recall was at a Caltech seminar. He and Paul Thomas had driven up from Irvine to attend a talk, a trip I would later repeat many times with Myron after Paul moved on and I became his replacement as an Irvine postdoc.

Myron and Paul were seated against the outside wall in the 4th floor Lauritsen conference room, opposite from me across the conference table in the center of the room. The talk (I cannot remember who gave it, or what it was about!) was well-attended. All the seats were taken except for one chair at the head of the table closest to the screen. This was invariably taken by Murray Gell-Mann.

Richard Feynman came in just before Gell-Mann, and seeing no available seats, except for the one reserved for Murray, he dramatically lay down on the floor between the table and the row of seats where Myron and Paul sat. Thus he was out of sight to anyone subsequently entering the room.

Murray came in last, and seated himself with a groan and a displeased look. This caught everyone's attention, perhaps intentionally, and in explanation Murray said, "I'm suffering from rich man's disease." Immediately from the floor behind the table came a loud remark, "You mean a fat ass?" Completely surprised, Murray attempted to regain his composure and said something like, "Well ... no ... I think that would facilitate sitting." But the room had already erupted with laughter, including Myron's.

However, as I looked across the room, I could see Paul was not so much amused as he was perplexed. I suppose he did not believe that anyone would have the nerve to say something so coarse to Gell-Mann. It was immediately clear why. He turned to Myron, pointed his finger toward the unknown (to him, apparently) guy on the floor and mouthed the words "Who's that?!" I could see Myron respond, "That is Feynman."

2)

Later, I came to Irvine in the fall of 1976 as a postdoc in the particle theory group (with Myron, Gordon Shaw, and Dennis Silverman). In addition to carrying out research, the tasks given to me included organizing the preprint collection (articles arrived regularly in the US mail, printed on real paper!) and arranging some seminars. I can remember especially one talk attended by Professor John Pellam, in a wheel chair, during the later stages of his battle with ALS. I was seated behind him in the back of the room. John had been wheeled into the room by an attendant, and he wore weighted wrist bands to help hold his arms in their correct position on the armrests of the wheel chair. During the talk his right arm slipped off the armrest and fell down against the wheel of the chair. He couldn't raise it up.

I was distracted by this for the rest of the talk, concerned that Pellam's attendant would wheel him away at the end of the seminar without seeing that his arm was out of position, thereby injuring John as the chair's wheel turned against his skin. Even so, I waited until the end of the talk to offer any help. Then when the talk concluded, before I could do or say anything, Pellam quickly spoke to Myron who was seated next to him but who had not yet noticed John's predicament.

John said in a remarkably clear, firm voice, "Myron, would you please give me my arm?" Myron did so, gently but deliberately, with such great care and grace that the images of this simple kindness are something I have not forgotten. In a way that I cannot clearly express, this brief incident summarizes for me Myron's thoughtfulness and concern for others. From that moment onward I felt very comfortable working under his supervision.

(John Pellam died on July 23, 1977. There was a Pellam Symposium here at UC Irvine, 12 May 1978, in memoriam of the eminent low-temperature physicist. It was one of the last events I participated in at the Physics Department here, before I left to take up my second postdoctoral postion at the EFI, University of Chicago. As I recall, it is the only other such memorial symposium that I have attended. So now, some 35 years later, by participating in this symposium I feel like I have completed a long and circuitous journey back to its starting point.)

3)

When I came here in 1976 I knew many of the people at Caltech (having been a student there for so long!) and this furnished me with a ready pool of seminar speakers to invite to Irvine. At that time the theory group here provided a whopping $15 to entertain each speaker. This was usually enough to cover the cost of an abalone dinner for both me and the visitor at Woody's Wharf. Those were the days! I see that an abalone starter plate now goes for $22 at the Cannery (that is, if they have any).

Jonathan Schonfeld was then a new postdoc at Caltech, so I invited him to Irvine to talk about instantons in November 1976. He gave us an excellent seminar. Early in his talk, I asked an uninspired question about the axial anomaly. Myron joined in with Jonathan to provide a rapid response. I cannot remember the details, but their answers gave me the impression that they both thought I had asked a "dumb" question. Feeling somewhat uneducated as a consequence of this, I began to think about how to compute anomalies, more to make up for my own shortcomings than to do anything new.

Oddly enough, although I had worked on the Coleman-Weinberg mechanism of symmetry breaking for my thesis, I had never worked out from scratch the anomalous trace of the energy-momentum tensor. After dinner with Jonathan, I went back to the office before going home. I attempted a quick back-of-the-envelope calculation of the trace anomaly. Somewhat surprisingly, I got the known answer spot-on after a single line of algebra. That was fun. Now, it was well-known even at that time that the anomalous trace of the energy-momentum tensor caused scale and conformal invariance to be broken. Having also worked on radiative corrections for supersymmetric theories for my thesis, another obvious exercise suggested itself to me: Compute the spinor trace of the supercurrent which I knew to be related to broken super conformal invariance. As before, a similar one line calculation gave an anomalous answer. But this time, it was not only fun, it turned out to be a new result.

The next day I discussed this with Myron. He was enthusiastic, and he encouraged me to work through the problem in complete detail to justify my quick calculation. So I plunged into that analysis. Looking back, I am surprised that he did not join in the calculation, but he did not. Perhaps it was because supersymmetry was such a conceptual stretch given the physics known at that time. Unlike now.

In the meanwhile, I also had to finish up the last chapter of my thesis and return to Caltech to defend it all. (How I was a "postdoc" before defending my thesis is a long story for another venue. Suffice it to say Myron made it possible.) So there was a bit of delay before I completed all the details and checks for the supercurrent calculation. In fact it was about 6 months later that I recall asking Myron if I should write up the result and send it out as a preprint. I don't remember his exact words, but the upshot was, why yes, of course ... why as in why do you need to ask?! After all, wasn't that what I was hired to do? Find and publish new results! But Myron was more specific. He encouraged me to write a short version, and submit it quickly, and leave the bulk of the details until later. This was prescient.

I wrote it up and submitted it to the various preprint libraries, and to Physics Letters. Soon thereafter the members of the Irvine particle theory group dispersed for the summer. I went to SLAC. Myron went to Fermilab. A few days after going to Fermilab, he called me at SLAC to say that another preprint by Abbott, Grisaru, and Schnitzer on supercurrent anomalies had arrived at Fermilab. "Don't worry," Myron said, "your paper got here first." I think there was some institutional as well as individual competition underlying that statement.

4)

The superanomaly study led to some interesting developments. Regrettably, I never did write up all the details of my own calculations on the subject. I let it go at just that one Physics Letters article. Still, I did reap some benefits from the work. In particular, as a result I became friends with Dan Freedman and subsequently worked with him on another interesting problem.

Dan was at Caltech on sabbatical for the 1977-78 academic year. I met him there on one of my many trips back to Pasadena from Irvine. As I recall, I was introduced to Dan by John Schwarz, who described me as "yet another superman" or words to that effect. I replied that I was somewhat of an anomaly among supermen, a remark appreciated by Freedman. Later that fall I invited Dan to Irvine to give a seminar.

That seminar offered an opportunity for Myron to display his excellent intuition for theoretical physics. He did.

Freedman spoke about his research on antisymmetric tensor fields. The results had just appeared as a Caltech preprint and were very interesting.

==================== (my emphasis)

In his Irvine talk, Dan repeated his assertion that the 2 → 2 amplitude was of opposite sign to that in the nonlinear sigma model, and he therefore concluded that these two models were not equivalent.

It did not take long for Myron to object. He argued by an analogy with φ 4 theory, suggesting that only one sign of the amplitude would be consistent with a positive energy condition. And it was more or less clear from Dan's formalism that the energy density for his antisymmetric tensor model was positive. So, Myron concluded by suggesting that there must be a sign mistake in Dan's calculation of the 2 → 2 scattering amplitude.

Freedman said he would make some checks after he returned to Caltech. Eventually he did and, yes indeed, he found that his original sign was in error.

Kudos to Myron!

Let me re-emphasize that I describe this vignette as an example of Myron's excellent intuition, and not as criticism of Dan's original and very interesting preliminary paper. As every physicist knows, it is always possible to inadvertently flip a sign, especially when working alone without collaborative cross-checks. So far as I can tell, Dan's 1977 paper on the subject only appeared as a Caltech preprint. However, in a later final published paper co-authored with Paul Townsend, Freedman established the equivalence of the antisymmetric tensor theory to the sigma model. As is evident in the following excerpts from that later paper, it also helps to have good students.

5)

Now permit me to fast-forward 35 years to Myron's talk about de Sitter space at Miami 2012, the conference where he passed away the final night of the meeting, just after Carol had arrived to join him for a cruise of the Caribbean.

The slides from that talk are available online. Myron considered the supposedly conserved currents that generate the isometries of de Sitter space in light of quantum effects. The main point of his talk is summarized on this slide, with his usual humor.

In other words, he claimed these currents are anomalous. So, 35 years after showing interest and offering critical advice to me on related matters, I learned that his enthusiasm for this subject was undiminished. In fact, it had increased.

Myron backed up his claim about quantum corrections by considering the effects of a scalar field, a subject somewhat in vogue during the last year or so.

I will refer you to his talk and his published work cited therein for details.

Let me just use this to segue to some of my recent work, also on scalar fields, which perhaps may have some relevance to the other technical talks to follow at this symposium. Allow me to offer some brief remarks about scale versus conformal invariance in the context of galileon models, and to discuss some features of classical field configurations for these systems. This is something I have worked on during the last year or so with David Fairlie, one of the originators of galileon models.

I had hoped to interest Myron in this subject. But, sadly, circumstances did not permit that. Physics is not going to be as much fun without him.

Galileon theories are a class of models for hypothetical scalar …elds whose Lagrangians involve multilinears of …rst and second derivatives, but whose nonlinear …eld equations are nonetheless still only second order. These have been studied by numerous authors (please see references in the written version). The simplest example involves a single scalar …eld. This galileon …eld is usually coupled to all other matter through the trace of the energymomentum tensor, (matter) , and is thus gravitation-like by virtue of the similarity between this universal coupling and that of the metric g to (matter) in general relativity.

But surely, in a self-consistent theory, for the coupling to be truly universal, the galileon should also be coupled to its own energy-momentum trace, even in the ‡at spacetime limit.

Some consequences of this additional self-coupling are discussed in my papers with David Fairlie.

The action for the lowest non-trivial member of the galileon hierarchy can be written in various ways upon integrating by parts. Perhaps the most elegant of these is

A 2 = 1 2 Z d n x : (1)
where is the scalar galileon …eld, = @ (x) =@x , etc., and where repeated indices are summed using the Lorentz metric = diag (1; 1; 1;

). The …eld equation that follows from locally extremizing A 2 is

0 = A 2 = = E 2 [ ] ;
(2)

E 2 [ ] : (3)
As previously stated, the …eld equation is only second-order, albeit nonlinear.

It is straightforward to include in A 2 a covariant minimal coupling to a spacetime metric and hence to deduce a symmetric energy-momentum tensor. In the ‡at-space limit, the result is

(2) = + : (4)
This is veri…ed to be conserved upon using the …eld equation,

@ (2) = E 2 [ ] : (5)
An interesting wrinkle now appears:

(2) is not traceless.

Consequently, the usual form of the scale current, x 2) , is not conserved. On the other hand, the action (1) is homogeneous in and its derivatives, and is clearly invariant in n dimensions under the scale transformation (x) ! s (n 4)=3 (sx) for …eld con…gurations su¢ ciently localized so that boundary e¤ects are insigni…cant. Hence the corresponding Noether current must be conserved. This current is readily found, especially for n = 4, so let us restrict our attention to four spacetime dimensions in the following.

In that case the trace is easily identi…ed to be a total divergence:

(2)

(2) = @ ( ) :

That is to say, for n = 4 the virial is the trilinear V = . So a conserved scale current is given by the combination,

S = x (2) : (7)
Interestingly, this virial is not a divergence modulo a conserved current, so this model is not conformally invariant despite being scale invariant.

But our main interest was in the fact that the nonzero trace suggests an additional interaction where couples directly to its own (2) . This is similar to coupling a conventional massive scalar to the trace of its own energy-momentum tensor. In that previously considered example, however, the consistent coupling of the …eld to its trace required an iteration to all orders in the coupling. Upon summing the iteration and making a …eld rede…nition, the Nambu-Goldstone model emerged. But, for the simplest galileon model in four spacetime dimensions a consistent coupling of …eld and trace is much easier to implement. No iteration is required.

The …rst-order coupling alone is consistent, after integrating by parts and ignoring boundary contributions, so that

1 4 Z @ ( ) d 4 x = 1 4 Z d 4 x : (8)
(Similar quadrilinear terms have appeared previously, only multiplied by scalar curvature R so that they would drop out in the ‡at spacetime limit that we have emphasized here.)

Consistency follows because (8) gives an additional contribution to the energy-momentum tensor which is traceless, in 4D spacetime:

(3) = 1 4

;

(3) = 0 :

Of course, were quantum e¤ects taken into consideration, say by dimensional regularization, it is clear that this trace may well be anomalous. We leave such quantum e¤ects as (backof-the-envelope?) exercises for the reader, and proceed with a classical …eld perspective.

Based on these elementary observations, we considered a model with action

A = Z 1 2 1 2 1 4 d 4 x ; (10)
where for the Lagrangian L we take a mixture of three terms: the standard bilinear, the trilinear galileon, and its corresponding quadrilinear trace-coupling.

The …eld equation of motion for the model is now

0 = A= = E [ ],
where

E [ ] ( ) ( ) : (11)
As expected, this …eld equation is again second-order, albeit nonlinear.

Following the lead of many other studies, we discussed solutions of these classical equations in our papers, including some novel general relativistic e¤ects that arise from minimal covariant coupling to gravity. We were surprised to discover the model readily exhibits naked singularities, unlike the usual situation where only a bilinear scalar …eld action is coupled minimally to gravity. With the static metric expressed in Schwarzschild coordinates as (ds) 2 = e N (r) (dt) 2 e L(r) (dr) 2 r 2 (d ) 2 r 2 sin 2 (d') 2 ;

(12) we de…ne the constants M and C in the usual way in terms of r ! 1 data for static solutions in asymptotically ‡at spacetime:

e L=2 r!1 1 + M r + O 1 r 2 ; (13) e N=2 r!1 1 M r + O 1 r 2 ; (14) r!1 C r + O 1 r 2 : (15)
We found the following curve separating solutions with naked singularities from those with event horizons. This plot is based on a numerical study. Computed r ! 1 asymptotic data points (small circles) and an interpolating curve separating solutions with naked singularities from solutions with event horizons.

  the course of his talk, Dan exhibited a lowest order calculation of 2 → 2 scattering, as described in the preprint.

Acknowledgments

I would like to thank the organizers of this symposium for the opportunity to speak here.

It was an honor to be able to remember Myron in this forum. This work was supported in part by NSF Award PHY-1214521 and by a University of Miami Cooper Fellowship.