Bez Nicolas
email: nicolas.bez@ird.fr

Renard Didier
email: didier.renard@mines-paristech.fr

Ahmed-Babou Dedah
email: dedah.ahmed-babou@ird.fr

Nicolas Bez

Empirical Orthogonal Maps (EOM) and Principal Spatial Patterns: Illustration for Octopus Distribution Off Mauritania Over the Period 1987-2017

Keywords: Spatiotemporal data, Factorization, Principal maps, Distance between maps, Dimension reduction

Analysis of spatiotemporal observations often leads to decomposition of the problem into a spatial part multiplied by a temporal part (factorization). Principal component analyses produce factors that are temporally uncorrelated but that remain spatially correlated, leading to incomplete factorization. Min-max autocorrelation factors developed many years ago are adapted here to ecological applications, leading to empirical orthogonal maps (EOMs). EOMs owe their name to the fact that they are indeed an enhancement of empirical orthogonal functions which extract the spatial patterns that explain most of the variability of a set of spatiotemporal observations indexed by time. Application on a time series of 61 scientific monitoring surveys targeting octopus distribution off the Mauritanian coast indicates that ten basic maps are sufficient to recover 68% of the total variability, and that the first two EOMs explain 38.4% of this variability. This manuscript clarifies the concept of orthogonality between factors in a spatial context. This provides the conditions for using Euclidean distance between spatial distributions, which in turn supports the reduction of a large set of spatial distributions into a small subset of basic spatial distributions explaining most of the variability within the set of input maps.

Several approaches, stochastic or deterministic, parametric or empirical, can be explored to analyze data with dual spatial and temporal dimensions. The complexity of the landscape of possible approaches is complicated by the semantic that, sometimes, introduces ambiguity. For instance, it is common to say that some methods are spatial while, in fact, they are not. Our criterion is to consider that a method is spatial (respectively temporal), if the exchange of two observations in space (respectively in time) modifies the result. Spatial representation of outputs is not sufficient to claim that the method that generated the outputs is spatial. For instance, empirical orthogonal functions -EOFs [START_REF] Lorenz | Empirical orthogonal functions and statistical weather prediction[END_REF] produce outputs that can be represented on a geographical scale. However, EOFs do not account for spatial structures and spatial auto-correlation that exist in the observations and are not, per se, spatial statistics. This was precisely the improvement brought in by [START_REF] Switzer | Min-max autocorrelation factors for multivariate spatial imagery[END_REF] with the min-max auto-correlation factors -MAF that explicitly account for spatial correlation.

