M Melich

J Dupont-Roc

Ph Jacquier

Search for anisotropic effects of hcp solid helium on optical lines of cesium impurities

Keywords: numbers: 67.80.Mg, 78.30.-j, 32.70.Jz, 32.80.Ys

The anisotropic effect of a hcp 4 He solid matrix on cesium atoms has been proposed as a tool to reveal the parity violating anapole moment of its nucleus. It should also result in splitting the D2 optical excitation line in a way depending on the light polarization. An experimental investigation has been set up using oriented hcp helium crystals in which cesium metal grains are embedded. Atoms are created by laser sputtering from this grains. Optical absorption spectra of the D2 line have been recorded in the temperature range of 1.0 to 1.4 K at liquid/solid coexistence pressure by monitoring the fluorescence on the D2 line at 950 nm. No significant effect of the light polarization has been found, suggesting a statistically isotropic disordered solid environment for the cesium atoms.

I. INTRODUCTION

Optical properties of heavy alkali metal atoms (Cs, Rb) embedded in liquid or solid helium have been studied since 1990 [1] and reviewed in several articles [2]. Liquid helium and b.c.c solid behave as isotropic matrices for the atomic impurities. On the contrary, uniaxial h.c.p. solid is expected to perturb the atomic levels as a quadrupolar hamiltonian, lifting the degeneracy between Zeeman sublevels with different values of |m|. Such a feature was indeed observed experimentally in the two ground state hyperfine sublevels [3] and was interpreted as an anisotropic hyperfine interaction induced by the crystal. This effect was later proposed as a tool to reveal the parity violating anapole moment of the cesium nucleus [4]. An electric and a magnetic field with well defined orientations with respect to the helium crystal axis should produce a linear Stark effect, proportional to this anapole moment. This experiment requires a uniform crystal orientation throughout the sample. This remains to be achieved, since the experiments in the Cs ground state quoted above were interpreted as resulting from a statistically isotropic distribution of the local crystal c-axis. Also unknown is the perturbation of the local crystal order around each Cs atom. Optical absorption spectra are a sensitive tool to test the isotropy of the local environment of the alkali atoms in the solid matrix. We show that for an anisotropic matrix the shape of the D 2 absorption line is expected to depend on the excitation light polarization, as well as that of the detected fluorescent light. In addition we report the result of an experimental investigation on cesium atoms sputtered in an initially mono-domain hcp crystal.

II. CESIUM IMPURITY ABSORPTION SPECTRUM IN LIQUID AND SOLID HELIUM

Optical transitions of a cesium atom embedded in a dense phase of helium have been described using various models [5]. Neglecting the hyperfine structure, the electronic energy levels reduce to the 8 sublevels shown in Fig. 1. The D 1 line corresponding to the 6S 1/2 →6P 1/2 transition is simply broaden by the interaction with the surrounding helium atoms. The properties of the D 2 line (6S 1/2 → 6P 3/2 transition) is due to the removal of the 4-fold degeneracy of the 6P 3/2 level by any anisotropic environment, resulting into two Kramers doublets. Let us first ignore fluctuations

1/3 2/3 1/3 2/3 -1/2 +1/2 m= 6P1/2 6S1/2 - - 2/3 1/3 1/3 2/3 6S1/2 6P3/2 -1/2 -3/2 +1/2 +3/2 1 1
FIG. 1: Energy levels of cesium with a quadrupolar perturbation. Hyperfine structure is ignored. Numbers shown are proportional to the electric dipole matrix elements. The energy gaps between the various energy levels are not to the same scale. Ground state is duplicated for clarity.

and assume that the environment has an axial symmetry, as the hcp crystal does. Then taking the z-axis along this direction, its action on the 6P 3/2 level may be written as H = A + B 3J 2 z -J(J + 1) /6. The energy sublevels are the m = ±3/2 and m = ±1/2 magnetic sublevels with respect to the z-axis, with energy A + B/2 and A -B/2 respectively. The matrix elements of the electric dipole moment responsible for the optical transitions are reminded in Fig. 1. They imply that the weights of the two absorption line components strongly depend on the light polarization with respect to the z-axis. For z-polarized light, only the |1/2| component is excited, while for x-or y-polarization, the 3/2 and 1/2 components are excited with the probability 1/2 and 1/6 respectively. Experimentally, the absorption of light is very weak. Since radiative decay is the dominant decay process for the 6P excited state, one measures instead the intensity in the atomic fluorescence band as a function of the excitation wavelength. This intensity depends also on fluorescence polarization. Table I gives relative line intensities for various polarizations of fluorescence light. It is implicitly assumed that relaxation in the excited state prior to fluorescence preserve the populations of the m sublevels. This amounts to neglect effects of the hyperfine coupling in the excited state. A fourth case is also considered which assumes on the contrary a complete depolarization of the excited state before fluorescence. It appears clearly that in any case the absorption spectrum strongly depends on excitation and detection polarizations. Consider now the effect of dynamic fluctuations. The effect of such a perturbation has been considered by Kinoshita et al. in the case of liquid helium [6], with no static anisotropy. Because of the three dimensions of space, an isotropic perturbation is less probable than anisotropic ones, and the line is split into two broad components of equal weights. Because of the statistical isotropy, no effect of polarization on the two component weights is expected. The more general case with a remaining static anisotropy is more complicated. It can be shown to reduce the contrast between the two components.

If the static anisotropy is not small with respect to the fluctuations, some effect of light polarization should remain. This motivated the following experimental investigation.

III. EXPERIMENTAL ARRANGEMENT AND PROCEDURE

The experimental cell is a 4 cm stainless steel cube, with five silica windows (diam. 2.5 cm), cooled from the top by a pumped helium-4 fridge in the 1.0-1.4 K temperature range. To achieve single crystal growth, the nucleation and growth is made at constant temperature and pressure (1.2 K and 25.5 bar). An electro-crystallization device [7] provides a unique seed which falls on the bottom of the cell. The crystal is subsequently grown from this seed. The orientation of the c-facet (most often horizontal) is easily monitored visually during the growth process which takes place below the corresponding roughning transition 1.3 K [8]. Cesium atoms are produced in two steps by laser sputtering from a metal target [1] located in the upper part of the cell filled with superfluid. The laser pulses (1-5 mJ, 8 ns, λ=532 nm) are provided by a Nd-YAG laser. The first step produces metallic grains, about 100 µm in size, which fall on the crystal surface, when the interface is located in the middle of the windows. These grains are then incorporated into solid helium by subsequent growth. In the second step, the laser is focused on these grains, and atoms are produced in their immediate vicinity. This production phase lasts between a few minutes to half an hour at a rate of 2 pulses/second with the aim of maximizing the fluorescence signal on the D 1 line. The D 1 and D 2 lines are excited by a tunable Ti:Sa CW laser scanned between 770 and 870 nm, with a z-or y-polarization, the z-axis being along the vertical c-axis of the grown crystal. The fluorescence light in the 950 nm band [9,10] is detected by a cooled APD and selected along the y-axis by a 900 nm-edge filter which blocks the stray light from the laser. The other 880 nm band is also occasionally used to monitor the D 1 fluorescence.

IV. RESULTS

A first observation is the time T a over which cesium atoms can be observed after their initial preparation. A very significant extension of T a was obtained when lowering temperature from 1.4 K to 1.0 K. While at 1.4 K, fluorescence from Cs atoms disappears in a few hours, as observed in Weis' group [11], at low temperatures it can be kept for days and even weeks. Because the diffusion of impurities in solid helium is expected to be thermally activated in this temperature range, such a behaviour is not unexpected. It is however of primary importance for future experiments with Cs atoms which will need long integration times.

Absorption spectra of the D 2 line are shown in Fig. 2 for various polarizations of excitation and detection. A splitting of the line into two components is clearly visible. Its magnitude corresponds to an energy gap of 200±10 cm -1 . The line shape is similar to what has been reported for liquid helium [6] and solid helium at higher pressure [12]. Quantitatively, line widths are about 15% larger than those in the liquid at 20 bar while splittings are 20% larger. The effect of the excitation polarization on the line shape is much smaller than expected. More precisely, the component weight ratio is 0.7 for z-excitation, z-analyzer and 0.85 for y-excitation, x-analyzer, to be compared with 0 and 9 expected respectively (see Table I).

V. DISCUSSION

There are two ways for explaining this result. The first one is to conclude that there is no permanent anisotropy in the environment of the atoms and that the splitting observed is mostly due to dynamic anisotropic fluctuations similar to those observed in liquid helium [6] as already discussed in reference [12]. The similarity of line shapes is indeed striking. Whether such an explanation is compatible with the hyperfine and magnetic resonance spectra observed in the ground state remains to be investigated. The second explanation is the one put forward in reference [3] : the sputtering process destroys the global order of the crystal which is transformed in a mixture of randomly oriented small hcp domains. Then no effect of the light polarization remains. If so, improvement of the crystal quality is to be searched for. Partial fusion, and re-crystallization from an ordered part of the solid is to be tried. It is unclear whether atoms will survive to such a process. Less damaging sputtering process using a femtosecond laser could also be attempted.

FIG. 2 :

 2 FIG. 2: Absorption spectra of the D2 line for different polarizations. (a) z-excitation, z-analyzer on detection, (b) y-excitation, x-analyzer on detection. Fitting the spectra with a sum of two gaussian lines yields the ones shown as plain curves.

TABLE I :

 I Intensities of the two components I 3/2 and I 1/2 of the D2 absorption line excited by a beam on the x-axis and detected on the fluorescence intensity in the y-direction for various polarizations.

Acknowledgments

We acknowledge support from ANR, grant 05-BLAN-0084-01.

[1] A. Fujisaki, K. Sano, T. Kinoshita, Y. Takahashi, T. Yabuzaki, Phys. Rev. Lett. 71, 1039 (1993) ; M. Arndt, S.I. Kanorsky, A. Weis, T.W. Hänsch, Phys. Lett. A 174, 298 (1993).


Antoine Boubault
email: antoine.boubault@mines-paristech.fr

Seungwoo Kang

Nadia Maïzi

Closing the TIMES integrated assessment model (TIAM-FR) raw materials gap with life-cycle inventories

Keywords: Raw materials, Integrated assessment models, Prospective scenario, TIMES, TIAM-FR, ecoinvent, Life-cycle inventories, Integrated assessment, Material flow analysis

Introduction

Following the recent development of large-scale integrated assessment models (IAMs), it is possible to investigate the complexity of anthropogenic energy systems, their interactions with Earth systems, and their multiscale evolution. Among the variety of existing models, global models stand out due to their wide scope, comprehensive overview and finite nature (resources come from the environment only, not as imports from other regions). The most authoritative projections, such as the IEA's World Energy Outlook (International Energy Agency 2016), the EIA's International Energy Outlook (U.S. Energy Information Administration 2017) and the WEC's World Energy Scenarios (World Energy Council 2016) provide global energy scenarios constrained by physical, technological, and economic drivers, as well as environmental policies. These models result from bottom-up, top-down or hybrid descriptions of the global energy system following optimized (function minimization) or simulated development trajectories that satisfy a set of constraints. Although these have become increasingly integrated and complex using large quantities of data, they typically do not consider nonfuel material availability, stocks, flows, transformations, functions and services.

