G Parent

P Dular

F Piriou

A Abakar

Accurate Projection Method of Source Quantities in Coupled Finite Element Problems

Keywords: Coupled Problems, Finite Element Method, Projection Method

Projecting discrete quantities from one mesh to another one is a critical operation which is usually required in multi-physical problems. The efficiency of a Galerkin projection method is studied for magneto-thermal and magneto-elasticity couplings. As examples of application, two electromagnetic devices are studied: a three-phase busbar system as a magneto-thermal problem and an electromagnet as a magneto-mechanical problem.

I. INTRODUCTION

In many studies, such as multi-physical coupled problems, one can need to use different meshes to take into account some regions or not according to their influence on the computed field. In addition, each problem can ask for a particular fineness of mesh. Hence, the main issue about coupled problems using dedicated meshes is to properly transfer each discrete field obtained on one mesh onto another one.

The finite element (FE) method allows evaluating a field everywhere in a studied domain via the direct interpolation of basis functions. This method can be used to perform mesh to mesh interpolation. Another approach for mesh to mesh interpolation consists in defining the discrete expression of a field on the two considered meshes and then minimize the gap between them through a weighted residual method to build an easy to solve linear matrix system. In a previous study [START_REF] Parent | Using a Galerkin Projection Method for Coupled Problems[END_REF], the theory of the method has been presented for the projection of magnetic, thermal and mechanical quantities. In this paper, this theory is recalled and the method is used to perform the transfer of source quantities in the cases of magneto-thermal and magneto-mechanical problems using dedicated meshes for each physical phenomenon. As examples of applications, two electromagnetic devices are studied: a three-phase busbar system as a magneto-thermal problem and an electro-magnet as a magneto-mechanical problem. The influence of the projected source quantities on the temperature -in the first case -and the displacement -in the second case -is studied.

II. FIELDS COMPUTATION

A. Magnetic Field Computation

Maxwell's equations and material constitutive laws are written in a studied domain Ω:

curl , curl t = = -∂ h j e b (1)
,

σ µ = = j e b h ( 2
)
where h is the magnetic field, b is the magnetic flux density, e is the electric field, j is the current density, µ is the magnetic permeability and σ is the electric conductivity. Boundary conditions are defined on two complementary parts Γ b and Γ h of the boundary of the studied domain:

0 on and 0 on b h ⋅ = Γ × = Γ b n h n (3)
with n the unit normal exterior to Ω. By combining equations ( 1) with ( 2) one can obtain either the magnetodynamic h-formulation [START_REF] Bossavit | Force-related Nuts and Bolts in the Discretization Toolkit for Electromagnetic[END_REF] or a-v formulation [START_REF] Bossavit | Edge-Element Computation of the Force Field in Deformable Bodies[END_REF].

1 curl( curl ) ( ) 0 t - + ∂ = h h σ µ (4) 1 curl( curl )
( grad ) 0 t v - + ∂ + = a a µ σ (5)
where a is the magnetic vector potential and v is the scalar electric potential, with b = curl a and e = -∂ t agrad v.

B. Thermal Field Computation

The heat transfer in a studied domain is described by the Fourier's law:

( )

div grad p t C T T p ∂ - = ρ λ (6)
where ρ is the mass density, C p is the specific heat and λ is the thermal conductivity. In addition to (6) boundary conditionsideal insulation, isotherm wall, radiation and convection -also have to be taken into account.

C. Elastic Field Computation

In the case of small displacements u, the equilibrium equation is given by:

0 B ∇ + = σ F (7)
where F B represents body forces and σ σ σ σ is the stress. The latter is linked to the strain ε ε ε ε via Hook's law = ⋅ σ E ε where E is Young modulus.

III. COUPLING STRATEGIES AND QUANTITIES

For coupled problems, two different coupling strategies, either direct (strong) or indirect (weak), can be used. The most natural and obvious approach from a mathematical point of view is the direct solution since it allows taking into account every single field interaction between each other. Nevertheless, such a strategy can be computationally expensive and is only justified in particular configurations.

On the other hand, under certain conditions, one can use indirect coupling strategy. For magneto-thermal coupled problems, using such an approach is justified by the fact that each field has its own specific time scale. For magnetomechanical coupled problems, indirect coupling can be used in the case of small deformations.

With such an approach, the interaction between the different fields is taken into account through coupling quantities, which are the heat sources and the temperature (which influences the electrical conductivity) in a magneto-thermal study and the magnetic force in a magneto-mechanical study.

In our case, using dedicated grids for each physical phenomenon requires projecting those coupling quantities onto different meshes (Fig. 1).

IV. PROJECTION METHOD

A. Galerkin Projection Method General Formulation

One considers a studied domain Ω and two function spaces

( ) 2 U L ∈ Ω and
( )

2

V L ∈ Ω with a known field u U ∈ . The aim of the method is to find v V ∈ by using a weighted residual method, with test functions ψ, so that:

( ) 0 , v u d V Ω - = ∀ ∈ ∫ ψ τ ψ (8)
In the discrete domain, a field can be interpolated via the basis functions defined on the mesh used to discretize Ω. Hence, if the field u is discretized on a source mesh and v on a target mesh, writing (8) in the discrete domain leads to: ,

j j i i t s v u t t s s j i v w d u w d V ∈ ∈ Ω Ω = ∀ ∈ ∑ ∑ ∫ ∫ ε ε ψ τ ψ τ ψ (9)
where w are the basis functions associated to the entity i of the source mesh and j of the target mesh respectively. Hence, the Galerkin method is applied with each basis function j t w as a test function ψ so (9) leads to a linear matrix system where the discrete coefficients

j t v
are the unknown entities. The numerical integration has to be done using an appropriate number of Gauss points to accurately sample to source field variations.

B. Projection of Magneto-Thermal Coupling Quantities

