Yalchin R Efendiev
email: efendiev@math.tamu.edu

Lev Truskinovsky

a driven bi-stable FPU chain

Keywords: Fermi-Pasta-Ulam system, Bi-stable elements, Thermolization, Ergodicity, Statistical mechanics, Latent heat

come

Introduction

About 60 years ago Fermi et al. [START_REF] Fermi | Studies of Nonlinear Problems[END_REF] asked whether a nonlinear chain with fixed ends and nonzero initial velocities can reach a state of thermalization and over the years this problem attracted a lot of attention (see reviews [START_REF] Berdichevsky | Thermodynamics of Chaos and Order[END_REF][START_REF] Campbell | The Fermi-Pasta-Ulam problem-the first fifty years[END_REF][START_REF] Ford | The Fermi-Pasta-Ulam problem: paradox turns discovery[END_REF][START_REF] Gallavotti | The Fermi-Pasta-Ulam Problem[END_REF][START_REF] Jackson | Perspectives of Nonlinear Dynamics[END_REF][START_REF] Weissert | The Genesis of Simulations in Dynamics: Pursuing the Fermi-Pasta-Ulam Problem[END_REF]). It was soon realized that for springs with convex energies the thermalization can be reached only if the initial kinetic energy is above a certain threshold. Below this threshold the nonlinearity is weak and the behavior of the system remains close to the quasi-periodic behavior of the integrable linear system.

The important aspect of the classical FPU nonlinearity was its linearizability. In this paper we consider non-convex springs which implies that there exists at least one state (spinodal) where the linearized model degenerates. In this state one can expect non-integrability even at infinitesimal energy of initial excitations.

To support this intuition we provide numerical and analytical evidence that for the springs with double-well energy the state of thermalization or, at least of equipartition, can be reached with zero kinetic energy input. To reach the spinodal state we subjected the chain to infinitely slow driving in a hard loading device. In this case the 'cold' system will first accumulate a certain amount of potential energy while remaining in a metastable state, and then release it through a dynamic spinodal instability. Reaching thermalization in this way would then mean that the system was brought exactly to the boundary of the KAM island of regular behavior.

To be more specific, consider a mass-spring chain governed by a dimensionless system of Newton equations. Denote by x i (t), i = 0, . . . , N the coordinates of the particles and by e( ) the elastic energy of the springs. Then

δ d 2 x i dt 2 = σ x i+1 -x i δ -σ x i -x i-1 δ (1)
where σ ( ) = e ( ). We assume zero initial velocities and apply the following time-dependent boundary conditions:

x 0 = 0, x N (t) = x 0 N + vt. ( 2
)
The system (1) depends on the two dimensionless parameters. The first parameter, measuring the degree of discreteness of the system, is δ = h/L = N -1 , where h is the interatomic distance and L is the macroscopic length of the system. The second parameter, measuring the slowness of the driving, is v = V /( √ E/ρ), where ρ is mass density, E is the characteristic elastic modulus and V is the velocity of driving. We are interested in the behavior of the system (1) at δ → 0 (thermodynamic limit) and v → 0 (quasi-static driving).

Non-convex energies at a micro-scale appear often in the modeling of materials and distributed structures with internal instabilities. A typical example is a shape memory alloy where the presence of a phase transition evokes description in terms multi-well snap-springs [START_REF] Muller | A model for an elastoplastic body[END_REF]. The energy landscape for such chains was found to be extremely rugged (e.g. [START_REF] Puglisi | Mechanics of a discrete chain with bi-stable elements[END_REF]). First numerical experiments with a fully inertial bi-stable chain [START_REF] Balk | Dynamics of chains with non-monotone stress-strain relations. I. Model and numerical experiments[END_REF] showed that the system becomes highly oscillatory after reaching the spinodal threshold. The use of a micro-canonical ensemble for a thermodynamical description of a minimal chain with two snap-springs was rationalized in [START_REF] Cardin | On the Helmholtz-Boltzmann thermodynamics of mechanical systems[END_REF]. A symmetric bi-stable chain in a soft device was studied numerically in [START_REF] Yoshikawa | Transition from mechanical equilibrium to thermal equilibrium in a chain of unstable springs[END_REF] and a canonical distribution was found to be adequate for the prediction of the macroscopic behavior. Different numerical methods of finding equilibrium properties for generic 1D lattices were recently reviewed in [START_REF] Blanc | Finite-temperature coarse-graining of one-dimensional models: mathematical analysis and computational approaches[END_REF]. Other contributions to the thermalization problem in discrete systems with non-convexity can be found in [START_REF] Dreyer | Cold, thermal and oscillator closure of the atomic chain[END_REF][START_REF] Basile | Energy transport in stochastically perturbed lattice dynamics[END_REF][START_REF] Kastner | Molecular-dynamics of a 2D model of the shape memory effect, part II: thermodynamics of a small system[END_REF][START_REF] Müller | Thermodynamic aspects of shape memory alloys[END_REF][START_REF] Weinan | Cauchy-born rule and the stability of crystalline solids: dynamic problems[END_REF].

In the present work, most of which was completed in 2001 (see [START_REF] Yoshikawa | Transition from mechanical equilibrium to thermal equilibrium in a chain of unstable springs[END_REF]), we study the emergence of ergodicity in an asymmetric bi-stable FPU chain driven quasi-statically in a hard device. In our numerical experiment we first bring the 'cold' system to the point of instability and study the equipartition of the energy. We then perform a cyclic quasi-static loading of the 'self-heated' chain through the phase transition. This leads to reversible and well reproducible temperature oscillations which we associate with the release (or absorbtion) of the latent heat. To justify the ergodicity hypothesis we compare the time-averaged behavior of the system with the statistical computations based on canonical distribution. For the case of a bi-quadratic springs we obtain explicit expression for the free energy and show that the time-averaged behavior of the chain agrees almost perfectly with the computed adiabats. We then study the emergence of hysteretic behavior in the damped chain; the overdamped limit was analyzed in [START_REF] Puglisi | Thermodynamics of rate independent plasticity[END_REF].

The paper is organized as follows: In Sect. 2 we report the results of our numerical experiments, demonstrate the equipartition of the energy, and clarify distinction between reversible and irreversible heating. In Sect. 3 we use the classical methods of statistical mechanics to derive the thermodynamic response of the system. We show that in case of bi-quadratic energy the statistical problem can be reduced to a 1D Ising problem for a special diamagnetic lattice. In Sect. 4 we compare the results of the thermodynamic analysis with direct numerical simulations. The final Sect. 5 contains a brief discussion of the under-damped dynamics.

Numerical studies

Assume that δ = 1 and denote by j = x j+1x j the strain in the jth spring, j = 0, . . . , N -1. To make analytical computations possible, we assume that the function e( ) is bi-quadratic

e( ) = μ 1 2 ( -a 1 ) 2 if ≤ a 0 μ 2 2 ( -a 2 ) 2 if > a 0 , (3)
where the spinodal strain is a 0 = (1 -κ)/(1 + κ), κ = √ μ 2 /μ 1 (see Fig. 1). While one of the parameters μ 1 or μ 2 can be chosen to be equal to 1, we preserve both of them for the symmetry of the analytical expressions. We say that the interval < a 0 belongs to Phase I and the interval > a 0 -to Phase II. The spinodal region is represented in this model by a single-point = a 0 .

We first prepare a "cold" chain with all springs initially at rest (in phase I) and then drive it quasi-statically with a constant strain rate v. In the actual numerical experiments we diminished the value of v until the averaged behavior stabilized. For numerical integration of the system (1) we used a symplectic numerical method widely known as Stormer-Verlet algorithm (see [START_REF] Stuart | Dynamical Systems and Numerical Analysis[END_REF] for more details). To quantify the observed phenomena, we need to introduce a measure of intensity of the micro-oscillations. Consider, for instance, the stage of the stretching process where all springs oscillate in a single quadratic energy well. Due to the linearity of the system we can decompose the motion into the average and the fluctuations. The average deformation is both homogeneous and quasi-static

x i (t) = (i/N ) x 0 N + vt .
The general motion can then be represented in the form

x i (t) = x i (t) + ω i (t),
where ω i (t) is the oscillatory contribution with zero average over a time interval which is much larger than the period of the microscopic oscillations ∼1 and much smaller than the characteristic time of the loading ∼ v -1 . Next, we can decompose the oscillatory motion into normal modes in the quadratic regions. Recall that ω k (t) is the solution of the linear system

d 2 dt 2 ω(t) = Aω,
where A is a tri-diagonal matrix. To define the modes we need to diagonalize the matrix A, which can be done explicitly (e.g. [START_REF] Berdichevsky | Thermodynamics of Chaos and Order[END_REF]). If we assume μ = 1 (otherwise we can re-scale the time) we obtain that the eigenvalues of A are λ l = -2 + 2 cos lπ N +1 and the corresponding eigenvectors are

sin lπ N +1 , sin 2lπ N +1 , . . . , sin l Nπ N +1 t .
Once the eigenvectors and eigenvalues are known, the matrix A can be represented as

D = U T AU, ( 4
)
where D is a diagonal matrix consisting of eigenvalues and the columns of U are corresponding eigenvectors.

The the system can be represented as a set of noninteracting harmonic oscillators d dt 2 φ = Dφ, where φ = U ω. The energy of each mode (each oscillator) is then

E l = 1 2 d dt φ l 2 + λ 2 l φ 2 l .
Consider now the symmetric bi-quadratic energy with μ 1 = μ 2 = μ. In this case the definition of the modes can be formally extended into the regime where the springs change phase and where oscillators are 2 Total energy of the modes for the system of five springs during quasistatic stretching of an initially 'cold' chain no longer independent. In Fig. 2 we show the evolution of the total energy of the modes for the system of five springs which is plotted as a function of the total elongation. As we have already mentioned, the chain is initially "cold" and the energy of the modes is equal to zero. Before the system reaches the spinodal point the springs remain in Phase I and the energy of oscillations remains equal to zero which is expected in view of the quasi-static nature of the driving. At large elongations all springs are again oscillating in a single well (Phase II) and the energy of oscillations again remains constant because the average motion decouples. At this stage the system regains its linearity and can be represented by a set of noninteracting harmonic oscillators with a randomness in initial conditions inherited from the mixed-phase dynamics. As the chain enters the nonlinear, interphase regime the springs start to oscillate between phases I and II. In this regime the micro and macro motion become coupled and, in particular, during the spinodal instability a finite amount of energy is converted from the average motion into the oscillatory motion. This energy is marked as 'available energy' in Fig. 1. One half of the available energy goes into kinetic energy of oscillations and another half goes into potential energy.

