Pascal Delange
email: pascal.delange@polytechnique.edu

Steffen Backes
email: steffen.backes@polytechnique.edu

Ambroise Van Roekeghem
email: ambroise.vanroekeghem@gmail.com

Leonid Pourovskii
email: leonid@cpht.polytechnique.fr

Hong Jiang
email: jianghchem@pku.edu.cn

Silke Biermann
email: biermann@cpht.polytechnique.fr

Novel Approaches to Spectral Properties of Correlated Electron Materials: From Generalized Kohn-Sham Theory to Screened Exchange Dynamical Mean Field Theory

The most intriguing properties of emergent materials are typically consequences of highly correlated quantum states of their electronic degrees of freedom. Describing those materials from first principles remains a challenge for modern condensed matter theory. Here, we review, apply and discuss novel approaches to spectral properties of correlated electron materials, assessing current day predictive capabilities of electronic structure calculations. In particular, we focus on the recent Screened Exchange Dynamical Mean-Field Theory scheme and its relation to generalized Kohn-Sham Theory. These concepts are illustrated on the transition metal pnictide BaCo 2 As 2 and elemental zinc and cadmium.

Introduction

Technological progress has been intimately related with progress in materials science since its very early days. The last century has seen the development of refined capabilities for materials elaboration, characterisation and control of properties, culminating among others in the unprecedented possibilities of the digital age: materials properties have become an object of theoretical simulations, stimulating systematic searches for systems with desired characteristics. This branch of condensed matter physics is continuously evolving into a true new pillar of materials science, where the theoretical assessment of solid state systems becomes a fundamental tool for materials screening.

Obviously, the success of this program hinges on the degree of predictive power of modern simulation techniques, which have to allow for performing calculations without introducing adjustable parameters. This requirement is even more severe since the most popular branch of ab initio calculations, the Density Functional Theory (DFT) approach, is restricted to ground state properties of the solid, while most properties of potential technological interest stem from excited states: the calculation of any type of transport phenomenon, for example electric, thermal, thermoelectric or magnetoelectric transport, of optical or spectroscopic properties, or of magnetic, charge-or orbital susceptibilites requires the accurate assessment of highly non-trivial response functions within a finite-temperature description.

It is not a coincidence that the materials with the most exotic features are also most challenging to get to grips with from a theoretical point of view: both the complexity of their properties and the difficulty of their description are in fact a consequence of the very nature of the underlying electronic structure: typically, one is dealing with transition metal, lanthanide or actinide compounds, where the electrons in partially filled d-or f -shells strongly interact with each other through electron-electron Coulomb interactions, leading to highly entangled many-body quantum states.

Even when using a simplified description of the solid in terms of an effective lattice model, assessing these quantum correlated states is a tremendous challenge. If local quantum fluctuations on a given atomic site are dominant, dynamical mean-field theory (DMFT) 1,2) yields an accurate description of the system. In particular, DMFT captures both the strong coupling Mott-insulating limit and the weakly interacting band limit of the Hubbard model, and is able to describe the salient spectral features of a correlated metal with coexisting quasi-particle and Hubbard peaks even in the intermediate correlation regime. In the multi-orbital case, additional degrees of freedom can lead to even richer physics with unconventional (e.g. orbital-selective) behavior, 3) or complex ordering phenomena. [4][5][6][7][8][9] In order to recover the material-specific character of the calculations, DMFT has been combined with DFT 10,11) into the so-called "DFT+DMFT" scheme, which is nowadays one of the most popular workhorses of electronic structure theory for correlated electron materials. Successful applications include transition metals, [12][13][14] transition metal oxides, 4,[15][16][17][18][19][20][21][22][23][24][25][26] lanthanide 27,28) or actinide [29][30][31][32] systems. Comparisons be-tween calculated and measured spectral functions have sometimes led to impressive agreement. 33,34) However, despite these successes the construction of the Hamiltonian used in these calculations is rather ad hoc and remains a limitation to the predictive power of the approach. Therefore, the elaboration of more systematic interfaces between electronic structure and many-body theory has become an active area of research. 35) In this work, we review a recent scheme combining screened exchange and DMFT. [36][37][38] We discuss its relation to many-body perturbation theory and generalized Kohn-Sham (KS) Theory and analyse the effects included at the different levels of the theory. As an illustration we describe its application to the low-energy spectral properties of the cobalt pnictide BaCo 2 As 2 .

How does a theory thas has been designed for strongly correlated materials reduce to an a priori simpler version in the case of a weakly correlated system ? For methods based on DFT that include only the static on-site Coulomb interaction between localized states, like DFT+U 39) or DFT+DMFT, the spectrum in the limit of vanishing Hubbard interactions reduces trivially to the Kohn-Sham spectrum of DFT. For more complex interfaces of electronic structure and many-body theory which we will discuss here, however, it is a non-trivial question and allows for interesting possibilities to check their consistency. Here, we will discuss this issue on the example of the transition metals zinc and cadmium. This also allows us to comment on the challenge of including states in a wider energy range, identifying a challenging problem on these seemingly "simple" systems.

The paper is organised as follows: in section 2 we review how Screened Exchange Dynamical Mean Field Theory derives from the combined many-body perturbation theory + dynamical mean field scheme "GW+DMFT". In section 3 we analyse the relation of screened exchange schemes to generalized KS theory. Section 4 provides an example of the application of such methods to BaCo 2 As 2 , while section 5 discusses the electronic structure of elemental zinc and cadmium within a simplified approach. Finally, we conclude in section 6.

From GW+DMFT to Screened Exchange Dynamical Mean Field Theory

Improving the predictive power of methods that combine electronic structure and many-body theory poses the challenge of properly connecting the two worlds, without double counting of interactions or screening. At the heart of this challenge lies the mismatch between the density-based description of DFT and the Green's function formalism used at the many-body level, as well as the difficulty of incorporating the feedback of high-energy screening processes governed by the unscreened Coulomb interaction onto the low-energy electronic structure. Conceptually speaking, these difficulties can be avoided by working on a large energy scale in the continuum with the full long-range Coulomb interactions, and a Green's function-based formalism even at the level of the weakly correlated states. These features are realised within the combined many-body perturbation theory and dynamical mean field theory scheme "GW+DMFT": [40][41][42][43][44][45][46][47][48][49] screening is assessed by the random phase approximation in the continuum, augmented by a local vertex correction, while the starting electronic structure for the DMFT calculation can be roughly interpreted as a "non-local GW" calculation. 50) Screened Exchange Dynamical Mean Field Theory [36][37][38] can be understood as an approximation to this full GW+DMFT scheme. It is based on the recent observation 36,50) that within the GW approximation the correction to LDA can be split into two contributions: a local dynamical self-energy Σ loc (ω) and a k-dependent but static self-energy Σ nloc (k), which does not contain any local component. If such a separation was strictly valid in the full energy range, the non-local static part of the full self-energy would be given by the non-local Hartree-Fock contribution, since the dynamical part vanishes at high frequency. In many realistic systems such a decomposition holds to a good approximation in the low-energy regime that we are interested in, where the static part Σ nloc (k) is quite different from the Fock exchange term. It is approximately given by a screened exchange selfenergy, leading to a decomposition of the GW self-energy into

Σ GW = [GW(ν = 0)] nonloc + [GW] loc .
Here, the first term is a screened exchange contribution arising from the screened interaction W(ν) evaluated at zero frequency W(ν = 0). The second term is the local projection of the GW self-energy. It is simply given by the GW self-energy evaluated using a local propagator G loc and the local screened Coulomb interaction W loc . Exactly as in GW+DMFT, Screened Exchange Dynamical Mean Field Theory replaces this term by a nonperturbative one: it is calculated from an effective local impurity problem with dynamical interactions. In current practical applications of Screened Exchange + DMFT the RPAscreened Coulomb potential W has been replaced by its longwavelength limit, which reduces to a simple Yukawa-form. 36) Quite generally the dynamical character of the interactions results in an additional renormalization Z B of the hopping amplitudes, which can be understood as an electronic polaron effect: the coupling of the electrons to plasmonic screening degrees of freedom leads to an effective mass enhancement corresponding to the hopping reduction, manifesting itself as a narrowing of the band. This effect can be estimated from the plasmon density of modes as given by the imaginary part of the frequency dependent interaction W. An explicit expression for Z B has been derived in Ref. 51.

Relation to Generalized Kohn-Sham Theory

Screened Exchange Dynamical Mean Field Theory can also be viewed as a specific approximation to a spectral density functional theory based on the Generalized Kohn-Sham (GKS) scheme of Seidl et al. 52) In GKS theory, alternative choices for the reference system that are different than the familiar Kohn-Sham system of DFT are explored. In particular, a generalized Kohn-Sham scheme where the reference system is a screened exchange Hamiltonian can be constructed. The main motivation for the inclusion of screened exchange in the literature has been to improve upon the band gap problem in semiconductors. Indeed, it can be shown that the screened exchange contribution, which corresponds to a non-local potential, effectively reintroduces to some degree the derivative discontinuity that is missing in the pure DFT description based on local exchange-correlation potentials. 53) Since the derivative discontinuity corresponds to the discrepancy between the true gap and the Kohn-Sham gap in exact DFT, a substantial improvement of the theoretical estimate for band gaps can be expected on physical grounds and has indeed been found. Here, our goal is somewhat different: motivated by the analysis of the role of screened exchange in GW+DMFT described above, we would like to connect the Screened Exchange DMFT scheme introduced above to generalized KS schemes making direct use of the non-local screened exchange potential.

With this aim in mind, we briefly review the generalized Kohn-Sham construction in the case an effective Kohn-Sham system including screened exchange. Hereby, we follow closely Seidl et al., 52) both in notation and presentation. First, one defines a functional

S [Φ] = Φ T Φ + U H {φ i } + E sx x {φ i } (1)
that includes, in addition to the familiar kinetic energy term Φ| T |Φ and the Hartree energy U H {φ i } also the screened Fock term

E sx x {φ i } = - N i< j drdr ′ × φ * i (r)φ * j (r ′ )e -k TF |r-r ′ | φ j (r)φ i (r ′ ) |r -r ′ | (2)
Here, Φ are Slater determinants of single-particle states φ i . k TF is the Thomas Fermi wave vector. In order to derive a functional of the density Seidl et al. define a functional F s via the minimisation

F S ρ = min Φ→ρ(r) S [Φ] = min {φ i }→ρ(r) S {φ i } (3)
Next we define the energy functional

