Bilel Bousselmi
email: bilel.bousselmi@insa-rennes.fr

Jean-François Dupuy
email: jean-francois.dupuy@insa-rennes.fr

Abderrazek Karoui
email: abderrazek.karoui@fsb.rnu.tn

Censored count data regression with missing censoring information

Keywords: Poisson regression, Asymptotic properties, Missing data, Regression calibration, Multiple imputation, Augmented inverse probability weighting, Simulations

We investigate estimation in the Poisson regression model when the count response is right-censored and the censoring indicators are missing at random. We propose several estimators based on the regression calibration, multiple imputation and augmented inverse probability weighting methods. Under appropriate regularity conditions, we prove the consistency of our estimators and we derive their asymptotic distributions. Simulation experiments are carried out to investigate the finite sample behaviour and relative performance of the proposed estimates.

Introduction

Poisson regression is a popular tool for modeling the relationship between a count response (such as the number of cases of a specific disease in epidemiology, or the number of insurance claims within a given period of time) and a set of predictors or covariates. Over the past years, Poisson regression has been extended to accommodate censored count data. Although censoring is usually associated to lifetime data analysis, count data can also be censored, the most common type being right-censoring, which occurs when it is only known that the true count is higher than the observed one. For example, consider a study investigating the smoking habits of some population, where people report their number of cigarettes smoked per day. If one possible answer is "20 cigarettes or more", all cigarettes counts greater than 20 are right-censored at 20. Ignoring censoring is known to yield biased estimates and thus, incorrect inferences. Statistical inference in the censored Poisson regression model and its extensions was therefore addressed by several authors ; see, for example, [START_REF] Terza | A Tobit-type estimator for the censored Poisson regression model[END_REF], [START_REF] Caudill | Modeling household fertility decisions: Estimation and testing of censored regression models for count data[END_REF], [START_REF] Famoye | Censored generalized Poisson regression model[END_REF], [START_REF] Xie | Diagnostics analysis in censored generalized Poisson regression model[END_REF], [START_REF] Mahmoud | On estimating parameters of censored generalized Poisson regression model[END_REF] for censored generalized Poisson regression, [START_REF] Karlis | Finite mixtures of censored Poisson regression models[END_REF] for finite mixtures of censored Poisson regressions and [START_REF] Saffari | Zero-inflated Poisson regression models with right censored count data[END_REF] and Nguyen andDupuy (2019, 2020) for zero-inflated censored Poisson regression. Censored models for count data can be conveniently specified by introducing a censoring indicator which is set to 1 if the observed count is not censored and 0 otherwise. In this paper, we consider the situation where the censoring indicator is missing for some sample individuals. In the context of survival analysis, this issue has been considered by several authors. For example, [START_REF] Van Der Laan | Efficient estimation from right-censored data when failure indicators are missing at random[END_REF] and [START_REF] Subramanian | Survival analysis for the missing censoring indicator model using kernel density estimation techniques[END_REF][START_REF] Subramanian | Multiple imputations and the missing censoring indicator model[END_REF] address estimation of the survival function of a random survival time with missing censoring indicators. [START_REF] Mckeague | Product-limit estimators and Cox regression with missing censoring information[END_REF] consider estimation in Cox proportional hazards regression model with missing completely at random censoring indicators. [START_REF] Wang | Hazard function estimation with cause-of-death data missing at random[END_REF] and [START_REF] Brunel | Nonparametric estimation for survival data with censoring indicators missing at random[END_REF] propose various nonparametric estimates of the hazard and conditional hazard functions with censoring indicators missing at random. [START_REF] Wang | Linear regression analysis of survival data with missing censoring indicators[END_REF] investigate estimation in the linear regression model for survival data with missing censoring indicators. Estimation in censored Poisson regression with missing censoring information is still an open problem. Our aim in this paper is to provide and compare several estimates adapted to this setting.

Missing data problems have given rise to a rich literature and several estimation methods have been proposed for this setting. A common and simple approach is to exclude individuals with missing data. This is the so-called complete-case analysis. This method can induce bias and substantial variance increase. Two alternatives are regression calibration and multiple imputation. In regression calibration, missing data are replaced by their conditional expectation given the observed data. In multiple imputation, missing data are replaced by data generated from an imputation model. This imputation is repeated M times, generating M completed data sets. Each of them is analysed and an overall estimator is obtained by combining the estimates of the M completed samples. Both methods require a model for the missing data given the observed data. Inverse probability weighting constitutes another alternative method for dealing with missing data (see for example [START_REF] Seaman | Review of inverse probability weighting for dealing with missing data[END_REF] for a review of this method). Similarly to the complete-case analysis, inverse probability weighting only uses complete cases, but weights are used to rebalance the set of complete cases. Calculating these weights requires a model for the probability that an individual has complete data. Augmented inverse probability weighting was then proposed to ensure robustness against misspecification of the missingness model (see, for example, [START_REF] Tsiatis | Semiparametric Theory and Missing Data[END_REF] for a detailed account on the method).

In this paper, we investigate, both theoretically and numerically, the regression calibration, multiple imputation and augmented inverse probability weighting estimators of the regression parameter in the censored Poisson regression model with missing censoring indicators. Our analysis of these estimates will be based on parametric assumptions for the conditional models for missing data and the missingness mechanism. The plan of the paper is as follows. In Section 2, we describe the model setup and we introduce the notations that will be used throughout the paper. In Section 3, we introduce our regression calibration estimator and we establish its consistency and asymptotic normality. In Sections 4 and 5, we propose our multiple imputation and augmented inverse probability weighted estimators, and we derive their asymptotic properties. All our theoretical derivations are based on an incomplete gamma function formulation of the distribution function of the Poisson regression model. Consistent asymptotic variance estimates are also proposed for the regression calibration, multiple imputation and augmented inverse probability weighted estimators. In Section 6, we conduct a simulation study to assess the finite sample performance and robustness to parametric assumptions of the proposed estimates. Discussion and perspectives are given in Section 7. All proofs are deferred to appendices.

Model, data, notations

Let Y denote the count of interest and X = (1, X 2 , . . . , X p ) be a p-vector of covariates ( denotes the transpose operator). We assume that the conditional distribution of Y given X is given by a Poisson regression model with parameter λ = exp(β X), where β ∈ R p is a vector of unknown parameters.

We consider the situation where Y can be right-censored, that is, instead of the true Y , we eventually observe a value which is smaller than Y . This can be formalised by introducing a finite random variable C such that we observe either Assume that n independent individuals are available and that for each of them, we observe the triplet (Y * i , X i , δ i ) (with i ∈ {1, . . . , n}). Based on these observations, the likelihood of β is calculated as:

Y if Y < C or C if Y ≥ C,
L n (β) = n i=1 P(Y i = Y * i |X i ) δ i P(Y i ≥ Y * i |X i ) 1-δ i ,
from which we easily deduce the loglikelihood n (β) = log L n (β):

n (β) = n i=1    δ i Y * i β X i -e β X i -log(Y * i !) + (1 -δ i ) log   1 - Y * i -1 k=0 e -exp(β X i )+kβ X i k!      (2.1)
By standard asymptotic theory, the maximum likelihood estimator βn = arg max β n (β) is consistent and asymptotically normal with variance -E[∂ 2 1 (β)/∂β∂β ]. Now, we consider the situation where some additional uncertainty can arise in the observations. Precisely, we consider the situation where the censoring indicator δ i is missing for some individuals. Let ξ be a missingness indicator, that is, ξ = 1 if δ is observed and ξ = 0 otherwise. Then, for individual i ∈ {1, . . . , n}, the observed data are

(Y * i , X i , δ i , ξ i = 1) or (Y * i , X i , ξ i = 0). (2.2)
We consider a missing at random (MAR) mechanism, which means that ξ and δ are independent given all other observed variables (a more restrictive assumption is that ξ and δ are independent, which is called "missing completely at random"). In the next sections, we propose, investigate and compare several estimators of β in this context.

Regression calibration estimation

The proposed estimator

Our first estimator is based on the regression calibration idea. It consists in replacing any missing δ i in (2.1) by its conditional expectation E(δ i |W i ), where W i contains the observed variables Y * i and X i and eventually (if available) some observed surrogate variables V i for δ i . Thus, we let

W i = (Y * i , X i , V i ) (we denote by q the dimension of W i
). An approximated version of δ i can then be defined as:

δi = ξ i δ i + (1 -ξ i )E(δ i |W i ).
The conditional expectation E(δ i |W i ) (or conditional probability P(δ i = 1|W i )) will generally be unknown and will have to be estimated. As is usual with the regression calibration approach, we assume that E(δ i |W i ) can be specified by a parametric model m(W i , θ), where θ is an unknown q-dimensional parameter with true value θ 0 . Remark 1. A convenient candidate for m(•, •) is the logistic regression model m(W i , θ) = logit -1 (θ W i ) but other choices, such as the probit, are possible. One may also allow for polynomial, spline and interaction terms in these models, in order to make them as flexible as desired. In what follows, we assume a general model m(W i , θ) with some regularity conditions stated in section 3.2.

At a first stage, we estimate θ 0 by maximizing a likelihood based on complete cases i ∈ {1, . . . , n|ξ i = 1} only:

θn = arg max θ n i=1 m(W i , θ) ξ i δ i (1 -m(W i , θ)) ξ i (1-δ i ) . (3.3) Let ṁ(W i , θ) = ∂m(W i , θ) ∂θ , m i (θ) = ṁ(W i , θ) m(W i , θ)(1 -m(W i , θ)) ,
and

Θ(θ) = E ṁ⊗2 (W, θ) m(W, θ)(1 -m(W, θ)) ξ ,
where for any column vector u, u ⊗2 = uu . Then it is rather straightforward to see that θn is asymptotically linear with influence function

Θ -1 (θ 0 ) m i (θ 0 )ξ i (δ i -m(W i , θ 0 )), that is : √ n( θn -θ 0 ) = 1 √ n n i=1 Θ -1 (θ 0 ) m i (θ 0 )ξ i (δ i -m(W i , θ 0 )) + o P (1). (3.4)
Finally, it will be useful to note that if Y is distributed as Poisson with parameter λ, then for any u ∈ N, P(Y ≤ u) = u k=0 exp(-λ)λ k /k! = Γ(u + 1, λ)/u! where Γ(u, λ) = ∞ λ t u-1 exp(-t)dt is the incomplete gamma function, whose derivative with respect to λ is given by ∂Γ(u, λ)/∂λ = -exp(-λ)λ u-1 . Now, letting δi (θ) = ξ i δ i +(1-ξ i )m(W i , θ) be the approximation of δ i based on model m(W i , θ), we define our regression calibration estimator of β as

βn = arg max β n (β, θn ), where n (β, θn ) = n i=1    δi ( θn ) Y * i β X i -e β X i -log(Y * i !) + (1 -δi ( θn )) log   1 - Γ Y * i , e β X i (Y * i -1)!     
is an approximated version of (2.1).