Observe that due to simplicity of the model where only the nearest neighbors interact and the elastic energy is permutationally invariant we do not see phase boundaries propagating in an organized fashion. Adding linear interaction of next to nearest neighbors of the ferromagnetic type (see [START_REF] Truskinovsky | About the origin of the nucleation peak in transformational plasticity[END_REF]) fixes this problem and localizes transformation in several phase boundaries, whose dynamics is expected to remain complex.

The next experiment is intended to show the details of the self-heating process. We prepare a 'cold' chain with i = a 0 (for all i) and then we slightly perturb one element. Note that we do not load the chain in this case. In Fig. 3 we show the ensuing motion. One can see that after the regular nucleation stage described in [START_REF] Balk | Dynamics of chains with non-monotone stress-strain relations. I. Model and numerical experiments[END_REF][START_REF] Ngan | Thermoelastic aspects of nucleation in solids[END_REF][START_REF] Truskinovsky | Nucleation and growth in elastodynamics[END_REF] the motion becomes increasingly complex due to multiple interaction of the waves reflected from the boundaries. We demonstrate the equipartition by plotting in Fig. 4 the time dependence of the averaged kinetic energies of the particles

1 2 kT i = 1 τ t 0 +τ t 0 dx i (t) dt 2 dt,
where k is a dimensional (Boltzmann) constant and τ is the averaging window. One can see that after a finite time the system achieves the state of equipartition with

T 1 = • • • = T N -1 = T.
The same conclusion can be reached when one considers the energy of the modes. Therefore, our numerical experiment indicates an important sign of ergodicity of the infinitesimally perturbed system with regular behavior. We can now study the time evolution of the temperature during several cycles of loading and unloading through the phase transformation (several subsequent elongations and contractions). For the case of equal moduli of the phases the total energy of the linear modes is shown in Fig. 5. Note that in both phases, the coordinate transformation U (see ( 4)) is the same and this transformation is used in the mixed region to compute the total energy of the linear modes. Notice again that after the initial irreversible "self heating" the system does not heat or cool any more, and the total temperature of the modes remains practically constant except in the mixed-phase regions. This is consistent with an assumption that the quasi-static driving in the hard device can be represented at the macroscale as the adiabatic evolution of the quasi-static thermodynamical system. The continuity of the temperature in the pure phases can then be viewed as the consequence of the fact that the entropies of the phases are equal. In mechanical terms one can argue that in the "hot" chain instead of massive spinodal instability we see fluctuations-induced tunneling through the barrier. This process does not require overstepping and is not accompanied by the irreversible macroscopic energy conversion into microscopic vibrations.

To confirm our thermodynamical interpretations, we studied the case when the elastic moduli in the two energy wells are different. Then the phase with more shallow energy minimum will have a higher entropy than the phase with the more narrow minimum, and transition from former to latter in adiabatic conditions would require a decrease of temperature. Indeed, the experiments for the chains with μ 1 = μ 2 show that the value of the temperature (represented by the total energy of the modes defined in each of the phase) oscillates in the predicted way as one switch between the phases (see Fig. 6). More precisely, our numerical experiments indicate that the ratio

E μ 1 /E μ 2 ≈ √ μ 1 /μ 2 ,
where E μ i is the total energy of the modes in the single phase stage of deformation with stiffness μ i , (i = 1, 2). In terms of the temperature of the particles this relation translates into

T μ 1 /T μ 2 ≈ μ 1 /μ 2 .
In the next section we show that this numerical result is fully supported by the computations of the time-averaged behavior based on ergodicity hypothesis.

To check the ergodicity of the system one would like to compare statistical mechanical computations with the macroscopic stress-strain and energy-strain relations obtained from the direct numerical simulations for the self-heated chain driven quasi-statically through the phase transition. To perform time averaging and filter oscillations, we used the window size equal to about 30 oscillations. The results for the system with symmetric energy wells are presented in Fig. 7, where we show the macroscopic stress and potential energy for the system of 15 springs. According to our assumption those curves describe the adiabats. For comparison, we plot them against the corresponding curves obtained by the global minimization of the energy (e.g. [START_REF] Fedelich | Hysteresis in discrete systems of possibly interacting elements with a double-well energy[END_REF]). Notice that in the two-phase region, where the purely mechanical description predicts zero elastic modulus, our numerical data suggest nonzero "entropic" elasticity. The numerically computed macroscopic potential energy is close to being convex (as can be expected in a 1D system) and is always above the energy of the mechanical ground state due to the potential energy of oscillations (half of the 'available energy'). The goal of the statistical theory presented in the next section is to qualitatively reproduce these numerical results.

Statistical modeling

Under the assumption of ergodicity at the microscale one can use the techniques of statistical mechanics to compute the macroscopic constitutive response of the chain. For our finite 'self-heated' system the exact nature of the invariant measure remains to be established; however, the presence of equipartition suggests the applicability of the canonical distribution at least when the number of elements is sufficiently large (e.g [START_REF] Berdichevsky | Thermodynamics of Chaos and Order[END_REF]). To define the ensemble we fix the total strain = N i=1 i and denote the mean strain by ¯ = /N . We introduce a canonical measure by assuming that a microstate with the energy E = N i=1

1 2 p 2 i + e( i ) has a probability f ( p i , i ) = 1 Z exp (-β E)
, where β = (kT ) -1 is the new macroscopic parameter (k is the Boltzmann constant). The statistical sum Z is defined as a normalization of the probability distribution

Z (¯ , β) = exp(-β E)d p 1 , . . . , d p N -1 d 1 , . . . , d N -1
and the integration is carried over all micro-configurations characterized by N -1 strains i and N -1 momenta p i (we assume that p 0 = 0 and p N = 0).

The partition function can be immediately factorized into kinetic and configurational parts Z = Z k Z p . Due to the independence of momenta, the computation of the kinetic energy contribution is straightforward

Z k (β) = (2π/β) N -1/2 . The configurational contribution Z p is given by Z p (¯ , β) = 1 ,..., N , i i =N ¯ exp - 1 2 β E p d 1 , . . . , N . ( 5
)
where E p = N i=1 e( i ) is the total potential energy of the system. After the statistical sum is known, the free energy can be computed as follows:

F(¯ , β) = - 1 β log(Z ) = F k + F p , (6)
where

F k = - 1 β (N -1) 2 log (2π) β , F p = - 1 β log(Z p ).
The expression for the entropy takes the form

S(¯ , β) = k N -1 2 1 + log 2π β + k log(Z p ) -β ∂ log(Z p (β, N )) ∂β . ( 7
)

Single phase system

To illustrate the above formulas we first consider quadratic potential with e( ) = 1 2 μ 2 . In this case the problem reduces to the computation of a Gaussian integral

Z p = 1 ,..., N , i i =N ¯ exp - 1 2 βμ i 2 i d 1 , . . . , N . ( 8
)
As is common in such calculations we introduce the new variables r 0 = 0, r k = k i=1 i , and then use the known formulas

r i exp -β μ 2 N i=1 (r i -r i-1 ) 2 dr 1 , . . . , dr N = N 1/2 2π βμ N -1 2 r i exp -β μ 2N (r N -r 0 ) 2 dr 1 . . . , dr N . ( 9
)
to obtain

Z p = N 1/2 2π βμ N -1 2 exp - 1 2 βμN 2 .
The average potential energy takes the form E p = (N -1) 2β + 1 2 N μ¯ 2 . At large N the energy density E p /N splits into the macroscopic elastic energy 1 2 μ 2 and the microscopic energy 1/(2β) = kT /2. The microscopic energy is due to fluctuations which behave as an ideal gas. The specific heat at constant strain is equal to c = (∂ E/∂ T ) = k N = R, with half of it due to kinetic energy of oscillations and half due to their potential energy. For the free energy function we obtain

F p (¯ , β) = 1 (2β) log(N ) - N -1 (2β) log 2π (βμ) + 1 2 N μ¯ 2 .
One can see that the stress-strain relation at constant temperature

σ = ∂(F p /N )/∂ ¯ = μ
remains purely mechanical. The expression for the finite size entropy takes the form

S = k (N -1) 2 log (4π 2 ) (β 2 μ) - 1 2 log(N ) + N -1 . ( 10
)
We can use the last formula to explain the equilibrium temperature jump as the "hot" asymmetric chain is driven quasi-statically through the phase transition. Indeed, due to adiabatic nature of the loading the entropies of the two limiting linear regimes must be identical. If the corresponding elastic moduli are different the equilibrium temperatures will also be different due to the release or absorption of the latent heat. From (10) we obtain

β 2 1 μ 1 = β 2 2 μ 2 .
which agrees with the relation between the limiting temperatures found in our numerical experiments.

Two-phase system

For the two-phase system with a bi-quadratic energy we can write

E p = N i=1 μ i 2 ( i + α i ) 2 .
Here the discrete parameter α which locates the bottoms of the energy wells takes only two values -a 1 and -a 2 . Similarly, μ i takes values μ 1 or μ 2 depending on the value of i . Note that both μ i and α i depend on i .

For the sake of analytical transparency we can start with a simplified model where we expand the quadratic energy wells beyond the spinodal point into the whole real line and neglect the nonlinear switching condition at i = 0 . Such approximation implying the multivalued energy can be justified in the limit of small temperatures where the contribution from the non single-valued response is negligible. Our numerical results confirm that the simplified model is accurate for the values of temperature generated during the 'self-heating' of the 'cold' chain loaded quasi-statically through the phase transition.

More specifically, the simplified model adopts the following approximate expression for the configurational partition function:

Z p (¯ , β) = α i 1 ,..., N , i i =N ¯ exp - 1 2 β E p d 1 , . . . , N . (11)
To compute the integrals, we can again introduce the variables r k

r 0 = 0, r k = k i=1 i + k i=1 α i . ( 12
)
In terms of these new variables, the function Z p can be written in the form

Z p = α i r 1 . . . r N -1 dr 1 , . . . , dr N -1 exp - β 2 N i=1 μ i (r i -r i-1 ) 2 . ( 13
)
To proceed further, we denote