As we mentioned it earlier, coupling quantities in a magneto-thermal problem are heat sources and temperature. According to Tonti diagram (Fig. 2) those quantities are discretized in the volume FE space 3 E and the nodal FE space 0 E respectively [2].

g q p t T ∂ t ∂ T p C ρ λ - div grad 0 E 1 E 2 E 3 E Fig. 2. Tonti diagram applied to thermal problem
Hence, the projection of the heat source and of the temperature is performed by replacing j t w , and so the test function ψ, by a volume basis function and a nodal basis function respectively. However, the projection of the heat sources can be tricky. Indeed, the fact that it is discretized in the volume FE space means that the heat sources are constant per FE. Hence, in some cases, such as devices with skin effects for example, a step effect can appear on the projected quantity. Moreover, most of thermal codes historically use nodal element space to discretize the heat sources. Hence, for both reasons, one can use nodal basis functions for both the projected heat source and its associated test function ψ in the projection formulation.

C. Projection of the Force

The proposed projection method (9) is made for the projection of discrete quantities in the sense of FE. That means that one has to be careful about the space in which quantities to be projected take place, especially about the force since there are various methods to determine it. For example, in the case of Laplace's force × j b , the right function to use to replace the test function ψ is nodal basis function [START_REF] Bossavit | Force-related Nuts and Bolts in the Discretization Toolkit for Electromagnetic[END_REF].

In the case of nodal forces given by the virtual work method, the projection cannot be done directly by using (9). Indeed, those forces do not define a field [START_REF] Bossavit | Edge-Element Computation of the Force Field in Deformable Bodies[END_REF], but rather a set of discrete force components acting on each node. Hence, once the nodal forces F n given by the virtual work method are computed for each node n of the mesh supporting the magnetic computation, one needs to determine the vectorial values f n of the corresponding force field. Those f n then define a force density and are, by definition, related to F n via the expression ,

s i i s s n s s n s s i N w w d n N τ ∈ Ω = ∀ ∈ ∑ ∫ f F (10)
with N s the set of nodes belonging to the source mesh. Hence, once the force density is evaluated on the source mesh (supporting the magnetic computation), it can be projected onto the target one (supporting the mechanical computation) by replacing the test function ψ by a nodal basis function.

V. APPLICATIONS

A. Magneto-Thermal Coupled Problem

A three-phase busbar system is taken as an example for the application of the coupling strategy presented in Fig. 1 (left). It is composed of five copper bars included in a concrete bloc (Fig. 3). Dimensions of each bar and of the concrete bloc are 100 8 mm² × and 104 200 mm² × respectively.

Fig. 3. Three-phase busbar system Three different meshes are generated for the study. M magsupporting the electromagnetic computation -is fine enough to take into account the skin effect that appears in the conducting parts with good accuracy, whereas M ther1 and M ther2 are coarser. The number of elements of these three meshes are 28791, 3763 and 1144, respectively.

The study is performed in two steps. M mag is first used to perform a magneto-thermal computation via the formulations mentioned in section II. Due to the high number of elementsand so high number of degrees of freedom -belonging to M mag , the so-computed FE solution is considered as a reference one for both magnetic and thermal problems. On a second step, M mag is used for the magnetodynamic computation only. The Joule losses are then projected successively onto M ther1 and M ther2 using (3). For that computation, both nodal and volume basis functions are used. Results obtained for that operation are compared to the reference solution in Fig. 4. Moreover, total Joule losses evaluated on target meshes are compared in table I.

The good accuracy between the projected Joule losses evaluated thermal meshes and the reference one is pointed out in Fig. 4 andtable I. Moreover, even if local error due to "step effect" appear when using volume basis functions for the projection, the global Joule losses do not suffer from any major loss of accuracy, which is an advantage of the proposed method in comparison with the classical interpolation technique.

Then, the so-computed heat sources are used to perform a thermal computation via (6) on M ther1 and M ther2 . Results of that computation are shown in Fig. 5. The results obtained with projected heat sources are in very good agreement with those obtained by using one single mesh.

B. Magneto-Elastic Coupled Problem

An electro-magnet above a plate is taken as an example for the application of the coupling strategy presented in Fig. 1 (right). For the mechanical computation, the left side of the plate is clamped and the studied quantity is the displacement of the right side of the plate. Since the mechanical studied domain is the plate, we generate three different meshes -M mag , M mech1 and M mech2differing from each other on the number of elements inside that region, which are 2000, 1200 and 800 respectively. As in the previous section, the study is performed in two steps. A magneto-mechanical computation is first performed on M mag to get a reference solution in term of force densities and displacement at the free side of the plate. Hence, in a second step, the electromagnetic and mechanical problems are solved on dedicated meshes, requiring the projection of the force densities from M mag onto M mech1 and M mech2 . As an illustration of the fact that the forces given by the virtual work method are strongly linked to the mesh on which they are computed, those latter are shown in Fig. 7 for both M mag (left) and M mech1 (right). The corresponding force density distributions are shown in Fig. 8. Results of the force density projection operation are shown in Fig. 9. A good accuracy between the force density computed on M mag and the projected ones is obtained. Then, the projected force densities are used for the magneto-elastic computation. The displacement of the free side of the plate for both computations -using M mech1 and M mech2 -are presented in table II. The mechanical computation using dedicated meshes -and so projected source terms -gives very accurate results in comparison with the reference solution.

VI. CONCLUSION

A Galerkin projection method has been used for the projection of source quantities in both magneto-thermal and magneto-mechanical problems using dedicated meshes for each physical field computation. It has been shown that the projected quantities are in good agreement with the original ones and can be used as source terms for multi-physical problems using dedicated meshes.

Fig. 1 .

 1 Fig. 1. Magneto-thermal (left) and magneto-mechanical (right) indirect coupled problems using dedicated meshes

ε

  are sets of geometrical entities belonging to the source and the target mesh respectively. Moreover,

Fig. 4

 4 Fig. 4 Joule losses densities versus y-position in a bar; projected quantities with volume (top) or nodal (bottom) basis functions

