Umedjon Khalilov

Annemie Bogaerts

Shahzad Hussain

Eva Kovacevic

Pascal Brault

Chantal Boulmer-Leborgne

Erik C Neyts

Nanoscale mechanisms of CNT growth and etching in plasma environment

Keywords: carbon nanotube, Ni nanocatalyst, PECVD, computer simulations, electron microscopy, H 2 plasma, hydrocarbons, etching

Plasma-enhanced chemical deposition (PECVD) of carbon nanotubes has already been shown to allow chirality control to some extent. In PECVD, however, etching may occur simultaneously with the growth, and the occurrence of intermediate processes further significantly complicates the growth process.

We here employ a computational approach with experimental support to study the plasmabased formation of Ni nanoclusters, Ni-catalyzed CNT growth and subsequent etching processes, in order to understand the underpinning nanoscale mechanisms. We find that hydrogen is the dominant factor in both the re-structuring of a Ni film and the subsequent appearance of Ni nanoclusters, as well as in the CNT nucleation and etching processes. The obtained results are compared with available theoretical and experimental studies and provide a deeper understanding of the occurring nanoscale mechanisms in plasma-assisted CNT nucleation and growth.

Introduction

Ever since the discovery of carbon nanotubes (CNT) [START_REF] Iijima | Helical microtubules of graphitic carbon[END_REF][START_REF] Iijima | Single-shell carbon nanotubes of 1-nm diameter[END_REF][START_REF] Bethune | Cobalt-catalysed growth of carbon nanotubes with single-atomic-layer walls[END_REF], obtaining control over its chirality-dependent properties has been a major topic of research. Such control is needed to allow their commercial application in full [START_REF] Baughman | Carbon Nanotubes-the Route Toward Applications[END_REF]. While the precise control of the obtained structure is still a daunting task, significant achievements have already been attained both experimentally and computationally [START_REF] Tessonnier | Recent Progress on the Growth Mechanism of Carbon Nanotubes: A Review[END_REF][START_REF] Jourdain | Current understanding of the growth of carbon nanotubes in catalytic chemical vapour deposition[END_REF][START_REF] Neyts | PECVD growth of carbon nanotubes: From experiment to simulation[END_REF].

Among other growth techniques [START_REF] Iijima | Single-shell carbon nanotubes of 1-nm diameter[END_REF][START_REF] Bethune | Cobalt-catalysed growth of carbon nanotubes with single-atomic-layer walls[END_REF][START_REF] Thess | Crystalline Ropes of Metallic Carbon Nanotubes[END_REF][START_REF] Li | Large-Scale Synthesisi of Aligned Carbon Nanotubes[END_REF], thermal chemical vapor deposition (TCVD) and plasma enhanced chemical vapor deposition (PECVD) show several advantages [START_REF] Jourdain | Current understanding of the growth of carbon nanotubes in catalytic chemical vapour deposition[END_REF][START_REF] Kato | Growth of Single-Walled Carbon Nanotubes by Plasma CVD[END_REF].

Especially, PECVD allows low-temperature growth, which is envisaged as a means towards chirality control and a reduced chirality distribution with respect to thermal growth [START_REF] Neyts | PECVD growth of carbon nanotubes: From experiment to simulation[END_REF][START_REF] Hatakeyama | Plasma-synthesized single-walled carbon nanotubes and their applications[END_REF][START_REF] Ghorannevis | Narrow-Chirality Distributed Single-Walled Carbon Nanotube Growth from Nonmagnetic Catalyst[END_REF]. Furthermore, plasma pre-treatment of the catalyst can greatly modify the resulting nanotubes [START_REF] Cantoro | Plasma restructuring of catalysts for chemical vapor deposition of carbon nanotubes[END_REF][START_REF] Esconjauregui | Growth of high-density vertically aligned arrays of carbon nanotubes by plasma-assisted catalyst pretreatment[END_REF]. So far, however, the nature of the catalyst formation process by annealing metal thin films (such as Ni, Co or Fe) with or without plasma is not fully understood.

Particularly, the reason and mechanism of the appearance of small and uniform nanoparticles in H 2 plasma is still not clear [START_REF] Robertson | Controlling the Catalyst During Carbon Nanotube Growth[END_REF].

Besides the size control of the catalyst nanoparticles, several other factors are also of considerable importance, including the carbon chemical potential [START_REF] Amara | A tight-binding potential for atomistic simulations of carbon interacting with transition metals: Application to the Ni-C system[END_REF][START_REF] Amara | Understanding the nucleation mechanisms of carbon nanotubes in catalytic chemical vapor deposition[END_REF], the necessity of a carbide phase [START_REF] Gómez-Gualdrón | Characterization of carbon atomistic pathways during single-walled carbon nanotube growth on supported metal nanoparticles[END_REF][START_REF] Ohta | Rapid Growth of a Single-Walled Carbon Nanotube on an Iron Cluster: Density-Functional Tight-Binding Molecular Dynamics Simulations[END_REF][START_REF] Page | QM/MD Simulation of SWNT Nucleation on Transition-Metal Carbide Nanoparticles[END_REF][START_REF] Börjesson | First Principles Studies of the Effect of Nickel Carbide Catalyst Composition on Carbon Nanotube Growth[END_REF], the importance of metal-mediated defect healing [START_REF] Diarra | Role of defect healing on the chirality of single-wall carbon nanotubes[END_REF][START_REF] Neyts | Changing Chirality during Single-Walled Carbon Nanotube Growth: A Reactive Molecular Dynamics/Monte Carlo Study[END_REF], the influence of growth precursor, etc. All of these factors have been investigated to some extent, but a clear and consistent picture has not yet emerged. Also, in PECVD-based CNT growth, the role of C atoms or C 2 dimers [START_REF] Gómez-Gualdrón | Effects of Precursor Type on the CVD Growth of Single-Walled Carbon Nanotubes[END_REF] or energetic C-ions [START_REF] Shariat | Plasma enhanced growth of single walled carbon nanotubes at low temperature: A reactive molecular dynamics simulation[END_REF][START_REF] Shariat | On the low-temperature growth mechanism of single walled carbon nanotubes in plasma enhanced chemical vapor deposition[END_REF] as input gas, as well as the influence of Ar + bombardment [START_REF] Shariat | Plasma enhanced growth of single walled carbon nanotubes at low temperature: A reactive molecular dynamics simulation[END_REF][START_REF] Neyts | Defect Healing and Enhanced Nucleation of Carbon Nanotubes by Low-Energy Ion Bombardment[END_REF], have been computationally studied. Nevertheless, the role of the plasma-generated growth precursor in CNT growth, and in particular the hydrocarbon C n H m species in the incubation stage of CNT nucleation [START_REF] Khalilov | Atomic scale simulation of carbon nanotube nucleation from hydrocarbon precursors[END_REF], is not yet properly understood.

Under plasma conditions, the growth kinetics of CNTs can be expressed as a balance between the growth and etching processes [START_REF] Kato | Growth of Single-Walled Carbon Nanotubes by Plasma CVD[END_REF]. Besides its damaging effect on a growing CNT structure, plasma-assisted etching can also be used to separate semiconductor CNTs from their metallic counterparts [START_REF] Zhang | Selective Etching of Metallic Carbon Nanotubes by Gas-Phase Reaction[END_REF]. Overall, however, the etching onset is poorly studied during the CNT nucleation from plasma species, including H-radicals and ions. However, the ions are not direct sources in both CNT growth/etching, in contrast to the H-radical and hydrocarbon species [START_REF] Kato | Formation of Freestanding Single-Walled Carbon Nanotubes by Plasma-Enhanced CVD[END_REF].

Overall understanding of these processes, i.e., re-structuring of the metal catalyst, and CNT nucleation and subsequent etching, require a careful investigation of the intermediate processes taking place in nanoseconds, which is difficult to attain experimentally. On the other hand, simulations are typically computationally expensive if they are to address the inherent time and length scales. Nevertheless, both obstacles may be overcome by pursuing a combined modeling/experimental approach addressing the aforementioned issues. In the current contribution, we investigate both subjects by the modeling of the nanocatalyst formation and subsequent CNT nucleation/etching, and by performing experiments to obtain nanoclusters and CNT carpets in both thermal and plasma mode.