Regularity conditions and asymptotic results

The following regularity conditions are needed to establish the asymptotic properties of the regression calibration estimator. We assume: C1 The covariates vectors X i and V i are bounded, for every i = 1, 2, . . . Before stating the asymptotics of βn , we introduce some further notations. Let h β be the function defined by:

h β (y, x) =
e -e β x +β xy (y -1)! -Γ(y, e β x ) (3.5)

for any β ∈ R p , x ∈ R p and y ∈ N\{0}. Let also π(W) = P(ξ = 1|W) and define the matrices

Σ 1 (β) = E XX δe β X + (δ -1) Y * -e β X -h β (Y * , X) h β (Y * , X) , Σ 2 (β, θ) = E X ṁ (W, θ) Y * -e β X -h β (Y * , X) (1 -π(W)) , Σ 3 (β, θ) = E X ṁ (W, θ) Y * -e β X -h β (Y * , X) .
We are now in position to state our first theorem. The proof is given in Appendix A.

Theorem 3.1. Assume that conditions C1-C4 hold. Then βn P -→ β 0 as n → ∞ and √ n( βn -β 0 ) is asymptotically normal with mean zero and variance Σ, where

Σ = Σ -1 1 (β 0 ) Σ 1 (β 0 ) + (2Σ 3 (β 0 , θ 0 ) -Σ 2 (β 0 , θ 0 )) Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ) Σ -1 1 (β 0 ).
Remark 3. If π(W) is identically equal to 1 (that is, if there is no missing data), Σ reduces to the asymptotic variance of the maximum likelihood estimator βn in (2.1), which in turn reduces to the usual asymptotic variance (E[XX e β 0 X ]) -1 in Poisson regression if m(W, θ 0 ) is identically equal to 1 (that is, no censoring can affect the data).

A consistent estimator of Σ is given by

Σ -1 1,n ( βn , θn ) Σ 1,n ( βn , θn ) + 2Σ 3,n ( βn , θn ) -Σ 2,n ( βn , θn ) Θ -1 n ( θn )Σ 2,n ( βn , θn ) Σ -1 1,n ( βn , θn ),
where

Σ 1,n (β, θ) = 1 n n i=1 X i X i δi (θ)e β X i + ( δi (θ) -1) Y * i -e β X i -h β (Y * i , X i ) h β (Y * i , X i ) , Σ 2,n (β, θ) = 1 n n i=1 X i ṁ (W i , θ) Y * i -e β X i -h β (Y * i , X i ) (1 -ξ i ), Σ 3,n (β, θ) = 1 n n i=1 X i ṁ (W i , θ) Y * i -e β X i -h β (Y * i , X i ) , Θ n (θ) = 1 n n i=1 ṁ⊗2 (W i , θ) m(W i , θ)(1 -m(W i , θ)) ξ i .
The consistency proof of the variance estimator uses similar arguments as the proof of consistency of βn , it is thus omitted. The estimator βn will be evaluated in the simulation study of Section 6.

Several methods have been proposed to address missing data problems in regression. Among them is the multiple imputation, which provides an alternative, popular and widely-used approach. The basic idea is to create several (say M ) completed data sets, by filling in plausible values for the missing data. Then, each filled sample is analysed as if it were the complete data set. Finally, the M imputed-samples inferences are combined into a single overall inference. In the next section, we investigate this approach for estimating β in our problem.

Multiple imputation

In this section, we assume, as in Section 3, that the conditional expectation E(δ i |W i ) can be specified by a parametric model m(W i , θ 0 ), and we denote by θn the maximum likelihood estimator of θ 0 based on the complete cases i ∈ {1, . . . , n|ξ i = 1}.

The imputation procedure is as follows. Each missing δ i is replaced by a random draw from the Bernoulli distribution B(m(W i , θn )). We obtain a completed data set. This procedure is repeated M times to form M imputed data sets. For a given θ, let D i,j (θ) ∼ B(m(W i , θ)) denote the imputation of δ i in the j-th completed data set (j = 1, . . . , M ). Let also

δ * i,j (θ) = ξ i δ i + (1 -ξ i )D i,j (θ)
be the random variable which is equal to

δ i is ξ i = 1 (that is, if δ i is observed) and to D i,j (θ) if ξ i = 0 (that is, if δ i is missing) (note
the difference between the imputation method, where δ * i,j (θ) ∈ {0, 1}, and the regression calibration approach, where δi (θ)

∈ [0, 1]). A single-imputation estimator β * n,j of β 0 is obtained by maximizing the imputed log-likelihood * n,j (β, θn ) = n i=1    δ * i,j ( θn ) Y * i β X i -e β X i -log(Y * i !) + (1 -δ * i,j ( θn )) log   1 - Γ Y * i , e β X i (Y * i -1)!      .
The final multiple imputation estimator β * n is obtained by averaging the M estimators β * n,j as:

β * n = 1 M M j=1 β * n,j .
The next theorem gives the asymptotic properties of β * n . Its proof is given in Appendix B.

Theorem 4.1.

For j = 1, . . . , M , let f β,θ,j (O i ) = X i {δ * i,j (θ)[Y * i -e β X i -h β (Y * i , X i )]+h β (Y * i , X i )}, where O i denotes the observation (2.2). Let also Σ * 1 (β, θ) = var( 1 M M j=1 f β,θ,j (O 1 ). If conditions C1-C4 hold, then β * n P -→ β 0 as n → ∞ and √ n( β * n -β 0 )
is asymptotically normal with mean zero and variance Σ * , where

Σ * = Σ -1 1 (β 0 ) Σ * 1 (β 0 , θ 0 ) + (2Σ 3 (β 0 , θ 0 ) -Σ 2 (β 0 , θ 0 )) Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ) Σ -1 1 (β 0 ).
A consistent estimator of Σ * can be obtained as

Σ-1 1,n ( β * n , θn ) Σ * 1,n ( β * n , θn ) + 2Σ 3,n ( β * n , θn ) -Σ 2,n ( β * n , θn ) Θ -1 n ( θn )Σ 2,n ( β * n , θn ) Σ-1 1,n ( β * n , θn ),
where

Σ * 1,n (β, θ) is the empirical covariance of the vectors 1 M M j=1 f β,θ,j (O i ) (i = 1, . . . , n), Σ1,n (β, θ) is the average 1 M M j=1 Σ 1,n,j (β, θ), with Σ 1,n,j (β, θ) = 1 n n i=1 X i X i δ * i,j (θ)e β X i + (δ * i,j (θ) -1) Y * i -e β X i -h β (Y * i , X i ) h β (Y * i , X i ) ,
and Σ 2,n , Σ 3,n and Θ n are as given in Section 3.

Regression calibration and multiple imputation rely on the ability of the investigator to formulate an appropriate model for E(δ|W). Misspecifying this model is likely to yield biased estimates of the parameters of interest. An alternative approach is to specify the selection probabilities π(W i ) = P(ξ i = 1|W i ) and to use the inverse probability weighting (IPW) of complete-case technique of [START_REF] Horvitz | A generalization of sampling without replacement from a finite universe[END_REF]. The basic idea of IPW is to adjust a complete-case analysis by weighting individuals with no missing data by the inverse of their selection probability. Selection probabilities are generally unknown and have to be estimated. Again, misspecifying the π(W i ), i = 1, . . . , n is likely to yield biased inference. Moreover, by discarding individuals with missing data, IPW is also known to yield loss of efficiency. For these reasons, the augmented IPW approach (AIPW henceforth, see [START_REF] Robins | Estimation of regression coefficients when some regressors are not always observed[END_REF] was proposed to improve the basic IPW. Since its introduction, the method has been shown to be doubly robust in several models, such as the proportional hazards model [START_REF] Wang | Augmented inverse probability weighted estimator for Cox missing covariate regression[END_REF]), the single-index model [START_REF] Guo | Empirical likelihood for single index model with missing covariates at random[END_REF], the additive hazards model [START_REF] Sun | Analysis of two-phase sampling data with semiparametric additive hazards models[END_REF] and the accelerated failure time model [START_REF] Steingrimsson | Estimation in the semiparametric accelerated failure time model with missing covariates: improving efficiency through augmentation[END_REF]. Double robustness refers to the fact that the AIPW estimates are consistent as long as either the selection probability model or the conditional expectation of the missing data is correctly specified. In the next section, we propose an augmented IPW estimating equation adapted to our problem, and we investigate the asymptotic properties of the resulting estimator.

Augmented inverse probability weighted estimation

Inspired by [START_REF] Horvitz | A generalization of sampling without replacement from a finite universe[END_REF], the inverse probability weighting of complete cases has become a classical estimation method in missing data problems. One drawback of the method is that the observed variables of subjects with missing data are not fully used, except through the estimation of the unknown selection probabilities. The AIPW method improves IPW by introducing an additional term involving contributions from individuals with some missing data (we refer to [START_REF] Tsiatis | Semiparametric Theory and Missing Data[END_REF] for a detailed account on the method and numerous references). Adapting this idea, we propose the following augmented IPW estimating equation for β:

n i=1 X i ξ i δ i π(W i ) + 1 - ξ i π(W i ) E(δ i |W i ) Y * i -e β X i -h β (Y * i , X i ) + h β (Y * i , X i ) .
The quantities E(δ i |W i ) and π(W i ) are unknown and have to be estimated. We assume that they can be specified by some parametric models m(W i , θ) and π(W i , γ) respectively, where θ and γ are unknown q-dimensional parameters with true values θ 0 and γ 0 . Let θn and γn be the maximum likelihood estimates of θ 0 and γ 0 . θn is given by (3.3). Similarly, γn can be obtained as

γn = arg max γ n i=1 π(W i , γ) ξ i (1 -π(W i , γ)) 1-ξ i .
Finally, our AIPW estimator βn of β solves the estimating equation ˘ n (β, θn , γn ) = 0, where

˘ n (β, θn , γn ) = n i=1 X i ξ i δ i π(W i , γn ) + 1 - ξ i π(W i , γn ) m(W i , θn ) Y * i -e β X i -h β (Y * i , X i ) + h β (Y * i , X i ) .
Before stating the asymptotic properties of βn , we introduce some further notations and regularity conditions. For any θ, γ ∈ R q , we let

δi (θ, γ) = ξ i δ i π(W i , γ) + 1 - ξ i π(W i , γ) m(W i , θ).
Assuming the parametric model π(W i , γ) for the selection probabilities, the maximum likelihood estimator γn is asymptotically linear with influence function

Σ -1 4 (γ 0 ) π i (γ 0 )(ξ i -π(W i , γ 0 )), where π(W i , γ) = ∂π(W i , γ) ∂γ , π i (γ) = π(W i , γ) π(W i , γ)(1 -π(W i , γ)) ,
and

Σ 4 (γ) = E π⊗2 (W, γ) π(W, γ)(1 -π(W, γ))
.