P(μ) = ⎡ ⎢ ⎢ ⎢ ⎣ μ 1 + μ 2 μ 2 0 • • • 0 μ 2 μ 2 + μ 3 μ 3 • • • 0 0 μ 3 μ 3 + μ 4 • • • 0 • • • • • • • • • • • • • • • 0 0 0 μ N -1 μ N -1 + μ N ⎤ ⎥ ⎥ ⎥ ⎦ , φ(μ) = (0, . . . , 0, μ N r N ). ( 14
)
Then ( 13) can be rewritten as

Z p = α i r 1 . . . r N -1 exp - β 2 ((P(μ)r, r )) + β(φ(μ), r ) dr 1 , . . . , dr N -1 exp - β 2 μ N r 2 N , (15)
where (•, •) denotes scalar product and r = (r 1 , . . . , r N -1 ). Next, we use again the well known property of Gaussian integrals exp

⎛ ⎝ - 1 2 i, j P i, j x i x j + k φ k x k ⎞ ⎠ dx 1 , . . . , dx N = (2π) N /2 √ det(P) exp ⎛ ⎝ 1 2 i, j P -1 i, j φ i φ j ⎞ ⎠ , ( 16
)
to obtain

Z p = α i (2π) N -1 2 β N -1 2 √ det(P) exp 1 2 β(μ N P -1 N -1,N -1 -1)μ N r 2 N , ( 17
)
where P -1 N -1,N -1 is the last entry (last column and last row) of the inverse of P(μ). Now, noting that r N = + i α i , the expression for Z p can be rewritten as

Z p = α i (2π) N -1 2 β N -1 2 √ det(P) exp 1 2 β(μ N P -1 N -1,N -1 -1)μ N ( + i α i ) 2 . ( 18
)
Every term that involves P depends on μ i and consequently on α i . Thus det(P) can be evaluated explicitly det(P) = N i=1 D/μ i , where D = μ 1 . . .μ N . The term P -1 N -1,N -1 can be written in the closed form (see [START_REF] Da Fonseca | Explicit inverse of a tridiagonal k-toeplitz matrix[END_REF])

P -1 N -1,N -1 = μ N -1 θ N -1 /θ N ,
where

θ i = (μ i + μ i+1 )θ i-1 -μ 2 i-1 θ i-2 , θ 0 = 1, θ -1 = 0. Denoting η i = θ i /θ i-1
, we obtain the following recursive relation

η i = 1 μ i + μ i+1 -μ 2 i-1 η i-1 , η 1 = 1 μ 1 + μ 2
We then note that

P N -1,N -1 = μ N -1 η N ,
and by induction argument show that

P N -1,N -1 = N -1 i=1 μ i N -1 i=1 1 μ i N i=1 μ i N i=1 1 μ i .
Using these expressions we finally obtain

Z p = α i (2π) N -1 2 β N -1 2 N i=1 μ i N i=1 1 μ i 1/2 exp ⎛ ⎝ - 1 2 β N i=1 1 μ i -1 ( + i α i ) 2 ⎞ ⎠ . ( 19
)
Notice that we were able to eliminate all linear variables (strains) and reduce the problem to solving a particular (diamagnetic) one-dimensional Ising model. The mean field nature of this model is due to absence of interactions beyond nearest neighbors, when the discrete statistical sum ( 19) is invariant with respect to the permutation of μ i and can be rewritten in terms of the number of springs in the first phase p. First we observe that

N i=1 μ i = μ p 1 μ N -p 2 , N i=1 1 μ i = p μ 1 + N -p μ 2 , α i = pa 1 + (N -p)a 2 .
Then, we can write the following explicit expression capturing the size effect:

Z p ( , β) = N -1/2 (2π) N -1 2 β N -1 2 N p=0 N ! p!(N -p)! (μ p/N 1 μ 1-p/N 2 ) -N /2 p/N μ 1 + 1 -p/N μ 2 -1/2 × exp - β 2 N p/N μ 1 + 1 -p/N μ 2 -1 + p N a 1 + 1 - p N a 2 2 . ( 20
)

Mechanical limit

To make connection with the purely mechanical theory based on the global energy minimization of the elastic energy of the chain (ground state response) we need to consider the limit β → ∞ when the simplified model becomes asymptotically exact. In this limit the sum (20) can be computed analytically using Laplace method.

More precisely, we can use the fact that in the limit of large β one can approximate the partition function Z p by exp(-

β P(x))dx = exp(-β P(x 0 )) (2π) N det(A) , ( 21
)
where the matrix A is defined by

A i j = -β ∂ 2 ∂ x i ∂ x j P(x)| x=x 0 .
and x 0 is the point of the minimum of P(x). Now, suppose for simplicity that we are dealing with equal moduli case and assume that E( 1 , . . . , N ) has minimum E min . In the neighborhood of the minimum E has quadratic behavior with

∂ 2 ∂ i ∂ j E = μδ i j ,
therefore, by using the Laplace formula ( 21) we obtain

Z p = 2π βμ N 2 α i exp (-β E min ). ( 22
)
By using the explicit expression for the ground state configurations [START_REF] Fedelich | Hysteresis in discrete systems of possibly interacting elements with a double-well energy[END_REF][START_REF] Puglisi | Mechanics of a discrete chain with bi-stable elements[END_REF], we can rewrite [START_REF] Truskinovsky | About the origin of the nucleation peak in transformational plasticity[END_REF] in the form

Z p = 2π βμ N 2 α i exp ⎛ ⎝ -β μ 2N + N i=1 α i 2 ⎞ ⎠ .
In terms of the average strain = /N we obtain

Z p ( , β) = 2π βμ N 2 α i exp - Nβμ 2 ( + (2 p/N -1)) 2 . ( 23
)
The corresponding stress-strain relation reads

σ = μ -μ(1 -2 p/N ),
where p is a piece wise constant function of (see [START_REF] Fedelich | Hysteresis in discrete systems of possibly interacting elements with a double-well energy[END_REF][START_REF] Puglisi | Mechanics of a discrete chain with bi-stable elements[END_REF] for more detail).

Thermodynamic limit

To study the limiting case N → ∞, we introduce the new variable x = p/N and write

Z p = N -1/2 (2π) N -1 2 β N -1 2 1 x=0 x μ 1 + 1 -x μ 2 -1/2 × exp -N β 2 x μ 1 + 1 -x μ 2 -1 ( + xa 1 + (1 -x)a 2 ) 2 +x log(x) + (1 -x) log(1 -x) + x 2 log(μ 1 ) + 1 -x 2 log(μ 2 ) . ( 24
)
If we now denote

ψ(x) = β 2 x μ 1 + 1 -x μ 2 -1 ( + xa 1 + (1 -x)a 2 ) 2 + x log(x) + (1 -x) log(1 -x) + x 2 log(μ 1 ) + 1 -x 2 log(μ 2 ).
and replace the sum by an integral, we obtain

Z p ≈ N 1/2 (2π) N -1 2 β N -1 2 1 0 x μ 1 + 1 -x μ 2 -1/2 exp(-N ψ(x))dx. ( 25
)
Since N is large we can use Laplace method which gives

Z p = (2π) N 2 β N -1 2 1 ψ (x * ) 1 2 x * μ 1 + 1 -x * μ 2 -1/2 exp(-N ψ(x * )) ( 26
)
where x * is a minimum of ψ(x).

For large | |, the minimum of ψ(x) is either at x = 0 or at x = 1. More precisely, when ψ(x) is increasing we have

Z p ≈ N 1/2 (2π) N -1 2 β N -1 2 1 N ψ (0) μ 2 (N -1)/2 exp -N β 2 μ 2 ( + a 2 ) 2 , ( 27
)
while for decreasing ψ(x), we obtain

Z p ≈ N 1/2 (2π) N -1 2 β N -1 2 1 N ψ (1) μ 1 (N -1)/2 exp -N β 2 μ 1 ( + a 1 ) 2 .
Here ψ (0) and ψ (1) need to be determined by computing the corresponding derivatives with respect to (at constant β). One can see that we recover the thermodynamics of the single phases.

For the intermediate values of , we can write

σ = 1 β ∂ψ(x * ) ∂ | β = x * μ 1 + 1 -x * μ 2 -1 ( + x * a 1 + (1 -x * )a 2 ) + 1 β ∂ψ ∂ x * ∂ x * ∂ . Denote Ψ = ∂ψ ∂ x = μ 1 -μ 2 2μ 1 μ 2 β (x/μ 1 + (1 -x)/μ 2 ) 2 ( + xa 1 + (1 -x)a 2 ) 2 + 2 β (x/μ 1 + (1 -x)/μ 2 ) ( + xa 1 + (1 -x)a 2 ) + log(x) -log(1 -x) + 1/2 log(μ 1 ) -1/2 log(μ 2 ). (28) If x * ( ) is the solution of Ψ (x, ) = 0, we obtain σ = x * μ 1 + 1 -x * μ 2 -1 ( + x * a 1 + (1 -x * )a 2 ). ( 29
)
Using (29) we can rewrite equation for x * in the following way:

μ 1 -μ 2 2μ 1 μ 2 βσ 2 + 2βσ log(x * ) -log(1 -x * ) + 1/2 log(μ 1 ) -1/2 log(μ 2 ) = 0. ( 30
)
From here we obtain

x * = 1 1 + μ 1 μ 2 exp μ 1 -μ 2 2μ 1 μ 2 βσ 2 + 2βσ
.

This allows us to write the expression for the function (σ ) in the closed form

= σ ⎛ ⎜ ⎝ 1 μ 1 1 1 + μ 1 μ 2 exp μ 1 -μ 2 2μ 1 μ 2 βσ 2 + 2βσ + 1 μ 2 μ 1 μ 2 exp μ 1 -μ 2 2μ 1 μ 2 βσ 2 + 2βσ 1 + μ 1 μ 2 exp μ 1 -μ 2 2μ 1 μ 2 βσ 2 + 2βσ ⎞ ⎟ ⎠ - ⎛ ⎜ ⎝ 1 1 + μ 1 μ 2 exp μ 1 -μ 2 2μ 1 μ 2 βσ 2 + 2βσ a 1 + μ 1 μ 2 exp μ 1 -μ 2 2μ 1 μ 2 βσ 2 + 2βσ 1 + μ 1 μ 2 exp μ 1 -μ 2 2μ 1 μ 2 βσ 2 + 2βσ a 2 ⎞ ⎟ ⎠. (31)
In the case of symmetric double well energy with equal moduli μ 1 = μ 2 = 1 and α 1 = -1, α 2 = 1 we obtain