Methods

Simulation details

The simulated deposition of Ni atoms onto a semiconductor substrate is performed using the LAMMPS open source code [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF]. In the simulations, two Si crystals with a size of 5.43 nm × 5.43 nm × 2.6 nm and 8.15 nm × 8.15 nm × 2.6 nm are used as a substrate, containing 3995 and 8888 silicon atoms, respectively. In both cases, the substrate temperature is equilibrated to 300 K, controlled by a Berendsen thermostat [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] with a damping constant of 0.1 ps. The atoms of the bottom layer are fixed. During the deposition, a Ni atom is randomly positioned at least 1 nm above the Si(100) surface every 0.5 ps : this value prevents thermal relaxation interactions between consecutive Ni impacts on the Si surface [START_REF] Xie | Molecular Dynamics simulations of clusters and thin film growth in the context of plasma sputtering deposition[END_REF]. The Ni initial velocities are randomly and uniformly chosen in the range 0.05-0.3 nm/ps towards the surface. Tersoff [START_REF] Tersoff | New empirical approach for the structure and energy of covalent systems[END_REF] and MEAM potentials [START_REF] Baskes | Modified embedded-atom potentials for cubic materials and impurities[END_REF] are used for the description of Si-Si and Ni-Ni interactions, respectively. Interactions between Ni and Si atoms are described using Lennard-Jones (LJ) potentials [START_REF] Xie | AlCoCrCuFeNi high entropy alloy cluster growth and annealing on silicon: A classical molecular dynamics simulation study[END_REF].

For the LAMMPS-based simulations of the re-structuring of the obtained thin Ni-film, a wider Si structure with a surface area of 8.15 nm × 8.15 nm and with a reduced thickness of 1.69 nm is chosen, to avoid the cell size effect and shorten the simulation time [START_REF] Maekawa | Dewetting dynamics of nickel thin film on alpha-quartz substrate: A molecular dynamics study[END_REF]. The Ni thin film, with a thickness in the range of 0.3-1.2 nm, is annealed at temperatures of 600 K, 800 K, 1000 K and 1200 K. Ni-Ni, Si-Si and Ni-Si interactions are described using the Tersoff [START_REF] Tersoff | New empirical approach for the structure and energy of covalent systems[END_REF], MEAM [START_REF] Baskes | Modified embedded-atom potentials for cubic materials and impurities[END_REF] and LJ potentials [START_REF] Xie | AlCoCrCuFeNi high entropy alloy cluster growth and annealing on silicon: A classical molecular dynamics simulation study[END_REF], respectively. Note that the LJ pair potential does not allow to correctly describe physical / chemical processes such as Ni-atom diffusion into the Si substrate or the formation of Ni-silicide [START_REF] Alberti | Role of the early stages of Ni-Si interaction on the structural properties of the reaction products[END_REF][START_REF] Thron | Structural changes during the reaction of Ni thin films with (100) silicon substrates[END_REF][START_REF] Bakovets | Mechanism of Ni Film CVD with a Ni(Ktfaa) 2 precursor on a Silicon Substrate[END_REF]. Therefore, to avoid an unphysical diffusion of Ni atoms and an evaporation of Si atoms at high temperatures, the substrate atoms are immobilized during the simulation.

Combined molecular dynamics / time-stamped force-bias Monte Carlo (MD/tfMC) simulations are used for simulation of the catalyzed CNT growth and hydrogenation/etching processes [START_REF] Bal | On the time scale associated with Monte Carlo simulations[END_REF][START_REF] Neyts | Catalyzed Growth of Carbon Nanotube with Definable Chirality by Hybrid Molecular Dynamics Force Biased Monte Carlo Simulations[END_REF][START_REF] Neyts | Establishing Uniform Acceptance in Force Biased Monte Carlo Simulations[END_REF]. In tfMC calculations, all the atoms in the system are displaced at once in every time step with unit probability, and thus generate a system evolution in a MD-like fashion. This technique can therefore very efficiently be coupled to canonical MD simulations [START_REF] Bal | On the time scale associated with Monte Carlo simulations[END_REF]. Both during the MD and the tfMC cycles, ReaxFF is applied to describe the bond dissociation and formation processes, i.e., the C-C bond dissociation, dehydrogenation, rehydrogenation and H 2 formation, as well as other reactions [START_REF] Van Duin | ReaxFF: A Reactive Force Field for Hydrocarbons[END_REF]45]. During the simulations, the Ni 55 nanocatalyst is physisorbed on a virtual Al or Si substrate, employing a z-integrated Lennard-Jones potential [START_REF] Khalilov | Microscopic mechanisms of vertical graphene and carbon nanotube cap nucleation from hydrocarbon growth precursors[END_REF]. Subsequently, the plasma-generated carbon precursor (C, C 2 , C 2 H, CH 3 , CH 4 , C 2 H 2 , C 6 H 6 , C 7 H 8 , and C 8 H 10 ) is allowed to impinge on the cluster and its total density is kept constant in the gas-phase. In the case of CNT etching, also H-radicals are allowed to impinge on the tube wall. While these plasma species can have a wide energy range from thermal to hyperthermal [START_REF] Hash | Model based comparison of thermal and plasma chemical vapor deposition of carbon nanotubes[END_REF], we choose thermal hydrocarbon species and Hradicals as plasma species in the simulations. When a gas-phase particle adsorbs on the nanocluster (or the nanotube), the resulting structure is allowed to relax by application of tfMC [START_REF] Bal | On the time scale associated with Monte Carlo simulations[END_REF][START_REF] Neyts | Establishing Uniform Acceptance in Force Biased Monte Carlo Simulations[END_REF]. During the relaxation, no new plasma species are allowed to impinge on the structure. All simulations are performed at a temperature in the range of 1000 K -2000 K controlled by the canonical Bussi thermostat [START_REF] Bussi | Canonical sampling through velocity-rescaling[END_REF].

Finally, in order to clarify the role of H 2 -plasma on the annealing of Ni-catalyst clusters, the interactions of two Ni nanoclusters are simulated using the combined reactive MD and tfMC technique in two regimes: [START_REF] Iijima | Helical microtubules of graphitic carbon[END_REF] annealing without H 2 plasma and (2) annealing with H 2 plasma. In the simulations, Ni-Ni and Ni-H interactions are described using the ReaxFF potential, as parametrized by Mueller et al. [START_REF] Van Duin | ReaxFF: A Reactive Force Field for Hydrocarbons[END_REF]45]. Prior to the annealing, two adjacent Ni 55 nanoclusters are initially positioned at least 1 nm from each other. Both nanoclusters are physisorbed on a virtual Si substrate, employing a z-integrated Lennard-Jones potential [START_REF] Khalilov | Microscopic mechanisms of vertical graphene and carbon nanotube cap nucleation from hydrocarbon growth precursors[END_REF].

Structures in both regimes are thermalized between 300 K and 800 K using the canonical Bussi thermostat [START_REF] Bussi | Canonical sampling through velocity-rescaling[END_REF].

Experimental procedure

The detailed description of the PECVD reactor and the growth of multi-wall CNTs (MWNTs) was reported in our previous studies [START_REF] Labbaye | In situ Raman spectroscopy for growth monitoring of vertical aligned multiwall carbon nanotubes in plasma reactor[END_REF][START_REF] Gaillard | High-resolution transmission electron microscope observations of multiwalled carbon nanotube microstructures grown by plasma enhanced chemical vapor deposition[END_REF]. In these experimental studies, the growth of MWNTs is achieved on a Si substrate using Ni nanocatalyst by PECVD in three steps. In the first step, a 10 nm thin film of Ni is deposited onto the substrate by using a precision etching and coating system (PECS). Afterwards, the sample is transferred into the PECVD reactor (radio frequency capacitively coupled plasma (RF CCP) working at 13.56 MHz) and mounted on the grounded electrode. The 2 nd step started once the base pressure inside the reactor is below 4•10 -4 Pa. In this step annealing with and without RF CCP treatment of the catalyst thin film is done at 923 K with a ramp time of 1200 s in a reducing atmosphere of hydrogen (10 sccm) at a pressure of 0.25 mbar and maintaining these conditions for 600 s. During 2 nd step, the electrode distance is of about 2cm. In the 3 rd step, the H 2 flow was increased to 40 sccm, the growth precursor C 2 H 4 (20 sccm) is introduced, the pressure is increased up to 1 mbar, electrode distance was increased to 6 cm, and growth of CNTs is performed during 1800 s. In the initial stage, incoming Ni atoms weakly adsorb on the Si surface, freely move over the surface and eventually find other Ni adatoms to form a small Ni cluster (Fig. 1b), in agreement with experimental evidence [START_REF] Niwa | Deposition Mechanism of Ni on Si(100) Surfaces in Aqueous Alkaline Solution[END_REF]. During the deposition process, small clusters quickly sinter due to the weak interaction between the Ni-cluster and the Si surface (Fig. 1c) [START_REF] Ghosh | Ultrathin Metal Transparent Electrodes for the Optoelectronics Industry, Chapter 2 Basics of Ultrathin Metal Films and Their Use as Transparent Electrodes[END_REF]. As time advances, large clusters and subsequently islands gradually appear on the substrate (Fig. 1d-h) before Ni completely covers the whole Si surface as a thin film (Fig. 1i).