Fig. 5 .

 5 Fig. 5. Evolution of the temperature versus x-position on top of bars; projected Joule losses with volume (top) or nodal (bottom) basis functions

Fig. 6 .

 6 Fig. 6. Electro-magnet above a plate

Fig. 7 .

 7 Fig. 7. Force given by virtual work method under a tooth using Mmag (left) and Mmech1 (right)

Fig. 9 .

 9 Fig. 9. Force density versus x-position along the top of the plate, under a tooth

TABLE I VALUE

 I

	MESH USED FOR THERMAL COMPUTATION	OF THE TOTAL LOSSES IN THE DEVICE AFTER PROJECTION	VALUE OF THE RELATIVE ERROR COMPARED WITH THE REFERENCE SOLUTION (%)
		v w	n w	v w	n w
	Mther1	347.76	347.76	0.028	0.028
	Mther2	346.98	346.98	0.252	0.252

ACKNOWLEDGMENT This work was supported in part by the French CNRT project (Centre National de Recherche Technologique) and by the Belgian Science Policy (IAP P6/21).


Hélène Morlon
email: helene.morlon@bio.ens.psl.eu

Stéphane Robin

Florian Hartig

Studying speciation and extinction dynamics from phylogenies: addressing identifiability issues

Keywords: Speciation, extinction, diversification, phylogenies, parameter identifiability, model congruency

A lot of what we know about past speciation and extinction dynamics is based on statistically fitting birth-death processes to phylogenies of extant species. Despite their wide use, the reliability of these tools is regularly questioned. It was recently demonstrated that vast 'congruent' sets of alternative diversification histories cannot be distinguished (i.e. are not identifiable) using extant phylogenies alone, reanimating the debate about the limits of phylogenetic diversification analysis. Here, we summarize what we know about the identifiability of the birth-death process and how identifiability issues can be addressed. We

conclude that extant phylogenies, when combined with appropriate prior hypotheses and regularization techniques, can still tell us a lot about past diversification dynamics.

Glossary

Asymptotic (or theoretical) unidentifiability: situation when there are distinct combinations of the model parameters that cannot be told apart even in the limit of an infinite number of

Molecular phylogenies and diversification dynamics

The diversity of life on Earth has arisen from a succession of speciation and extinction events (see Glossary). The rates at which ancestral species give rise to new daughter species (the speciation rate, 𝜆) or go extinct (the extinction rate, 𝜇) reflect underlying ecological and evolutionary processes, and shape species richness over geological timescales. Understanding how these rates have changed through time has long been of interest to evolutionary biologists [START_REF] Van Valen | A new evolutionary law[END_REF][START_REF] Gould | The shape of evolution: a comparison of real and random clades[END_REF][START_REF] Raup | Mathematical models of cladogenesis[END_REF][START_REF] Alroy | Dynamics of origination and extinction in the marine fossil record[END_REF][START_REF] Silvestro | Bayesian Estimation of Speciation and Extinction from Incomplete Fossil Occurrence Data[END_REF][START_REF] Simpson | Tempo and mode in Evolution[END_REF][START_REF] Stanley | Macroevolution: Pattern and Process[END_REF][START_REF] Quental | How the Red Queen Drives Terrestrial Mammals to Extinction[END_REF]. While the first estimates of speciation and extinction rates were derived from the fossil record, researchers now also widely use dated phylogenies of present-day species (so-called 'reconstructed (or extant) phylogenies', thereafter referred to as 'phylogenies' for simplicity) to study past speciation and extinction dynamics [START_REF] Ricklefs | Estimating diversification rates from phylogenetic information[END_REF][10][START_REF] Stadler | Recovering speciation and extinction dynamics based on phylogenies[END_REF][START_REF] Pennell | An integrative view of phylogenetic comparative methods: connections to population genetics, community ecology, and paleobiology[END_REF].

Nee et al. [START_REF] Nee | The reconstructed evolutionary process[END_REF] showed, using the homogeneous birth-death (BD) process, that despite extinct species being absent from a phylogeny, extinctions leave a distinctive signal in the timing of branching patterns, known as the 'pull of the present'. Under the assumption of homogeneous and constant speciation and extinction rates, it is therefore possible to estimate these rates from phylogenies. A wide range of more complex models grounded on the homogeneous BD process have now been developed, and used to test hypotheses about past diversification dynamics [START_REF] Rabosky | Explosive Evolutionary Radiations: Decreasing Speciation or Increasing Extinction Through Time?[END_REF][15][START_REF] Stadler | Mammalian phylogeny reveals recent diversification rate shifts[END_REF][START_REF] Etienne | Diversity-dependence brings molecular phylogenies closer to agreement with the fossil record[END_REF][18][START_REF] May | A Bayesian approach for detecting the impact of massextinction events on molecular phylogenies when rates of lineage diversification may vary[END_REF][START_REF] Valente | A simple dynamic model explains the diversity of island birds worldwide[END_REF][21][START_REF] Jetz | The global diversity of birds in space and time[END_REF].

Increasing flexibility of the models brings new issues, however, such as parameters that may not be identifiable. Here, we discuss the identifiability of speciation and extinction rates in a variety of homogeneous BD models, and clarify the theoretical limits that non identifiability imposes on phylogenetic diversification analysis. We conclude that although speciation and extinction histories are statistically unidentifiable if the underlying functions are completely unconstrained [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF], this does not imply that phylogenies can't reveal speciation and extinction dynamics [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF][START_REF] Pagel | Evolutionary trees can't reveal speciation and extinction rates[END_REF]. We hold that in most practical scenarios, a priori hypotheses, biological knowledge or statistical regularization can make the problem identifiable.

Identifiability of speciation and extinction rates

To clarify the issue of identifiability, it is useful to make a distinction between asymptotic (or fundamental) and practical unidentifiability. Asymptotic unidentifiability corresponds to the case when distinct parameter combinations cannot be told apart, even in the limit of an infinite number of observations; practical unidentifiability corresponds to the case when parameters cannot be told apart from the limited number of observations available in practice.

