Mixture experimental design applied to gallium-rich GaO3/2-GeO2-NaO1/2 glasses

A mixture-design method was used as a multivariate data analysis tool to investigate properties of gallate glasses in the GaO3/2-GeO2-NaO1/2 system. The objective was to design compositions having gallium oxide as the main component of the glass expected for further applications like wave guided transmission in the near-and mid-infrared optical domain. Accordingly to this tool a set of glass compositions were defined in a constrained domain of the ternary glass system rich in gallium oxide and prepared by the traditional melt-quenching technique. Sample density, refractive index, multiphonon cut-off wavelength as well as thermal characteristics were measured as representative multiple responses and a linear model fitting those physical properties vs glass composition was established. The obtained model allows to predict and optimize the different responses for any composition in the studied domain and a better understanding of the influence of each glass component on its physical properties. This study demonstrates that gallium-richest compositions lead to an increase in the polarizability and contribute to a slight enlargement of the transmission window which are the primary conditions requested for intended applications. Complementary structural investigation allows to confirm the influence of the component content in the glass composition revealed by the model.

Introduction

The demand for adapted glass compositions implemented in optical applications is continuously rising in the near and mid-infrared (up to 5 μm) where most of the molecules present a fingerprint absorption (harmonic vibrations). Silicate, phosphate and aluminate glasses have been widely investigated for NIR applications [START_REF] Shah | Preparation and studies of some thermal, mechanical and optical properties of xAl2O3(1-x) NaPO3 glass system[END_REF] but little attention has been devoted to gallate glass matrices, which are potentially excellent candidates. Some studies on binary and ternary gallium oxide containing glasses have been reported in the literature [START_REF] Lapp | Gallium oxide glasses[END_REF], [START_REF] Nishida | A linear relationship between the glass transition temperature and local distortion of calcium gallate barium gallate and calcium alluminate glasses[END_REF], [START_REF] Hannon | The effect of composition in lead gallate glasses: a structural study[END_REF], [START_REF] Dumbaugh | Heavy-metal oxide glasses[END_REF], [START_REF] Hwa | Polarized Raman scattering in lanthanum gallogermanate glasses[END_REF], [START_REF] Ilieva | Structural study on Ga2O3-TeO2 glasses[END_REF], [START_REF] Jewell | Properties of BaO-R2O3-Ga2O3-GeO2 (R = Y, Al, La, and Gd) glasses[END_REF], [START_REF] Schwarz | Thermal properties of Ga2O3-PbO-P2O5 glass system[END_REF], [START_REF] Bayya | BaO-Ga2O3-GeO2 glasses with enhanced properties[END_REF], [START_REF] Miyaji | Neutron and X-ray diffraction studies of PbO-Ga2 O3 and Bi2O3-Ga2O3 glasses[END_REF], [START_REF] Miyaji | Coordination of Ga 3 + ions in PbO-Ga2O3 glasses as determined by 71 Ga NMR[END_REF], [START_REF] Sakka | Structures of gallate, aluminate and titanate glasses[END_REF] and [START_REF] Murthy | Properties and structure of glasses int he systeme M2O-Ga2O3-GeO2 (M = Li, Na, K)[END_REF] but only a few present a gallium oxide rich content.

From a structural point of view, gallium is expected to act like aluminum in the vitreous oxide network. When the latter is introduced to an alkali-silica glass system, Si 4 + ions can be substituted by Al 3 + ions, which are then compensated by alkali ions, thus reducing the number of non-bridging oxygens (NBOs) [START_REF] Angot | A high-temperature Raman scattering study of the phase transitions in GaPO4 and in the AlPO4-GaPO4 system[END_REF] and [START_REF] Barton | Le verre: Science et Technologie[END_REF]. Regarding optical properties, gallate glasses are expected to exhibit a wider transparency window in the middle infrared compared to classical oxide glasses [START_REF] Bayya | BaO-Ga2O3-GeO2 glasses with enhanced properties[END_REF], [START_REF] Lapp | Gallium oxide glasses[END_REF], [START_REF] Higby | Properties of barium gallium germanate glasses[END_REF] and [START_REF] Angot | A high-temperature Raman scattering study of the phase transitions in GaPO4 and in the AlPO4-GaPO4 system[END_REF], such as borates, phosphates or silicates, because of the higher atomic weight of gallium. Also they exhibit significant polarizability and hyperpolarizability which makes them promising materials for nonlinear optical properties [START_REF] Lapp | Gallium oxide glasses[END_REF]. Their low phonon energy (~ 850 cm -1 for germanogallate glasses, ~ 675 cm -1 for some alkali-earth-gallate glasses [START_REF] Sakka | Structures of gallate, aluminate and titanate glasses[END_REF] and [START_REF] Fukumi | Raman spectroscopic study of the structural role of alkaline earth ions in alkaline earth gallate glasses[END_REF]) makes them of interest for hosting rare-earth ions in luminescent and laser materials [START_REF] Wen | Tm 3 + doped barium gallo-germanate glass single-mode fibers for 2.0 μm laser[END_REF]. From the foregoing gallate-based glasses would offer an excellent compromise for mid-infrared transmission applications with the required thermal stability and mechanical resistance for fiber waveguide shaping [START_REF] Fan | Structure and physical properties of Bi2O3-B2O3-Ga2O3 glasses[END_REF]. Recently, studies on optical waveguiding with gallophosphate glasses have been reported by Vangheluwe et al. [START_REF] Vangheluwe | Enhancement of nanograting formation assisted by silver ions in a sodium gallophosphate glass[END_REF] but none with compositions rich in gallium oxide. Most of the previous investigations conducted on germanogallate glass systems, which are expected to fulfill these properties, have been focused on the germanium oxide (GeO2) rich portion of the glass system [START_REF] Dumbaugh | Heavy-metal oxide glasses[END_REF], [START_REF] Jewell | Properties of BaO-R2O3-Ga2O3-GeO2 (R = Y, Al, La, and Gd) glasses[END_REF], [START_REF] Bayya | BaO-Ga2O3-GeO2 glasses with enhanced properties[END_REF], [START_REF] Higby | Properties of barium gallium germanate glasses[END_REF], [START_REF] Wen | Tm 3 + doped barium gallo-germanate glass single-mode fibers for 2.0 μm laser[END_REF], [START_REF] Jewell | Optical properties of BaO-Ga2O3-GeO2 glasses for fiber and bulk optical applications[END_REF], [START_REF] Zhao | Effect of alkali metal oxides R2O (R = Na, K) on 1.53 μm luminescence of Er 3 + -doped Ga2O3-GeO2 glasses for optical amplification[END_REF], [START_REF] Padlyak | Nature of intrinsic luminescence in the glasses of CaO-Ga2O3-GeO2 system[END_REF], [START_REF] Padlyak | Optical spectroscopy of the Er-doped glasses with 3CaO-Ga2O3-3GeO2 composition[END_REF], [START_REF] Bufetov | Bi-doped fiber lasers[END_REF], [START_REF] Cao | A new fluorogermanate glass[END_REF] and [START_REF] Fan | Dependence of fluorescence properties on substitution of BaF2 for BaO in barium gallo-germanate glass[END_REF].

On the contrary this study will focus on the investigation of gallium-rich glass compositions belonging to the sodogermanogallate ternary system with the perspective of applications in the mid-infrared range (3 to 5 μm). For this purpose, a first modelling of the glass properties (thermal stability toward crystallization, IR transmission, linear refractive indices and density) is explored with the mixture design technique in the gallium oxide rich system. In our knowledge, this represent the first study in the literature with this system which does not permit an effective use of properties calculations features from a glass property database as the SciGlass information system [30]. Properties linear modelling of the measured properties is expected to predict the influence of each component in the studied domain. Complementary structural investigations will be conducted using Raman and 71 Ga NMR spectroscopy to confirm, precise and support the model predictions.

Experimental details 2.1 Glass preparation

Here, the used components for the ternary composition studied system are GaO3/2, GeO2 and NaO1/2 which are more relevant to highlight the role of each cation. Glasses in the GaO3/2-GeO2-NaO1/2 system were then prepared by traditional melt-casting technique from gallium oxide Ga2O3 (99.99%), germanium oxide GeO2 (99.99%) and sodium carbonate Na2CO3 (99.95%). The weighed and mixed powders were pre-sintered at 950 °C for 1 h and then melted in platinum crucible in ambient atmosphere for 30 min at 1300-1550 °C depending on the compositions. The melted glasses were poured on a stainless steel mold preheated at 40 °C below the glass transition temperature and annealed for 6 h at the same temperature before being slowly cooled to room temperature. All the obtained glasses were colorless, transparent and bubble-free. The glass samples were finally cut and polished on both parallel faces for optical characterization.