Results and discussion

The simulations thus show Ni-thin film growth through the Volmer-Weber (island) growth mode [START_REF] Volmer | Nucleus formation in supersaturated systems[END_REF]. The obtained Ni film (with a thickness about 1.2 nm in Fig. 1i) is amorphous and the surface roughness is low, in agreement with experimental observations [START_REF] Cantoro | Plasma restructuring of catalysts for chemical vapor deposition of carbon nanotubes[END_REF][START_REF] Ghosh | Ultrathin Metal Transparent Electrodes for the Optoelectronics Industry, Chapter 2 Basics of Ultrathin Metal Films and Their Use as Transparent Electrodes[END_REF]. In Fig. 2, snapshots during the restructuring of the Ni film with thickness of 0.3 nm at a temperature of 1000 K and 1200 K are shown. When the annealing begins, the Ni thin film breaks due to its dewetting [START_REF] Maekawa | Dewetting dynamics of nickel thin film on alpha-quartz substrate: A molecular dynamics study[END_REF][START_REF] Krishna | Self-Organization of Nanoscale Multilayer Liquid metal Films: Experiment and Theory[END_REF][START_REF] Peláez | Dynamics of laser induced metal nanoparticle and pattern formation[END_REF] and vacant holes start to appear on the Si surface after 8 ps and 5 ps at 1000 K and 1200 K, respectively. The area of the bare Si surface subsequently expands (see Fig. 2, 16 ps -72 ps at 1000 K and 30 ps -65 ps at 1200 K) and vacant holes coalesce (see Fig. 2, 96 ps and 85 ps at 1000 K and 1200 K, respectively). Subsequently, the Ni thin film eventually converts to a Ni nanocluster (Fig. 2, last panels for both temperatures). The snapshots demonstrate that the dewetting process strongly depends on the annealing temperature and the thickness of the Ni thin film, as well as the adhesion of the film-surface interface, which is in qualitative agreement with experimental observations [START_REF] Cantoro | Plasma restructuring of catalysts for chemical vapor deposition of carbon nanotubes[END_REF]56].

Due to the Gibbs-Tomson effect, i.e., the melting temperature depression in nanoscale systems [START_REF] Buffat | Size effect on the melting temperature of gold particles[END_REF][START_REF] Shibuta | A molecular dynamics study of the phase transition in bcc metal nanoparticles[END_REF], the physical state of small cylindrical structures (indicated by the arrow in Fig. 3) changes (i.e., its fluidity increases). During the dewetting process, such small structures break up or collapse easily (Fig. 3, cf. 72 ps and 78 ps) due to the capillary effect, which can be explained by the Plateau-Rayleigh instability [START_REF] Maekawa | Dewetting dynamics of nickel thin film on alpha-quartz substrate: A molecular dynamics study[END_REF][START_REF] Plateau | Statique Experimentale et Theorique des Liquides Soumis aux Seules Forces Moleculaires[END_REF][START_REF] Rayleigh | [END_REF]. According to this theory, the cylindrical fluid breaks up when its length (l) satisfies the condition l > 2πr, where r is the radius of the cylindrical fluid. In our case, the radius and length of the small cylindrical (fluid-like) structure shown in Fig. 3 Indeed, the Plateau-Rayleigh instability explains the breakup of small structures observed in the simulations (see e.g. Fig. 3, 78 ps). Furthermore, our simulations reveal that the breakup speed depends on the surface temperature and inversely depends on the film-surface adhesion.

In our experiments, both small and large Ni clusters as well as Ni islands appear on the surface during the restructuring process (compare Fig. 2 and Fig. 4a). The number of larger clusters is seen to increase and the small clusters are seen to simultaneously disappear during the annealing process, which can be explained by Ostwald ripening [61]. According to this theory, the small clusters dissolve due to the minimization of the surface energy and their powders (or nanoclusters) redeposit onto the large cluster, which is energetically (thermodynamically) favored over the small clusters [START_REF] Voorhees | The Theory of Ostwald Ripening[END_REF]. Consequently, only small and large clusters remain on the surface due to fast incorporation of the kinetically favored powders (nanoclusters), as shown in Fig. 4a. Indeed, the Ostwald ripening explains well these processes [START_REF] Hansen | Sintering of Catalytic Nanoparticles: Particle Migration or Ostwald Ripening?[END_REF]. However, this theory fails to explain the remaining nanoclusters in the case of annealing with H 2 plasma as shown in Fig. 4b. To understand this phenomenon, MD/tfMC simulations were performed, taking the following two observations into account: (i) the concentration of H radicals in a H 2 plasma is significantly higher than their concentration in TCVD [START_REF] Hash | Model based comparison of thermal and plasma chemical vapor deposition of carbon nanotubes[END_REF][START_REF] Behr | Hydrogen etching and cutting of multiwall carbon nanotubes[END_REF], and (ii) reactive H radicals/atoms bind strongly with Ni atoms and therefore nickel hydride is stable at low temperatures [START_REF] Mueller | Structures, Energetics, and Reaction for CH x Bound to the nickel (111) Surface[END_REF]. In the simulations, H atoms are added to the gas phase and their number is kept constant during the annealing in the plasma case (Fig. 5b). Furthermore, hydrogenated Ni 55 (h-Ni) nanoclusters are used in the annealing process, in agreement with the second observation.

Structures in both regimes are thermalized at a temperature from 300 K up to 800 K. At the end of the annealing process, the nanoclusters coalesce in the first (thermal) regime (Fig. 5a), while the h-Ni nanoclusters do not sinter on the surface in the second (plasma) regime. In the annealing-only regime, the clusters are connected through surface Ni atoms, which have dangling bonds. Such bonds, however, are terminated by reactive H atoms (radicals) in the "annealing + plasma" regime even at low temperature (Fig. 5b). As a result, the interaction between adjacent Ni hydride nanoclusters becomes very weak and they cannot sinter. The results show that at low temperature (300 K), the nanoclusters do not sinter/coalesce, regardless of the presence or absence of H 2 , due to their low mobility and weak mutual interaction. At high temperature (800 K), however, the coalescence depends on the presence of H 2 , as at this temperature their mobility is much higher. Note also that when the temperature continues to increase, the sintering speed gradually enhances due to the temperature-dependent H 2 desorption and consequently the removal of the hydrogen shell from these nanoclusters.

In general, the obtained results indicate that the Ostwald ripening discontinues or retards due to the H-effect in the second regime. Due to this effect, Ni nanoparticles (powders) can remain stable for a longer time on the surface in the plasma regime. Indeed, the obtained results allow us a better understanding the difference of the re-structuring nature of Ni thin film in the annealing with H 2 plasma. Finally, the obtained Ni-nanoclusters are used as a catalyst in our simulations for the CNT nucleation from plasma species.

CNT nucleation and etching

2.1 CNT nucleation from plasma-generated species. From the simulations, we distinguish three main stages in the catalyzed nucleation of carbon nanotubes from hydrocarbon C n H m plasma species: (i) incubation, (ii) cap formation and (iii) continued tube growth. Also, these stages are discussed in our previous work, devoted to the CNT nucleation from hydrocarbon molecules in a TCVD setup [START_REF] Khalilov | Atomic scale simulation of carbon nanotube nucleation from hydrocarbon precursors[END_REF].