World energy scenarios usually assume growing economic demand. This results in the development of industrial supply chains that require increasing primary energy and materials. In 2015, about 81% of world's primary energy production consisted in fossil energy carriers. More than a quarter (26.5%) of the world's primary energy supply (13.6 Gtoe) was directed to electricity and combined heat-and-power (CHP) plants in the form of fossil energy carriers (International Energy Agency 2017), as shown by Figure 1. This share increased steadily for decades (dotted line) until stabilizing recently.

In the context of sustainable development, removing fossil fuels from electricity-generating (EG) technologies will require the widescale development of renewable energy technologies. Multiple prospective studies deduce that more materials -especially metals -could be required in the coming decades following the projected development of low-carbon and energy-efficient technologies [START_REF] Van Vuuren | Long-term perspectives on world metal use -a system-dynamics model[END_REF][START_REF] Halada | Forecasting of the consumption of metals up to 2050[END_REF][START_REF] Gerst | Linking Material Flow Analysis and Resource Policy via Future Scenarios of In-Use Stock: An Example for Copper[END_REF][START_REF] Graedel | On the Future Availability of the Energy Metals[END_REF][START_REF] Kleijn | Metal requirements of low-carbon power generation[END_REF][START_REF] Keramidas | Future availability and demand for oil gas and key minerals[END_REF][START_REF] Elshkaki | Dynamic analysis of the global metals flows and stocks in electricity generation technologies[END_REF][START_REF] Vidal | Metals for a low-carbon society[END_REF]). These suggest that available material reserves and processing capabilities might not be sufficient to cover future material needs, driven by economic growth and environmental constraints. This perspective leads to two complementary questions: 1) What are the implications of long-term prospective scenarios for mineral resources? and 2) What are the impacts of a constrained mineral resource supply on long-term prospective scenarios?

Determining the answer to the first question would bring valuable knowledge about potential highlydemanded mineral resources in the future, and it could be used to achieve more efficient resource management [START_REF] Ali | Mineral supply for sustainable development requires resource governance[END_REF]) and more resilient material supply chains, thus mitigating stress periods.

Addressing the second question would provide elements to determine whether or not mineral resources should be considered in a given scenario. For instance, scenarios not taking into account adequate mineral resource supply could be irrelevant or might need to be reassessed in light of new resource supply constraints. Ideally, the two questions above are addressed together using a systems approach. From this perspective, a study of the mutual interaction between raw materials and prospective energy scenarios is relevant and may be analyzed using different industrial ecology tools like material flow analysis (MFA) and life-cycle assessment (LCA). In this direction, Hertwich et al.

performed a prospective integrated life-cycle assessment of global electricity generation to 2050, based on IEA's 2015 Energy Technology Perspective report [START_REF] Hertwich | Integrated life-cycle assessment of electricity-supply scenarios confirms global environmental benefit of low-carbon technologies[END_REF]. Their work revealed that LCA metrics carry crucial information that needs to be considered by both energy system modelers and climate policymakers. Currently, life-cycle metrics are rarely considered in integrated assessment models [START_REF] Pauliuk | Prospective Models of Society's Future Metabolism: What Industrial Ecology Has to Contribute[END_REF][START_REF] Pauliuk | Industrial ecology in integrated assessment models[END_REF], although some studies aim to cover this gap for greenhouse gas emissions (GHG) [START_REF] Daly | Indirect CO2 emission implications of energy system pathways: linking IO and TIMES models for the UK[END_REF][START_REF] Pehl | Understanding future emissions from low-carbon power systems by integration of life-cycle assessment and integrated energy modelling[END_REF][START_REF] Mcdowall | Is the optimal decarbonization pathway influenced by indirect emissions? Incorporating indirect life-cycle carbon dioxide emissions into a European TIMES model[END_REF]) and for other life-cycle impacts [START_REF] Arvesen | Deriving life cycle assessment coefficients for application in integrated assessment modelling[END_REF].

In this paper, we address question 1) above and make a first step towards representing mineral resource use in the TIMES integrated assessment model (TIAM-FR version). TIAM-FR is a bottomup integrated assessment model that was developed by the Energy Technology Systems Analysis Program (ETSAP) group of the IEA [START_REF] Loulou | ETSAP-TIAM: the TIMES integrated assessment model Part I: Model structure[END_REF]. It provides decision-makers with prospective multiregional energy scenarios. Its bottom-up description of the energy system allows us to access the life-cycle material footprint of each technology using life-cycle inventories (LCIs). At that point, we introduce a novel but straightforward life-cycle material accounting method to estimate the life-cycle material use of EG technologies. First, we describe the method that attributes LCIs to the TIAM-FR EG processes. In particular, the technologies' LCIs are disaggregated into three phases: construction, operation and decommissioning, in order to represent each technology's use of fossil fuels and mineral resources over time. Second, using TIAM-FR we determine a scenario following a baseline second shared socioeconomic pathway (SSP2) as described by [START_REF] O'neill | A new scenario framework for climate change research: the concept of shared socioeconomic pathways[END_REF]. Third, the prospective requirements of fossil fuels and metallic and nonmetallic mineral resources in EG technologies are estimated on the 2010-2100 horizon. Results at the regional scale are provided.

Fourth, we investigate how sensitive the total material use is to a change in the allocation of infrastructure LCIs between the construction and the decommissioning phases. Finally, we provide perspectives relating to the full integration of nonfuel raw materials into TIAM-FR.

Method

The TIAM-FR model TIAM is a bottom-up model in the TIMES family and was created in the framework of the ETSAP group at the IEA. TIAM-FR is the TIAM version that has been developed at MINES Paristech Center for Applied Mathematics to address energy and climate-related questions [START_REF] Dubreuil | Water modeling in an energy optimization framework -The water-scarce middle east context[END_REF][START_REF] Selosse | Achieving negative emissions with BECCS (bioenergy with carbon capture and storage) in the power sector: New insights from the TIAM-FR (TIMES Integrated Assessment Model France) model[END_REF][START_REF] Kang | Is GHG mitigation policy enough to develop bioenergy in Asia: A long-term analysis with TIAM-FR[END_REF]. It depicts the global energy system with a high level of disaggregation on energy sources, technologies, and end uses. A Reference Energy System (RES) connects commodities to current and future technologies in different energy supply sectors (upstream [fuel processing], heat, electricity) and end-use sectors (agriculture, commercial, industry, residential and transportation). The RES includes the collection, transformation, distribution, trade and end-use of main energy carriers including fossil energy carriers, bioenergy and uranium. These material commodities have defined resource supply curves (fossil resources, uranium) or potentials (solar, bioenergy, hydro and wind energy). They are converted into energy-equivalent quantities to be monitored and processed in a common accounting unit. Demand for steel, iron, pulp, and non-ferrous metals is also defined in the model, but only as an exogenous driver for the industrial sectors concerned. The individual supply of mineral resources is not modeled. In each of the 15 world regions represented (Africa, Australia and New Zealand, Canada, China, Central and South America, Eastern Europe, Former Soviet Union, India, Japan, Middle East, Mexico, Other Developing Asian countries, South Korea, United States of America, Western Europe), TIAM-FR establishes a balance between supply technologies and a set of forty-two demand drivers including gross domestic product (GDP) growth (supp. Figure 1), population growth (supp. Figure 2), demand for electricity, heat, transportation, lighting, hot water, etc. [START_REF] Loulou | ETSAP-TIAM: the TIMES integrated assessment model Part I: Model structure[END_REF]. The model computes the total net present value of the total annual cost of all activities, discounted at 5% from the selected reference year 2010, and then it optimizes the technological landscape until a least-cost solution is found. The optimization is performed on 11 periods the lengths of which are provided in supplementary Table 1.