That is :

√ n(γ n -γ 0 ) = 1 √ n n i=1 Σ -1 4 (γ 0 ) π i (γ 0 )(ξ i -π(W i , γ 0 )) + o P (1). (5.6)
If the models m(W i , θ) and π(W i , γ) are misspecified, then by [START_REF] White | Maximum Likelihood Estimation of Misspecified Models[END_REF], there exists θ * and γ * such that θn P → θ * and γn P → γ * . Moreover, the asymptotic linear expansions for θn and γn are given by (3.4) and (5.6), with θ 0 and γ 0 replaced by θ * and γ * respectively. If the model m(W i , θ)

(respectively π(W i , γ)) is correctly specified, then θ * = θ 0 (respectively γ * = γ 0 ). Finally, let Σ 5 (β, θ, γ) = E X Y * -e β X -h β (Y * , X) 1 - ξ π(W, γ) ṁ (W, θ) , Σ 6 (β, θ, γ) = E X Y * -e β X -h β (Y * , X) ξ π (W, γ) π 2 (W, γ) (m(W, θ) -δ) , Σ 7 (β, θ, γ) = Σ 1 (β) + (2Σ 3 (β, θ) -Σ 5 (β, θ, γ)) Θ -1 (θ)Σ 5 (β, θ, γ),
and

Σ 8 (β, θ, γ) = Σ 1 (β) -Σ 6 (β, θ, γ)Σ -1 4 (γ)Σ 6 (β, θ, γ),
We assume the following additional regularity conditions:

C5 The parameter space for γ is a bounded set G ⊂ R q and the true parameter value γ 0 lies in the interior of G.

C6

The function π(w, γ) is strictly greater than 0 for all value of w in the support of W and all γ ∈ G.

C7

The function π(w, γ) is differentiable with respect to γ, for every w.

For every γ, γ ∈ G, |π(w, γ) -π(w, γ)| ≤ g(w) γ -γ for some bounded function g with E[g(W)] = u.
Conditions C5 and C7 for γ and π(•, •) are similar to conditions C2 and C4 for θ and m(•, •). We are now in position to state the asymptotic properties of our AIPW estimator of β. Proofs are given in Appendix C (consistency) and D (asymptotic normality).

Theorem 5.1. Assume that conditions C1-C7 hold. If either or both of the models m(W i , θ) and π(W i , γ) are well specified, then βn

P -→ β 0 as n → ∞.
From this result, the proposed estimator βn is doubly robust, in the sense that it estimates consistently β 0 as long as one of m(W i , θ) and π(W i , γ) is correctly modeled. The next theorem describes the asymptotic distribution of βn .

Theorem 5.2. Assume that conditions C1-C7 hold. Then, as n → ∞, √ n( βn -β 0 ) converges in distribution to the Gaussian random vector N (0, J), where

J =            Σ -1 1 (β 0 )Σ 7 (β 0 , θ 0 , γ * )Σ -1 1 (β 0 ) if m(W i , θ) is correctly specified, Σ -1 1 (β 0 )Σ 8 (β 0 , θ * , γ 0 )Σ -1 1 (β 0 ) if π(W i , γ) is correctly specified, Σ -1 1 (β 0 ) if both m(W i , θ
) and π(W i , γ) are correctly specified.

In order to estimate the asymptotic variance of βn , let:

Σ1,n (β, θ, γ) = 1 n n i=1 X i X i δi (θ, γ)e β X i + ( δi (θ, γ) -1) Y * i -e β X i -h β (Y * i , X i ) h β (Y * i , X i ) , Σ 4,n (γ) = 1 n n i=1 π⊗2 (W i , γ) π(W i , γ)(1 -π(W i , γ)) , Σ 5,n (β, θ, γ) = 1 n n i=1 X i Y * i -e β X i -h β (Y * i , X i ) 1 - ξ i π(W i , γ) ṁ (W i , θ), Σ 6,n (β, θ, γ) = 1 n n i=1 X i Y * i -e β X i -h β (Y * i , X i ) ξ i π (W i , γ) π 2 (W i , γ) (m(W i , θ) -δ i ) , Σ 7,n (β, θ, γ) = Σ1,n (β, θ, γ) + (2Σ 3,n (β, θ) -Σ 5,n (β, θ, γ)) Θ -1 n (θ)Σ 5,n (β, θ, γ),
and

Σ 8,n (β, θ, γ) = Σ1,n (β, θ, γ) -Σ 6,n (β, θ, γ)Σ -1 4,n (γ)Σ 6,n (β, θ, γ),
where Σ 3,n and Θ n are as given in Section 3. Then a consistent estimator of J is given by:

J n =              Σ-1 1,n ( βn , θn , γn )Σ 7,n ( βn , θn , γn ) Σ-1 1,n ( βn , θn , γn ) if m(W i , θ) is correctly specified, Σ-1 1,n ( βn , θn , γn )Σ 8,n ( βn , θn , γn ) Σ-1 1,n ( βn , θn , γn ) if π(W i , γ) is correctly specified, Σ-1 1,n ( βn , θn , γn ) if both m(W i , θ) and π(W i , γ) are correctly specified.
The proof is omitted.

A simulation study

In this section, we evaluate and compare the finite sample performance of the regression calibration (RC), multiple imputation (MI) and AIPW estimators. The simulation design is as follows.

For each of n individuals, the count response Y is simulated from a Poisson regression model with parameter

λ = exp(β 1 + β 2 X 2 + β 3 X 3 + β 4 X 4 + β 5 X 5 ),
where β = (0.2, -0.1, 0.4, 0.3, 0.5), X 2 ∼ N (0, 1), X 3 ∼ Bernoulli(0.3), X 4 ∼ N (0, 1.5) and X 5 ∼ uniform [2,5]. The censoring and missingness mechanisms are set to be logit(m(W, θ)

) = θ 1 +θ 2 X 2 + θ 3 X 3 +θ 4 X 4 +θ 5 X 5 +θ 6 Y and logit(π(W, γ)) = γ 1 +γ 2 X 2 +γ 3 X 3 +γ 4 X 4 +γ 5 X 5 +γ 6 Y * respectively.
In a first experiment, θ and γ are chosen so that approximately 20% of the subjects have a censored count response and 20% have a missing censoring indicator. In a second experiment, approximately 20% of the subjects have a censored count and 40% have a missing censoring indicator. We take a sample size n = 500.

We compare the RC, MI and AIPW estimates under three scenario : (i) only m(W, θ) is correctly modeled, (ii) only π(W, γ) is correctly modeled, (iii) both m(W, θ) and π(W, γ) are correctly modeled. In the first scenario, π(W, γ) is incorrectly modeled as logit(π(W, γ)

) = γ 1 + γ 2 X 2 + γ 3 X 3 + γ 4 Y * . In the second scenario, m(W, θ) is incorrectly modeled as logit(m(W, θ)) = θ 1 + θ 2 X 2 + θ 3 X 3 + θ 4 Y * .
Our simulation results are based on N = 1000 simulated samples. For each estimator, we report the average bias, average standard error (SE), empirical root mean square error (RMSE) and empirical coverage probability (CP) of 95%-level confidence intervals. MI estimates are obtained with M = 50 (from our numerical experiments, this is large enough to ensure stability of the estimates). To establish a benchmark for comparisons, we also include an estimator based on the full data set with no missing censoring indicators and the complete-case (CC) estimator which maximizes the log-likelihood (2.1) on the subsample of complete cases only.

The results are summarized in Table 1 (experiment 1) andTable 2 (experiment 2). In the first scenario, the RC, MI and AIPW methods appear to have similar performance, with slightly smaller average SE and RMSE for the AIPW estimates. Coverage probabilities are close to the nominal confidence level, indicating that the asymptotic variances are appropriately estimated. In the second scenario, the bias of the RC and MI estimates increase substantially, resulting in coverage probabilities smaller than desired. This was expected due to misspecification of model m(W, θ). On the other hand, the bias of the AIPW estimate stays moderate, and is similar to the first scenario, which is also expected due to the double robustness property stated in Theorem 5.1. The AIPW method generally achieves the smallest SE and RMSE but when the proportion of missing data increases, its coverage probabilities decrease around 85%-90%, which might indicate that the asymptotic variance is slightly under-estimated. Finally, when both models are correct, all three methods perform similarly (results for the RC and MI methods are the same as for the first scenario). Overall, this simulation study confirms the theoretical results stated in the previous sections. The regression calibration, multiple imputation and robust IPW methods provide similar results when either m(W, θ) or both m(W, θ) and π(W, γ) are correctly specified. When m(W, θ) is misspecified, the AIPW approach performs better than RC and MI in terms of point estimation (with substantially smaller bias for AIPW) but MI seems to provide better coverage probabilities (at the price of larger variance estimates). The CC estimates are outperformed by the three methods in all scenarios. Unreported simulation results with other sample sizes and censoring fractions provide similar observations.

Discussion

In this article, we have investigated several estimators of the regression parameter of the censored Poisson regression model when censoring indicators are partially missing. The regression calibration and multiple imputation estimates and their asymptotic variance estimators lead to reliable inferences when the model for the missing data given the observed variables is correctly specified, while the augmented inverse probability weighted estimator is robust against misspecification of either the model for the missing data or the missingness mechanism. Our estimators rely on parametric assumptions for the conditional models for the missing data and missingness mechanism. It is now important to assess the sensitivity of the statistical inference based on these estimates, with regard to deviations to the model assumptions. An alternative estimation strategy may rely on semiparametric or nonparametric estimation of the models for the missing data and mechanism. This stimulating topic is the subject for our next research.

Appendix A: Proof of Theorem 3.1

Consistency. The consistency of βn can be proved by verifying the conditions of the inverse function theorem [START_REF] Foutz | On the unique consistent solution to the likelihood equations[END_REF]. We describe the main steps of the proof and omit calculation details.

Let ˙ n (β, θ) := ∂ n (β, θ)/∂β. Straightforward calculations yield:

˙ n (β, θ) = n i=1 X i δi (θ) Y * i -e β X i + (1 -δi (θ))h β (Y * i , X i ) ,
where h β (y, x) is given by (3.5). We first need to show that ∂ ˙ n (β, θn )/∂β exists and is continuous in a neighborhood of β 0 . The map β → ˙ n (β, θn ) is trivially differentiable with respect to β and its derivative is given by:

∂ ˙ n (β, θn ) ∂β = n i=1 X i X i -δi ( θn )e β X i + (1 -δi ( θn )) Y * i -e β X i -h β (Y * i , X i ) h β (Y * i , X i ) ,
which is continuous in β.

Secondly, we need to show that n -1 ˙ n (β 0 , θn ) = o P (1). To see this, we decompose n -1 ˙ n (β 0 , θn ):