E S {φ i }; v e f f = S {φ i } + drv e f f (r)ρ(r) (4)
where now the potential v e f f does not only include the external potential v as in usual DFT, but also a contribution by the exchange-correlation part

v e f f = v + v sx xc ρ . ( 5
)
The additional contribution, the generalized (local) exchangecorrelation potential

v sx xc = ∂E sx xc ρ ∂ρ . ( 6
)
is the functional derivative of the generalized (local) exchange-correlation functional

E sx xc ρ = E xc ρ -E sx x ρ + T ρ -T sx ρ (7)
which comprises the difference between the exchangecorrelation potential of standard Kohn-Sham DFT and the non-local exchange energy defined above, as well as the difference between the kinetic energies of the standard and generalized Kohn-Sham systems. The functional derivative will eventually have to be evaluated self-consistently at the converged density. This construction leads to the generalized Kohn-Sham equations

-∇ 2 φ i (r) + v(r)φ i (r) + u ([ρ]; r) φ i (r) -dr ′ v sx,NL x (r, r ′ )φ i (r ′ ) + v sx xc ([ρ]; r) φ i (r) = ε i φ i (8)
with the Hartree potential u and the non-local screened Fock potential

v sx,NL x (r, r ′ ) = - N j=1 φ i (r)e -k TF |r-r ′ | φ * j (r ′ ) |r -r ′ | . ( 9
)
and the effective (local) generalized Kohn Sham potential v sx xc defined above. The generalized Kohn-Sham equations have the form

Ô {φ i } φ j + ve f f φ j = ǫ j φ j ( 10
)
where Ô is a non-local operator, generalizing the standard Kohn-Sham operator consisting solely of kinetic energy and Hartree potential. The ground state energy for a system in the external potential v is then given by the expression

E SEx-DFT [v] = F S ρ S 0 v e f f + E sx xc ρ S 0 v e f f + drv(r)ρ S 0 v e f f (11)
Now, the relation to Screened Exchange DMFT is becoming clear: one may construct a spectral density functional in the same spirit as in DFT+DMFT, 54) but starting from the generalized KS functional. In 55, the expression for the total energy within the standard DFT+DMFT case was derived to be

E = E DFT - l ǫ KS l + H KS + (H int -H dc ) ( 12
)
where l ǫ KS l is the sum of the occupied Kohn-Sham eigenvalues, H KS = tr [H KS Ĝ], and H int and H dc denote the local interaction part of the Hamiltonian and the corresponding double counting term, respectively.

Instead of using the usual Kohn-Sham Hamiltonian for the construction of the one-body part, Screened Exchange DMFT relies on the generalized Kohn-Sham reference system that includes the screened exchange potential. The generalization of (12) to the present case thus replaces the Kohn-Sham Hamiltonian H KS in the expression for the energy by its non-local form, keeping track of the effective potential part:

E = E SEx-DFT - l ǫ SEx-KS l + Ô + ve f f -6 -5 -4 -3 -2 -1 0 1 2 3 Γ X M Γ E -E F (eV) -6 -5 -4 -3 -2 -1 0 1 2 3 -6 -5 -4 -3 -2 -1 0 1 2 3 -6 -5 -4 -3 -2 -
+ (H int (V ee , λ s , ω s ) -H dc ) (13)
Furthermore, the local interaction term is taken in the more general form of a dynamical interaction, thus corresponding to a local Hubbard term with unscreened interactions and local Einstein plasmons of energy ω s coupling to the electrons via coupling strength λ s . This concludes our description of the generalized Kohn-Sham interpretation of Screened Exchange DMFT, resulting in particular in an energy functional expression. However, in the practical calculations presented in the following we do not minimise the full energy expression as given above, but rather work at the converged DFT density and then investigate spectral properties using the Screened Exchange DMFT formalism. This amounts to a one-shot Screened Exchange-DFT+DMFT calculation that uses the DFT density as a starting point. The advantage of such an approach is obvious: Numerically, this procedure allows us to avoid the expensive evaluation of non-local exchange terms within the selfconsistency cycle of GKS theory. Moreover, as is well-known, while severe deviations of the true spectrum from the Kohn-Sham spectrum are quite common, the ground state density obtained even from approximate DFT functionals is often a good representation of the true one. In the case of the exact DFT functional, our approach would also lead to the exact ground state density and energy, with additional improvements of the spectrum over standard Kohn-Sham DFT. 3d states changes from d 6 to d 7 , which strongly reduces the degree of Coulomb correlations, making BaCo 2 As 2 a moderately correlated compound. 36,56) Indeed, the power-law deviations from Fermi liquid behavior above an extremely low coherence temperature discussed in the Fe-based compounds 57) are a consequence of the d 6 configuration and strong intraatomic exchange interactions. This is no longer the case in the cobalt pnictides, where angle-resolved photoemission spectroscopy (ARPES) identifies well-defined and long-lived quasiparticle excitations with relatively weak mass renormalization. Nevertheless, the DFT-LDA derived Fermi surface differs from experiment. 36,56,58,59) Therefore, BaCo 2 As 2 provides an ideal testing ground for new approaches to spectroscopic properties. We do not discuss here the details of yet another interesting question, which is the absence of ferromagnetism despite a high value of the DFT density of states at the Fermi level, suggestive of Stoner ferromagnetism, but refer the reader to Ref. 36, where the solution to this puzzle was discussed in detail. Fig. 1 shows the DFT band structure of BaCo 2 As 2 , in comparison to a screened exchange calculation. As in the ironbased pnictides, the dominantly 3d-derived states are located around the Fermi level; in this case in a window of about -3 eV to 2 eV. As compared to the parent iron pnictides with 3d 6 configuration of the Fe shell, the Fermi surface topology is modified due to the larger 3d 7 filling. The hole pocket at the Γ point that is present in most Fe-based pnictide compounds is pushed below the Fermi level, as well as the band forming the electron pocket at M, which is now fully filled. In standard DFT-LDA (see Fig. 2), a characteristic flat band of dominant x 2y 2 character lies directly on the Fermi level around the M point, giving rise to a huge peak in the density of states. The close proximity of this band to the Fermi level will render its energy highly sensitive to the details of the calculation and thus provides a perfect benchmark for improved computational techniques.

Comparison to angle-resolved photoemission experiments, 36) as reproduced in Fig. 3, reveals that the overall bandwidth of the LDA band structure is too wide by roughly a factor of 1.5. In a DFT+DMFT calculation, as shown in Fig. 4, this band renormalisation is reproduced, giving an overall occupied bandwidth of about 1.5 eV. The fine details of the Fermi surface, and in particular the position of the flat x 2y 2 band are however not well described. For a detailed comparison we refer the reader to Ref. 36.

Including screened exchange in the form of a Yukawa potential on top of DFT, as shown in Fig. 1 widens the band by a considerable amount as compared to DFT-LDA. Even more striking are the modifications at the Fermi level: the x 2y 2 band has been shifted above the Fermi level, with an energy at the M point of about 0.15 eV above E F . While the lowering of the filling of this band improves the agreement between theory and experiment, the shift is too large to reproduce the experimental Fermi surface. However, when applying DMFT with dynamical interactions on top of the screened exchange Hamiltonian, equivalent to the one-shot Screened Exchange DMFT procedure described above, this band is renormalized by the electronic interactions and ends up again close (but above!) the Fermi level. Its energy is slightly higher than in DFT-LDA and is now in excellent agreement with the experimental Fermi surface. A more detailed analysis presented in Ref. 36) reveals that also the higher energy features such as the bands within the range of up to 2 eV are well reproduced in this one-shot Screened Exchange DMFT approach.

The comparison of the band structures and spectral func- tions within the different computational schemes illustrates the effects of the different terms in an instructive way: Improving on the description of exchange by replacing the local exchange as contained in DFT-LDA by a non-local screened Fock exchange widens the bandwidth and strongly "overcorrects" the Fermi surface. Improving at the same time on the correlation part by applying DMFT with frequencydependent interactions then leads to a partial cancellation of the band widening, giving an overall bandwidth surprisingly close to the LDA+DMFT one. This suggests that the good description of the overall bandwidth in DFT+DMFT is the result of an error cancellation between the local approximation to exchange and partial neglect of correlations. However, the Fermi surface is strongly modified by the non-local corrections, leading to a significant improvement over LDA+DMFT and resulting in good agreement with experiment. These findings might suggest that possible inconsistencies between theoretical and experimental Fermi surfaces that are observed in many pnictide materials could be corrected by adding simple non-local self-energy corrections stemming from screened exchange effects. Whether this treatment can fully account for the "red-blue shift" of Fermi surface pockets found in the experimental-theoretical comparisons 60) is a most interesting topic for future research. The current example may give rise to optimism.

How does Screened Exchange Dynamical Mean Field

Theory behave for weakly correlated materials?

Hamiltonians built as combinations of a DFT part and local Hubbard-type interaction terms trivially reduce to the DFT Kohn-Sham electronic structure when assuming that in weakly correlated materials the effective local interactions become negligible. The question of the recovery of the weakly interacting limit is, however, more interesting in the case of Screened Exchange DMFT. While the static part of the effective local interaction may be assumed to lose its importance, band widening by the replacement of the DFT exchange correlation potential by the non-local exchange-correlation GW self-energy persists. On the other hand, plasmonic effects are also present in weakly correlated materials and continue to renormalize the low-energy band structure through electronplasmon coupling. This raises the question of what the resulting spectra for weakly correlated materials look like in screened exchange + DMFT. In Ref. 38, this question has been studied for early transition metal perovskites, where it was found that the band widening effect induced by non-local exchange and the electronic polaron effect counteract each other and tend to approximately cancel, thus resulting in a low-energy electronic structure close to the DFT Kohn-Sham band structure as long as static Hubbard interactions are disregarded. Here, we address this question in the case of the seemingly "simple" transition metals zinc and cadmium.

Both elements nominally display a d 10 configuration, with fully occupied 3d orbitals in the case of Zn and 4d in the case of Cd, the dominantly d-derived bands being located several eV below the Fermi level. In Fig. 7 and Fig. 8 we show the band structure calculated within DFT for both materials. Here and in the following we use the experimental crystal structure. DFT puts the occupied d states at around -8 eV in Zn and -9 eV in Cd. The conducting states of this transition metal are formed by dispersive 4s(5s) in Zn(Cd) states around the Fermi level, that hybridize with the p-manifold.