Asymptotic identifiability of the homogeneous BD model

Nee et al. [START_REF] Nee | The reconstructed evolutionary process[END_REF] showed that the homogeneous constant rate BD model with complete sampling (i.e. all present-day species are represented in the phylogeny) is asymptotically identifiable.

Incomplete sampling can be accounted for by assuming that each extant species is sampled with the same probability 𝜌 (𝜌 < 1), but already in this simple extension of the model, if 𝜌 is a parameter to be estimated, 𝜆, 𝜇 and 𝜌 are not asymptotically identifiable [START_REF] Stadler | On incomplete sampling under birth-death models and connections to the sampling-based coalescent[END_REF]. To solve this identifiability problem, the fraction of present-day species represented in the phylogeny is often included as prior information on 𝜌, which renders 𝜆 and 𝜇 asymptotically identifiable.

Extending the work of Nee, Stadler [START_REF] Stadler | Mammalian phylogeny reveals recent diversification rate shifts[END_REF] showed that the « episodic » birth-death model (EBD, also called birth-death-shift, BDS), where diversification rates are piecewise constant (i.e. constant on successive time intervals, or epochs) is asymptotically identifiable. More recently, Legried & Terhorst [26] confirmed this result and showed that it holds even if the epochs are not fixed. However, the BDS model with mass extinction events, i.e. including the possibility that sudden (simultaneous) extinction events can occur at the end of each epoch (equivalent to sampling each species with an epoch-specific probability 𝜌), is not identifiable [START_REF] Stadler | Mammalian phylogeny reveals recent diversification rate shifts[END_REF].

In the case when 𝜆(𝑡) (or 𝜇(𝑡)) are smooth functions of time and are not constrained to follow specific functional forms such as the exponential or any other biologically-motivated function, Louca & Pennell [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF] showed that there is an infinity of 'congruent' functions that yield the same likelihood, meaning that this process is not asymptotically identifiable (Box 1).

Practical identifiability of the homogeneous BD model

