J Braga
email: braga@cict.fr

Y Heuze
email: heuze@cict.fr

Quantifying variation in human dental development sequences: An EVO-DEVO perspective

Keywords: dental development, dental evolution, EVO-DEVO, modularity, human evolution

The present paper describes a novel analytical approach to provide a comprehensive description of the complex interactions that exist between the growing permanent mandibular teeth (excluding the third molars), and to quantify variability in sequences of key events during crown and root formation, independent of chronological age. Importantly, our method integrates the fundamental concept of modularity and rejects the old statistical fallacy of analyzing data on the assumption that it contains no information beyond that revealed on a tooth-bytooth analysis. Indeed, interactions between growing teeth may also contain some information, which enables developmental or evolutionary information to be uncovered. Our training sample is based upon cross-sectional standardized panoramic radiographs of the teeth of a total of 2089 children (1206 girls and 883 boys) of different geographic origins (mainly Western Europe, Southern Iran, and Ivory Coast). We observe that, in extant humans sampled so far, the relative development of the permanent incisors is more plastic and varies more than for other teeth. Therefore, we consider that the quantification of possible variations between onsets, durations and rates of development of different teeth in any given child, within a large sample, is a prerequisite to the analysis of fossil hominids. In particular, we seriously question the assumption that the anterior teeth can serve as a reliable substitute for the other permanent teeth, and in particular for interpretations on somatic maturity and brain size. Our hypothesis of modularity in dental development and our method derived from this concept can serve as a basis for identifying and studying patterns of dental growth and, importantly, for comparisons between extant populations, and/or fossil species. These studies do not need to be hedged with age assessments of unknown accuracy and reliability levels (particularly in fossils), or the assumption of independence between growing teeth.

Introduction