1 n ˙ n (β 0 , θn ) = 1 n ˙ n (β 0 , θn ) - ˙ n (β 0 , θ 0 ) + 1 n ˙ n (β 0 , θ 0 ).
By the weak law of large numbers, n -1 ˙ n (β 0 , θ 0 ) converges in probability to

E X δ(θ 0 )(Y * -e β 0 X ) + (1 -δ(θ 0 ))h β 0 (Y * , X) = E X E( δ(θ 0 )|W)(Y * -e β 0 X ) + (1 -E( δ(θ 0 )|W))h β 0 (Y * , X) , (7.7)
where the second line follows by taking the conditional expectation given W. Under the missing at random assumption,

E( δ(θ 0 )|W) = E(ξδ + (1 -ξ)E(δ|W)|W) = E(ξ|W)E(δ|W) + (1 -E(ξ|W))E(δ|W) = E(δ|W).
Therefore, (7.7) is equal to

E X E(δ|W)(Y * -e β 0 X ) + (1 -E(δ|W))h β 0 (Y * , X) = E X δ(Y * -e β 0 X ) + (1 -δ)h β 0 (Y * , X) ,
which is equal to 0 (this can be seen by taking successively the conditional expectations given {δ = 1} and X). Convergence to 0 of n -1 ( ˙ n (β 0 , θn ) -

˙ n (β 0 , θ 0 )
) is a consequence of the consistency of θn and of assumptions C1, C2, C4. Details are omitted.

Thirdly, we need to show that n -1 ∂ ˙ n (β, θn )/∂β converges in probability to a fixed matrix, uniformly in an open neighborhood of β 0 . We have:

1 n ∂ ˙ n (β, θ) ∂β = 1 n n i=1 X i X i -δi (θ)e β X i + (1 -δi (θ)) Y * i -e β X i -h β (Y * i , X i ) h β (Y * i , X i ) .
We proceed as above and decompose n -1 ∂ ˙ n (β, θn )/∂β as

1 n ∂ ˙ n (β, θn ) ∂β = 1 n ∂ ˙ n (β, θn ) ∂β - ∂ ˙ n (β, θ 0 ) ∂β + 1 n ∂ ˙ n (β, θ 0 ) ∂β .
The first term converges to 0 (by the consistency of θn and assumptions C1, C2, C4) and n -1 ∂ ˙ n (β, θ 0 )/ ∂β converges in probability to -Σ 1 (β) (by the weak law of large numbers). Therefore, n -1 ∂ ˙ n (β, θn )/ ∂β converges in probability to -Σ 1 (β). Under conditions C1 and C2, the derivative of n -1 ∂ ˙ n (β, θn )/ ∂β with respect to β is bounded, for every n. Hence the sequence (n

-1 ∂ ˙ n (β, θn )/∂β ) n is equicon- tinuous. It follows from Ascoli theorem that the convergence of n -1 ∂ ˙ n (β, θn )/∂β to -Σ 1 (β) is uniform around β 0 .
Having proved the conditions of the inverse function theorem, we conclude that βn converges in probability to β 0 . Asymptotic normality. A Taylor's expansion of ˙ n ( βn , θn ) around (β 0 , θ 0 ) yields

√ n( βn -β 0 ) = - 1 n ∂ ˙ n (β 0 , θ 0 ) ∂β -1 1 √ n ˙ n (β 0 , θ 0 ) + 1 n ∂ ˙ n (β 0 , θ 0 ) ∂θ √ n( θn -θ 0 ) + o P (1).
We have:

1 n ∂ ˙ n (β, θ) ∂θ = 1 n n i=1 X i Y * i -e β X i -h β (Y * i , X i ) (1 -ξ i ) ṁ (W i , θ) = Σ 2 (β, θ) + o P (1).
Combining this and (3.4), we can write:

√ n( βn -β 0 ) = - 1 n ∂ ˙ n (β 0 , θ 0 ) ∂β -1 1 √ n n i=1 X i δi (θ 0 )(Y * i -e β 0 X i ) + (1 -δi (θ 0 ))h β 0 (Y * i , X i ) + Σ 2 (β 0 , θ 0 )Θ -1 (θ 0 ) m i (θ 0 )ξ i (δ i -m(W i , θ 0 )) + o P (1) := - 1 n ∂ ˙ n (β 0 , θ 0 ) ∂β -1 1 √ n n i=1 U i + o P (1). Now, note that var X i δi (θ 0 )(Y * i -e β 0 X i ) + (1 -δi (θ 0 ))h β 0 (Y * i , X i ) = Σ 1 (β 0 ), and var Σ 2 (β 0 , θ 0 )Θ -1 (θ 0 ) m i (θ 0 )ξ i (δ i -m(W i , θ 0 )) = Σ 2 (β 0 , θ 0 )Θ -1 (θ 0 )E m i ⊗2 (θ 0 )ξ i (δ i -m(W i , θ 0 )) 2 ×Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ) = Σ 2 (β 0 , θ 0 )Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ),
since under the missing at random assumption, we have :

E m i ⊗2 (θ 0 )ξ i (δ i -m(W i , θ 0 )) 2 = E ṁ⊗2 (W i , θ 0 ) {m(W i , θ 0 )(1 -m(W i , θ 0 ))} 2 E ξ i (δ i -m(W i , θ 0 )) 2 |W i = E ṁ⊗2 (W i , θ 0 ) {m(W i , θ 0 )(1 -m(W i , θ 0 ))} 2 E [ξ i |W i ] E [δ i -2δ i m(W i , θ 0 ) +m 2 (W i , θ 0 )|W i = E ṁ⊗2 (W i , θ 0 ) m(W i , θ 0 )(1 -m(W i , θ 0 )) π(W i ) = Θ(θ 0 ).
We consider now the covariance structure of U i . We have

cov X i δi (θ 0 )(Y * i -e β 0 X i ) + (1 -δi (θ 0 ))h β 0 (Y * i , X i ) , Σ 2 (β 0 , θ 0 )Θ -1 (θ 0 ) m i (θ 0 )ξ i (δ i -m(W i , θ 0 )) = E X i m i (θ 0 )E δi (θ 0 )(Y * i -e β 0 X i -h β 0 (Y * i , X i )) + h β 0 (Y * i , X i ) ξ i (δ i -m(W i , θ 0 ))|W i ×Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ), and
E δi (θ 0 )(Y * i -e β 0 X i -h β 0 (Y * i , X i )) + h β 0 (Y * i , X i ) ξ i (δ i -m(W i , θ 0 ))|W i = E ξ i δ i (1 -m(W i , θ 0 ))(Y * i -e β 0 X i -h β 0 (Y * i , X i )) + ξ i h β 0 (Y * i , X i )(δ i -m(W i , θ 0 ))|W i = (1 -m(W i , θ 0 ))(Y * i -e β 0 X i -h β 0 (Y * i , X i ))m(W i , θ 0 )π(W i ), therefore, cov X i δi (θ 0 )(Y * i -e β 0 X i ) + (1 -δi (θ 0 ))h β 0 (Y * i , X i ) , Σ 2 (β 0 , θ 0 )Θ -1 (θ 0 ) m i (θ 0 )ξ i (δ i -m(W i , θ 0 )) = E X i m i (θ 0 )(1 -m(W i , θ 0 ))(Y * i -e β 0 X i -h β 0 (Y * i , X i ))m(W i , θ 0 )π(W i ) Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ) = E X i ṁi (W i , θ 0 )(Y * i -e β 0 X i -h β 0 (Y * i , X i ))π(W i ) Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ) = (Σ 3 (β 0 , θ 0 ) -Σ 2 (β 0 , θ 0 )) Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ). It follows that var(U i ) = Σ 1 (β 0 ) + (2Σ 3 (β 0 , θ 0 ) -Σ 2 (β 0 , θ 0 )) Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ).
Finally, Theorem 3.1 follows from the multivariate central limit theorem and Slutsky's theorem.

Appendix B: Proof of Theorem 4.1

Consistency can be proved in much the same way as βn ; the proof is therefore omitted. We turn to asymptotic normality. A technical lemma is needed. For j = 1, . . . , M , let

˙ * n,j (β, θ) = ∂ * n,j (β, θ) ∂β = n i=1 X i δ * i,j (θ) Y * i -e β X i -h β (Y * i , X i ) + h β (Y * i , X i ) (7.8) := n i=1 f β,θ,j (O i ).
Then the following holds:

Lemma 1. Under conditions C1, C2 and C4 :

1 √ n ˙ * n,j (β 0 , θn ) -nE[ ˙ * 1,j (β 0 , θn )] -˙ * n,j (β 0 , θ 0 ) -nE[ ˙ * 1,j (β 0 , θ 0 )] P -→ 0 (7.9)
as n → ∞.

Proof of Lemma 1. In this proof, for notational simplicity, we will write f θ instead of f β 0 ,θ,j . First, note that

1 √ n ˙ * n,j (β 0 , θ) -nE[ ˙ * 1,j (β 0 , θ)] = 1 √ n n i=1 f θ (O i ) -nE[f θ (O 1 )] = G n f θ ,
where G n f θ denotes the empirical process evaluated at f θ . To prove the lemma, we first prove that the class of functions {f θ : θ ∈ Θ} is Donsker (see, for example, [START_REF] Van Der Vaart | Asymptotic Statistics[END_REF] for a detailed account on empirical processes and Donsker classes). For that purpose, we decompose f θ in (7.8) as

f θ (O i ) = X i (f 1,θ (O i )+f 2,θ (O i )+f 3,θ (O i )), where f 1,θ (O i ) = -δ * i,j (θ)e β 0 X i +h β 0 (Y * i , X i ), f 2,θ (O i ) = δ * i,j (θ)Y * i and f 3,θ (O i ) = -δ * i,j (θ)h β 0 (Y * i , X i
) and we show that the classes F 1 := {f 1,θ : θ ∈ Θ}, F 2 := {f 2,θ : θ ∈ Θ} and F 3 := {f 3,θ : θ ∈ Θ} are Donsker.

For illustration purpose, we show that F 1 is Donsker. Here, it is useful to see

D i,j (θ) ∼ B(m(W i , θ)) as the random variable 1 {U i ≤m(W i ,θ)} , where U i is a uniform random variable on [0, 1], independent of O i . Let d := diam(Θ) denote the diameter of Θ ⊂ R q .
Then the size of Θ in every direction is at most d and thus, we can cover Θ with fewer than (d/κ) q cubes of length κ. The circumscribed balls have radius a multiple κ * := ακ of κ (α > 0) and these balls also cover Θ. Now, for a given θ ∈ Θ, consider the set

{f 1, θ : θ ∈ Θ ∩ B(θ, κ * )}, where B(θ, κ * ) = { θ ∈ R q : θ -θ ≤ κ * } is the ball of radius κ * and center θ. If θ ∈ B(θ, κ * ), condition C4 implies that |m(w, θ) -m(w, θ)| ≤ h(w)κ * , hence m(w, θ)-h(w)κ * ≤ m(w, θ) ≤ m(w, θ)+h(w)κ * and thus 1 {U i ≤m(W i ,θ)-h(W i )κ * } ≤ 1 {U i ≤m(w, θ)} ≤ 1 {U i ≤m(w,θ)+h(w)κ * } . From this, we can see that f L θ (O i ) ≤ f 1, θ(O i ) ≤ f U θ (O i ),
where