In the beginning of stage I, gas-phase hydrocarbon species impinge on the cluster and they gradually dissociate by virtue of the Ni nanocatalyst. Subsequently, free H-adatoms move around the nanocluster surface and find another H-adatom, resulting in H 2 desorption. The remaining C adatoms subsequently dissolve into the Ni cluster. This "C n H m adsorption/dissociation, C dissolution and H 2 desorption" scenario continues until the nanocluster is saturated (Fig. 6a, frame 1). Subsequently, segregated C-atoms diffuse over and through the catalyst (sub)surface until recombining with other carbon atom(s) in order to form initial surface carbon structures, i.e., C 2 dimers, sp 2 and sp 3 C atoms and short C y polyyne chains without or with H-terminations (Fig. 6a, frame 2). When rehydrogenation occurs, the surface carbon structures, including polyyne chains and the graphene-like network (Fig. 6a, frame 3) become partially hydrogenated and vertically freestanding graphene sheets can be formed on the catalyst surface (Fig. 6a, frame 4). Due to competition between dehydrogenation and rehydrogenation, both horizontally and vertically oriented carbon patches can be found on the Ni surface (cf. Fig. 6a, frames 3 and 4) [START_REF] Khalilov | Microscopic mechanisms of vertical graphene and carbon nanotube cap nucleation from hydrocarbon growth precursors[END_REF]. The results show that this competition depends on the growth temperature, type of hydrocarbon species and their gas-phase pressure (or flux), as well as the catalyst-surface adhesion [START_REF] Khalilov | Atomic scale simulation of carbon nanotube nucleation from hydrocarbon precursors[END_REF].

Such competition can explain the appearance of unstable carbon structures in the incubation stage of CNT growth, as observed in the in-situ transmission electron microscopy (TEM) observations by Yoshida et al. [START_REF] Yoshida | Atomic-Scale In-situ Observation of Carbon Nanotube Growth from Solid State Iron Carbide Nanoparticles[END_REF]. As time advances, the vertical free-standing networks invariably and irreversibly transform to horizontal networks covering the surface. Also, we found in previous investigations that such vertical carbon networks can remain stable on the surface, and they stack up parallel at low temperatures [START_REF] Khalilov | Atomic scale simulation of carbon nanotube nucleation from hydrocarbon precursors[END_REF][START_REF] Khalilov | Microscopic mechanisms of vertical graphene and carbon nanotube cap nucleation from hydrocarbon growth precursors[END_REF]. The appearance of such structures could also explain the onset of MWNT formation [START_REF] Khalilov | Atomic scale simulation of carbon nanotube nucleation from hydrocarbon precursors[END_REF], as suggested earlier by Hoffman et al. [START_REF] Hofmann | In situ Observations of Catalyst Dynamics during Surface-Bound Carbon Nanotube Nucleation[END_REF]. At the end of the first growth stage, the graphitic network(s) start(s) to gradually cover the nanocluster surface.

The second stage of the CNT nucleation is characterized by the appearance of a carbon cap on the cluster (Fig. 6a, frame 5), which is essential for the nucleation of a CNT. The simulation results confirm our previous conclusions that cap nucleation in CNT growth from hydrocarbon species does not require the formation of long polyyne chains [START_REF] Page | Insights into carbon nanotube and graphene formation mechanisms from molecular simulations: a review[END_REF], which is in contrast to the many simulations using hydrogen-free carbon as growth precursor [START_REF] Khalilov | Atomic scale simulation of carbon nanotube nucleation from hydrocarbon precursors[END_REF].

In the third stage, the carbon nanotube continues to grow on the catalyst surface (Fig. 6a, frames 6 and 7) due to the incorporation of the segregated C atoms, diffusing C adatoms or through direct adsorption of partially dehydrogenated C n-x H m-y fragments near the catalyst/cap interface followed by dehydrogenation. Overall, two scenarios turn out to be possible after the saturation sub-stage: (1) either C atoms continue to diffuse into the cluster, provided they completely lose their hydrogens, or (2) they diffuse over the catalyst surface, in case they keep or only partially lose their hydrogen. In either case, the C-atoms eventually incorporate into the growing network [START_REF] Khalilov | Atomic scale simulation of carbon nanotube nucleation from hydrocarbon precursors[END_REF]. Also, our simulations show that the C adsorption/incorporation rate during the CNT nucleation depends on the number of C and H atoms as well as the H/C ratio in each plasma species (Fig 6b). Specifically, if the ratio tends to zero, the number of system C atoms rises faster, i.e., the adsorption rate of C atoms increases in the following order for the precursors:

CH 4 → CH 3 → C or C 2 H 2 → C 2 H → C 2 .
In the case where the H/C ratio is (nearly) equal for two different precursors, the adsorption rate depends on the number of C atoms in the precursor (Fig. 6b, cf. C 2 H 2 and C 6 H 6 ). Indeed, the C-addition rate decreases and consequently the C adsorption/incorporation can be delayed, if the hydrogen concentration increases in the system. Due to this effect, the CNT nucleation slows down and the grown structures can additionally be damaged in the plasma-assisted CNT growth [START_REF] Kato | Growth of Single-Walled Carbon Nanotubes by Plasma CVD[END_REF][START_REF] Zhang | Selective Etching of Metallic Carbon Nanotubes by Gas-Phase Reaction[END_REF]. Note also that µm-length CNTs can catalyze some carbon nanostructures, including carbon nanowalls, when the concentration of C radicals exceeds that in the plasma environment [START_REF] Malesevic | Combined growth of carbon nanotubes and carbon nanowalls by plasma-enhanced chemical vapor deposition[END_REF][START_REF] Morjan | High growth rates and wall decoration of carbon nanotubes grown by plasma-enhanced chemical vapor deposition[END_REF].

However, no experimental observation of this phenomenon in the CNT nucleation stage has been reported.

Plasma-assisted etching of carbon nanotubes.

Our simulations show that the incipient or small carbon nanotubes during their plasma-assisted growth can be eventually removed from the catalyst surface, if the plasma contains a high concentration of H-radicals. Fig. 7a shows that the ideal nanotube (frame 1) is gradually damaged and etched (frames 2-8) before its total elimination (frame 9) upon hydrogenation. Incoming H-atoms preferable stick on the catalyst particle, as well as on the top of the tube (Fig. 7a, frame 2) due to its relatively high curvature. The latter is explained by the increased reactivity due to the increase of the pyramidalization and misalignment of the π-orbitals of the C-atoms of the tube cap [START_REF] Niyogi | Chemistry of Single-Walled carbon Nanotubes[END_REF].