Perfect competition and full foresight are assumed in the TIAM-FR model. These assumptions lead to economically ideal situations rather than forecasts, from where real phenomena and future pathways can be discussed and compared. A shared socioeconomic pathway 2 (O'Neill et al. 2014) baseline scenario is investigated where no climate constraints are set. This means the world follows a "business-as-usual" trajectory leading to a 2100 temperature rise ranging from between 2.5°C and 3°C compared to the preindustrial era.

The full model is used, but only the EG sector is analyzed due to a wide coverage of main EG technologies in the ecoinvent LCI database, and a consensus about the increasing need for electricity in the future, whether or not climate targets are achieved. The 103 current and future electricitygenerating processes that constitute the EG sector in TIAM-FR are grouped into ten categories: bioenergy, coal, geothermal, hydropower, oil, natural gas, nuclear energy, marine, solar photovoltaics, solar thermal, and wind. Individual technologies, classification, and service lifetimes can be found in the supplementary material (supp. Table 2). Most of them are documented with their current and future costs and performance on the ETSAP-TIAM community's website (Energy Technology Systems Analysis Program). Three electricity-related outputs are reported: the newly-installed electric capacities in GW, the electricity production in PW h, and the end-of-life capacities in GW. These quantities are annualized over each period by dividing by the period's length.

Life-cycle inventories

The ecoinvent 3.3 database is currently one of the largest process-based life-cycle inventory databases [START_REF] Wernet | The ecoinvent database version 3 (part I): overview and methodology[END_REF]. In ecoinvent, available activities, products and technologies are built from elementary processes. This structure makes it a consistent framework for life-cycle assessment (LCA) and product comparisons. Simapro 8.3 software is used to extract the LCIs.

LCIs rely on the definition of system boundaries, which describe what flows of energy, mass, and other substances enter and exit the system. ecoinvent's EG technologies datasets generally include a foreground list including all activities directly related to the construction, operation (comprising fuel extraction and transport), maintenance, and decommissioning of the power plants, and a background system comprising "upstream" and "downstream" activities such as material extraction, processing, transportation, and waste treatment. In this study, we aim at quantifying the life-cycle material requirements (or footprint) of EG technologies in order to determine what categories of materials will be most sensitive to electricity production. This is why we consider direct and indirect use of material flows and stocks but do not differentiate them. Thus, accumulation of stocks and outflows are not investigated. Several calculation methods are available depending on how foreground and background activities are linked to form a system (European Commission et al. 2010, 70-72). We opt for the default "allocation at the point of substitution" method and discuss the choice of an allocation method rather than a consequential method in the supplementary material (supp. section 4). Other databases and literature that could complement missing data or provide more "balanced" estimates were intentionally avoided for consistency, transparency and reproducibility purposes. The only external LCI we use is a 100 MW solar tower dataset that was provided by [START_REF] San Miguel | D6.3 Full Life Cycle Analysis of HYSOL Technology: Environmental and socio-economic assessment[END_REF]Corona et al. 2016) to represent concentrated solar power technologies for which no dataset was found in ecoinvent. Regular updates, verified data and a large number of users make the ecoinvent database suitable for comparing datasets, filling data gaps and spotting unrealistic or inconsistent assumptions [START_REF] Steubing | The ecoinvent database version 3 (part II): analyzing LCA results and comparison to version 2[END_REF].

Fossil fuels and mineral resources

Each LCI comprises a list of substances that are involved in the technology's activity within the selected boundaries. Substances may be energy, materials (solids, liquids, gases) or other LCA-pertinent indicators, such as land occupation and transformation. Categorizing raw materials can be difficult and requires establishing definitions. In this study, we use the Committee for Mineral Reserves International Reporting Standards' definition of mineral resources, which describe the latter as "a concentration or occurrence of solid material of economic interest in or on the Earth's crust in such form, grade or quality, and quantity that there are reasonable prospects for eventual economic extraction" [START_REF] Drielsma | Mineral resources in life cycle impact assessment-defining the path forward[END_REF]. Following this definition, mineral resources are divided into two categories: metallic and nonmetallic mineral resources. A specificity of ecoinvent is that some materials such as metals are reported not as refined products nor as ores, but as resource contents in crude ore, distinguishing all individual valuable resources. Thus, a metallic mineral resource refers to the metal content in metal ore. Fossil fuels (including peat) are considered in a separate category. We define raw materials as a more general category of natural resources of economic interest including mineral resources, fossil fuels, biomass, water, and other liquids and gases. We choose to restrict the scope of this study to fossil fuels and metallic and nonmetallic mineral resources that are found in the 'resource' compartment and 'in ground' sub-compartment as defined in ecoinvent (peat, located the 'biogenic' sub-compartment, is also considered a fossil fuel). A list of mineral resources and their categories can be found in the supplementary material (supp. Table 3).

Combining LCIs with TIAM-FR outputs

Our objective is to quantify and compare the amount of fossil fuels and metallic and nonmetallic mineral resources extracted and consumed by the EG technologies, not only in direct processes but also in upstream activities indirectly contributing to these processes. To do so, we attribute each TIAM-FR process with an LCI.

A variety of electricity-producing technologies are documented in ecoinvent, but not each of the 103 TIAM-FR technology processes has a corresponding LCI (e.g. concentrating solar power, marine and carbon capture and sequestration technologies). This is why the same LCI is attributed to all processes in a similar technical family (e.g. marine and wind offshore technologies are considered in the same technical family, and thus, are attributed with the same LCI dataset. However, wind offshore is grouped with wind onshore in the "wind" technological category, and not the "marine"). Furthermore, the availability of ecoinvent's LCI datasets for different world regions allows us to select a corresponding dataset for each of the 15 world regions in TIAM-FR. A nearby region, or the "Rest-of-World" dataset is used when a given region is not available in the database. As the technologies' service lifetimes are approximately equivalent in ecoinvent and TIAM-FR, no corrections are made to the LCIs or to the TIAM-FR processes.

LCIs of EG processes list all substances according to the functional unit kg/(kW h). This format is ideal for rating technologies that have short lifetimes and invariable life-cycle phases. However, it is not convenient for analyzing long-lifetime technologies with a high variability of material use throughout their construction, operation, and decommissioning phases. Furthermore, a single massper-electricity-output figure is often not adequate to represent the intensity of the technology's impacts (e.g. material use, CO 2 emissions, etc.) when they occur for a relatively short period of time (e.g. construction of a hydropower dam). In the TIAM-FR model, installed capacities may produce electricity intermittently or not at all, depending on the model constraints. For these technologies, a minimum of three phases -construction, operation, and decommissioning -seems reasonable to capture the differences in their material use over their life cycle. In this case, a convenient feature is the introduction of an additional functional unit: kg/W, for the construction and decommissioning phases where, instead of electricity, a capacity is produced (or removed). We propose to follow the evolution of life-cycle material use evolution in the three phases, multiplying each functional unit by their corresponding electric output in the scenario. A similar approach, where material consumption is assumed to be linearly driven by the production of a unit service -in our case a kilowatt or a kilowatthour -is common in the material flow-stock literature, e.g. [START_REF] Müller | Stock dynamics for forecasting material flows-Case study for housing in The Netherlands[END_REF] who used floor space per capita as a driver for housing material stock dynamics. Such an approach requires identifying the life phase in which each activity occurs, which is not always documented in ecoinvent.

For instance, construction, repairs and decommissioning activities are usually merged in an "infrastructure" category, while other activities are available in a "transformation" category. In order to overcome the life-phase separation issue, we make the following approximation: first, "transformation" activities are allocated to the operation phase; second, activities related to the power plant's infrastructure are allocated between the construction and decommissioning phases. We use a ratio  to set a percentage of the infrastructure LCI that is attributed to the decommissioning phase, while the complementary part, 1-, is attributed to the construction phase. Practically, this coarse approximation can be interpreted as if part of the construction of new power plant was instead performed during the decommissioning phase of an end-of-life power plant of the same type being replaced. This enables us to evaluate the sensitivity of the total material use to more or less materialintensive construction and decommissioning phases, while maintaining consistency and mass balance with the original infrastructure LCI. A more refined life-phase separation method was recently proposed by [START_REF] Arvesen | Deriving life cycle assessment coefficients for application in integrated assessment modelling[END_REF]) in order to combine LCIs with IAMs.

In this study, we use "system" LCIs, which means direct and upstream material requirements appear in an aggregated form. The available infrastructure LCIs (construction + decommissioning phases) are first expressed in kg of materials required to produce one power plant. These materials are then divided by their corresponding power capacities to yield kg/W figures. The transformation LCIs (operation phase) are shown by the unit g/(kW h). Figure 2 shows the schematic disaggregation method and coupling with the TIAM-FR outputs. The construction LCI is multiplied by the new capacities, the operation LCI is multiplied by the electricity production, and the decommissioning LCI is multiplied by the retiring (or end-of-life) capacities. The LCI datasets and their attribution to the TIAM-FR technology processes are detailed in the supplementary material (supp. Table 4).

Results

TIAM-FR scenario

Figure 3 shows the TIAM-FR scenario outputs. New and end-of-life power plant capacities are relative to the reference year (2010). Global electricity production grows from 21 PW h in 2010 to 45 PW h in 2050, and 61 PW h in 2100. In 2050, it is mostly generated by coal power plants with 20.0 PW h, followed by natural gas (8.5 PW h), hydropower (5.6 PW h), and nuclear (4.1 PW h). Solar and wind technologies become key electricity producers from 2040 (8.6 PW h and 5.1 PW h in 2100, respectively). Other sources of electricity production like oil, marine, bioenergy, and geothermal energies together contribute less than 10% of total generation. New power capacities of gas and coal increase until 2050 when they fluctuate in the range of 65-120 GW/yr. Wind capacity increases rapidly but is caught up by solar photovoltaics (PV) in 2040, which reach a global average installation rate of 90 GW/yr in the second half of century. Nuclear, marine and hydro power plants develop slightly and reach annual new capacities of 5 GW/yr to 25 GW/yr. Bioenergy electricity production (2.7 PW h in 2050) is supported by new methane abatement technologies (e.g. biogas). However, these technologies have no defined capacities and, thus, their infrastructure cannot be considered (only material use during operations is considered). End-of-life capacities consist mostly of coal, natural gas, and nuclear power plants before 2040. An increasing number of wind farms and solar power plants are decommissioned in the second half of the century. Periods with increased end-of-life capacities are generally compensated by new capacities to support the growing electricity production.

Life-cycle material use of electricity-generating technologies

Figure 4 shows the sum of the masses of all fossil fuels and mineral resources used in the infrastructure and operation phases of each power-generating technology, averaged over the time periods and regions. There is a clear distinction between fossil and bioenergy power plants versus alternative technologies. Fossil fuel and bioenergy technologies have lower material requirements for their infrastructure -below 0.7 kg/W-but high material consumption during operation -120 g/ (kW h) (bioenergy) to 664 g/(kW h) (coal)-due to the large quantities of fuels consumed. Of the fossil energy technologies, natural gas demonstrates the lowest material needs for both infrastructure and operation due to its higher efficiency and lighter infrastructure. Coal and oil power plants have additional mechanical and chemical processing steps for fuel and wastes management. In contrast with fossil fuel technologies, alternative technologies have higher material needs for their infrastructurefrom 1 kg/W (wind, solar, and marine technologies) to 132 kg/W (hydropower) -but lower needs when operating, at about 510 5 g/(kW h) (solar PV) to 0.1 g/(kW h) (wind).

Nuclear plants have intermediate material requirements for both the infrastructure and operation phases at 0.9 kg/W and 7 g/(kW h), respectively. Solar thermal (concentrating solar power) has higher infrastructure requirements than nuclear but use fewer materials during operations. The large variations observed among the technologies reflect the very different ways materials are involved in the construction, operation and decommissioning life phases across the technologies.

Life-cycle material use of EG technologies is generally very sensitive to location, size, and technology vintage [START_REF] Gibon | Life cycle assessment demonstrates environmental co-benefits and trade-offs of low-carbon electricity supply options[END_REF]. This is especially the case for hydropower dams, for which large amounts of materials can be displaced and huge quantities of cement and steel can be required depending on the construction site. This is confirmed by IPCC AR5 (fig. 7.6, Bruckner et al. 2014), where hydropower has the greatest uncertainty and the largest GHG emission range . Figure 5 shows the life-cycle use of aluminum, copper and iron by several EG technologies, according to this study (ecoinvent 3.3), and studies by [START_REF] Hertwich | Integrated life-cycle assessment of electricity-supply scenarios confirms global environmental benefit of low-carbon technologies[END_REF] and (de [START_REF] Koning | Metal supply constraints for a low-carbon economy?[END_REF].

technologies and materials, except for the use of aluminum in the solar tower dataset (112%).

Observed differences are due to the consideration of different technologies, especially for coal and natural gas. Another reason is that [START_REF] Hertwich | Integrated life-cycle assessment of electricity-supply scenarios confirms global environmental benefit of low-carbon technologies[END_REF]) used custom LCIs with a hybrid approach combining process-based LCIs and a multiregional environmentally extended input-output database to cover economic activities that are not addressed in process-based LCIs [START_REF] Gibon | A Methodology for Integrated, Multiregional Life Cycle Assessment Scenarios under Large-Scale Technological Change[END_REF]. Therefore, some activities may be considered in their analysis that are truncated in ecoinvent 3.3. Nevertheless, the relative uses of the three materials are consistent. A relevant comparison between the technologies and the time periods is thus possible, while the individual value of material use should be interpreted inseparably from the system boundaries (Majeau-Bettez et al. 2011).

We deliberately choose not to adjust nor calibrate the ecoinvent datasets to maintain a high level of reproducibility. As the ecoinvent database is often used for technology comparisons, our approach enables one to identify potential gaps and investigate the convergence between prospective LCI, MFA and integrated assessment modeling.

Material use dynamics

requirements for each technology, by material type. The results are calculated at  = 0.1, meaning 10% of each substance in the infrastructure LCI is attributed to the decommissioning phase, and 90% to the construction phase. This value is inspired by [START_REF] Hertwich | Integrated life-cycle assessment of electricity-supply scenarios confirms global environmental benefit of low-carbon technologies[END_REF], which assume energy requirements for decommissioning power plants are 10% of those for their construction.

Figure 6 shows the life-cycle use of fossil fuels and metallic materials and nonmetallic mineral resources generated by the developing EG technologies (see full data in supp. Table 5). In the construction phase (6a), fossil fuels are mostly used by hydropower in the first half of century (76 Tg/yr in 2050) and by solar PV in the second half (94 Tg/yr in 2100). At the end of horizon, the development of geothermal facilities requires a substantial amount of fossil fuels (46 Tg/yr in 2100).

Fossil fuel consumption in the construction phase of other technologies (marine, natural gas and nuclear) is 33 Tg/yr in 2050 and 43 Tg/yr in 2100. Looking at material use in the decommissioning phase (6c), about 0.8 Tg yr -1 of fossil fuels are used in the coal and nuclear capacity decommissioning phases. After 2030, decommissioning of solar power plants is responsible for the largest share of fossil fuel use in this phase. However, use of fossil fuels in both the construction and decommissioning phases is negligible compared to the life-cycle amount of fossils burned in coal and natural gas power plants during operations (6b) (14,380 Tg for coal and 895 Tg for natural gas in 2100). Total annual consumption of fossil fuels (6d) due to EG technologies plateaus in 2090 at 15,648 Tg/yr.

Comparing our results with historical statistics, 2.1 Gtoe of coal and 1.0 Gtoe of natural gas were consumed directly by electricity and CHP plants in 2010 according to (International Energy Agency 2017). This is 3627 Tg of coal and 856 Tg of natural gas, assuming IEA's average net calorific values of 24.05 kg/MJ for coal and 50.03 kg/MJ for natural gas (International Energy Agency 2005). In contrast, our model gives a life-cycle consumption of 5314 Tg of coal and of 1349 Tg of natural gas in 2010, which is 47% and 58% higher than the direct consumption, respectively. However, IEA statistics do not include the life-cycle effects associated with each sector, which partly explains the difference observed. Another reason is the use of a limited number of LCIs for coal and natural gas technologies (two LCI datasets each) in our analysis.

The use of metallic resources in the construction phase (6a) follows a similar evolution to that of fossil fuels, except for a relatively larger share of metals in coal and wind generation plants, and a lower share in solar power plants (compared to their share of fossil fuels). Metal use due to decommissioning activities (6c) is dominated by coal power plants. Solar and wind end-of-life capacities generate substantial use of metallic resources in the second half of century. The consumption of metals is mostly driven by fossil and bioenergy plants during the operation phase (6b). Total use of metallic mineral resources (6d) is 23 Tg/yr in 2010 which is in agreement with [START_REF] Hertwich | Integrated life-cycle assessment of electricity-supply scenarios confirms global environmental benefit of low-carbon technologies[END_REF] which estimated the sum of life-cycle consumption of iron, aluminum and copper by EG technologies at 23 Tg/yr in 2010, in an IEA baseline scenario. However, use of metallic mineral resources reaches 80-90 Tg/yr by 2050 while it stays constant in their study due to decreasing material intensities.

Most of the nonmetallic mineral resources used are attributable to the new hydropower infrastructure (6a) for which outstanding material use is assumed in the ecoinvent LCI. The total use of nonmetallic materials (6d) tends to decrease after a surge in 2040 resulting from the construction of large hydropower plants. , consumption of nonmetallic resources is also substantial due to coal, gas, and bioenergy power generation, which displace large amounts of raw materials during mining activities and for the production of biofuels. In the second half of century, decommissioning activities (6c) in solar, wind, coal, marine and nuclear plants generate most of the use of nonmetallic mineral resources. The total material use generated in the decommissioning phase (6c) ranges from 1% to 3% of the total material use observed in the construction phase (6a) (at  = 0.1).

The totals plot (6d) shows that the amounts of mineral resources used increase rapidly. The total amounts of mineral resources required to construct power plants (6a) overtakes the amount of mineral resources needed for operations (6b) as early as 2020, due to the development of hydropower.

However, the amount of fossil fuels and biofuels burnt during the operation phase (6b) by far exceeds the total consumption of mineral resources.

Looking at the total material use in 2100, China is the highest consumer with about half (8,385 Tg/yr) of global use of fossil and mineral resources (17,659 Tg/yr ), followed by India (2,879 Tg/yr ) and the USA (2,221 Tg/yr ), while Western Europe (764 Gg/yr ) and other developing Asian countries (632 Gg/yr ) come in 4th and 5th positions. These results feature in more detail in the supplementary material (section 7).

The 2010 global extraction of fossil fuels, metal ores, and nonmetallic minerals were estimated by [START_REF] Schandl | Global Material Flows and Resource Productivity: Forty Years of Evidence: Global Material Flows and Resource Productivity[END_REF]) at 13,300 Tg, 7,200 Tg, 30,700 Tg, respectively. Comparing these estimates with our results, we find electricity and CHP plants were responsible for the life-cycle use of 51.3% of global fossil fuels (compared to 36% [4,732 Tg] according to IEA's direct consumption data and previously assumed calorific values), 0.3% of global metal ores, and 1.3% of global nonmetallic minerals in 2010. However, metallic mineral resources only include the metal content in ore, while the nonmetal part is considered as nonmetallic mineral resource. Summing all mineral resources -thus overcoming this counting issue-we find 1.1% of global mineral resources extracted were used in the life-cycle of EG technologies in 2010.

Sensitivity to a material use impact shift towards the decommissioning phase

When  varies, the material consumption dynamics are modified. Figure 7 shows the total annual use of fossil fuels and metallic and nonmetallic mineral resources for three values of :  = 0,  = 0.5, and  = 1. As expected, fossil fuel consumption is not significantly affected by  as most of the fossil fuels are used during the operation phase. When  is higher, material use is lower as there are fewer end-oflife capacities than new ones. The raw mineral resources that were assumed in the construction phase of hydropower are partly shifted to the decommissioning phase (which occurs after 2100 due to an assumed service lifetime of 100 years). Compared to the case where  = 0, the material requirements at  = 0.5 in 2050 are 0.5% for fossil fuels, 20.8% for metallic mineral resources, and 32.4% for nonmetallic mineral resources.

Discussion

The TIAM-FR model -as well as most IAMs -does not consider the mineral resources needed for products and infrastructures in the energy system. This results in prospective scenarios that assume unlimited mineral resource supply. Not only does it prevent IAMs from assessing the impacts of mineral resource cycles and circular economy policies, but it can also generate inconsistencies when energy and climate policies are in fact associated with unrealistic material availability.

The purpose of this prospective exercise is not to predict future material consumption to 2100 but to demonstrate the possibility to combine the bottom-up structure of TIAM-FR with life-cycle inventories to extract valuable information on the relations between energy scenarios and material resource use. Linking the TIAM-FR model with LCIs, we show that the SSP2 scenario investigated has strong implications for fossil fuel and mineral resource use, with a high increase in the use of metallic and nonmetallic mineral resources due to the development of low-carbon EG technologies. A specific feature of our approach is that all fossil fuels and mineral resources considered in ecoinvent are analyzed, providing an overview of material consumption in the EG sector. The very different dynamics of material use in the construction, operation and decommissioning phases could not be represented using aggregated LCIs, especially for large and long-lived facilities. Using disaggregated LCIs, we showed that periods of investment in new capacities are associated with a strong increase in the resource requirements. This may in turn generate increasing GHG emissions. Thus, low-carbon solutions must be found at every stages of the life-cycle of renewable energy technologies, from material extraction to recycling and waste disposal. In the same time, a zero carbon emission world is likely to have a much larger mineral resource footprint if economic growth is maintained. Therefore, it is of crucial importance to investigate the material limits to the global energy system's development, and establish circular economy scenarios.

From this perspective, soft-linking IAMs with LCIs can be a convenient and useful solution to energy modelers and decision-makers in performing a first screening of energy and environmental impacts of material resources on various socioeconomic pathways.

Based on the current results, it is not possible to say that EG technologies will face material supply shortages. To answer this question, the present method could be adapted to cover other resourceintensive sectors of the energy system (e.g. tranportation, electricity consumption, or building) using adequate functional units. This would require separating direct from indirect material use for each process LCI to avoid double-counting due to overlapping supply chains. Another important task would be to analyze the material use dynamics of individual resources regarding their specific availability and environmental impacts (e.g. GHG emissions, land use, abiotic resource depletion, etc.).

Eventually, such an approach could lead to a full description of material flows in the energy system, and help determine resource constraints not currently considered in IAMs.

Our approach does not circumvent errors in ecoinvent and TIAM-FR, nor does it palliate missing or outdated data. On the contrary, it makes it easier to identify them, and if relevant, to correct them, or submit new datasets. Thus, we would like to invite other integrated assessment modelers to look at the nonfuel material issue and develop scenarios that are compatible with industrial ecology principles.

Further work and recommendations

Life-cycle inventory is one of the most effective tools to compare different technologies or products based on their life-cycle material use. However, datasets often come from different sources, use different data, at different times, and assume different system boundaries and linking models. Expert judgment is thus generally required to combine, harmonize, correct or update existing LCIs, depending on the question addressed. This leads to seemingly more accurate estimates, but also to a profusion of LCI datasets that may be reused in situations where initial data sources are unclear and where their applicability is questionable [START_REF] Masanet | Life-Cycle Assessment of Electric Power Systems[END_REF]. In this study, we took an alternative route, relying on the ecoinvent database only. This choice also comes with its own issues and limitations. Even though ecoinvent is one of the most complete LCI databases, only a limited number of technologies are represented (e.g. no CCS technologies, no CSP plants). Each has its own data gaps and assumptions, which might not be fully consistent with others, or with their associated TIAM-FR technologies. The LCI system boundaries may generate discrepancies in total material requirements, by including or excluding economic sectors and activities that might -or might not-be included in other LCIs. For example, the hydropower dam infrastructure LCI include large earth movements that end up as wastes, while these wastes could be used for other functions (building, roads, etc.). In addition, the dataset is calculated as an average technology (of Swiss dams) for which no power capacity is explicitly provided.

It was shown previously that process-based LCIs cannot cover all economic sectors and may underestimate material use and environmental impacts by 20% -60% following truncation [START_REF] Majeau-Bettez | Evaluation of Process-and Input-Output-based Life Cycle Inventory Data with Regard to Truncation and Aggregation Issues[END_REF]. From this perspective, life-cycle inventory faces the same "bottom-up vs. topdown" issue that took place in the integrated assessment community and led to the development of hybrid IAMs [START_REF] Hourcade | Hybrid Modeling: New Answers to Old Challenges Introduction to the Special Issue of The Energy Journal[END_REF]. Hybrid LCI methods represent the "natural" evolution towards analyzing more and more complex systems with a consequential insight. [START_REF] Majeau-Bettez | Evaluation of Process-and Input-Output-based Life Cycle Inventory Data with Regard to Truncation and Aggregation Issues[END_REF] writes: "the development of generic life cycle data sets should be seen as something more than the mere compilation of individual inventories, but rather as an opportunity to model the metabolism of society". In this direction, segmenting the life cycle of products into multiple phases could provide clearer information for decision-making, and a higher flexibility for modeling time-dependent consequential LCIs.

Ideally, information loss between LCI providers and interpreters should be minimized, and the life phase during which each activity occurs should be preserved, especially for large technologies or products that have heterogeneous material consumption intensities and extended service lifetimes.

In this analysis, the LCIs are static and identical by scenario. Using invariable allocation LCIs may create deviations in the material use when the surrounding environment has changed and no longer corresponds to the initial system. The further the horizon, the higher the uncertainty and bias generated by use of static LCIs. To partially address this issue, LCIs can be adjusted according to the evolution of the background system in order to reflect, for example, a more or less carbon-intensive power generation sector [START_REF] Hertwich | Integrated life-cycle assessment of electricity-supply scenarios confirms global environmental benefit of low-carbon technologies[END_REF][START_REF] Pehl | Understanding future emissions from low-carbon power systems by integration of life-cycle assessment and integrated energy modelling[END_REF]). However, several key aspects of existing prospective and consequential approaches are left unaddressed, such as the increasing material use of decarboning resource extraction treatment processes, as well as the non-linearities involved in consequential LCIs of technologies, the service lifetime of which may range from 15 years to more than 100 years. Hybrid or consequential LCIs that are tied to a specific scenario may not be more applicable to other scenarios than an allocation approach. This is why consequential LCI needs improvement to make more room for scenario customization [START_REF] Masanet | Life-Cycle Assessment of Electric Power Systems[END_REF][START_REF] Pauliuk | Lifting Industrial Ecology Modeling to a New Level of Quality and Transparency: A Call for More Transparent Publications and a Collaborative Open Source Software Framework: Open Source Software for Industrial Ecology[END_REF], while keeping some representative invariants allowing relevant comparisons.