f L θ (O i ) = h β 0 (Y * i , X i ) -(ξ i δ i + (1 -ξ i )1 {U i ≤m(W i ,θ)+h(W i )κ * } )e β 0 X i , f U θ (O i ) = h β 0 (Y * i , X i ) -(ξ i δ i + (1 -ξ i )1 {U i ≤m(W i ,θ)-h(W i )κ * } )e β 0 X i .
Moreover, under conditions C1, C2 and C4, there exists a finite positive constant c 1 such that

E f U θ (O i ) -f L θ (O i ) 2 ≤ 2c 1 κ * v. Therefore, [f L θ , f U θ ] is an ε-bracket for {f 1, θ : θ ∈ Θ ∩ B(θ, κ * )}, with ε 2 = 2c 1 κ * v.
Since we can cover Θ with fewer than (d/κ) q balls of radius κ * , we can cover

F 1 = {f 1, θ : θ ∈ Θ} with fewer than (d/κ) q ε-brackets [f L θ , f U θ ], with ε = √ 2c 1 κ * v.
The number of such ε-brackets is thus bounded by (αd/κ * ) q = (2αc 1 dv/ε 2 ) q , which is order ε -2q . Hence, the bracketing integral is of order 1 0

√

-2q log ε dε, which is finite. Therefore, the class of functions F 1 is Donsker, by Theorem 19.5 of [START_REF] Van Der Vaart | Asymptotic Statistics[END_REF].

By using similar arguments, we can prove that F 2 and F 3 are also Donsker classes. It follows that the class of functions {f 1,θ + f 2,θ + f 3,θ : θ ∈ Θ} is Donsker (sums of Donsker classes are Donsker). Finally, X is bounded (by condition C1), thus the class of functions {f θ : θ ∈ Θ} is Donsker.

It follows that the sequence of processes {G n f θ : θ ∈ Θ} converges in distribution to a tight limit process, and as such, is stochastically equicontinuous. Thus, Lemma 14.3 of [START_REF] Tsiatis | Semiparametric Theory and Missing Data[END_REF] and the consistency of θn imply that G n f θn -G n f θ 0 P -→ 0, which is exactly (7.9). This concludes the proof.

We come back to the proof of asymptotic normality. By a Taylor expansion of ˙ * n,j ( β * n,j , θn ) around β 0 (for j = 1, . . . , M ), we have:

0 = 1 √ n ˙ * n,j ( β * n,j , θn ) = 1 √ n ˙ * n,j (β 0 , θn ) + 1 n ∂ ˙ * n,j (β 0 , θn ) ∂β √ n( β * n,j -β 0 ) + o P (1).
Then, using Lemma 1, we obtain: (7.10) where the second line follows from a Taylor expansion of E[ ˙ * 1,j (β 0 , θn )] around θ 0 . Two technical lemmas are now needed : Lemma 2. For j = 1, . . . , M , we have

0 = 1 √ n ˙ * n,j (β 0 , θ 0 ) - √ nE[ ˙ * 1,j (β 0 , θ 0 )] + √ nE[ ˙ * 1,j (β 0 , θn )] + 1 n ∂ ˙ * n,j (β 0 , θn ) ∂β √ n( β * n,j -β 0 ) + o P (1) = 1 √ n ˙ * n,j (β 0 , θ 0 ) + √ n ∂E[ ˙ * 1,j (β 0 , θ 0 )] ∂θ ( θn -θ 0 ) + o P ( θn -θ 0 ) + 1 n ∂ ˙ * n,j (β 0 , θn ) ∂β √ n( β * n,j -β 0 ) +o P (1),
∂E[ ˙ * 1,j (β, θ)] ∂θ = Σ 2 (β, θ).
Proof of Lemma 2. First, we note that

E[δ * 1,j (θ)|W 1 ] = E[ξ 1 δ 1 + (1 -ξ 1 )D 1,j (θ)] = π(W 1 )m(W 1 , θ 0 ) + (1 -π(W 1 ))m(W 1 , θ). (7.11)
Hence, using (7.8) and iterating the expectation with conditioning on W 1 , we obtain:

E[ ˙ * 1,j (β, θ)] = E X 1 δ * 1,j (θ) Y * 1 -e β X 1 -h β (Y * 1 , X 1 ) + h β (Y * 1 , X 1 ) = E X 1 (π(W 1 )m(W 1 , θ 0 ) + (1 -π(W 1 ))m(W 1 , θ)) Y * 1 -e β X 1 -h β (Y * 1 , X 1 ) +h β (Y * 1 , X 1 ))] .
Finally, straightforward calculations yield

∂E[ ˙ * 1,j (β, θ)] ∂θ = E X 1 (1 -π(W 1 )) ṁ (W 1 , θ) Y * 1 -e β X 1 -h β (Y * 1 , X 1 ) = Σ 2 (β, θ).
Lemma 3. For j = 1, . . . , M ,

1 n ∂ ˙ * n,j (β 0 , θn ) ∂β P -→ -Σ 1 (β 0 ).
Proof of Lemma 3. Let j = 1, . . . , M . Straightforward calculations yield:

∂ ˙ * n,j (β 0 , θn ) ∂β = n i=1 X i X i -δ * i,j ( θn )e β 0 X i + (1 -δ * i,j ( θn ))h β 0 (Y * i , X i ) Y * i -e β 0 X i -h β 0 (Y * i , X i ) .
Then we decompose n -1 ∂ ˙ * n,j (β 0 , θn )/∂β as:

1 n ∂ ˙ * n,j (β 0 , θn ) ∂β = 1 n ∂ ˙ * n,j (β 0 , θn ) ∂β - 1 n ∂ ˙ * n,j (β 0 , θ 0 ) ∂β + 1 n ∂ ˙ * n,j (β 0 , θ 0 ) ∂β = 1 n n i=1 X i X i -e β 0 X i (1 -ξ i )(D i,j ( θn ) -D i,j (θ 0 )) +h β 0 (Y * i , X i ) Y * i -e β 0 X i -h β 0 (Y * i , X i ) (1 -ξ i )(D i,j (θ 0 ) -D i,j ( θn )) + 1 n ∂ ˙ * n,j (β 0 , θ 0 ) ∂β = 1 n n i=1 X i X i (1 -ξ i ) e β 0 X i + h β 0 (Y * i , X i ) Y * i -e β 0 X i -h β 0 (Y * i , X i ) ×(D i,j (θ 0 ) -D i,j ( θn )) + 1 n ∂ ˙ * n,j (β 0 , θ 0 ) ∂β = 1 n n i=1 Z i (1 {U i,j ≤m(W i ,θ 0 )} -1 {U i,j ≤m(W i , θn)} ) + 1 n ∂ ˙ * n,j (β 0 , θ 0 ) ∂β , (7.12)
where

Z i := X i X i (1 -ξ i )[e β 0 X i + h β 0 (Y * i , X i )(Y * i -e β 0 X i -h β 0 (Y * i , X i ))]
(in what follows, we will denote by Z i,( ,k) the ( , k)-th element of Z i ) and U i,j is a uniform random variable on [0, 1], independent of all other random variables.

Consider the first term in the right-hand side of (7.12). The random variable |1 {U i,j ≤m(W i ,θ 0 )} -1 {U i,j ≤m(W i , θn)} | is equal to 0 or 1 and takes the value 1 with probability |m(W i , θ 0 ) -m(W i , θn )|. Let ε > 0. Then, for , k ∈ {1, . . . , p}, Markov's inequality implies that

P 1 n n i=1 Z i,( ,k) (1 {U i,j ≤m(W i ,θ 0 )} -1 {U i,j ≤m(W i , θn)} ) > ε ≤ 1 ε E 1 n n i=1 Z i,( ,k) (1 {U i,j ≤m(W i ,θ 0 )} -1 {U i,j ≤m(W i , θn)} ) .
Under conditions C1 and C2, there exists a finite positive constant c 2 such that |Z i,( ,k) | ≤ c 2 . Thus,

P 1 n n i=1 Z i,( ,k) (1 {U i,j ≤m(W i ,θ 0 )} -1 {U i,j ≤m(W i , θn)} ) > ε ≤ c 2 εn n i=1 |m(W i , θ 0 ) -m(W i , θn )| ≤ c 2 εn n i=1 h(W i ) θ 0 -θn ≤ c 2 ε θ 0 -θn (v + o P (1)),
where the last two lines follow from the condition C4. Finally, consistency of θn implies that

1 n n i=1 Z i,( ,k) (1 {U i,j ≤m(W i ,θ 0 )} -1 {U i,j ≤m(W i , θn)}
) converges in probability to 0, and the first term in the right-hand side of (7.12) also converges to 0.

We consider now the second term in the right-hand side of (7.12). By the weak law of large numbers, n -1 ∂ ˙ * n,j (β 0 , θ 0 )/∂β converges in probability to

E X 1 X 1 -δ * 1,j (θ 0 )e β 0 X 1 + (1 -δ * 1,j (θ 0 ))h β 0 (Y * 1 , X 1 ) Y * 1 -e β 0 X 1 -h β 0 (Y * 1 , X 1 ) . (7.13) Using the fact that E[δ * 1,j (θ 0 )|W 1 ] = m(W 1 , θ 0 ) (see (7. 11
)), and iterating the expectation in (7.13) with conditioning on W 1 , we easily show that (7.13) is equal to -Σ 1 (β 0 ).

Thus, we have shown that n -1 ∂ ˙ * n,j (β 0 , θn )/∂β converges in probability to -Σ 1 (β 0 ), which concludes the proof.

By combining (7.10) with Lemmas 2 and 3, we obtain the following approximation of β * n,j :

√ n( β * n,j -β 0 ) = Σ -1 1 (β 0 ) 1 √ n n i=1 X i δ * i,j (θ 0 )[Y * i -e β 0 X i -h β 0 (Y * i , X i )] + h β 0 (Y * i , X i ) + Σ 2 (β 0 , θ 0 )Θ -1 (θ 0 ) m i (θ 0 )ξ i (δ i -m(W i , θ 0 )) + o P (1),
which in turn implies the approximation of the multiple imputation estimator β * n :

√ n( β * n -β 0 ) = 1 M M j=1 √ n( β * n,j -β 0 ) = Σ -1 1 (β 0 ) 1 √ n n i=1   1 M M j=1 X i δ * i,j (θ 0 )[Y * i -e β 0 X i -h β 0 (Y * i , X i )] + h β 0 (Y * i , X i ) + Σ 2 (β 0 , θ 0 )Θ -1 (θ 0 ) m i (θ 0 )ξ i (δ i -m(W i , θ 0 )) + o P (1) := Σ -1 1 (β 0 ) 1 √ n n i=1   1 M M j=1 f β 0 ,θ 0 ,j (O i ) + V i   + o P (1), (7.14)
where