When applying birth-death models to real data, a further issue arises: the size of phylogenies is typically not huge. Finite data sizes impose limits to the identifiability of any given model, as the confidence in the parameter estimates decreases with decreasing sample sizes. This is well illustrated by estimates of the extinction rate and the extinction fraction ( * + ), which typically have wide confidence intervals even for asymptotically identifiable models (see e.g. Table S9 in [START_REF] Stadler | Mammalian phylogeny reveals recent diversification rate shifts[END_REF]), such that accurate estimates often require sample sizes that are not achieved in practice. Speciation rates, on the other hand, can be estimated with good accuracy on phylogenies of moderate size for the constant-rate BD model [START_REF] Stadler | How Can We Improve Accuracy of Macroevolutionary Rate Estimates?[END_REF], as well as for the BDS model if the number of epochs is kept small [START_REF] Stadler | Mammalian phylogeny reveals recent diversification rate shifts[END_REF]. Similarly, in BD models with rates that are constrained to follow a specific and simple functional dependency (such as the exponential) to time [START_REF] Rabosky | Explosive Evolutionary Radiations: Decreasing Speciation or Increasing Extinction Through Time?[END_REF]15] or the environment [START_REF] Lewitus | Detecting Environment-Dependent Diversification From Phylogenies: A Simulation Study and Some Empirical Illustrations[END_REF], parameters determining the time-or environmentdependency of the extinction rate have wide confidence intervals, while those associated with the speciation rate can be estimated with good accuracy [15,[START_REF] Lewitus | Detecting Environment-Dependent Diversification From Phylogenies: A Simulation Study and Some Empirical Illustrations[END_REF]. However, by the usual arguments about degrees of freedom, the functional complexity that can be supported by a typically-sized phylogeny of a few hundred tips is probably in the order of a few parameters. Thus, practical identifiability alone dictates that we must put constraints on the flexibility of the models used to infer diversification dynamics.

Dealing with practical versus asymptotic identifiability issues

Asymptotic and practical identifiability issues are common in science, and a large set of ideas has emerged to address such problems. Practical identifiability issues are commonly understood as manifestations of the bias-variance trade-off, which states that model complexity must be adjusted to the data size to minimize the total error (bias + variance) of the inference (Box 2). This can be achieved by a variety of statistical model selection or regularization techniques (Box 2). For example, the practical identifiability of the asymptotically identifiable BDS model (without mass extinctions) can be improved by introducing temporally-autocorrelated rates drawn from a Bayesian prior, rendering parameter estimates with time divided in hundreds of epochs identifiable on relatively small phylogenies (200 tips) [START_REF] Magee | Locally adaptive Bayesian birth-death model successfully detects slow and rapid rate shifts[END_REF].

Addressing asymptotic identifiability issues, such as the non-identifiability of the BD model with unconstrained 𝜆 and 𝜇 highlighted by Louca & Pennell [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF], is a different problem, as the error of our inference does not decrease with increasing data size. Yet there are approaches for dealing with asymptotic identifiability as well, that we detail below.

Reparametrization

A solution to asymptotic identifiability issues is to reparameterize the model with identifiable quantities. For example, in the BD model with incomplete sampling and free 𝜌 (which needs to be considered when total diversity is unknown, which is the case of most microbial and insect groups), the net diversification rate 𝜆 -𝜇 and 𝜆𝜌 are identifiable. The drawback of this approach, however, is that the reparameterized quantities are often scientifically less interesting. For example, Louca & Pennell (2020) [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF] suggest estimating the pulled speciation and diversification rates 𝜆 D and 𝑟 D instead of 𝜆(𝑡) and 𝜇(𝑡) (Box 1), but these pulled rates are difficult to interpret biologically (see [30] for an attempt), which considerably limits their practical utility.

Independent data sources

Another approach to dealing with asymptotic identifiability issues is to add additional, independent data sources. Considerable progress has been made in recent years to use both phylogenetic and fossil data, which is achieved by adding fossil sampling processes to the BD process [31][START_REF] Gupta | The probability distribution of the reconstructed phylogenetic tree with occurrence data[END_REF][START_REF] Manceau | The probability distribution of the ancestral population size conditioned on the reconstructed phylogenetic tree with occurrence data[END_REF][START_REF] Gavryushkina | Bayesian Inference of Sampled Ancestor Trees for Epidemiology and Fossil Calibration[END_REF][35][START_REF] Stadler | Sampling-through-time in birth-death trees[END_REF][START_REF] Didier | The reconstructed evolutionary process with the fossil record[END_REF][38]. In the most elaborate versions of these "fossilized" Birth-Death (FBD) models, two distinct sampling processes are considered: one with rate 𝜓 for fossils with character (or molecular) data, which are included in the tree, and one for simple fossil occurences without character data. The former process is asymptotically identifiable when 𝜆, 𝜇, and 𝜓 are constant [START_REF] Gavryushkina | Bayesian Inference of Sampled Ancestor Trees for Epidemiology and Fossil Calibration[END_REF], unless samples are removed upon sampling [START_REF] Gavryushkina | Bayesian Inference of Sampled Ancestor Trees for Epidemiology and Fossil Calibration[END_REF][START_REF] Louca | Fundamental Identifiability Limits in Molecular Epidemiology[END_REF]. The latter, however, is irrelevant in the case of modeling diversification dynamics, as extinctions and fossilizations are independent processes. As long as samples are not removed upon sampling, the process remains identifiable even if the sampling probability at present 𝜌 is unknown (a case when the process is not identifiable from extant species alone), which illustrates that fossils can alleviate identifiability issues [START_REF] Gavryushkina | Bayesian Inference of Sampled Ancestor Trees for Epidemiology and Fossil Calibration[END_REF].

Despite these encouraging results, more work is needed to determine if and under which circumstances the FBD process is identifiable when 𝜆, 𝜇, and 𝜓 vary as piecewise constant or continuous functions of time, to assemble empirical datasets on which to apply FBD models for diversification analyses (the FBD has so far mainly been applied to improve divergence times rather than diversification rate estimates, but see e.g. [35,[START_REF] May | Inferring the Total-Evidence Timescale of Marattialean Fern Evolution in the Face of Model Sensitivity[END_REF]), to improve their computational efficiency (current implementations limit the applicability of the model to small datasets), as well as to assess whether the inclusion of fossils provide realistic estimates of extinction rates [START_REF] Marshall | Five palaeobiological laws needed to understand the evolution of the living biota[END_REF] (see Outstanding questions).

Constraints from a priori hypotheses

Identifiability issues are more likely to arise the more flexible our models are. Flexibility is put to the extreme by Louca & Pennell [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF], who set the task to be able to identify any possible functional forms 𝜆(𝑡), 𝜇(𝑡) from extant phylogenies. A hypothesis-driven research framework limits this complexity by comparing only a small number of alternative a priori ideas about the underlying process [42]. Such a priori hypotheses will usually constrain the functional forms of 𝜆 and 𝜇 and thus render the corresponding BD models identifiable.

The foundational study of Nee et al. [START_REF] Nee | Tempo and mode of evolution revealed from molecular phylogenies[END_REF] followed such a hypothesis-driven philosophy. After demonstrating that their bird phylogeny was incompatible with a constant-rate diversification model and grounded in Simpson's evolutionary theory of adaptive radiations [START_REF] Simpson | The Major Features of Evolution[END_REF], they hypothesized that rates of cladogenesis might be affected by niche-filling processes. Finding that a diversity-dependent model indeed fitted their data better, they concluded that diversitydependent cladogenesis was a more plausible scenario to explain the diversification of birds. This hypothesis-driven approach has inspired more than 30 years of research in phylogenetic diversification analyses [10]. Exponential time-dependencies have been used, for example, to mimic early burst patterns expected from adaptive radiation theory [START_REF] Simpson | The Major Features of Evolution[END_REF], or as an approximation to diversity-dependent cladogenesis [START_REF] Rabosky | Density-dependent diversification in North American wood warblers[END_REF] (see Box 3 for an illustration with the Madagascan vangas, Vangidae). In the context of environment-dependent models, functional hypotheses have often been derived from foundational theories of biodiversity, such as the metabolic theory of biodiversity [18] and MacArthur & Wilson's theory of island biogeography [START_REF] Valente | A simple dynamic model explains the diversity of island birds worldwide[END_REF]. Phenomenological models, such as simple linear time-or environmentaldependencies, have also been used, but typically either as null models [START_REF] Rabosky | Density-dependent diversification in North American wood warblers[END_REF] or as the simplest way to model the effect of an explanatory environmental variable on evolutionary rates [18].

The primary goal of this research has been to fit, test and compare diversification scenarios that were defined a priori to reflect different evolutionary hypotheses. Louca & Pennell's congruent models do not correspond a priori to any evolutionary hypotheses, and would never be considered in a hypothesis-driven model selection procedure in the first place [42] (Box 3).

A drawback of hypothesis-driven research is that the biological conclusions we draw are contingent on the a priori hypotheses we formulate. In particular, our hypotheses typically do not correspond completely to the process underlying the empirical data ("the truth"). Still, it is usually assumed that if a given hypothesis is statistically supported within a well-chosen set of alternatives, it is likely that this hypothesis is the closest to the truth. Whether this is the case for BD models, considering the existence of a large number of congruent models, remains an open question to be explored in more details (see The future of phylogenetic-based diversification research and Outstanding questions).

Constraints on complexity and statistical regularization techniques

Even in the absence of additional data or a priori hypotheses, there are certain philosophical, statistical or information-theoretic principles that may allow us to prefer some congruent solutions over others.

For example, a widely accepted scientific method of deciding between alternative explanations is the principle of parsimony (or Occam's razor, Box 2). If we follow this traditional thinking in science, when several explanations with different degrees of complexity are asymptotically unidentifiable, we should prefer the simplest, which is most probably true, all other things equal. A possible solution to the identifiability issue highlighted by Louca & Pennell [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF] consists then in selecting the simplest diversification scenario in a congruence class. This preference for simplicity is distinct from the problem of optimizing complexity to avoid overfitting in the case of finite data, and applies to the case of infinite data as well. Quantifying and penalizing complexity can be challenging, but it is a classical problem that can be addressed with a variety of statistical regularization techniques (Box 2).

Penalizing complexity is just one example of a more general class of regularization techniques that add additional constraints to solve an ill-posed (for example asymptotically unidentifiable) problem [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF]. Constraints can also come from prior biological knowledge, information theory or model selection principles, added in the statistical inference in the form of shrinkage estimators [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], or as priors in the case of Bayesian inference (Box 2). For example, as shown by May et al. [START_REF] May | A Bayesian approach for detecting the impact of massextinction events on molecular phylogenies when rates of lineage diversification may vary[END_REF], using Bayesian priors that represent the prior belief that on average 10% of species survive a mass-extinction event in the BDS model with mass extinction events (an asymptotically unidentifiable model) allows distinguishing rate shifts from mass extinction events. This example provides a clear counter-example to the conclusion of Louca & Pennell that regularization cannot solve asymptotic identifiability issues ([39], S2.2). Another well-known example in phylogenetics is the dating of divergence times: substitution rates and time are unidentifiable with only sequence data from extant species, but Bayesian priors on divergence times (e.g. informed by fossils) combined with relaxed clock models solve this issue (see, e.g. Fig. 1 in [START_REF] Dos Reis | Bayesian molecular clock dating of species divergences in the genomics era[END_REF]).

The future of phylogenetic-based diversification research

The asymptotic non-identifiability of the homogeneous BD process led Louca & Pennell [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF] to conclude that phylogenetic-based diversification research should switch from a focus on speciation and extinction rates to a focus on the identifiable pulled rates. Yet, scientists interested in testing specific evolutionary hypotheses would have a hard time formulating their hypotheses in terms of these quantities, which do not correspond to a particular biological mechanism. Moreover, estimating these rates from limited-size phylogenies is still a challenging task (SI Text S2 & S3).

Instead of abandoning the goal of developing models with explicit hypotheses on speciation and extinction rates, we argue to put more efforts in using all available data (including fossil data), and testing how robust the inference from these models really is in practice, when using either a hypothesis-driven research approach, or appropriate statistical regularization techniques (Fig. 1). In this area, two key questions remain: how robust are biological conclusions in practice, when we use a hypothesis-driven research framework, given the existence of congruence classes? And can parsimony considerations or other regularizing techniques successfully shrink solutions in the congruence class towards the truth? The answer to these questions depends on the nature of congruence classes, for example on whether congruence classes typically contain a wide range of disjunct models that all correspond to reasonable biological hypotheses, or that have similar parsimony/regularization properties, which remains to be explored by future research.

We can think of several ways to explore these questions, such as: i) Studying the geometric properties of congruence classes mathematically, as L&P have started to do but without definitive conclusions (their S. 1.8). This would help make the regularization choices most likely to render the models identifiable. ii) Simulating phylogenies under general ecoevolutionary models [49-51] and checking whether the application of a hypothesis-driven framework (with well-chosen a priori hypotheses) selects the hypothesis that best captures a given simulated scenario; in comparison to the simulation analyses that are already usually performed to evaluate the power and type I error rates of newly-developed methods, in which simulations correspond exactly to one of the fitted models, this requires using less idealized simulation models representing the eco-evolutionary processes that shape diversification dynamics. iii) Pursuing current efforts to develop regularized models, as detailed in the following paragraph, and use eco-evolutionary simulations (as in ii) to check whether these models provide estimates of speciation and extinction rates that approach simulated rates. Moreover, in real applications, practical identifiability is often as much a problem as asymptotic identifiability. Given that regularization can solve practical as well as asymptotic identifiability issues, developing suitable and biologically motivated regularization approaches that act directly on speciation and extinction rates seems more promising to us.