According to our previous calculations, the energy barrier for H-diffusion (0.57 eV) on the Ni-catalyst surface is about three times lower than the barrier for H 2 -formation and desorption (1.63 eV) [START_REF] Somers | Temperature influence on the reactivity of plasma species on a nickel catalyst surface: An atomic scale study[END_REF][START_REF] Khalilov | How the alignment of adsorbed ortho H pairs determines the onset of selective carbon nanotube etching[END_REF]. Therefore, H-adatoms can easily move on the tube wall before they either locally associate (i.e., H-clustering effect [START_REF] Andreoni | Exohedral Hydrogen Chemisorption on a Carbon Nanotube: The Clustering Effect[END_REF]) or desorb as H 2 -molecules, according to the Langmuir-Hinshelwood (LH) recombinative desorption or the Eley-Rideal (ER) abstraction mechanisms [START_REF] Sha | Quantum studies of Eley-Rideal reactions between H atoms on a graphite surface[END_REF][START_REF] Zecho | Abstraction of D chemisorbed on graphite (0001) with gaseous H atoms[END_REF][START_REF] Neyts | Plasma-Surface Interactions in Plasma Catalysis[END_REF]. As a result, the nanotube surface is not fully covered by a "hydrogen shell" and the coverage does not reach one-third of the surface [START_REF] Khalilov | How the alignment of adsorbed ortho H pairs determines the onset of selective carbon nanotube etching[END_REF] due to simultaneous dehydrogenation. At the same time, the tube surface is locally amorphized and the tube curvature is locally increased due to the H-coalescence [START_REF] Andreoni | Exohedral Hydrogen Chemisorption on a Carbon Nanotube: The Clustering Effect[END_REF]. In this amorphized and high-curvature site of the tube, three different H-pairs, i.e., para, ortho and meta pairs can be found. Among them, in particular the ortho H pairs seriously affect the carbon network. Specifically, they significantly elongate the C-C bond length which eventually leads to the initiation of the etching process (Fig. 7a, frame 3). Recently, we found that the onset of selective CNT etching also depends on the angle between the HC-CH bond and the tube chirality vector, in addition to the tube metallicity and curvature [START_REF] Khalilov | How the alignment of adsorbed ortho H pairs determines the onset of selective carbon nanotube etching[END_REF]. After the first C-C bond breaking, the etching area propagates around this defected site due to fast H-occupation. Also, adsorbed H atoms have the tendency to move and cluster on the opposite part of the defected tube and begin the etching at that site as well (Fig. 7a,frame 4). Such etching behavior leads to the formation of parallel nanotube ribbons, which has already been experimentally observed [START_REF] Chopra | Fully collapsed carbon nanotubes[END_REF]. The defected tube finally collapses and partially covers the cluster surface (Fig. 7a, frame 5). In this tube/sheet transition, the removal rate of hexagon rings in the carbon network significantly decreases, as shown in Fig. 7b (after 0.7⋅10 8 MD(+MC) steps). Also, the number of pentagon rings considerably decreases. This clearly indicates that the obtained nanosheet becomes a nearly perfect C network due to the etching or catalyst-assisted conversion of pentagons, heptagons, and other C rings into hexagons. Such conversion eventually allows the appearance of graphene patches on the cluster [START_REF] Khalilov | Microscopic mechanisms of vertical graphene and carbon nanotube cap nucleation from hydrocarbon growth precursors[END_REF]. Due to hydrogenation or rehydrogenation, their edges are partially terminated by H atoms and thus free-standing graphene patches or graphene nanowalls can form on the catalyst surface (Fig. 7a, frame 6). The formation of such nanowalls is also in agreement with experimental evidence that no carbon nanowalls are formed without H 2 gas flow or without irradiation by H radicals [START_REF] Kondo | Nucleation and Vertical Growth of Nano Graphene Sheets, Graphene -Synthesis, Characterization, Properties and Applications[END_REF]. Subsequently, sheet/wall or wall/sheet transformations occur, due to the competition between (re)hydrogenation and dehydrogenation as well as H-etching. As mentioned in the previous section, such unstable carbon structures were also detected during in situ TEM observations before the CNT nucleation from C 2 H 2 [START_REF] Yoshida | Atomic-Scale In-situ Observation of Carbon Nanotube Growth from Solid State Iron Carbide Nanoparticles[END_REF]. The size of these unstable nanostructures continuously decreases during the hydrogenation/etching process (Fig. 7a, frame 7). Finally, we can find ordinary carbon structures, i.e., simple carbon rings and carbon polyyne chains on the catalyst surface (Fig. 7a, frame 8) before the complete removal of the CNT. Note that while increasing flow of hydrogen can speed up the degradation/etching process, it does not change the etching scenario.

During the entire hydrogenation and subsequent etching process, the number of H adatoms first increases until the etching onset and their number gradually reduces afterwards due to both the H 2 formation and the etching, while the surface C atoms are continuously removed. Finally, we observed mostly methane and ethylene molecules, as well as methyl radicals besides other etched hydrocarbon species in the gas phase, in addition to a high number of H radicals and H 2 molecules during the hydrogenation and subsequent etching process.

Conclusion

We study the nanoscale mechanisms for the Ni-catalyzed CNT growth in a hydrogencontaining plasma environment using our hybrid MD/MC simulation technique and TEMbased experiments.

We conclude that the role of hydrogen in the re-structuring of a Ni film and consequent appearance of Ni nanoclusters, as well as in the CNT nucleation/etching from plasma species, is very important and thus its careful control during aforementioned processes is essential. Specifically, the H 2 plasma can delay the sintering of the Ni powders and this leads to maintain nanoparticles longer on the semiconductor surface. Also, in CNT nucleation, the C adsorption rate is found to inversely depend on the H content in the plasma precursor.

Moreover, due to the H-effect, i.e., hydrogenation, dehydrogenation and H-etching, different C structures, including nanosheets, nanowalls, and free-standing polyyne chains can be found during both CNT nucleation and etching processes. All obtained results are in agreement with available literature data, including theoretical and experimental evidences.

In general, this study clarifies the nanoscale mechanisms in the nanocatalyst formation and subsequent carbon nanotube nucleation/etching in a plasma environment which is important for the eventual goal of understanding the selective synthesis of CNTs with a definable chirality.

1 .

 1 Ni nanocluster formation 1.1 Deposition of Ni thin film. While Ni deposition has been experimentally studied [13, 40, 51], atomic-scale simulation results can provide details at the nanoscale of intermediate processes.

Figure 1

 1 Figure 1 Snapshots of top and side views of the simulation cell during deposition of Ni atoms (green) on the Si(100) substrate (brown) at room temperature.

1. 2

 2 Re-structuring of Ni thin film. After the deposition process, the obtained Ni thin film is annealed in order to obtain Ni nanoclusters on the Si surface (see Fig. 2, 0 ps, 1000 K and 1200 K).

Figure 2

 2 Figure 2 Evolution of the restructuring of a Ni film to the Ni nanocluster (or nanocatalyst) at temperatures of 1000 K and 1200 K during the annealing.

  (72 ps, indicated by the arrow) are approximately 0.3 nm and 3.6 nm, respectively, and the break up condition is thus satisfied.

Figure 3

 3 Figure 3 Breaking of a cylindrical Ni structure (indicated by the arrow) during the Ni re-structuring process at 1000K.

Figure 4

 4 Figure 4 SEM images after the annealing of Ni film (white) with initial thickness of 6 nm on Si|SiO 2 substrate (black) at 650°C, (a) without and (b) with applying a H 2 -plasma.

Figure 5

 5 Figure 5 Top views of the interaction of two Ni 55 nanoparticles on a virtual Si surface at temperatures of 300 K and 800 K in the regimes of annealing, (a) without and (b) with H 2plasma. Here, Ni and H atoms are shown in green and blue colors, respectively.

Figure 6

 6 Figure 6 (a) Nucleation process of a CNT from C n H m species on a Ni nanocatalyst, which is physisorbed on a virtual Al or Si substrate. Here, Ni, C and H atoms are in green, red and blue colors, respectively; (b) evolution of adsorbed C atoms on/in the Ni nanocluster during the CNT nucleation from different hydrocarbon species. For the sake of clarity, the different cases are presented in different colors.

Figure 7

 7 Figure 7 (a) Hydrogenation and etching process for a base-grown (5,5) tube; and (b) number of carbon rings (hexagon, pentagon, etc) as a function of the number of MD(+MC) steps.

Acknowledgments U. K. gratefully acknowledges financial support from the Fund of Scientific Research Flanders (FWO), Belgium (Grant No. 12M1315N). The work was carried out in part using the Turing HPC infrastructure of the CalcUA core facility of the Universiteit Antwerpen, a division of the Flemish Supercomputer Centre VSC, funded by the Hercules Foundation, the Flemish Government (department EWI) and the Universiteit Antwerpen. The authors also thank Prof. A. C. T. van Duin for sharing the ReaxFF code.


Clément Chatras
email: clement.chatras@renault.com

Vincent Giard
email: vincent.giard@dauphine.fr

Mustapha Sali
email: mustapha.sali@dauphine.fr

High variety impacts on Master Production Schedule: a case study from the automotive industry

Keywords: mass customization, planning Bill of Materials, generic Bill of Materials, alternative services, master production schedule, forecast, automotive industry