While LCIs may evolve depending on the scenario, the scenario itself may also be affected by LCIs, especially if technological deployment is constrained by material supply stresses. Raw material reserves, life-cycle impacts, and criticality could be used to adjust the technologies' costs in the TIAM-FR model. This could be done by integrating raw materials and LCA metrics directly into the TIAM-FR modeling structure. In this case, LCA, MFA and prospective scenario analysis could be achieved simultaneously to inform energy, resource, and environment decision-making regarding life-cycle aspects [START_REF] Masanet | Life-Cycle Assessment of Electric Power Systems[END_REF]. TIAM-FR is already capable of dealing with material flows and stocks in technology processes using specific parameters. However, achieving such an integration into the TIAM-FR model would necessitate considerable amounts of dynamic MFA data defining how material use might evolve over time, driven by efficiency changes, material substitution, economic reserves, prices, and many other parameters.

Evolution towards a circular economy is now recognized as a way to contribute to achieving the United Nations Sustainable Development Goals (SDGs) (European Commission 2015). As secondary material resources (i.e. recycled or reused) are expected to play an increasing role in the future [START_REF] Ekins | Resource Efficiency: Potential and Economic Implications[END_REF], investigating their implications in the TIAM-FR model could help determine more effective resource efficiency policies, especially in the context of climate change where environmental constraints are likely to exacerbate the stress on raw materials. Thus, an interesting extension to this work would be to analyze the build-up of in-use stocks and outflows of materials by distinguishing flows and stocks in the LCI datasets. Thermoecological approaches based on exergy and entropy [START_REF] Stanek | Fundamentals of Exergy Analysis[END_REF][START_REF] Valero | Exergy as a resource efficiency indicator for industries[END_REF] appear to us promising ways to describe material cycles more physically.