When dealing with spatio-temporal observations, an objective shared by a large part of the different approaches is to factorize the spatio-temporal problem into a spatial part multiplied by a temporal part. Based on the stochastic partial derivative equation -SPDE framework recently developed by [START_REF] Lindgren | An explicit link between gaussian fields and gaussian markov random fields: the stochastic partial differential equation approach[END_REF], [START_REF] Thorson | Spatial factor analysis: a new tool for estimating joint species distributions and correlations in species range[END_REF] conducted a spatial factor analysis -SFA. This approach uses a reduced number of orthogonal random fields with SPDE characteristics. Each latent random field is a stationary random Gaussian field characterized by a Matérn spatial covariance function (related to the distance between observations). In such a model, spatial independence between factors is guaranteed by construction and their number is constrained to be small for inference and parsimony considerations. This framework is thus quite appealing in that it allows identifying few factorial maps summarizing the input spatio-temporal signal. A possible drawback of this approach is that the factors, both their numbers and their shapes, are model based. Once the model types are chosen (i.e., Matérn spatial auto-correlations functions), the parameters are inferred using integrated nested Laplace approximation -INLA algorithms making the overall approach efficient. [START_REF] Switzer | Min-max autocorrelation factors for multivariate spatial imagery[END_REF] proposed a factorization based on so-called min-max auto-correlation factors -MAFs to filter out the noise of a series of multi-channel spatial imagery data. Even though their approach uses variogram computations, no variogram model is required and their approach is model-free. The MAF procedures is a sequence of two principal component analyses -PCA, and can be fully developed empirically, that is, without any parametric assumption. Since their seminal work, several studies have been developed using MAF decompositions (e.g., [START_REF] Shapiro | Extracting time trends from multiple monitoring sites[END_REF] notably in marine ecology [START_REF] Fujiwara | Identifying interactions among salmon populations from observed dynamics[END_REF][START_REF] Wackernagel | Using min/max autocorrelation factors of survey-based indicators to follow the evolution of fish stocks in time[END_REF][START_REF] Petitgas | Analysing temporal variability in spatial distributions using min-max autocorrelation factors: sardine eggs in the bay of biscay[END_REF], and several papers have contributed to further describe the connection between MAF and the linear model of coregionalization (LMC) [START_REF] Desbarats | Geostatistical simulation of regionalized pore-size distributions using min/max autocorrelation factors[END_REF][START_REF] Vargas-Guzman | Computational properties of min/max autocorrelation factors[END_REF]. These latter references established also that for particular situations like the intrinsic correlation models, the absence of cross-correlation at some particular distance propagates to all spatial distances insuring full orthogonality between factors.

When applied to multivariate time series, the different variables, taken in columns of the database, are often systematically measured over time taken as the rows (preferably regularly, but not always). In these cases, temporal MAFs can be readily applied, as in [START_REF] Solow | Detecting change in the composition of a multispecies community[END_REF] or in [START_REF] Wackernagel | Using min/max autocorrelation factors of survey-based indicators to follow the evolution of fish stocks in time[END_REF]. The situation is more problematic considering time series of spatial distributions. In these cases, each column is a spatial distribution ; the rows corresponding the geographical sites. Precisely because it is a model-free approach, the main drawback of MAF is that they need isotopic observations, that is that data should be observed at the same spatial points over time [START_REF] Wackernagel | Using min/max autocorrelation factors of survey-based indicators to follow the evolution of fish stocks in time[END_REF]. This is seldom achieved for spatial survey following a scientific random design. In these cases, observations must be transformed into isotopic sets of data prior to factorization by ad hoc methods. For instance, [START_REF] Petitgas | Analysing temporal variability in spatial distributions using min-max autocorrelation factors: sardine eggs in the bay of biscay[END_REF] used migration to the nearest point of a common grid prior to the computation of MAFs. However, when the sampling protocol varies from one survey to the other, or when the number of samples fluctuates between surveys, heterotopy cannot be solved by a simple migration; instead, a spatial interpolation (e.g., kriging) is required to return to isotopy. Regular and isotopic grids correspond to particularly suitable situations. In these cases, as we will see in this manuscript, the first PCA of the MAFs is nothing but an EOF. This paves the way for a new nomenclature that is justified below.

MAFs are ordered by say, descending, spatial auto-correlation intensity (the first MAFs are, by construction, those with the strongest auto-correlation at small geographical distance). This is the reason why the first MAFs suppress noise components. However, in the ecological context for instance, what matters is rather the relative importance of the different spatial factors with regards to the variability of the observations, which corresponds to a different perspective. The objective is no longer to filter out the noise but to build the spatio-temporal patterns that best explain the time series of distribution maps, with the possible consequence that the most important part of the spatio-temporal signal could be the noise (nevertheless, this is important to know from an ecological point of view). In this paper, MAFs are thus reformulated so that the first factors of the decomposition, hereafter called empirical orthogonal maps -EOMs, make it possible to reproduce a reasonable portion of the variability of the initial data, which is not necessarily ensured by MAF.

The objectives of this manuscript are then threefold. First, the calculation of EOMs are detailed, clarifying the standardization steps which is important for the interpretation of the outputs. Second, the percentage of variance explained is defined which gives the basis for the distinction between MAF and EOM. Third, the orthogonality between factors is enlarged to encompass mean orthogonality for spatial scales not restricted to the shortest ones. This opens the possibility to compute Euclidean distance between spatial distributions and to map them into small dimensional spaces.

The developments are illustrated by an application on a time series of sixty one scientific monitoring surveys of the octopus spatial distribution off the Mauritanian coast.

METHOD

The whole framework can be developed outside of any probabilistic framework even though the use of random function formalism could be relevant. In this context, capital letters are used to denote matrices. The factors are obtained by linear combination of a set of spatially regular distribution maps represented by the matrix Z of dimension S × T where S indicates the spatial dimension (e.g., the number of geographical locations sampled each time) and T the temporal dimension (e.g., the number of times the spatial distribution is observed). The matrix notation provides a connection with the continuous space-time framework where Z[s,t] represents the value at site s, s = 1, . . . , S and at time t,t = 1, . . . , T .

Without loss of generality, Z is supposed to be centered. Its standardized version is denoted Ż = Z • D s -1 where D s -1 is the T × T diagonal matrix of the inverse standard deviations of the column of Z.

SPATIALLY NON-CORRELATED FACTORS

The factorization consists of a sequence of two PCAs. When dealing with regular gridded data, the first PCA is nothing but an EOF. It is based on the eigen elements of Γ 0 , the correlation matrix of Z, that is, the variance-covariance matrix of Ż. The correlation considered here is the correlation between observations made at the same site but at different time, that is for a 0 distance in the geographical space, hence the notations used (subscript 0 means "related to distance equals to 0"). If we denote λ 0 the vector of the eigenvalues and P 0 the T × T matrix of the corresponding eigenvectors of Γ 0 , the EOFs are given by

F 0 = ŻP 0 , (1)
where F 0 is a matrix with the same dimension as Z (S × T ). The variances of the factors are given by the corresponding eigenvalues. Their sum equals T , that is the number of input spatial distributions. This first PCA is a projection of the standardized raw data in an orthogonal (but not orthonormal) base formed by the empirical orthogonal factors. In the spatio-temporal representation, the orthogonality between the factors refers to non-correlation of the factor values at the same geographical points. However, these factors may exhibit spatial correlations in the sense that their spatial covariances may be different from zero. The factors of an EOF are statistically, but not spatially, orthogonal.

The second PCA thus aims at constructing new factors with no spatial correlation for a given spatial distance (which usually corresponds to the distance to the nearest neighbor in the case of systematic spatial sampling or to the average distance to the nearest neighbor in irregular cases). The factors of the first PCA are first standardized, that is divided by the square root of their eigenvalues

Ḟ0 = F 0 D λ -1/2 0 .
(2)

The second PCA is then based on the diagonalization of the variance-covariance matrix of their spatial increments for a given reference distance h = r, that is on the eigen decomposition of (twice) the matrix of variogram and cross-variogram values between standardized EOFs for the distance h = r. Denoting Γ r this T × T matrix, λ r the vector of its eigenvalues and P r the matrix of its eigenvectors, the factors associated to the reference distance h = r are

F r = Ḟ0 P r = Z D s -1 P 0 D λ -1/2 0 P r = ZΦ r , (3)
where F r is an S × T matrix, and where the linear operator Φ r to transform directly the input raw data into factors is defined by

Φ r = D s -1 P 0 D λ -1/2 0 P r . (4)
In this expression, the matrix Φ r explicitly describes the sequence of operations required to build the factors:

1. first, the data are standardized (D s -1 ); 2. then, their statistical variance-covariance matrix is diagonalized and the data are projected in this orthogonal space (P 0 ); 3. then, the factors constituting this orthogonal space are normalized (D λ -1/2 0 ); 4. and, finally, (twice) the variogram-cross-variogram matrix at distance h = r of the projected data is diagonalized and the projected data are projected in this new orthogonal space (P r ).

By construction, the final factors F r , have a unit variance, and are uncorrelated locally (h = 0) and at distance h = r (the proof is given in [START_REF] Switzer | Min-max autocorrelation factors for multivariate spatial imagery[END_REF]. Each factor (F r [.,t],t = 1, . . . , T ) represents indeed a spatial distribution, which is a linear combination of the input spatial distributions. As they are obtained by linear combinations of the input spatial distributions, this can be reversed given that Φ r is invertible with Φ -1 r = Ψ r (see the section "Practical considerations" below). Each input spatial distribution can, in turn, be expressed as a linear combination of the full set of the factors without any approximation (back-transformation) by

Z = F r Ψ r .
(5)

This can be summarized by the following way and back transformation scheme

Z Φ r -⇀ ↽ - Ψ r F r . (6)

ORDERING THE FACTORS AND DIMENSION REDUCTION: MAF VERSUS EOM

The back-transformation equation Eq. ( 5) offers the possibility to use only a subset of the factors when back-computing the spatial distributions from the factors. Approximation or filtration of the input spatial distributions are then obtained by the linear combination of, say, the n-first factors

Ẑ1:n = F r [., 1 : n]Ψ r [1 : n, .]. ( 7
)
The question is thus to order the factors in accordance with the objectives of the dimension reduction.

MAF

While the primary objective of [START_REF] Switzer | Min-max autocorrelation factors for multivariate spatial imagery[END_REF] was to remove the noise part of a set of images, they ordered the factors, called MAFs, by increasing variograms values at distance h = r in order to favor the spatially most regular factors and to remove the factors with pure or strong nugget structure. The eigenvalues λ r equal (twice) the variogram of the final factors at distance h = r, so that for two factors ranked i and j, one gets

γ i, j (r) = 1 2 λ r [i] if i = j (simple variogram) 0 if i ̸ = j (cross variogram) . ( 8
)
This is the basis for ordering MAF which are organized in increasing order of eigenvalues λ r so that the first MAF gets the strongest spatial structures.

EOM

An alternative objective is to select as few factors as possible that explain as much variability of the observations as possible, further called EOMs. Obviously, the larger the n parameter, the better the approximation. In standard PCA, the total variance or the total inertia I is equal to the sum of the variances of the input variables, given by the trace of the variance-covariance matrix of I = tr(C Z ), and the percentage of variance explained by the first factors is equal to the sum of their eigenvalues over the total inertia. This cannot be directly transposed to the spatial factors that are produced by a sequence of two PCAs. Alternatively, one can consider, one by one, the sets of T distributions obtained by the back-transformation of a single factor, say factor i

Ẑi = F r [., i]Ψ r [i, .]. (9)
The percentage of variance explained by this factor is hereafter defined by

p i = tr(C Ẑi ) tr(C Z ) = tr(Ψ r [i, .] t Ψ r [i, .]) tr(C Z ) = ∑ T k=1 Ψ r [i, k] 2 tr(C Z ) , (10)
This makes it possible to rank and re-arrange the factors (eigenvalues and eigenvectors) according to the percentage of variance that they explain and select the ones that reproduce the largest part of the input variance. This ordering leads to empirical orthogonal maps -EOMs.

There is a priori no reason for the two arrangements, MAF and EOM, to coincide.

MATHEMATICAL ORTHOGONALITY AND STATISTICAL INDEPENDENCE

In spatial statistics, the orthogonality implies the absence of correlation between factors at any given possible distance (i.e., not only for zero distance). Strictly speaking, this means that the set of EOMs forms an orthogonal basis if and only if they are spatially uncorrelated. However, by construction, EOMs get zero covariance only at distance 0 and at distance h = r. While it is a step ahead towards spatial non-correlation compared to traditional EOFs, the EOMs are not fully spatially orthogonal and do not form a basis sensus stricto. In some cases, notably in case of intrinsic correlations models, absence of cross-correlation at some distance between principal components propagates to all distances [START_REF] Goovaerts | Spatial orthogonality of the principal components computed from coregionalized variables[END_REF][START_REF] Desbarats | Geostatistical simulation of regionalized pore-size distributions using min/max autocorrelation factors[END_REF][START_REF] Rondon | Teaching aid: minimum/maximum autocorrelation factors for joint simulation of attributes[END_REF]). In these particular cases, statistical non correlation is equivalent to full spatial orthogonality. However, this relies on particular models of coregionalization.

In an empirical approach which MAF are in essence, one can not rely on model's characteristics. An empirical alternative is however suggested below.

In their seminal paper, [START_REF] Switzer | Min-max autocorrelation factors for multivariate spatial imagery[END_REF] computed MAF through a sequence of two PCA, the second one using spatial increments of a fixed spatial distance (e.g., the pixel size for grid data). A natural extension is to consider a range of distances instead of precise distance for this second PCA. A weak definition of orthogonality is thus envisaged here. Hereafter, regionalized variables are said to be weakly orthogonal if the mean value of all their cross-variograms for all possible distances is 0, that is if

γ i̸ = j (h) = 0, ∀h. (11)
We thus suggest to compute EOMs using a range of spatial distances (h ∈ [0, R]) rather than a reference distance (h = r) as in the original paper of [START_REF] Switzer | Min-max autocorrelation factors for multivariate spatial imagery[END_REF]. To do so, the second PCA relies on Γ [0,R] the T × T matrix of (twice) the mean variogram and mean cross-variogram values for all possible distances between 0 and R. Using the above notation, this would lead to the following factorization

F [0,R] = ZΦ [0,R] with Φ [0,R] = D s -1 P 0 D λ -1/2 0 P [0,R] . (12)
The only difference between this decomposition (Φ [0,R] ) and the former one (Φ r ) relies on the second PCA (P [0,R] instead of P r ); the first PCA, non spatial, remains unchanged. Weak non-correlation can be applied to the full set of the EOMs or to the n-first ones in case of dimension reduction, in order to diagnose if the space where the input spatial distributions are projected is more or less spatially orthogonal.

DISTANCE BETWEEN SPATIAL DISTRIBUTIONS

In the EOMs framework, each input spatial distribution is represented by the vector of coefficients of their EOMs decompositions. If the n-first EOMs are uncorrelated, these coefficients give the coordinates of the input spatial distributions in the orthonormal space defined by the n-first EOMs. Each dimension of this space is defined by an EOM, that is, a basic spatial distribution. The first factorial plan, which corresponds to the first two EOMs, is the two-dimensional space defined by the two principal spatial distributions that explain most of the variability of the input spatial distributions. Each input distribution can be represented by a point in this two-dimensional space.

This can be generalized to any dimensional space and opens the use of a distancebased metric to compare distribution maps or clustering techniques. For instance, hierarchical ascending classification or k-means can be used to group spatial distributions whose decompositions in the basis of EOMs are similar. Given the above discussion on orthogonality (in the mathematical and statistical senses), the use of Euclidean distance between spatial distributions is as relevant as the orthogonality of the EOMs is effective. In this context, there is a strict equivalence in considering:

1. the coefficients of the decomposition of an input spatial distribution on n basic spatial distributions/EOMs, 2. the coefficients of the approximation of an input spatial distribution by the linear combination of n basic spatial distributions/EOMs and, 3. the coordinates of an input spatial distribution in an orthogonal space made of n basic spatial distributions/EOMs.

STANDARDIZED AND NON-STANDARDIZED EOMS

Back-transformation, dimension reduction and approximation can refer to the standardized input data. This allows analyzing the shape of the input spatial patterns irrespective to their level of variability. In this case, spatial distributions are similar if they have the same patterns, up to a multiplicative value. The EOMs are the same as their computation is based on standardized data, but the transformation and back-transformation matrices differ slightly as matrix D s -1 must be removed.

Ż = F r Ψr , (13)
with

Φr = P 0 D Λ -1/2 0 P r = Ψ -1 r . ( 14
)
The percentage of correlation explained by the n-first EOMs may however differ from the percentage of variance explained.

Practical considerations

EOMs are built for a given reference distance h = r. However one often needs to use some tolerance around the reference distance to account for sampling sites that are not regularly spaced, or a given reference distance lag h ∈ [0, R] to ensure weak non-correlation.

The signs of the eigen elements are purely conventional but coherent between eigenvectors and their eigenvalues. Therefore, their interpretation must be established jointly.

Empirical variance-covariance matrices are not always positive definite. In particular, when T ≥ S, that is when the number of surveys is larger than or equal to the number of sampling sites per survey, the variance-covariance matrices is singular and cannot be inverted. In this case, its eigen elements do not exist and the EOM decomposition is not possible.

The time series of input spatial distributions may not be regular in time. It is also worth mentioning that indices t could be interchanged without modifying the EOM outputs. EOM is not a method that is explicitly temporal. However, unlike EOF that are not a spatial method (sites could be interchanged without modifying the EOF outputs), EOM are spatially explicit.

Each term is important and meaningful: "empirical" indicates that no parametric assumption is required, "orthogonal" refers to non correlation, and "map" specifies that the factors of the decomposition are maps or spatial distributions (contrary to PCA where the factors are variables).

All the computations were performed under R using the package RGeostats (MINES-ParisTech/ARMINES 2021). The scripts and the data required to reproduce the analysis are available here: https://github.com/abambad/EOM_Proj.

APPLICATION

DATA

EOMs are used to analyze the time series of sixty-one (T = 61) octopus surveys made with the research vessels N'Diago and Al-Awam during the periods 1987-1996 and 1997-2017 respectively. Each survey follows a stratified random sampling based on three latitudinal strata (Fig. 1). The average number of samples is 102 samples per survey with some surveys having only few tens of samples. In each sampling site, the density of octopus is provided in number of individuals per swept area (on average 0.055 km²). An inter-calibration experiment between the two research vessels was carried out for the period 1987-89 in order to make the data series homogeneous by taking into account the change of fishing gear that took place in 1989 [START_REF] Gascuel | Decline of demersal resources in north-west africa: an analysis of mauritanian trawl-survey data over the past 25 years[END_REF]. The timing of the surveys is not regular over years but there is at least one survey per year.

The geographical locations of the sampling sites of each survey are drawn at random and are thus different from one survey to another, with surveys with a low spatial coverage. So, prior to EOMs computations, survey data were interpolated by ordinary block kriging (Chilès and Delfiner 2012) on a regular 0.1°x 0.1°grid restricted to the polygon of presence of octopus with a kriging neighborhood of 0.75°. Given the small latitudes of the study area no projection was required (no significant space deformation at these latitudes). The omni-directional experimental variograms got reasonable spatial structures (see Supplementary information). The number of active grid cells is S = 341 (T < S).

As all regression techniques, kriging is smoothing. This means that the kriging maps do not have the same level of variability than the raw input observations. This favors using standardized EOMs that allow comparing and grouping the surveys only based on the shapes of their spatial distributions.

RESULTS

The first EOM alone explains 28.5% of the overall variability (Fig. 2) and four basic EOMs were enough to recover more than half of the input variability. The ranking based on the percentage of explained correlation is not fully consistent with the ranking based on the variogram value (see for instance, the sixth EOM whose spatial regularity is not intermediate between that of the fifth and the seventh ones). The local variability, that is the variogram value at the reference distance of 0.1°, starts increasing after the tenth EOMs. Moreover the ten first EOMs carry 68% of the overall variability. This is considered as a good compromise for dimension reduction (n = 10).

The cross-variogram values of the ten selected EOMs are equal to 0 for the reference distance (i.e., h=0.1°) and are equal to zero on average when mixing all pairs of EOMs together (Fig. 3a) which is consistent with a weak non-correlation. However, their fluctuations clearly increase for larger distances. When EOMs are computed for a large interval of reference distance, say for distances between 0°and 1°, their cross-variogram values are strongly reduced around 0 on average for all distance classes between 0°and 1°(Fig. 3b) indicating that they form an orthogonal basis in the strong sense. However, this produces EOMs whose spatial structures are rapidly become pure noise with low descriptive power (Fig. 2). The first factorial space, that is the space formed by the two first EOMs, explains 38.4% of the overall input correlation. In this factorial space, the surveys display two clear groups with little overlap (Fig. 4). This is further investigated through a hierarchical ascending classification (HAC) based on the 10 first EOMs using the Ward distance. The HAC underlines the existence of two groups of surveys with similar spatial patterns (Fig. 5) that strongly matches the climatic seasons (Fig. 5; accuracy = 68%). While the first three EOMs get clearly and statistically (ANOVA with p.value < 10 -3 ; Fig. 6) different scores for the two clusters, ANOVA diagnostics are also statistically significant for EOMs ranked 6, 8, 9 and 10. Finally, the vectors of the average scores per cluster are used to estimate the mean spatial distribution of each cluster (Fig. 4).

DISCUSSION

Being defined by the sequence of two PCAs, EOM is fundamentally an empirical approach. The fact that it refers to variogram in the second PCA does not change this fact. Incidentally, what is minimized is not the value of a variogram model at a given distance, but the value of empirical covariance between the increments of two EOFs for a given reference geographical distance. The recourse to random function is thus external to the EOM: either before, in order to estimate or to simulate values on fixed sampling sites over the study period to get isotopic data base, or after, to map the EOMs if needed. The recourse to random functions and to model of coregionalization help understanding when the factors are spatial factors sensus stricto [START_REF] Goovaerts | Spatial orthogonality of the principal components computed from coregionalized variables[END_REF][START_REF] Desbarats | Geostatistical simulation of regionalized pore-size distributions using min/max autocorrelation factors[END_REF]. Being empirical, EOM has the drawback that it can only be performed for sets of sampling sites (spatially regular or not) that are systematically observed over time (isotopy). Factorizing a set of het-Fig. 1 Sampling protocol. Typical survey data (March, 2015). Octopus densities are in number of individuals per swept area (on average 0.055km 2 ) erotopic variables (i.e., variables that are not observed at the set sets of geographical points) by linear combinations remains an open question. Spatial covariance (simple and cross covariance) can deal with heterotopic variables [START_REF] Wackernagel | Using min/max autocorrelation factors of survey-based indicators to follow the evolution of fish stocks in time[END_REF]. However, the fundamental nature of factorization being to make a linear combination of the variables at the same points, the problem cannot be fully solved by the recourse to covariance.

Similar to EOFs that are not spatial in their construction but that can be represented spatially, EOMs are not temporal but can be represented temporally. It is thus an abuse of language to say that EOM offers a spatio-temporal approach.

In statistics, the orthogonality of the factors (PCA, EOF, etc) refers to the absence of mutual correlations. In a spatial context, this means the absence of correlation at the same geographical points. After the second PCA, the orthogonality is extended to a given geographical distance but not to all possible distances. So EOMs do not reach a full orthogonality by construction. The only known case were orthogonal factors in the statistical sense are also fully spatially uncorrelated is when all covariances are proportional, also called intrinsic correlation model [START_REF] Chilès | Geostatistics, modeling spatial uncertainties[END_REF][START_REF] Goovaerts | Spatial orthogonality of the principal components computed from coregionalized variables[END_REF][START_REF] Desbarats | Geostatistical simulation of regionalized pore-size distributions using min/max autocorrelation factors[END_REF]. This model is very peculiar and very specific. In the present analysis, we have shown that a weak orthogonality can be considered when mutual correlations are not null one by one, but on average. We also have indicated that enlarging the distance interval used in the EOM computation can help reaching orthogonality. However, the EOMs that are obtained in Fig. 2 Values of the omnidirectional variogram of the EOMs for the grid cell distance (0.1°) as a function of the percentage of correlation explained by the r-first EOMs. In black, the EOMs corresponding to a reference distance equal to the grid cell size (h=0.1+/-0.05). In red, the EOMs obtained when h=0.5+/-0.5. The first ten EOMs are depicted by plain circles Fig. 3 Cross variogram values for the first distance lags for the first ten EOMs. Left: EOMs used in the study corresponding to a reference distance equal to the grid cell size (h=0.1+/-0.05). Right, EOMs obtained when using a large interval reference distance (h=0.5+/-0.5) Fig. 4 Factorial space made of the first two EOMs. In this space, the survey are located according to their standardized scores in the EOMs factorization. The two EOMs associated to each dimension of the factorial space are represented. The mean distribution of each group is also represented. Surveys that happened during the hot season are denoted "H" and "C" for the cold season. Polygons are the convex hull encompassing the surveys that belong to the same cluster of the hierarchical ascending analysis this case lack of ecological interest. There is thus a compromise between orthogonality, that makes orthogonal representations and Euclidean distances relevant, and ecological interest when EOMs get meaningful spatial patterns.

MAFs were initially developed to eliminate noise in a set of images and extract their common signal. They are thus naturally ordered by decreasing auto-correlation at the reference distance, that is, by increasing order of the eigenvalue of the second PCA. Considering all MAFs in the analysis and interpretation allows to restore all the initial information and to describe perfectly the spatio-temporal variability of the observed data. However, there is no reason why the most important factors should be the most spatially regular. Indeed, the most important EOMs are rather those that endorse most of the variability. In our study case, the percentages of variance explained by the EOMs indicate that the two rankings are not similar, even though the first ten EOMs are the same.

Depending on the EOMs' algorithm, the approach refers either to the raw data or to their standardized version. The empirical Taylor's power law [START_REF] Taylor | Aggregation, variance and the mean[END_REF]) has been widely established in ecology. It can be summarized as the relationship between mean and variance (of count data). In this context, looking at spatial patterns rela- tively to the standard deviations amounts to study the spatial patterns relatively to the biomass. Therefore it is crucial to know precisely what kind of PCAs is used and which EOMs are generated. Analyses reported here concern the spatial patterns in relative terms, in order to compare and to group the surveys considering only the shapes of their spatial distributions. This also weights down the impact of the kriging performed prior to the analyses. An EOM decomposition together a the dimension reduction make it possible to summarize a large set of distribution maps onto a single factorial plan and to extract meaningful information. When the dispersion of a subgroup of maps is small, this opens the possibility to built a relevant mean distribution map which is of particular importance in ecology. In the present study, this allowed grouping surveys that are characteristic of each ecological season, and allowed building the expected spatial distribution for each ecological season. This means that the climatic season coincided with an ecological season. This corresponds to an a priori knowledge however with no strong empirical foundations. The present comparison and clusering of the spatial distributions strengthens the status of such a knowledge: it is no longer an a priori, but it is now deduced from a long series of spatial distributions.