In mass customization production systems, end-products diversity proceeds from the combination of hundreds of alternative components (ACs) involved in the final assembly stage. Since diversity is very important, forecasts cannot be done directly at the end-product level for production planning purposes. Thus, it is generally agreed in a mass customization context that Master Production Schedules (MPSs) have to be defined for each AC rather than each end-product. Nevertheless, some issues may hamper this approach. First, the MPSs have to be consistent with the technical and commercial constraints that prevent some ACs combinations. Second, due to longer lead times resulting from globalization, production planning must accommodate with uncertainty over the frozen horizon. To address these issues in the automotive sector, a new way to represent product through the concept of the Alternative Services (ASs), which refers to the alternative commercial features offered to the customer, has been introduced by some carmakers. This paper proposes a description of this approach and discusses its relevance regarding the issues of the classical production planning process. Some suggestions are proposed to overcome the limitations of the approach based on the ASs.

INTRODUCTION

Over the last five decades, many organizations have switched from mass production to mass customization in order to better meet customer expectations. This has resulted in an explosion of endproduct diversity that far outstrips the annual production capacity [START_REF] Pil | Linking Product Variety to Order-Fulfillment Strategies[END_REF]. The automotive industry is a perfect example of this context since a car results from the combination of thousands of alternative components. According to [START_REF] Anderson | Agile product development for mass customization: how to develop and deliver products for mass customization[END_REF], the automotive industry is one of the few industries that have pushed the envelope further in implementing mass customization. Accordingly, all the examples referred to in our paper are drawn from a simplified but real automotive case study.

End-product diversity is at the origin of production and material requirements planning issues especially for components whose lead time exceeds the frozen horizon. The classical Bill of Material (BOM) explosion mechanism that relies on end-product forecasts provided by the sales department is not manageable in a high diversity context. Indeed, due to the millions of potentially manufacturable products, it is impossible to obtain directly reliable forecasts at end-product level. This paper proposes to describe the approach used in the automotive sector for sales and production planning based on an original representation of product diversity. This approach is discussed and compared to another alternative that deals with demand uncertainty over the frozen horizon.

From a production standpoint, millions of end products may be obtained on an assembly line through the combination of alternative components and a few physical postponement operations (such as painting…). We assume that an alternative component (AC) is a variant chosen from a set of alternative components (SAC) that are used in the same assembly line workstation. Obviously, end products assembled on a given line share many components but those components that are systematic do not impact the end-product variety. For instance, a variety analysis of a Renault assembly line shows that vehicle production involves around 700 systematic components and as many SACs. The 700 SACs account for a total of more than 1,900 ACs.

Because of huge product variety, forecasting is difficult in practice either at end product level or at alternative component level. Thus, the particular context of mass customization has led automakers to devise a new approach for production planning and diversity representation. This solution hinges on functional product description that reflects the customer perception through what we call "alternative services" (ASs) rather than the organic perception directly based on components. In this new perception, a car can be completely defined by the choice of one AS in each different Set of Alternative Services (SASs) just like what customers do through the online configurators. The use of each AC on the assembly line is related to a certain combination of ASs called "predicates". The link between ASs and ACs, which remains needed for operational matters, is done through the predicates. The transition from a functional perception to an organic one is described in [START_REF] Chatras | High variety impacts on Bill of Materials Structure: Carmakers case study[END_REF].

The sales departments can only make their forecasts at this functional level, which is the only one that can be understood by customers. Based on the definition of predicates for each AC, the determination of the list of all the ACs from the list of all the ASs of a car is straightforward within the frozen horizon. Unfortunately, the calculation of AC forecasts from AS forecasts over the frozen horizon is not so easy because of the numerous AS combinations and incompatibilities. In this paper we describe and analyze two possible ways to solve the methodological problem of drawing MPS that arises from the mass customization. Starting from the AS forecasts, the first proposal aims at drawing up MPSs independently for each level 1 item of the BOM (i.e. alternative components level) while the second proposal aims at drawing up MPSs for level 0 items of the BOM (i.e. end product level). The second solution is the current solution adopted by several carmakers to deal with this issue.

In the next section (Section 2), we describe the MPS definition problem within the context of mass customization. Sections 3 and 4 are dedicated to the description and analysis of MPSs construction respectively at the AC and end-product levels. This paper ends with a conclusion (Section 5) that provides some strategic perspectives to overcome the current issues of production planning encountered in the mass customization sectors.

THE MPS DEFINITION PROBLEM IN THE CONTEXT OF MASS CUSTOMIZATION

In a mass customization context, defining the MPS for the assembly line requires a prior definition of overall production volume for each period of the planning horizon. Because this type of assembly line is characterized by fixed cycle time, the periodic production volume is directly deducted from overall production time derived from the number of shifts for the period. In order to manage teams and also part supplies, this information must be planed over the planning horizon. We will suppose this to be known for purposes of this paper.

The large diversity obtained by combining ACs from different SACs makes it difficult to define the MPS at the end-product level. This difficulty can be lifted since it is not mandatory to use the finest level of combination (i.e. end product) for procurement planning. Indeed, the diversity induced by some cheap and easily obtainable ACs pertaining to SACs used for postponement (stickers, paintings…) can be ignored in the planning process. Despite the huge remaining endproduct diversity, some companies still use an MPS at this level. This solution is discussed under Section 4.

The alternative solution that consists in defining MPSs at the AC level is easy to implement when the information (BOM coefficients and total production volume per period) used for those components is reliable. Each MPS is related to a SAC and defined through the combination of the planning BOM (Stonebraker 1986) and a periodic production volume considered as known. Beyond the frozen horizon, the planning BOM related to a SAC are considered as average usage rates related to each AC. Assuming these BOM coefficients (i.e. average usage rates) are reliable, two MPS calculation approaches beyond the frozen horizon are possible.

-In the first approach, the requirements for a period beyond the frozen horizon are calculated as the total production volume for the period multiplied by the AC BOM coefficients. This calculation is based on the assumption of a certain demand. In order to limit the propagation of disruptions in the upstream part of the supply chain due to demand uncertainty that leads to MPS rolling changes, it is important to dispose of safety stocks for those ACs. -In the second approach, AC demand beyond the frozen horizon is considered known in probability. The planning BOM coefficients for a SAC are used as a vector of a multinomial distribution in which the number of trials is the total number of cars produced for the period [START_REF] Camisullis | The information to share in upstream supply chains dedicated to mass production of customized products for allowing a decentralized management[END_REF]Giard 2010, Sali andGiard 2015). AC demand in a given period beyond the frozen horizon can be seen as a random variable following a binomial distribution and the safety stock calculation can be done by defining a replenishment level associated with a predetermined shortage risk. In this context, it can be shown that the gross requirements of any component in the BOM are random variables defined as a weighted sum of binomial random variables [START_REF] Sali | Monitoring the production of a supply chain with a revisited MRP approach, Production Planning & Control[END_REF].

The relevance of those two approaches is based on the quality of the planning BOM coefficients' forecasts. When the number of SACs is low and there are no combination constraints between ACs, one can easily define planning BOM coefficients' forecasts from their demand records using for instance the exponential smoothing technique when structural changes in demand are slow. If there are constraints that limit the free combination between ACs, the use of these approaches is possible only if the structural characteristics of demand are stable.

Because of the product diversity and the multiple combinatory constraints between ACs, many companies of the automotive sector propose a new representation of end products through a list of ASs instead of a list of ACs [START_REF] Chatras | High variety impacts on Bill of Materials Structure: Carmakers case study[END_REF][START_REF] Astesana | Constraint-Based Modeling and Exploitation of a Vehicle Range at Renault's: New Requests for the CSP Formalism[END_REF]. This representation hinges on functional product description that reflects customer perception rather than the classical organic description that reflects the engineering perception. Limited in number, the ASs allow an easier management of the combinatory constraints and thus facilitate the planning process at this level.

To move from a functional to an organic perception, which is essential for scheduling, production planning and procurement purposes, a system of predicates serves to identify a list of compatible ACs from a list of compatible ASs. A predicate is a logical expression associated to an AC that indicates the specific combinations of ASs implying its use (i.e. the AC). The ASs involved in a predicate determine the use of the associated AC. Since it is possible to move from a functional to an organic description of the product through the predicates, let us see how the planning process can be performed at the level of the ASs. The planning BOM method, initially applied at the AC level, can be used at the AS level. In fact, this method is more adapted to SASs since they are fewer than SACs and correspond to a customer and sales standpoint whereas SACs correspond to a physical perception. For production planning and procurement purposes, several scenarios illustrated in table 1 must be distinguished:

-If AC procurement lead time is shorter than the frozen horizon, no MPS needs to be drawn up beyond the frozen horizon. Accordingly, the definition of systematic MPS component is purposeless since it uses known information (i.e.: production volume and planning BOM coefficient). -If AC procurement lead time is beyond the frozen horizon, an MPS beyond the frozen horizon is required. In that case, two scenarios are possible:

-AC is determined by a single AS. In this situation, the forecast is straightforward based on the planning BOM of the AS (i.e. the periodic production associated with the average usage rate of the AS within the SAS to which it belongs). The planning BOM of the AC is the same as the planning BOM of the AS. -AC is determined by several ASs from different SASs. In this situation, the definition of the AC planning BOM coefficient is a hard problem, even if all the AS planning BOM coefficients have been forecasted. Indeed, the presence of several ASs that belong to different SASs in the logical expression of a predicate is a real issue we address below through a simple example inspired from a real case in the automotive industry.

Table 2 highlights how each alternator (i.e. A1, A2 and A3 that constitute a SAC) is defined by a combination of ASs that belong to the motorization and cooling system services. The three predicates associated to the three ACs are given under the table.

Table 1 : Typology of problems for defining level 1component MPSs Table 2 : Determination of alternators via the motorization and cooling system services

Using a standard probability notation, table 3 describes the calculation of the alternator coefficients based on the motorization and the cooling system services. The calculations using the twelve variables of the table (i.e. colored cells in the table) are given under the table grouped by a hug. As the marginal probabilities P(MOi) and P(CSj) are assumed to be known, defining the planning BOM coefficients consists in finding a solution for the following system with 9 equations and 12 variables.

We are therefore faced with an underdetermined system that accepts an infinity of solutions, which doesn't allow the calculation of the alternator planning BOM coefficients. The under-determination is systematic when an AC is determined by two or more ASs belonging to different SASs. In our example, in order to define the alternator MPS beyond the frozen horizon, one needs to add three additional constraints (i.e. three equations to the previous system) or add an optimization function. Both solutions imply a measure of arbitrary decision to solve the problem.

As mentioned above, in a mass customization context the MPS could be defined for components of level 1 in the BOM (i.e. ACs level) (Section 3) or for end products (Section 4). We identified two possible uses of planning BOM coefficients to define an MPS: a deterministic approach and a stochastic one. In the MPS definition for components of level 1, the stochastic approach beyond the frozen horizon is adopted. To our knowledge, this approach is a novelty since it is not yet used by practitioners. In the end product MPS definition, which is the method used by some carmakers, the deterministic approach is adopted.

DEFINING MPSs FOR LEVEL 1 OF THE BOM

Table 1 shows that the MPS definition at the level 1 of the BOM is simplified by the setting aside of the systematic components and those whose lead time does not exceed the frozen horizon. Thus, the problem arises for the ACs not included in these two categories.

The probabilistic approach of MPS definition beyond the frozen horizon described in the previous section can be used in the presence of AC BOM coefficient forecasts. Two cases can be distinguished:

-Case 1. AC BOM coefficients are calculated from the ASs which determine them. -Case 2. AC BOM coefficients are calculated on the basis of their time series without reference to ASs. In this case, it is useless to make forecasts for ASs that determine these ACs.

Case 1 -SAC BOM coefficients calculated from AS coefficients

The first case requires knowing BOM coefficients associated to the ASs that determine the relevant ACs through the predicates. Two possibilities have to be distinguished:

-AC is determined by a single AS. The BOM coefficient forecast for this AC is determined by those of the relevant AS. One may then use the probabilistic approach referred to under Section 2, using an arbitrarily weak risk of stockout, to determine component production launches, required to fully or partly cover random requirements. This approach has been implemented on cases drawn from real situations in the automotive industry [START_REF] Sali | Monitoring the production of a supply chain with a revisited MRP approach, Production Planning & Control[END_REF]. The stochastic approach yielded clearly better results than the deterministic approach based on end-product MPS definition. -AC is determined by several ASs belonging to different SAS. In this situation, the BOM coefficient forecast problem is more complex as shown at the end of Section 2. As mentioned, one possible solution to calculate the AC BOM coefficients consists in introducing an objective function to optimize under constraints. In the example of table 3, the values of the variables 1 P (MO CS )

j i j i j ∧ - ∧     ∑
. The left-hand part of table 4 corresponds to the "historical" table of the AC BOM coefficients (i.e. joint coefficients of combined ASs). The right-hand part provides the forecasted AS BOM coefficients (i.e. marginal distributions of the motorization and cooling system services) and the solution which minimizes the selected criterion, under constraints of respect of the margins (see equations under table 3). Note that this method is used here locally and thus differently from what will be described under Section 4. Indeed, for the definition of AC BOM coefficients within each SAC, only the BOM coefficients of the ASs involved in the relevant predicate are used. The AS BOM coefficients are independently reused as many times as necessary for the different SACs.

Table 3: planning BOM coefficient definition for alternators based on motorization and cooling system services

Case 2 -Calculation of SAC BOM coefficients based on time series

Where ACs are determined by more than two ASs, the above approach becomes unwieldy. The MPS should then be prepared at AC level and the solution consists in performing direct BOM coefficient forecasts on the ACs. Rather than using the last observed BOM coefficients values for a given AC, one may use for instance a simple exponential smoothing technique, well adapted in the case of slow changing demand structures.

AC P( )

i is thus replaced by

 AC P ( ) t t i ′> =  1
AC AC P( ) (1 ) P ( )

t t i i α α - ⋅ + -⋅
for the periods beyond t. In these calculations, the constraint

 AC P ( ) 1 t t i i ′> = ∑
guarantees that the sum of BOM coefficients is equal to 1 within the same SAC. This is guaranteed at the initialization of the calculation process of the exponential smoothing that uses BOM coefficients whose sum must be equal to 1. This condition is fulfilled since the same smoothing coefficient is used for all the ACs of the same SAC. The exponential smoothing technique is adapted to slow demand structure changes and short procurement lead times, failing which, the sales department should adopt a different forecasting technique.

Two complementary remarks can be made. The value of the smoothing coefficient can be adjusted periodically to minimize the forecasting errors. In addition, the use of a linear filter generates oscillations known as the Slutsky-Yule effect, which induce a skew in the forecasts [START_REF] Kendall | Time series[END_REF]. This skew can be easily neutralized by increasing forecasts by the quantity corresponding to the standard deviation of forecast errors. This avoids undervaluation if the forecast is on the high part of an oscillation (which can be known only latter), at the cost of building excess safety stock in the contrary case. Direct BOM coefficients forecasts using a simple exponential smoothing were subject to a real time three-month benchmark in an automotive assembly plant, on the supply of wiring harnesses which are determined by several ASs. A comparison with the system in place based on the deterministic approach described in section 4, showed a clear superiority of the probabilistic approach based on direct forecasts [START_REF] Camisullis | The information to share in upstream supply chains dedicated to mass production of customized products for allowing a decentralized management[END_REF].

DEFINING MPS FOR LEVEL 0 OF THE BOM

Some carmakers (Renault and PSA…) propose a process that gradually transforms the sales information in stages, expressed at SAS level, to draw up an MPS at end-product level. To move from SAS level forecasts to end-product forecast level, the ASs from different SASs are iteratively combined to gradually expand the definition of an object which, ultimately, becomes a finished product.

The combination of the ASs representing the different SASs gives rise to successive resolutions of equation systems that are generally under-determined as previously shown. The final result of this iterative process, called "enrichment process" by one of the manufacturers that uses this approach, is a forecast expressed in Completely Defined Vehicles (CDV). Let us briefly examine this iterative process.

At each step of the enrichment process, the resolution of the equation system consists in choosing a particular solution among the infinite number of possible solutions. To select a particular solution, one must define a selection criterion that may be for example, as proposed in the previous section, the minimization of a difference between the forecasts (i.e. decision variables) and the last observed values (i.e. historical values). To pursue the logic of calculating the coefficients at a higher level of AS combination (eg. combination involving three SASs), MO CS i j ∧ combinations such as P(MO CS ) 0

i j ∧ > are
considered as ASs that pertain to a set that can be qualified as a Meta-SAS. This iterative process therefore aims to gradually expand, at each period of the planning horizon, the sales forecasts expressed at the ASs level.