These facts raise the immediate expectation of negligibly small correlation effects on the occupied d shells. An effective Hubbard interaction calculated for the d-manifold within the constrained random phase approximation coincides with the fully screened interaction, since as a consequence of the complete filling of the d-shell there are no intra-d transitions to be cut out as opposed to an open shell system, where transitions inside the shell contribute to screening effects. Figs. 9 and 10 display the local component of this fully screened interaction projected on the d-manifold. The lowfrequency limit approaches a value of 4.8 eV and 3.1 eV for Zn and Cd respectively. Even though this value is similar to their oxides, which are open shell systems, where correlations on the d states are significant, the high binding energy of these states far away from the Fermi level effectively prevents dynamical fluctuations. This suggests that screened exchange + DMFT should in fact reduce to screened exchange renormalized by the bosonic factor Z B discussed above. Nevertheless, this does not mean that static effects of the interaction are properly treated in DFT. Even if this were the case, there is no reason that the DFT Kohn-Sham spectrum, being derived from an effective non-interacting system, provides an accurate description of the experimental situation.

Figs. 11 and 12 compare the experimental photoemission spectra 61,62) from the literature to the density of states (DOS) derived from DFT and Hartree-Fock (HF) theory. The resulting discrepancy in terms of an underestimation of the binding energy of the d states in DFT of several eV had been noted in the literature before: [63][64][65] Norman et a.l 63) discussed it in terms of a self-interaction error, proposing a correction in terms of an approximate substraction of self-interaction contributions contained in DFT. 66) Hartree-Fock calculations are, on the one hand, self-interaction free, but on the other handdue to the absence of screening -widen all bands and place the d-bands far too low in energy, as can be seen in Figs. 11 and12.

Since local dynamical correlations can assumed to be small as discussed before, an improved treatment of the screened interaction and the self-interaction correction at the same time is likely to improve the shortcomings of Hartree-Fock, resp. DFT. Here we will discuss the two possible extensions of Screened Exchange plus a bosonic renormalization factor Z B and the GW approximation. The computationally cheaper option of screened exchange is including only static exchange contributions with a Yukawa-type interaction potential, and an effective renormalization Z B which originates from the spectral weight transfer to plasmonic excitations. GW is computationally more demanding, but has the advantage of treating the dynamical part of the screening and correlation. Even though both methods include a self-interaction correction, the selfinteraction contained in the Hartree term is not completely cancelled since the exchange contributions are derived from a screened interaction and not the bare one, as opposed to the Hartree term.

In Figs. 13 and 14 we show comparisons of the DOS of zinc and cadmium, calculated within DFT, screened Exchange(+Z B ) and GW 67) to photoemission spectra. In the GW calculation we used 7 × 7 × 3 k-points and 5 additional high-energy local orbitals. Interestingly, while in both systems the GW approximation provides a significant correction of the DFT Kohn-Sham spectrum in the right direction, but still underestimates the binding energy, the screened exchange scheme places the d-states too low in energy for Zn while providing a slightly better estimate than GW in Cd. The addition of the bosonic renormalization factor Z B merely renormalizes the d bandwidth but keeps the average level position constant.

This raises the interesting question of which effects are missing in screened exchange and GW? The incomplete cancellation of the self-interaction in both approaches is expected to lead to an overall underestimation of the binding energy, since the additional unphysical interaction increases the energy of the d states. A more accurate estimate of this term would lead to an improvement of GW in both systems, but an even larger error of screened exchange in Zn. Another effect neglected in screened exchange is the Coulomb hole contribution: this term, discussed by Hedin as part of the "Coulomb hole screened exchange (COHSEX)" approximation translates the fact that the presence of an electron at a position r pushes away charge at r (in the language of a lattice model, the charge-charge correlation function exhibits a reduction of the double occupancy), and the interaction of this effective positive charge with the electron presents an energy gain expressed in the form of an interaction of the electron with a "Coulomb hole". This term, contained in GW but not in screened exchange, also increases the binding energy.

This leads to the overall picture that screened exchange with the inclusion of the static corrections just discussed has a tendency to overestimate the binding energy in general, while GW underestimates it. The obvious difference between the two methods is the dynamical treatment of the screened interaction, which is treated more appropriately in GW, but it is not clear a priori whether the static approximation of the screened interaction or the approximated form of the screened interaction in terms of a Yukawa potential gives rise to the difference between screened exchange and GW. The GW description of zinc and cadmium is close to the experiment. The remaining discrepancy to experiment is likely explained by remaining self-interaction contributions and/or missing self-consistency. Self-consistency (or quasiparticle self-consistency) has been investigated in the homogeneous electron gas 68) and various solid state systems. 69) These questions are left for future work.

Conclusions

In this work we have reviewed and applied existing as well as novel approaches to obtain spectral properties of correlated electron materials. Guided by the need of a proper treatment of the long-range Coulomb interaction and non-local exchange effects we presented a lightweight version of the general GW+DMFT approach, the so-called Screened Exchange DMFT. It can be derived as a simplication to GW+DMFT in terms of a generalized screened exchange DFT scheme where local interactions are treated by dynamical DMFT.

Analysis and the application of a simplifed form of this scheme to BaCo 2 As 2 indeed showed that non-local exchange and electronic screening lead to significant corrections to the electronic spectrum which are necessary to obtain a proper description of the experimental observations. Furthermore, we discussed the case of the elemental transition metals Zn and Cd, where strong local correlations are unimportant but the position of the occupied d manifold is very sensitive to a proper treatment of screened exchange effects.

SB and AvR thank their collaborators of Refs.36 and 37 for the fruitful collaborations and discussions. This work was supported by IDRIS/GENCI Orsay under project t2017091393, the European Research Council under its Consolidator Grant scheme (project 617196), and the ECOS-Sud grant A13E04. HJ acknowledges financial support from National Natural Science Foundation (21373017, 21621061) and the National Basic Research Program of China (2013CB933400). E -E F (eV)

Fig. 1 .

 1 Fig. 1. (Color online) Kohn-Sham band structure of BaCo 2 As 2 within DFT-LDA (red lines) and the screened exchange approximation (black lines).

Fig. 2 .

 2 Fig. 2. (Color online) Kohn-Sham band structure of BaCo 2 As 2 within DFT-LDA, projected on the d x 2 -y 2 orbital.

Fig. 3 .

 3 Fig. 3. (Color online) Angle-resolved photoemission spectrum of BaCo 2 As 2 . Adapted from Ref. 36.

Fig. 4 .

 4 Fig. 4. (Color online) k-resolved spectral function of BaCo 2 As 2 within LDA+DMFT. Adapted from Ref. 36.

Fig. 5 .

 5 Fig. 5. (Color online) k-resolved spectral function of BaCo 2 As 2 within Screened Exchange DMFT (with dynamical interactions). Adapted from Ref. 36.

Fig. 6 .

 6 Fig. 6. (Color online) k-resolved spectral function of BaCo 2 As 2 within Screened Exchange DMFT (with dynamical interactions), projected on the d x 2 -y 2 orbital. Adapted from Ref. 36.

Fig. 7 .

 7 Fig. 7. (Color online) The band structure of elemental Zn calculated within DFT. The orbital character is indicated by the intensity of the different colors.

Fig. 8 .

 8 Fig. 8. (Color online) The band structure of elemental Cd calculated within DFT. The orbital character is indicated by the intensity of the different colors.

Fig. 9 .

 9 Fig. 9. (Color online) The fully screened effective local Hubbard interaction on the 3d manifold for Zn.

Fig. 10 .

 10 Fig. 10. (Color online) The fully screened effective local Hubbard interaction on the 4d manifold for Cd.

Fig. 11 .

 11 Fig. 11. (Color online) The Density of States of elemental Zn calculated within Density Functional Theory (DFT, black solid line) and Hartree-Fock (HF, red solid line) in comparison with Photoemission experiments61,62) (dashed line, symbols). Density Functional Theory (DFT) calculations underestimate the binding energy of the Zn 3d states, while the HF overestimates the binding energy significantly (see explanation in the text).

Fig. 12 .

 12 Fig. 12. (Color online) The Density of States of elemental Cd calculated within Density Functional Theory (DFT, black solid line) and Hartree-Fock (HF, red solid line) in comparison with Photoemission experiments 61)(dashed line, symbols). Density Functional Theory (DFT) calculations underestimate the binding energy of the Cd 3d states, while the HF overestimates the binding energy significantly (see explanation in the text).

Fig. 13 .Fig. 14 .

 1314 Fig. 13. (Color online) The Density of States of elemental Zn calculated within different theoretical methods (solid lines), in comparison with Photoemission experiments61,62) (dashed line, symbols). Density Functional Theory (DFT) calculations significantly underestimate the binding energy of the Zn 3d states, while the GW approximation obtains a much better agreement. Screened exchange overestimates the binding energy significantly.

As a first illustration of Screened Exchange DMFT, we review calculations on BaCo 2 As 2 , which is the fully Cosubstituted representative of the so-called "122" family of the iron-based superconductors, isostructural to the prototypical parent compound BaFe 2 As 2 . The Fe → Co substitution has however important consequences: the nominal filling of theΓ X M Γ


Strain tensor evaluation in polycrystalline materials by scanning high energy X-ray diffraction

András Borbély