x * = 1 1 + exp(2βσ ) ,
and

σ = + 2 1 + exp(2βσ ) -1.
We now illustrate the obtained formulas. For the case of equal moduli, we take a 1 = 1 and a 2 = -1, and for the case of different moduli, we take a 1 = 1/ √ μ 1 and a 2 = -1/ √ μ 2 . In all cases, the discontinuity point for the stress-strain relationship is at = 0. In Fig. 8, we plot the isothermal stress-strain curves at different values of β. At large β (low temperatures), our results are close to those obtained through global energy minimization. At finite temperatures the isotherms exhibit finite slope in the transition region, which is a manifestation of "rubber elasticity". The finite slope is expected because in 1D the Ising model is known not to show phase transition in the thermodynamic limit.

Adiabats in the thermodynamic limit

The general expression for the entropy is given in [START_REF] Dreyer | Cold, thermal and oscillator closure of the atomic chain[END_REF]. For large N , we can write

log(Z p ) ≈ N 2 log(2π) - N -1 2 log(β) -N ψ(x * ),
and

β ∂ log(Z ) ∂β ≈ - N -1 2 -Nβ ∂ψ(x * ) ∂β .
Noting that ∂ψ(x * )/∂ x * = 0, we obtain Then the expression for the entropy takes the form

∂ψ(x * ) ∂β = 1 2 x * μ 1 + 1 -x * μ 2 -1 ( + x * a 1 + (1 -x * )a 2 ) 2 .
S = k N 1 2 + log(2π) -log(β) -x * log(x * ) -(1 -x * ) log(1 -x * ) - x * 2 log(μ 1 ) - (1 -x * ) 2 log(μ 2 ) . ( 32
)
This formula is illustrated in Fig. 9 where we plot the stress-strain curves σ ( , S) for different values of entropy. To this end, the equation S = const is first solved and the value of x * is expressed via β. This value is then substituted into the stress-strain relation.

To compare our thermodynamical theory with the results of the direct numerical simulations we would need to know the temperature-strain distribution along an adiabat β( , S). This relation is shown in Fig. 10 for different values of entropy. Notice that in the mixed-phase region the temperature decreases (see also [START_REF] Yoshikawa | Transition from mechanical equilibrium to thermal equilibrium in a chain of unstable springs[END_REF]), which can be explained by the high entropy of the corresponding two-phase inhomogeneous configuration. Indeed, one lattice "melts" before it becomes another lattice and it is then expected that the entropy of the "liquid" phase is above the entropies of the two coexisting solid phases.

A simple formula for the maximum temperature drop associated with this "solid to liquid" transition can be obtained in the case μ 1 = μ 2 . We observe that in this case x * = 0.5 at = 0 independently of β and write On the right plot one can easily check that for the limiting values

S = 1 2 + log(2π) -log(β t ) -x * log(x * ) -(1 -x * ) log(1 -x * ) = 1 2 + log(2π) -log(β t ) + log(2),
β 2 1 μ 1 = β 2 2 μ 2
where β t corresponds to the value of temperature at = 0. Once the chain reaches the state where all the springs are in one phase, then x * = 0 or x * = 1, and

S = 1 2 + log(2π) -log(β t ) -x * log(x * ) -(1 -x * ) log(1 -x * ) = 1 2 + log(2π) -log(β ± ),
where β ± characterize the temperature in the single phase states. Since the entropy of the chain does not change along the adiabat, we obtain β t /β ± = 2. In the general case of asymmetric wells, μ 1 = μ 2 , we can obtain β t /β ± = 1 1-x * (1/x * -1) x * , where x * must be computed at = 0.

Thermodynamic theory versus direct numerical simulations

We can now compare the thermodynamic theory with the results of our direct numerical simulations. We recall that the theory is based on the plausible but heuristic assumption that the motion induced by the irreversible self-heating is ergodic. We also recall that the obtained explicit formulas are approximate in the sense that we implicitly deal with a double-valued energy obtained by the continuation of each of the two parabolic wells beyond the spinodal point into the whole real axis. The last assumption can be easily dropped at the expense of analytical transparency, but we found that our simplified theory is sufficient to show the validity of the thermodynamic approach.

The first problem is to specify the specific adiabat β( , S) resulting from the irreversible self-heating of a 'cold' chain. As we have already noted earlier, this transition is irreversible and is accompanied by the entropy increase from S = 0 to, say, S = S 0 . To find the value of S 0 , we need to know the value of temperature at one value of . Observe that at the spinodal instability threshold the energy which transforms from macro to micro level is known exactly (see Fig. 1). More specifically, half of this (available) energy is transformed into the kinetic energy and half into the elastic energy. The kinetic energy per mode is equal to kT from where the temperature and therefore the value of β can be easily computed.

Our Fig. 12 and Fig. 11 compare the adiabats σ ( ) and β( ) obtained using our analytical model with the corresponding results of the direct numerical simulations for the quasi-statically driven chain with 30 springs. In simulations, the numerical results were smoothed using window averaging. We have tested various window sizes. The window size away from the mixed-phase region is taken to be five times larger than that in the mixed-phase region. To compute the stress-strain relation from the statistical model, we obtained the value of S 0 , generated associated adiabat β( , S 0 ) and substituted it into the stress strain relation σ (β, ).

The stress-strain curves presented in Fig. 11 show a very good agreement between the statistical model and the direct numerical simulations. The temperature-strain distributions β( , S 0 ) shown in Fig. 12 exhibit a slight mismatch in mixed-phase region, which we associate with the finite averaging window.

The overall agreement between the theory and the numerical experiment can be viewed as a strong argument in favor of the assumption of ergodicity achieved during spinodal instability. It remains to be understood at the formal level why the available potential energy at the instability point is exactly sufficient to thermalize the whole system.

Miscellaneous

In this section we discuss two issues tangentially related to the main subject of the paper. One concerns the dependence of the results on the type of the loading device and the other one concerns the role of small damping.

We begin with the comparison of the thermodynamical models for a chain with N elements loaded in soft and hard devices. Numerical experiments and statistical modeling for the case of soft device was done in [START_REF] Yoshikawa | Transition from mechanical equilibrium to thermal equilibrium in a chain of unstable springs[END_REF] for symmetric energy, and here we generalize these results for the case of asymmetric energy wells.

To analyze the chain in a soft device we need to modify the energy by adding the work of the loading device. The controlling parameter now is the applied stress σ and the configurational statistical sum takes the form One can see that in the case of soft device we deal with noninteracting elements (vs elements interacting through the mean field in the case of hard device). The statistical sum is then a product of the contributions due to independent elements, and we do not even need our simplifying assumption to obtain an analytic result