Such approaches have already started to be developed (e.g. [START_REF] May | A Bayesian approach for detecting the impact of massextinction events on molecular phylogenies when rates of lineage diversification may vary[END_REF][START_REF] Magee | Locally adaptive Bayesian birth-death model successfully detects slow and rapid rate shifts[END_REF]), and including further general ideas from statistics and machine learning, for example the fused lasso [START_REF] Tibshirani | Regression Shrinkage and Selection Via the Lasso[END_REF] or generalized additive models (GAMs, [START_REF] Hastie | Generalized additive models[END_REF][START_REF] Wood | Generalized Additive Models: An Introduction with R, Second Edition[END_REF]) could lead to further advances (Box 2).

The problems as well as their solutions discussed here are likely not limited to homogeneous BD models. In recent years, models with diversification rates that vary across lineages have been developed to understand why some groups of organisms are richer than others and to avoid biased inferences linked to model misspecification [15,[START_REF] Alfaro | Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates[END_REF][START_REF] Rabosky | Automatic Detection of Key Innovations, Rate Shifts, and Diversity-Dependence on Phylogenetic Trees[END_REF][57][START_REF] Ronquist | Universal probabilistic programming offers a powerful approach to statistical phylogenetics[END_REF][59]. Unlike for the homogeneous BD model, for which all topologies are equally likely and therefore only branching times are informative, both branching times and topology are informative in the case of heterogeneous BD models. Despite this additional source of information, it is very likely that models with heterogeneous rates are asymptotically unidentifiable in the absence of any constraint. Working with biologically interpretable speciation and extinction rates has helped regularizing this problem, for example by favoring rare rate shifts with large effects corresponding to the invasion of new ecological space [START_REF] Alfaro | Nine exceptional radiations plus high turnover explain species diversity in jawed vertebrates[END_REF][START_REF] Rabosky | Automatic Detection of Key Innovations, Rate Shifts, and Diversity-Dependence on Phylogenetic Trees[END_REF][57] or by favoring frequent shifts with small effects corresponding to heritable rates, formalized by regularization in the form of autocorrelated Bayesian priors [59,[START_REF] Maliet | Fast and Accurate Estimation of Species-Specific Diversification Rates Using Data Augmentation[END_REF].

Concluding Remarks

Identifiability issues naturally arise in approaches that try to infer the potentially unlimited complexity of historical processes from limited contemporary data, and inference of past diversification history from phylogenies of present-day species is no exception. These identifiability issues are one of the reasons why scientists adhere to hypothesis-driven research, use parsimony or regularization principles, or integrate multiple data types.

Phylogenetic-based diversification analyses have already adopted these methods in the past, and need to pursue this effort to provide increasingly robust tools for understanding past diversification histories from the data that is available today (see Outstanding Questions).

Box 1: Model congruency and pulled diversification rates

Louca & Pennell [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF] consider the homogeneous (i.e. lineage-independent) stochastic birthdeath process of cladogenesis with rates of speciation (birth, 𝜆) and extinction (death, 𝜇) that can change arbitrarily over time 𝑡. They show that for any given derivable (and therefore continuous) time-dependent speciation function 𝜆 > 0 and extinction function 𝜇 ≥ 0, there exists an infinite set of alternative functions 𝜆 * > 0 and 𝜇 * ≥ 0 such that the probability distribution of extant trees under the corresponding birth-death processes M and M * is identical. Consequently, M or M * yield identical likelihood values for any given empirical tree, which implies that 𝜆(𝑡) and 𝜇(𝑡) are not uniquely identifiable unless further constraints are imposed on their functional form.

Louca & Pennell then re-parameterize the problem to have only identifiable quantities, which they call the pulled rates. The pulled speciation rate is given by:

𝜆 D = 𝜆(1 -𝜙)
where 𝜙 is a function of time that denotes the probability that a lineage alive at time 𝑡 has no descendant in the phylogeny, and which analytical expression is given, for example, by Eq.2 in [15]. The pulled diversification rate is given by:

𝑟 D = 𝜆 -𝜇 + Q + R+ RS
Congruent models are found by solving Eq. 2 in [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF]: where the penalty term controls the "strength" of the preference for simplicity.

𝑑𝜆 * 𝑑𝑡 =
In statistics, the traditional motivation to favor simplicity is based on the bias-variance trade-off, which posits that increasing model complexity reduces the systematic misfit (bias), but at the cost of increasing variance (uncertainty) of the parameter estimates. One can prove that, with limited data, inducing a bias towards simpler models decreases total estimation error (bias + variance), even if the true underlying model is more complex. The complexity penalty is selected to optimize the total error. This logic underlies most frequentist regularization and model selection approaches.

There is a second argument for constraining model complexity, which is independent of the data size and the bias-variance trade-off. This argument, known as the law of parsimony or Occam's razor, relies on an a priori assumption that natural processes tend to be simple and smooth. The principle of parsimony is not a mathematically provable law, but it underlies centuries of thinking and experience from physics to machine learning, and from philosophy as well (see [START_REF] Coelho | A parsimonious view of the parsimony principle in ecology and evolution[END_REF] for a discussion).

When implementing preferences for simplicity, it typically makes no difference if they originate from bias-variance or parsimony principles. The main difference is that in the former the penalty is chosen from the data, such that more complex models are preferred as the data size increases, whereas in the later the penalty is chosen independent of the data, based on prior beliefs. How to best define complexity is a question of constant debate and development in statistics: we may for example decide that a model is simple if it is interpretable, if it involves less parameters, if it prevents fast variations, or yet other criteria.

Various statistical regularization techniques implementing these criteria exist. For example, information-theoretical measures (e.g. the AIC or BIC, [42,[START_REF] Aho | Model selection for ecologists: the worldviews of AIC and BIC[END_REF]) add a direct penalty for the number of parameters, shrinkage estimators such as lasso or ridge or their corresponding Bayesian priors add a penalty on the deviation of model parameters from zero [START_REF] Tibshirani | Regression Shrinkage and Selection Via the Lasso[END_REF] and statistical smoothers [START_REF] Wood | Smoothing Parameter and Model Selection for General Smooth Models[END_REF] penalize the roughness of the fitted model (as in generalized additive models GAMs, see [START_REF] Hastie | Generalized additive models[END_REF][START_REF] Wood | Generalized Additive Models: An Introduction with R, Second Edition[END_REF]).

Box 3: Diversification of the Madagascan vangas

We illustrate hypothesis-driven research by performing an analysis of the diversification of the Madagascan vangas (Vangidae) using the logic that would be applied in the field [START_REF] Jønsson | Ecological and evolutionary determinants for the adaptive radiation of the Madagascan vangas[END_REF], but simplified for illustrative purposes. We hypothesize that diversification followed an 'Early Burst' pattern [START_REF] Harmon | Early Bursts of Body Size and Shape Evolution Are Rare in Comparative Data[END_REF], with fast speciation at the origin of the group and subsequent slowdown, rather than constant-rate diversification. The Early Burst pattern, related to the idea of adaptive radiations [START_REF] Simpson | The Major Features of Evolution[END_REF], is modeled by an exponential decay of the speciation rates through time, used as an approximation of diversity-dependence. We also consider the hypothesis that a substantial number of extinction events occurred during the diversification of this group.

Among the four corresponding models, the model with an exponentially declining speciation rate 𝜆(𝑡) = 𝜆 X 𝑒 aS (time 𝑡 is measured from the present to the past), with speciation rate at present 𝜆 X = 0.018, rate of decline 𝛼 = 0.1 and no extinction 𝜇(𝑡) = 0, noted M, is best supported by the data (see SI Table S1). We conclude that the hypothesis of Early Burst diversification with negligible extinctions is the most likely of the four hypotheses we considered.

In order to better grasp the nature of congruent models, we explore models congruent to our best model M (see SI Text S1). First, we choose the extinction function to be a constant 𝜇 Q * (𝑡) = 𝜇 X and compute 𝜆 Q * (𝑡). Second, we choose the speciation function to be a constant The black curves represent the dynamics of speciation (solid line) and extinction (dashed line) corresponding to our best-fit model M (exponential decline in speciation rate, non-significant extinctions). The colored curves illustrate the rate dynamics of congruent models that were obtained by (A) fixing increasing values of a constant extinction rate (M1 * ) and (B) fixing the speciation rate to 𝜆 X (M2 * ). In the case of constant extinction (A), we can choose any value for 𝜇 X and find 𝜆 Q * (𝑡) (so there is an infinity of congruent models), while in the case of constant speciation (B), 𝜆 X needs to be taken equal to the 𝜆 X of model M , as two congruent models necessarily have the same speciation rate at present if 𝜌 is fixed [START_REF] Louca | Extant timetrees are consistent with a myriad of diversification histories[END_REF] (so there is only one congruent model). Note that M1 * infers a speciation rate decline regardless of the assumed extinction rate and that M2 * infers biologically implausible negative extinction rates. See also Online Supplemental Information Figure S1 & S2. small number of biologically motivated hypotheses (H1 to H4), the models will usually be identifiable, meaning that the optimum solution under a given hypothesis is unique (one black circle per hypothesis), and we will select the hypothesis that comes closest to the congruence class (here, H1, dashed lines convey the distance to the congruence class). This hypothesis, which is the one with highest likelihood, is traditionally assumed to be the closest to the true process. B: Parsimony and regularization assumptions constrain the congruence class (grey circle). From the experience in other fields, we would expect the congruence class to be constrained towards the true process. These two expectations are likely to be met if biologically and statistically (i.e. with respect to parsimony and regularity properties) reasonable models within the congruence class cluster around the true process. Whether this assumption holds in reality is a question for future research.

observations

  Bias-variance trade-off: trade-off between systematic model error due to limited flexibility (bias) and uncertainty of the parameter estimates (variance) Extinction: disappearance of a species, corresponding to the death of its last individual Homogeneous birth-death (BD) process: birth-death process where speciation and extinction rates are identical across lineages at any time. Rates may vary in time, but not across lineages Identifiability: when fitting statistical models, identifiability means that any two combinations of parameter values can be distinguished Likelihood: function of the parameters of a given model that measures the probability of the observations given the model and its parameter values Model misspecification: situation when the distribution of data implied by the model (under best possible parameterization) differs from the distribution of data under the true generating process Net diversification rate: speciation rate minus extinction rate Practical unidentifiability: situation when there are distinct combinations of the model parameters that cannot be told apart from the limited number of observations available in practice Reconstructed phylogeny: estimated phylogenetic tree for present-day species (missing lineages that have gone extinct and are thus unsampled) Regularization: set of statistical techniques that consist in adding a regularization term (or penalty) to the optimization function (typically the likelihood) to solve an ill-posed problem or avoid overfitting Speciation: process by which two populations of the same ancestral species give rise to two distinct descendant species Extinction fraction: extinction rate divided by speciation rate

Box 2 -

 2 𝜆 * U𝑟 D -𝜆 * + 𝜇 * V Given any 𝜇 * , we can compute 𝜆 * using the solution to this equation, provided in Louca & Pennell [23]'s SI (Eq. 39 & 40, 𝜂 X = 𝜌𝜆(0), 𝜇 X = 𝜇(0)): any 𝜆 * , we can compute 𝜇 * as: Reasons and approaches to select simple models Deciding between alternative hypotheses through a preference for simplicity is ubiquitous in statistics and the sciences. Mathematically, this is expressed by viewing the evidence in favor of a respective hypothesis (or model, denoted by M) as a combination of: Evidence = Likelihood(M) -Penalty*Complexity(M)

  ) = 𝜆 X and compute 𝜇 d * (𝑡). We find (SI Text S1; Fig.I; here we take 𝜌 = 1 as the tree of the Madagascan vangas is complete[START_REF] Jønsson | Ecological and evolutionary determinants for the adaptive radiation of the Madagascan vangas[END_REF]): ) = 𝜆 X -𝛼 -𝜆 X 𝑒 aSThe biological interpretation of these models and of their parameters is not obvious. The equation for 𝜇 d * looks more interpretable at first, but it expresses the temporal change and the extinction rate at present through the same parameter 𝛼, which means that a positive extinction rate at present (𝛼 < 0) will force extinction rates to decline over time. Here M2 * infers negative extinction rates, and is therefore not plausible (Fig.I). M1 * infers a decline in speciation rate from the origin of the group to the present for extinction rates 𝜇 X ranging from at least 0.05 to 0.3, consistent with our previous results (Fig.I). While rate estimates do vary substantially, the general temporal trend is preserved.

Fig I .

 . Fig I. Diversification of the Madagascan vangas as inferred from congruent models.

Fig. 1 .

 1 Fig. 1. Conceptual figure illustrating how constraints imposed by prior hypotheses and regularization may help to approach the true process. Following Fig. 3 in [23], the pink area represents the congruence class of the true process (red circle). A: When considering a

Acknowledgments

We thank Amaury Lambert and the Morlon (current and past) research team, in particular Jeremy Andréoletti, Joelle Barido-Sottani, Julien Clavel, Jonathan P. Drury, Isaac Overcast, Odile Maliet, and Ignacio Quintero for discussions. We also thank Knud Jønsson for providing us with the Magagascan vanga tree. H. Morlon acknowledges funding from the ERC (grant ERC-CoG PANDA).