Characterization techniques

Physico-chemical and optical response measurements

Two physico-chemical experimental responses, density and thermal stability, are studied here. The density of the glass, ρ, is related to composition and structural evolution, as well as the thermal stability against crystallization, corresponding to the difference, ΔT, between the onset temperature of crystallization and the glass transition temperature. The latter must be as large as possible (≥ 100-120 °C) to facilitate any future technological transfer. Two optical parameters are also studied: the IR cut-off wavelength, which is crucial for mid-infrared applications, and the refractive index, n, which is expected to be the signature of high third-order optical susceptibilities which can be involved in higher level applications like Raman gain or frequency mixing [START_REF] Rivero | Quantifying Raman gain coefficients in tellurite glasses[END_REF] and [START_REF] Boling | Empirical relationships for predicting nonlinear refractive index changes in optical solids[END_REF].

Confirmation of nominal glass stoichiometry was verified with a deviation of ± 2 mol% by using Electron Probe Micro Analysis on a VG Microlab 310F apparatus. The glass transition temperature, Tg, and the onset temperature of crystallization, Tx, were determined by differential scanning calorimetric (DSC), using a Netzsch DSC Pegasus 404F3 apparatus on glass pieces in Pt pans at a heating rate of 10 °C/min and with a precision of ± 2 °C. The glass density was determined at room temperature by the Archimedes' method with an Alfa-Mirage MD-300S densimeter using deionized water as immersion liquid. The measurement precision was estimated to be ± 0.01 g/cm 3 .

The infrared transmission spectra were recorded on glass samples of about 2 mm thickness on a FTIR Perkin Elmer Frontier spectrometer. The linear refractive index was measured by employing the M-lines prism coupling technique (using a Metricon 2010) at 532 nm with an accuracy of ± 0.005.

Complementary structural characterization

The Raman spectra were obtained using a Horiba Jobin-Yvon LabRAM HR800 micro-Raman spectrometer equipped with a 532 nm laser. 71 Ga magic-angle spinning NMR spectra were recorded in a magnetic field of 21.1 T using a Bruker Avance II 900 MHz spectrometer at the Canadian National Ultrahigh-field NMR Facility for Solids (Ottawa, Ontario). Samples were ground to a fine powder, packed into 1.3 mm (outer diameter) rotors and spun at 60 kHz. Single-pulse excitation with a 15 degree tip angle was used to acquire 8000 transients, separated by an optimized relaxation delay of 5 s. The chemical shift reference was 1 M Ga2(SO4)3 (aq) at 0ppm. Spectra were fit using a Czjzek model in the software package dmfit [START_REF] Massiot | Modelling one-and twodimensional solid-state NMR spectra[END_REF], which allows parameter distributions suitable for network-forming oxide glasses [START_REF] D'espinose De Lacaillerie | MAS NMR spectra of quadrupolar nuclei in disordered solids: the Czjzek model[END_REF].

3. Design of experiments: glass composition study in the GaO3/2-GeO2-NaO1/2 system

Definition of the experimental domain

The realization of an optimized design of experiments provides a maximum of information with a minimum of experiments. Within the context of the selection of compositions belonging to the ternary GaO3/2-GeO2-NaO1/2 system, designing experiments consists in polynomial modelling of the properties of interest defined in terms of responses of the system as a function of the composition. Here, the cationic molar percentages of GaO3/2, GeO2 and NaO1/2 are noted xGa; xGe; xNa respectively. In a first step, a preliminary exploratory experiment was performed using the selected preparation method to reach a first piece of information and permit the definition of the domain of study. In this case, the central point in the ternary diagram (domain of compositions) where xGa = xGe = xNa = 33.3% was prepared and characterized (composition N°2 in Fig. 1, Table 1 andTable 2). In order to define our experimental domain, the range of variations of the xGa, xGe and xNa compositions and the relational constraints between variables were fixed [START_REF] Fleury | Development of an experimental design to investigate the effects of R7T7 glass composition on the residual rate of alteration[END_REF]. The sum of xGa, xGe and xNa, i.e. of all the oxides percentages, must equal to 100%. The gallium oxide GaO3/2 content was selected to be larger than that of germanium oxide GeO2, to focus on rich-gallate glasses and investigate the role played by the gallium cation as a main constituent of the glass network. This leads to the following constraints equations:

Determination and quality of the design of experiments

The optimal design of glass compositions was obtained by using the 'optimal-D' criterion of optimization from the "Design Expert®" version 8 software [START_REF] Goupy | Plans d'expériences: les mélanges[END_REF]. Eight experiments, i.e. glass compositions, were selected, including the exploratory point. Adjustments were necessary for some compositions which proved to be outside the glass-forming region after preparation, using the conditions described above. In order to construct the model, compositions rich in information and showing glass-forming ability were therefore provided by the software. The final compositions are given in Table 1 and are represented in the ternary diagram shown in Fig. 1. Note that we are considering cationic relative concentrations in the present study to facilitate our investigation of the role played by gallium as a main glass constituent.

After measurements of the four selected physico-chemical and optical responses for the eight designed experiments, each response can be predicted in the domain of study by the use of a linear model with respect to wGa, wGe, wNa representing the reduced molar fractions in the reduced ternary diagram shown on Fig. 1, as expressed in Eq. ( 1) [START_REF] Rivero | Quantifying Raman gain coefficients in tellurite glasses[END_REF].

With wGa, wGe, wNa comprised between 0 and 1 and wGa + wGe + wNa = 1.

Reduced variables being comprised between 0 and 1 allow further direct comparison between a, b and c coefficients showing the estimated response for each corner of the reduced triangle. For example, if wGa = 1 and wGe = wNa = 0, for point A of the reduced ternary diagram in the direction of the richest compositions in gallium oxide, then the response is the coefficient a. The model coefficients a, b and c have then comparable values and they are estimated with their uncertainty depending on the remaining degrees of freedom for a confidence factor of 95%, considering a normal distribution of the response measurement [START_REF] Rivero | Quantifying Raman gain coefficients in tellurite glasses[END_REF]. The coefficient's estimated values give directly the possibility to evaluate the relative weights of the compositions on the studied response. They can be considered to be significant if the relative uncertainties for all estimations remain below 10%. The variance inflation factors are calculated in order to quantify the severity of multicollinearity in the least squares regression analysis leading to the coefficient's estimations [START_REF] Rivero | Quantifying Raman gain coefficients in tellurite glasses[END_REF]. For a good quality regression, these factors would not exceed 5 [START_REF] Rivero | Quantifying Raman gain coefficients in tellurite glasses[END_REF].

Results

Physico-chemical characterization

The measured glass transition temperature Tg, the onset temperature of crystallization Tx, the difference between Tx and Tg (ΔT = Tx -Tc) and the density ρ, are listed in Table 2. The thermal stability against crystallization corresponding to a large ΔT is considered as an important parameter for various practical applications [START_REF] Wen | Tm 3 + doped barium gallo-germanate glass single-mode fibers for 2.0 μm laser[END_REF]. The multiphonon absorption edge λIR, calculated here for a linear absorption coefficient of 10 cm -1 , as well as the refractive index n (measured at 532 nm), are presented in Table 2.

Modelling of physico-chemical properties

The coefficients a, b, and c, calculated by the software, for the density ρ, the thermal stability ΔT, the multiphonon cutoff wavelength λIR and the refractive index n are presented in Table 3 with their respective uncertainties and variance inflation factors. The coefficients, determined for the refractive index, the cut-off wavelength and the density, show excellent agreement with the experimental data (uncertainty below 1%). On the contrary, the coefficients estimated for the thermal stability (uncertainty higher than 12%) can be considered as less significant. All the inflation factors remain below 5, ensuring independent calculations of the coefficients estimations. Those reliable estimated values show equivalent weight of the three oxides normalized relative concentrations on the studied responses.

To predict the glass properties directly on the ternary diagram of composition the four linear models deduced from estimations in Table 3 are used coming back to the linear model given in Eq. ( 2) with real compositions xGa, xGe and xNa.