Z p ( σ , β) = i ∈R exp -β N i=1 μ i 2 ( i + α i ) 2 - σ N i=1 i d 1 , . . . , d N . ( 33
Z p ( σ , β) = √ π/2 √ βμ 1 exp( σ β) exp β σ 2 2μ 1 ERFC(-βμ 1 (1 + σ/μ 1 )) + √ π/2 √ βμ 2 exp(-σβ) exp βσ 2 2μ 2 (1 -ERFC(-βμ 2 (1 + σ/μ 2 ))) N , ( 34
)
where ERFC is the complementary error function. The stress-strain isotherms for the case of soft device are compared with their hard device analogs (computed using our analytical model) in Fig. 13. One can see that at small N and large β 1/μ the predictions of the two thermodynamical models can be rather different which is well known for the purely mechanical problem β → ∞ where in the soft device the system 'snaps' and in the hard device it 'pops' [START_REF] Puglisi | Mechanics of a discrete chain with bi-stable elements[END_REF]. As one can expect, in the thermodynamical limit the mode of loading becomes irrelevant.

The next issue which is appropriate to discuss in this Section concerns the role of small dissipation. We recall that the discrete system studied in the previous sections was Hamiltonian and therefore after initial 'self-heating' which erased the memory of the special preparation of the system, the effective macroscopic behavior under quasi-static cyclic loading was reversible. In many realistic systems with bi-or multi-stability at the micro-scale such reversibility, which can be associated with rubber-type behavior, is replaced by rate independent hysteresis.

To establish a link between the behavior of the undamped system and the overdamped behavior of a hysteretic system [START_REF] Puglisi | Thermodynamics of rate independent plasticity[END_REF] we consider a bi-stable FPU system with finite visco-elasticity at the micro-scale(different type of dissipation was introduced in [START_REF] Balk | Dynamics of chains with non-monotone stress-strain relations. I. Model and numerical experiments[END_REF][START_REF] Yoshikawa | Transition from mechanical equilibrium to thermal equilibrium in a chain of unstable springs[END_REF]). More precisely, we consider springs with viscosity D and write the system of governing equations in the form

ρh d 2 x i dt 2 -D d dt x i+1 -2x i + x i-1 h = Eσ x i+1 -x i h -Eσ x i -x i-1 h x 0 = 0, x N (t) = x 0 N + V t. ( 35
)
Next we introduce nondimensional variables

M = V 2 ρ/E, W = DV /(E L), δ = h/L = 1/N
d 2 x i dt 2 -W d dt x i+1 -2x i + x i-1 δ = σ x i+1 -x i δ -σ x i -x i-1 δ x 0 = 0, x N (t) = x 0 N + t. (36)
In the previous sections we studied the limit W = 0, M → 0, while the asymptotics considered in [START_REF] Puglisi | Thermodynamics of rate independent plasticity[END_REF] corresponded to the limit M = 0, W → 0. It is then natural to perform numerical simulations for the case when both M and W are small but finite. We assume that δ is finite, leaving for future studies the comprehensive analysis of the triple asymptotic structure for the singular point M = 0, W = 0, δ = 0. In Fig. 14 we show the time dependence of the kinetic energy for the initially 'cold' chain with 14 springs loaded quasi-statically in a hard device. We show one loading path (solid line) and one unloading path (dotted line). One can see that the elements start transforming from one phase to another at the spinodal point. The transition takes place in the form of avalanches accompanied by the emergence of oscillations. The oscillations quickly disappear due to dissipation before the arrival of the next avalanche. During unloading similar behavior of kinetic energy can be observed again with three major transformation events. In Fig. 15, we present the corresponding stress-strain diagram. Here also the solid line corresponds to elongation, and the dotted line corresponds to contraction. The overall behavior of the system is clearly hysteretic which is due to finite dissipation.

As the value of the parameter W increases, the individual avalanches get smaller and occur in shorter time intervals. The picture asymptotically approaches the one predicted for the overdamped system in [START_REF] Puglisi | Thermodynamics of rate independent plasticity[END_REF].

Conclusions

Our highly idealized model suggests that even in quasi-static conditions the homogenized description of the mechanical structures with internal instability (internal buckling) must be necessarily thermodynamical and that the canonical distribution may already be adequate for systems with relatively small number of degrees of freedom.

As we showed, in order to reach the thermodynamic level of description for our bi-stable FPU chain one needs to consider a double limit: δ → 0 (discrete to continuum) and v → 0 (quasi-static driving). The naive long-wave homogenization of the FPU chain suggests that the limiting problem will be of the form σ (ε) x = 0, (37) where σ (ε) = e (ε) is the reversible stress strain relation for the elastic material. A more careful line of reasoning going back to von Neumann [START_REF] Von Neumann | Proposal and analysis of a numerical method for the treatment of hydrodynamical shock problem[END_REF] suggests that the limiting equations will be of thermo-elasticity σ (ε, T ) x = 0, S(ε, T ) t = 0 (38)

where σ (ε, T ) = f ε (ε, T ) is the stress, S(ε, T ) =f T (ε, T ) is the entropy, T is the temperature and f (ε, T ) is the free energy. In this paper we derive an explicit relation for the f (ε, T ) for the bi-stable FPU chain. Our results suggest that the choice between (37) and (38) depends on the convexity of the function e(ε): if it is convex, we expect (37) to be true for the case of a quasi-static driving of initially 'cold' chain, while if it is non-convex, one eventually would have to use (38). Despite its simplicity our model shows both reversible heat exchange associated with release and absorbtion of the latent heat and irreversible heating due to energy dissipation. It is also a first example when a system can self-thermalize due to microscopic instability. One interesting extension of the model would be a chain with bi-stable nearest and linear ferromagnetic next to nearest neighbor interactions (e.g. [START_REF] Truskinovsky | About the origin of the nucleation peak in transformational plasticity[END_REF]). Such model is expected to exhibit more realistic transformation mechanism with localized phase boundaries separating the two phases. Introducing nonlinearity inside single phases will contribute to thermalization and show a rudiments of diffusive heat transfer. On the mathematical side, the formal justification of the use of canonical measure for the simplest bi-stable FPU model remains a challenge.

Fig. 1

 1 Fig.1Stress-strain and energy-strain relations for a bi-linear snap-spring. Available energy is released during spinodal instability

  Fig.2Total energy of the modes for the system of five springs during quasistatic stretching of an initially 'cold' chain

Fig. 3 Fig. 4

 34 Fig.3The motion of the chain with 80 particles after spinodal instability

Fig. 5

 5 Fig. 5 Total energy of linear modes versus time for the system with μ 1 = μ 2 = 1

Fig. 6 Fig. 7

 67 Fig. 6 Time dependence of the total energy of the linear modes for a 'self-heated' chain driven cyclically (and quasi-statically) through the phase transition. Left figure μ 1 = 9, μ 2 = 1. Right figure μ 1 = 0.25, μ 2 = 1

Fig. 8 Fig. 9

 89 Fig. 8 Stress-strain isotherms for different values of β. Left μ 1 = μ 2 . Right μ 2 = 2μ 1

Fig. 10

 10 Fig. 10 Temperature-strain relation along the adiabats. Left μ 1 = μ 2 . Right μ 2 = 2μ 1 . On the right plot one can easily check that for the limiting values β 2 1 μ 1 = β 2 2 μ 2

Fig. 11 Fig. 12

 1112 Fig. 11 Comparison of the stress-strain curves obtained in direct numerical simulations (dotted line) with the adiabats computed in the statistical model (solid line). σ versus . Left μ 1 = μ 2 . Right μ 2 = 2μ 1

)Fig. 13

 13 Fig. 13 Stress-strain isotherms for a chain with 30 springs driven in hard and soft devices. Left figure high temperature. Right figure low temperature

Fig. 14

 14 Fig. 14 Kinetic energy versus time for a damped chain during one cycle of loading and unloading

Fig. 15

 15 Fig. 15 Stress-strain relation for a damped chain during one cycle of loading and unloading

Acknowledgments

We would like to acknowledge useful discussions with Pierre Collet, Michael Ortiz and Henning Struchtrup.


Lei Pan

Clément Zanolli

Comparative observations on the premolar root and pulp canal configurations of Middle Pleistocene Homo in China

Objectives: The aim of this study is to explore the root and root canal morphology of Homo fossil occupying China during the Middle Pleistocene period. Human occupation and evolutionary dynamics in East Asia during the Middle Pleistocene period is one of the most intriguing issues in paleoanthropology, with the coexistence of multiple lineages and regional morphs suggesting a complex population interaction scenario. Although premolar root and canal morphology has certain phylogenetic, taxonomic, and functional implications, its morphological diversity, possible evolutionary trend and characteristics regarding Middle Pleistocene hominins inhabiting East Asia are still insufficiently understood; where these populations fits within the Homo lineage (with respect to root and pulp canal structure) needs to be explored. Materials and methods: Using microtomography, we directly observed and assessed the nonmetric variability of root and canal forms in maxillary and mandibular premolars of Chinese Middle Pleistocene Homo (N = 19), and compared our observed variations with Eurasian Early Pleistocene specimens from the Asia continent (N = 1) and Java (N = 2), as well as with Neanderthals (N = 28) and recent modern humans (N = 67). Results: A total number of nine types of root-canal forms were recorded. As a whole, the Chinese Middle Pleistocene record shows an evolutionary trend toward a modern human-like condition (a reduction of root/canal number and a simplification of root surface structure). We documented primitive signals like high percentage of Tomes' root in lower premolars. A considerable occurrence of incompletely separated root branches and bifid root and canal apices, representing evolutionary transformation from multi-root to single-root condition was also noticed. The results were compared with previous publications on Early and Middle Pleistocene Homo in East Africa, North Africa, and Eurasia. Conclusion: This work provides new original data, incorporates the latest human fossil discoveries and suggests that analyzing the variation of premolar root structural organization, notably integrating together root/canal form and number, could possibly contribute to taxonomic and phylogenetic assessments. The mid-Middle Pleistocene populations, or "classic" Homo erectus, in our study show closer affinity to Early and Middle Pleistocene hominins in Eurasia, than to East African early Homo, which supports the suggestion that at least some of the Early Pleistocene hominin groups in Eurasia contribute to the later population; on the other hand, it is still difficult to clearly trace the evolutionary fate of those late Middle Pleistocene populations (roughly assigned as archaic Homo sapiens through a craniodental perspective). More comparable materials from the Early to Middle Pleistocene period as well as precise chronological framework is needed to further explore the evolutionary trends of archaic hominins in the Asian continent before the arrival of modern humans.

Introduction