Uncertainties were not addressed in this paper due to the fact that many sources of uncertainty are not identified (in LCIs) or currently impossible to quantify (in the TIAM-FR model). Further work is needed to understand how uncertainty can be better dealt with in more and more complex and integrated decision-making tools.

Transparency and reproducibility were motivations for using the ecoinvent database. At IEA, similar motivations led to the decentralization of the TIAM model, which is now used and developed by research teams in nearly seventy countries within the ETSAP community. However, it is clear that IAMs are still perceived as black boxes not only by the public but also by other research communities. This is why we subscribe to the ongoing discussions on open science and the need for more transparency in energy research [START_REF] Pfenninger | Opening the black box of energy modelling: Strategies and lessons learned[END_REF].

Conclusion

Combining the TIAM-FR global energy system model with disaggregated life-cycle inventories enabled us to compute the prospective life-cycle requirements of fossil fuels and mineral resources due to electricity-generating technologies to 2100. Fossil fuels and metallic and nonmetallic mineral resource requirements were computed in the construction, operation and decommissioning phases of the technologies, providing the first prospective assessment of mineral resource requirements by life phase and region. The results show that the development of hydro, solar and wind power generate increasing global demand for metallic and nonmetallic minerals within the next few decades in the scenario investigated, mainly due to their infrastructure. End-of-life capacities rapidly increase after 2040, which means more fossil fuels and mineral resources could be needed for decommissioning activities. The presented method can be adapted to investigate other prospective scenarios, and possibly other energy sectors, such as transportation, building construction and information technology. Total integration of raw materials into the TIAM-FR model still needs to be developed, as the current method does not allow endogenous interactions with the scenarios. This will require identifying those material-related parameters that may have a significant influence on the scenario (supply availability, material criticality, etc.). Eventually, LCA, MFA, and IAM may converge into a single tool to devise LCA-oriented energy and resource policies.


Tecla Bonci

Valentina Camomilla

Raphaël Dumas

Laurence Cheze

Aurelio Cappozzo

A modal approach for the soft

INTRODUCTION

Optimal bone pose and joint kinematics estimators, that use data provided by stereophotogrammetry, aim at compensating for soft tissue artefacts (STA). We suggest that they embed a mathematical representation of the STA that adequately approximates it and uses a sustainable number of parameters. To this purpose we propose to represent STAs with a modal approach using different STA definitions: individual marker displacements (MrkD), marker-cluster geometrical transformations (GT), or skin envelope shape variations (SV) [1]. Modes can be ranked and a number of them selected leading to different degrees of approximation of the STA reconstruction and, therefore, compensation. This study investigated how, for each STA definition, the selected number of modes, the resulting STA approximation and accuracy of joint kinematics estimates relate to find the most appropriate STA mathematical representation.

MATERIALS and METHODS

Experimental data provided by Reinschmidt et al., (Journal of Biomechanics, 1997:729-732) were used. STA vectors and reference knee kinematics were determined using intracortical pin data during running trials. For each segment i, a STA field V i (k) was generated. For each STA definition, V i (k) was split into additive modes (amplitude a i l (k) and direction Φ i l ) that were ranked. As a paradigmatic case, for each definition, a number of modes (r i ) was selected that represented the 95% of the STA deformation energy and the residual STA (5%) fields computed Ṽ i (k). Knee kinematics was calculated, through an SVD approach, using skin marker positions embedding V i (k) (raw data: RD) or Ṽ i (k) (STA virtually compensated data) associated with each STA definition (MrkD, GT, SV). For each of the four kinematics estimates thus obtained, the rms difference with respect to the reference kinematics was calculated (rmse).

RESULTS

DISCUSSION

The selection of a STA definition and its modal representation allows for a trade off between virtual STA compensation effectiveness and number of modes, that is number of parameters involved, and, therefore, mathematical tractability of the bone-pose estimator. The GT definition grants for the best knee kinematics estimate accuracy, but uses 6-10 modes; while the SV definition allows for a slightly worst accuracy, but uses only 2-3 modes.

KEY REFERENCE

[1] Dumas, R., Camomilla, V., Bonci, T., Cheze, L., Cappozzo, A., 2014. Generalized mathematical representation of the soft tissue artefact. Journal of Biomechanics 47,476-481


E De Prunelé

The hard sphere quantum propagator: exact results via partial wave analysis

An exact expression for the partial wave hard sphere quantum propagator in threedimensional space is given in terms of the Fourier transform of an expression involving spherical Bessel functions. For the s wave case, the Fourier transform is calculated analytically. This Fourier transform is also evaluated for the general case by contour integration in terms of a Laplace transform and residues contributions. The accuracy of previous approximations for a hard sphere quantum propagator can thus be evaluated. In particular, the Van Vleck-Gutzwiller classical approximation is found to be very accurate provided that the classical action is greater than only a few units of reduced Planck constant.

Introduction

Propagators are important conceptual tools in quantum mechanics. For a given physical system, propagators are matrix elements of the time evolution operator U(t f , t i ) which can be used to obtain the state |ψ(t f ) of the system at time t f if its state at time t i is |ψ(t i ) . In the Schrödinger picture : |ψ(t f ) = U(t f , t i )|ψ(t i ) . For example, they are studied and discussed in the non-relativistic context in [START_REF] Feynman | Quantum Mechanics and Path Integrals[END_REF] and in both non-relativistic and relativistic contexts in [START_REF] Bjorken | Relativistische Quantenmechanik Band[END_REF]. Propagators are particularly important for Feynman formulation of quantum mechanics in terms of path integrals.

In the non-relativistic context to be considered here, analytic exact expressions for propagators are known for only a few simple systems, for example, free particles or harmonic oscillators [START_REF] Feynman | Quantum Mechanics and Path Integrals[END_REF][START_REF] Beauregard | [END_REF]. An exact analytic expression for the propagator of a particle interacting only with a hard sphere is not known. This has motivated the search for approximate expressions. The so-called image approximation propagators [START_REF] Uhlenbeck | [END_REF][5][6] have been the first ones. They have now been superseded by a more recent attempt of Cao and Berne [7], to be called hereafter the CB propagator or the CB approximation. The hard sphere propagator, for real or imaginary time, is a basic ingredient for the study of thermodynamic properties of hard sphere systems (fluid or solid) and Bose-Einstein condensation of gases [8][9][10][11][12][13][14][15][16]. The CB approximation is based on a purely quantal approach to scattering, and has largely been used in these applications. But there is also another approximation for the hard sphere quantum propagator [17][18][START_REF] Gutzwiller | Chaos in Classical and Quantum Mechanics[END_REF], involving only classical quantities as inputs. We shall call it the Van Vleck-Gutzwiller propagator, to be called hereafter the VG propagator or the VG approximation, which has been used for example in [START_REF] Goussev | [END_REF].

In the present paper, an exact expression for the partial wave hard sphere quantum propagator is given in terms of the Fourier transform of an expression involving spherical Bessel functions. This Fourier transform is also evaluated by contour integration in terms of a Laplace transform and residues contribution (see section 4.3). The accuracy of previous attempts at approximations for the hard sphere quantum propagator can thus be evaluated by comparison with exact results obtained by numerical integration of the exact expressions based on Fourier or Laplace transforms. It is also shown that for the s wave hard sphere quantum propagator, the Fourier transform can be calculated analytically, leading to a simple expression (see equation (15)).