The dependence of the properties on the components can be mapped by contour line in Fig. 2. In all cases, in the gallium-rich compositions the density and the refractive index increase and the transmission window in the mid-infrared enlarges. For the density, the sodium oxide concentration is clearly unfavorable, the higher density requiring the lower NaO1/2 content as shown in Fig. 2(a). One can observe for a high xGa/xGe ratio the highest density indicating that the gallium-rich compositions have a more compact glass network.

Regarding the refractive index, one can clearly observe in Fig. 2(b) that the polarizability increase is mainly driven by the gallium oxide content, which is in agreement with the data collected by Murthy et al. [START_REF] Murthy | Properties and structure of glasses int he systeme M2O-Ga2O3-GeO2 (M = Li, Na, K)[END_REF]. Those observations are in accordance with a larger polarizability of the gallium oxide as compared to the germanium oxide [START_REF] He | First-principles study of the structural, electronic, and optical properties of Ga2O3 in its monoclinic and hexagonal phases[END_REF] and [START_REF] Weber | Handbook of Optical Materials[END_REF]. The increase of gallium oxide, to the detriment of germanium oxide content is the most important parameter for expanding the IR window (Fig. 2(c)). One can observe that the increase of gallium oxide tends to decrease the ΔT as shown in Fig. 2(d).

Then except for thermal properties, the ratio of gallium oxide appears to be the factor with the most important weight in the increase of density and optical properties. Nevertheless the ratio of germanium oxide and sodium oxide play a less important and different role when comparing the density and refractive indices on one side and the cut-off wavelength on another side. Complementary structural characterizations are expected to bring information to explain this interesting observation deduced from the mixture-design model.

Complementary structural characterizations

The Raman spectra recorded on four glass compositions with various xNa/xGa ratios and chosen for location on different edges of the experimental domain of study (= 1 for N°4: Ga38Ge25Na37 and N°2: Ga33Ge23Na33; and = 0.71 for N°3: Ga42Ge28Na30) are presented in Fig. 3.

The Raman bands are assigned according to McKeown et al. [START_REF] Mckeown | Raman spectroscopic studies of BaO-Ga2O3-GeO2 glasses[END_REF]. The band at 300 cm -1 is related to network-modifying cation vibrations in large interstitial sites. The principal band at about 500-550 cm -1 is attributed to vibrations of X O X (X = Ga or Ge) bonds. These bonds are in a polymerized network built with GeO4 and GaO4 tetrahedra, which are essentially located in rings made of 5 or 6 tetrahedra. The last two bands located at about 750 cm -1 and 870 cm -1 are respectively attributed to GeO4 and GaO4 tetrahedra with two non-bridging bonds (chain-like Q 2 units) and with one nonbridging bond (sheet-like Q 3 units). One can note that both bands are slightly shifted toward lower frequencies for N°4: Ga38Ge25Na37.

71

Ga MAS NMR experiments recorded for three selected samples with different xNa/xGa ratios are presented in Fig. 4. A single broad asymmetric peak around 180 ppm is observed for all samples, which is typical for quadrupolar nuclei in glassy environments [START_REF] Kroeker | NMR spectroscopy of glasses[END_REF]. The charge-balanced spectra (xNa/xGa = 1) can be modeled using a single site with an isotropic chemical shift of 207 ppm, quadrupole coupling constant of 11 MHz and modest parameter distributions typical of fourcoordinate gallium in oxide glasses [START_REF] Massiot | 71 Ga NMR of reference Ga IV, Ga V, and Ga VI compounds by MAS and QPASS, extension of gallium/aluminum NMR parameter correlation[END_REF]. With a reduction of sodium content with respect to gallium, additional spectral intensity can be detected on the low-frequency side of the main peak (ca. 100 ppm) which cannot be modeled using a single tetrahedral Ga site, even with more generous parameter distributions. Fig. 4 shows a single-site fit to sample 1, illustrating that additional coordination environments are required to satisfactorily fit the data. Although there is no clear spectral resolution in that region, it may be estimated by difference that approximately 15% of the Ga is found in highcoordinate (5-or 6-fold) sites.

Discussion

From the proposed models, it can be deduced that increasing the gallium oxide content produces an increase in density, linear refractive index and multiphonon cut-off wavelength. On the contrary, the thermal stability is reduced with the addition of gallium oxide to the matrix. However, the thermal stability of the prepared glasses remains excellent, i.e. larger than 100 °C, over the whole compositional region studied.

The role of the glass network structure on the physical properties has been already investigated from Raman and NMR spectroscopies [START_REF] Mckeown | Raman spectroscopic studies of BaO-Ga2O3-GeO2 glasses[END_REF] in the gallium germanate glass system to precise the effects observed from the mixture-design model on these physical properties. The bands observed above 600 cm -1 in the Raman spectra show similar depolymerization of the glass network due to the introduction of sodium oxide content comparable to that of germanium oxide (N°2: Ga33Ge33Na33 compared to N°3: Ga42Ge28Na30). For the lowest quantity of germanium oxide and an xNa/xGa ratio around 1 (N°4: Ga38Ge25Na37), one can see from Fig. 3 a shift in wavenumber of the band around 870 cm -1 attributed to Q 3 units. Henderson et al. [START_REF] Henderson | Q speciation in alkali germanate glasses[END_REF] have previously attributed in alkali germanate systems the vibration at 835 cm -1 to modified Q 3 species. In the Na2O-GeO2 system such a shift appears for the ratio xNa/xGe > 1. This was attributed to the formation of sodium digermanate-like structures. In our case it would correspond to sodium ions acting only as modifiers for the germanium. Such a hypothesis requires a compensation mechanism such as interconnected (mixed corner/edge) gallium oxide entities.

The 71 Ga NMR spectra indicate that nearly all of the gallium in these glasses is present in anionic four-fold coordination, as expected to balance the sodium cationic charge. Only when there is insufficient sodium to balance these tetrahedral gallium units is there spectral evidence for high-coordinated gallium. From these observations, the gallium environment can be tentatively described. When xNa/xGa = 1, for xGe ≤ xGa, each gallium ion occupies an anionic tetrahedral site, the charge being essentially balanced by one sodium cation as there is the same amount of gallium and sodium ions. In such a configuration, the role played by NBOs to compensate the charge of Na + cations is negligible. This insertion mode reaches its limit when xGa = xGe = xNa (glass composition N°2: Ga33Ge33Na33). Indeed, whereas the Na + cation charge is already compensated in the vitreous network, additional incorporation of gallium oxide substituted to germanium oxide would preferentially result in the formation of gallium cations with higher (5 or 6) coordination numbers, necessitating new connections between the rings made of GeO4 and GaO4 tetrahedra [START_REF] Mckeown | Raman spectroscopic studies of BaO-Ga2O3-GeO2 glasses[END_REF]. The formation of such 5 or 6-fold gallium ions can be directly related to the so-called gallate glass, with xGa > xGe and xGa > xNa (as for glass composition N°1: Ga44Ge35Na31). The shoulder observed in the NMR spectra attributed to high coordinated sites in Fig. 4 is effectively observed for the glass composition N°1: Ga44Ge35Na31. The addition of high coordinated gallium ions interconnecting the glass network is expected to favor an increase of the density and consequently the refractive index [START_REF] Higby | Properties of barium gallium germanate glasses[END_REF]. A more comprehensive NMR characterization will be necessary to prove the existence of high coordinated gallium ions for different compositions where xNa/xGa < 1 and to accurately characterize the evolution of the glass-network connectivity vs the cationic ratio xNa/xGa. On the contrary the introduction of excess sodium oxide substituted to the germanium oxide from the central point N°2 is expected to break connections in the network leading to new Q 3 species, as observed in the Raman spectra in Fig. 3. The transformation is expected when xGa = xNa for the glass N°4 (Ga38Ge25Na37), for which a shift to lower energy of the vibrational bands is observed. For this composition, no shoulder attributed to high coordinated sites can be seen in the NMR spectra, as shown in Fig. 4.

In the studied compositions the formation of a gallate network with the formation of high-coordinated sites due to an increase of GaO3/2 content with a decrease of the ratio xNa/xGa can explain the direction toward high refractive index and density expected from the model (Fig. 2). Regarding the infrared cut-off, Fig. 2 indicates that the polymerization is not a significant parameter, the main parameter remains the lowering of the xGe/xNa ratio. The evolution of the multiphonon cut-off wavelength observed from the model could be related to the shift to lower energy of the vibrational bands (Fig. 3) and the appearance of a new Q 3 unit. One must be cautious in drawing definitive conclusions since the multiphonon edge can correspond to diverse intrinsic vibration combinations.