Recent discoveries revealed different evolutionary lineages and a considerable degree of isolation among regional populations in Middle Pleistocene Asia (rev. in [START_REF] Kaifu | Archaic hominin populations in Asia before the arrival of modern humans: their phylogeny and implications for the "Southern Denisovans[END_REF].

Marked differences were reported for the internal tooth crown morphology of Chinese and Indonesian Middle Pleistocene humans, possibly suggesting that distinct groups may have populated continental and (pen)insular Asia [START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF][START_REF] Liu | Geometric morphometric analysis of the early Pleistocene hominin teeth from Jianshi, Hubei Province, China[END_REF][START_REF] Liu | Morphological diversities of the late-Middle and Late Pleistocene human fossils in China (in Chinese)[END_REF][START_REF] Xing | Morphological variation of Zhoukoudian H. erectus teeth[END_REF][START_REF] Xing | Hominin teeth from the Middle Pleistocene site of Yiyuan, Eastern China[END_REF]. In particular, complex occlusal crown surface and enamel-dentine junction morphology (represented by additive traits like interconnected ridges, bifurcated essential crests, and accessory ridges and cusps) is found in the Chinese dental assemblages across the Middle Pleistocene, including the earliest specimens from Zhoukoudian [START_REF] Xing | The fossil teeth of the Peking Man[END_REF], the mid-Middle Pleistocene teeth from Chaoxian, Hexian, and Yiyuan [START_REF] Bailey | A comparative dental metrical and morphological analysis of a Middle Pleistocene hominin maxilla from Chaoxian (Chaohu), China[END_REF][START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF][START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF][START_REF] Xing | Hominin teeth from the Middle Pleistocene site of Yiyuan, Eastern China[END_REF], and the late Middle to early Late Pleistocene material from Panxian Dadong and Xujiayao [START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Liu | Morphological diversities of the late-Middle and Late Pleistocene human fossils in China (in Chinese)[END_REF]Xing et al., 2015).

These findings revealed characteristic patterns for the populations that inhabited East Asia during the Middle Pleistocene period, like highly crenulated enamel-dentine junction (EDJ) and its imprint on the roof of the pulp cavity [START_REF] Xing | Hominin teeth from the Middle Pleistocene site of Yiyuan, Eastern China[END_REF][START_REF] Xing | The fossil teeth of the Peking Man[END_REF], and highlighted a combination of primitive and derived features that cannot be organized along a temporal scale [START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF]. Based on dental evidence, those early-and mid-Middle Pleistocene specimens are aligned to the "classic" Homo erectus hypodigm and it is proposed that Early Pleistocene hominins in Eurasia contributed to at least some of the later H. erectus lineages, and that these populations experienced relatively independent evolutionary courses compared with H. erectus in the African continent [START_REF] Kaifu | Archaic hominin populations in Asia before the arrival of modern humans: their phylogeny and implications for the "Southern Denisovans[END_REF][START_REF] Kaifu | Taxonomic affinities and evolutionary history of the Early Pleistocene hominids of Java: dentognathic evidence[END_REF]Martinón-Torres et al., 2007[START_REF] Peiris | Root and canal morphology of human permanent teeth in a Sri Lankan and Japanese population[END_REF]. On the other hand, craniodental evidence is still insufficient to resolve the taxonomic and phylogenetic identity of late Middle Pleistocene hominins and how they relate to modern Homo sapiens. Except for a retention of primitive traits, like large and robust roots, and a trend toward root number reduction [START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF][START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF][START_REF] Xing | The fossil teeth of the Peking Man[END_REF], the morphological diversity and phylogenetic implications of the root-canal structure of these hominin populations are still unclear; comparative study is needed to fully reveal the root and canal structure for the Chinese materials from this period. The variability of form, number and position of fossil hominin and extant human premolar roots and/or pulp canals has been described in previous research as reflecting functional/biomechanical demand, developmental aspects, sexual dimorphism, and taxonomy [START_REF] Abbott | A comparative study of tooth root morphology in the great apes, modern man and early hominids[END_REF]Kottoor, Albuquerque, Velmurugan, & Kuruvilla, 2013;[START_REF] Kupczik | Comparative observations on the tooth root morphology of Gigantopithecus blacki[END_REF]Kupczik, Spoor, Pommert, & Dean, 2005;[START_REF] Moore | Premolar root and canal variation in extant non-human hominoidea[END_REF][START_REF] Moore | Premolar root morphology and metric variation in Pan troglodytes verus[END_REF]Moore, Thackeray, Hublin, & Skinner, 2016;[START_REF] Spencer | Tooth-root form and function in platyrrhine seed-eaters[END_REF][START_REF] Wood | Analysis of the dental morphology of Plio-Pleistocene hominids. V. Maxillary postcanine tooth morphology[END_REF]. Root morphology and pulp canal form, number, and size are the result of variations in the development of Hertwig's epithelial root sheath that determines the location, number and expression of inter-radicular processes (Kovacs, 1967(Kovacs, , 1971;;[START_REF] Shields | Mandibular premolar and second molar root morphological variation in modern humans: what root number can tell us about tooth morphogenesis[END_REF][START_REF] Moore | Premolar root morphology and metric variation in Pan troglodytes verus[END_REF]. Hominid root complexity (single vs. multiple roots) and length are generally interpreted to play a role in the dispersion of occlusal loading during mastication and paramasticatory activities, but root number alone does not reflect a functional adaptation related to area of root attachment [START_REF] Kupczik | Premolar root number variation in hominoids: genetic polymorphism vs. functional significance[END_REF]. Root and canal configurations are also associated with complex genetic, and epigenetic mechanisms as well as environmental influences [START_REF] Brook | Multilevel complex interactions between genetic, epigenetic and environmental factors in the aetiology of anomalies of dental development[END_REF][START_REF] Shields | Mandibular premolar and second molar root morphological variation in modern humans: what root number can tell us about tooth morphogenesis[END_REF][START_REF] Tobias | Root number in the maxillary third premolars: a very ancient polymorphism[END_REF]. But tooth root structural organization also holds critical taxonomic information and represents a phylogenetic indicator for distinguishing closely related extinct hominid species [START_REF] Moore | Premolar root morphology and metric variation in Pan troglodytes verus[END_REF][START_REF] Moore | Premolar root and canal variation in South African Plio-Pleistocene specimens attributed to Australopithecus africanus and Paranthropus robustus[END_REF]. Within the Pan-modern human lineage, while variations exist, the plesiomorphic great-ape pattern is represented by (a) three-rooted maxillary premolars with three canals, arranged as two buccal and one lingual root, and (b) double-rooted mandibular premolars with three (P3) or four (P4) canals [START_REF] Abbott | A comparative study of tooth root morphology in the great apes, modern man and early hominids[END_REF]. In contrast, modern humans show a series of apomorphic condition: modern human upper premolars are usually single-rooted [START_REF] Abbott | A comparative study of tooth root morphology in the great apes, modern man and early hominids[END_REF][START_REF] Peiris | Root and canal morphology of human permanent teeth in a Sri Lankan and Japanese population[END_REF][START_REF] Wood | Analysis of the dental morphology of Plio-Pleistocene hominids. V. Maxillary postcanine tooth morphology[END_REF], but in certain populations, predominant doublerooted P3s have been recorded [START_REF] Shaw | The teeth, the bony palate and the mandible in Bantu races of South Africa[END_REF], and sub-Saharan populations even show a wider range of variability [START_REF] Shields | Mandibular premolar and second molar root morphological variation in modern humans: what root number can tell us about tooth morphogenesis[END_REF]. Regarding lower premolars, the modern human condition is a single root, although variations such as Tomes' root have been reported (see [START_REF] Loh | Root morphology of the maxillary first premolar in Singaporeans[END_REF][START_REF] Peiris | Root and canal morphology of human permanent teeth in a Sri Lankan and Japanese population[END_REF][START_REF] Sert | Evaluation of the root canal configurations of the mandibular and maxillary permanent teeth by gender in the Turkish population[END_REF][START_REF] Shields | Mandibular premolar and second molar root morphological variation in modern humans: what root number can tell us about tooth morphogenesis[END_REF][START_REF] Tomes | A manual of dental anatomy: human and comparative[END_REF]. Among fossil taxa, elaborated or "molariform" lower premolar roots (each plate-like and contains two canals) are seen in Australopithecus and Paranthropus (Moggi-Cecchi, Menter, Boccone, & Keyser, 2010;[START_REF] Moore | Premolar root and canal variation in South African Plio-Pleistocene specimens attributed to Australopithecus africanus and Paranthropus robustus[END_REF]White, Suwa, Simpson, & Asfaw, 2000;[START_REF] Wood | Analysis of the dental morphology of Plio-Pleistocene hominids. IV. Mandibular postcanine root morphology[END_REF]. Regarding Homo, an apomorphic condition with respect to australopiths and a trend towards the reduction in root radicule and pulp canal number is expressed. Although doublerooted premolars exist, the East African Homo exhibits predominantly single or Tomes' root conditions, whereas primitive double-rooted lower premolars can be observed in some other members of Early Pleistocene H. erectus s.l. from Dmanisi (D2600) and Sangiran (Sangiran 5, 8, and 9; [START_REF] Kaifu | Taxonomic affinities and evolutionary history of the Early Pleistocene hominids of Java: dentognathic evidence[END_REF][START_REF] Martinón-Torres | Dental remains from Dmanisi (Republic of Georgia): morphological analysis and comparative study[END_REF], as well as in H. antecessor and H. naledi (Berger et al., 2015;Bermúdez de Castro, Martinón-Torres, Sier, & Martín-Francés, 2014;[START_REF] Bermúdez De Castro | Dental remains from Atapuerca-TD6 (Gran Dolina site, Burgos, Spain)[END_REF]. In addition, the mandibular premolars of H. naledi are somewhat molarized (Berger et al., 2015).

Goals of this study

To date, root and pulp canal structure has been reported for a number of archaic Homo in East Asia [START_REF] Chang | The first archaic Homo from Taiwan[END_REF][START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF][START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF][START_REF] Xing | The fossil teeth of the Peking Man[END_REF]Zanolli et al., 2018); previous investigations [START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Xing | The fossil teeth of the Peking Man[END_REF] have revealed a general trend toward root number reduction, but a retention of primitive traits is still represented in particular cases like Hexian H. erectus [START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF][START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF] and Penghu 1 mandible [START_REF] Chang | The first archaic Homo from Taiwan[END_REF]. However, the root and pulp canal structure of Chinese Middle Pleistocene hominins has not been sufficiently recorded (especially late Middle Pleistocene specimens that are generally assigned to archaic H. sapiens). Moreover, comparative investigations focusing on the root and canal morphology are needed, to better understand where these populations fit along the Homo lineage. We incorporated the latest discoveries and characterized here the extent of root number and pulp canal variation in an extensive sample of Chinese Middle Pleistocene premolars. More specifically, we used micro-CT scans of the majority of the Chinese Middle Pleistocene Homo specimens suitable for rootcanal investigation to document and assess the variability of premolar rootcanal system. We compared results on this Chinese fossil sample with Early Pleistocene Homo from the Asian continent and Java, as well as with Late Pleistocene Neanderthals and recent modern humans.

Materials and method

Study sample

The samples comprise a total of 117 permanent maxillary and mandibular premolars from fossil and extant individuals from a variety of chronological and geographical origins. The Chinese Middle Pleistocene Homo sample (CMPH) includes 19 teeth, with a chronological age ranging from mid-Middle to late-Middle Pleistocene, scattering amongst 8 archaeological localities. It should be noted that the specimens from Xichuan lack accurate stratigraphic context, but their morphological signature could be compared to Zhoukoudian H. erectus and therefore it is acknowledged that Xichuan hominins are H. erectus [START_REF] Wu | Human fossil teeth from Xichuan, Henan (In Chinese with English abstract)[END_REF]. We compared the results of Chinese Middle Pleistocene Homo with Neanderthals (NEA, N=28), recent modern humans (RMH, N=67) and two Early Pleistocene H. erectus individuals from East Asian continent and Indonesia. Our modern human reference material comes from clinical extractions and anatomical collections, sampling Asian, European and South African individuals; our Neanderthal sample comes from Krapina site, in Croatia, and La Chaise, France (NESPOS database, 2017) (SOM Table 1). We selected specimens with well-preserved roots, when antimeres exist, the better-preserved side is chosen; when both antimeres are well preserved, the left side is chosen (regarding our sample, well-preserved antimeres exhibit the same root pattern).

Except for PMU M 3887 housed at Uppsala University (which was scanned using X-ray microtomography (µCT) at the Multidisciplinary Laboratory of the ICTP, Trieste [START_REF] Tuniz | The ICTP-Elettra X-ray laboratory for cultural heritage and archaeology. A facility for training and education in the developing world[END_REF]Zanolli et al., 2018a), all other Chinese specimens (fossil and recent humans) were scanned using a 225 kV-µXCT scanner at the Institute of Vertebrate Paleontology and Paleoanthropology (IVPP, Chinese Academy of Sciences); micro-CT data or STL mesh of the Neanderthals were downloaded from the NESPOS Database (NESPOS database, 2017); data for South African modern human specimens were produced on an X-Tek (Metris) XT H225L industrial XCT system at the South African Nuclear Energy Corporation (Necsa) [START_REF] Hoffman | Characteristics of the micro-focus X-ray tomography facility[END_REF]; European modern human specimens were acquired using a Phoenix Nanotom 180 scanner from the FERMAT Federation from the Inter-university Material Research and Engineering Centre (CIRIMAT, UMR 5085 CNRS).

Isometric voxel size ranged from 40 to 70 µm. Each specimen was segmented in Avizo 8.0 (Visualization Sciences Group, www.vsg3d.com). After the segmentation, the root and canal system was generated using the "unconstrained smoothing" parameter. (Wood et al., 1988): 1R, 2T, 2R MB + D (mesiobuccal and distal) and 2R M + D (mesial and distal). 1R was defined as a single root with a single (main) canal (single-rooted Tomes' root with one apex is included in this category); 2T was defined as the Tomes' root or C-shaped root with bifid apex [START_REF] Tomes | A manual of dental anatomy: human and comparative[END_REF], where each root has a separate pulp canal. This classification system was later extended by [START_REF] Moore | Premolar root morphology and metric variation in Pan troglodytes verus[END_REF][START_REF] Moore | Premolar root and canal variation in extant non-human hominoidea[END_REF][START_REF] Moore | Premolar root and canal variation in South African Plio-Pleistocene specimens attributed to Australopithecus africanus and Paranthropus robustus[END_REF] to combine both the external root traits and internal canal parameters.

The classification typology we used here was adapted from formula set by [START_REF] Moore | Premolar root morphology and metric variation in Pan troglodytes verus[END_REF][START_REF] Moore | Premolar root and canal variation in extant non-human hominoidea[END_REF][START_REF] Moore | Premolar root and canal variation in South African Plio-Pleistocene specimens attributed to Australopithecus africanus and Paranthropus robustus[END_REF], slight modifications were made in order to highlight number of roots.

The formula and description of each root/canal type is listed in Table 2. Multi-rooted premolars were defined on the basis of furcation exceeding more than 25% of the total root length, from cervix to root apex, and the classification of Tomes' root was done following ASUDAS [START_REF] Turner | Scoring procedures for key morphological traits of the permanent dentition: the Arizona State University dental anthropology system[END_REF]. The evaluation of root and canal types were performed three times by the same observer (LP).

Results

Among the examined groups, we identified a total of nine types of root and pulp canal configurations (Table 2). Upper premolars show single, double or triple roots with one, two or three canals; lower premolars are single-or Tomes'rooted and have one or two canals.

The Figures 1 and2 show the percentage of root-canal type in each group with reference to dental position, as well as a representative cross-sectional image of each root-canal type, while the Figure 3 shows the percentage of premolar root/canal type in the three groups investigated in the present study. Root type of each fossil specimen was listed in Supplementary Online Material (SOM) Table 1.

Variation of root and canal types among dental position and groups

Maxillary premolars

Two-thirds of Chinese fossil P 3 s shows double roots with two canals (Figs. 1 and3), one specimen (PA 832) has triple roots with two buccal and one lingual canals (SOM Table 1).

The root configuration of this specimen has been previously reported [START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF].

Zhoukoudian specimen PA 67 shows incompletely separated B+L roots, connected by cementum; P 3 specimens sh.y.003 as well as Chaoxian have two independent, B+L roots, and PA 76 has single root with B+L canals (SOM Table 1). The Early Pleistocene H.

erectus representatives Sangiran 4 and PA 105-4 also show two separate roots with a pulp canal in each branch (SOM Table 1), because of a small sample size, these early specimens were not introduced in the frequency counts or in Figures 1, 2 and 3, but micro-CT sections are presented in SOM. The frequency of double-rooted P 3 s is also high in recent modern humans and Neanderthals (68.8% and 100%, respectively) (Figs. 1 and3). Modern human P 3 s present two cases of single joined root with two canals, and one specimen has apical bifurcation (Fig. 1).

P 4 s vary in root number between one and two. One half of Chinese Middle Pleistocene sample show double roots and two canals (Fig. 1). Specimen PA 68 has a bifid root apex with a bifid canal (SOM Table 1); sh.y.007 has single, plate-like root with two canals; double roots and two B+L canals form is expressed for PA 524 (SOM Table 1), the buccal and lingual components are separated by a deep groove at the mesial side (n.b., this root is damaged, but a sign of two branches can be seen). The P 4 of the Javanese specimen Sangiran 4 shows a deep furrow at the distal face of the root with two pulp canals almost linked together by thin buccolingual elongations, similar as the configuration of Chaoxian UP4 (SOM Table 1). Double-rooted P 4 s is observed in modern humans (38.5%), and more frequently, in Neanderthals (87.5%) (Figs 1 and3).

Mandibular premolars

Across the sample, root type varies between single circular, bifid, and Tomes' root, canal number varies between one and two (Figs. 2 and3). A key difference between Chinese Mid-Pleistocene fossil and the other two groups is the relatively high prevalence of Tomes' root (Figs. 2 and3). In P3s, a two-canaled Tomes' type with C-shaped buccal face and deep MB-DL groove is expressed for Chinese fossils (60%), followed by single circular root with one canal (40%); in P4s, 2 fossil specimens were discovered with Tomes' root (50%).

The two-canal, Tomes' rooted P4 is reported in single case (PA 831), at the external surface of the root (SOM Table 1), a mesiobuccal and a distolingual radical grooves are separated by a deep lingual invagination [START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF]. Note that only the apical third of this P4 is preserved, embedded in the mandibular corpus. Three moderate radicular grooves (M, D and L) are seen at the root surface of PA 81, together with a crescent-shaped canal in cross-section, these characteristics thusly differentiate this root from a typical single-rooted premolar, and categorize itself as a non-typical Tomes' root (SOM Table 1). In modern human lower premolars, the prevalence of Tomes' root is 27.8% for P3 and 15% for P4 (Figs. 2 and3). It is also interesting to note that compared with the other two groups, Neanderthals P3s show a relatively high percentage of single root with bifurcation at the root and canal apex (37.5%) (Figs. 2 and3).

Additional morphological traits and variants within root/canal formula

As noticed by [START_REF] Moore | Premolar root and canal variation in extant non-human hominoidea[END_REF][START_REF] Moore | Premolar root and canal variation in South African Plio-Pleistocene specimens attributed to Australopithecus africanus and Paranthropus robustus[END_REF], double-rooted upper premolars (referred as formula 11B + 11 L or root type 2-A-2 by Moore and colleagues) have a number of variants within hominoids, and this is in accordance with our observations. In P 3 of Zhoukoudian H. erectus, two root branches are presented but not completely separated: cementum tissue connects root branches at the mesial face, and a deep and wide furrow is seen at the distal face, located at the middle third and extending to the apex of the root, with two, clearly separated B + L canals. This root configuration is also seen in the P4s of Javanese H.

erectus Sangiran 4, Chaoxian P4, and also in some of the Neanderthal P4s (Supporting Information Table 1), but is absent in early H. erectus specimen PA 105-4, not frequently observed in our modern humans.

Discussion

Inter-group comparisons of root and canal form

Recent studies on the dentition of Chinese hominin fossils from Middle Pleistocene period (including a number of fossils involved in this study) have reported a combination of derived and primitive traits expressed at premolar roots and/or root canals. Compared with early members of H. erectus s.l. (e.g., specimens from Dmanisi, East Africa, and Sangiran), upper premolars from Zhoukoudian and Yiyuan have less complex root and canal structures, whereas the lower premolars are closer to the early Homo condition [START_REF] Martinón-Torres | Dental remains from Dmanisi (Republic of Georgia): morphological analysis and comparative study[END_REF][START_REF] Xing | The fossil teeth of the Peking Man[END_REF][START_REF] Zanolli | Inner tooth morphology of Homo erectus from Zhoukoudian. New evidence from an old collection housed at Uppsala University, Sweden[END_REF]. On the other hand, Hexian H. erectus differs from northern Chinese H. erectus and European Middle Pleistocene hominins by showing primitive affinities (three-rooted upper premolar, robust mandible and premolar roots) with the Early Pleistocene specimens from Africa (H. ergaster) and Java (H. erectus; [START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF][START_REF] Xing | Middle Pleistocene hominin teeth from Longtan Cave, Hexian, China[END_REF]. Another study states that P3 from Panxian Dadong shows a slightly robust root compared with modern H. sapiens, but the archaic features are very weakly expressed [START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF]; our observation on the root and canal form confirms that this late Middle Pleistocene hominin specimen has a simple root system (single root with one canal), a condition that is commonly found in recent humans. Besides these above-mentioned, previous records, the root and canal system of Middle Pleistocene hominins from China is still insufficiently documented: little is known about their morphological variability, leaving inconclusive discussions about their taxonomical and phylogenetic places in relation to contemporary hominins from Africa and Eurasia. We extend the knowledge on fossil human tooth variability by providing the first detailed document of the root-pulp canal system variation of Asian hominins living in this particular time period. The primitive hominin root morphology, as represented by the maxillary premolars of South African australopiths, are double-or triple-rooted [START_REF] Moore | Premolar root and canal variation in South African Plio-Pleistocene specimens attributed to Australopithecus africanus and Paranthropus robustus[END_REF]. Triple-rooted P 3 s are also observed in Early Pleistocene Homo from Sangiran (older deposits; [START_REF] Kaifu | Taxonomic affinities and evolutionary history of the Early Pleistocene hominids of Java: dentognathic evidence[END_REF][START_REF] Von Koenigswald | Fossil hominids from the Lower Pleistocene of Java[END_REF]. Contrasting this "molarized" primitive condition, our representatives of Asian Early Pleistocene Homo (Sangiran 4 and PA 105-4), as well as Chinese Middle Pleistocene sample, show a clear trend of root reduction (toward modern human condition). Chinese specimens in general exhibit a derived condition of predominantly double-rooted P 3 s and single-rooted P 4 s (see also [START_REF] Xing | The fossil teeth of the Peking Man[END_REF], close to the frequency seen in modern humans (Figure 1), although variations like triple-and single-rooted P 3 s do exist. It has to be noted that [START_REF] Xing | The fossil teeth of the Peking Man[END_REF] considered Zhoukoudian H. erectus specimen PA 67 as three-rooted because of the presence of a small and unattached pulp diverticulum between the two main canals [START_REF] Xing | The fossil teeth of the Peking Man[END_REF] Figure 2). However, because of its position in the inter-radicular blade and flatten shape, we consider this to represent a side-effect of the partial fusing of the two radicles and not a third root canal (see Supporting Information Table 1). If we enlarge our sample to other contemporaneous Zhoukoudian H. erectus materials (which were lost during WWII), [START_REF] Weidenreich | The dentition of Sinanthropus pekinensis: a comparative odontography of the hominids[END_REF] described a few double-rooted premolars with incompletely separated root branches; in fact, 2/3 of P 3 s and 2/2 of P 4 s recorded by [START_REF] Weidenreich | The dentition of Sinanthropus pekinensis: a comparative odontography of the hominids[END_REF] show this feature. Partly divided upper premolar roots are also frequently found in Neanderthals (this study) and is seen in one Middle Pleistocene H. floresiensis-like specimen (SOAMM5) from Mata Menge, Flores (van den [START_REF] Van Den Bergh | Homo floresiensis-like fossils from the early Middle Pleistocene of Flores[END_REF]. While the premolar root morphology of the early modern human assemblage from Jebel Irhoud, in Morocco, is still unreported (Hublin et al., 2017), the recently published late Middle Pleistocene modern human maxilla from Misliya, in Israel, showed singlerooted upper premolars [START_REF] Hershkovitz | The earliest modern humans outside Africa[END_REF]. Double-rooted lower premolars are frequently reported for Early Pleistocene hominins found in Eurasian sites like Atapureca TD6 (Bermúdez de [START_REF] Bermúdez De Castro | Dental remains from Atapuerca-TD6 (Gran Dolina site, Burgos, Spain)[END_REF] and Sangiran [START_REF] Kaifu | Taxonomic affinities and evolutionary history of the Early Pleistocene hominids of Java: dentognathic evidence[END_REF], but East African early Homo has predominantly single-or Tomes' rooted lower premolars (Wood et al., 1988) [START_REF] Chang | The first archaic Homo from Taiwan[END_REF]. A three-grooved P4 root (with one circular canal) is seen in late Middle Pleistocene specimen from Changyang (this study), representing a fusion of three root branches-a rather primitive feature. In all, our Middle Pleistocene record from China fits this evolutionary scenario, with a retention of primitive and derived features associated with root and canal simplification. Single-rooted lower premolars are also frequently reported for the Middle Pleistocene hominins from Eurasia, including specimens from Zhoukoudian [START_REF] Weidenreich | The dentition of Sinanthropus pekinensis: a comparative odontography of the hominids[END_REF][START_REF] Xing | The fossil teeth of the Peking Man[END_REF][START_REF] Zanolli | Inner tooth morphology of Homo erectus from Zhoukoudian. New evidence from an old collection housed at Uppsala University, Sweden[END_REF] and Chenjiawo (this study), the P4 of Penghu 1 [START_REF] Chang | The first archaic Homo from Taiwan[END_REF], Qesem Cave (Hershkovitz et al., 2011), Atapuerca Sima de los Huesos (Martinón-Torres et al., 2012), and Visogliano (Zanolli et al., 2018).

In the Late Pleistocene, Neanderthal premolars shows more homogeneity in their rootcanal form. All the P3s have double roots and all the P3s have single root. Moreover, all the Neanderthal P4s show single root (1R) with one or two canals except for one specimen (KRD 35) that exhibits Tomes' root with deep MB groove and two MB + DL canals (Supporting Information Table 1). These results based on a larger sample are comparable to those from Prado-Simón, Martinón-Torres, Baca, Olejniczak, et al. (2012). However, the majority of our Neanderthal sample comes from Krapina and further investigation is needed to test whether this low variability on the premolar rootcanal form exists in other Neanderthal populations. The modern human lower premolars show derived root forms of predominantly, if not all, single root with one canal. Despite the tendency towards a simplification of the root system across the Middle Pleistocene in Eurasian human groups, the ancestral pattern is found to be retained in multiple human groups through time and space, probably indicating the survival of primitive lineages, but more evidence is needed to make any phylogenetic interpretations. In all, premolar root and canal form of mid-Middle Pleistocene hominins inhabiting China (or "classic" H. erectus) show closer affinities with their contemporaneous or older counterparts in Eurasia, than with East African early Homo. This evidence is consistent with the hypothesis that at least some of the Early Pleistocene hominins occupying Eurasia contributed to the later "classic" H. erectus population [START_REF] Kaifu | Taxonomic affinities and evolutionary history of the Early Pleistocene hominids of Java: dentognathic evidence[END_REF]Martinón-Torres et al., 2007[START_REF] Peiris | Root and canal morphology of human permanent teeth in a Sri Lankan and Japanese population[END_REF]. As for the late Middle Pleistocene specimens, although they are often assigned to archaic H. sapiens and display a series of derived features, we record here a number of primitive features that cannot be strictly aligned across a chronological sequence. Together with root-canal configurations reported here, current detal evidence cannot confidently link these Middle Pleistocene population to the H. sapiens lineage, although an evolutionary trend toward this condition is clearly observed [START_REF] Bailey | A comparative dental metrical and morphological analysis of a Middle Pleistocene hominin maxilla from Chaoxian (Chaohu), China[END_REF][START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF].

Conclusions

Previous studies have reported a high degree of morphological complexity in tooth crown and enamel-dental junction [START_REF] Liu | Late Middle Pleistocene hominin teeth from Panxian Dadong, South China[END_REF][START_REF] Liu | Morphological diversities of the late-Middle and Late Pleistocene human fossils in China (in Chinese)[END_REF][START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF][START_REF] Xing | Hominin teeth from the Middle Pleistocene site of Yiyuan, Eastern China[END_REF][START_REF] Xing | Crown morphology and variation of the lower premolars of Zhoukoudian Homo erectus[END_REF], compatible with our results on the root and pulp canal. Together with a clear trend of root number reduction during the Middle and Late Pleistocene, our study revealed a previously undocumented level of diversity in root-canal configuration within the hominins that occupied East Asia in Middle Pleistocene period. In all, premolar root and canal form we observed in "classic" H. erectus population is comparable to Early and Middle Pleistocene hominins in Eurasia, but show distant morphological affinities to East African early hominins. This work shows that analyzing the variation of premolar root structural organization, and notably integrating together root/canal form and number, could possibly contribute to taxonomic and phylogenetic discussions. Future investigations is needed, for example, by characterizing the root-canal morphology of the Indonesian Early-Middle Pleistocene hypodigm in order to better understand the still obscure evolutionary relationships and population dynamics with the Asian continental groups.

  : (A) Huang et al. (1993) ; (B) Grün et al. (1997) ; (C) Shen et al. (2001) ; (D) Grün et al. (1998) ; (E) Shen et al. (2010); (F) Yuan et al. (1986) ; (G) Han et al. (2016) ; (H) An and Kun (1989) ; (I) Schepartz and Miller-Antonio (2004); (J) Zhu et al. (2015) ; (K) Huffman (2001) ; (L) Rink et al. (1995) ; (M)[START_REF] Macchiarelli | How Neanderthal molar teeth grew[END_REF]

Figure 1 .

 1 Figure 1. The frequency of root-canal types in upper premolar. Examples of each root type in cross-section image are shown.

Figure 2 .

 2 Figure 2. The frequency of root-canal types in upper premolar. Examples of each root type in cross-section image are shown.

Figure 3 .

 3 Figure 3. Root and canal type in Chinese Middle Pleistocene Homo (CMPH), Neanderthals (NEA) and recent modern humans (MH).

Table 1 .

 1 Composition of the study sample.

Table 2 .

 2 List of premolar root and canal formulae identified in our sample

	Formula	Root	Canal	Description
		number	number
	1R1	1	1	No bifurcation; 1 root/1 canal
	1R(Bf)1Bf	1	1	1 bifid root/1 bifid canal
	1R1Bf	1	1	1 single root/1 bifid canal
	1R2(1B+1L)	1	2	1 root (joined 1B+1L roots)/2 circular canals (1B+1L)
	1T1	1	1	Tomes' root; 1 C-shaped root body (deep ML
				groove, round DB face)/1 C-shaped canal
	1T2(1B+1L)	1	2	Tomes' root; 1 C-shaped root body/2 circular canals
				(1B+1L)
	1T2(1MB+1DL)	1	2	Tomes' root; 1 C-shaped root body /2 circular canals
				(1MB+1DL)
	2R: 11B+11L 2	2	2 roots (1B+1L circular)/2 circular canals (1B+1L).
				Double roots but incompletely separated, fully-grown
				root branches are included
	3R: 22B+11L 3	3	3 roots (2B+1L)/3 circular canals (2B+1L)

  , high percentage of Tomes' root is also observed in Early Pleistocene specimens from Atapuerca TD6 and Sima del Elefante[START_REF] Prado-Simón | Three-dimensional evaluation of root canal morphology in lower second premolars of early and middle pleistocene human populations from atapuerca (Burgos, Spain)[END_REF],Dmanisi and Trinil (Martinón-Torres et al., 2008). As for the Middle Pleistocene period, Tomes' root has been described for the specimens found in Tighenif (accompanied with two or three canals;[START_REF] Zanolli | Endostructural characterization of the H. heidelbergensis dental remains from the early Middle Pleistocene site of Tighenif, Algeria[END_REF], mid-Middle Pleistocene specimens like H. erectus P4 from Hexian[START_REF] Liu | A mandible from the Middle Pleistocene Hexian site and its significance in relation to the variability of Asian Homo erectus[END_REF] and Zhoukoudian (3/5 P3s and 1/2 P4s;[START_REF] Weidenreich | The dentition of Sinanthropus pekinensis: a comparative odontography of the hominids[END_REF], P3 from Chenjiawo (this study) and a late Middle Pleistocene P3 of Penghu 1 (

Acknowledgements

For access to comparative materials we are grateful to P. Bayle, G. Kruger, O. Kullmer, E. L'Abbé, D. S. Li, A. Oettlé, F. Schrenk, and R. Machiarelli. We also acknowledge the Nespos Society (www. nespos.org). For the elaboration of micro-CT record of the comparative material we thank F. Bernardini, B. Duployer, Y. M. Hou, and C. Tenailleau. We are also indebted to W. Liu and S. Xing for sharing micro-CT record of some Chinese fossil specimens, and also for scientific discussion and comments. Research supported by Strategic Priority Research Program of Chinese Academy of Sciences (No. XDB26000000), State Key Laboratory of Palaeobiology and Stratigraphy (Nanjing Institute of Geology and Palaeontology, CAS; No. 173119), National Natural Science Foundation of China (No. 41702026, No. 41872030).