V i := Σ 2 (β 0 , θ 0 )Θ -1 (θ 0 ) m i (θ 0 )ξ i (δ i -m(W i , θ 0 )).
We have already shown (see proof of Theorem 3.1) that var(

V i ) = Σ 2 (β 0 , θ 0 )Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ).
Similar calculations as in the proof of Theorem 3.1 yield:

cov(f β 0 ,θ 0 ,j (O i ), V i ) = (Σ 3 (β 0 , θ 0 ) -Σ 2 (β 0 , θ 0 )) Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ). Therefore, var   1 M M j=1 f β 0 ,θ 0 ,j (O i ) + V i   = var   1 M M j=1 f β 0 ,θ 0 ,j (O i )   + var(V i ) + 2 M M j=1 cov(f β 0 ,θ 0 ,j (O i ), V i ) = Σ * 1 (β 0 , θ 0 ) + Σ 2 (β 0 , θ 0 )Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ) +2 (Σ 3 (β 0 , θ 0 ) -Σ 2 (β 0 , θ 0 )) Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ) = Σ * 1 (β 0 , θ 0 ) + (2Σ 3 (β 0 , θ 0 ) -Σ 2 (β 0 , θ 0 )) Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ). (7.15)
as n → ∞ (calculations are omitted). Moreover, if the model m(W i , θ) is correctly specified (that is, θ * = θ 0 ), then Σ 6 (β 0 , θ 0 , γ * ) = 0. Similarly, if model π(W i , γ) is correctly specified (and thus, γ * = γ 0 ), then Σ 5 (β 0 , θ 0 , γ * ) = 0. Now, taking Taylor's expansion of ˘ n ( βn , θn , γn ) around (β 0 , θ * , γ * ) gives

√ n( βn -β 0 ) = - 1 n ∂ ˘ n (β 0 , θ * , γ * ) ∂β -1 1 √ n ˘ n (β 0 , θ * , γ * ) + 1 n ∂ ˘ n (β 0 , θ * , γ * ) ∂θ √ n( θn -θ * ) + 1 n ∂ ˘ n (β 0 , θ * , γ * ) ∂γ √ n(γ n -γ * ) + o P (1). (7.17)
Finally, combining (3.4), (5.6), (7.16) and (7.17) and using the limit central theorem yield the asymptotic distribution of √ n( βn -β 0 ) when either m(W i , θ) or π(W i , γ) is correctly specified. Formulas for the asymptotic variance follow from easy albeit tedious calculations.

If both m(W i , θ) and π(W i , γ) are correctly specified, Σ 7 (β 0 , θ 0 , γ 0 ) = Σ 8 (β 0 , θ 0 , γ 0 ) = Σ 1 (β 0 ) and the asymptotic variance of βn reduces to Σ -1 1 (β 0 ), which concludes the proof. -0.0022 -0.0006 0.0008 0.0000 0.0005 0.0311 0.0021 -0.0065 -0.0048 -0.0073 -0.0022 -0.0006 0.0008 0.0000 0.0005 RC SE 0.0780 0.0153 0.0329 0.0117 0.0192 0.0769 0.0153 0.0330 0.0115 0.0188 0.0780 0.0153 0.0329 0.0117 0.0192 RMSE 0.1114 0.0218 0.0463 0.0163 0.0273 0.1177 0.0220 0.0476 0.0174 0.0286 0.1114 0.0218 0.0463 0.0163 0.0273 CP 0.9490 0.9490 0.9520 0.9570 0.9430 0.9100 0.9450 0.9330 0.9160 0.9120 0.9490 0.9490 0.9520 0.9570 0.9430 bias -0.0078 -0.0009 0.0020 0.0007 0.0017 -0.0069 -0.0008 0.0018 0.0009 0.0015 -0.0060 -0.0008 0.0017 0.0006 0.0013 AIPW SE 0.0765 0.0150 0.0322 0.0114 0.0187 0.0747 0.0149 0.0319 0.0112 0.0183 0.0766 0.0150 0.0323 0.0115 0.0188 RMSE 0.1106 0.0217 0.0459 0.0161 0.0271 0.1151 0.0221 0.0474 0.0170 0.0280 0.1106 0.0217 0.0459 0.0162 0.0271 CP 0.9410 0.9450 0.9430 0.9500 0.9390 0.8990 0.9280 0.9290 0.9150 0.9070 0.9420 0.9420 0.9450 0.9520 0.9420 bias -0.0026 -0.0006 0.0009 0.0000 0.0006 0.0301 0.0019 -0.0062 -0.0047 -0.0071 -0.0026 -0.0006 0.0009 0.0000 0.0006 MI SE 0.0772 0.0150 0.0324 0.0115 0.0189 0.0805 0.0155 0.0340 0.0120 0.0196 0.0772 0.0150 0.0324 0.0115 0.0189 RMSE 0.1109 0.0216 0.0460 0.0162 0.0272 0.1199 0.0222 0.0482 0.0177 0.0290 0.1109 0.0216 0.0460 0.0162 0.0272 CP 0.9410 0.9450 0.9430 0.9460 0.9380 0.9340 0.9500 0.9470 0.9330 0.9300 0.9410 0.9450 0.9430 0.9460 0.9380 0.7500 0.8940 0.9400 0.7800 0.8920 0.7500 0.8940 0.9400 0.7800 0.8920 0.7500 0.8940 0.9400 0.7800 0.8920 bias 0.0065 0.0004 -0.0021 -0.0006 -0.0013 0.0779 0.0010 -0.0173 -0.0113 -0.0181 0.0065 0.0004 -0.0021 -0.0006 -0.0013 RC SE 0.0793 0.0154 0.0336 0.0119 0.0196 0.0764 0.0153 0.0335 0.0115 0.0187 0.0793 0.0154 0.0336 0.0119 0.0196 RMSE 0.1134 0.0217 0.0465 0.0167 0.0281 0.1387 0.0220 0.0504 0.0204 0.0335 0.1134 0.0217 0.0465 0.0167 0.0281 CP 0.9430 0.9480 0.9580 0.9490 0.9460 0.7750 0.9420 0.9200 0.8080 0.7920 0.9430 0.9480 0.9580 0.9490 0.9460 bias -0.0052 0.0000 0.0004 0.0010 0.0012 -0.0061 -0.0005 0.0008 0.0014 0.0013 -0.0017 0.0000 -0.0002 0.0007 0.0003 AIPW SE 0.0765 0.0150 0.0322 0.0112 0.0188 0.0698 0.0149 0.0310 0.0105 0.0173 0.0765 0.0150 0.0323 0.0115 0.0188 RMSE 0.1115 0.0215 0.0455 0.0163 0.0276 0.1185 0.0225 0.0480 0.0177 0.0291 0.1113 0.0216 0.0455 0.0164 0.0275 CP 0.9390 0.9400 0.9410 0.9470 0.9370 0.8501 0.9080 0.8925 0.8273 0.8635 0.9410 0.9410 0.9430 0.9550 0.9420 bias 0.0056 0.0004 -0.0018 -0.0005 -0.0011 0.0763 0.0010 -0.0168 -0.0110 -0.0177 0.0056 0.0004 -0.0018 -0.0005 -0.0011 MI SE 0.0777 0.0150 0.0326 0.0116 0.0191 0.0822 0.0156 0.0352 0.0125 0.0200 0.0777 0.0150 0.0326 0.0116 0.0191 RMSE 0.1124 0.0215 0.0458 0.0165 0.0278 0.1412 0.0222 0.0514 0.0209 0.0340 0.1124 0.0215 0.0458 0.0165 0.0278 CP 0.9440 0.9420 0.9440 0.9470 0.9420 0.8330 0.9450 0.9390 0.8690 0.8440 0.9440 0.9420 0.9440 0.9470 0.9420

correct m(W, θ) / incorrect π(W, γ) incorrect m(W, θ) / correct π(W, γ) both models correct estimator β 1 β 2 β 3 β 4 β 5 β 1 β 2 β 3 β 4 β 5 β 1 β 2 β 3 β 4 β 5 bias -0.

  and an indicator δ (called censoring indicator thereafter) which is equal to 1 if Y < C and 0 if Y ≥ C. Finally, we denote by Y * the observed count value (that is, Y * = min(Y, C)).

  .9388 0.9519 0.8656 0.9188 0.8676 0.9388 0.9519 0.8656 0.9188 0.8676 0.9388 0.9519 0.8656 0.9188 bias

			0060 -0.0008 0.0014 0.0005 0.0013 -0.0060 -0.0008 0.0014 0.0005 0.0013 -0.0060 -0.0008 0.0014 0.0005 0.0013
	full data SE	0.0766 0.0150 0.0323 0.0115 0.0188 0.0766 0.0150 0.0323 0.0115 0.0188 0.0766 0.0150 0.0323 0.0115 0.0188
		RMSE 0.1099 0.0215 0.0453 0.0162 0.0269 0.1099 0.0215 0.0453 0.0162 0.0269 0.1099 0.0215 0.0453 0.0162 0.0269
		CP	0.9460 0.9490 0.9500 0.9560 0.9420 0.9460 0.9490 0.9500 0.9560 0.9420 0.9460 0.9490 0.9500 0.9560 0.9420
		bias	0.0688 -0.0038 -0.0051 -0.0109 -0.0092 0.0688 -0.0038 -0.0051 -0.0109 -0.0092 0.0688 -0.0038 -0.0051 -0.0109 -0.0092
	CC	SE	0.0866 0.0163 0.0350 0.0131 0.0207 0.0866 0.0163 0.0350 0.0131 0.0207 0.0866 0.0163 0.0350 0.0131 0.0207
		RMSE 0.1421 0.0238 0.0493 0.0213 0.0310 0.1421 0.0238 0.0493 0.0213 0.0310 0.1421 0.0238 0.0493 0.0213 0.0310
		CP	0.8676 0