Conclusion

A mixture experimental design, using the optimal-D based design from the "Design Expert®" software, was implemented on glasses from the GaO3/2-GeO2-NaO1/2 ternary system in order to map the evolution of physical properties (thermal stability and density) and optical properties (refractive index and multiphonon cut-off wavelength). A linear model relating each property with the glass composition has been established and represented in a contour line. Correlation using Raman and NMR spectroscopies indicates that replacing germanium oxide by gallium oxide tends to improve the network polymerization through the creation of bridging Ge O Ga bonds and probably leads to the insertion into the network of high coordinated gallium ions. The increase in network connectivity results in an increase of the density and the refractive index while maintaining a good thermal stability against crystallization. The multiphonon cut-off wavelength also increases with the xGa/xGe ratio but for comparable gallium and sodium contents.

  xGa + xGe + xNa  100 mol% xGa≥xGe The range of variation for each parameter xGa, xGe and xNa was fixed around 10 mol% to stay in the homogeneous glassforming region and to build up a linear model of the experimental responses: GaO3/2(mol%) : 33≤xGa≤44 GeO2(mol%) : 24≤xGe≤34 NaO1/2(mol%) : 29≤xNa≤38 These ranges of variations were chosen so as to inscribe the domain of study inside a reduced equilateral triangle moved from the central preliminary test to the richest part of gallium oxide in the ternary diagram. The reduced triangle and domain of composition are depicted in the diagram presented in Fig. 1.

Fig. 1 .

 1 Fig. 1. Ternary GaO3/2-GeO2-NaO1/2 diagram (cationic mol%) representing the determined experimental domain of study (in orange), the vitreous domain determined by Murthy et al. [14] (in blue) and the 8 glass compositions prepared and characterized to build up the model.

Fig. 2 .

 2 Fig. 2. Contour line representation of the properties in the studied domain. (a) Density, ρ -(b) refractive index, n -(c) thermal stability, ΔT -(d) IR cut-off, λc.

Fig. 3 .

 3 Fig. 3. (a) Raman spectra for GaGeNa glass system with different xNa/xGa ratios -(b) Q3 band normalized Raman spectra for the GaGeNa glass system.

Fig. 4 .

 4 Fig. 4. 71 Ga MAS NMR spectra of sodium gallogermanate glasses recorded at 21.1 T with a sample rotation speed of 60 kHz. Red lines represent the best-fit spectrum using a single site (see text).

Acknowledgements

This research was supported by the Canadian Excellence Research Chair program (CERC) in Photonics Innovations and the Aquitaine Region. The authors are also grateful to the Fonds de Recherche Québecois sur la Nature et les Technologies (FRQNT) (097991) and the Canadian Foundation for Innovation (CFI) (29191 and 32532) for the financial support. Access to the 21.1 T NMR spectrometer was provided by the National Ultrahigh-Field NMR Facility for Solids (Ottawa, Canada), a national research facility funded by a consortium of Canadian universities, supported by the National Research Council of Canada and Bruker Biospin, and managed by the University of Ottawa. We are grateful to Dr Victor Terskikh, Dr. John Wren and Mr. Kirill Levin for valuable assistance with the NMR experiments. Mobilities were supported by grants of the association of Campus France and Mitacs Globalink as well as grants of the Agence Nationale de la Recherche (ANR) with the program "Investissement d'avenir" number ANR-10-IDEX-03-02.


Claude Yugma
email: claude.yugma@emse.fr

Christian Artigues

Stéphane Dauzère Pérès

Alexandre Derreumaux

Olivier Sibille

A Batch Optimization Solver for diffusion area scheduling in semiconductor manufacturing

Keywords: Batch scheduling, disjunctive graph, local search, simulated annealing, wafer fabrication

This paper presents a method and a software for solving a batching and scheduling problem in the diffusion area of a semiconductor plant, the ATMEL fabrication unit in Rousset, France. The diffusion area is one of the most complex area in the fab. A significant number of lots has to be processed while satisfying complex equipment process and line management constraints. The purpose of this study is to investigate approaches to group lots in batches, to assign the batches on the equipment and to sequence these batches. Three indicators are used to evaluate the quality of a solution: the total number of moves, the batching coefficient and the X-factor. The problem is modeled through the disjunctive graph formulation. A constructive algorithm is proposed and improvement procedures based on iterative sampling and simulated annealing are developed. Computational experiments, carried out on actual industrial problem instances, show the ability of iterative sampling to improve significantly the initial solution. The proposed simulated annealing method brings in turn important enhancements to the results of iterative sampling. The software based on these methods, named Batch Optimization Solver (BOS), is currently used in the diffusion area of ATMEL. The disjunctive graph model allows in addition a high level of interactivity with the decision makers.

Introduction

Semiconductor wafer fabrication can be described as a multistage process with re-entrant flows. The processing is done layer by layer. Each layer requires several steps of processing such as chemical-mechanical polishing, diffusion, film deposition, photolithography, implant (doping) and etching. For each of the product types, and depending on the technology, a wafer goes through more than 400 process steps over a period of a few weeks. Wafer fabrication planning and scheduling is a complex task due to the large number of products and machines involved. It is further complicated by additional constraints such as re-entrant flow of operations [START_REF] Kumar | Scheduling semiconductor manufacturing plants[END_REF], setup issues, preventive maintenances, random machine breakdowns, etc.

The importance of scheduling on the performance of semiconductor wafer fabrication facilities is known for many years [START_REF] Wein | Scheduling semiconductor wafer fabrication[END_REF][START_REF] Varadarajan | A survey of dispatching rules for operational control in wafer fabrication[END_REF]. Recently [START_REF] Toba | Dynamic load balancing among multiple fabrication lines through estimation of minimum inter-operation time[END_REF] have also shown the importance of using predictive scheduling in load balancing methods among multiple semiconductor wafer fabrication lines in order to make a better estimation of the lot cycle times. The cycle time is the total time necessary to manufacture a lot of wafers. It measures the duration time spent since the lot enters the fab and the moment it leaves it. According to experimental studies, decreasing the cycle time results in an increase of the productivity in a fab [START_REF] Jr | Competing against time: how time-based competition is reshapping global markets[END_REF].

The complexity of the fabrication process is such that most previous existing work on scheduling in semiconductor manufacturing focus on dispatching rules [START_REF] Wein | Scheduling semiconductor wafer fabrication[END_REF][START_REF] Dabbas | A new scheduling approach using combined dispatching criteria in wafer fabs[END_REF][START_REF] Varadarajan | A survey of dispatching rules for operational control in wafer fabrication[END_REF]. Other approaches deal with Integer Linear Programming formulations of the problem and apply lagrangean relaxation of the equipment limitation constraints [START_REF] Liao | Daily scheduling for r&d semiconductor fabrication[END_REF][START_REF] Hwang | Design of a lagrangian relaxation-based hierarchical production scheduling environment for semiconductor wafer fabrication[END_REF]. For a given set of lagrangean multipliers, the relaxation can be solved with minimal cost network flow methods. Subgradient optimization is used to solve the dual linear program. A feasible solution is obtained by adjusting heuristically the solution of the relaxation.

In this paper, instead of scheduling the entire fab, we propose to focus on a bottleneck part of the manufacturing process. Among the complex operations involved in the fabrication of a wafer, the diffusion phase is of critical importance since the batching decisions that are involved may affect the performance of the entire wafer fab [ICN + 03, MH03]. This is also the case in the ATMEL fabrication unit of Rousset (France). The diffusion phase is used primarily to alter the type and level of conductivity of semiconductor materials. It is used to form bases, emitters and resistors in bipolar devices, and also to dope polysilicon layers. The Batch Optimization Software (BOS) project aims at developing a method and a software for solving efficiently the batching and scheduling problem of the diffusion area while taking into account numerous constraints and optimizing various criteria. Another objective of the BOS is to provide an easy way to build, modify and use the results of the software. As for the global problem, batching and scheduling in diffusion areas is mainly performed in the literature by dispatching rules. In [ICS + 03] and [START_REF] Monch | Simulation-based assessment of batching heuristics in semiconductor manufacturing[END_REF], the authors propose to use a dispatching rule to batch the operations and a genetic algorithm to assign the batches to equipment. In [START_REF] Mason | A modified shifting bottleneck heuristic for minimizing total weighted tardiness in complex job shops[END_REF], a scheduling problem in semiconductor manufacturing related to the one tackled in this paper is solved by a modified shifting bottleneck heuristic [START_REF] Adams | The shifting bottleneck procedure for job shop scheduling[END_REF].