∧

are not included in the Meta-SAS combining the Motorization and Cooling system SASs).

The existence of alternative components determined by several SASs involves multiple steps in the enrichment process. The extreme case is that of an AC determined by ASs from all the SASs. In practice, this method leads to forecasts expressed at end-product level (i.e. level 0 of the BOM) excluding postponed components. Without any improvement, the described enrichment process would be difficult to implement in practice. Indeed, at the last stage of enrichment, the number of decision variables reaches the diversity of end products.

To work around this problem, ASs are combined into groups of strongly dependent SASs. The notion of interdependency between SASs refers to the existence of technical and/or commercial constraints between ASs belonging to these sets (see Chart 1). Meta-SASs formed from interdependent SASs include combinations of ASs whose coefficients are calculated using the previously illustrated optimization method. By grouping SASs into Meta-SASs, which by construction have more limited dependency between them, it becomes possible to calculate coefficients by combining two by two the elements of the Meta-SASs. The coefficients obtained are then used to generate, for each period of the planning horizon, CDVs that are characterized over all the SASs. The total number of CDVs per period is equal to the total periodic production of the assembly line.

Table 4: example of determination of alternator BOM coefficients based on motorization and cooling system services using the minimization criterion

In accordance with the time-driven technical and commercial constraints, a CDV is one of the possible ASs combinations representative of all the SASs. A CDV has no physical existence, but is a reflection of what could be a real finished product assembled at a given time. It may be considered as a virtual product. By construction, a CDV respects the commercial and the technical constraints. This property gives the CDV the same consistency as marketable products.

The ACs making up a CDV can thus easily be identified from predicates.

To illustrate how this approach works, take the example of a vehicle that is configurable from the SASs specified in Chart 1 (Commercial Region, Motorization, Equipment level, Transmission, Driver Side, Wheel rim, Cooling System, Embedded Computer). To simplify matters, let's assume a stable commercial offering over the planning period and a production line dedicated to the assembly of the relevant vehicle model.

In the approach described in this section and inspired from the process currently used by several car manufacturers, the sales department expresses its demand in the form of forecasts for each AS for all the SASs. These forecasts are expressed as a percentage of total production per period. The sales forecast enrichment process can be described in four steps. The role of each step is illustrated through the previous example.

The first step in the enrichment process is construction of several Meta-SASs as shown below. It may be noted that in this arbitrary construction, SASs belonging to a Meta-SAS are interdependent and that the Meta-SASs, taken in pairs, are independent except where one of them is that formed from the SASs: Equipment Level, Commercial Region and Motorization. This particular Meta-SAS, called "principal" and highlighted in dark grey in Chart 1, constitutes a base from which sales forecasts will be enriched by combining their elements with those of the secondary Meta-SASs highlighted in light grey in Chart 1.

The second stage of the enrichment process consists in applying the optimization method referred to above (see beginning of the section) to each Meta-SAS in order to calculate coefficients for AS combinations. This yields marginal coefficients related to new synthetic ASs, called Meta-ASs, resulting from the combination of ASs. The coefficients thus obtained for the secondary Meta-ASs from different secondary Meta-SASs are considered as independent probabilities.

The third step also uses the above optimization method to calculate the coefficients associated with the combination of two Meta-ASs, one of which belongs to the principal Meta-SAS. CR1 0,020 0,050 0,070 0,090 0,080 0,060 0,010 0,015 0,060 0,070 0,005 0,070 0,600 CR2 0,050 0,050 0,040 0,020 0,040 0,080 0,025 0,010 0,080 0,005 0,400 CRSet 0,020 0,100 0,120 0,090 0,120 0,080 0,050 0,095 0,085 0,080 0,085 0,075 1,000 CR1 0,027 0,000 0,000 0,073 0,000 0,000 0,000 0,017 0,083 0,127 0,107 0,117 0,550 CR2 0,067 0,107 0,103 0,053 0,043 0,077 0,000 0,000 0,000 0,000 0,450 4)

Chart 1 : Meta-SAS Definition

SAS

Meta-SAS

The fourth and final step of the enrichment process is to build, using the coefficients calculated in the third step, a number of CDVs that correspond to assembly line total production for the period. To this end, the coefficients calculated in the third step are transformed into conditional probabilities (probability of having a CDV with a Meta-AS belonging to a secondary Meta-SAS knowing that this CDV is characterized by a Meta-AS of the principal Meta-SAS). These conditional probabilities are used to generate a random integer number of different CDVs. The total number of generated CDVs must be equal to total assembly line production for the period. The independence of the secondary Meta-SASs serves to characterize independently a CDV for each one of them. The conditional probabilities are updated after each draw (draw without replacement) to reflect the changes in the properties of the population. The products generated and characterized for all SASs correspond to the MPS at level 0 of the BOM. The quantities generated for the MPS are considered as reliable over the non-frozen horizon. It is then possible to calculate AC requirements from these quantities through their predicates.

Because they match the physical features of CDVs, MPS defined at BOM level 0 are comfortable for upstream supply chain management purposes. Indeed, their use for medium-term planning is treated as a mere extension of the short-term schedule. By construction, the use of finished product level MPS guarantees consistent component requisitions in terms of technical and commercial constraints. While finished product level production scheduling is mandatory for the short-term, conducting this exercise beyond the frozen horizon produces the impression of having reliable information for the medium-term.

In practice, the quest for strict consistency through compliance with the technical and commercial constraints, by means of finished product level MPSs has substantial consequences on the performance of the upstream supply chain. The rules inherent in the enrichment process can generate serious disturbances that propagate along the upstream supply chain, generating a bullwhip effect [START_REF] Childerhouse | On the impact of order volatility in the European automotive sector[END_REF][START_REF] Niranjan | Measuring information distortion in real-world supply chains[END_REF].

The detailed analysis of the enrichment process used by Renault, briefly summarized in this section, revealed the existence of this effect [START_REF] Sali | Exploitation de la demande prévisionnelle pour le pilotage des flux amont d'une chaîne logistique dédiée à la production de masse de produits fortement diversifiés[END_REF]. From AS planning BOM coefficients provided by the sales department, the solution selected by the automaker to generate MPSs at end-product level uses mathematical programming techniques to solve a series of cascading optimization problems for each period of the planning horizon. By propagating errors in each iteration, this solution ultimately produced unreliable forecasts leading to shortages and costly emergency supplies, or overstock that impacted cash flow.

CONCLUSION

In the context of mass customization, characterized by multiple combinatory restrictions of the ACs, organizations represent end products through combinations of alternative services that correspond to functional features understandable by both customers and sales departments.

One must rely on SASs to prepare MPSs beyond the frozen horizon since sales department forecasts are based on a functional product description. Switching from AS-to AC-based forecasts is however not straightforward when ACs are determined by several ASs. Two MPS formulation options are available and described in our paper: one, used by several carmakers, boils down to defining MPSs at level-0 of the BOM, the other is based on level-1. Both solutions present advantages and drawbacks, which calls for further analysis and research.

Whatever the preferred solution, both increasing diversity and longer lead times in the supply chain severely impact forecasting reliability beyond the frozen horizon. This induces unwarranted supply chain costs resulting from stockout and associated remedial measures. From a strategic standpoint one should focus on the root cause of both these approaches implemented in the hope of increasing margin often without proper evaluation of their potential impact on supply chain management cost.

In order to improve their planning process sustainably in a context of strong diversity, carmakers have strategic levers on product and upstream supply chain design at their disposal. Streamlining diversity to limit the cases of determination involving a multiplicity of steps requires standardization and modularization measures at the product design stage including ACs. To address technical unfeasibility of perfect modularity involving fully standardized interfaces, the supply chain design lever may be relied on to switch back to local procurement to cut lead times beyond the frozen horizon for the ACs determined by more than one AS.

  of BOM coefficient determination for a combination of three ASs belonging respectively to the Motorization, Cooling System and Commercial Region SASs *Forecast coefficients obtained at the last iteration (see. Table

  Table 5 illustrates the calculation of coefficients associated with combinations of ASs from the SASs Motorization, Cooling system and Commercial Region (prohibited combinations MO CL

	i	j