Table 1 :

 1 Simulation results for experiment 1. SE: average standard error. RMSE : root mean square error. CP: empirical coverage probability of 95%-level confidence intervals. .0001 -0.0001 0.0006 0.0003 -0.0015 -0.0001 -0.0001 0.0006 0.0003 -0.0015 -0.0001 -0.0001 0.0006 0.0003 full data SE 0.0765 0.0150 0.0323 0.0115 0.0188 0.0765 0.0150 0.0323 0.0115 0.0188 0.0765 0.0150 0.0323 0.0115 0.0188 RMSE 0.1100 0.0213 0.0447 0.0162 0.0270 0.1100 0.0213 0.0447 0.0162 0.0270 0.1100 0.0213 0.0447 0.0162 0.0270 CP 0.9370 0.9450 0.9520 0.9540 0.9440 0.9370 0.9450 0.9520 0.9540 0.9440 0.9370 0.9450 0.9520 0.9540 0.9440 bias 0.1224 0.0121 -0.0128 -0.0177 -0.0152 0.1224 0.0121 -0.0128 -0.0177 -0.0152 0.1224 0.0121 -0.0128 -0.0177 -0.0152 CC SE 0.0993 0.0186 0.0391 0.0151 0.0233 0.0993 0.0186 0.0391 0.0151 0.0233 0.0993 0.0186 0.0391 0.0151 0.0233 RMSE 0.1864 0.0289 0.0562 0.0275 0.0363 0.1864 0.0289 0.0562 0.0275 0.0363 0.1864 0.0289 0.0562 0.0275 0.0363 CP

		correct m(W, θ) / incorrect π(W, γ)	incorrect m(W, θ) / correct π(W, γ)		both models correct
	estimator	β 1	β 2	β 3	β 4	β 5	β 1	β 2	β 3	β 4	β 5	β 1	β 2	β 3	β 4	β 5
	bias	-0.0015 -0

Table 2 :

 2 Simulation results for experiment 2. SE: average standard error. RMSE : root mean square error. CP: empirical coverage probability of 95%-level confidence intervals.

Acknowledgements

Authors acknowledge financial support from the Hubert Curien "PHC-Utique" program (CMCU number : 20G1503 -Campus France number : 44172SL), implemented by Campus France.

Finally, it follows from (7.14), (7.15) and the multivariate central limit theorem that √ n( β * n -β 0 ) converges in distribution to a Gaussian vector with mean zero and variance Σ -1 1 (β 0 ) Σ * 1 (β 0 , θ 0 ) + (2Σ 3 (β 0 , θ 0 ) -Σ 2 (β 0 , θ 0 )) Θ -1 (θ 0 )Σ 2 (β 0 , θ 0 ) Σ -1 1 (β 0 ), which concludes the proof.

Appendix C: Proof of Theorem 5.1

Assume that the model m(W i , θ) is correctly specified. It is straightforward to check that the map β → ∂ ˘ n (β, θn , γn )/∂β exists and is continuous in a neighborhood of β 0 (condition i).

Now, we show that n -1 ˘ n (β 0 , θn , γn ) = o P (1) (condition ii). To see this, decompose n -1 ˘ n (β 0 , θn , γn ) as:

First, we consider the term Q

(1)

We have:

n,1, and Q i, denote the -th component of the vectors

n,1 and Q i respectively (for = 1, . . . , p), we have:

Conditions C1 and C6 ensure that there exists a finite positive constant c 3 such that

and the condition C7 implies that

Finally, the convergence of γn to γ * implies that Q

(1)

n,1, ( = 1, . . . , p), and thus

n,1 , converge to 0 as n → ∞. Similarly, under conditions C1 and C6, there exists a finite positive constant c 4 such that

and condition C4 implies

n,3, ( = 1, . . . , p), and thus Q

(1) n,3 , converge to 0 as n → ∞. Finally, by the law of large numbers,

With obvious notations, we have Q

(2)

n (θ 0 ). By the law of large numbers, Q

(2)

which is equal to 0 (see proof of Theorem 3.1). We also have

and using similar arguments as for Q

(1)

n,3 , we can show that this converges to 0 if model m(W i , θ) is correctly specified. Finally, Q

(2) n ( θn ) = o P (1), which concludes the proof of condition ii. Now, we prove that n

, some easy calculations yield:

Using similar arguments as for Q

(1) n,3 (respectively Q

(1)

n and Q

(1) n,2 ), we can show that T

(1)