Our starting point is the well-known expression (see, e.g., [START_REF] Bjorken | Relativistische Quantenmechanik Band[END_REF][START_REF] Cohen-Tannoudji | Processus D'interaction Entre Photons et Atomes[END_REF]) of the retarded time evolution operator as a Fourier transform of the resolvent G(z) = (z -H ) -1 for a system with time-independent Hamiltonian:

θ(t f -t i )U (t f -t i ) = - 1 2π i +∞ -∞ dE exp -i E(t f -t i ) h G(E + ). ( 1
)
The Heaviside function θ is zero when its argument is negative, and unity when its argument is positive. The z-plane of G(z) is cut along the real positive axis, and the notation G(E + ) means lim ε->0 G(E + iε) with ε > 0. Otherwise stated, on the integration interval [0, +∞], the cut has to be approached from above, i.e. with a vanishing positive imaginary part. Equation ( 1) can be verified by considering the inverse Fourier transform [START_REF] Cohen-Tannoudji | Processus D'interaction Entre Photons et Atomes[END_REF].

From now on, only the positive time interval t = t ft i will be considered, so that the Heaviside function will not appear in all the following.

An exact expression in terms of Fourier transform

What is the expression of the resolvent for a particle of mass M, in the three-dimensional space, interacting only with a hard sphere of radius a ? To answer this question, let us first consider the Hamiltonian H (λ):

H (λ) = H 0 + λV H 0 = p 2 2M V = a 3 +∞ =0 m=-
|a, , m a, , m|.

The vector |r, , m is a eigenvector of the squared orbital angular momentum with eigenvalue ( +1)h 2 , an eigenvector of the component z of the orbital angular momentum with eigenvalue mh, and a generalized eigenvector of the radial position operator with generalized eigenvalue r. All operators are relative to a given frame O, x, y, z. The normalization is [START_REF] De Prunelé | [END_REF] r , , m |r, , m = δ(rr) r 2 δ δ mm .

It has been shown (see equation (13) of [START_REF] De Prunelé | [END_REF]) that the stationary scattering partial wave functions for H (λ) are r, , m|p, , m, + = i 2 π j (pr) -j (pa)

1 λ2Ma 3 + ph + (pa)j (pa) ph + (pr)j (pa) (2)
with r a. It is easy to verify that for the limiting case λ → +∞, these stationary scattering partial wavefunctions are eigenfunctions of H 0 with vanishing boundary condition on the sphere of radius a (see also [START_REF] Morse | Methods of Theoretical Physics[END_REF]). It can also be noted that an explicit expression for the partial wave-scattering amplitude for H (λ) can be obtained (see equation (12) of [START_REF] De Prunelé | [END_REF]). For the limiting case λ → +∞, these partial wave-scattering amplitudes then lead exactly to the phase shift for hard sphere scattering (compare with equation (50) of section X -13 of [START_REF] Messiah | [END_REF]). Now we are interested in the partial-wave hard sphere propagators:

r 2 , , m|U(t)|r 1 , , m = - 1 2πi ∞ -∞ dE exp -i Et h r 2 , , m|G(E + )|r 1 , , m (3)
and we first consider a partial-wave resolvent for the Hamiltonian H (λ). From the algebraic identity

G λ (z) = G 0 (z) + G 0 (z)λV G λ (z) one obtains r 2 , , m| G λ (z) |r 1 , , m = r 2 , , m| G 0 (z) |r 1 , , m + λa 3 r 2 , , m| G 0 (z) |a, , m a, , m| G λ (z) |r 1 , , m a, , m|G λ (z)|r 1 , , m = a, , m|G 0 (z)|r 1 , , m + λa 3 a, , m|G 0 (z)|a, , m a, , m|G λ (z)|r 1 , , m a, , m|G λ (z)|r 1 , , m = a, , m|G 0 (z)|r 1 , , m 1 -λa 3 a, , m|G 0 (z)
This equation could also be obtained more directly by noting that the right member of equation ( 4) satisfies the correct differential equations with respect to the variable r 1 or r 2 with vanishing boundary conditions at r 1 = a and at r 2 = a, and correct asymptotic behavior (see equation (5)). An essentially analogous expression for the two-dimensional resolvent for a hard disk is given in [START_REF] Brack | Semiclassical Physics[END_REF] with a reference to an unpublished thesis [START_REF] Snaith | [END_REF]. In [START_REF] Hill | [END_REF], the expression for the partial wave hard sphere resolvents (equations ( 4) and ( 5)) is given for negative values of the energy (see equation ( 21) of [START_REF] Hill | [END_REF]). In [START_REF] Hill | [END_REF], the starting point is different from that of the present paper in two respects. Firstly, it is the resolvent that is expressed as an integral transform of the evolution operator (the first equation of [START_REF] Hill | [END_REF]), whereas it is the evolution operator that is expressed as an integral transform of the resolvent in our equation [START_REF] Feynman | Quantum Mechanics and Path Integrals[END_REF]. Secondly, imaginary time for the evolution operator is considered in [START_REF] Hill | [END_REF] in order to consider the Boltzmann operator exp(-βH ). An asymptotic series for the resolvent is then inverse-Laplace transformed term-by-term in [START_REF] Hill | [END_REF] to obtain high-temperature successive quantum corrections for hard spheres second virial coefficients. The interest in the present method is that it gives the partial wave resolvent not only for λ → ∞ (hard sphere), but also for arbitrary λ values. The well-known expression for the free resolvent (see, e.g., [START_REF] Taylor | Scattering Theory[END_REF]) is:

r 2 , , m |G 0 (E + )|r 1 , , m = -δ δ m m 2Mph + (pr > )j (pr < ) (5)
with

p = (2ME) 1/2 . ( 6
)
r > is the largest value and r < is the smallest value of r 2 , r 1 . The spherical Bessel functions are defined as in [START_REF] Messiah | [END_REF]. Specifically, the spherical Bessel functions j and h + are related to the usual Bessel function J +1/2 regular at origin (see e.g [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]) by the equation

j (z) ≡ π 2z 1/2 J +1/2 (z) h ± (z) ≡ (-1) π 2z 1/2 J --1/2 (z) ± ij (z).
It is stressed that (2ME) 1/2 is defined on the z plane with a cut along the positive real axis. This means that on the first Riemann sheet, the argument of z to be denoted by arg 1(z) satisfies 0 arg 1(z) < 2π and z 1/2 = |z| 1/2 exp(i arg 1(z)/2). Thus z 1/2 should not be confused with √ z which is traditionally defined on the z plane cut along the negative real axis, i.e with an argument arg(z) which satisfies -π < arg(z) π, so that √ z = |z| 1/2 exp(i arg(z)/2). On the second Riemann sheet of the complex E plane, the argument of z to be denoted by arg 2(z) satisfies -2π arg 2(z) < 0 and z 1/2 = |z| 1/2 exp(i arg 2(z)/2).

One finally obtains

r 2 , , m|G(E + )|r 1 , , m = r 2 , , m|G 0 (E + )|r 1 , , m + r 2 , , m|G n (E + )|r 1 , , m (7)
r 2 , , m|G 0 (E + )|r 1 , , m = -2Mph + (pr > )j (pr < ) (8) r 2 , , m|G n (E + )|r 1 , , m = 2Mph + (pr 2 )h + (pr 1 )j (pa) h + (pa) , ( 9
)
where G n (E + ) represents the difference between the full resolvent of the particle in the presence of a hard sphere and the free particle resolvent. From now on, G n will be called the non-free part of the hard sphere resolvent. To avoid any confusion on the vocabulary, we write now several equations as a single index equation:

r 2 , , m|U k (t)|r 1 , , m = - 1 2πi ∞ -∞ dE exp -i Et h r 2 , , m|G k (E + )|r 1 , , m . ( 10
)
If the index k is absent, the left-hand side is the partial wave hard sphere propagator, and r 2 , , m|G(E + )|r 1 , , m on the right-hand side is the partial wave hard sphere resolvent. If the index k is 0, the left-hand side is the free partial wave propagator, and r 2 , , m|G 0 (E + )|r 1 , , m on the right-hand side is the free partial wave resolvent. If the index k is n, the left-hand side is the non-free part of the partial wave propagator, and r 2 , , m|G n (E + )|r 1 , , m on the right-hand side is the non-free part of the partial wave resolvent. Finally, the transition from the partial wave propagator to the propagator is achieved through a series involving Legendre polynomial P :

r 2 |U k (t)|r 1 = ∞ =0 2 + 1 4π P (cos(α)) r 2 , , m|U k (t)|r 1 , , m (11)
with α being the angle between r 2 and r 1 . If the index k is absent, the left-hand side will be called the hard sphere propagator, if the index k is 0, it will be called the free propagator and if the index k is n, it will be called the non-free part of the hard sphere propagator. Alternatively, the partial wave propagator is obtained from the propagator through the relation

r 2 , , m|U k (t)|r 1 , , m = d 2 d 1 Y m ( 2 ) r 2 |U k (t)|r 1 Y m ( 1 ) (12)
with Y m a spherical harmonic, and a solid angle.

A well-known exact result (see, e.g., [START_REF] Feynman | Quantum Mechanics and Path Integrals[END_REF][START_REF] Gutzwiller | Chaos in Classical and Quantum Mechanics[END_REF][START_REF] Durand | Mécanique Quantique Tome[END_REF]) for the free propagator is

r 2 |U 0 (t)|r 1 = M 2πhit 3 exp i M[r 2 -r 1 ] 2 2th . ( 13
)
It is worth noting that the modulus of the free propagator is independent of positions and depends only on the length √ ht/M. From now on, we shall work with units where h = 1, and we do no longer write the reduced Plank constant in the formulae. From this result, and the relations

r 1 |r, , m = δ(r 1 -r) r 2 Y m (r 1 ) exp(iβ[r 2 -r 1 ] 2 ) = exp iβ r 2 2 + r 2 1 4π ∞ =0 m =- (-i) j (2βr 2 r 1 )Y m (r 2 )Y m (r 1 )
one obtains the partial wave free propagator:

r 2 , , m|U 0 (t)|r 1 , , m = 4π M 2πit 3 exp i M 2t r 2 2 + r 2 1 (-i) j M t r 2 r 1 . ( 14
)

An exact analytical result for the s wave hard sphere propagator

For the s wave, = 0, the spherical Bessel functions are very simple

j 0 (z) = sin(z) z h + 0 (z) = exp(iz) z
and equations ( 10), ( 5), (9) give for r 2 r 1

π ir 2 r 1 M r 2 , 0, 0|U 0 (t)|r 1 , 0, 0 = ∞ -∞ dE exp(-iEt) exp(ipr 2 ) sin(pr 1 ) p - πir 2 r 1 M r 2 , 0, 0|U n (t)|r 1 , 0, 0 = ∞ -∞ dE exp(-iEt) exp(ip(r 2 + r 1 -a)) sin(pa) p .
Since the expression for r 2 , 0, 0|U 0 (t)|r 1 , 0, 0 is known from equation ( 14)

πir 2 r 1 M r 2 , 0, 0|U 0 (t)|r 1 , 0, 0 = π ir 2 r 1 M 4π M 2π it 3 exp i M 2t r 2 2 + r 2 1 j 0 M t r 2 r 1 = 2π Mit exp i M 2t r 2 2 + r 2 1 sin M t r 2 r 1 .
one obtains the non-free part of the s wave propagator by the substitutions r 2 → r 2 + r 1a and r 1 → a:

r 2 , 0, 0| U n (t) |r 1 , 0, 0 = - 1 r 2 r 1 2iM πt exp i M 2t ((r 2 + r 1 -a) 2 + a 2 ) × sin M t (r 2 + r 1 -a)a . ( 15
)
4. An exact expression for the partial wave hard sphere propagator in terms of a Laplace transform

For the = 0 partial wave non-free propagators, we resort to numerical computation of the Fourier transform (equations ( 3), ( 9)). The integrand is however of doubly oscillating nature on the real positive axis, and numerical integration can fail, or requires too many digits of precision or computation time. This, of course, depends on the numerical values of the parameters t, a, r 1 , r 2 . For small t values, the difficulty as well as the computation time increases. It is therefore of interest to have another method for comparing numerical results. This other method is based on contour integration and residues theorem. It leads essentially to the numerical computation of a Laplace transform. This numerical method also requires highdigit numbers and a large computation time when t decreases. It is particularly well suitable for large-t values. Anyway, it is useful to have two methods when numerical integration is difficult in order to test accuracy. The expression of the propagator as a Laplace transform is also interesting for study of asymptotic behaviors with respect to different parameters. The Fourier transform is an integral over the interval ]-∞, +∞[. The integration over the interval ]-∞, 0[ will be called the left part of the Fourier transform to be denoted by F L , the integration over the interval ]0, ∞[ will be called its right part to be denoted by F R :

F L = - 1 2πi 0 -∞ dE exp(-iEt) r 2 , , m|G n (E + )|r 1 , , m F R = - 1 2πi +∞ 0 dE exp(-iEt) r 2 , , m|G n (E + )|r 1 , , m .

Contour in the complex energy plane for the left part of the Fourier transform

Figure 1 shows the complex energy plane with its cut along the positive real axis. On the first Riemann sheet, 0 arg(E) < 2π . The close contour AOBA on the left part of figure 1 is on the first Riemann sheet. It consists of the interval [-R, 0] on the negative real axis, the interval [0, -iR] on the negative imaginary axis and a quarter circle of radius R starting at the point -iR and ending at the point -R. As the spherical Bessel functions are entire functions, the only singularities of the integrand ( 9) are possibly at the zeros of h + ((2ME) 1/2 a). The zeros of h + (z) in the z plane have negative imaginary part (see [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]). For = 1, one zero at z = -i. For = 2, 3, . . . , 10, they are shown in figure 2.

With the change of variable given by equation ( 6), the close contour AOBA in the E plane becomes the close contour aOba in the p plane (see figure 3).

It consists of the interval [i √ 2MR, 0] on the positive imaginary axis, the interval 0, √ 2MR exp i 3π contribution of the integration along the circular part goes to zero when R → +∞, the left part of the Fourier transform is equal to the integral along the negative imaginary axis ]-i∞, 0[ on the first Riemann sheet of the E plane, where E = |E| exp i 3π 2 , and therefore

p = √ 2M |E| exp i 3π 4 = - √ 2M |E| exp -i π 4 .
From the symmetry properties of the spherical Bessel function

j (-z) = (-1) j (z) h + (-z) = (-1) +1 h -(z).
With the change of variable E = -iu, one obtains finally an integral on the real positive axis which is a Laplace transform: with

F L = M π ∞ 0 du exp(-ut)z(u)j (z(u)a) - h -(z(u)r 2 )h -(z(u)r 1 ) h -(z(u)
z(u) = √ 2Mu exp -i π 4 . ( 16
)

Contour in the complex energy plane for the right part of the Fourier transform

In order to have a close contour involving the integral on the real positive axis (first sheet of the E plane) with a vanishing contribution of a quarter circle, the imaginary part on the quarter circle must be negative. Since the E plane is cut along the real positive axis, the integrand must be analytically continued into the second Riemann sheet. The closed contour OCDO in figure 1 is the following: OC is the interval [0, R] on the real positive axis (first sheet of the E plane), CB is a quarter circle of radius R in the second Riemann sheet and DO is the purely imaginary interval [-iR, 0] in the second Riemann sheet. The analytical continuation of the integrand into the second Riemann sheet is simply achieved by the equation E = |E| exp(iα) with now -π 2 α 0. With the change of variable given by equation ( 6), the close contour OCDO in the E plane becomes the close contour OcdO in the p plane (see figure 3). It is seen in figures 2 and 3 that for 4, the integrand has poles inside the contour with respect to the variable p. The number of poles increases slowly with : one pole for from 4 to 7, two poles for from 8 to 12, etc. The right part of the Fourier transform is finally equal to minus the integral along the negative imaginary axis ]-i∞, 0[ on the second Riemann sheet of the E plane, where E = |E| exp -i π 2 , and therefore p = √ 2M |E| exp -i π 4 , minus the contribution of the poles. With the change of variable E = -iu, one finally obtains an integral on the real positive axis which is a Laplace transform:

F R = M π ∞ 0 du exp(-ut)z(u)j (z(u)a) h + (z(u)r 2 )h + (z(u)r 1 ) h + (z(u)a) + p 0 residues(f (p))
with z(u) again defined by equation ( 16), and

f (p) = 2 exp -i p 2 2M t p 2 h + (pr 2 )j (pa)h + (pr 1 ) h + (pa) . ( 17
)
The summation over p 0 is over all zeros of h + (pa) with respect to the variable p between the real positive axis and the second bisector (see figure 2). For < 4 , this summation over p 0 is empty and therefore yields 0.

Final results

Let us summarize the results. The partial wave hard sphere propagator for a sphere of radius a is given by r 2 , , m|U(t)|r 1 , , m = r 2 , , m|U 0 (t)|r 1 , , m + r 2 , , m|U n (t)|r 1 , , m (18) where the free part r 2 , , m|U 0 (t)|r 1 , , m is given by equation ( 14) and where the non-free part of the partial wave propagator r 2 , , m|U n (t)|r 1 , , m is given by the Fourier transform:

r 2 , , m|U n (t)|r 1 , , m = - M πi ∞ -∞ dE exp(-iEt) ph + (pr 2 )h + (pr 1 )j (pa) h + (pa) ( 19
)
with p defined by equation ( 6). Another expression for the non-free part of the partial wave propagator is given by a Laplace transform (with residues contributions if 4):

r 2 , , m|U n (t)|r 1 , , m = M π ⎧ ⎪ ⎨ ⎪ ⎩ ∞ 0 du exp(-ut)z(u)j (z(u)a) h + (z(u)r 2 )h + (z(u)r 1 ) h + (z(u)a) - h -(z(u)r 2 )h -(z(u)r 1 ) h -(z(u)a) ⎫ ⎪ ⎬ ⎪ ⎭ + p 0 residues(f (p)) (20)
with z(u) defined by equation ( 16) and f (p) defined by equation (17). The summation over p 0 is over all zeros of h + (pa) whose arguments are in the interval ]-π/4, 0[. For s wave ( = 0), the integration has been done analytically and yields the result given by equation (15).

Numerical results

Comparison with a quantal approximation

Several attempts to approximate analytical expressions for the hard sphere propagator has been made in the past [START_REF] Uhlenbeck | [END_REF]5]. These approximations have been called image approximations. More recently, the CB [7] approximation has superseded image approximations. The CB propagator, noted r 2 | U CB (t) |r 1 is [7]:

r 2 | U CB (t) |r 1 = M i2πt 3/2 exp i M (r 2 -r 1 ) 2 2t - a(u -a) r 2 r 1 exp i MD 2t ( 21
)
with u = r 2 + r 1 D = u 2 -2au + 2a 2 -2a(u -a) cos(α).
The first term involving exp i M(r 2 -r 1 ) 2 2t corresponds to the free propagator, and the non-free part of the CB propagator, to be denoted by r 2 | U CBn (t) |r 1 is thus One deduces from equation ( 12) the non-free part of the partial wave CB propagator:

r 2 | U CBn (t) |r 1 = - M i2πt 3/2 a(u -a) r 2 r 1 exp i MD 2t .
r 2 , , m|U CBn (t)|r 1 , , m = -4π M 2πit 3 a(u -a) r 2 r 1 exp i M(u 2 -2au + 2a 2 ) 2t × (-i) j Ma(u -a) t . ( 22
)
Comparison with equation (15) shows that the CB propagator gives the exact result for the s wave = 0. However, the CB partial wave propagator rapidly becomes inaccurate when increases, as shown in table 1. The modulus of the CB partial wave propagator does not decrease sufficiently rapidly for large values. For example, for = 25, the exact result is

(-3.02 + 1.98i) × 10 -30 whereas the CB result is (2.4 + 17.1i) × 10 -20 .
It is of interest to see how the non-free part of the CB partial wave propagator can be obtained from the exact result (equation ( 19)). One first replaces all h + functions of the integrand by their asymptotic form for large argument: h + (z) exp iz -π 2 /z. This gives a pr 2 r 1 exp i p(r 2 + r 1a) -π 2 . Now one considers the limit form of all h + functions of the integrand for small argument: h + (z) z --1 (2 + 1)!!/(2 + 1). This gives p --1 (a/(r 1 r 2 )) +1 (2 + 1)!!/(2 + 1). Then the replacement of all h + functions of the integrand by a r 2 r 1 (r 2 + r 1a)h + (p(r 2 + r 1a)) ensures the correct asymptotic behavior for large argument and a correct behavior with respect to p, namely proportional to p --1 in the limit of small p values. The Fourier transform then can be evaluated analytically and one obtains equation [START_REF] De Prunelé | [END_REF]. Since h + 0 (z) coincides with its asymptotic form, one understands why the Cao-Berne partial wave propagator is exact for = 0. One understands also why the disagreement increases with since the asymptotic behavior of h + (z) is attained for increasing values of z when increases.

Summation over partial waves (see equation ( 11)) gives the propagator as a function of the angle α between the two vectors r 2 , r 1 (see figure 4). Specifically, the hard sphere propagator has been numerically calculated by adding its free part computed according to equation ( 13), and its non-free part computed by summing the partial wave propagators (equation (11) with index k equal to n). The critical values α c for which the vector r 2 -r 1 is tangent to the hard sphere of radius a are given by the geometrical relation 21)). Points: Van Vleck-Gutzwiller approximation (see equation ( 24)). The vertical line indicates the value of the critical angle α c . The horizontal lines indicate the values of the modulus of the free propagator (see equation ( 13)).