In this paper, we propose to model solutions of the considered batching and scheduling problem through a variant of the disjunctive graph proposed in [START_REF] Mason | A modified shifting bottleneck heuristic for minimizing total weighted tardiness in complex job shops[END_REF]. A constructive heuristic is designed that builds an initial solution as well as improvement procedures based on local search and simulated annealing.

This paper is organized as follows. In Section 2, we give the problem description. In Section 3, we present the proposed disjunctive graph representation. The method for computing an initial solution and improvement procedures based on iterative sampling and simulated annealing are described in Section 4. Experimental results are given on real problem instances issued from the Rousset ATMEL fabrication unit in Section 5. The BOS software is briefly presented in Section 6. The Section 7 concludes the paper with recommendations for further research.

Problem description

The diffusion area defines a batching and scheduling problem of wafer lots on two types of equipment: furnaces and cleaning machines. These resources are able to perform several lots simultaneously. Each lot requires one or more consecutive operations on the equipments and each operation has a recipe which determines its duration and the set of machines that are able to process it. On a 24-hour basis, each operation has to be assigned to an equipment, and included into a batch, i.e. a set of operations of the same recipe that are processed simultaneously by the equipment. The constraints in the diffusion area are divided into 3 types: equipment constraints, process constraints and line management constraints.

Constraints

Some of these constraints are common to the two types of resources while others are dedicated to furnaces.

Equipment constraints

• Common constraints

Dedicated equipment: any equipment is able to process a limited set of recipes.

Maximum batch size: any equipment defines a maximum batch size corresponding to the capacity of that equipment.

Loading and unloading times: a time may be needed to load and unload a batch on the equipment.

Unavailability periods: the equipment may be unavailable during some periods (defined by time windows) due to qualification, repair, maintenance etc.

In process jobs: the equipment may be occupied at the beginning of the time horizon by in-process operations that have to be completed before the equipment becomes available.

• Specific furnace constraints

Minimum time between two batches on an furnace: some furnaces must be inspected at the end of each processed batch.

Process constraints

• Common constraints Precedence: lot operations must be performed following the manufacturing process of the lot. The operations of a lot are chained and no operation can start before the end of its predecessor, except for the first operation of each lot. Minimum time lag: there is a fixed handling and transport time between every two successive operations of a lot.

Release dates: they correspond to the arrival times of lots at the cleaning machines and furnaces. It means that a lot cannot be scheduled before its release date. Because the diffusion area is a stage of the complex global production process, the release date are precomputed by a simulation package designed by ATMEL. Fixed recipe: each operation of a lot is associated with a recipe, i.e. the lot should be processed on resources that are qualified for the corresponding recipe. This implies that all lots in the same batch must be processed with the same recipe. Process time: the process time of a batch on an equipment depends on the recipe. Maximum time lag: a time limit is given for two successive operations x and y of a lot. The difference between the starting time of y and the completion time of x cannot exceed this limit. Maximum time lags depend on the operations x and y. They can correspond to the maximum time limit between cleaning and furnace equipment.

Line management constraints

• Specific furnace constraints Minimum distance between two batches of same recipe: there is a minimum time between the beginning of two batches of the same recipe on two different furnaces.

Objectives

In semiconductor factories, several indicators are used to measure the performance of a fab. We describe thereafter the indicators that are relevant in our study. The interested reader can refer to [START_REF] Montoya-Torres | Manufacturing performance evaluation in wafer semiconductor factories[END_REF] for more details.

• The number of moves (to maximize): It corresponds to the number of completed operations on the planning horizon, which can be compared to the target number fixed by the production managers.

• Batching coefficient (to maximize): Defined on the planning horizon, it is calculated as the number of moves divided by the sum of the number of batches performed on each machine, times the maximum capacity of that machine. Note that the denominator is the number of lots that could be performed if the equipment was loaded up to their maximum capacities.

• The X-factor (to minimize): This indicator is used to evaluate the waiting times of lots in the diffusion area in order to reduce the cycle times. For a given lot, this factor is calculated as the total staying time of the lot in the diffusion area (i.e. waiting time from its arrival in the diffusion area and processing time) over its processing time.

These three indicators will be used to evaluate solutions in the fab.

Others indicators not considered in this paper are of interest. Among them, the Work-In-Process (WIP), is defined as the number of wafers being in the fab at a given period, either in a production state or in a non-production state (e.g. transport and waiting). The law of Little [START_REF] Little | A proof for the queuing formula l = λw[END_REF] establishes that if a system is stable and stationary then the average WIP is proportional to the average cycle time. The throughput is defined as the outgoing number of wafers of the fab per unit of time. The evolution of the throughput in time makes it possible to know if the system is stable i.e. to know if there is no accumulation of product in the fab or if the system evolves/moves according to forecasts. Authors, such as [START_REF] Glassey | Closed-loop job release control for vlsi circuit manufacturing[END_REF], observed that there is a relation between the increasing of the throughput and the output of a fab. The cycle time drastically increases when throughput is close to the maximal capacity of the fab. Hence for both WIP and throughput indicators, the X-factor indicator is adequate.

The goal is to optimize the various performance measures, while taking into account the numerous complex constraints. Another objective is to make the proposed schedule easily modifiable by the decision makers at the fabrication level. Machines may be declared as down, operations may be moved, lots may be removed or inserted. Such requirements can be achieved thanks to the flexibility of the disjunctive graph representation. The next subsection described the mathematical formulation used to model the problem.

Problem Formulation

The scheduling problem can be formulated as follows. For sake of clarity, the above-described Unavailability periods, In-process jobs and Minimum distance between batches of the same recipe constraints are not included. In Section 3, we describe how we tackle these additional characteristics.

A set of jobs (lots) J = {J i |i = 1, . . . , n} has to be processed during a period T by a set of machines M = {M k |k = 1, . . . , m}. Each job J i is made of n i operations such that each operation O ij has a duration p ij > 0 and a set M ij ⊆ M of machines (the furnaces or the cleaning equipments of the diffusion area) able to process it.

Let O k = {O ij ∈ O|M k ∈ M ij }
denote the set of operations that can be assigned to machine M k . The value of p ij and the elements of the set M ij are determined by the recipe of operation O ij denoted ρ ij . In general we have M ij ⊂ M since each machine cannot be configured for all recipes. Each operation O ij has to be included in a batch on a resource k ∈ M ij . Each machine has a finite capacity R k which gives the maximal number of lots in the same batch. On each machine k, S k denote the setup time needed before starting a new batch, D k denote the removal time needed after the completion of a batch and s k denote the constant setup time needed between two different batches. s 0 k denotes the initial setup time on machine k, depending on the state of the resource at time 0. Two consecutive operations O ij and O i(j+1) of the same job are linked by minimal and maximal time lags. Once the batch of O ij is completed and removed from k, the setup for the batch of O i(j+1) cannot start before a minimal time lag τ min ij and has to start before a maximal time lag τ max ij . Let O = {O ij |i = 1, . . . , n; j = 1, . . . , n i } denote the set of all operations. Each job J i has a relative priority c i (c i < c j =⇒ J i is more urgent than J j ). Each job corresponds to a number w i of wafers produced when the job is completed. Finding a feasible solution for the problem lies in making four types of decisions: D1 -Partition the operations into batches, D2 -Select a resource to process each batch, D3 -Order the batches on each resource, D4 -And assign a start time to each batch. Decisions D1-D3 can all be represented by a family of batches B = {B kq } k∈{1,...,m},q∈{1,...,ν k } where B kq is the batch sequenced at position q on machine M k . ν k ∈ {0, . . . , |O k |} denote the number of batches assigned to machine M k . Decision D4 lead to a family of start times

T = {t ij } O ij ∈O assigned to the operations. Once a solution {B, T } is determined, we can deduce a machine assignment {m ij } O ij ∈O where m ij denote the machine O ij is assigned to, i.e. verifying ∃q ∈ {1, . . . , ν k } such that O ij ∈ B m ij q .
To be feasible, a solution {B, T } and its corresponding assignment {m ij } O ij ∈O have to satisfy the following constraints. The operations of the same batch must have the same recipe, i.e.

ρ ij = ρ xy ∀B ∈ B, ∀O ij , O xy ∈ B (1)
Each operation must be assigned to a machine able to process its recipe.