Mines Saint-Etienne, Univ. Lyon, CNRS, UMR 5307 LGF, Centre SMS, F -42023 Saint-Etienne France 3 rd generation high-energy synchrotron sources made possible during the last 20 years the development of several diffraction-imaging methods, which led to meaningful physical insights into the real structure and dynamics of bulk polyrcystalline materials. A not exhaustive list could include the work on stress corrosion cracking (King et al., 2008), predictions of crystal plasticity (Pokharel et al., 2014), stress variation in copper through-Si vias (Levien et al., 2015), martensitic transformation [START_REF] Sedmák | [END_REF] or on long-range symmetry breaking in embedded ferroelectrics (Simons et al., 2018). The pioneering work of [START_REF] Poulsen | Three-Dimensional X-ray Diffraction Microscopy: Mapping Polycrystals and Their Dynamics[END_REF] paved the way for the majority of the techniques using monochromatic radiation, commonly called Three Dimensional X-Ray Diffraction (3DXRD) or High Energy Diffraction Microscopy (HEDM) (Suter et al., 2006). Since then, a wide variety of techniques have been developed, which are mostly classified based on the distance between the detector and the specimen in a) far-field or b) near-field techniques. However, from the point of view of the final reconstruction describing the spatial distribution of the crystallographic orientation (implicitly the grain shape) and strain/stress, a classification based on their spatial and angular resolutions seems to be more appropriate. According to this two-resolution criterion all existing techniques can be included in one of the following three categories: a) High-spatial and low-angular resolution methods use a near-field setup with a high resolution detector (effective pixel size of about 1 µm) closely placed downstream the specimen at a distance of a few millimeters. Their angular resolution is relatively low (~0.1°), but they can resolve the structure at micrometer/sub-micrometer length scale over millimeter sample sizes. Grain and intra-grain orientation is obtained via scans done either with a broad (Ludwig, 2008) or a planar beam (Li et al., 2013), the beam size being usually adapted to the grain size to avoid peak overlap. b) Low-spatial and high-angular resolution methods involve a low resolution detector with pixel size of about 50 to 200 µm placed in the far-field (~0.2 to 1 m as a function of the pixel size and beam energy). Since diffraction peak positions can be determined with sub-pixel accuracy (Borbély et al., 2014) this setup has high angular resolution (< 0.01°) that allows the elastic strain to be determined, however, the grain shape remains unknown. Measurements are usually performed with a broad beam. c) Finally, high-spatial and high-angular resolution methods enable characterization of both the crystallographic orientation and strain inside single grains at sub-micrometer length scale, even for deformed materials. It is expected that these "holy grail" methods will provide the missing local experimental evidence for understanding unsolved problems in materials science such as polycrystal plasticity, recrystallization or damage initiation. There are already two methods achieving high-spatial and high-angular resolution, the Differential-Aperture X-ray Microscopy (DAXM, Larson et al. 2002) and Dark Field X-ray Microscopy (DFXM, Simons et al. 2015). The first uses a polychromatic pencil beam with sub-micron cross section for identifying the local orientation from the resulting Laue pattern and an additional energy scan (Chung & Ice, 1999) for strain determination. DFXM uses monochromatic radiation and compound refractive lenses in the diffracted beam, which magnify a small diffractive volume of the analyzed grain/subgrain.

To improve the spatial resolution of the low-spatial & high-angular resolution methods Hayashi et al., (2015) proposed a "scanning 3DXRD" approach, using a monochromatic pencil beam, where the spatial resolution could be controlled by the beam size. Combining a lateral scanning of the sample and tomographic approach (rotation around an axis perpendicular to the beam) the authors were able to reconstruct the variation of the crystallographic orientation inside single grains. Very recently, the method was also applied for determining local strain and stress tensors (Hektor et al., 2019, Hayshi et al., 2019). These developments, however, considered a simplified reconstruction procedure based on the assumption that the crystallographic structure in a given voxel can be refined from the subset of diffraction peaks generated by the beams that intersect its volume. Evidently, this assumption, which could be referred to as "single crystal reconstruction" (SCR), does not obey the tomographic principle, when the recorded signal is an average response over the X-ray beam path. The paper of Henningsson et al., (2020) in this journal issue suggests two solutions to the tomographic problem of lattice strain, both considering all voxels in the beam path. The first called "polycrystal reconstruction" (PCR) is based on the minimization of the Euclidean distance between the experimental and simulated peak centers measured on the detector. The simulation approximates the scattered intensity by assuming that the contribution of a given voxel to the measured peak is proportional to the fraction of the voxel illuminated by the beam. The second solution is inspired by the common algebraic reconstruction (ART) method, as used, for example, by Alpers et al. (2006) to reconstruct grain shapes from 3DXRD data, extended to reconstruct tensor fields of smoothly varying intragranular strain. Both the PCR and ART methods were checked against simulated data and performed better than the simple SCR method. Finally, the strain state in a columnar tin grain was determined from experimental data measured at the ID11 beam line of the ESRF using a pencil beam cross with section of 0.25 µm x 0.25 µm. The reconstructions obtained with all three methods (SCR, PCR & ART) evidence the existence of a strain gradient in the grain.

By establishing the correct tomographic approach, the paper of Henningsson et al., (2020) upgraded "scanning 3DXRD" to a viable alternative to DAXM and DFXM. Since it only involves lateral sample scanning and rotation around a single axis the technique can be easily implemented in practice, which helps its wide use. However, to gain more confidence in the results, checking of real reconstructions against other techniques (such as near-field HEDM, similarly to Renversade et al., (2016) would be valuable. Today, scan duration for representative volumes is somewhat long (12h for scanning a volume of 37³ µm³ with a beam size of 1 µm x 1 µm in steps of 1 µm, Hayasi et al., 2019), however, with the advent of 4 th generation synchrotron sources and fast detectors the scanning time can be reduced by one order of magnitude, which already allows in situ studies to be performed during common beam time duration.


Yann Heuzé

Simeon A Boyadjiev

Jeffrey L Marsh

Alex A Kane

Elijah Cherkez

James E Boggan

Joan T Richtsmeier

New insights into the relationship between suture closure and craniofacial dysmorphology in sagittal nonsyndromic craniosynostosis

Keywords: Isolated craniosynostosis, geometric morphometrics, suture closure pattern, suture closure sequence, dysplastic and compensatory growth, scaphocephaly, skull shape

Premature closure of the sagittal suture occurs as an isolated (nonsyndromic) birth defect, or as a syndromic anomaly in combination with other congenital dysmorphologies. The genetic causes of sagittal nonsyndromic craniosynostosis (NSC) remain unknown. Though variation of the dysmorphic (scaphocephaly) skull shape of sagittal NSC cases has been acknowledged, this variation has not been quantitatively studied in 3D. We have analyzed the computed tomography skull images of 43 infants (aged from 0.9 to 9 months) with sagittal NSC using anatomical landmarks and semilandmarks to quantify and characterize the within-sample phenotypic variation. Suture closure patterns were defined by dividing the sagittal suture into three sections (anterior, central, posterior) and coding each section as "closed" or "fused". Principal components analysis of the Procrustes shape coordinates representing skull shape of 43 cases of NSC did not separate individuals by sex, chronological age, or dental stages of deciduous maxillary first molar. However, analysis of suture closure pattern allowed separation of these data. The central section of the sagittal suture appears to be the first to fuse. Then, at least two different developmental paths towards complete fusion of the sagittal suture exist; either the anterior section or the posterior section is the second to fuse. Results indicate that according to the sequence of sagittal suture closure patterns, different craniofacial complex shapes are observed. The relationship between craniofacial shape and suture closure indicates not only which suture fused prematurely (in our case the sagittal suture), but also the pattern in which the suture closes. Whether these patterns indicate differences in aetiology cannot be determined with our data and requires analysis of longitudinal data, most appropriately of animal models where prenatal conditions can be monitored.

Introduction

Craniosynostosis is defined as the premature fusion of one or more cranial suture and is a common malformation occurring in 3-5 per 10,000 live births in all ethnic groups [START_REF] Cohen | Epidemiology of craniosynostosis[END_REF]. Craniosynostosis of the sagittal suture is the most common form of craniosynostosis accounting for 55-60% of all craniosynostoses [START_REF] Anderson | Craniosynostosis: a survey of 204 cases[END_REF][START_REF] Shilito | Craniosynostosis. A review of 519 surgical patients[END_REF][START_REF] Hunter | Craniosynostosis. I. Sagittal synostosis: its genetics and associated clinical findings in 214 patients who lacked involvement of the coronal suture(s)[END_REF][START_REF] Lajeunie | Genetic study of scaphocephaly[END_REF].

Nonsyndromic (isolated) craniosynostosis (NSC) of the sagittal suture represents the majority of all cases of sagittal synostosis with a birth prevalence of 1.4-1.5 per 10,000 live births [START_REF] Lajeunie | Genetic study of scaphocephaly[END_REF][START_REF] Boulet | A population-based study of craniosynostosis in Metropolitan Atlanta[END_REF]. Though environmental and genetic risk factors for craniosynostosis have been identified, no single risk factor has emerged as being necessary or sufficient to cause premature fusion of any of the calvarial sutures. In utero head constraint has historically been considered as an environmental component of craniosynostosis [START_REF] Koskinen-Moffett | Cranial Synostosis and intra-uterine compression: a developmental study of human sutures[END_REF][START_REF] Graham | Craniofacial deformation[END_REF][START_REF] Cohen | Epidemiology of craniosynostosis[END_REF]. Twin studies reported by [START_REF] Lajeunie | Genetic considerations in nonsyndromic midline craniosynostoses: a study of twins and their families[END_REF] suggested a multifactorial pattern of inheritance, supporting both an environmental (in utero head constraint due to multiple births) and genetic origin in midline craniosynostosis (i.e. metopic and sagittal). Recently, an in vitro model of forceinduced sagittal craniosynostosis has corroborated the importance of intrauterine constraint in initiation of midline craniosynostosis [START_REF] Oppenheimer | Force-induced craniosynostosis in the murine sagittal suture[END_REF].

Advanced maternal age (35 years and older) and birth weight (>4,000 g) are additional environmental components of sagittal NSC [START_REF] Boulet | A population-based study of craniosynostosis in Metropolitan Atlanta[END_REF].

Little is known about the genetic basis of NSC. FGFR 1, 2, 3 and TWIST are among the genes commonly identified for their involvement in craniosynostosis syndromes (Apert, Crouzon, Pfeiffer, and Saethre-Chotzen syndromes) [START_REF] Muenke | A common mutation in the fibroblast growth factor receptor 1 gene in Pfeiffer syndrome[END_REF][START_REF] Hollway | Mutation detection in FGFR2 craniosynostosis syndromes[END_REF]), but a role for these genes in NSC has not yet been identified [START_REF] Sakai | Sequence analysis of fibroblast growth factor receptor 2 (FGFR2) in Japanese patients with craniosynostosis[END_REF][START_REF] Boyadjiev | Analysis of candidate genes for non-syndromic craniosynostosis[END_REF][START_REF] Zeiger | Genetic and environmental risk factors for sagittal craniosynostosis[END_REF][START_REF] Morriss-Kay | Growth of the normal skull vault and its alteration in craniosynostosis: insights from human genetics and experimental studies[END_REF][START_REF] Boyadjiev | Genetic analysis of non-syndromic craniosynostosis[END_REF].

Sagittal NSC is thought to represent a consistent phenotypic group of unknown aetiology (see for example the 2D quantitative study of [START_REF] Guimarães -Ferreira | Sagittal synostosis: I. Preoperative morphology of the skull[END_REF]) that shares a common skull phenotype described as scaphocephalic even if variation in calvarial shape has been reported [START_REF] Schmelzer | Identifying reproducible patterns of calvarial dysmorphology in nonsyndromic sagittal craniosynostosis may affect operative intervention and outcomes assessment[END_REF].