Fig. 5

 5 Fig. 5 Dendrogram (Ward distance) of the 61 surveys based on the first ten EOMs. The color bar is defined by the season: in blue, the months corresponding to the cold season (light blue for the two extreme months of the cold season) and, in red, the months corresponding to the hot season (light red for the two extreme months of the hot season)

Fig. 6

 6 Fig. 6 Boxplots per cluster of the scores of the decompositions of the 61 input spatial distributions on the ten first EOMs. The color's transparency is proportional to the p-value of the difference between means per cluster: the more transparent, the less significant the difference between the two seasons

Acknowledgements Authors would like to give deep recognition to the people working at the Institut Mauritanien de Recherche Océanographique et des Pêches -IMROP in Mauritania that organized the surveys, collected the raw data and made the application possible.


Philippe Robutel

Adrien Leleu

Exoplanets Hans Trojan

Deeg

Hans Edition

Juan Antonio Deeg

Belmonte

Trojan Exoplanets

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

Introduction

The existence of Trojans in extrasolar systems has already been discussed in the chapter "Special Cases: Moons, Rings, Comets, and Trojans". The present chapter generalizes this concept by considering co-orbital planets. This kind of configuration consists of a planet pair trapped in the 1:1 mean motion resonance. Thus, con-Philippe Robutel IMCCE, Observatoire de Paris, PSL University, CNRS, Sorbonne Université, 77 avenue Denfert-Rochereau, 75014 Paris, France, e-mail: philippe.robutel@obspm.fr Adrien Leleu Observatoire de Genève, Université de Genève, Chemin Pegasi, 51, 1290 Versoix, Switzerland, e-mail: adrien.leleu@unige.ch trary to the Trojan case, discussed in the above-mentioned chapter, for which the mass is supposed to be small enough not to perturb the motion of the two other bodies (restricted case), it is possible to have co-orbital planets of comparable masses provided their sum remains small with respect to the mass of the star (less than 25 times smaller).