m ij ∈ M ij ∀O ij ∈ O (2)
The batch capacity cannot be exceeded and each batch includes at least one operation.

1 ≤ |B kq | ≤ R k ∀B qk ∈ B (3)
An operation appears in only one batch i.e.

B ∩ B ′ = ∅ B, B ′ ∈ B, B = B ′ (4)
All operations are included in a batch

∪ B∈B B = O (5)
The start time of the first operation of each lot cannot exceed the lot release date:

t i1 ≥ r i ∀J i ∈ J (6)
Each operation O ij , j > 1, cannot start before a minimal time lag after the end of preceding operation O i(j-1) , which takes account of the removal time of the batch of O i(j-1) , the minimal time lag τ min i(j-1) and the setup time of batch of O ij :

t ij -t i(j-1) ≥ D m i(j-1) + p i(j-1) + τ min i(j-1) + S m ij ∀i ∈ {1, . . . , n}, ∀j ∈ {2, . . . , n i } (7)
Each operation O ij , j > 1, has to start before a maximal time lag after the end of preceding operation O i(j-1) , which takes account of the removal time of the batch of O i(j-1) , the maximal time lag τ max i(j-1) and the setup time of batch of O ij :

t ij -t i(j-1) ≤ D m i(j-1) + p i(j-1) + τ max i(j-1) + S m ij ∀i ∈ {1, . . . , n}, ∀j ∈ {2, . . . , n i } (8)
The start times of two operations of the same batch must be equal:

t ij = t xy ∀B ∈ B, ∀O ij , O xy ∈ B (9)
An operation of a batch which is not at the first position on its machine cannot start before the end of the preceding batch on the machine, plus the necessary removal time of the preceding batch, plus the minimal setup time on the machine between two batches, plus the necessary setup time for the next batch.