Though current research efforts are focused on identifying subtle associated phenotypes in other developmental systems [START_REF] Boyadjiev | Genetic analysis of non-syndromic craniosynostosis[END_REF], NSC is still thought of as primarily a craniofacial disease.

Since the advent of CT imaging, detailed analysis of the three-dimensional morphology of both the ecto-and endocranial surfaces of the skull has been possible and much has been learned about the facial and cranial base correlates of premature neurocranial suture closure. However, most quantitative analyses have been limited to traditional anthropometric measures or linear measures between anatomical landmarks. Anatomical landmarks are precise locations on biological forms that hold some biological meaning or significance and ensure correspondence among the forms being compared [START_REF] Richtsmeier | Precision, repeatability and validation of the localization of cranial landmarks using computed tomography scans[END_REF]. Landmark locations are recorded in 2-or 3D space with an estimated degree of precision and accuracy. Because the cranial vault is relatively devoid of anatomical landmarks, especially in comparison to the density of landmarks on the cranial base and the facial skeleton [START_REF] Valeri | Capturing data from threedimensional surfaces using fuzzy landmarks[END_REF], information pertaining to shape change of the cranial vault is limited when only anatomical landmarks are used. Here, we introduce semilandmarks to the investigation of craniofacial shape of infants with sagittal craniosynostosis.

Semilandmarks were introduced to adapt landmark-based statistics to smooth curves in 2D [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF][START_REF] Bookstein | Landmark methods for forms without landmarks: Morphometrics of group differences in outline shape[END_REF], and the algebra was extended to accommodate the analysis of 3D data [START_REF] Gunz | Statistical and geometric reconstruction of hominid crania: Reconstructing australopithecine ontogeny[END_REF]Gunz, et al. 2005). Semilandmarks present "deficient" coordinates and require specific computational steps to define their final location on their defined curves or surfaces. After this computational stage, semilandmarks acquire a geometric correspondence across individuals so that comparative analyses can be conducted. The use of semilandmarks allows a more complete characterization of the overall morphology of the cranial vault, a region of high interest in the study of NSC.

Here we investigate skull shape in a sample of infants with sagittal NSC using anatomical landmarks and semilandmarks and appropriate methods of analysis to identify and quantify phenotypic variation in sagittal NSC. Using these data we sequentially test three different hypotheses.

Research Design

Since [START_REF] Virchow | Ueber den cretinismus, namentlich in Franken, und ueber pathologische Schädelformen[END_REF] study, an implicit assumption of a direct causal relationship between suture closure and cranial vault shape has underlain most craniosynostosis research and treatment protocols. Here we use cross sectional data of a sample of individuals with sagittal NSC who had not yet undergone any surgery or manipulation to test hypotheses about the relationship between suture closure and cranial shape.

A generally accepted idea is that the earlier the synostosis takes place, the greater the effect on skull shape and, conversely, the later that synostosis occurs, the less effect on skull shape [START_REF] Cohen | Perspectives on craniosynostosis[END_REF]. In the present study, no control sample is considered since our aim is to better characterize sagittal NSC and not to compare its phenotype with the unaffected skull phenotype. Because we did not know the timing of onset of craniosynostosis for individuals in our sample and we did not use a control sample, the relationship between timing of craniosynostosis and effect on skull shape, as well as the notion of severity, could not be tested in the present manuscript. However, the relationship between shape variation and age and sex represent key issues that warrant study. Consequently we test the hypothesis that cranial shape is associated with sex, age, or level of maturity (H1). Support for H1 will be characterized by the identification of different groups of individuals that share common shape features according to their sex, age or level of maturity. [START_REF] Albright | Suture pathology in craniosynostosis[END_REF] studied the timing of sagittal suture closure in 14 NSC patients. Their results revealed a heterogeneous sequence of closure patterns suggesting variability in fusion of the sagittal suture in sagittal NSC cases. In the present study, we focus on the pattern of closure of the sagittal suture of 43 patients with sagittal NSC. We divide the sagittal suture into three anatomic sections (anterior, central, posterior) and each of these sections is coded to characterize the nature of suture closure (i.e. closed, fused). In the present study a stepwise and irreversible progression is assumed between the different stages of fusion: initially, the suture is patent, then it is closing, and finally it becomes fused. When the sagittal suture is "fused", no evidence of a suture between the two parietal bones is visible on 3D CT reconstructions. When the sagittal suture (or part of it) is "closed" but not yet fused, the suture remain visible as a fine line. Note that a single sagittal suture can present a closure pattern involving different stages (e.g. a sagittal suture fused in its central and anterior sections but closed in its posterior section). We refer to the resulting set of observations as the "closure pattern". During craniofacial development, the sagittal suture will present different closure patterns. It is then possible to study the sequence of these closure patterns. Sequence as used in this paper is different from timing and does not imply any reference to absolute time. The timing of a specific event (e.g. sagittal suture fusion) includes the moment when this event starts and ends, and its duration. Even if results of cross sectional studies based on large sample can be extrapolated to describe the timing of an event, only longitudinal studies provide the information required to define the timing of an event. In the present crosssectional study we only describe different sequences of suture closure patterns and do not focus on timing. From our data on closure pattern, we test a second hypothesis that a unique sequence of sagittal suture closure patterns is associated with sagittal NSC (H2). [START_REF] Virchow | Ueber den cretinismus, namentlich in Franken, und ueber pathologische Schädelformen[END_REF] proposed that cranial deformity found in subjects with craniosynostosis could be explained by growth inhibition at right angles to the fused suture, with compensatory overexpansion of the cranium at patent sutural sites accommodating normal growth of the brain. Since then, this concept of dysplastic and compensatory growth has been shown to involve not only the calvaria, but all regions of the craniofacial complex in patients with different types of craniosynostosis [START_REF] Kreiborg | Postnatal growth and development of the craniofacial complex in premature craniosynostosis[END_REF]. Though a heterogeneous aetiology has been suggested, [START_REF] Boyadjiev | Genetic analysis of non-syndromic craniosynostosis[END_REF] samples of sagittal NSC individuals have implicitly been considered as representing a single, albeit variable, calvarial phenotype (scaphocephaly). Some have tried to quantitatively address the phenotypic variation of NSC [START_REF] Richtsmeier | Precision, repeatability and validation of the localization of cranial landmarks using computed tomography scans[END_REF] with limited success. Difficulties stem from the relatively small sample sizes, the varied quality of image data, and the specific quantitative approaches used in analysis. Here we take full advantage of available data and recently developed tools to quantify shape variation and investigate in 3D detail the potential relationship between the closure patterns sequence of the sagittal suture and corresponding skull shape. We test a third hypothesis that different sagittal suture closure pattern sequences generate different craniofacial shapes (H3).

Materials and methods

Images

Computed tomography (CT) images of children diagnosed with sagittal NSC were acquired from several medical centers (John Hopkins Hospital; Children's Hospital -Boston; University of Oklahoma Medical Center; Washington University), some of which participate in the International Craniosynostosis Consortium (genetics.ucdmc.ucdavis.edu). These images were deposited into the Image Analysis and Morphometrics Laboratory archive housed at the Pennsylvania State University (getahead.psu.edu). Here we analyze data from the images of 43 individuals diagnosed with sagittal NSC. Individuals diagnosed with a known genetic syndrome were excluded and none of the individuals included in our study sample had extracranial anomalies. For 22 out of 43 individuals, a detailed clinical genetics evaluation was done by hot-spot mutation analysis of the areas associated with known craniosynostosis syndromes (FGFR1 exon IIIa, FGFR2 exons IIIa and IIIc, FGFR3 exon IIIa, and the entire coding sequence of TWIST). None of these 22 individuals carried any of the known mutations. Our sample consequently consists of 22 individuals who did not express any of these known mutations and 21 individuals with a clinical diagnosis of sagittal NSC for who molecular data was not available. At the time of CT exam, infants had not undergone any surgical procedure and ranged in age from 0.9 to 9.0 months (mean age: 4.7 months, sd: 2.2 months) (Table 1).

Infants were classified into three-month age intervals for subsequent analyses. The sample includes a majority of males (i.e. 35 males, 8 females), which reflects the male to female ratio reported for sagittal synostosis (e.g., 7:3 to 4:1 [START_REF] Shilito | Craniosynostosis. A review of 519 surgical patients[END_REF][START_REF] Hunter | Craniosynostosis. I. Sagittal synostosis: its genetics and associated clinical findings in 214 patients who lacked involvement of the coronal suture(s)[END_REF][START_REF] Lajeunie | Genetic study of scaphocephaly[END_REF].

Closure of the sagittal suture

Because individuals within this sample were medically diagnosed with sagittal NSC (criteria for inclusion in the present study), all have premature fusion of the sagittal suture. However, we found that the exact condition of the sagittal suture is not identical for all individuals. We scored the condition of the sagittal suture as either totally or partially fused (Fig. 1). When the sagittal suture is 'fused', no evidence of a suture between the two parietal bones is visible on 3D CT reconstructions. When the sagittal suture (or part of it) is 'closed' but not yet fused, the suture remains visible as a fine line. Note that "closed" as used in the present study is sufficient to diagnose sagittal craniosynostosis.

In the present study, we divided the sagittal suture into three sections based on anatomical location: anterior, central, posterior, and each section was scored as fused (F) or closed (C). Infants with a completely fused sagittal suture were coded as "FFF" (i.e. the anterior, central, and posterior sections of the sagittal suture are fused). Infants whose sagittal suture was fused anteriorly but closed posteriorly were coded as "FFC" or "FCC" according to the status of the central section of the suture.

Finally, infants whose sagittal suture was closed anteriorly and fused posteriorly were coded as "CFF" or "CCF" according to the status of the central section of the suture (Table 1).

Dental development stage assessment

Chronological age of each individual at CT exam was computed by determining the number of days between the date of birth and the date of the CT exam. However, since all of our subjects are infants, it is possible that chronological age may be biased because our sample might include infants born prematurely or who experience some sort of developmental or growth delay. Unfortunately we do not have access to data regarding gestational age for all the individuals in our sample.

To address this potential bias and to determine a degree of maturity for each individual in our sample, dental development stage of each infant was assessed using the deciduous maxillary first molar (m 1 ). We adopted the radiographic method originally proposed by [START_REF] Demirjian | A new system of dental age assessment[END_REF] for permanent teeth and transposed to deciduous teeth by [START_REF] Liversidge | Variation in crown and root formation and eruption of human deciduous teeth[END_REF] as adapted to 3D CT reconstructed deciduous teeth by [START_REF] Bayle | Dental developmental pattern of the Neanderthal child from Roc de Marsal: a high-resolution 3D analysis[END_REF]. Although eight stages are defined from the first signs of mineralisation of the cuspids (stage A) to the completion of the root(s) (stage H1), the relatively young age range of our sample limits us to the crown formation stages (A to D) (Fig. 2). We used thresholding tools to segment tooth anatomy and therefore most likely underestimated the degree of dental mineralisation. However since we used the same tool for every individual, this potential underestimation is consistent across individuals in our sample. In the present study, no attempt was made to link the degree of maturity assessed by means of dental stages to chronological age. Consequently, we do not refer to any biological or developmental age, but only to a certain degree of maturity.

Landmark data collection

Images were reconstructed using a threshold that enabled visualization of bone. A set of 33 anatomical landmarks (LM) were defined and located on the 3D reconstruction of the CT images of each infant and the corresponding x,y,z coordinates were recorded (Table 2, Fig. 3). The 3D reconstructions of the CT images were manually closed using a segmentation tool at the region of the patent fontanelles to allow the placement of landmarks in these regions. In addition to the anatomical landmarks, a total of 189 semilandmarks were defined on four predefined curves (53 curve semilandmarks; CL) and four surface patches (136 surface semilandmarks; SL) on each skull. In order to gain geometric correspondence, semilandmarks were "slid" or moved along these curves and patches according to a sliding algorithm that minimizes the bending energy [START_REF] Bookstein | Landmark methods for forms without landmarks: Morphometrics of group differences in outline shape[END_REF]; for a more detailed description of the definition, recording and treatment of the semilandmarks see supplementary data 1).

Superimposition and exploration of the shape space

The 43 individuals defined on the basis of 222 landmarks (anatomical and semilandmarks combined) were superimposed by translating, rotating and scaling all forms with the aim of reducing the sum of the squared distances between homologous landmarks by means of a general Procrustes analysis [START_REF] Rohlf | Extensions of the Procrustes method for the optimal superimposition of landmarks[END_REF]. The coordinates of the resulting centered, scaled, and rotated landmarks are called the Procrustes shape coordinates. A measure of overall size called centroid size (CS; the square root of the sum of the squared distances of the landmarks to the centroid) was estimated for each individual and used as a proxy for overall cranial size in subsequent analyses. Centroid size is a measure of scale for landmark configurations, which has been shown to be approximatively uncorrelated with shape for small isotropic variation [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF]Dryden & Mardia, 1998). A Procrustes average shape (PAS) is computed as the coordinate-wise average of the Procrustes shape coordinates.