Co-orbital configurations are well known in the solar system. More than a century after the pioneer works of [START_REF] Euler | Considerationes de motu corporum coelestium[END_REF], who discovered the aligned configurations in the general 3-body problem, and [START_REF] Lagrange | OEuvres complètes[END_REF] who demonstrated the existence of the famous equilateral configurations, was discovered Achilles, the first Trojan asteroid [START_REF] Wolf | Photographische Aufnahmen von kleinen Planeten[END_REF]. Since then, thousands of small bodies are observed in the Jovian Trojan swarms L 4 and L 5 as well as some Neptune, Uranus, Mars and Earth Trojans (see https://www.minorplanetcenter.net/iau/lists/ Trojans.html). Saturn's satellite system also hosts co-orbital satellites, in particular the famous Janus and Epimetheus pair which, unlike most Trojans in the solar system, follow horseshoe trajectories of very large amplitude [START_REF] Dermott | The dynamics of tadpole and horseshoe orbits II. The coorbital satellites of saturn[END_REF][START_REF] Thomas | The inner small satellites of saturn: A variety of worlds[END_REF].

Although many theoretical works suggest that co-orbital exoplanets may also exist, such a pair of planets has yet to be observed. In particular, out of the hundreds of multi-planetary systems discovered by the Kepler mission, none were found in co-orbital configuration thus far. Although this result points to a relative rarity of the configuration, part of this absence might be due to observational biases.

After reviewing the stability issues of co-orbital systems, the main methods dedicated to their detection will be outlined.

Co-orbital configurations and their stability

The first stability result was established by [START_REF] Gascheau | Examen d'une classe d'équations différentielles et application à un cas particulier du problème des trois corps[END_REF] who proved that the Lagrange equilateral configuration was stable, in the coplanar circular case, providing the masses m j satisfy the following relation

(m 0 + m 1 + m 2 ) 2 ≥ 27(m 0 m 1 + m 0 m 2 + m 1 m 2 ).
In order to be satisfied, this condition requires that the mass of one body, for instance m 0 , dominates the others. Thus, it can be verified that, the mass ratio µ = (m 1 + m 2 )/m 0 must remain smaller than 0.04 to ensure stability. Therefore, it is not necessary for one of the masses to be negligible in front of the other two, as is the case with the Jovian Trojans. In particular, there is no reason why a Lagrange configuration consisting of one solar mass star and two Jovian planets should not be stable. These particular periodic solutions being linearly stable, their neighborhoods are filled with quasi-periodic orbits, known as tadpole-orbits, or also trojan-orbits by analogy with the Jovian Trojans , that correspond to small to moderate deformations of the equilateral configurations (Fig. 1, solid lines). Although the Lagrangian triangular equilibria are linearly stable, the overall stability of the tadpole region is not guaranteed. Indeed, for some particular values of the mass ratio µ associated with resonances, the stability domain shrinks, or even disappears entirely. (e.g. [START_REF] Danby | Stability of the triangular points in the elliptic restricted problem of three bodies[END_REF][START_REF] Giorgilli | Effective stability for a hamiltonian system near an elliptic equilibrium point, with an application to the restricted three body problem[END_REF][START_REF] Roberts | Linear stability of the elliptic Lagrangian triangle solutions in thethree-body problem[END_REF][START_REF] Nauenberg | Stability and Eccentricity for Two Planets in a 1:1 Resonance, and Their Possible Occurrence in Extrasolar Planetary Systems[END_REF][START_REF] Érdi | Secondary resonances of co-orbital motions[END_REF]. For a given tadpole orbit, the libration angle ζ , equal to the difference of the planetary mean-longitude λ 1 and λ 2 , oscillates around a value close to ±60 • (depending on whether one considers the L 4 or L 5 neighborhood) with a frequency ν equal to n √ 27µ/2 very close to Lagrange equilibria, n being the orbital frequency common to the two planets, and remains of the same order throughout the tadpole areas. In case of significant eccentricities, perihelia slowly precess with small secular frequency of order µ, which is therefore small compared to the libration frequency. Tadpole (solid lines) and horseshoe orbits (dashed lines) seen in a rotating frame with the common orbital frequency n of both planets. The coordinates of the central body are (0, 0) while while the co-orbital's orbits are plotted in (a j cos(λ jnt), a j sin(λ jnt)). The purple orbits correspond to the planet 1 whose mass satisfies m 1 /m 0 = 10 -4 while those of the second planet with m 2 /m 0 = 3 10 -5 are plotted in purple. The grey equilateral triangle indicates the Lagrange configuration. The initial conditions of the horseshoe orbit satisfy When µ is lower than a specific value close to 0.0004 [START_REF] Laughlin | Extrasolar Trojans: The Viability and Detectability of Planets in the 1:1 Resonance[END_REF], the libration angle ζ can undergo variations of much greater amplitude given birth to horseshoe orbits, which become stable below this critical value (Fig. 1, dashed lines). Along horseshoe orbits, the libration angle ζ evolves with an amplitude which is greater that 312 • , while for tadpole orbits this amplitude is bounded by 156 • . As with the tadpole orbits, the coplanar dynamic is driven by three distinct time scales: the fast one, related to the orbital period; the semi-fast, associated with the libration inside the 1:1 mean-motion resonance, and the secular related to orbital precession. For large eccentricities, the width of the tadpole and horseshoe domains shrink to make room for quasi-satellites orbits [START_REF] Giuppone | Dynamics of two planets in co-orbital motion[END_REF][START_REF] Pousse | On the co-orbital motion in the planar restricted three-body problem: the quasi-satellite motion revisited[END_REF]. Other families of orbits also exist for highly eccentric orbits [START_REF] Leleu | On the coplanar eccentric non-restricted co-orbital dynamics[END_REF]. As it is the case in the solar system, it can also be assumed that the orbits of the two co-orbital are not coplanar. If the mutual inclination is low, the main dynamical features are very similar to those described above. Higher inclinations allow transitions between theses different types of orbits and merging between them [START_REF] Namouni | Secular Interactions of Coorbiting Objects[END_REF], but in any case, the mutual inclination of stable co-orbital tadpoleorbits may not exceed about 60 • degrees (see [START_REF] Qi | Phase Structure of Co-orbital Motion with Jupiter[END_REF]. Other kinds of trajectories, which are not discussed here, also exist within the co-orbital resonance such as retrograde motions [START_REF] Morais | A numerical investigation of coorbital stability and libration in three dimensions[END_REF][START_REF] Sidorenko | Quasi-satellite orbits in the general context of dynamics in the 1:1 mean motion resonance: perturbative treatment[END_REF].

a 1 = a 2 = 1, e 1 = e 2 = 0, λ 1 - λ 2 = ϖ 1 -ϖ 2 = 15 • , while λ 1 -λ 2 = ϖ 1 -ϖ 2 = 45 • is set for tadpole.
When the co-orbital bodies are embedded in a multi-planetary system, the dynamics become much more complex and the size of the stability regions may decrease significantly. Because of the additional fundamental frequencies associated with the whole planetary system (three frequencies per planet) a wide range of resonances can destabilize the coorbital system [START_REF] Marzari | Stability of jupiter trojans investigated using frequency map analysis: the matros project[END_REF][START_REF] Robutel | The resonant structure of Jupiter's Trojan asteroids I: Long-term stability and diffusion[END_REF][START_REF] Leleu | Co-orbital exoplanets from close-period candidates: the TOI-178 case[END_REF]). These phenomena have many consequences. They are the cause of the slow erosion of the Jovian Trojan swarms [START_REF] Levison | The long-term dynamical of Jupiter's Trojan asteroids[END_REF][START_REF] Robutel | The resonant structure of Jupiter's Trojan asteroids I: Long-term stability and diffusion[END_REF], they introduce unstable areas in the tadpole and horseshoe regions [START_REF] Leleu | Co-orbital exoplanets from close-period candidates: the TOI-178 case[END_REF], and can go so far as to make the whole coorbital region unstable, as is the case with Saturn [START_REF] Holman | Dynamical stability in the outer solar system and the delivery of short period comets[END_REF][START_REF] Nesvorny | How long-live are the hypothetical Trojan populations of Saturn, Uranus, and Neptune?[END_REF][START_REF] Hou | Saturn Trojans: a dynamical point of view[END_REF].

Formation and evolution of co-orbital exoplanets

To understand the particularities of the co-orbital resonance, it is useful to recall the current understanding for the formation and evolution of other resonant configurations.

First and second-order MMRs are natural and common outcomes of models of formation of planetary systems: as planets form in the protoplanetary disc, they tend to migrate toward the star. Depending on their relative migration speed and eccentricities, they have a probability of being captured in each of the MMRs they cross (e.g. [START_REF] Lee | Dynamics and Origin of the 2:1 Orbital Resonances of the GJ 876 Planets[END_REF], the most commonly observed being the 3:2 MMR. Models hence predict that toward the end of the protoplanetary disc phase, numerous planetary systems are in close-in, compact resonant state (e.g. [START_REF] Cresswell | Three-dimensional simulations of multiple protoplanets embedded in a protostellar disc[END_REF][START_REF] Coleman | Pebbles versus planetesimals: the case of Trappist-1[END_REF][START_REF] Emsenhuber | The New Generation Planetary Population Synthesis (NGPPS). I. Bern global model of planet formation and evolution, model tests, and emerging planetary systems[END_REF]. As the protoplanetary disc dissipates, the eccentricity damping lessens, which can lead to instabilities (e.g. [START_REF] Terquem | Migration and the Formation of Systems of Hot Super-Earths and Neptunes[END_REF][START_REF] Pu | Spacing of Kepler Planets: Sculpting by Dynamical Instability[END_REF][START_REF] Izidoro | Breaking the chains: hot super-Earth systems from migration and disruption of compact resonant chains[END_REF]. For planets that remained in resonance, and are close enough to their host star, tides become the dominant force that affect the architecture of the systems, which can then lead to a departure of the period ratio from the exact resonance [START_REF] Henrard | A second fundamental model for resonance[END_REF][START_REF] Papaloizou | On the dynamics of multiple systems of hot super-Earths and Neptunes: tidal circularization, resonance and the HD 40307 system[END_REF][START_REF] Delisle | Dissipation in planar resonant planetary systems[END_REF]. This succession of dissipative processes partly explain why the observed close-in multi-planetary systems are often just outside MMRs, but rarely inside [START_REF] Fabrycky | Architecture of Kepler's Multi-transiting Systems. II. New Investigations with Twice as Many Candidates[END_REF].

Both of the aforementioned phenomenon, i.e. the smooth convergent migration in resonance during the protoplanetary disc, and the slow departure from resonance due to tides, cannot be applied to the co-orbital resonance. Indeed, this resonance is surrounded by a chaotic area due to the overlap of first-order MMRs [START_REF] Wisdom | The resonance overlap criterion and the onset of stochastic behavior in the restricted three-body problem[END_REF][START_REF] Deck | First-order Resonance Overlap and the Stability of Close Two-planet Systems[END_REF][START_REF] Petit | Hill stability in the AMD framework[END_REF], and a slow crossing of this area would generally result in the excitation of the bodies' eccentricities, leading to collision or scattering instead of the capture into the resonance.

Nevertheless, several formation scenarios of co-orbital configurations can be found in the literature. [START_REF] Laughlin | Extrasolar Trojans: The Viability and Detectability of Planets in the 1:1 Resonance[END_REF] suggested planet-planet gravitational scattering followed by strong eccentricity damping due to the protoplanetary disc, or accretion in situ at the L 4 /L 5 points of a primary. Alternatively, the presence of an external body such as a massive planet can help the capture of coorbital bodies through resonant crossing (see [START_REF] Nesvorny | Capture of Trojans by Jumping Jupiter[END_REF], for the case of the solar system ).

In the in situ scenario, different studies yielded different upper limit to the mass that can form at the L 4 /L 5 equilibrium point of a giant planet: [START_REF] Beaugé | Co-orbital terrestrial planets in exoplanetary systems: a formation scenario[END_REF] obtained a maximum mass of ∼ 0.6M ⊕ , while [START_REF] Lyra | Standing on the shoulders of giants. Trojan Earths and vortex trapping in low mass self-gravitating protoplanetary disks of gas and solids[END_REF] obtained 5 -15M ⊕ planets in the tadpole area of a Jupiter-mass primary.

Co-orbital pairs form naturally in models of formation of planetary systems (e.g. [START_REF] Cresswell | On the growth and stability of Trojan planets[END_REF][START_REF] Coleman | Pebbles versus planetesimals: the case of Trappist-1[END_REF][START_REF] Emsenhuber | The New Generation Planetary Population Synthesis (NGPPS). I. Bern global model of planet formation and evolution, model tests, and emerging planetary systems[END_REF]) in typically a few percent of the systems, through scattering, but also during the migration and encounter of resonant chains: in the simulation of the formation of Trappist-1 by [START_REF] Coleman | Pebbles versus planetesimals: the case of Trappist-1[END_REF], co-orbitals were mostly found as a part of a resonant chain, which appears to have a stabilising effect (Leleu et al. 2019a). However, models of the formation of planetary systems model the disc-planet interactions using analytical expression computed for single planets. The presence of a co-orbital configuration has been shown to impact significantly the local proto-planetary disc structure, sometime changing the stability of the co-orbital pair [START_REF] Pierens | Disruption of co-orbital (1:1) planetary resonances during gas-driven orbital migration[END_REF]Leleu et al. 2019a).

Nonetheless, co-orbitals can also form in hydro-dynamical simulations of the proto-planetary disc with two or more planets embedded (e.g. [START_REF] Crida | Minimum Mass Solar Nebulae and Planetary Migration[END_REF][START_REF] Penzlin | 1:1 orbital resonance of circumbinary planets[END_REF]. In several simulations of the formation of the outer part of the solar system, [START_REF] Crida | Minimum Mass Solar Nebulae and Planetary Migration[END_REF] found that Uranus and Neptune ended up in co-orbital configuration, while being both trapped in the same MMR with Saturn. This reinforces the idea that resonant chains might be favorable to the formation of co-orbital configuration. As long as the planets do not fully open a gap in the proto-planetary disc (typically for planet smaller than Neptune), co-orbitals evolving in the protoplanetary disc appear to be more stable when the leading planet is more massive [START_REF] Pierens | Disruption of co-orbital (1:1) planetary resonances during gas-driven orbital migration[END_REF]Leleu et al. 2019a), while trailing more massive planet tends to disrupt the configuration. In the case of a terrestrial planet co-orbital with a Jupiter-sized planet, the stability seemed to mainly depend on the proto planetary disc parameters.

For close-in systems, the effect of tides raised by the star onto co-orbital planets appears to always be disruptive [START_REF] Rodríguez | Tidal evolution of close-in exoplanets in co-orbital configurations[END_REF][START_REF] Couturier | An analytical model for tidal evolution in co-orbital systems[END_REF][START_REF] Dobrovolskis | Do tides destabilize Trojan exoplanets?[END_REF]. However this phenomenon only impacts planets close to their host star, and most co-orbitals that are at more than 10 days of orbital period at the disc dispersal should survive for billions of years [START_REF] Couturier | An analytical model for tidal evolution in co-orbital systems[END_REF]).

Detection of co-orbital planets

Stacking the phase-folded lightcurve of thousands of KOIs, [START_REF] Hippke | Photometrys Bright Future: Detecting Solar System Analogs with Future Space Telescopes[END_REF] found statistical hints of dips near the L 4 and L 5 equilibria of known planets. However, despite the discovery of thousands of exoplanets in hundreds of multiplanetary systems, no extra-solar co-orbital configurations has been confirmed thus far. The subsequent sections discuss the signature of co-orbital configurations in photometric (transit) and spectroscopic (radial velocity) measurements.

Transit

Transit surveys such as Kepler, TESS and the upcoming PLATO mission have discovered the bulk of known exo-planetary systems. As for any other pair of planets, the detection of co-orbitals by transit would require both planets to orbit roughly in the same orbital plane. In addition, their detection requires that each of the co-orbital is large enough for its individual transits to be detected.

For smaller planets, the data-reduction pipelines of transit surveys use an approach similar to the box least square algorithm (BLS, see [START_REF] Kovács | A box-fitting algorithm in the search for periodic transits[END_REF][START_REF] Jenkins | Overview of the Kepler Science Processing Pipeline[END_REF][START_REF] Jenkins | The TESS science processing operations center[END_REF]. These algorithms fold each light curve over a large number of different periods and look for transits in the folded data.

This method is adapted to planets on unperturbed or moderately perturbed orbits for which the transits occur periodically. However, gravitational interactions between planets in resonant systems produce transit timing variations (TTVs). Figure 3 shows examples of TTVs for different co-orbital configurations. These TTVs can prevent the detection of small planets if their amplitude σ T TV (amplitude of the oscillation of the curves along the x-axis) is larger than the transit duration D of the planet (width of the curves), and their period is shorter than the mission duration (y-axis). In this case, the signal-to-noise ratio of the transit of the planet recovered by BLS-like algorithms decreases by a factor 1/(1 + σ T TV /D) [START_REF] Leleu | Alleviating the transit timing variation bias in transit surveys. I. RIVERS: Method and detection of a pair of resonant super-Earths around Kepler-1705[END_REF]). In addition, the recovered transit shape is altered by the TTVs, and therefore the planet can be missed, or mistaken for a false positive [START_REF] Leleu | II. RIVERS: Twin resonant Earth-sized planets around Kepler-1972 recovered from Kepler's false positive[END_REF].

Assuming that the masses m p of the two planets are comparable and that it is the same for their eccentricities e p , a rough calculation (see [START_REF] Gallardo | Semianalytical model for planetary resonances. Application to planets around single and binary stars[END_REF][START_REF] Nesvorný | Dynamics and Transit Variations of Resonant Exoplanets[END_REF], and references therein for accurate expressions) shows that, for a MMR of order q, the TTVs occur on a timescale proportional to:

P T TV ∝ (m p /m ⋆ ) -1/2 e -q/2 p P orb .
(1)

Note that this formula is not singular because the eccentricities cannot be cancelled simultaneously if q ≥ 1. With q = 0, co-orbitals hence generate faster TTVs than other MMRs. In addition their amplitude can be a significant fraction of the orbital period of the planets (e.g. [START_REF] Nesvorný | The Effect of Conjunctions on the Transit Timing Variations of Exoplanets[END_REF]Leleu et al. 2019b, and Fig. 3). Co-orbitals are hence especially subject to the TTV bias, and dedicated methods are required to recover potential pairs of small co-orbitals in transit surveys data.

Assuming a large mass discrepancy between co-orbitals, Janson (2013) looked for dips near the L 4 and L 5 equilibrium of known planets, using an analytical model for small libration of the resonant angle and low eccentricities. No significant signal was found, despite a sensitivity level down to Earth-size objects. However, such a signal could still exist for larger amplitudes of libration or eccentricities.

If the co-orbitals are similar in size, both of them could be affected by the TTV bias, in which case both would have been missed by BLS-like approaches, and all orbital periods have to be probed to recover them. The QATS algorithm is adapted to the research of planets with TTVs: it allows the timing of each transit to vary at the same time that it adjusts the best transit depth and duration [START_REF] Carter | The Quasiperiodic Automated Transit Search Algorithm[END_REF]. Alternatively, neural networks trained to recover the track of planets in river diagrams were shown to be capable of detecting planets with large TTVs and low S/N [START_REF] Leleu | Alleviating the transit timing variation bias in transit surveys. I. RIVERS: Method and detection of a pair of resonant super-Earths around Kepler-1705[END_REF][START_REF] Leleu | II. RIVERS: Twin resonant Earth-sized planets around Kepler-1972 recovered from Kepler's false positive[END_REF].

Co-orbitals might also be detected even if only one of them transit: with a good phase coverage and high S/N transits, it is possible, based on the TTVs of a single planet, to determine in which MMR that planet and its perturber are (see the case of KOI-142, Nesvorný et al. 2013). When co-orbitals librate with a significant amplitude the peculiar shape of their TTVs (arches for tadpole-orbits as shows in Fig. 3

Radial velocity

Co-orbitals librating in Trojan or horseshoe configurations modulate the radial velocity (RV) signal with a period equal to their libration period. The amplitude of this modulation increases when the mass of the planets tend toward equal value and the amplitude of libration increases [START_REF] Laughlin | Extrasolar Trojans: The Viability and Detectability of Planets in the 1:1 Resonance[END_REF][START_REF] Ford | Using Transit Timing Observations to Search for Trojans of Transiting Extrasolar Planets[END_REF]. If this modulation is observed, the estimation of the Keplerian signal and the first modulation harmonics is sufficient to recover the orbit and masses of the planets [START_REF] Leleu | Detectability of quasi-circular co-orbital planets. Application to the radial velocity technique[END_REF].

For a pair of co-orbital at the exact Lagrangian equilibrium, the modulation vanishes and the generated radial velocities are identical from those induced by a single planet located at their longitudinal barycenter with a mass of m 2 1 + m 2 2 + m 1 m 2 , (see [START_REF] Giuppone | Origin and detectability of coorbital planets from radial velocity data[END_REF][START_REF] Leleu | Detection of co-orbital planets by combining transit and radial-velocity measurements[END_REF], and hence the co-orbital configuration cannot be identified. However, in the best configuration (equal mass planets, large libration amplitude), the maximum modulation amplitude can be up to ∼ 30% for tadpole orbit and 100% for horseshoe (the RV signal vanishes when the planets are in opposition with respect to the star). In this case, the modulation would be detectable with current radial-velocity instruments. It however requires a long observational baseline to cover the libration period, and might be difficult to disentangle from other dynamical effects or noise. In the case of quasi-satellite orbits, [START_REF] Laughlin | Extrasolar Trojans: The Viability and Detectability of Planets in the 1:1 Resonance[END_REF] remark that the RV signal is also similar to one induced by a single planet, however a small fraction of the orbital phase shows a strong variation in RV.

The co-orbital signature is similar in astrometry. However, the system observed in astrometry tend to be at larger orbital period, making it unlikely to observe a system long enough to see the resonant modulation of the signal.

Radial velocity and transit

Combining RVs and transits allows observers to detect co-orbitals when only one of them transits and without the need to observe the resonant libration. For a single planet in circular orbit, the time of mid-transit coincides with the instant where the radial-velocity reaches its mean value. If the transiting planet has a co-orbital companion, there is a time shift ∆t between these two phenomena that depends on the mass distribution between the co-orbitals and their angular separation ζ [START_REF] Ford | Observational Constraints on Trojans of Transiting Extrasolar Planets[END_REF][START_REF] Leleu | Detection of co-orbital planets by combining transit and radial-velocity measurements[END_REF]). If P is the common orbital period, this time shift reads:

∆t = P 2π ζ 2 -arctan (1 -2δ ) tan ζ 2 with δ = m 2 m 1 + m 2 , (2)
m 1 being the transiting planet mass and m 2 the co-orbital one. But this time shift also exits with a single planet on an elliptic orbit of eccentricity e and argument of pericentre ω. In this case, the time shift reads: ∆t = -e cos ω P 2π .

(3)

The sensibility of this approach hence depends on the precision reached on the quantity e cos ω of the transiting planet. So far, this approach could only be used on hot Jupiter, for which the observation of the secondary eclipse can constrain e cos ω [START_REF] Madhusudhan | Empirical Constraints on Trojan Companions and Orbital Eccentricities in 25 Transiting Exoplanetary Systems[END_REF]Lillo-Box et al. 2018a,b). No co-orbital companion were found in these studies, which is consistent with their instability to tides [START_REF] Couturier | An analytical model for tidal evolution in co-orbital systems[END_REF].

Conclusion

Co-orbital exoplanets can be stable for duration longer than the life of their host star if they are sufficiently far away from the star to prevent large tidal disruption, and in the absence of large perturbations by additional planets. In addition, models of formation of planetary systems tend to create pairs of co-orbital exoplanets. However this resonance has several particularities, like that of being surrounded by large chaotic zones, or to modify the local proto-planetary disc structure because of the presence of two planets at the same distance to the star. This require extra care when studying to formation and evolution of co-orbital systems.

Co-orbital exoplanets have not be found yet. However, the methods used to detect the bulk of the exoplanets have biases that impede their capability to recover exoplanets in this kind of configuration. The development of dedicated data analysis methods is required, and the continuous increase in measurement precision and observation baseline could allow us to identify a co-orbital signature in the future. Co-orbitals could also be found by other methods such as direct imaging which, combined with stability studies, has been shown to be able to identify the resonant state of multi-planetary systems [START_REF] Goździewski | An exact, generalized laplace resonance in the hr 8799 planetary system[END_REF].

  Fig.1Tadpole (solid lines) and horseshoe orbits (dashed lines) seen in a rotating frame with the common orbital frequency n of both planets. The coordinates of the central body are (0, 0) while while the co-orbital's orbits are plotted in (a j cos(λ jnt), a j sin(λ jnt)). The purple orbits correspond to the planet 1 whose mass satisfies m 1 /m 0 = 10 -4 while those of the second planet with m 2 /m 0 = 3 10 -5 are plotted in purple. The grey equilateral triangle indicates the Lagrange configuration. The initial conditions of the horseshoe orbit satisfy a 1 = a 2 = 1, e 1 = e 2 = 0, λ 1λ 2 = ϖ 1ϖ 2 = 15 • , while λ 1λ 2 = ϖ 1ϖ 2 = 45 • is set for tadpole.

Fig. 2

 2 Fig.2Quasi-satellite orbit. Trajectory of the first planet (purple) and of the second one (black) seen in a rotating frame with the common orbital frequency n of the both planets. The masses are the same than in Fig.1, while the initial elliptic elements satisfy a 1 = a 2 = 1, e 1 = 0.05, e 2 = 0.6.

Fig. 3

 3 Fig.3River diagrams[START_REF] Carter | Kepler-36: A Pair of Planets with Neighboring Orbits and Dissimilar Densities[END_REF]) for co-orbitals in tadpole orbit (a.), horseshoe orbit (b.) and quasi-satellite (c.). These diagrams show the transit timings of each planet as function of the phase of the average period of the planets (here 10 days, x axis), for each epoch (y axis). In each case, the yellow track shows the transit of the most massive of the two planets, and the blue track shows the transits of the least massive. Assuming a sun-like star, m 1 = 2m 2 = 10m ⊕ for a. and c.; and m 1 = 2m 2 = 2m ⊕ for b.

  -a and triangle for horseshoe as shown in Fig. 3-b and in the figure 8 by Vokrouhlický and Nesvorný 2014) could suffice to confirm the existence of a co-orbital configuration.