n (respectively T

(2)

It follows that n -1 ∂ ˘ n (β, θn , γn )/∂β converges in probability to -Σ 1 (β). Uniformity of the convergence follows by the same arguments as in the proof of Theorem 3.1. Finally, having proved conditions i, ii and iii, we apply the inverse function theorem of [START_REF] Foutz | On the unique consistent solution to the likelihood equations[END_REF] and conclude that βn converges in probability to β 0 if m(W i , θ) is correctly specified. The consistency proof of βn when model π(W i , γ) is correctly specified proceeds along the same lines and is omitted.


domain: phys chim spi

text:
K Chevrier

J M Benoit

C Symonds

S K Saikin

J Yuen-Zhou

J Bellessa

Anisotropy and controllable band structure in supra-wavelength polaritonic metasurfaces

Keywords: numbers: 71.36.+c, 42.25.Kb, 73.20.Mf, 78.55.Kz

In this letter we exploit the extended coherence length of mixed plasmon/exciton states to generate active metasurfaces. For this purpose, periodic stripes of organic dye are deposited on a continuous silver film. Typical metasurface effects, such as effective behavior and geometry sensitivity, are measured for periods exceeding the polaritonic wavelength by more than one order of magnitude. By adjusting the metasurface geometry, anisotropy, modified band structure and unidimensional polaritons are computationally simulated and experimentally observed in reflectometry as well as in emission.

The strong-coupling phenomenon that occurs when light-matter interaction prevails over damping [START_REF] Kimble | Structure and Dynamics in Cavity Quantum Electrodynamics[END_REF] leads to remarkable changes in the optical properties of matter, modifying both the spectrum [START_REF] Weisbuch | [END_REF] and dynamics of excitations [3]. In molecular systems [4], this regime can be reached with a reduced number of molecules, opening the way to the realization of single photon emitters and quantum information processing architectures [5]. When a large number of emitters is involved, a coherent coupling of the emitters through the electromagnetic mode occurs [6,7]. In metal/semiconductors systems, the strongly-coupled state (polariton) arising between molecules and a propagating plasmon has an extension ranging from a few [7,8] to several tens of microns [9]. The conductivity of molecular films can be dramatically modified through this coherent coupling [10,11]. The strong-coupling regime also affects other properties such as chemical reactivity [12] or efficiency of second harmonic generation [13]. Within this context, the control of polaritonic modes is of crucial interest. Polaritons can be tailored by shaping their electromagnetic part: mesa [14], pillar structuration [15] or photonic metasurfaces [16,17] have, for example, allowed for an efficient control of their spectra. Various types of plasmonic modes have also been exploited in strong coupling to reduce volume and losses [18,19]. This control can also be achieved by tailoring the excitonic part: a recent demonstration of polarization control and energy tunability has been performed with aligned carbon nanotubes coupled to a microcavity [20]. Beside the properties inherited from their light or matter components, the hybrid nature of polaritonic states leads to unique properties which could be obtained neither by the photon nor exciton alone. This has yield for instance the first observation of Bose condensation in condensed matter physics, with exci-ton/photon polaritons [21][22][23]. This hybrid nature can also be exploited to tailor the polaritons properties in the linear regime.

In the present work, we exploit the mixing of propagating plasmons with a set of localized emitters to control their band structure and generate anisotropy. More specifically, by patterning an organic film on silver, we investigate a new type of metasurface that builds on the polariton extended coherence. As the coherence length is one order of magnitude larger than the wavelength, an efficient geometry control can be achieved in the visible range with micrometric structuration. These findings pave the way towards a novel type of artificial materials which require micro-instead of nanostructuration. The developments of new functionalities based on polaritonic metasurfaces could also build on the extraordinary achievements of metamaterials [START_REF] Engheta | Metamaterials: physics and engineering explorations[END_REF][START_REF] Monticone | [END_REF][26] and metasurfaces [27,28], with potential electrical or optical control enabled by the excitonic part of the polaritons [29,30].

The sample consists of stripes of active dye, periodically deposited on a 50 nm continuous silver film thermally evaporated on a glass substrate. For this purpose a solution of J-aggregated dye (TDBC) [7] is spin coated on top of the silver film. The spin coated layer has a thickness of 17 nm, measured by plasmon fitting, with a sharp absorption peak at 2.1 eV. The inactive stripes have been created in the homogeneous cyanine dye film by photobleaching the dye molecules using local UV-irradiation performed with a 355 nm solid state laser. By running an adjustable mechanical slit at an intermediate image plane, a micrometric precision is reachable for the size of the patterned structures.The UV-irradiated material has the same optical index as the cyanine dye film, except in the 2.1 eV region where dielectric variation related to the sharp absorption is eliminated [40]. Optical properties of the samples are studied using leakage radiation microscopy [31]through a Nikon oil immersion objective of 1.49 numerical aperture [7]. Reflectometry images in the Fourier plane, which represent the dispersion relation (energy as a function of wavevector), are collected by imaging the microscope objective Fourier plane by a Nikon CCD camera or an Andor spectrometer, coupled to a CCD camera. Luminescence is triggered with a 532 nm laser diode beam focused onto the upper face (dye layer) of the sample. Direct space images, with selected region of Fourier space can be performed by mounting an iris (selecting the central part of the Fourier plane) or a circular beam blocker in an intermediate Fourier plane.

Before addressing the patterned surfaces, the properties of uniform layers of UV-irradiated and of active dye on silver were measured. Figure 1(a) shows the dispersion relation of irradiated dye on silver. The line corresponding to the surface plasmon is uniformly broadened showing that the dye absorption resonance has been completely removed by the irradiation. On the other hand, the dispersion relation of the active dye film on silver (Fig. 1(b)) shows a typical polariton anticrossing associated with the exciton/plasmon hybridization in a strong coupling regime. The Rabi splitting at resonance is 155 meV with a bare dye broadening of 30 meV and a plasmon linewidth of 40 meV, fulfilling the strong coupling conditions [32]. The polariton coherence length deduced from its wavevector width at resonance is about 7 µm [9] [40].

To demonstrate averaging of optical properties within the coherence length, characteristic of metasurfaces, we compared patterned films with two periods: one larger (25 µm) and the other smaller (5 µm) than the coherence length. In both cases the filling factor was 1/2. The dispersion relations recorded for propagation direction (wavevector) perpendicular to the stripes are shown in Figs. 1(c) and 1(d). For the larger period (Fig. 1(c)), the dispersion shows the coexistence of a continuous line and an anticrossing between two branches. The Rabi splitting associated with the anticrossing is 145 meV, which is close to the one obtained with the homogeneous dye film (Fig. 1(b)), while the continuous line has essentially the same position as the bare plasmon line (Fig. 1(a)). This behavior can be viewed as the response of two distinct regions with strong and weak plasmon-molecule coupling, and negligible interaction between them. For the smaller period (Fig. 1(d)), the measured dispersion relation is drastically different: only an anticrossing is visible, showing no regular plasmon propagation anymore. Furthermore, the Rabi splitting at resonance ( Ω= 98 meV) is reduced compared to the continuous dye film, confirming the effective behavior of the metal/dielectric metasurface, in that it differs from the basis constituents. It has to be noticed that another anticrossing, translated by k =1.31µm -1 , can be seen in Fig. 1(d). This additional line corresponds to the first diffraction order of the periodic structure (for the 5µm period, 2π/period=1.26µm -1 ) and is also present for larger periods.

The first qualitative interpretation for these observations is that, for small periods, the system averages properties of the materials with and without the dye resonance absorption. Since the average number of emitters interacting with the plasmon mode is now reduced by a factor of 2 compared to the continuous dye layer, the effective interaction energy should be reduced by a factor √ 2 [33]. Patterns of different periods were studied and the extracted Rabi splittings at resonance are summarized in Fig. 1(e). The splittings have been normalized with respect to the Rabi splitting of the continuous dye film. For small periods up to 15 µm, a quasi-plateau can be measured, while for large periods the splitting reaches the continuous film value. The striking point lies in the large value of the period at which the transition between the mean-field (effective) and multicomponent (polariton and plasmon) behavior occurs, more than one order of magnitude greater than the wavelength.

The properties of the system experimentally studied above can be theoretically rationalized by invoking a quantum mechanical framework for strong plasmonexciton coupling. The resulting Hamiltonian [40] is given by H = K H K ( = 1), where

H K = κ ω K+κ a † K+κ a K+κ + ω e ′ κ ′ σ † K+κ ′ σ K+κ ′ + κ ′ κ ′ J K+κ K κ-κ ′ σ † K+κ ′ a K+κ + h.c. ...(1)
Here, the two-dimensional wavevector K belongs to the Brillouin zone 0, 2π L × -π a , π a defined by the unit cell of the simulation, taken to be L × a, where a ≪ L is the smallest spatial resolution of the simulation and L is the patterning period. a † k (a k ) creates (annihilates) a plasmon with energy ω k , at wavevector k. σ † k (σ k ) is the active dye Frenkel exciton in Fourier space, which creates (annihilates) a delocalized collective exciton with energy ω e at wavevector k [40]. J k is a collective plasmonexciton coupling for the continuous dye film (with no patterning). Notice that the dye patterning enables mixing between plasmons and excitons (created by the operator σ † k ) as long as their wavevectors differ by a Bragg vector κ = κ x x, where κ x = 2πm L (with m integer). The efficiency with which this mixing occurs is dictated by the Fourier transform

F ∆κ = 1 ζ Nuc-1 ν=0
p ν e i∆κνa of the patterning function p ν , where νa = 0, a, • • • , (N uc -1)a are the x-positions of the localized molecular excitations (the stripes are along the y axis, the patterning along x ) and N uc a = L. Given a filling factor f , p ν = 1 for ν = 0, • • • , f N uc -1 and p ν = 0 otherwise, while the normalization is given by ζ = √ f N uc . Finally, ′ κ ′ indicates a constrained summation over Bragg vectors associated with an integer number of wavelengths per active stripe. It is clear that in the absence of patterning, F κ-κ ′ = δ κ,κ ′ and Eq. (1) reduces to the standard polaritonic dispersion Hamiltonian.

Plasmon-weighed dispersion plots are obtained by numerical diagonalization of H K [40]. Results for the continuous inactive and active dye films are plotted in Figs. 1(f) and 1(g), respectively. For a patterned film of small period (5 µm, Fig. 1(i)), a very clear anticrossing with 1/ √ 2 Rabi splitting of the continuous film can be observed, verifying our effective medium picture. For a large period of 25 µm (Fig. 1(h)), the stripes support a large number of modes, upon which the calculation shows a clear partitioning between those which are confined in the active and inactive dye stripes, respectively. The former feature strong coupling with the dyes, while the latter remains as plasmons localized in the inactive dye stripes. As no losses are included in the calculation, the calculated lines have no width. The losses, which limit the coherence length, broaden the experimental lines and explain the merging of the lines which is observed in the experiments: the plasmonic discrete modes appearing in the anticrossing gap form the propagating plasmon band, and the different diffracted polaritonic lines tend to merge to form a new anticrossing.

Another feature of the metasurface geometry studied here is its anisotropy, which has been investigated for a small pattern period (5 µm) by recording the dispersion relations for different propagation directions. For propagation along the stripes, the dispersion presented in Fig. 2(a) can be understood as the coexistence of bare plasmon and polariton lines, in contrast to that recorded perpendicular to the stripes (Fig. 2(b)) which only shows polariton lines. The value of the Rabi splitting parallel to the stripes is 160 meV, almost equal to that of the continuous dye film. A first interpretation based on this energy is that the effective medium effect is not present for propagation parallel to the stripe: polaritons confined along the active dye stripes coexist with propagating plasmons confined along the bare metal regions. Our numerical calculations unambiguously reproduce the anisotropic dispersion of the system as shown in Figs. 2(c) and 2(d). Figure 2(e) shows the full experimental two-dimensional dispersion for a 5 µm period. This surface is deduced by tracking reflectivity minima as a function of azimuthal angle. A propagating plasmon is present only for angles between ±20 • about the direction of the stripes. For larger angles, the bare plasmon line is not detected anymore, and only the polariton lines remain.

The described metasurfaces are also accompanied by substantial modifications in the emission properties of the constitutive dye. The hybridization of a large number of molecular excited states with an extended plasmonic mode results in two types of states: extended-polaritonic states (with a well-defined wavevector) and localized incoherent-dispersionless purely molecular states also known as dark or reservoir states [6,34]. The polaritonic states and the plasmon modes can be respectively selected by choosing wavevector ranges under the light cone (non-radiative modes). Figure 3(a) presents a Fourier image of the emission of a polaritonic metasurface, with the horizontal axis x perpendicular to the stripes, nonresonantly excited by a continuous 532 nm laser. Figure 3(b) shows the direct luminescence imaging of the incoherent states, where the emission coupled to the plasmon is removed by selecting the central part of the Fourier plane, i.e. the small wavevectors. In this case, the dye stripes appear clearly, and no emission is observed between them in the inactive dye areas.

When the plasmon/polariton emission is selected (large wavevectors), a strong polarization dependence can be observed. Selecting the polarization in real space is equivalent to choosing the propagation direction in re- ciprocal (Fourier) space, as the plasmon and thus the polaritons are TM polarized. Strikingly, for a horizontal polarization (i.e. propagation direction k x , Fig. 3(c)) the stripes are hardly visible. This uniform luminescence confirms the interpretation of an effective medium behavior for periods smaller than the coherence length. For a vertical polarization (i.e. propagation direction k y , Fig. 3(d)) the stripes appear clearly although with a lower contrast than the incoherent luminescence. To analyze the luminescence images, the low energy polariton wavefunction at the resonance has been calculated. Figures 3(e) and 3(f) shows the spatial partition of the plasmonic and excitonic densities associated with this wavefunction. An essentially uniform plasmonic density is obtained for a propagation perpendicular to the stripes whereas large maxima within the stripes can be seen for a propagation parallel to the stripes. As the intensity collected in the experiments is associated with the plasmonic part [3] of the wavefunction, the calculated densities and measured luminescence are fully consistent.

The hybridization with the plasmon renders luminescence dispersive and modifies emission energies. The luminescence as a function of wavevector for the 5 µm period array (dispersion of Fig. 2) is shown in Figs. 4(a) and 4(b). The polariton appears for both propagation directions, while an emission through the plasmon can be detected only along k y . By changing the filling factor, the effective concentration of active dye molecules within the plasmon coherence length can be modified, allowing a transition from strong to weak coupling. With a filling factor of 1/10, an emission without any polaritonic contribution can be observed on Fig. 4(c) for a detection perpendicular to the stripes. This emission corresponds to the typical weak coupling regime for molecules close to a metal surface: an emission at the dye transition energy into the plasmon. For propagation along k y (Fig. 4(d)), both emission into plasmon and a polariton can be detected. These observations agree with the previous interpretation of two coexisting modes, spatially localized in the regions of the active (polariton) and inactive (plasmon) dyes. We have thus demonstrated unidirectional polaritons and anisotropic luminescence starting from an isotropic emitting material. It has to be noted that the hybrid state emission can be observed at low energies (1.95 eV) [40] where the emission of the bare dye vanishes. This shows that the observed emission cannot be interpreted in term of filtering, but that the metasurface really alters the energy of the molecules.

In conclusion, we have shown that the properties of molecules strongly coupled to surface plasmons can be drastically modified by structuring the sample on distances which are large compared to wavelengths. The modified material band structure can find applications in polarization controlled sources, as well as in topological photonics based on anisotropic dispersions and Dirac cones [35] or exceptional points [20]. The anisotropic behavior demonstrated here could also open the way to anisotropic conductivity based on strong coupling [10]. The active polaritonic metasurfaces we propose are not restricted to plasmonic modes nor organic dyes alone, but can be extended to a large variety of structures in strong coupling such as cavity polaritons which have similar extended coherence lengths [36]. Finally, building on the excitonic part of the hybrid states, these ideas can also be harnessed to design novel ultrafast reconfigurable metasurfaces which can be optically [37][38][39] or electrically [29] switched upon pumping.

FIG. 1 :

 1 FIG. 1: (a-d) Experimental reflectometry dispersion images of a 50 nm silver film covered by a 17 nm thick dye layer, measured along the propagation direction perpendicular to the stripes: (a) homogeneously UV-irradiated film; (b) homogeneous active dye film; (c) patterned layer with a large period (25 µm); (d) patterned layer with small period (5 µm). (e) Rabi splitting recorded as a function of the array period and normalized to the Rabi splitting of a homogeneous active dye layer. The dashed lines correspond to a Rabi splitting decreased by a factor of √ 2 compared to a full active dye film. (f-i) Computed dispersion curves for continuous and patterned films: (f) bare plasmon dispersion; (g) plasmon coupled to a continuous active dye film; (h) patterns of 25 µm period; (i) patterns of 5 µm period. All the patterned samples have a filling factor of 1/2.

FIG. 2 :

 2 FIG.2: Dispersion relations measured for a 5 µm period pattern with propagation direction (a) parallel to the stripes and (b) perpendicular to the stripes, as well as corresponding calculated dispersion relations (c),(d). (e) Two-dimensional experimental dispersion surface reconstructed across all propagation directions (from directions orthonormal to the stripes kx to directions parallel to the stripes ky). For clarity, the region of wavevectors with magnitudes lower than 10 µm -1 is removed.

FIG. 4 :

 4 FIG. 4: Anisotropic emission dispersion properties of polaritonic metasurfaces with two different filling factors. (a-b) Filling factor 1/2 (stripe width is 2.5 µm, period is 5 µm). (c-d) Filling factor 1/10 (stripe width 1 µm).

  C.S., J.M.B, and J.B. acknowledge financial support from ANR on AAPG 2018 project PlasHybrid, S.K.S. thanks the Center for Excitonics, an Energy Frontier Research Center funded by the U.S. Department of Energy, under Award No. DE-SC0001088 and also the Ministry of Education and Science of the Russian Federation, Award No. 3.2166.2017/4.6. J.Y.Z. acknowledges support from NSF CAREER Award CHE-164732.