α c = arccos a r 1 + arccos a r 2 . ( 23
We call the domain α < α c the classical domain, and the domain α > α c the classically forbidden domain.

Exact numerical results for the modulus of the hard sphere propagator are compared with the CB approximation in figures 5 and 6. In these figures, the CB approximation reproduces the main qualitative behavior of the exact results. In the classical domain, the oscillations have a very large amplitude and the phase may be incorrect, although the number of oscillations is correct for the examples considered. The exact results exhibit a significant decrease before coming into the classically forbidden domain. This decrease is physically quite natural. The CB approximation does not always decrease sufficiently, and, for example, an unphysical maximum just after α c for t = 1/5 appears. To summarize, the CB approximation has been found globally qualitatively correct but quantitatively inaccurate. An interesting phenomenon exhibited by the exact results is a local maximum for α = π . This effect is particularly visible in figure 6. This counter intuitive effect from a classical view point is a purely quantal reminiscent of a similar effect in optics: central bright spot (called Maraldi or Poisson or Arago spot) at the centre of the geometrical shadow in Fresnel diffraction by a circular obstacle [START_REF] Ditchburn | Light[END_REF][START_REF] Born | Principles of Optics[END_REF].

Classical arguments and limit cases

How rapidly does the summation over partial wave converge (see equation (11))? An estimation for the maximum value M to be used for truncating the series (11) is obtained by looking at a majorant for the classical angular momentum compatible with M, a, r 1 , r 2 , t values. For fixed values of the velocity v, the angular momentum is maximum at maximal impact parameter a, when the trajectory is tangent to the hard sphere. The modulus of the angular momentum in this configuration is the product aMv. For a given time t, the velocity is maximum for the maximum distance along trajectories compatible with r 1 and r 2 . The trajectory distance between r and the hard sphere is maximum for the tangential trajectory, and is then √ r 2a 2 . A classical majorant for angular momentum is then:

c = aM r 2 1 -a 2 + r 2 2 -a 2 t .
One therefore expects that the contribution to the sum of the series (11) should be negligible for greater than about

M = α aM(r 1 + r 2 ) t + β.
The choice for the values of the coefficients α, β depends on the precision required. The choice α = 3/2 and β = 5 has been found largely sufficient for the precision of drawing figures 5 and 6. The reason for the addition of the number β is that for small quantum numbers, classical argument may fail. For the cases considered in figure 5, M is respectively equal to 5.3, 6.5, 12.5 and 42.5, for t respectively equal to 25, 5, 1 and 1/5. The values of above which the visual appearance of the graph does not change are 2, 2, 7, 30, for t respectively equal to 25, 5, 1 and 1/5. The calculations have been made up to equal to 7, 15, 25 and 50 for t respectively equal to 25, 5, 1 and 1/5 and it has been verified that the decrease is very rapid above M . For the case considered in figure 6, M is equal to 20 and it has been verified that the visual appearance of the graph does not change as long as the truncation value is equal to or above 13. If one requires increasing the relative accuracy, the majorant must of course also be increased.

Physically, an important result is that in the limit a → 0 or t → ∞ or M → 0, only the s and p waves contribute significantly. (Of course the limit M → 0 is not a physical one.) What about the opposite cases, the short time limit ? It will be seen in the next section that the classical result of Van Vleck-Gutzwiller provides an excellent approximation in the classical domain.

Comparison with a classical approximation

In the classical region α α c (see equation ( 23)), or otherwise stated, for r 2 outside the geometric shadow of the hard sphere of radius a, the VG approximation for the hard sphere quantum propagator, to be denoted by r 2 | U V G (t) |r 1 is [18,[START_REF] Gutzwiller | Chaos in Classical and Quantum Mechanics[END_REF]

r 2 | U V G (t) |r 1 = r 2 | U 0 (t) |r 1 - 1 (2π i) 3/2 det - ∂ 2 S ∂r 2 ∂r 1 exp(iS) (24)
In equation ( 24), r 2 | U 0 (t) |r 1 is the free propagator (see equation ( 13)). S is the classical action, the integral over time of the Lagrangian L:

r 2 r 1 Ldt = r 2 r 1 p 2
2M dt for the real classical trajectory starting at time 0 at the point r 1 and ending at time t at the point r 2 after one collision with the hard sphere of radius a. (The present action S is called the Hamilton's principal function in [START_REF] Gutzwiller | Chaos in Classical and Quantum Mechanics[END_REF][START_REF] Brack | Semiclassical Physics[END_REF] , where it is noted R. In these references, the name action and the letter S refer to another function.) The minus sign on the right-hand side of equation [START_REF] Messiah | [END_REF] originates from the reflection at the point of collision with the hard sphere. The action S is thus a function of the variables M, t, a, r 2 , r 1 . Finally ∂ 2 S ∂r 2 ∂r 1 denotes the matrix of order 3 whose ij element is the second derivative

∂ 2 S ∂(x 2 ) i ∂(x 1
) j of the action with respect to the coordinates of the end points r 2 , r 1 of the classical trajectories. As emphasized in [START_REF] Gutzwiller | Chaos in Classical and Quantum Mechanics[END_REF], the burden of computing explicitly r 2 | U V G (t) |r 1 lies completely within classical mechanics, and the VG approximation can appropriately also be called the classical propagator. Specifically

S = M 2t l 2 (25)
with l the length of the classical trajectory connecting the points r 1 , r 2 with one collision with the hard sphere. In the applications to be discussed now, the length l has been determined numerically by minimization, and the second derivative of S has been determined numerically according to the relation [29]

∂ 2 f ∂x∂y (x, y) f (x + b, y + b) -f (x -b, y + b) -f (x + b, y -b) + f (x -b, y -b) 4b 2
The accuracy has been tested by changing the values of b. For example, for b = 1/100 and b = 1/10, the relative change of all points in figures 5 and 6 (to be defined as the ratio between the result with b = 1/10 minus the result with b = 1/100 and the result with b = 1/100) is less than 5 × 10 -4 . It is seen in figure 5 that the accuracy of the VG approximation increases when t decreases. This accuracy is excellent for t = 1/5, except when approaching the classically forbidden region (α α c ). This also holds for figure 6. This behavior with respect to time was expected since the classical action increases when t decreases or when l increases (see equation [START_REF] Brack | Semiclassical Physics[END_REF]). The action S increases with α. For figure 5, the minimum values of S, for α = 0, are 0.18h, 0.9h, 4.5h, 22.5h for t = 25, 5, 1, 1/5, respectively. For figure 6, the minimum value of S, for α = 0, is 32h. Path integral formulation of quantum mechanics [START_REF] Feynman | [END_REF][START_REF] Feynman | Quantum Mechanics and Path Integrals[END_REF][START_REF] Zinn-Justin | Intégrale de Chemin en Mé canique Quantique: Introduction[END_REF] tells us that for a very large value of the classical action, the sum over all paths (classical and non-classical) is dominated by the classical path due to the stationary value of the action for this path. It is therefore natural that the accuracy of the VG propagator increases when S increases (at a fixed value of α). What is remarkable is the degree of accuracy for such low values of the action S, only a few units of h. This is however not quite a surprise since the exact value of the free propagator (see equation ( 13)) coincides with the free VG propagator, and the real classical trajectory with one collision with the hard sphere is nothing but two free trajectories (one before the collision, one after).

Only the moduli of propagators are reported in figures 5 and 6, but when the modulus of the approximate propagators agrees with the modulus of the exact result, the real and imaginary parts also agree.

Conclusion

By an explicit determination of the matrix elements of the hard sphere resolvent (see equation ( 4)), an exact expression for the partial wave hard sphere quantum propagator in a three-dimensional space has been given in terms of a Fourier transform of an expression involving spherical Bessel functions. For the s wave case, the Fourier transform has been calculated analytically. This Fourier transform has also been evaluated for the general case by contour integration in terms of a Laplace transform and residues contribution. This Laplace transform is particularly suitable for numerical computation for the large time value t (small values of the classical action S (see equation ( 25)). For small t values (more generally for large values of the classical action S in units of h), numerical calculations by both methods (Fourier and Laplace) are difficult, and require a large number of digits for precision and a large computation time.

The calculation of the hard sphere quantum propagator by summation over partial wave contribution has shown an excellent accuracy of the Van Vleck-Gutzwiller approximation in the classical domain α < α c for α not in the vicinity of α c , and when the classical action S (see equation [START_REF] Brack | Semiclassical Physics[END_REF]) is greater than only a few units of h. In the classically forbidden domain (α > α c ), the exact result rapidly decreases as expected, but exhibits a small local maximum for α = π , i.e. for the point r 2 on the other side of the sphere, so that the three points r 1 ,

4 on the second bisector and part of a circle ofFigure 1 .Figure 2 .

 412 Figure 1. Contours in the cut energy plane E. The dashed line curve lies on the second Riemann sheet.

Figure 3 .

 3 Figure 3. Contours in the p complex plane, corresponding to the contours of figure 1.

) cFigure 4 .Figure 5 .

 45 Figure 4. The origin of the two vectors r 1 , r 2 is the center of a hard sphere of radius a = 1. α is the angle between r 1 , r 2 . The norms of the two vectors are constant, r 1 = 2, r 2 = 3. The extremal values of |r 2r 1 | are 1 (parallel case) and 5 (antiparallel case). When r 2r 1 becomes tangent to the hard sphere of radius a, the angle α has the critical value α c above which no classical trajectory exists.

1 Figure 6 .

 16 Figure 6. Modulus of the hard sphere quantum propagator as a function of the angle α between the two vectors r 1 , r 2 for a = 1, r 1 = r 2 = 5, M = 1, t = 1. Full curve: exact result. Dashed curve: Cao-Berne approximation (See equation (21)). Points: Van Vleck-Gutzwiller approximation (See equation (24)). The vertical line indicates the value of the critical angle α c . The horizontal lines indicate the values of the modulus of the free propagator (see equation (13)).

Table 1 .

 1 First line: , exact non-free part of the partial wave propagator r 2 , , m|U n (t)|r 1 , , m (see equations[START_REF] Gutzwiller | Chaos in Classical and Quantum Mechanics[END_REF], (20)), non-free part of the CB propagator r 2 , , m|U CBn (t)|r 1 , , m (see equation[START_REF] De Prunelé | [END_REF]), for M = 1, t = 1, a = 1, r 1 = 2, r 2 = 3.

	3, , m|U n (1)|2, , m	3, , m|U CBn (1)|2, , m
	0 0.0997 -0.0140 i	0.0997 -0.0140 i
	1 0.0122 + 0.0565 i	0.0086 + 0.0612 i
	2 0.1216 -0.0560 i	0.146 -0.020 i
	3 -0.0746 -0.0764 i	-0.0169 -0.1208 i
	4 -0.0165 + 0.0536 i	-0.0658 + 0.0092 i
	5 0.0206 -0.0085 i	0.0038 + 0.0273 i
	6 -0.00676 -0.00168 i	0.0162 -0.0023 i
	7 0.00140 + 0.00095 i	-0.00037 -0.00263 i
	8 -0.000239 -0.000174 i	-0.000652 + 0.000091 i
	9 0.0000307 + 0.0000116 i	0.000020 + 0.000143 i

Acknowledgments

The author thanks the referees for drawing his attention to reference [27].

the center of the sphere and r 2 are aligned. This effect is reminiscent of Arago bright spot in Fresnel diffraction [START_REF] Ditchburn | Light[END_REF][START_REF] Born | Principles of Optics[END_REF].