The aim of our cross-sectional analysis of subadult radiographs is to examine, using an appropriate method and samples from different extant human populations, the variation in the relative sequences (or patterns) of key events during crown and root formation. To date, investigators have focused on the timescale (chronology) of dental development with published median or mean ages of emergence for permanent teeth in different modern human populations (e.g., [START_REF] Kuykendall | Dental development in chimpanzees (Pan troglodytes) and implications for dental development patterns in fossil hominids[END_REF][START_REF] Liversidge | Variation in modern human dental development[END_REF] (in some studies, emergence is often confused with eruption; see Marks and Cahill, 1993 for a definition of the latter term). The debates have focused on the existence of real differences between samples, or on biases due to sampling, rating effects, and statistical procedures [START_REF] Smith | Standards of human tooth formation and dental age assessment[END_REF]. However, sequence and timing represent two separate aspects of variation. For example, the first permanent molars of two individuals can mineralize at the same time but in different sequences relative to the other teeth. Compared to variation in timing, variation in patterns is much less documented in the literature. The question of variability in patterns is fundamental to the study of dental developmental processes as a possible foundation for morphological changes during human evolution. For permanent teeth, data are available for variation in sequences of emergence (e.g., [START_REF] Garn | Negro-Caucasoid differences in permanent tooth emergence at a constant income level[END_REF][START_REF] Smith | Polymorphisms in eruption sequence of permanent teeth of American children[END_REF][START_REF] Nonaka | Changes in the eruption order of the first permanent tooth and their relation to season of birth in Japan[END_REF]) and for differences in sequences of development revealed only between pairs of teeth (e.g., FanningandMoorrees,1969;Tompkins,1996;[START_REF] Liversidge | Growth of permanent mandibular teeth of British children aged 4 to 9 years[END_REF]. Variability in sequences of dental development in extant humans is acknowledged by some scholars [START_REF] Mann | Patterns of ontogeny in human evolution: evidence from dental development[END_REF][START_REF] Smith | Growth and development and its significance for early hominid behaviour[END_REF] but the problem has not properly addressed so far, due to lack of definition, concepts and methods (see below). Little is known of possible shifts in the sequences of formation of teeth that might be due to sex, age, jaw, and geographic origin, among other factors. It is generally assumed that "the pattern of development should be more robust than is age of stage appearance" (Smith, 1989: 77) and that there is little variability in relative sequences of key events during tooth growth. The comparative analyses of the sequences of dental development observed in fossil hominids have led to opposite views regarding their "modern human-like" or "apelike" status [START_REF] Smith | Dental development in Australopithecus and early Homo[END_REF][START_REF] Smith | Growth and development and its significance for early hominid behaviour[END_REF][START_REF] Smith | Patterns of dental development in Homo, Australopithecus, Pan, and Gorilla[END_REF][START_REF] Mann | Patterns of ontogeny in human evolution: evidence from dental development[END_REF][START_REF] Lampl | Furtherobservations on a method for estimating hominoid dentaldevelopmentalpatterns[END_REF][START_REF] Conroy | Paleopediatrics: or when did human infants really become human?[END_REF]. Discrepancies are due to the lack of studies of variability in patterns of dental development both within extant humans and between the two chimpanzee species. In addition to the problem of quantifying variability and the controversy over which standard is the most appropriate, other limitations on the study of dental developmental patterns are mainly methodological: (i) first, the understanding of patterns of developmental relationships should not be age related because non-adult dental age assessment is another complex problem with factors influencing its quality (accuracy and reliability) [START_REF] Ritz-Timme | Age estimation: the state of the art in relation to the specific demands of forensic practice[END_REF][START_REF] Braga | Non-adult dental age assessment: correspondence analysis and linear regression versus Bayesian predictions[END_REF]; (ii) second, developing teeth do not mineralize randomly with respect to one another. We do not expect partial correlations between the rates of tooth formation to approach zero when controlled for age. Indeed, teeth are topographically, developmentally and functionally associated with each other. Teeth essentially grow as a unit (see below) and should be statistically considered as dependent units which grow within a developmental module (Figure 1). Any assumption of independence (in a statistical sense) between developing teeth is a fallacy. Age should be estimated not only from as many available teeth as is possible, but also by considering teeth as dependent units [START_REF] Braga | Non-adult dental age assessment: correspondence analysis and linear regression versus Bayesian predictions[END_REF]; this can be done even in the case of fragmentary remains [START_REF] Heuzé | Chronologie et étiologie de la maturation macrostructurale des dents définitives[END_REF]; and (iii) third, usually more than two teeth are represented in sequences of dental development. The intricate relationships existing in the set of all teeth need to be examined for an accurate, reliable and comprehensive outline of dental developmental patterns.

The dental developmental status of an individual corresponds to the degree of various growth processes: the formation of crypts, the mineralization of crowns and roots, and their eruption status. Any combination of offset, duration and rate between developing teeth results in a sequence of dental development. A sequence of dental development is a series of discrete events, from the first evidence of mineralization to the closure of the root apex in all teeth (i.e., based on ratings of tooth formation made from panoramic radiographs) (Figure 1). It affects the growing teeth during pre-eruptive (movements of the growing tooth germs within the alveolar process before root formation), pre-functional eruptive (pre-emergent and post-emergent) and functional eruptive phases, when teeth move relative to each other and relative to the developing jaws. A dental mineralization sequence then corresponds to growing elementary units (teeth) embedded within the jaws, observed at a given developmental stage of an individual. Teeth are connected by many interactions that channel variation in sequences of dental development in nonrandom ways, and that lead to morphological patterns of co-variation. Moreover, even if minor, teeth also exhibit some degree of interaction with their respective jaws. We still have not clearly identified these interactions; however, experimental and clinical observations have demonstrated that tooth development and alveolar bone growth are interdependent [START_REF] Marks | Regional control by the dental follicle of alterations in alveolar bone metabolism during tooth eruption[END_REF][START_REF] Wise | Cellular, molecular, and genetic determinants of tooth eruption[END_REF]. There is no unifying analytical approach that provides a description of the dynamics of dental development. It is, therefore, important to identify within and between species both the more conservative and more plastic patterns of co-variation during dental growth. Indeed, conservative interactions among developing teeth arise primarily by strong, direct connections between their developmental pathways. In this case, developing teeth will be considered to co-vary as a unit during evolution, and represent an integrated module [START_REF] Olson | Morphological Integration[END_REF]. On the contrary, plastic patterns of interactions arise primarily in parallel but separate developmental pathways. Therefore, analysis of the complex interactions existing between growing teeth may well represent a new field in evolutionary developmental biology (EVO-DEVO). Are developmental changes in a set of teeth accompanied by correlative changes in other teeth in extant humans and chimpanzees? How did morphological patterns of covariation evolve in humans, and which changes in interactions between growing teeth may have produced evolutionary changes? To contribute to answering this question, we hypothesize that the developing dentition represents a developmental module with its several discernible properties defined by [START_REF] Raff | The Shape of Life: Genes, Development and the Evolution of Animal Form[END_REF]: (i) it has "an autonomous, genetically discrete organization"; (ii) it is composed of "hierarchical units and may in turn be parts of larger hierarchical entities"; (iii) it has a physical location within the developing system; (iv) it exhibits varying degrees of connectivity to other modules; (v) it undergoes temporal transformations. Indeed, a sequence of dental development in an individual is not necessarily stable through time. For example, in a given individual examined at two successive stages A and B, of his somatic development; teeth A1 and A2 may show, at stage A, the same degree of mineralization whereas, at stage B, the same teeth, B1 and B2, may be shifted and then show two different degrees of mineralization. As a cross-sectional assessment, this study does not examine the temporal variability of a sequence in an individual. This is also the case in studies applied to fossil hominids and those based on radiographic stages.

This paper describes a novel analytical approach: (i) to provide a comprehensive description of the complex interactions existing between the growing permanent teeth (third molars excepted); and (ii) to quantify variability in sequences of dental development among extant humans at the gross anatomical level and independently of chronological age. Importantly, our method integrates the fundamental concept of modularity by rejecting the old fallacy of analyzing data on the assumption that it contains no information beyond that revealed on a toothby-tooth analysis. On the contrary, we believe that the phylogenetic and taxonomic significance of dental developmental patterns arises from the interactions between growing teeth, for it is these that may contain the information which enables developmental or evolutionary information to be uncovered.

Materials and Methods

Materials

The training sample is based on crosssectional standardized panoramic radiographs of the teeth of 2089 children (1206 girls and 883 boys) aged 36 to 192 months. The geographic composition of the sample is presented in Table 1. These radiographs were taken when indicated by treatment need, and later digitally photographed for subsequent studies in our department (University of Bordeaux 1). Importantly, because magnification may vary and distortions are not uncommon, and also because of the absence of growth standards for tooth length in children from various geographic populations, we decided to assess mineralization radiographically by using fractions of crown and root formed (see above). We selected the children who were clinically free of anomalies in tooth numbers, size or shape. None of the children had undergone orthodontic treatment at the time of the selected radiograph. Our study was ethically approved by a French institutional board (Commission Nationale de l'Informatique et des Libertés); no regulation existed in Iran or in Ivory Coast. Our study includes the left mandibular permanent teeth except third molars, owing to the well-known practical difficulties in the radiographic visualization of the maxillary teeth. We exclude the third molar from our study for several reasons. First, this tooth is often not observed in our sample, and, second, the absence of this developing tooth may or may not correspond to its agenesis. Owing to the high variability (still to be measured across extant human populations) in the timing of the third molar development [START_REF] Mincer | study of third molar development and its use as an estimator of chronological age[END_REF], we cannot determine whether an absence of the germ corresponds to an agenesis or to a "delay" in the mineralization of this tooth.

Radiology has the disadvantage of detecting mineralization with some delay as compared to histology. However, the principal strength of radiology is that it allows a much better assessment of developmental variability: it allows large samples of individuals of varying geographic origin to be observed and any given child will be represented by each tooth category (i.e., incisors, canine, premolars, and molars). Indeed, when microstrutural studies have been conducted on relatively large samples (e.g., [START_REF] Yuan | Perikymata counts in two modern human sample populations[END_REF] (still smaller samples than in radiographic studies), all the teeth of a dentition have never been analysed.

Rating Dental Development

Each tooth was assigned a developmental status from the first evidence of crown mineralization to the closure of the root apex, by using the height-stages (from A to H) system described by [START_REF] Demirjian | A new system of dental age assessment[END_REF] with its precise descriptive points, radiographs and line drawings. We believe that this rating method, with its relatively small number of formation stages, is more accurate than the 14-stage system developed by Moorrees et al. (1963a) for the permanent teeth. Reducing the number of stages decreases the possibilities of inter-and intra-observer errors. These errors also decrease with increased training in scoring. Some scholars argue that rating dental development by using fractions of crown and root formed is largely subjective, because we do not know the total length of the future completed crown or root. However, when quantifying mineralization from linear measurements taken from radiographs (even if not considering the problem of magnification and distortion; see above), we also do not know the total length of the future completed crown and root. So how can we handle the problem of variability in crown and root length in extant human populations, and how do we take it into account to assess dental development? Moreover, inter-and intra-observer errors also exist when taking measurements (only automatic, operator independent, methods are purely nonsubjective). One of us (YH) rated the crown and root mineralization in the total training sample. The other author rated dental development in a fraction of our sample. Disagreement occurred in less than 5% and, at the most, by one stage (usually, at the end of the incisor root formation) (4.2%).

Decomposing a Sequence Into Combinations

Our method aims to quantify patterns of dental development and is based on the hypothesis that any sequence represents a developmental module composed of its hierarchical units (i.e., isolated teeth or groups of teeth) (Figure 1), which exhibit varying degrees of interactions. Indeed, any method aiming to measure patterns yet ignoring this aspect in statistical analyses by assuming independence between teeth is rooted on a circular fallacy. For example, in [START_REF] Smith | Growth and development and its significance for early hominid behaviour[END_REF] method, a dental mean ape is assigned to each developing tooth of an individual from average standards of human and great ape development. This tooth-by-tooth age assessment method obviously assumes independence between teeth because average standards for each tooth are derived from different subsamples of individuals with differing age distributions (e.g., the average standard for the stage observed on the canine is not derived from the same individuals than the average standard observed on the first molar). Indeed, Moorrees et al. (1963a) developed the human standards used by [START_REF] Smith | Growth and development and its significance for early hominid behaviour[END_REF] and mentioned this important limitation on the assessment of dental maturation by acknowledging "possible variation between rates of development of different teeth in a given child". As previously stated by [START_REF] Braga | Non-adult dental age assessment: correspondence analysis and linear regression versus Bayesian predictions[END_REF], the variability across age classes is more likely higher for a given radiographic stage observed on an isolated tooth than for a given mineralization sequence from which this isolated tooth is derived. In other words, individuals will be distributed in a larger age range when considering a single tooth with its radiographic stage than when considering the same tooth within a given sequence. As a consequence, we simply argue that Smith's (1989) age based method produces results which largely depend on the independence assumption. However, this method has the advantage that it is easy to use with a simple program and with only stages of dental development as input. Our aim is also to devise an easily used program, but one that is not based on age assessments and an assumption of independence between growing teeth. Inversely, we use the concept of modularity. Degrees of interactions between isolated teeth, or groups of teeth, can be measured by decomposing any sequence into two subsets with no element in common (Figure 1). Decomposition begins because any sequence comprising 7 developing permanent teeth can be represented by a rearrangement of its seven elements, where none are lost, added, or changed, but distributed into two subsets (each comprising one tooth or more, and noted "sub.k" and "sub.q"; where k + q = n and n = 7), which are combined by a conditional probability. Any sequence of dental development with seven teeth can then be decomposed in a finite number of combinations (as opposed to permutations because in combinations the order of the teeth in each subset is not important, e.g., a subset comprising the canine and the lateral incisor will be identical to the subset comprising the lateral incisor and the canine). Such combinations correspond to conditional probabilities (Figure 1) of observing one subset at radiographic stages (a prior event) as we condition on the attained developmental status of the other subset (a posterior event). The total number of combinations derived from 7 teeth by grouping into two subsets, is 126 (Figure 1). This number is given by the so-called binomial coefficient. Importantly, this decomposition approach is suitable for use with fragmented remains (e.g., fossil hominids). In this case, the total number of combinations will decrease according to the number of preserved teeth.

Calculating Conditional Probabilities

The probability of observing, in our training sample, each of the 126 combinations derived from any complete (7 teeth) sequence of dental development, is calculated using Bayes's rule of conditional probability, with teeth being considered as statistically dependent units (as opposed to "Naive Bayes", with units are considered as independent). The Bayesian approach has recently been used in non-adult dental age assessment from ratings of crown and root formation, considered as statistically dependent units [START_REF] Braga | Non-adult dental age assessment: correspondence analysis and linear regression versus Bayesian predictions[END_REF]. The conditional probability of observing the development status of "subset q" (noted "sub.q"), as we condition on the developmental status of "subset k" (noted "sub.k"), is written as follows: P(subq/subk) = P(sub.k / sub.q) × Pprior (sub.q) ______________________________________________________ [P(sub.k / sub.q) × Pprior (sub.q)] + [P(sub.k / sub.q) × Pprior (sub.q)] The term P(sub.k/sub.q) represents, for subset q, the observed proportion of individuals evincing, in the training sample, subset k. The term P(sub.k/sub.q) represents, for all subsets differing from subset q, the observed proportion of individuals evincing subset k (see [START_REF] Braga | Non-adult dental age assessment: correspondence analysis and linear regression versus Bayesian predictions[END_REF] for more details about the computations). Crucial to this approach is the selection of appropriate prior probabilities ("priors"). A comparison of conditional probabilities from different trials requires the assumption that individuals are derived from samples of balanced distributions of radiographic stages. Training samples never exhibit unbiased distributions of radiographic stages, because they are constructed by "availability" sampling. Consequently, in absence of any specific prior knowledge, we assume that each mineralization stage has the same prior probability to be observed (a so-called "uninformative prior"). The Bayesian procedure is then effective in removing part of the factors of confusion due to sampling bias. As the uniform priors appear in both numerator and denominator, by factorization, they have no influence on the calculation of the conditional probabilities.

Bootstrapping and Mean Conditional Probabilities

Because of finite training sample size, it is possible that some dental mineralization sequences in extant human populations will not be sampled. Bootstrapping is a way of testing the reliability of our training sample, i.e., the degree to which adding or removing individuals in the training sample changes the estimate of variability in sequences of dental development. Therefore, we generated 1000 pseudoreplicate datasets by randomly resampling with replacement 1000 individuals from our training sample (bootstrapping). Consequently, instead of deriving one conditional probability from our total training sample, we obtained distribution of 1000 randomized conditional probabilities for each of the 126 combinations. The mean value of each distribution is then used as a measure of the likelihood of each of the 126 combinations (Figure 2). If this mean probability is higher or equal to 0.75 it is considered to be likely. Inversely, if the mean probability is lower or equal to 0.25 the combination is considered to be unlikely (Figure 2). A mean probability ranging between 0.25 and 0.75 is considered to be not informative (Figure 2).

Let us now consider the example given in Figure 2. Using Demirjian et al.'s (1973) rating system for all permanent teeth (except the third molar), from mesial to distal, the sequence tested is written as follows: G, F, F, E, E, G, E. In this example, 4 combinations (numbered 22, 58, 85 and 121) out of 126, are unlikely, with their mean probabilities falling under 0.25. As shown in Figure 1, the combinations numbered 22, 58, 85 and 121 correspond, respectively, to the following conditional probabilities: CM2 if I1I2P3P4M1, CP4M2 if I1I2P3M1, I1I2P3P4M1 if CM2, I1I2P3M1 if CP4M2. Therefore, from these combinations, we can conclude that, in this sequence, the developmental rating(s) of the canine and/or the second premolar and/or the second molar deviate significantly from what we observe in our training sample. When we now examine the combinations corresponding to the developmental status of each of these three teeth taken separately, versus of all other permanent teeth ( P4,and I1I2CP3P4M1 if M2;corresponding,respectively,to numbers 3,5,7,66,68 and 70), we observe that the lowest probabilities are for the relative dental development of the second molar. Therefore, we conclude that this tooth shows an abnormal developmental status in this sequence.

i.e., C if I1I2P3P4M1M2, P4 if I1I2CP3M1M2, M2 if I1I2CP3P4M1, I1I2P3P4M1M2 if C, I1I2CP3M1M2 if
Figure 2. Example of a sequence of dental development with each of its 126 combinations corresponding to a conditional probability. For each combination, we calculate the mean of 1000 probabilities obtained after 1000 resamplings with replacement.

In order to automate the calculations, Dr. F. Houët (from our department) created a Microsoft Excel 2000® Visual Basic macro to determine the relative frequencies of subsets q and k, and to apply Bayes's rule of conditional probability. This program is easy to use, with only stages of dental development as input. We will make it accessible through a web-site in the future.

Testing Sequences and Cross-Validation

After devising our method, which decomposes any sequence into 126 combinations (Figure 1), and which estimates the generalization error (using bootstrapping) for each of the corresponding 126 conditional probabilities (Figure 2), we need to test each sequence of dental development by using a leave-one-out cross-validation. We test sequences of dental development rather than individuals because most individuals share a common sequence, as illustrated in Figure 2 (see above). After cross-validation, the generalization error of our model can be assessed with a representation of the likelihood and the variability of each of the 126 combinations. For each combination, the likelihood and its variability are represented, respectively, by the mean conditional probability and its coefficient of variation, calculated across all sequences that have been tested (Figure 3). This method can then help us to identify the more conservative and the more plastic combinations, or morphological patterns of covariation during dental growth, in our training sample.

Results

In our training sample, 1926 individuals out of 2089 had their first permanent molars at stages G or H (roots completed or almost completed; see [START_REF] Demirjian | A new system of dental age assessment[END_REF]. These 1926 individuals do not share the same developmental sequence. Instead, they are distributed across 263 different sequences that need to be tested. As shown in Figure 2, some of these 263 sequences with the M1 at stages G or H (117), are represented only once in our training sample. This is the lowest occurrence level for a sequence. These sequences are called "single" as opposed to "multiple". Each single sequence corresponds only to one individual. Each multiple sequence corresponds to two or more individuals. Multiple sequences demonstrate different occurrence levels. Nevertheless, as the occurrence level increases, the number of distinct sequences showing this occurrence level, decreases rapidly in the training sample (Figure 2). In our training sample, when we consider only the 1926 individuals with their first permanent molars at stages G or H, there are 117 distinct single sequences (sequences occurring once, each representing only one individual), 38 distinct double sequences (sequences occurring twice, each representing two individuals, 16 distinct triple sequences, and so on (see Figure 2). When we now consider, in our training sample, the 163 individuals with their first permanent molars before or at stage F, we observe the same trend: as the occurrence level increases, the number of distinct sequences showing this occurrence level decreases rapidly. We now examine which of the 126 combinations are plastic or conserved across the 182 multiple sequences that were tested. In Figure 3, the combinations are numbered and ordered as in Figures 1 and2, with their corresponding mean probabilities (ranging from 0.58 to 0.91) and coefficients of variation (ranging from 20.8 to 43.9). Assuming a normal distribution of all 126 mean probabilities, 95% are within two standard deviations (within 0.62 and 0.90). The following six combinations fall outside this range:

I2CP3P4M1M2 if I1, I1CP3P4M1M2 if I2, I1I2P3P4M1M2 if C, I1I2CP3P4M2 if M1, I1I2CP3P4M1 if M2, CP3P4M1M2 if I1I2 (
respectively, numbers 64 to 66 and 69 to 71) (Figure 3). Five out of six of these combinations represent the developmental status of six permanent teeth given the developmental rating of one isolated tooth. From these results, the premolars do not seem to deviate significantly from the other teeth during dental growth. Interestingly, the sixth combination falling outside the 95% confidence interval of the distribution of all mean probability values corresponds to the combination numbered 71: CP3P4M1M2 if I1I2 (Figure 1). This means that, in any sequence that was tested, the probability of observing the permanent incisors at their relative developmental status, in the training sample, is low. The corresponding low mean probability is associated with the second higher coefficient of variation (41.9) (the highest -43.9is for combination number 6 which is associated with a mean probability falling inside the 95% confidence interval). In other words, these results mean that the mandibular incisors are generally developmentally more plastic compared to other permanent mandibular teeth. Discussion Robert (2004: 129) proposes "…a broad interpretation of evo-devo, one according to which the development of whole organisms and the evolution of their modular parts are deemed the primary analysands, rather than as secondary to the epigenetic expression of purely or primarily genetic potential." Robert (2004: 129) goes on to write: "Such is the force of taking development seriously. A developing system is clearly organized, but in a systemic way; that is, the interrelations between its parts are structured into causal, generative systems with ontogenetic control dispersed throughout these systems". In this evo-devoist perspective, we believe that our results are important. Using radiographic studies [START_REF] Smith | Growth and development and its significance for early hominid behaviour[END_REF][START_REF] Tompkins | Relative dental development in Upper Pleistocene hominids compared to human population variation[END_REF], scholars have asserted that there is "a modern human-like sequence of dental development" that was already present in Homo erectus and Neanderthals [START_REF] Dean | Growth processes in teeth distinguish modern humans from Homo erectus and earlier hominins[END_REF]RamirezRozziandBermúdezdeCastro,2004). It is assumed that "a modern human-like sequence of dental development", as opposed to various human-like sequences, characterizes extant humans and their Upper and Middle Pleistocene fossil relatives. This assumed 'modern humanlike sequence" is then used to compare fossil hominids represented by single, or few, isolated teeth and to understand tooth (especially enamel) growth processes [START_REF] Dean | Growth processes in teeth distinguish modern humans from Homo erectus and earlier hominins[END_REF][START_REF] Ramirez Rozzi | Surprisingly rapid growth in Neanderthals[END_REF]. In particular, this assumption has been used in microanatomical (or histological) analyses to argue that, the anterior teeth can serve as a reliable substitute for the other permanent teeth (e.g., the molars). It has also been used for interpretations of somatic maturation and brain size. However, even if microanatomical studies now allow the study of small samples of children and the quantification of variation in rates of enamel formation of a given tooth (usually, the permanent incisors), we still have not qualified possible variability between rates of development of different teeth in a given child.

The results of our basic analytical approach (Figure 3) suggest that the developmental relationships between, on the one hand, the permanent incisors (one subset), and on the other hand, the canine, premolars and molars (another subset), are the most statistically plastic in extant humans. Therefore, these two subsets of teeth may well represent two hierarchical units within the same developmental module. Obviously, before applying our method to fossil hominids (preferably in single fossil hominids rather than in fossil samples mixing different periods, populations, or even species), our results need to be confirmed by further testing in extant human populations. We should add that the incisor developmental instability and variation demonstrated in this study may arise primarily by parallel variation of separate developmental pathways and simultaneous effects by an extrinsic source of variation. This source could be the developmental constraints due to the growing premaxilla or incisive bone, which may represent another modular unit. In this case, the developing canine, premolars and molars could be considered to co-vary as a hierarchical unit during evolution, arising primarily by strong connections between developmental pathways.

This study aimed to assess variation in the relative sequences (or patterns) of key events during formation of crowns and roots from samples representing different extant human populations, by using a method that does not assume tooth independence. Our training sample is biased because early radiographic stages (from A to F), especially for early developing permanent teeth (incisors and first molar), are underrepresented (Figure 2). Indeed, a high proportion of individuals with their first permanent molars before or at stage F, is represented by single dental mineralization sequence (63 individuals out of 163; i.e., 38.7%) (Figure 2). For individuals with their first permanent molars at stages G or H, this proportion decreases considerably (117 individuals out of 1926; i.e., 6%) (Figure 2). A temporary limitation of our method, which does not assume independence, is that no computations can be made for those individuals represented by unique dental development sequences. For these unique sequences, a "Naïve Bayes" approach can easily be developed but, again, will not be biologically realistic (because of the independence assumption) [START_REF] Heuzé | Chronologie et étiologie de la maturation macrostructurale des dents définitives[END_REF][START_REF] Braga | Non-adult dental age assessment: correspondence analysis and linear regression versus Bayesian predictions[END_REF]. Therefore, in future studies, we need to sample more individuals with their first molars before or at stage F. More sampling will also allow us to quantify possible sexual dimorphism and/or geographic variation of sequences of dental development across extant human populations.

So far, our results are limited to extant humans with their very limited lower facial prognathism coupled with reduced jaw elevator muscles (when compared with chimpanzees). An interpretation of our results in a totally different context of facial morphology and growth would be misleading. Importantly, we still need data on dental development in the pygmy chimpanzee and close attention to the interrelationships between jaw and lower/upper tooth growth in both chimpanzee species in which the patterns of closure of the premaxillary (or incisive) suture are totally different [START_REF] Braga | Chimpanzee variation facilitates the interpretation of the incisive suture closure in South African Plio-Pleistocene hominids[END_REF] but also in extant humans. In this context, we suggest that it would be premature to make direct comparisons, using our method, between extant humans and chimpanzees, or even, between fossil hominids and chimpanzees, in which the premaxilla (containing the alveolar ridge of the four upper incisors), acting as a stabilizing element within the facial skeleton, plays a totally different role in both craniofacial growth and morphology.

Finally, more testing is needed to determine how stable a crown and root-formation sequence is through ontogeny. In other words, would our results be the same using a different sample with a much younger age distribution? Indeed, Moorrees et al. (1963a, b) have demonstrated that later formation stages vary more in their timing than early stages and later-forming teeth vary more in their timing than early forming teeth. However, timing and pattern are two separate things and this trend for timing and its implications for developmental patterns need to be confirmed. This question cannot be addressed simply by examining pairs of developing teeth and testing whether the stage of development of one given tooth increases when another given tooth may increase its range of calcification. This kind of simplistic observation may represent a false correlation caused by the influence of a third single tooth or group of teeth. Therefore, as previously stated, all possible combinations between teeth need to be examined. Nevertheless, our results suggest a developmental instability of the incisors and integration of the developing canine, premolars and molars. We regard this result as consistent with findings suggesting that the first permanent molar provides indications of overall somatic maturity and brain size in living primates [START_REF] Smith | Growth and development and its significance for early hominid behaviour[END_REF].

Conclusions

Dental development studies fall mainly into two methods of observation, which now co-exist but which succeeded each other historically, the first method leading to the hypothesis of an evolutionarily early appearance of a human-like dental development and the later method giving birth its negation. Indeed, scholars first used radiographs and attempted to focus on dental macrostructural development. Pioneering comparative studies [START_REF] Mann | Some Paleodemographic Aspects of the South African Australopithecines[END_REF][START_REF] Skinner | Dental maturation, dental attrition and growth of the skull in fossil hominidae[END_REF][START_REF] Skinner | Atlas of Radiographs of Early Man[END_REF]) suggested that early hominids followed a human-like schedule of dental growth. This hypothesis of an evolutionarily early appearance of a human-like dental development was established on the principal strength of developmental variability in both extant humans and the great apes [START_REF] Mann | Patterns of ontogeny in human evolution: evidence from dental development[END_REF]. The opposing hypothesis was built up with a different method of observation, which was based on the study of incremental growth features in enamel to estimate age at death of fossil hominids [START_REF] Bromage | Re-evaluation of the age at death of immature fossil hominids[END_REF]. These microanatomical observations are based on the record of a precise timescale for enamel developmental events on growth increments. A large number of studies (e.g., [START_REF] Dean | Growth layers and incremental markings in hard tissues: a review of the literature and some preliminary observations about enamel structure in Paranthropus boisei[END_REF][START_REF] Beynon | Radiographic and histological methodologies in estimating the chronology of crown development in modern human and great apes: a review, with some applications for studies on juvenile hominids[END_REF][START_REF] Dean | Growth processes in teeth distinguish modern humans from Homo erectus and earlier hominins[END_REF] now consider this latter approach as a more secure model on which to base studies of the paleodemography, growth and maturation of early hominids. This led to the more general consensus that truly human-like dental development emerged relatively late in human evolution because all early hominids documented thus far, including some earliest fossils attributed to Homo, did not resemble modern humans in their development. These two views can now be reconciled because the conceptual framework of modularity enables us to build new methods to achieve an agreement. More precisely, modularity and EVO-DEVO [START_REF] Raff | The Shape of Life: Genes, Development and the Evolution of Animal Form[END_REF][START_REF] Bolker | Modularity in development and why it matters to evo-devo[END_REF] offer a new conceptual framework and new perspectives to provide a picture of the variability in sequences of dental development, as well as in enamel and dentine growth processes, as possible foundations for morphological changes during human evolution. Indeed, our hypothesis of modularity in dental development and our method derived from this concept can serve as a basis for identifying and studying patterns of dental growth and, importantly, for comparisons populations or species that do not need to be hedged with age assessments based on a tooth-by-tooth analysis and its fallacies. At this stage, we suggest that our first comprehensive quantification of variability in patterns of development of permanent teeth in three extant human populations opens an EVO-DEVO perspective. Therefore, we aim to contribute to EVO-DEVO in the future by examining modularity/patterning in fossil hominid samples and by further testing (e.g., influence of sex, age and geographic origin) of the conclusions presented here.

Figure 1 .

 1 Figure 1.

Figure 3 .

 3 Figure 3. Leave-one-out cross validation in our training sample: 182 multiple sequences (see Figure 2) have been tested.

Table 1 .

 1 Geographic composition of the sample. ( * some of them being, at least in part, of European origin)

	Girls	Boys		Girls	Boys
	242	214	South of France	54 to 192	47 to 189	Four grand parents from Western Europe
	392	291	South of France	42 to 192	43 to 192	Unknown, or mixed, geographic origin
	249	135	Southern Iran	69 to 191	77 to 190	Iran
	323	244	Ivory Coast	36 to 181	45 to 190	Ivory Coast *

Acknowledgments

We would like to express our gratitude to Shara Bailey and Jean-Jacques Hublin for inviting us to Leipzig to contribute to this special edition and their help throughout its preparation. This study was funded by the French Ministry of Research ("Action Concertée Incitative") and the GDR 2152 of the CNRS. We thank Drs. Chabadel, de Brondeau, Garde, de la Talonnière, Roux, Boileau, Chouvin, Hassid, Blanca, Thévenet, Frapier (France), Sonan (Ivory Coast) and Gueramy (Iran), for access to the orthopantomographs. We also thank F. Houët and B. Dutailly for their help in developing the program for the calculations and an anonymous referee who helped us to improve our discussion.


Ilya Karuseichyk

Giacomo Sorelli

Vyacheslav Shatokhin

Mattia Walschaers

Nicolas Treps

Exploiting separation-dependent coherence to boost optical resolution

The problem of resolving point-like light sources not only serves as a benchmark for optical resolution but also holds various practical applications ranging from microscopy to astronomy. In this research, we aim to resolve two thermal sources sharing arbitrary mutual coherence using the spatial mode demultiplexing technique. Our analytical study includes scenarios where the coherence and the emission rate depend on the separation between the sources, and is not limited to the faint sources limit. We consider the fluorescence of two interacting dipoles to demonstrate that the dependence of emission characteristics on the parameter of interest can boost the sensitivity of the estimation and noticeably prolong the duration of information decay.

I. INTRODUCTION

The role of coherence in imaging problems has been intensively studied and discussed for decades [START_REF] Considine | Effects of Coherence on Imaging Systems[END_REF][START_REF] Nayyar | Two-point resolution of Gaussian aperture operating in partially coherent light using various resolution criteria[END_REF][START_REF] Aert | Resolution of coherent and incoherent imaging systems reconsidered -Classical criteria and a statistical alternative[END_REF]. This topic became especially interesting with the formulation of the imaging problem in terms of parameters estimation [START_REF] Helstrom | Image Restoration by the Method of Least Squares[END_REF][START_REF] Helstrom | Resolvability of Objects from the Standpoint of Statistical Parameter Estimation[END_REF]. A widely studied form of this problem is the estimation of the separation between two point sources of light, which not only serves as a benchmark for optical resolution but also has numerous practical applications ranging from microscopy to astronomy [START_REF] Baddeley | Biological Insight from Super-Resolution Microscopy: What We Can Learn from Localization-Based Images[END_REF][START_REF] Tsang | Resolving starlight: a quantum perspective[END_REF]. The parameter estimation approach to this simple imaging problem allows to find an ultimate quantum limit for the resolution [START_REF] Tsang | Quantum Theory of Superresolution for Two Incoherent Optical Point Sources[END_REF][START_REF] Nair | Far-Field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit[END_REF][START_REF] Lupo | Ultimate Precision Bound of Quantum and Subwavelength Imaging[END_REF][START_REF] Sorelli | Quantum limits for resolving Gaussian sources[END_REF]. The Fisher information of the spatial-mode demultiplexing (SPADE) technique was shown to saturate the quantum limit for incoherent [START_REF] Tsang | Quantum Theory of Superresolution for Two Incoherent Optical Point Sources[END_REF][START_REF] Nair | Far-Field Superresolution of Thermal Electromagnetic Sources at the Quantum Limit[END_REF][START_REF] Rehacek | Optimal measurements for resolution beyond the Rayleigh limit[END_REF] and partially coherent [START_REF] Tsang | Resurgence of Rayleigh's curse in the presence of partial coherence: comment[END_REF] faint sources. The sensitivity measure based on the method of moments [START_REF] Gessner | Metrological Nonlinear Squeezing Parameter[END_REF] allowed to study a more general case of bright incoherent [START_REF] Sorelli | Moment-based superresolution: Formalism and applications[END_REF][START_REF] Sorelli | Optimal Observables and Estimators for Practical Superresolution Imaging[END_REF], fully coherent [START_REF] Karuseichyk | Resolving mutually-coherent point sources of light with arbitrary statistics[END_REF] and entangled [START_REF] Zhang | Quantum super-resolution for imaging two pointlike entangled photon sources[END_REF] sources, demonstrating non-trivial scaling of the sensitivity with the brightness of the thermal sources and showing that the SPADE approach remains quantum optimal in this more general scenario. Moreover, the method of moments provides a simple estimation strategy for the parameter of interest that relies solely on the measured mean intensities and does not require maximum likelihood estimation or other heavy computations. At the same time, the sensitivity obtained from the method of moments saturates the Cramer-Rao bound in the limit of faint sources.

In this study we analyze the case of partially coherent sources that sparked hot debates in the community [START_REF] Tsang | Resurgence of Rayleigh's curse in the presence of partial coherence: comment[END_REF][START_REF] Larson | Resurgence of Rayleigh's curse in the presence of partial coherence[END_REF][START_REF] Larson | Resurgence of Rayleigh's curse in the presence of partial coherence: reply[END_REF][START_REF] Hradil | Quantum Fisher information with coherence[END_REF][START_REF] Hradil | Exploring the ultimate limits: super-resolution enhanced by partial coherence[END_REF][START_REF] Liang | Coherence effects on estimating two-point separation[END_REF][START_REF] Liang | Quantum Fisher information for estimating N partially coherent point sources[END_REF][START_REF] Kurdzialek | Back to sources -the role of losses and coherence in super-resolution imaging revisited[END_REF][START_REF] Wang | Quantum-limited resolution of partially coherent sources[END_REF]. We aim to resolve two bright thermal sources * ilya.karuseichyk@lkb.upmc.fr sharing arbitrary mutual coherence, including cases when the latter depends on the separation between the sources. This situation was partially covered in Ref. [START_REF] Wang | Quantum-limited resolution of partially coherent sources[END_REF] using the approximation of faint sources and assuming a lossless optical system. In this contribution, we instead consider sources with arbitrary brightness, since the approximation of faint sources is not always applicable in practical situations. Furthermore, previous findings indicate that the lossless model tends to underestimate sensitivity in cases where the initial brightness of the sources is known [START_REF] Tsang | Resurgence of Rayleigh's curse in the presence of partial coherence: comment[END_REF][START_REF] Kurdzialek | Back to sources -the role of losses and coherence in super-resolution imaging revisited[END_REF]. To address this limitation, our model explicitly incorporates losses in the imaging system.

Our approach, employing the method of moments, extends to scenarios where both mutual coherence and the emission rate of the sources are dependent on the separation. This method allows us to address a practically significant issue of resolving two reflectors illuminated by a source of light with finite coherence width, where a faint source limit is not always applicable. We show that the finite coherence width of the illumination can significantly enhance the optical resolution, however, this enhancement becomes less pronounced with the growing number of coherently emitted photons, i.e. for a more extended coherence time of the illumination source.

In the second example, we resolve the separation between two interacting dipole emitters. We show that the dependence of emission characteristics on the separation can drastically boost the estimation sensitivity. As a result, one can obtain useful information about the dipoles' separation on a timescale significantly larger than the radiative lifetime of the excited state of the dipole. For the late stages of decay, even though the probability of detecting a photon is very small, the detection event contains a lot of information because of the strong interaction-induced correlations between the emitters. We show that this effect is resistant to the dephasing of the emitters and persists even under fast

Ƹ 𝑠 1 Ƹ 𝑠 2 𝑑 𝜅 ො 𝑎 1 ො 𝑎 𝑚 𝜎 𝑃𝑆𝐹 𝛾 FIG.
1. The scheme of the two point sources resolving. The sources emit light into orthogonal modes ŝ1,2 with a mutual coherence γ. The light propagates through the lossy imaging system with transmissivity κ and PSF width σP SF before being decomposed into spatial modes fm(x) with corresponding field operators âm and subsequently detected.

dephasing.

II. THE MODEL OF EMITTERS

We consider a traditional optical scheme for singleparameter imaging in which the only estimated parameter is the separation d between two point-like light sources (Fig. 1). All other parameters of the sources are assumed to be known, i.e. the imaging setup can be aligned in such a way that both sources lie on its X-axis and the centroid of the sources lies at the coordinate origin. Thus the problem can be treated in one dimension with coordinates of the sources being x j = (-1) j d/2. The sources emit light into orthogonal modes with field operators ŝ1,2 . Each of these modes describes a spherical wave emitted by the corresponding source during its coherence time. We only consider equally bright sources emitting thermal light, including cases of correlated emission. Since these stochastic sources do not have any external phase reference, the emitted state should be symmetric over the global phase shift. This implies that ⟨ŝ j ⟩ = 0 and ⟨ŝ j ŝk ⟩ = 0. Thus the state is fully characterized by the first-order coherency matrix

Γ (1)ŝ jk = ⟨ŝ † j ŝk ⟩ = N γN γ * N N , ( 1
)
where N is the average number of photons emitted by each source per its coherence time. Assuming the coherence time of the sources to be fixed we refer to N as brightness and define sources with N ≪ 1 as faint sources. Note, however, that the number of photons emitted by such sources per integration time of the detection system can still be large if the coherence time of the sources is short compared to the integration time. Unlike most research on this topic, our approach allows for studying bright thermal sources with N > 1, demonstrating the non-trivial scaling of the separation estimation sensitivity with this parameter. The complex number γ = γ 0 e iϕ with -1 ≤ γ 0 ≤ 1 and 0 ≤ ϕ < π stands for the mutual coherence between the sources. In the general case, γ can depend on the separation. The presence of a non-diagonal part in the coherency matrix corresponds to correlations of the quadratures

V jk = 1 2 { Qj , Qk } = =     N + 1 2 0 N Re γ N Im γ 0 N + 1 2 -N Im γ N Re γ N Re γ -N Im γ N + 1 2 0 N Im γ N Re γ 0 N + 1 2     , (2)
with quadrature vector ⃗

Q = (q 1 , p1 , q2 , p2 ) T , with qj = (ŝ † j + ŝj )/ √ 2 and pj = i(ŝ † j -ŝj )/ √ 2.
In this case mutual coherence originates from thermal correlations, as opposed to coherent states where mutual coherence comes from a nonzero mean-field.

III. FIELD TRANSFORMATION THROUGH THE IMAGING SYSTEM

Light emitted by the point sources propagates through a diffraction-limited imaging system with transmissivity κ (see Fig. 1). The corresponding transformation of the field operators in the paraxial approximation can be represented as [START_REF] Lupo | Ultimate Precision Bound of Quantum and Subwavelength Imaging[END_REF]

ĉj = √ κŝ j + √ 1 -κv j , (3)
where operators ĉj are associated with the images of the sources u 0 (x -x j ), and u 0 (x) is the point spread function (PSF) of the imaging system, which is assumed to be translationally invariant, and thus κ does not depend on the position of the source. Here we assume, without losing generality, that the magnification factor of the imaging system equals 1. Field operators vi are associated with auxiliary mutually non-orthogonal modes [v i , v † j ] ̸ = δ ij in the vacuum state. These auxiliary modes describe losses of light, that unavoidably appear due to the finite aperture of the imaging system and in the general case lead to modes non-orthogonality [ĉ j , ĉ † k ] ̸ = δ jk . The measurements are performed in the spatial modes f m (x) with the corresponding field operators âm , which can be expressed as [START_REF] Karuseichyk | Resolving mutually-coherent point sources of light with arbitrary statistics[END_REF]

âm = j f * m (x)u 0 (x -x j )dx ĉj + v′ m , (4)
where v′ m are non-normalized field operators of auxiliary vacuum modes.

In most cases, the PSF is symmetric, i.e. u 0 (-x) = u 0 (x). Thus it is often reasonable to use a symmetric mode basis for the measurement, i.e. perform SPADE in a basis of modes with well-defined parity

f m (-x) = (-1) m f m (x). ( 5
)
Photon counting in such symmetric bases was shown to be quantum optimal in other scenarios [START_REF] Rehacek | Optimal measurements for resolution beyond the Rayleigh limit[END_REF][START_REF] Sorelli | Moment-based superresolution: Formalism and applications[END_REF][START_REF] Karuseichyk | Resolving mutually-coherent point sources of light with arbitrary statistics[END_REF]. In this case Eq. ( 4) reduces to âm =

2 j=1 (-1) jm β m ĉj + v′ m , (6)
where

β m = dx f * m (x) u 0 x - d 2 . ( 7
)
The average photon numbers in the measurement modes reads

N m = ⟨ Nm ⟩ = ⟨â † m âm ⟩ = ξ m + ζ m γ 0 cos ϕ, (8)
where

ξ m = 2κN β 2 m , ζ m = 2κN β 2 m (-1) m . (9)
The second term in (8) accounts for light interference arising from partial coherence. The reader may notice that the average detected photon numbers depend solely on the real part of the degree of mutual coherence, Re γ = γ 0 cos ϕ. However, as we show in the following, the imaginary part Im γ does influence the photon counting statistics and, consequently, changes the sensitivity of the parameter estimation.

All the field characteristics after passing through the imaging system contain sources' original brightness N only in combination with the transmissivity factor κ. Thus important parameter of this problem is the combination κN , which defines the brightness of the sources imaged with a given optical system.

The total number of photons detected in K measurement modes equals

N D = K m=0 N m = 2κN (A + δγ 0 cos ϕ). (10)
Coefficients A and δ do not depend on the measurement basis if it is full (K → ∞) and symmetric (5):

A = K m=0 β 2 m ----→ K→∞ 1 (11) δ = K m=0 β 2 m (-1) m ----→ K→∞ u 0 x - d 2 u 0 x + d 2 dx.
(12) In this limit (K → ∞), the coefficient δ depends only on the shape of the PSF u 0 (x) and the value of the separation d, representing the overlap between images of the sources. In the case of a full measurement basis all the photons, that passed through the imaging system are being detected. Generally, the total number of detected photons N D [START_REF] Lupo | Ultimate Precision Bound of Quantum and Subwavelength Imaging[END_REF] depends on the separation d through the overlap δ. One can leverage this dependence to achieve a more accurate estimation of the separation d, only provided that the losses are correctly accounted for in the model [START_REF] Tsang | Resurgence of Rayleigh's curse in the presence of partial coherence: comment[END_REF][START_REF] Kurdzialek | Back to sources -the role of losses and coherence in super-resolution imaging revisited[END_REF]. If the theoretical model is defined with an arbitrary loss factor, measuring only one observable does not provide useful information.

IV. MEASUREMENT SENSITIVITY

To describe the sensitivity of the SPADE measurement, we use the method of moments [START_REF] Gessner | Metrological Nonlinear Squeezing Parameter[END_REF]. It sets the bound for estimators d(N m ) that are based on the mean values of the observables N m = ⟨ Nm ⟩:

∆ 2 d ≥ 1 µM , ( 13
)
where µ is the number of measurement repetitions (assumed to be large) and the sensitivity M is defined as

M = K m,n=0 ∂N m ∂d Γ -1 mn ∂N n ∂d , (14)
with Γ being the covariance matrix of the measured observables, with elements

Γ mn = ⟨ Nm Nn ⟩ -⟨ Nm ⟩⟨ Nn ⟩. (15)
Moreover, the method of moments presents an optimal estimator [START_REF] Gessner | Metrological Nonlinear Squeezing Parameter[END_REF]. This estimator relies on the linear combination of mean values N m , and asymptotically saturates the bound given by ( 13) for a large enough number of measurement repetitions µ.

Exploiting the property of non-displaced Gaussian states

⟨â † m âm â † n ân ⟩ = ⟨â † m âm ⟩⟨â † n ân ⟩+ + ⟨â † m â † n ⟩⟨â m ân ⟩ + ⟨â † m ân ⟩⟨â m â † n ⟩, (16)
and Eq. ( 6) we find the photon number covariance matrix

Γ = diag ⃗ Y 1 + ⃗ Y 2 + ⃗ Y 1 ⃗ Y T 1 + ⃗ Y 2 ⃗ Y T 2 + + 1 -γ 2 0 2 sin 2 ϕ ⃗ ζ ⃗ ζ T , ( 17
)
where diag(⃗ x) is a diagonal matrix with diagonal elements equals x n and

⃗ Y 1,2 = 1 ± γ 0 2 ⃗ ξ ± ⃗ ζ cos ϕ , (18)
with ⃗ ξ and ⃗ ζ defined in Eq. ( 9). In the limit of faint sources, κN ≪ 1, the covariance matrix [START_REF] Karuseichyk | Resolving mutually-coherent point sources of light with arbitrary statistics[END_REF] becomes diagonal Γ nm = δ nm N m , and the sensitivity M (14) coincides with the Fisher information. In this case only the real part of the mutual coherence Re γ = γ 0 cos ϕ defines the measurement statistics.

In the more general scenario, determining the sensitivity M requires inverting the non-diagonal photon number covariance matrix Γ [START_REF] Karuseichyk | Resolving mutually-coherent point sources of light with arbitrary statistics[END_REF]. It is possible to do analytically by applying the Sherman-Morrison formula [START_REF] Sherman | Adjustment of an Inverse Matrix Corresponding to a Change in One Element of a Given Matrix[END_REF] three times, but the resulting general expressions are very bulky. Therefore, we decided to focus on several special cases for which the inversion simplifies.

V. CONSTANT MUTUAL COHERENCE

A. Real-valued coherence (γ = γ0) First, we consider separation-independent real-valued coherence γ = γ 0 ∈ R, i.e. ϕ = 0. For this case, the photon number covariance matrix reduces to:

Γ = diag ⃗ Y 1 + ⃗ Y 2 + ⃗ Y 1 ⃗ Y T 1 + ⃗ Y 2 ⃗ Y T 2 . ( 19
)
The inversion of this matrix gives

Γ -1 mn = δ mn (ξ n + γ 0 ζ n ) -1 -t mn , (20)
where

t mn =    2/h 1 , if m and n are both even 2/h 2 , if m and n are both odd 0, otherwise, (21)
with

h 1,2 = 2 + 2κN (1 ± γ 0 )(A ± δ), (22)
and A and δ are defined in Eqs. [START_REF] Sorelli | Quantum limits for resolving Gaussian sources[END_REF] and [START_REF] Rehacek | Optimal measurements for resolution beyond the Rayleigh limit[END_REF].

If the mutual coherence γ and the brightness N are independent of the separation, i.e. ∂γ 0 /∂d = 0, and ∂N/∂d = 0, the derivative of the measured signal is

∂N m ∂d = 2 N m β m ∂β m ∂d . (23)
Then the normalised sensitivity of the separation estimation is

M Re 2κN = ∆k 2 -γ 0 β- -2κN (A ′ + δ ′ ) 2 (1 + γ 0 ) 2 2h 1 + (1 -γ 0 ) 2 2h 2 , (24)
where

∆k 2 = 4 K m=0 ∂β m ∂d 2 , (25)
β = -4 K m=0 ∂β m ∂d 2 (-1) m , (26)
δ ′ = ∂δ ∂d = 2 K m=0 (-1) m β m ∂β m ∂d , (27)
A ′ = ∂A ∂d = 2 K m=0 β m ∂β m ∂d . (28)
All these quantities are measurement basis independent for full symmetric bases:

∆k 2 ----→ K→∞ ∂u 0 (x) ∂x 2 dx, (29)
β ----→ K→∞ ∂u 0 (x -d/2) ∂d ∂u 0 (x + d/2) ∂d dx (30)
In this scenario, the sensitivity expressed in (24) coincides with the quantum Fisher information [START_REF] Lupo | Ultimate Precision Bound of Quantum and Subwavelength Imaging[END_REF][START_REF] Sorelli | Quantum limits for resolving Gaussian sources[END_REF]. This implies the possibility of constructing an estimator d(N m ) based solely on a linear combination of the mean photon count values in any basis with the parity defined in Eq. ( 5) and achieving the ultimate resolution limit dictated by the quantum Cramér-Rao bound.

To calculate numerical values of the sensitivity M (24) we consider the so-called soft-aperture imaging system with Gaussian PSF

u 0 (x) = 1 2πσ 2 P SF 1/4 exp - x 2 4σ 2 P SF , (31)
where σ P SF is the PSF width. Then for the full symmetric measurement basis, we have

δ = exp - d 2 8σ 2 P SF , ∆k 2 = 1 4σ 2 P SF (32) β = σ 2 P SF -(d/2) 2 4σ 4 P SF exp - d 2 8σ 2 P SF . (33)
The sensitivity calculated with Eq. ( 24) is shown in Fig. 2. The first plot corresponds to a small number of emitted photons and matches other results obtained in the same limit [START_REF] Tsang | Resurgence of Rayleigh's curse in the presence of partial coherence: comment[END_REF][START_REF] Kurdzialek | Back to sources -the role of losses and coherence in super-resolution imaging revisited[END_REF]. With an increasing number of photons, there is a reduction in the sensitivity for the separations d ≈ σ P SF . This reduction is also present in the QFI and occurs due to the quadratic term of the noise in the thermal statistics [START_REF] Lupo | Ultimate Precision Bound of Quantum and Subwavelength Imaging[END_REF]. However, as we show in the next section, this drop does not occur even for bright thermal states, if their mutual coherence degree γ is not a real number.

B. Imaginary coherence ( γ = iγ0)

Another interesting example that we consider is that of a purely imaginary degree of coherence γ = iγ 0 ∈ I, which corresponds to a relative phase ϕ = π/2 between the sources. In this case, the measured mean photon numbers N m do not contain an interference term, i.e. it is the same as for a pair of incoherent sources. However, the presence of coherence does influence the covariance matrix which only depends on γ 2 0 and thus is not affected by its sign. The inverse of this matrix reads

Γ = diag ⃗ ξ + 1 + γ 2 0 2 ⃗ ξ ⃗ ξ T + 1 -γ 2 0 2 ⃗ ζ ⃗ ζ T , (34)
Γ -1 mn = δ mn ξ n - (-1) m+n η 1 + η 2 -2κN δ((-1) n + (-1) m ) η 1 η 2 -(2κN δ) 2 , (35)
where η 1,2 = 2κN A + 2(1 ± γ 2 0 ) -1 . In this case, the normalized sensitivity takes the form

M Im 2κN = ∆k 2 -2κN η 1 (δ ′ ) 2 + η 2 (A ′ ) 2 -2κN δA ′ δ ′ η 1 η 2 -(2κN δ) 2 . ( 36
)
It is plotted in Fig. 3 for different values of the mutual coherence and sources' intensities.

One can see that in the limit of faint sources the sensitivity

M Im 2κN κN →0 ----→ ∆k 2 (37)
is the same for any imaginary degree of mutual coherence, including the case of incoherent sources. However, for bright incoherent thermal sources the sensitivity per photon drops significantly in the sub-Rayleigh region d < σ P SF , while for sources with imaginary mutual coherence, this drop is smaller or doesn't occur at all for perfect coherence γ = i. Generally, if imaginary coherence is close to the complex unity then the sensitivity can be expressed as

M Im 2κN = ∆k 2 -2κN (δ ′ ) 2 (1 -γ 0 ) + O((1 -γ 0 ) 2 ). ( 38
)
This means that the presence of correlations between sources can improve the sensitivity of the separation estimation even in cases when no interference is visible in the mean values N m . One can use a simple intuition to interpret this observation: two incoherent thermal sources do not share any correlations, while mutually coherent sources are correlated in the number of photons, even if they do not interfere due to the mutual phase ϕ = π/2. Utilizing these correlations allows to cancel part of the intensity noise in the measurement modes, resulting in a more precise estimation of the separation.

VI. PARAMETER-DEPENDENT COHERENCE

Now let us consider the situation when the mutual coherence γ depends on the separation d of the sources. Below we will discuss in details physical examples of such systems, but first, we analyze the general case of real-valued separation dependent mutual coherence

γ(d) = γ 0 (d) ∈ R.
This dependence leads to an additional term in the derivative

∂N m ∂d = 2 N m β m ∂β m ∂d + ∂γ 0 ∂d ζ m . ( 39
)
and in the sensitivity

M γ0(d) = M Re + ∆M γ0(d) , (40)
where

∆M γ0(d) 2κN = 2γ ′ 0 δ ′ 1 h 1 + 1 h 2 + A ′ 1 h 1 - 1 h 2 + + (γ ′ 0 ) 2 A -δγ 0 1 -γ 2 0 -2κN (A + δ) 2 2h 1 + (A -δ) 2 2h 2 , (41)
and M Re is defined in Eq. ( 24). The additional sensitivity ∆M γ0(d) originates from the fact that in the case of separation-dependent coherence, the measurement results change faster with changing separation, and one can do estimation with higher precision. The extra sensitivity ∆M γ0(d) typically has maxima around the maxima of γ ′ 0 , i.e. in the regions of fastchanging mutual coherence.

In the limit of faint sources, the additional sensitivity takes the form

∆M γ0(d) 2κN κN →0 ----→ 2γ ′ 0 δ ′ + (γ ′ 0 ) 2 A -γ 0 δ 1 -γ 2 0 . ( 42
)
This expression is valid for an arbitrary complex value of the mutual coherence γ if one replaces γ 0 with Re γ.

A. Finite coherence width of the illumination

In the first example, we investigate the estimation of separation between two reflective objects. When illuminated with light of finite coherence width, the mutual coherence of the reflected light becomes dependent on the separation between the reflectors (see Fig. 4). In many cases, the transverse coherence of the illuminating light can be effectively approximated by a Gaussian function. In such scenarios, the field reflected by two small reflectors exhibits mutual coherence

γ 0 (d) = exp - d 2 2ω 2 c , (43)
where the coherence width ω c can be adjusted by changing the optical parameters of the illumination system. Using the Gaussian PSF (31), we calculate the full sensitivity (40) of the reflectors' separation estimation (Fig. 5). The red lines on the plots correspond to the spatially coherent illumination. One can see that separationdependent coherence, compared to fully coherent case, results in significantly higher sensitivity for the separations close to the coherence width of the illumination source d ≈ ω c , thus it is desirable to use illumination with coherence width of the order of the measured separations. A similar increase in the resolution, when the coherence width of the source matches the size of the features of the studied object, was observed earlier experimentally and numerically for the case of quantum imaging with pseudo-thermal light [START_REF] Mikhalychev | Efficiently reconstructing compound objects by quantum imaging with higherorder correlation functions[END_REF].

The coherence width of the light expands as it passes through a diffractive imaging system. Thus, if the same optical system (with PSF width σ P SF ) is used for the imaging and the illumination of the object, the minimal achievable coherence width of the illumination is ω c ≥ 2σ P SF . This case was studied in Ref. [START_REF] Wang | Quantum-limited resolution of partially coherent sources[END_REF] in the faint source limit, where an analogous increase of the SPADE sensitivity was demonstrated for the separations around the Rayleigh limit (in comparison to spatially coherent illumination). However, as we show in Fig. 5, for brighter sources this increase is much less noticeable. In the extreme case of very bright sources, additional sensitivity per photon (41) vanishes

∆M γ0(d) 2κN ----→ N →∞ 0. ( 44
)
This is important to take into account for practical imaging with pseudo-thermal sources, since this scenario often does not fit into the faint sources approximation. Therefore, it is crucial to aim for the utilization of an illumination source with a short coherence time to ensure that the number of photons N reflected by each imaged object per coherence time is small. One could achieve better sensitivity in the sub-Rayleigh regime via a smaller coherence width, which is achievable in two possible scenarios: using an independent illumination scheme with narrower PSF (for example, if the illumination source is located closer to objects than the detection apparatus) or in the case of sources' mutual coherence originating from their interaction, which is considered in the next section.

B. Interacting emitters

Another possible scenario where the mutual coherence of the emitted light depends on the distance between the sources is the case of interacting emitters (see Fig. 6).

As an example, we consider two identical, dipole-dipole interacting two-level systems prepared initially in their excited states. In this scenario, each dipole emits precisely 1 photon during the decay of the excited state. For simplicity, we assume that the dipole moments are parallel to each other and the line connecting them is orthogonal to the main optical axes of the imaging system.

Ƹ 𝑠 1 Ƹ 𝑠 2 𝑑 ො 𝑎 1 ො 𝑎 𝑚 𝜎 𝑃𝑆𝐹 𝜅 𝑔 12 FIG. 6.
Dipoles emit partially coherent light due to the interaction.

The time evolution of the quantum state of the two coupled dipoles in the interaction picture can be described by the following master equation [START_REF] Lehmberg | Radiation from an $N$-Atom System. I. General Formalism[END_REF][START_REF] Agarwal | Quantum statistical theories of spontaneous emission and their relation to other approaches[END_REF]

d dt ρ = -if 12 [σ + 1 σ- 2 + h.c., ρ] + 1 2 2 j,l=1 g jl 2σ - j ρσ + l -σ+ j σ- l ρ -ρσ + j σ- l + 2 j=1 η j 2σ + j σ- j ρσ + j σ- j -σ+ j σ- j ρ -ρσ + j σ- j , ( 45
)
where σ± j are the dipoles' transition operators. The first term in Eq. ( 45) is proportional to the collective frequency shift f 12 of the energy levels [START_REF] Agarwal | Quantum statistical theories of spontaneous emission and their relation to other approaches[END_REF][START_REF] James | Frequency shifts in spontaneous emission from two interacting atoms[END_REF]. The second term describes the spontaneous emission of the individual emitters (j = l) and their collective radiative decay (j ̸ = l). In the case of identical dipoles,

g 11 = g 22 = 1/τ, ( 46
)
g 12 = 3 2 1 τ sin z z + cos z z 2 - sin z z 3 , ( 47
)
where τ is the natural lifetime of the excited state of an individual dipole, and z = 2πd/λ, with λ the wavelength of the dipole transition. The final term in the master equation accounts for the dephasing of the dipoles with the rates η i (i = 1, 2), due to their different local environments.

The radiation of the dipoles is described in the transient temporal modes ŝ1,2 (t), of duration ∆t ≪ τ , centered around time t. For these modes, the timedependent coherency matrix can be expressed through the atomic dipole correlators as

Γ (1)ŝ(t) jk = ⟨σ + j σ- k ⟩ ∆t τ . ( 48
)
The coherency matrix of the excited modes remains unaffected by the collective shifts f 12 . However, as discussed in Appendix A, the emission spectrum is influenced by f 12 . Consequently, the measurement results and the sensitivity of the measurement are independent of f 12 exclusively in the considered case of frequency-blind measurements.

Cross-correlation of the atomic operators ⟨σ + 1 σ-2 ⟩ vanishes if the collective decay rate g 12 = 0. The crucial role of these correlations in the emergence of quantum coherence was previously discussed in other contexts [START_REF] Shatokhin | Coherence turned on by incoherent light[END_REF][START_REF] Ames | Theory of multiple quantum coherence signals in dilute thermal gases[END_REF].

Both the unitless emission rate Ṅ = N τ /∆t = Γ depend on both time and separation, however, we will not explicitly indicate it in our notations.

1. Model without dephasing (η1 = η2 = 0) First, we examine the evolution of the dipoles' state without taking into account the dephasing process. The analytical solution for the dipole correlators in this scenario is presented in Appendix A. The associated emission rate Ṅ and the degree of mutual coherence γ for different time instances t are depicted in Fig. 7 as functions of the separation d normalized by the wavelength λ.

One can see that at an early stage of the emission (t = 0.02τ ) the dipoles fluoresce almost independently, with the emission rate close to that of an individual dipole, Ṅ ≈ e -t/τ , and with almost no coherence, γ ≈ 0. However, after some time, the dipole-dipole interaction creates correlations and the difference with the individual dipole emission becomes evident in both the emission rate [START_REF] Devoe | Observation of Superradiant and Subradiant Spontaneous Emission of Two Trapped Ions[END_REF] and the mutual coherence. Crucially, this difference depends on the separation between the dipoles, which allows for a measurement of the separation with higher accuracy. Without accounting for the decoherence process, the mutual coherence exhibits a particularly strong dependence on separation during the late emission stages (t > τ ) for separations d ≈ n λ 2 , n ∈ N. Therefore, we can expect heightened sensitivity in estimating separation around these points. Note that due to the geometrical symmetry of the problem the mutual coherence γ stays real at any moment, i.e. γ(t) ∈ R ∀t. However, all the conclusions below can be generalized to the case of the complex degree of mutual coherence by replacing γ 0 with Re γ.

The mean number of photons emitted by each dipole during a short time interval ∆t ≪ τ is low (N ≪ 1). Furthermore, the imaging system is characterized by high losses, κ ≪ 1. Therefore, irrespectively of the statistics of the sources, the detection statistics can be described by a Poisson distribution which also approximates that of faint thermal sources discussed earlier. Besides, the radiation emitted by the two dipoles in the far field represents a superposition of two spherical waves and therefore is locally equivalent to the radiation emitted by two-point sources (see Appendix A). Thus, in the paraxial approximation, the model developed in the previous sections can be applied to the problem of resolving dipoles.

The only modification required to Eq. ( 39) is the addition of an extra term, associated with the separationdependent emission rate: This extra-dependence results in an additional term ∆M N (d) in the sensitivity. In the limit κN = κ Ṅ ∆t/τ ≪ 1 it takes the form

∂N m ∂d = 2 N m β m ∂β m ∂d + ∂γ 0 ∂d ζ m cos ϕ + N m N ∂N ∂d . (49)
∆M N (d) = 2κN ′ N N ′ (A + γδ) + 2N (A + γδ) ′ , ( 50
)
where X ′ = ∂X/∂d. Combining all the contributions we get the normalized sensitivity per unit time in the form

Ṁ = 4σ 2 P SF 2κ τ ∆t M Re + ∆M γ0(d) + ∆M N (d) . ( 51
)
This quantity is plotted in Fig. 8 and9.

As one can expect, at an early stage of emission the problem resembles the resolving of two incoherent sources (i.e. non-interacting emitters). In the context of independent emitters, the sensitivity per unit time is directly proportional to the emission rate, exhibiting an exponential decrease, as indicated by Ṅ ≈ e -t/τ (depicted by the dashed black lines in Fig. 8 and9). On the contrary, in the case of interacting emitters, the sensitivity rate may even increase over time due to the cumulative effects of interaction. The specific values of the sensitivity depend on the ratio between the separation d, PSF width σ P SF and the wavelength λ. We present plots for the cases of σ P SF = λ and σ P SF = 5λ. In the paraxial approximation, the PSF width of the imaging system typically significantly exceeds the wavelength, i.e. σ P SF > λ. Therefore, the case σ P SF = λ stretches the boundaries of our model. Nevertheless, we investigate this case to explore the model's limits and offer an illustrative example that is easy to analyze. In Fig. 9, one can observe a pronounced enhancement in sensitivity within regions where the emission characteristics (the mutual coherence and the emission rate) exhibit rapid changes with a change in separation. Meanwhile, the sensitivity remains low around the extrema of the emission characteristics. In typical cases, the sensitivity for interacting dipoles is significantly larger and decaying significantly slower compared to non-interacting dipoles. Even after more than 15 lifetimes of the excited state, when the probability of photon detection is very low, the average sensitivity rate remains notably high due to the substantial amount of information in "late" detection events.

It is important to note that the moment-based sensitivity provides a limit for a local estimation strategy. The presence of multiple narrow peaks in the sensitivity plot signals the potential degeneracy of the estimator. Therefore, when dealing with a low-resolution optical system (σ P SF ≫ λ), one may require an increased number of measurement repetitions to reach the saturation of the bound [START_REF] Tsang | Resurgence of Rayleigh's curse in the presence of partial coherence: comment[END_REF].

In general, the interaction between the emitters induces an entanglement between them [START_REF] Tanaś | Entangling two atoms via spontaneous emission[END_REF] and creates temporal correlations of the emitted light [START_REF] Peshko | Quantum noise radar: superresolution with quantum antennas by accessing spatiotemporal correlations[END_REF]. However, if the losses in the imaging system are large κ ≪ 1, the detection events for different time intervals can be considered independent. The sensitivity of independent detection events is additive thus the total sensitivity can be calculated as

M T ot = 1 τ ∞ 0 Ṁ dt, ( 52
)
where Ṁ is defined in (51). The result of this calculation is shown in Fig. 10. One can see that the interaction effects can increase the total sensitivity of the separation estimation by several orders of magnitude around the points d ≈ n λ 2 , n ∈ N, where the mutual coherence of the emission strongly depends on the separation, if the dephasing effect is not present. In between these points, for the separations around d ≈ n -1 2 λ 2 , n ∈ N, both the emission rate and mutual coherence exhibit extrema. Consequently, the sensitivity is not heightened by the separation-dependent emission characteristics around these specific points. In realistic scenarios, the process of building up coherence through the dipole-dipole interactions competes with the coherence loss caused by dipoles' individual dephasing. The master equation (45) with dephasing can be solved analytically, however the solution is too long and cumbersome to be presented within this paper. Nevertheless, it demonstrates that the coherency matrix for our system only depends on the average dephasing rate η = (η 1 + η 2 )/2 of the dipoles. In the following, we discuss the impact of dephasing for different values of η. Figure 11 illustrates the emission characteristics (the emission rate and the mutual coherence) at the time instance t = 5τ . In the weak dephasing regime (η = 0.2/τ ), where the dephasing time is substantially longer than the natural lifetime of the excited state τ , the effect of dephasing has a limited influence on the emission process. On the other hand, when the dephasing process is faster than the emission (η = 2.0/τ ), one can observe a significant decrease in the mutual coherence γ.

As expected, the reduction in emission coherence due to the dephasing effect leads to a weaker sensitivity boost. Fig. 12 illustrates the total sensitivity of the separation estimation defined in (52). The top plot highlights the robustness of the considered scheme to weak dephasing, where the sensitivity is only minimally affected, compared to the dephasing-free case, shown in Fig. 10. Conversely, the bottom plot in Fig. 12 reveals that strong dephasing significantly diminishes the sensitivity boost arising from dipole interaction (pay attention to the different scale of the y-axis of this plot). However, even in this scenario, the sensitivity can be several times higher compared to non-interacting dipole emission.

Curiously, in all the considered cases the sensitivity boost is more pronounced for smaller separations, where the interaction between the emitters is stronger. As a result, in the context of single-parameter estimation, one can achieve higher precision in estimating smaller separations compared to larger ones.

Finally, we note that in the case of frequency-resolved measurements, the collective frequency shift f 12 may come into play, such that one may be able to locally estimate the dipoles' separation via the emission spectrum. Therefore, by integrating spectral and spatial approaches, the sensitivity could be further enhanced by introducing frequency-resolving detection into our scheme.

VII. CONCLUSION

We have analyzed the problem of resolving partially coherent thermal sources with SPADE measurement. To estimate the sensitivity of a separation estimation, we utilized the method of moments, which allowed us to make no assumptions about the brightness of the sources. We have found analytical expressions for the sensitivity including the cases of separation-dependent mutual coherence and emission rate. We studied two specific examples of separation-dependent coherence: a reflection of light coming from a finite-coherence-width illumination source and creating mutual coherence due to the interaction of the emitters. In both cases, we demonstrate the possibility of a significant boost in the separation estimation sensitivity due to the additional mechanism of the parameter encoding to the problem. Our analysis shows that for efficient resolving of the reflecting objects one needs to use an illumination source with a narrow coherence width (on the order of the reflectors' separation) and a short coherence time. This ensures the full advantage of the separation-dependent coherence of the reflected light. Examining the interacting emitters, we demonstrate that the sensitivity of separation estimation can be increased by several orders of magnitude compared to independent emitters. This enhancement arises from the separation-dependent mutual coherence and emission rates of the interacting dipoles. The effect remains robust in the presence of weak dephasing of dipoles. Although strong dephasing decreases the observed resolution boost, it does not eliminate it entirely.

FIG. 2 .FIG. 3 .

 23 FIG. 2. The normalized sensitivity MRe = MRe4σ 2 P SF /2κN for different real values of the mutual coherence γ and numbers of emitted photons N multiplied by the transmissivity κ vs. separation d.

FIG. 4 .

 4 FIG.[START_REF] Helstrom | Image Restoration by the Method of Least Squares[END_REF]. Objects reflect the light from the common thermal illumination source with finite coherence width ωc.

FIG. 5 .

 5 FIG. 5. Normalized sensitivity M γ 0 (d) = M γ 0 (d) 4σ 2 P SF /2κN in case of the separation-dependent mutual coherence with Gaussian profile (43).

( 1 )

 1 ŝ(t) 11 τ /∆t and the degree of mutual coherence γ = Γ

FIG. 7 .

 7 FIG. 7. Properties of the dipoles radiation: the normalised emission rate Ṅ and the degree of mutual coherence γ vs. separation d. The dashed line corresponds to the time t = 0.02τ , solid line t = 5τ . No dephasing (η = 0).

FIG. 8 .

 8 FIG.8. Normalized sensitivity rate Ṁ (51) of two dipoles' separation estimation vs. time t. The red line corresponds to σP SF = λ, the blue line to σP SF = 5λ, and the dashed line to non-interacting dipoles.

FIG. 9 .

 9 FIG. 9. Normalized sensitivity rate Ṁ (51) of two dipoles' separation estimation vs. separating d. The red line corresponds to σP SF = λ, the blue line to σP SF = 5λ, and the dashed line to non-interacting dipoles.

2 .

 2 Model with dephasing (η1, η2 ̸ = 0)

FIG. 10 .

 10 FIG. 10.Normalized total sensitivity MT ot (52) of two dipoles' separation estimation vs. separation d. The red line corresponds to σP SF = λ, the blue line to σP SF = 5λ, and the dashed line to non-interacting dipoles. No dephasing (η = 0).

FIG. 11 .

 11 FIG. 11. Emission characteristics: normalized emission rate Ṅ and degree of mutual coherence γ at the time moment t = 5τ . The red line corresponds to the dephasing-free model η = 0, the blue line -weak dephasing η = 0.2/τ , the grey line -strong dephasing η = 2.0/τ .

FIG. 12 .

 12 FIG.12. Normalized total sensitivity MT ot (52) for the model with dephasing. Top inset -weak dephasing η = 0.2/τ , bottom inset -strong dephasing η = 2.0/τ .

ACKNOWLEDGMENTS

I.K. acknowledges support from the PAUSE program. This project has received funding from the European Defense Fund (EDF) under grant agreement EDF-2021-DIS-RDIS-ADEQUADE (n°101103417). Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union. Neither the European Union nor the granting authority can be held responsible for them. This work was carried out during the tenure of an ERCIM 'Alain Bensoussan' Fellowship Programme.

Appendix A: Coherency matrix of dipole emission

In the far field of the dipole, the positive-frequency part of the emitted field operator in point ⃗ r reads [START_REF] Scully | Quantum Optics[END_REF] ⃗

where ⃗ p is a dipole moment. This emission mode does not have the spherical symmetry, however, in the far field, it is locally indistinguishable from a spherical wave. Thus in a paraxial ap-proximation, dipole emitters can be considered as point sources.

Due to the proportionality (A1) between the emitted field operator ⃗ E and the dipole transition operators σ-, the coherency matrix of the emitted field can be expressed as (48), where all the constant factors, except for ∆t, are omitted since they will later be included in the transmissivity factor κ. This factor represents the coupling between the emission modes and image modes of individual sources. Explicit calculation of (48), giving the evolution (45), result in

(g 11 + g 12 ) 2 e -g12t + (g 11 -g 12 ) 2 e g12t -4g 2 12 e -g11t , (A2)

where g jl is defined in (46). One can see, that in the non-interacting limit g 12 -→ 0 the emission rate decay exponentially Ṅ = e -t/τ and the mutual coherence does not appear γ = 0. Note that the unitary component of the master equation ( 45), featuring a collective frequency shift f 12 , does not impact the coherency matrix of the excited modes. Nevertheless, the spectra of the excited modes are influenced by f 12 . Consequently, in instances of frequency-resolving measurement, the dependence f 12 (d) introduces an additional mechanism for parameter encoding and, in general, enhances sensitivity.