t ij -t xy ≥ p xy + D k + s k + S k ∀B kq,q>1 ∈ B, ∀O ij ∈ B kq , ∀O xy ∈ B k(q-1) (10)
An operation of a batch in the first position on its machine cannot start before the initial setup time for this batch (we assume

S 0 k ≤ S k + D k + s k , ∀M k ∈ M): t ij ≥ s 0 m ij ∀O ij ∈ O (11)
By definition, a feasible solution includes each operation inside a batch. In the ATMEL Rousset production unit, the scheduling horizon is in fact limited to a scheduling period T . Hence only those batches released in the interval [0, T [ must be taken into account. Several criteria are used to measure the quality of a feasible solution. The number of moves is the number of wafers produced during the time period T:

f mov = J i ∈J w i θ i ( 12
)
where θ i = 1 when the job is completed before T and θ i ∈ [0, 1[ denotes the completion ratio of job i before time T otherwise.

θ i = O ij ,t ij +p ij <=T p ij + O ij ,t ij <T <t ij +p ij (T -t ij ) n i j=1 p ij
The batching coefficient is the average ratio of the actual size of each batch over to its maximal size. Let B T = {B ∈ B|t ij < T, ∀O ij ∈ B} denote the set of batches started before time period T .

batching coefficient = k B kq ∈B T |B kq |/R k |B T | (13) f batch = k B kq ∈B T |B kq |/(R k + Z k 100 ) |B T | (14)
where Z k means the number of total qualified recipe on the machine k.

The weighting of the coefficient of batching by the number of qualified recipes on the concerned machine makes it possible to use the least general-purpose machines initially. The average X-factor is the average of the X-factor of each job weighted by the job priority. All the jobs do not have the same priority. Certain jobs are more important than others (important customers, tests to be carried out quickly, delay to catch up with, etc). Thus, a weight is affected to the job to reflect the priority. Let J T = {J i ∈ J |t in i + p in i ≤ T } denote the set of jobs completed before T .

X-fac =

J i ∈J T (t in i -r i )/p in i |J T | (15) f X-fac = J i ∈J T c i (t in i -r i ) |J T | (16)
The choice made together with the decision makers in the ATMEL Rousset unit is to combine these different objectives into a single one by minimizing the weighted sum

f = αf mov + βf batch -γf X-fac (17)
where α, β and γ are adjustable weights allocated to each objective function. The objective functions have been designed to integrate some requests of ATMEL managers. On the one hand, in the calculation of f X-f ac , the delay of each lot is multiplied by its priority, in order to support the urgent lots. On the second hand, in the calculation of the total batching coefficient, the batching coefficient of each equipment is multiplied according to the number of qualified recipes on this equipment. This leads to occupy in priority the equipments able to treat less recipes, and thus to maintain the availability of the most flexible equipments. The objective of the problem is to search for a feasible selection {B, T } such that f is maximized. Note that given a (feasible) solution {B, T }, f can be computed in O(N ) time where N = n i=1 n i is the total number of operations.

The disjunctive graph representation

The considered problem can be seen as an extension of the job-shop problem and, consequently, the disjunctive graph model can be used for batching and scheduling problems as proposed by [START_REF] Mason | A modified shifting bottleneck heuristic for minimizing total weighted tardiness in complex job shops[END_REF]. In their article, the authors consider a different objective function (namely total weighted tardiness). Sequence-dependant setup times and reentrant flows are also tackled but there are no maximum time lags. Unfortunaly, these latter constraints increase considerably the difficulty of the problem [START_REF] Gentner | Batch production scheduling in the process industries[END_REF]. Let us explain how we modeled the problem using disjunctive graphs. We define the disjunctive graph G = (V, C, E) as follows.

• V is a set of nodes where there is one node per operation, denoted V ij , plus a dummy start node denoted 0.

• C is the set of conjunctive arcs representing the release dates and minimal and maximal time lags. There is an arc from node 0 to node V i1 of each job J i . There is an arc from V ij to V i(j+1) and an arc from V i(j+1) to V ij , for each consecutive operations O ij and O i(j+1) of each job J i . Let B denote a partial or complete batching for the problem satisfying at least constraints (1), (2), (3) and (4). B is a complete batching if constraint (5) is also verified, otherwise it is a partial batching . Recall B also defines the machine assignment m ij , for all O ij ∈ ∪ B∈B B. We assume

m ij = 0 if O ij is not batched in B, i.e. if we have O ij ∈ ∪ B∈B B.
B unambiguously defines a selection S as follows. For each distinct operations Once a selection is computed, we define a graph G(S) = (V, C ∪ S) where arcs C ∪ S are valuated as follows (we assume s 0 0 = s 0 = S 0 = D 0 = 0):

O ij and O xy such that m ij = m xy = k = 0, select arc (V xy , V ij ) k if O xy is in a batch sequenced before the batch of O ij , select arc (V ij , V xy ) k if O ij is in
• Each arc from 0 to the first operation O i1 is valuated by max(r i , S 0 m ij ), the maximal value between the release date of job i and the initial setup time for machine m i1 .

• Each arc from V ij and V i(j+1) is valuated by D m ij + p ij + τ min ij + S m i(j+1)
, the value of the minimal time lag between O ij and O i(j+1) plus the setup and removal times linked to the assignment of the operations.

• Each arc from O i(j+1) and O ij is valuated by -(D m ij +p ij +τ max ij +S m i(j+1)
), the (negative) value of the maximal time lag between O ij and O i(j+1) plus the necessary setup and removal times.

• Each arc (V ij , V xv ) k is valuated either by p ij + D k + s k + S k if the opposite arc is not
selected or by 0 if the opposite arc is selected. Indeed, in the first case this arc represent the decision to sequence O xy after O ij on machine k whereas in the second case, both arcs (V ij , V xv ) k and (V xv , V ij ) k represent the sychronisation of the operations included in the same batch.

We can state that the (partial) solution represented by the (partial) batching B and its selection S is feasible if and only if all longest path problems from node 0 to each node V ij in G(S) have a solution. In the positive case and if B is complete, a feasible schedule T can be obtained by setting t ij to the length of the longest path from 0 to V ij . Furthermore T is the best schedule compatible with B one can obtain when the objective is to maximize f . The problem can be formulated as follows: Find the batching B verifying constraints (1), (2), (3), (4) and ( 5) such that the corresponding selection S is feasible and maximizes f .

As stated in the previous sections, the actual data issued from the fab has other characteristics that have been tackled, such as in-process jobs and machine down times. We can model easily these characteristics thanks to the disjunctive graph formulation. They can be both represented by operations with fixed start times on the involved machine. A start time t can be fixed by linking the node with node 0 by two opposite arcs valued by t and -t. As stated in Section 2, the actual problem also includes an important line management constraint: a minimum distance between the start time of two batches of the same recipe scheduled on two different batches has to be respected. This can also be tackled through the disjunctive graph representation. A fictitious machine can be associated to each recipe and disjunctive arcs linking two operations of the same recipe can be defined. Then, whenever the two operations are batched on two different machines (furnaces), the disjunctive arc has to be orientated. Once orientated the arc is valuated by the minimum required distance.

Solution Methods

We propose a two-phase heuristic method to solve the problem. The first phase is a constructive heuristic based on successive job insertions. The second phase is a local search method which aims at improving the initial solution. Both phases are based on the evaluation of a (complete or partial) selection through longest path calculations.

Evaluation of a partial or complete batching

Any partial or complete batching B and its selection S can be evaluated through the calculation of start time t ij , equal to the longest path from 0 to V ij in G(S) of each operation O ij . To compute such longest paths, since the graph includes arcs with negative weights, we use the Bellman-Ford algorithm which has a O(N |S ∪ C|) time complexity. If the algorithms finds a path of positive length, then the partial or complete solution is unfeasible. Otherwise the algorithm returns start times t ij and the objective function value f can be determined.

Computing an initial solution by a priority rule-based constructive heuristic

The initial solution (selection) is computed by a simple job insertion method. The jobs are first sorted in a list L according to the order: jobs involving maximal time lags first, then increasing release dates, then job priority.

The method starts with an empty batching. Then the jobs are taken in the order given by the list and inserted one by one in the current batching. The insertion of J i is made as follows. Let B = {B kq } k∈M,q∈1,...,ν k denote the current batching including jobs located before J i in L.

For each operation O ij of J i and for each resource k ∈ M ij , there are 2ν k + 1 insertion positions of O ij in the batch sequence of machine M k : indeed, for each batch B kq , we may insert O ij inside batch B kq or create a new batch at any position. Each of these insertion positions is evaluated with the algorithm described in the previous section and the one that maximizes f is kept to update B. If none of the insertion positions is feasible for an operation O ij then all these partial solutions violate a maximal time lag and there exists insertion positions that violate only maximal time lag between O i(j-1) and O ij (the last position on each resource of M ij , for instance). Hence we have to remove O i(j-1) from B, the previous operation of job J i , and insert it at a later insertion position. If this is not feasible, they are not scheduled and delayed from the list of operations. Note that at most n i j=1 k∈M ij 2ν k + 1 insertion positions are tested per inserted lot.

Improving the initial solution by iterative sampling

Iterative random sampling consists in applying iteratively the constructive algorithm presented in Section 4.2 changing each time the input operation list randomly. Two different iterative sampling methods are proposed. The random sampling method initialize defines the input list of each iteration as a random permutation. The pseudo-random sampling method computes the list by making slight random modifications on the sorted list. It consists in dividing the list into intervals and, for each interval, in making a random modification following a uniform law.

Simulated annealing

Simulated annealing (SA) belongs to the class of randomized local search algorithms and has been developed by [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] to handle hard combinatorial problems. The use of simulated annealing supposes the definition of a combinatorial optimization problem and a neighborhood. We are here looking for a solution to the considered batching and a scheduling problem with the maximal objective function value. The neighborhood we consider is defined by the following moves:

1. "Batch move". A batch is moved randomly (both equipment and position can change).

2. "Operation move". An operation is moved randomly in an existing batch or a new batch is created (both equipment and position can change).

3. "Operation switch" Two batches with the same recipe are selected randomly. Then two randomly selected operations are switched.

50% of the considered moves are batch moves, 25% are operation moves and 25% are operation switches. According to the move, randomly generated, we proceed as follows. All random selections are made according to the uniform law. If a move is impossible due to a violation of a constraint, we restore the current solution, we try another move and we reiterate. Each neighbor solution is evaluated y means of the longest path computations presented in Section 4.1.

The algorithm starts with an initial solution s0 and then continually tries to find better solutions by searching in its neighborhood (obtained by the moves described above), and applying a stochastic acceptance criterion. When a neighbor (a solution s) of the current solution is selected we compare the difference ∆f = f (s)f (s0). If the difference ∆f is negative, the neighbor replaces the current solution. Otherwise, if the difference is positive the neighbor is accepted with a probability based upon the Boltzmann distribution P accept (∆f ) ≈ exp( -∆f kT ) where k is a constant and the temperature T is a control parameter. This temperature is gradually lowered following a geometrical law function g(T ) : g(T ) = αT with α < 1.

Computational experiments

We have coded the proposed algorithms in Java and we have tested them on a 1 gigaoctet RAM and 3.4 gigahertz processor computer. These tests have been performed on 8 actual instances of the problem issued for 8 production days in the ATMEL Rousset factory. There are 700 jobs yielding a total number of 1400 operations with about 50 different recipes to schedule on about 70 furnaces, and each of them has a capacity between 4 and 6 lots. There are 12 cleaning equipment with capacity about 2 to 4 lots. The time needed to load and unload a batch varies between 10 and 30 minutes. The minimum and maximum time lags varie from 10 minutes to 4 hours. The target time horizon is 24 hours. The weights of the components of the objective function have been computed both to normalize the different indicators and to take account of the user preferences. For the experiments, we have selected α = 601, β = 1500001, γ = 41.

In Table 1, we display the objective function value (experssion 17), the number of moves expression 12), the batching coefficient (expression 13), the X-factor (expression 15) and the required computational time of the 8 solutions obtained by the constructive method.

From Table 2 to Table 5, we display the results of 1000 runs of the random sampling method, the pseudo-random sampling method and the simulated annealing heuristic for the 4 indicators (objective function, number of moves, batching coefficient, X-factor) separately.

The percentage correspond for each instance to the relative improvement brought by 1000 runs of the considered method over the reference solution obtained by the constructive heuristic, i.e. Best solution-Sorted list Sorted list . We also display the standard deviation over the 1000 runs for each method and each instance.

The computational times are 1000 times greater than the constructive heuristic for the random sampling and pseudo-random sampling methods. The computational time of each run of the simulated annealing algorithm has been limited to 360 seconds.

The analysis of the Table 1 shows that we get the initial solution in a rather fast computation time: the smallest time is 9 seconds and the the largest time is 17 seconds. For the managers of the diffusion area, obtaining a solution in a fast time was a key issue. We underline also that the constructive method provides solutions which are evaluated as satisfactory by the ATMEL managers. They also point out that the results are consistant with the operator practices. This was also an important point for the ATMEL managers who would not accept that the proposed initial solutions change too deeply the organization of the diffusion area. Consequently the first results we obtained were well accepted in the fab and constituted good start points for further improvement.

Let us consider the results of the random and pseudo-random sampling method for the objective function (Table 2). The improvements brought by the random sampling method are better than the ones brought the pseudo random sampling method. The analysis of standard deviation allows us to state that the random list covers a greater solution space, meaning that there is a capability of reaching better solutions, but also worse. Indeed, we obtain an objective function increase of 6, 77% in the best case and of 0, 00% in the worst case for random sampling while for pseudo random sampling, the increase is of 6, 52% in the best case and of 0, 85% in the worst case. Among all instances, the average improvements are of 3,66% for random sampling and 3% for pseudo random sampling. Globally the random and pseudo-random sampling method succeed in improving the objective function at the expense of larger compuational times.

Looking separately at the objective function components, the move indicator, the batching coefficient, the X-factor are slightly improved in average showing the relative quality of the initial solution. The weighting of the delay by the priority of the batches makes it possible to schedule the most urgent batches first. Similarly, the weighting of the batching coefficient by the number of qualified recipes on the concerned equipment makes it possible to fill the least general-purpose machines first. For some instances (060317 for the random method and 060428 for the pseudo-random method), we can notice that the objective function is improved while the indicators are worsened which is a consequence of the selected weights.

For simulated annealing, the initial solution is taken according to the best solution obtained by the random sampling method and the pseudo random sampling method. First, it must be reported that there is no significant difference between starting from a high temperature or a low temperature. The displayed results correspond to a temperature equal to 32000 and the step parameter is set to 25. Considering we start from an initial solution which was already improved by (pseudo-)random sampling, the fact of starting from a high temperature does not have a significant influence on the objective function. Globally, the objective function can be improved up to 12,36 % (starting from the best random sampling solution) and 12,37% (starting from the best pseudo random sampling solution). The average objective function improvement is of 8,31% (random) and of 7,05% (pseudo random).

The "number of moves" indicator is rather difficult to improve. The majority of the tests slightly worsen it but we can reach an improvement on it up to 4,15% (random) and 4,28% (pseudo random). The average improvement for move indicator is of 0,40% (random ) and 0,77% (pseudo random) (see Table 3). The batching coefficient can be improved up to 3,78% (random list) and 4,04% (pseudo random list). For this indicator the average improvements are 1,24% (random list) and 2,16% (pseudo random list) (see Table 4). The X-factor is the indicator where the improvements are the most significant. There are improvements which go up to 5,22% (random) and 6,52% (pseudo random). The average improvement is of 2,81% (random) and 3,37% (pseudo random) (see Table 5). In comparison to the results of random and pseudo-random sampling, we improve in average a lot more each indicator. Again, at the expense of extra computational times, the simulated anealing algorithm succeeds in improving significantly the results of the sampling methods.

The BOS software

The use of a disjunctive graph brings significant improvements for interactive scheduling at the fab level. A prototype software (Figure 1 gives the graphical user interface of the software), the Batch Optimization Software, has been developed jointly by ATMEL, the University of Avignon (LIA) and the École des Mines de Saint-Etienne (CMP Gardanne). It includes the off-line batching and scheduling phase, as described above, and also an interactive module that allows the decision-makers to make modifications and to test options directly on the proposed plan. The longest path calculations in the disjunctive graph allow to almost immediately give the impact of a modification such as job removal, operation move, job or operation insertion. The software is currently in phase of intensive test in the ATMEL diffusion area.

Concluding remarks

We have proposed a model and a method based on a disjunctive graph for batching and scheduling problem in a semiconductor manufacturing factory while taking into account complex constraints and optimizing multiple measures. A constructive algorithm has been proposed to solve the problem, local search improvements based on the disjunctive graph representation have been defined and a simulated annealing algorithm has been developed. The computational tests made on real instances of ATMEL showed that we obtain quality solutions in fast calculation times. The simulated annealing procedure improve significantly the different criteria (number of moves, batching coefficient and X-factor) in comparison with a simple local search method.

This research can be extended in many ways. Other types of moves such as the simultaneous move of 2 jobs linked by a maximum time lag should be tested. The maximal time lags are currently hard constraints. However in practice some of them can be relaxed and treated as soft constraints or objectives. Another important issue can be in a multicriteria analysis study. Computing several Pareto solutions are useful particularly when the situation frequently changes as it is usually the case of semiconductor factories.

  a batch sequenced before the batch of O xy , and select both arcs (V ij , V xv ) k and (V xy , V ij ) k if O ij and O xy are assigned to the same batch.

Figure 1 :

 1 Figure 1: Interface of the prototype

  • E is the set of disjunctive arcs which represent the decisions of the problem. There are two opposite disjunctive arcs (V ij , V xv ) k and (V xy , V ij ) k for any machine k ∈ M and for any pair of operations O ij , O xy ∈ O k , O ij = O xy . These arcs represent the sequencing or the batching decision concerning O ij and O xy on machine k.

Table 1 :

 1 Results of the constructive method

		Random	Pseudo-random	SA (Random) 1	SA (Pseudo-random) 2
	Instance	Best	St. dev.	Best	St. dev.	Best	St. dev.	Best	St. dev.
	060303	0,00%	1,59%	6,52%	0,47%	8,32%	1,22%	12,37%	0,82%
	060310	2,64%	1,78%	1,20%	0,90%	11,32%	1,34%	8,57%	1,18%
	060317	1,53%	1,31%	0,85%	0,41%	9,96%	1,08%	4,69%	0,64%
	060331	1,56%	1,38%	1,56%	0,44%	4,12%	0,48%	4,36%	0,58%
	060407	6,77%	2,17%	6,18%	1,57%	3,45%	0,70%	7,40%	1,00%
	060414	5,38%	2,40%	5,10%	2,10%	6,08%	1,00%	6,23%	0,96%
	060421	5,21%	1,88%	0,90%	0,30%	12,36%	0,85%	5,72%	0,70%
	060428	6,17%	1,66%	1,65%	0,54%	10,91%	0,83%	7,05%	0,81%
	Average	3,66%	1,77%	3,00%	0,84%	8,31%	0,94%	7,05%	0,84%

Table 2 :

 2 Results of the improvement methods (Objective function)

		Random	Pseudo-random	SA (Random) 1	SA (Pseudo random) 2
	Instance	Best	St. dev.	Best	St. dev.	Best	St. dev.	Best	St. dev.
	060303	1,15%	1,92%	-0,66%	0,51%	2,15%	0,76%	0,98%	0,79%
	060310	2,96%	1,76%	1,49%	1,40%	0,48%	1,23%	3,20%	0,78%
	060317	-2,37%	1,69%	-1,43%	0,50%	-3,38%	0,95%	-2,96%	0,59%
	060331	-2,63%	1,71%	0,34%	0,41%	0,29%	0,45%	-0,38%	0,93%
	060407	4,49%	2,42%	1,89%	1,38%	1,59%	0,53%	4,28%	0,75%
	060414	0,29%	1,71%	1,00%	0,88%	-0,87%	0,90%	-1,65%	0,53%
	060421	1,99%	2,38%	0,53%	0,67%	4,15%	0,99%	2,98%	0,63%
	060428	-3,11%	1,70%	-0,99%	0,46%	-1,23%	1,18%	-0,28%	0,48%
	Average	0,35%	1,91%	0,27%	0,78%	0,40%	0,87%	0,77%	0,68%

Table 3 :

 3 Results of the improvement methods (Number of moves)

		Random	Pseudo-random	SA (Random) 1	SA (Pseudo random) 2
	Instance	Best	St. dev.	Best	St. dev.	Best	St. dev.	Best	St. dev.
	060303	2,73%	1,26%	1,36%	0,35%	3,67%	0,73%	2,98%	0,64%
	060310	-3,80%	1,06%	-2,53%	0,61%	-1,68%	0,66%	0,70%	0,51%
	060317	0,00%	1,19%	-1,33%	0,47%	0,16%	0,94%	0,59%	0,56%
	060331	-1,33%	1,16%	-1,33%	0,55%	1,71%	0,37%	2,36%	0,50%
	060407	1,43%	1,32%	1,43%	0,57%	3,78%	0,67%	4,04%	0,53%
	060414	-1,32%	1,29%	0,00%	0,75%	-0,80%	0,50%	2,88%	0,71%
	060421	1,45%	1,28%	1,45%	0,31%	1,02%	0,68%	2,03%	0,50%
	060428	1,43%	1,16%	0,00%	0,51%	2,02%	0,35%	1,72%	0,60%
	Average	0,07%	1,22%	-0,12%	0,52%	1,24%	0,61%	2,16%	0,57%

Table 4 :

 4 Results of the improvement methods (Batching coefficient)

		Random	Pseudo-random	SA(Random) 1	SA(Pseudo random) 2
	Instance	Best	St. dev.	Best	St. dev.	Best	St. dev.	Best	St. dev.
	060303	0,95%	1,62%	0,95%	0,52%	0,86%	0,99%	5,22%	0,16%
	060310	2,35%	1,53%	0,00%	1,56%	4,48%	1,12%	3,58%	0,14%
	060317	-1,22%	1,32%	0,61%	0,56%	1,23%	0,80%	1,24%	0,13%
	060331	0,82%	1,82%	2,45%	0,95%	2,00%	1,03%	4,40%	0,16%
	060407	-1,77%	1,42%	-0,29%	0,54%	5,00%	0,90%	6,52%	0,12%
	060414	-0,34%	1,44%	1,03%	0,86%	4,24%	1,08%	3,08%	0,20%
	060421	-3,73%	1,75%	-1,69%	0,35%	1,50%	1,25%	0,10%	0,12%
	060428	0,28%	1,12%	0,00%	0,28%	3,18%	0,93%	2,82%	0,12%
	Average	-0,33%	1,50%	0,38%	0,70%	2,81%	1,01%	3,37%	0,14%

Table 5 :

 5 Results of the improvement methods (X-factor)

Acknowledgements

This work is part of the MEDEA+ European project HYMNE (High Yield driven MaNufacturing Excellence in sub 65 nm CMOS), partly funded by the "Ministère de l' Économie, des Finances et de l'Industrie" (French Ministry of Economy, Finance and Industry).