The Procrustes shape coordinates were analyzed by a principal components analysis to reduce the dimensionality of the dataset and explore the placement of individuals within the shape space. The eventual goal is to study the specific combination of morphological variables that successfully separate individuals into groups of known membership by projecting them onto the shape space. PCA is computed by an eigendecomposition of the sample covariance matrix and is a rigid rotation of the data preserving the Procrustes distances among individuals. The eigenvectors (or principal components, PCs) contain the weightings for the linear combination of the original variables and can be visualized as actual shape deformation [START_REF] Mitteroecker | Advances in geometric morphometrics[END_REF]. A simulation of a continuous shape deformation based on the available data can be obtained by warping the PAS to the negative (or positive) direction of the principal component. This is done simply by adding a multiple of the eigenvector to the mean shape (PAS). A warping through the negative direction will produce a deformation that is the inverse of a warping through the positive direction. The warping was executed by thin-plate spline (TPS; [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF].

Software

Amira 4.0.1 (Mercury Computer Systems Inc.) was used to visualize the CT images, to segment and reconstruct the skulls, to measure anatomical landmarks, and to do surface warping using TPS [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF]. The surface semilandmarks were measured on the template in Landmarks v3.0.0.6 [START_REF] Wiley | Evolutionary morphing[END_REF].The semiautomatic procedure allowing the measurement of the semilandmarks on the target skulls using TPS was done in Viewbox 4.0.0.98 (dHAL software, Athens, Greece). The Procrustes superimposition and subsequent statistical analysis were done using MorphoJ 1.01a [START_REF] Klingenberg | MorphoJ. Faculty of Life Sciences[END_REF].

Results

The PCA of the Procrustes shape coordinates using all landmarks and semilandmarks of the 43 infants with sagittal NSC was computed. No separation of the data according to sex or age (grouped by 3-month intervals) was observed with the first 10 PCs accounting for a total of 82% of the total variance (see Fig. 4 for a plot of PC1 and PC2). A lack of separation among groups for the first 10 PCs persisted when the dental stage of the deciduous first molar (m 1 , Table 1) was substituted for the 3-month age groups in the PCA (not shown). No separation of individuals according to the availability of molecular data was observed with the first 10 PCs (not shown). Multivariate regression between CS (size; independent variable) and the first 10 PCs (shape; dependent variables) indicated that size only accounts for 4.5% of shape variation (percentage predicted = 1 -[sums of squares error/sums of squares total]). A permutation test (10,000 rounds) against the null hypothesis of independence between size (CS) and shape (first 10 PCs) gives a P-value of 0.053. Similar results were obtained when size was replaced by age (percentage predicted: In contrast to sex, chronological age and m 1 dental stages, the sagittal suture closure patterns allowed a separation of the 43 individuals according to skull shape as estimated by the PCA of the Procrustes shape coordinates (Fig. 6). Indeed, the individuals coded as FFC showed negative PC1 scores, while the individuals coded as CFC and CFF showed positive PC1 scores. Though the demarcation was not as clear as with PC1, individuals coded as CFC (negative scores) are separated from those coded as CFF (positive scores) along PC2.

Exploration of the shape continuum associated with the negative and positive PC scores allows the observer to better understand the information expressed by the PCA. The main differences between a warped skull corresponding to -0.05PC1 and a warped skull corresponding to +0.04PC1 were mainly located on the cranial vault and the cranial base. Few differences were noted for the face (Fig. 7). The anterior frontal of a warped skull corresponding to +0.04PC1 is supero-inferiorly elongated and protrudes anteriorly above glabella (Fig. 7,arrows A,E). The bregma is relatively "higher". The occipital is convex and more inferoposteriorly projected with a "lower" lambda (Fig. 7, arrow C). In lateral profile, the posterior sections of the parietals are more inferiorly positioned and formed a more concave profile (Fig. 7, arrow B). In a supero-inferior view, the skull appears more elongated in the anteroposterior direction and the medial and posterior crania fossae are narrower (Fig. 7, arrows F, G). The cranial base, especially the middle and posterior fossae, is more posteriorly shifted, as seen with the petrous temporal bones and the occipital (Fig. 7, arrows C, D, F).

Regarding PC2, the main differences between a warped skull corresponding to -0.04PC2 and a warped skull corresponding to +0.04PC2 were located on the face, the cranial base, and the posterior section of the cranial vault (parietal and occipital bones) (Fig. 7). A warped skull corresponding to +0.04PC2 demonstrates a wider upper face, wider orbits associated with a broader interorbital distance, and a relatively less anteriorly projected face. The cranial base, especially the anterior and middle cranial fossae, is more inferiorly projected, or depressed (Fig. 7, arrow J). The anterior and middle crania fossae are wider while the posterior crania fossa is narrower (Fig. 7, arrow K). The posterior section of the parietals and the occipital are less posteriorly projected but present a similar profile to that of a warped skull corresponding to -0.04PC2 (Fig. 7, arrows H, I, L).

Discussion

Deformation of the cranial vault has been observed as a defining characteristic in craniosynostosis since its definition. Qualitative assessment of calvarial dysmorphology has shown consistent variability in sagittal NSC [START_REF] Schmelzer | Identifying reproducible patterns of calvarial dysmorphology in nonsyndromic sagittal craniosynostosis may affect operative intervention and outcomes assessment[END_REF]). Quantitative analysis of the cranial vault has proven problematic due to its smooth contours and relative lack of characteristic locations that can be used to precisely place anatomical landmarks. Using geometric morphometric methods developed for the placement and analysis of semilandmarks, we have provided 3D information on cranial shape variation in sagittal craniosynostosis including detailed information about vault shape. As a consequence, we are able to study the relationship of premature fusion of the sagittal suture with the overall shape of the cranium.

H1: Cranial shape in sagittal NSC is associated with sex, age, or level of maturity

When sex, chronological age, and level of maturity assessed by means of m 1 were used as classifiers with the PCA of the Procrustes shape coordinates, no obvious separation of the data was seen. Multivariate regression between age and shape (first 10 PCs) indicated that age only accounts for 4.5% of shape variation.

Hypothesis 1 is not supported by these results which imply that shape variation as registered by the first 10 PCs was not associated with sex, age or level of maturity. In addition, multivariate regression between size and shape (10 first PCs) showed that size only accounts for 4.5% of shape variation, while a PCA of the residuals of the regression of shape with size produced almost identical placement of the individuals to the one obtained in the shape space. Consequently, allometry was not the main source of shape variation. Although we cannot address the "severity" of shape dysmorphology in this sample of sagittal NSC due to the lack of a control sample, our analysis demonstrates that neither chronological age nor level of maturity explain the shape variation quantified for this sample of 43 children with sagittal NSC. Our results may reflect the limited age range of this sample (i.e. 0.9 to 9 months).

Because the great majority of reconstructive surgeries are done before the age of one year, we were not able to study the influence of sagittal NSC on the shape of the craniofacial complex after the first year of life.

H2: A unique sequence of sagittal suture closure patterns is associated with sagittal NSC

By coding closure pattern of the sagittal suture we gain access to the sequence leading from a patent sagittal suture to a fused one. Our analysis depends on the principle that a stepwise and irreversible progression exists between the different stages of fusion defined in the present study: initially, the suture is open, then it is closing (C), and finally it becomes fused (F). Careful observations show that different stages of fusion can be found simultaneously across the sagittal suture. For this reason we divided the sagittal suture into three distinct sections and qualitatively assessed and scored each section. Since the central section of the sagittal suture is the region that was most often fused (95.6%, Table 1), we suggest that it is frequently the section of the suture that closes first. This confirms previous observations made by [START_REF] Vannier | Radiologic evaluation of craniosynostosis[END_REF]. The co-occurrence of the sagittal suture closure patterns CFF and FFC in roughly equal numbers (Table 1) does not support the hypothesis of a unique sequence of closure patterns since by definition closure (C) cannot follow fusion (F). The equal mean sizes of CFF and FFC (Fig. 5) also indicate the coexistence of multiple sequences of closure patterns. At least two different paths towards complete fusion (FFF) appear to exist in the sample. Either the anterior section will be the second part of the sagittal suture to fuse (FFC), or the posterior section will be the second part to fuse (CFF). A third sequence of closure patterns could exist if the anterior and posterior sections of the sagittal suture fused roughly in the same time (Fig. 8). The coexistence of at least two different sequences of closure patterns for the fusion of the sagittal suture in sagittal NSC does not support the hypothesis of a unique sequence of closure patterns (H2) and confirms results of a previous analysis of a small sample which showed various closure patterns of the sagittal suture [START_REF] Albright | Suture pathology in craniosynostosis[END_REF].

H3: Different pattern sequences of sagittal suture closure generate different craniofacial shapes

Of the variables considered (sex, chronological age, dental stages of m 1 and sagittal suture closure pattern), suture closure pattern is the only variable associated with placement of individuals on the plot of the PCA of the Procrustes shape coordinates of the 43 individuals (Fig. 6). Skull shapes that characterize the negative end of PC1 were associated with an expanded temporal width, a relatively "low" bregma, and a relatively flat occipital with a "high" lambda. In addition to the individuals showing the pattern FFF who were homogeneously distributed, the individuals plotted on the negative end of PC1 share the sagittal suture closure pattern FFC. Following the concept of dysplastic and compensatory growth of the craniofacial complex, we interpret the expanded temporal width, the superoposteriorly expanded parietals and the relatively "low" bregma and "high" lambda as corresponding with premature fusion of the median and anterior sections of the sagittal suture. In this interpretation, the posterior section of the sagittal suture is the last section to close, remaining patent and thereby experiencing the maximum of local compensatory growth. Individuals with positive PC1 scores are associated with a skull shape characterized by relatively superiorly placed bregma and a relatively inferiorly placed lambda, a relatively more convex and more inferoposteriorly projected occipital, a narrower posterior half of the skull, and a posteriorly and inferiorly displaced cranial base (Fig. 7). In addition to those individuals showing the pattern FFF who were homogeneously distributed, individuals situated on the right side of the PCA plot (i.e. those with positive PC1 scores) shared the patterns CFC and CFF. Finally, the posterior section of the metopic suture was largely patent in the majority of individuals (15 individuals out of 17) showing the closure patterns CFC and CFF (positive PC1 scores). In these cases, the anterior fontanelle took on a triangular or a more rectangular shape (as seen in Fig. 3). The individuals showing the closure patterns FFC (negative PC1 scores) presented a posterior section of the metopic suture that was largely patent in only 2 out of 7 cases. At least two different interpretations can be made about the presence of this large patent area: (1) the patent area is likely associated with the localized overgrowth of the frontal region; or (2) the patent area potentially reflects the compensation for increased intracranial pressure due to the constraint on calvarial width expansion.

Consideration of PC2 reveals a separation among the suture closure patterns CFC and CFF. Individuals with positive PC2 scores correspond to the pattern CFF and the corresponding warped skull (+0.04PC2) presented a relatively wider face and frontal portion of the skull (but a narrower occiput), a relatively reduced anteroposterior dimension, and a cranial base that is relatively more inferiorly placed (Fig. 7). Following the concept of dysplastic and compensatory growth of the craniofacial complex, we suggest that the individuals at the positive end of PC2 reveal a craniofacial shape that is a consequence of the premature fusion of the middle and posterior sections of the sagittal suture, which generates a physical constraint driving a compensatory response mainly located in the anterior part of the skull and along the supero-inferior axis.

Instead of using skulls warped according to the negative or positive direction of the eigenvector, it is also possible to compute and then compare Procrustes showing the closure pattern FFC. Because the PCA of the Procrustes shape coordinates of the 43 individuals allows for good separation between individuals according to their closure pattern, this alternative approach allows similar observations than those made with the warped skulls. Supplementary data 2 and 3 present this alternative approach based on the PASs.

Since [START_REF] Virchow | Ueber den cretinismus, namentlich in Franken, und ueber pathologische Schädelformen[END_REF], several studies on craniosynostosis have produced results supporting the concept of dysplastic and compensatory growth of the craniofacial complex (e.g. [START_REF] Kreiborg | Postnatal growth and development of the craniofacial complex in premature craniosynostosis[END_REF]. [START_REF] Virchow | Ueber den cretinismus, namentlich in Franken, und ueber pathologische Schädelformen[END_REF] thought the primary cause was the premature fusion of a cranial suture and that cranial base deformation followed. [START_REF] Moss | Growth of the calvaria in the rat: The determination of osseous morphology[END_REF][START_REF] Moss | Functional anatomy of cranial synostosis[END_REF] proposed the functional matrix hypothesis which introduced the role of soft tissues in cranial growth. He argued that cranial base dysmorphology result in aberrant tensile forces being communicated from the brain to the calvarium via dural attachments, ultimately causing craniosynostosis. Although our study does not support either of these theories in particular, it shows a strong relationship between the sequence of the closure pattern of the sagittal suture and deformation of the cranium. Our results add to the well recognized relationship between craniofacial shape and premature suture closure. For premature closure of the sagittal suture, our results suggest that the shape of the craniofacial complex not only identifies the sagittal suture as the suture that is closed prematurely but also the pattern in which the suture closes.

By analyzing variation in skull shape in a cross sectional sample of children with sagittal NSC we participate in a more general effort aimed at the understanding of the aetiology of this disorder and thereby adding to the formulation of clinical rules for craniofacial surgery and neurosurgical treatment that contribute to skull shape (e.g. [START_REF] Delashaw | Rules for cranial vault growth[END_REF][START_REF] Marsh | Surgical research on craniosynostosis[END_REF]. Multiple surgical techniques are available for the treatment of sagittal NSC. It is most important, that the surgical approach be suited to the age of the patient, the state of suture fusion and the particular deformity.

Besides the general scaphocephaly, anterior or posterior bulging of the skull and temporal indentation can be major problems that require appropriate variation in technique [START_REF] Jane | Neurosurgical treatment of craniosynostosis[END_REF]. Importantly, none of the 43 individuals (including our youngest individual aged about one month) showed a section of the sagittal suture still patent. Consequently, the closure of at least the central section of the sagittal suture most likely takes place prenatally, suggesting that sagittal NSC is a congenital malformation. We suggest that prenatal ultrasound could be used for third trimester prenatal diagnosis which could aid the family and the medical team for optimal monitoring of the patient.

Our results suggest a relationship between pattern and sequence of suture closure in sagittal NSC and calvarial dysmorphology. Tests of our results require different and possibly large samples where factors known to potentially influence skull shape are known (e.g. ethnicity, duration and condition of labour). Whether the patterns we have quantified indicate differences in aetiology cannot be determined with our data, but require longitudinal studies, most appropriately of animal models where prenatal conditions can be manipulated and monitored. semilandmark. By looking at the figure 3, one can easily imagine the grids used in this study. Because different patches share one or several edges, as a final step repetitive semilandmarks were manually deleted.

Although anatomical landmarks correspond across individuals due to commonalities of development and anatomy, this is not the case for semilandmarks.

Still, the defined curves and/or surfaces on which the semilandmarks lie correspond across individuals (e.g., the right orbit curve of a skull A corresponds to the right orbit curve of a skull B) making them useful in determining similarities and differences across shapes.

Contrary to classic anatomical landmarks, semilandmarks present "deficient" coordinates. For instance, a curve semilandmark will have well-defined coordinates on two axes but the coordinate on the third axis depends on the defined curve. The specific location of semilandmarks is adjusted by allowing them to "slide" along the defined curve (curve semilandmarks) or surface (surface semilandmarks) until a specified measure of shape difference among the configurations is minimized [START_REF] Bookstein | Landmark methods for forms without landmarks: Morphometrics of group differences in outline shape[END_REF][START_REF] Gunz | Semilandmarks in Three Dimensions[END_REF]. In actuality, to linearize the problem, the semilandmarks do not slide on the defined curve or surface, but instead on the tangent vector to the curve or the tangent plane to the surface. Semilandmarks "slide" along the defined curve or surface through an iterative process repeating the following three steps as outlined by [START_REF] Mitteroecker | Advances in geometric morphometrics[END_REF]: (1) computation of a tentative sample mean shape and of the tangent vectors for each semilandmark in each individual, (2) sliding of semilandmarks along the tangents to minimize "a specified measure of shape difference" to the mean, and (3) placing back the slid semilandmarks to the nearest point on the curve or the surface, as they may slip off the curvature.

The sliding process can be based on one of two minimization criteria: one algorithm is based on minimizing the "bending energy" [START_REF] Bookstein | Landmark methods for forms without landmarks: Morphometrics of group differences in outline shape[END_REF], while the other, based on a Procrustes algorithm minimizes Procrustes distances [START_REF] Bookstein | Corpus callosum shape and neuropsychological deficits in adult males with heavy fetal alcohol exposure[END_REF]. Application of either of these criteria causes the semilandmarks to slide along the defined curves or surfaces until they are placed in a location consistent with minimization of the chosen criteria. In the presence of limited shape variation the differences between the two minimization criteria have little effect.

The minimization calculation that determines the sliding of the semilandmarks to their final location confers a correspondence among the forms being analyzed that is both geometric due to the calculations used and biological due to the biological correspondence of the curves and surfaces under consideration. Once this correspondence is established, semilandmarks are analyzed using the same procedures adopted for anatomical landmarks in the subsequent statistical analyses.

Methods for determining the position of semilandmarks are designed to operate on smooth curves and surfaces. These methods become inappropriate as soon as the relative smoothness of the structure of interest diminishes.

Once the semilandmarks are located on the template skull, a semiautomatic procedure is utilised to place the semilandmarks on each target skull. This procedure relies on a warping of the template skull onto the target skull using TPS [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF]. As stated earlier, the 3D coordinates of the 33 anatomical landmarks as well as the location of the four curves are recorded by hand. A TPS is then computed between the template and the target using only the anatomical landmarks (Fig. 3 blue points) to warp the curve (Fig. 3 -red points) and surface semilandmarks (Fig. 3 -green points) from the template to the target. Finally, a new TPS based on all landmarks and semilandmarks is computed. At this point semilandmarks are slid according to the chosen sliding process. Here we used the sliding algorithm that minimizes the bending energy [START_REF] Bookstein | Landmark methods for forms without landmarks: Morphometrics of group differences in outline shape[END_REF]. This semiautomatic operation is repeated for each target skull of the sample.

Placement of the individuals on PC1 and PC2 in the shape space defined by the principal components analysis of the Procrustes shape coordinates using all landmarks and semilandmarks of the 43 infants with sagittal NSC. Individuals are coded for suture closure pattern. The key is the same as in Fig. 6

  4.5; permutation test P-value: 0.058). A PCA of the residuals of the regression of the PCs coefficient with CS resulted in almost identical placement of the individuals as the one obtained with the PCA of the Procrustes coordinates. In only two individuals out of 43 was the central section of the sagittal suture closed but not yet fused. The youngest child of the sample shows the pattern FCC and another individual shows the pattern CCC. Of the 41 remaining individuals, each child exhibited complete fusion (F) of the central region of the sagittal suture. However, we found the following variation in other sections of the suture: 10 individuals showed the anterior and posterior sections closed (i.e. CFC); seven displayed the anterior section closed and the posterior section fused (i.e. CFF); and eight individuals exhibited the anterior section fused and the posterior section closed (i.e. FFC). Finally, the sagittal suture was totally fused (i.e. FFF) in 16 individuals. Considering the relationship between overall size (CS being highly correlated with age; R²=0.7628, p<0.001) and closure patterns, individuals showing the CFC closure patterns were significantly smaller (mean CS: 91.1) than: (1) individuals showing the pattern CFF (mean CS: 96.2) (t test, p= 0.036); (2) individuals showing the pattern FFC (mean CS: 96.3) (t test, p=0.019); and (3) individuals showing the pattern FFF (mean CS: 98.7) (t test, p=0.001). Individuals showing the CFF and FFC closure patterns were of equal mean sizes (Fig. 5).

  average shapes (PAS) according to the sagittal suture closure pattern. The PAS of all individuals showing the closure pattern CFC was computed, as were the PAS of all the individuals showing the closure pattern CFF, and the PAS of all individuals

  . Procrustes average shape (PAS) of the skulls with the closure pattern FFC is represented in blue. PAS of the skulls with the closure pattern CFC is represented in green. PAS of the skulls with the closure pattern CFF is represented in orange. Arrows indicate direction of possible compensatory growth related to suture sections which are closed (dashed lines) opposed to suture sections which are fused (full lines). Because PC1 and PC2 allow the separation of the individuals according to their closure pattern (CFC, CFF, and FFC), one can observe similarities between Procrustes average shape skulls represented here and corresponding warped skull represented in Fig. 7.

Fig. 1

 1 Fig. 1 Superior view of a sagittal suture (anterior fontanelle at top) illustrating sagittal suture closure pattern. Anterior section of sagittal suture "closed" (C), central and posterior section "fused" (F). Consequently, the closure pattern of the sagittal suture is CFF. The black arrow indicates appearance of suture sections we coded as "closed".

Fig. 2

 2 Fig. 2 Radiograph and 3D CT surface reconstruction of crown stages showing deciduous first molar. Stage A: beginning of mineralization at cusp tips, Stage B: coalescence of cusp tips to form a regularly outlined occlusal surface, Stage C: occlusal surface is complete. Approximal edges of forming crown has reached future contact areas, Stage D: crown (enamel) is complete with full-thickness occlusal dentine present, and roof of the pulp chamber is mature. Beginning of root formation is seen (modified after Liversidge & Molleson, 2004).

Fig. 3

 3 Fig. 3 Illustration of the 222 points measured on 3D reconstruction of computed tomography scans of each individual in our sample. Anatomical landmarks (Table 2) are shown in blue, curve semilandmarks are shown in red, and surface semilandmarks are shown in green. Top left: endocranial surface of the cranial base, top right: anterior view, bottom left: lateral view, bottom right: inferior view.

Fig. 4

 4 Fig.4Placement of the individuals on PC1 and PC2 in the shape space (principal components analysis of the Procrustes shape coordinates using all landmarks and semilandmarks of the 43 infants with sagittal NSC). Male infants are denoted by triangles while circles denote female infants. Ages of individuals are indicated by colour: white = 0-3 months, grey = 3-6 months, black = 6-9 months.

Fig. 5

 5 Fig. 5 Mean and 0.95 confidence interval plot of centroid size (CS) grouped by sagittal suture closure pattern (C: closed, F: fused).

Fig. 6

 6 Fig.6Placement of the individuals on PC1 and PC2 in the shape space defined by the principal components analysis of the Procrustes shape coordinates using all landmarks and semilandmarks of the 43 infants with sagittal NSC. Individuals are coded for suture closure pattern (closed: C, fused: F). Convex hulls (transparent) are drawn for the CFC, CFF and FFC groups.

Fig. 7

 7 Fig.7Superimposed warped skulls according to PC1 (left) and PC2 (right). Lateral views of endocranial surface viewed through a mid-sagittal section (top) and superior view of endocranial surface of axial sections (bottom). The blue warped skull corresponds to -0.05PC1 and the red warped skull corresponds to +0.04PC1. The green warped skull corresponds to -0.04PC2 and the orange warped skull corresponds to +0.04PC2. The colours of the warped skulls (i.e. blue, orange and green) correspond to the colours used in fig.6. Arrows show the main shape differences (see text for more details).

Fig. 8

 8 Fig. 8 Illustrations of three potential closure pattern sequences of the sagittal suture in sagittal NSC cases as documented in our study. The sagittal suture is divided into three distinct sections: anterior, median, and posterior. Each section is coded as closed (C) or fused (F).

Table 1

 1 Sample sex and age description. The closure pattern of the sagittal suture is indicated (C: closed, F: fused) as well as deciduous maxillary first molar (m 1 ) dental stages. Observations are grouped by sagittal suture closure pattern. Relative proportions and mean age for each group are presented.

	Suture closure pattern	Age	AgeGroup
	(mean age in months)	(months)	(months)	Sex	m 1	Percentage
	CCC	6.0	6-9	M	C	2.3
	FCC	0.9	0-3	M	B	2.3
		1.1	0-3	M	C
		2.0	0-3	M
		2.2	0-3	M	C
		2.5	0-3	F	B
	CFC (3.0)	2.6 3.1	0-3 3-6	M M	C B	23.3
		3.2	3-6	M	B
		3.6	3-6	M	B
		4.7	3-6	M	C
		5.4	3-6	M	B
		1.8	0-3	M	B
		3.3	3-6	F	C
		4.1	3-6	M	B
	CFF (4.9)	4.7	3-6	M	C	16.3
		5.2	3-6	M	D
		7.2	6-9	M	C
		8.0	6-9	M
		2.3	0-3	M	B
		2.3	0-3	M
		2.5	0-3	M	B
	FFC (4.6)	3.6 4.6	3-6 3-6	M F	C B	18.6
		5.2	3-6	M	C
		7.2	6-9	M	C
		8.7	6-9	F	C
		3.3	3-6	M	B
		3.4	3-6	M	C
		3.5	3-6	M	C
		3.5	3-6	F	B
	FFF (5.8)	4.0	3-6	M	B	37.2
		4.8	3-6	M	C
		5.0	3-6	M	C
		5.7	3-6	F	C
		6.1	6-9	F	D

Acknowledgments

We are grateful to all study participants and to the following persons who participated in the data collection and management: Daniel Govier from Washington University; Jayesh Panchal from Oklahoma City Children's Hospital; Craig Senders, Travis Tollefson and Bridget Wilson from University of California Davis; Virginia Kimonis from University of California Irvine; Nisha Isaac, Doris Lin, George Jallo, Ben Carson, and Craig Vander Kolk from the Johns Hopkins Medical Institutions; Caroline D. Robson and Joan Stoler from the Children's Hospital Boston; Satama Sirivunnabood from Pennsylvania State University. We thank Neus Martínez-Abadías for early fruitful discussion of this manuscript. We thank Brenda C. Frazier, Chris Percival and John M. Starbuck for editorial suggestions. We thank the two anonymous referees for their constructive comments that helped to enhance the overall quality of this manuscript. This study was funded by Public Health Service grants R01 DE018500, 3R01DE18500-02S1, R01 DE016886, and CDC 5R01DD000350. SAB is partially supported by Children's Miracle Network (CMN) endowed chair in Pediatric Genetics.

Supplementary Material

Supplementary data 1:

Semilandmarks definition and recording

To create the model for placement of semilandmarks, one individual of the sample was arbitrarily selected as the template skull. The template skull is the skull on which all the landmarks and semilandmarks were "manually" placed. Consequently, even though the choice of template skull is arbitrary, this choice will have an effect on the placement of landmarks on all other skulls (Fig. 3). The template was subsequently used in a semiautomatic procedure to determine the location of semilandmarks on all the remaining skulls of the sample called the target skulls. For the target skulls, only the 33 anatomical landmarks and the 4 curves were "manually" measured and recorded. The four curves recorded on each skull were:

1. an open curve describing the mid-sagittal plane (from nasale to opisthion passing by bregma and lambda) 2. two closed curves describing the right and left orbits (each starting and ending at the fronto-zygomatic suture and following the orbital ridge) 3. a closed curve describing the nasal aperture (starting and ending at the nasale, passing by the anterior nasal spine and following the nasal aperture ridge)

On the template skull, the 53 curve semilandmarks were distributed at equal intervals along the defined curves. The positions of the semilandmarks of the target skulls were computed by means of thin-plate spline algorithm (TPS, [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF]. Seventeen CLs were distributed along the mid-sagittal plane, 14 CLs were distributed around each of the orbit curves, and 8 CLs were distributed around the nasal aperture curve.

The method for determining the location of surface semilandmarks on the template is more complex. First, four surfaces patches were defined: the frontal patch, the left and right parietal-temporal patches, and the occipital patch. These patches were anchored by controlling points which were, when it was possible, actual anatomical landmarks (not necessarily among the 33 included in this study). Once each patch was defined, a grid was created and superimposed on the skull surface.

The intersection of perpendicular lines of the grid defines the position of a

Supplementary data 3

Superimposed Procrustes average shape (PAS) of the skulls according to their sagittal suture closure pattern. PAS of the skulls with the closure pattern FFC is represented in blue. PAS of the skulls with the closure pattern CFC is represented in green. PAS of the skulls with the closure pattern CFF is represented in orange. Note that a single skull with the closure pattern FFF has been warped to create the PASs skulls represented here. Consequently, the sagittal suture closure pattern of the PASs' representations do not correspond to the reality. The visible parts of the PAS correspond to the ones that are the closest to the eye of the observer. For example the orange PAS presents a wider and more anteriorly projected frontal than the blue PAS.