Y J Dappe

M Andersen

R Balog

L Hornekaer

X Bouju

Adsorption and STM imaging of polycyclic aromatic hydrocarbons on graphene

Keywords: number(s): 68.37.Ef, 68.43.Fg, 71.15.Mb

The structural characterization of polycyclic aromatic hydrocarbon molecules adsorbed on graphene is of fundamental importance in view of the use of graphene or graphene nanoribbons for electronic applications. Before reaching this point, one has to determine the structure of the adsorbed molecules. Here, we study the case of benzene, coronene, and hexabenzocoronene on a graphene layer. First, the adsorption properties of single molecules are calculated using first-principles calculations at the level of density functional theory. We benefit from a recent scheme, particularly adapted for weakly adsorbed molecules, allowing us to precisely calculate the van der Waals contribution. Then, scanning tunneling microscopy (STM) is used to produce images of self-assembled molecules comparing different theoretical approaches to experimental observations. Finally, we consider the imaging of isolated molecules, and we show how the STM tip influences the molecule position by soft mechanical interaction during the scanning process.

I. INTRODUCTION

Polycyclic aromatic hydrocarbon (PAH) molecules such as benzene, coronene, or hexabenzocoronene (HBC), constitute building blocks for more complex molecules [START_REF] Peña | Ideas in Chemistry and Molecular Sciences: Advances in Synthetic Chemistry[END_REF][START_REF] Alonso | [END_REF]. The use of PAHs as molecular skeletons in combination with chemical functionalization leads to molecules of high interest for molecular electronics [3][4][5][6], for example. Consequently, the study of these building blocks and their interaction with metallic surfaces or graphitic materials is of fundamental importance. The characterization of their structural and electronic properties is a prerequisite for more complex studies such as electronic transport. In that manner, this work is focused on these two aspects, using density functional theory (DFT) for structural aspects, and scanning tunneling microscopy (STM) image calculation for the electronic characterization.

In addition, since the interaction between these molecules and graphene is dominated by weak covalent and van der Waals interactions [7,8], this constitutes also a model system for the study of dispersion interactions. Even though several approaches have been elaborated recently, these interactions remain complicated to handle, especially at the microscopic level. This is mainly due to the nonlocal character of the noncovalent bonding. For example, some methods have been proposed on the basis of DFT [9][10][11][12][13] with relatively successful results, but are limited by the size of the considered systems and the computational resources, depending on the considered calculation scheme. Additionally, inclusion of many-body contributions improves the accuracy of DFT calculations [14].

On the other hand, semiempirical methods remain important references in the field [15,16], but the parameters still have to be fitted for each new system to be studied. Regarding the specific adsorption of PAH molecules based on the DFT formalism, some studies have been achieved [17][18][19] on graphene or metallic surfaces, for structural characterization or influence on graphene electronic properties.

Here, our purpose is twofold. Following an intermolecular perturbation theory combined with DFT [20,21], we determine accurately the structural properties of PAH adsorbed on graphene. Second, we calculate STM images of benzene, coronene, and HBC and compare to experimental results. The influence of perturbations from the tip during the scanning process is investigated from a simple force-field parametrization of the van der Waals interactions between molecule and tip.

In Sec. II, we will present our calculation methodology. We will recall our intermolecular perturbation theory plus DFT to determine the structural configuration of the PAH molecules adsorbed on graphene. Then, we will present the adsorption results in Sec. III, and in Sec. IV A, we will show the STM characterization of these systems. Finally, in Sec. IV B we will discuss the tip influence on the molecule during the scanning process.

II. COMPUTATIONAL METHODS AND EXPERIMENTS

A. DFT calculations

In order to perform STM image calculations of PAHs adsorbed on graphene and to get a valid comparison with experiments, we need to accurately characterize the atomic structures of these systems. Therefore we have used DFT calculations to determine the equilibrium distances of the different PAHs with respect to graphene. Since the PAHgraphene interaction is dominated by π stacking and van der Waals (vdW) interactions, we have used the previously developed LCAO-S 2 + vdW formalism [20,21]. This approach is based on DFT in combination with an intermolecular perturbation theory to describe weak covalent and vdW interactions. Each interacting subsystem, namely in this case the PAH (benzene, coronene, or HBC) and the graphene plane, are calculated separately by DFT. The DFT computational scheme as well as the theoretical foundations underlying our calculations-a very efficient DFT localized orbital molecular dynamics technique (FIREBALL)-have been described in full details elsewhere [22][23][24][25].

We first analyze all the systems involved in this study by using a self-consistent version of the Harris-Foulkes LDA functional [26,27]-instead of the traditional Kohn-Sham (KS) functional based on the electronic density-where the KS potential is calculated by approximating the total charge by a superposition of spherical charges around each atom. The FIREBALL simulation package uses a localized optimized minimal basis set [28], and the self-consistency is achieved over the occupation numbers through the Harris functional [29]. Besides, the LDA exchange-correlation energy is calculated using the efficient multicenter weighted exchange-correlation density approximation (McWEDA) [23,24]. To these DFT calculations we add "weak interactions," which can be seen as two opposite contributions. The first one, named "weak chemical" interaction, is due to the small overlaps between electronic densities of the interacting subsystems. Therefore, this energy can be determined as an expansion of the wave functions and operators with respect to these overlaps. Typically, in graphitic or hydrocarbonated materials this contribution is repulsive [21].

The second contribution is the pure van der Waals interaction, which finds its origin in charge fluctuations in each subsystem, arising from oscillating dipoles whose interaction gives the attractive part of the cohesive energy. This interaction is treated in the dipolar approximation and added in perturbations to the total energy of the system.

The balance of the two contributions gives the equilibrium configuration of the system. This formalism [20,21] has already provided excellent results in the study of graphenegraphene or molecule-graphene interactions [30] or encapsulated molecules in nanotubes [31].

B. STM image calculations

For the monolayer of coronene molecules on graphene we compared two approaches to simulated STM images, namely DFT calculations coupled to the Tersoff-Hamann approximation [32] (DFT-TH) and the elastic scattering quantum chemistry (ESQC) method [33,34]. Whereas ab initio DFT calculations could potentially give a more accurate description of the electronic structure of the system [35][36][37][START_REF] Chen | Introduction to Scanning Tunneling Microscopy[END_REF][START_REF] Janta-Polczynski | [END_REF] compared to the semiempirical extended Hückel approximation implemented in ESQC, STM simulations in the Tersoff-Hamann approximation are rather crude since the STM tip is merely modeled as a sphere with a single s-wave state. The advantage of the ESQC method lies in allowing for the simultaneous description of the surface, the adsorbates in the tunneling junction, the tip apex, and the tip bulk as well as the explicit treatment of the tunneling current. Even though DFT-TH provides good simulated images of bare layered substrates such as graphite and graphene [35,[40][41][42], we will show in the following that its accuracy is inferior to ESQC for coronene on graphene. We have therefore preferred the ESQC method for all subsequent STM calculations of PAHs on graphene. The here reported STM simulations of monolayer coronene on graphene in the Tersoff-Hamann approximation are based on previous DFT calculations [43]. In the Tersoff-Hamann approximation the local density of states (LDOS) is integrated from the Fermi level to the bias voltage. Different tunneling currents can be modeled by considering different isocontour surfaces for the LDOS. In the ESQC method, the tunneling current is evaluated at low bias voltage from the full scattering matrix of the tunneling junction. The system is described by a monoelectronic Hamiltonian, where the matrix elements are evaluated within the semiempirical extended Hückel approximation. The valence orbitals of the atoms are described by a Slater-type basis set. Using standard Slater parameters for carbon, we found that the current perpendicular to the graphene layers of the graphite bulk is rather small [44]. To keep the electronic properties intact and to mimic as close as possible the experimental system, we therefore considered a decoupled graphene layer deposited on a Cu(111) substrate. It is well known that on the copper substrate the graphene layer is weakly bonded and preserves its gapless structure [45,46]. Moreover, as we will show, there is an excellent agreement with the experimental STM images of PAHs on graphite. Indeed, a graphene layer may adopt a particular structure once adsorbed, especially on a rough SiO 2 substrate, on SiC [47,48], on metallic surfaces [46,49] such as Ru(0001), Rh(111), Pt(111), Ir(111), or Ni(111), or on h-BN [50]. Even if a moiré pattern was observed on Cu(111) with STM [51], as well as preferential sticking of molecules on some parts of the moiré [52], one anticipates that the reconstruction of the graphene layer is small enough to consider an unperturbed substrate for molecular adsorption. For the Cu surface and the tip we limited the basis set to a single 4s orbital, except for the apex Cu atom which was modeled with two Slater exponents to allow for long-range interactions with the adsorbates in the tunneling junction. The carbon and hydrogen atoms in the tunneling junction were described by full standard 2s2p and 1s basis sets, respectively. Such a parametrization has been already used with success [53,54]. Notice that the tip apex is ended by a single atom supported by a perfect trigonal aggregate (its structure is threefold). The probe structure does not influence the image of molecular adsorbate because the tip-sample distance is such that the second plane of the tip aggregate does not contribute directly to the main tunneling current. The control of the tip apex structure and its orientation with respect to the surface is much less sensitive in STM than in the noncontact atomic force microscopy (AFM) [55][56][57][58], due to the z dependence of the tunnel current variation which is more acute for STM than for AFM.

C. Experimental setup

Measurements were performed using the home-built Aarhus STM [59] under ultrahigh vacuum (UHV) at a base pressure of 5 × 10 -10 mbar. The grade 1 HOPG substrate was obtained from SPI Supplies. Coronene samples of 99% purity were provided by Sigma-Aldrich. In order to remove the contamination in coronene samples introduced by handling and transfer procedures, the glass crucible filled with coronene was heated in vacuum to 120 • C for 2 hours. The ordered monolayer of coronene was then achieved by heating the crucible up to 180 • C for 120 s while keeping it 5 cm away from HOPG substrate. During deposition, the substrate was kept at room temperature but afterwards slowly (ramp 1 K/s) heated to 100 • C to desorb any coronene multilayers. This approach yields a densely packed and highly stable monolayer of coronene, which allows for performing STM measurements even at room temperature. However, in order to obtain the highest possible stability and resolution, all presented STM images are obtained while keeping the substrate at temperatures of about 140-150 K.

III. ADSORPTION

In this section, we present the main results of our structural calculations for the case of benzene, coronene, and HBC adsorption on graphene. The metallic substrate was not considered as it brings only a van der Waals background energy to the PAH molecule. All these calculations have been performed using the LCAO-S 2 + vdW as described previously. In each case, we have calculated the adsorption potential energy of each molecule as a function of the distance to the graphene plane for several (x,y) positions of the molecules in the plane, and we have kept the optimal energy for each (x,y) position in order to build the potential energy surface (PES). In Fig. 1, we present those PESs for the adsorption of benzene, coronene, and HBC on graphene, respectively. As a main result, we can observe that the three PESs share the same geometric characteristics; i.e., the most stable structure is obtained for an AB-like stacking, similarly to the interaction between two graphene planes. Hence, we obtain a minimum in the potential energy when the center of the molecule coincides with an atom of the graphene plane. Our calculations also indicate that the average distance of these graphitic molecules to the graphene plane is naturally close to the graphene-graphene distance (3.35 Å experimentally, 3.1 Å in the LCAO-S 2 + vdW formalism). More precisely, this distance lies between 2.9 and 3.0 Å for benzene on graphene and 3.0 and 3.1 Å for coronene on graphene, which means that the effect of corrugation is really small on the molecule-graphene distance.

However the energy variation is noticeable, between 0.54 and 0.62 eV, and 1.70 and 2.05 eV per benzene and coronene molecules, respectively. Regarding the HBC molecule, the distance range is more important, between 2.9 and 3.1 Å. This behavior can be explained by the bigger size of the molecule, which results in a larger sensitivity towards the graphene corrugation. The adsorption energy lies between 2.9 and 3.5 eV per molecule. In that case, one can also stress that the molecule-graphene distance is closer to the graphene-graphene distance, as another consequence of the bigger size of HBC. Finally, if we analyze the adsorption energy values with respect to the number of carbon atoms in the molecule, we obtain the following. For benzene, the energy per carbon atom is around 90-103 meV. This value should be compared to our previously published value for the graphenegraphene interaction (between 60-72 meV/atom [20]), which was obtained using the same methodology. While the value for the graphene-graphene interaction is in good agreement with results obtained using the vdW(TS) functional [12] and the optB88-vdW functional [60], our value for benzene is significantly higher. This can be explained by the lack of quadrupole-quadrupole interactions in our model (which remains at the dipole-dipole interaction level), which is important for small molecules, and decreases with the size of the molecule for symmetry reasons, like in the graphene-graphene configuration, for example. In the coronene case, the energy range is around 71-85 meV/atom, in good agreement with quantum chemical calculations at the PM6-DH2 level [61], and DFT calculations using the optB88-vdW functional [60].

As another comparison, the energy value found in this case using the vdW-DF functional gives a slightly lower result, around 63 meV/atom [60]. However this method is known to give results at the lower boundary of the energy range. Finally, the HBC adsorption is characterized by an energy range of about 69-83 meV/atom. As stated previously, the increasing size of the molecules stabilizes the adsorption energy per atom until reaching the value for the graphene-graphene interaction.

In the next section, we will use the above determined structures to model STM experiments on the PAH coronene adsorbed on graphite. Notice that adding more graphene layers to describe the graphite substrate does not modify the present results. Indeed, considering a graphite substrate just adds a background of van der Waals energy due to the large interlayer distance and does not modify significantly the corrugation energy of the PES.

IV. STM CHARACTERIZATION OF PAH ADSORPTION

A. STM images: Theory-experiments comparison

After determining the adsorption geometry of the molecules, we performed STM calculations of PAH molecules adsorbed on a supported graphene layer. In order to check the accuracy of the method, we first consider a monolayer of coronene molecules on graphene. In this case, the adsorption structure has been previously determined from DFT, employing the vdW-optB88 functional [43], and from a LEED study [62]. In Fig. 2 where the molecule appears as a dark spot and the substrate between the molecules is revealed by a bright and large signature, exhibiting a sixfold symmetry. Now, between these two ranges at V t = 1.74 V, we experimentally observe that this last signature shows a triangular shape with three distinct bumps. We recover such a behavior in the calculated image with a higher voltage, i.e., at ∼3.5 V. This discrepancy should be explained by the different substrates considered, but a qualitative agreement is revealed. Apart from this, the concordance between experimental and calculated images is very good. Based on calculated spectra, the change of image appearance is related to the particular molecular resonances probe by the STM bias voltage. The STM images simulated from DFT-TH show varying degrees of agreement with the measured images. At high bias voltage the agreement is quite good, and at low bias voltage there is a similarity, but the simulated image lacks the depression in the center of the molecule. However, for the medium bias voltage-in spite of testing many different bias voltages and LDOS isocontour values-the triangular shape observed in the experimental images and in the ESQC calculated images could not be observed. As mentioned previously, the Tersoff-Hamann approximation is rather crude, and the observed discrepancies could be related to the failure to take into account a realistic tip and the explicit treatment of all tunneling paths in the system, all of which is accounted for in the ESQC method. The good performance of ESQC allows us to have confidence in this method for the study of isolated molecules as described next.

In Fig. 3 we show ESQC calculated STM images of the isolated benzene, coronene, and HBC molecules obtained just below the Fermi level within the HOMO-LUMO gap. Since for benzene, the stacked [Fig. 3(a)] and the 30 • rotated [Fig. 3(b)] configurations are very close in adsorption energy we consider the two situations. At this energy the calculated images of both configurations are circular features almost identical in appearance. For coronene [Fig. 3(c)] and HBC [Fig. 3(d)] the calculated images show a sixfold symmetry with six bright protrusions on the edge of the molecule and a darker contrast towards the center. For all molecules, the calculated images at this energy are nearly independent of the orientation of the molecule with respect to the graphene layer (only shown for the two benzene configurations). At these imaging conditions, it is not possible to resolve the inner structure of the molecule; only an apparent shape can be obtained. In consequence, the determination of the precise position of the molecule with respect to the substrate is difficult without the help of calculations.

In Fig. 4, we show calculated images at 2.4 eV above the Fermi level for the benzene molecule in (a) the stacked configuration, (b) the 30 • rotated configuration, and (c) displaced away from the optimal adsorption site. At this energy there is a pronounced influence of the rotation or translation of the molecule with respect to the substrate.

As a final remark on those calculations, we want to stress here that the considered imaging conditions (low current and high tip-sample distance) generate low electric fields, which allows us to neglect shifts of the molecular levels in the calculations.

B. Influence of the STM tip: Molecular manipulation

As shown in a previous paper [43], the diffusion energy of each molecule on the graphene substrate is relatively low, about tens of meV (see in Fig. 1). This is the reason why it is not trivial to image single molecules without disturbing them during the scanning process. One has to use low tunneling currents to reach high enough tip-sample distance to not displace the adsorbate. We have checked this property by calculating STM images of a single HBC molecule in the constant-height mode: At a fixed tip height (z tip ) and for each (x,y) position of the probe during the scanning, we have calculated the total potential energy of the HBC molecule, that is to say the tipmolecule energy E MT + the molecule-substrate energy E MS . The E MT energy is evaluated by a standard MM4 force-field parametrization for van der Waals interactions [63][64][65]: (the slow scanning direction), until the top right point is reached. One can observe that the molecule is gently attracted and repelled by the tip around its central position. This induces an increase of the apparent shape of the imaged molecule if one compares with the case where the molecule is not free to move (right panels in Fig. 5). This increase reaches at least ∼31% of the apparent area. More interestingly, the shape corresponding to the molecule is largely distorted and it becomes difficult to recover the actual and fine structure of the molecule when the tip is close to the surface. Actually, the tip acts as a van der Waals trap whose depth is modulated by the tip height [66,67]. For lower tip-sample distances, a mechanical action occurs on the molecule with manipulation possibility during the scanning [68][69][70]. Even for z tip = 9 Å, i.e., for a tip-molecule distance greater than 6 Å in our case (taking into account the molecule-graphene distance), the molecule is sensitive to the interaction of the tip. This constitutes an important piece of information in the perspective of using single molecules to functionalize graphene layers or graphene nanoribbons. Indeed, their observation can be altered by the tip-molecule interaction, leading to incorrect interpretations of the experimental results. Nevertheless, it remains that local-probe based methods, such as STM and atomic force microscopy (AFM), are the most adapted techniques to characterize the structure of adsorbates [71], even for those with low adsorption and diffusion energies. On the other hand, this property can also be turned into an advantage, as manipulating molecules opens important potentialities for designing new electronic devices. Here, we have presented a specific case of weakly interacting molecules on graphene, as an academic case of study for molecular manipulation.

E MT = E ab (

V. CONCLUSION

In this work, we have investigated the interaction of PAH molecules with a graphitic substrate. Considering weak covalent and van der Waals interactions, we have determined the energetic landscape of benzene, coronene, and HBC adsorbed on graphene. We have observed that the graphene corrugation influences the adsorption energy very little, except for the HBC molecule due to its bigger size. Also, the molecular adsorption energy per atom tends to the expected graphene-graphene interaction energy per atom when the size of the PAH molecule is increasing, which is in good agreement with previous studies.

For the monolayer of coronene on graphene we have compared STM images calculated from ESQC and DFT-TH to experimental images. Whereas the ESQC method showed a very good agreement with the experiments, we could not recover all experimentally observed features in the images with the DFT-TH method. This could be caused by the very approximate treatment of the tunneling process in the latter method. For isolated PAHs, STM calculations at energies around specific molecular levels exhibited a significant influence of the molecular displacement with respect to the substrate, whereas nothing was observed for energies in the gap. This can be attributed to the symmetry of the molecular levels in interaction with the electronic states of the graphene plane. Finally, since an isolated molecule presents a high mobility with respect to a molecular self-assembly, we have modeled the influence of the STM tip on the molecular position during the scanning process. In particular we have shown a nonnegligible displacement, even at high tip height, leading to a deformation of the STM image. Nevertheless, this property can also be exploited for molecular manipulation to design electronic devices at a fundamental level. Consequently, PAHs on graphene constitute an academic case study of weakly interacting systems for molecular manipulation.

FIG. 1 .

 1 FIG. 1. (Color online) Contour energy plot of the adsorption potential energy on graphene, calculated using the LCAO-S 2 + vdW formalism for (left panel) benzene, (middle panel) coronene, and (right panel) HBC. The (x,y) position refers to the center of the molecule. The color scale refers to potential energy in eV. The graphene lattice is overlaid in gray.

  we compare experimental STM images [(a)-(c)] extracted from large images showing extensive selfassembled monolayer, ESQC-STM calculated images [(d)-(f)] and DFT-TH [(g)-(i)] simulated STM images of the coronene monolayer obtained at energies ranging from around 3 eV below the Fermi level of the substrate to around 4 eV above the Fermi level. The three selected bias voltages in the experimental images are representative for the different imaging modes observed when varying the voltage within the range accessible in STM experiments. For the calculated images many different bias voltages and currents were tested. Using the ESQC method, we were able to reproduce qualitatively very well the experimental images with current and bias conditions close to experimental ones. Here, one has to notice the change of the molecule images with a contrast reversal according to the bias voltage. For bias around -3.0 V [(a) and (d)], each coronene molecule looks like a donut with a circular protrusion for the peripheral rings of the coronene, and a small depression corresponding to the central ring. An opposite situation occurs for bias voltage around 4.2 eV [(c) and (f)],

FIG. 2 .

 2 FIG. 2. (Color online) (a)-(c) Experimental STM images of a monolayer of coronene on graphite. Tunneling conditions are (a) V t = -3.11 V, I t = 70 pA, (b) V t = 1.74 V, I t = 120 pA, and (c) V t = 4.14 V, I t = 110 pA. (d)-(f) ESQC calculated constant-current STM images of a monolayer of coronene adsorbed on graphene with the superimposed atomic structure. The graphene lattice is depicted in gray whereas the coronene is depicted with carbon atoms in dark blue and hydrogen atoms in white. Tunneling conditions and color scales (black to yellow) are (d) V t = -2.83 V, I t = 70 pA, [5.75 Å, 6.92 Å], (e) V t = 3.46 V, I t = 23 pA, [6.73 Å, 7.91 Å], (f) V t = 4.20 V, I t = 2.2 pA, [5.11 Å, 6.36 Å]. (g)-(i) DFT-TH simulated STM images. Tunneling conditions are (g) V t = -2.83 V, I t = 500 pA, (h) V t = 3.46 V, I t = 500 pA, (i) V t = 4.20 V, I t = 500 pA.

FIG. 3 .

 3 FIG. 3. (Color online) Constant-current calculated STM images of (a) benzene in stacked configuration, (b) benzene rotated 30 • , (c) coronene in stacked configuration, and (d) hexabenzocoronene in stacked configuration. The coloring of the superimposed atomic structures is the same as in Fig. 2. Tunneling conditions and color scales (black to yellow) are V t = -0.1 V, I t = 1 pA, [5.25 Å, 6.30 Å] for all images.

FIG. 4 .

 4 FIG. 4. (Color online) Constant-current calculated STM images of (a) benzene in stacked configuration, (b) benzene rotated 30 • , and (c) benzene displaced away from the optimal adsorption site. The coloring of the superimposed atomic structures is the same as in Fig. 2. Tunneling conditions are V t = 2.4 V, I t = 1 pA for all images. Color scales (black to yellow) are (a) [5.15 Å, 5.89 Å], (b) [5.15 Å, 5.97 Å], and (c) [5.15 Å, 7.02 Å].

  FIG. 5. (Color online) Comparison between the STM imaging of hexabenzocoronene in constant-height mode at (a) 9 Å, (b) 10 Å, and (c) 11 Å when the molecule is allowed to move under influence of the tip (left panels) and when the molecule is kept fixed (right images). The images are obtained at the Fermi energy.

ACKNOWLEDGMENT

Part of this work was performed using high-performance computing resources from the Calcul en Midi-Pyrénées (CALMIP) facilities (Grant No. 2011-[P0832]).


A S Albana1

K Aroui

G Alpan

Y Frein

MIXED MODEL ASSEMBLY LINE SEQUENCING TO MINIMIZE DELAYS USING META-HEURISTICS

Keywords: Mixed Model Assembly Line, Sequencing Problem, Meta-heuristics

Nowadays, the customers are increasingly demanding, pushing the companies to offer highly diversified products. This requires that different kinds of products be manufactured in intermixed product sequences on the same line. Such assembly lines are called the mixedmodel assembly lines (MMAL). Workers and machinery have to be flexible to reduce the setup times and costs. A good vehicle sequence in MMAL can have many positive effects on MMAL: It can permit producing more products in a shorter time period (providing cost reductions), it can also improve the work conditions by balancing the workload of the operators. In this article, we are interested in the sequencing of a mixed model assembly line for Truck Industry. In the literature, different objectives exist to solve the MMAL sequencing problem. In this article, we present methods to minimize the total work overload. In a previous work, a linear programming (LP) approach has been proposed for this problem. The MMAL is known to be an NP-hard problem. The exact methods such as LP can only handle small problems and their applications are limited in an industrial context. Therefore, we present here solutions based meta-heuristics. Three types of meta-heuristic algorithms are used in this research: Genetic Algorithm (GA), Simulated Annealing (SA), and finally a hybrid method based on both Genetic Algorithm and Simulated Annealing (GASA). Numerical tests are carried out to compare the performance of the proposed algorithms. For small instances, a benchmark data from the literature is used to compare the performance of the meta-heuristics versus the optimal solution found by the LP approach, based on the computational time and the quality of the solutions. The comparisons are also made for larger instances, for some generated data and the data from an industrial case study.

INTRODUCTION

Nowadays, the customer need in specific product increases. This forces the manufacturer to produce different types of products. High product variety makes the production management more difficult. The use of flexible workers and machinery is a common solution so that different kinds of products can be manufactured in intermixed product sequences on the same line, reducing the setup times and cost. Such production lines are called the mixed-model assembly lines (MMAL) [START_REF] Boysen | Sequencing mixed-model assembly lines: Survey, classification and model critique[END_REF]. This type of assembly line is generally used in automotive industry and consumer goods industries such as electronics, white goods, furniture and clothing.

There are two main problems in MMAL: assembly line balancing and product sequencing. Assembly line balancing problem deals with the assignment of tasks to workstations. The most common objective is to minimize the number of stations needed to manufacture a product in a line given a fixed cycle time, equivalent to a fixed production rate [START_REF] Bautista | Ant algorithms for assembly line balancing[END_REF]. And sequencing problem deals with sequencing different product models launched on an assembly line, so that work overload at the stations induced by direct succession of multiple labor-intensive models is avoided or part rate usage is minimized [START_REF] Boysen | Sequencing mixed-model assembly lines to minimise the number of work overload situations[END_REF]. Assembling various options leads to variations in processing times at work stations [START_REF] Boysen | Sequencing mixed-model assembly lines: Survey, classification and model critique[END_REF]. In automotive production, for instance, the installation of an electrical sunroof requires a different amount of time than that of a manual one. The work overload appears when the operator cannot finish the required tasks on a product within the predefined time window. The work overload refers to the remaining work. If several work intensive models follow each other at the same station, work overloads might occur, which need to be compensated, e.g., by additional utility workers. Work overloads can be avoided if a sequence of models is found, to alternate the high work-intensive models with the less work-intensive ones. This Work overload is referred as delay in our research.

Numerous researchers have worked on both problems. Amen [START_REF] Amen | Heuristic methods for cost-oriented assembly line balancing: A survey[END_REF], Chaves et al. [START_REF] Chaves | Hybrid metaheuristic for the assembly line worker assignment and balancing problem[END_REF], Bautista and Pereira [START_REF] Bautista | Ant algorithms for assembly line balancing[END_REF], and many others work on assembly line balancing. Hyun et al. [START_REF] Hyun | A genetic algorithm for multiple objective sequencing problems in mixed model assembly lines[END_REF], Ponnambalam et al. [START_REF] Ponnambalam | Genetic algorithms for sequencing problems in mixed model assembly lines[END_REF], Cano-Belmán et al. [START_REF] Cano-Belmán | A scatter search based hyper-heuristic for sequencing a mixed-model assembly line[END_REF], Rahimi-Vahed et al. [START_REF] Rahimi-Vahed | A new particle swarm algorithm for a multi-objective mixed-model assembly line sequencing problem[END_REF], Aroui et al. [START_REF] Aroui | Minimizing work overload in mixed model assembly lines: A case study from truck industry[END_REF] and many others work on product sequencing. In this article we will deal with the product sequencing problem.

On product sequencing problem, some of the researchersconsider multi-objectives, such as: minimizing total utility work, minimizing total setup cost and minimize total production rate variation [START_REF] Hyun | A genetic algorithm for multiple objective sequencing problems in mixed model assembly lines[END_REF], minimize total production rate variation and minimizing total setup cost [START_REF] Ponnambalam | Genetic algorithms for sequencing problems in mixed model assembly lines[END_REF], and the other work in single objective minimizing, such as: total utility work (Cano-Belmán et al. [START_REF] Cano-Belmán | A scatter search based hyper-heuristic for sequencing a mixed-model assembly line[END_REF]). Some of them solve it using exact methods, e.g. Mixed Integer Linear Programming [START_REF] Aroui | Minimisation des retards dans le séquencement des véhicules sur une ligne d'assemblage multi modèles[END_REF], Bounded Dynamic Programming [START_REF] Bautista | Heuristics and exact algorithms for solving the Monden problem[END_REF]. Some use meta-heuristic, e.g. Genetic Algorithm (GA) [START_REF] Hyun | A genetic algorithm for multiple objective sequencing problems in mixed model assembly lines[END_REF], [START_REF] Ponnambalam | Genetic algorithms for sequencing problems in mixed model assembly lines[END_REF], [START_REF] Rahimi-Vahed | A multi-objective scatter search for a mixed-model assembly line sequencing problem[END_REF], Particle Swarm Algorithm [START_REF] Rahimi-Vahed | A new particle swarm algorithm for a multi-objective mixed-model assembly line sequencing problem[END_REF], [START_REF] Mirghorbani | A multi-objective particle swarm for a mixed-model assembly line sequencing[END_REF], Greedy Randomized Adaptive Search Procedure (GRASP) [START_REF] Alpay | GRASP with path relinking for a multiple objective sequencing problem for a mixed-model assembly line[END_REF], and Ant Colony Optimization [START_REF] Zhu | Ant colony optimisation with elitist ant for sequencing problem in a mixed model assembly line[END_REF].

In 2013, Aroui et al. [START_REF] Aroui | Minimisation des retards dans le séquencement des véhicules sur une ligne d'assemblage multi modèles[END_REF] work on the sequencing problem of MMAL. They choose to act directly on the work overload to minimize the delay. This kind of approach is poorly developed in the literature [START_REF] Aroui | Minimisation des retards dans le séquencement des véhicules sur une ligne d'assemblage multi modèles[END_REF]. Their work is based on an industrial case of Bourg-en-Bresse plant of Renault Trucks. The truck assembly line is a typical MMAL with highly diversifying products compared to automotive industry. The number of vehicles to produce per day is much lower. They present a mixed integer linear programming (MILP). The main findings of Aroui et al. [START_REF] Aroui | Minimisation des retards dans le séquencement des véhicules sur une ligne d'assemblage multi modèles[END_REF] is that the sequence generated at the end of 2 hours is a better solution than the actual procedure at Renault Trucks, however, no optimal solution can be found.

The method presented in Aroui et al. [START_REF] Aroui | Minimisation des retards dans le séquencement des véhicules sur une ligne d'assemblage multi modèles[END_REF] gives good results, despite the long calculation time. In reality, managers need to make fast decisions, especially in production lines. To make a fast decision with a good solution quality, an approach less expensive in computation time is needed.

Based on problems from Aroui et al. [START_REF] Aroui | Minimisation des retards dans le séquencement des véhicules sur une ligne d'assemblage multi modèles[END_REF], this research aims to find another approach for solving thesame problem with faster calculation times but with the same solution quality. We use three types of heuristics, Genetic Algorithm (GA), Simulated Annealing (SA) and combined Genetic Algorithm -Simulated Annealing (GASA). Then, we compare the results of the algorithms (GA, SA and GASA) among themselves and with the results of Aroui et al. [START_REF] Aroui | Minimisation des retards dans le séquencement des véhicules sur une ligne d'assemblage multi modèles[END_REF] based on the performance and computational time.

DESCRIPTION OF THE PROBLEM

Mixed Model Assembly Line (MMAL) is a special assembly line where there is a single production line used for assembling multiple type of products. Normally this type of production system can be found in the automotive industry [START_REF] Boysen | Sequencing mixed-model assembly lines: Survey, classification and model critique[END_REF]. An illustration of a typical MMAL is given in Figure 1. The products in MMAL move on a continuous transportation system such as a belt conveyor. Each product has its own processing times, however, the cycle time at the workstations is constant.

Figure 1 Illustration of MMAL

Based on Boysen et al. [START_REF] Boysen | Sequencing mixed-model assembly lines: Survey, classification and model critique[END_REF],sequencing problemscan be categorized as three types:

1. Mixed-model sequencing: aims at minimizing sequence-dependent work overload based on a detailed scheduling which explicitly takes operation times, worker movements, station borders and other operational characteristics of the line into account.

Car sequencing:

To avoid the significant effort of data collection that accompanies mixed-model sequencing, car sequencing attempts to minimize sequence-dependent work overload in an implicit manner. 3. Level scheduling: While the first two approaches aim at minimizing violations of capacity constraints, level scheduling seeks to find sequences that are in line with the JIT-philosophy. For this purpose "ideal" production rates are defined and models are sequenced in such a manner that deviations between actual and ideal rates are minimized.

Each of those types has more detailed sub types; based on the type of the production line, station boundary, objective functions, etc. (for further detail see [START_REF] Boysen | Sequencing mixed-model assembly lines: Survey, classification and model critique[END_REF]. In this research, wefocus on the Mixed Model Sequencing Problem where the objective function is to minimize the work overload (MMSP-W). For our type of problem, we have several hypotheses:

-Line balancing has already been done. All tasks are assigned to different workstations, -Products move on the line at a constant speed. There are no inventory buffers between workstations, -Products require a specific amount of work (tasks with operation time) for each position. The operation times are deterministic. Setup time and the time required by operators to return to his initial position at the end of a task are included in the operating time.

The operators here areregular operators who areassigned tasks on every vehiclem i (see figure 2). Every vehicle requires a specific amount of work (tasks with operation time) for each position. And the cycle time is the time between two successive tasks. We denote cycle time as  . A delay happens when the task on a vehicle cannot be completed within the cycle time (see figure 2).

Figure 2 Illustration of delay in operator (Aroui et al. [11])

As mentioned before, Aroui et al. [START_REF] Aroui | Minimisation des retards dans le séquencement des véhicules sur une ligne d'assemblage multi modèles[END_REF] formulate this problem by MILP. For the sake of completeness, we recall their model below:

𝑀𝑖𝑛 𝑟 𝑗𝑝 𝑛 𝑗 =1 𝑝 𝑝=1 (1)
Subject to

𝑐 𝑗𝑝 = 𝑟 𝑗 -1,𝑝 + 𝑑 𝑖𝑝 𝑥 𝑖𝑗 𝑖=1 ∀𝑗∀𝑝 (2)
𝑟 𝑗𝑝 ≥ 0 ∀𝑗∀𝑝 (3) 𝑟 𝑗𝑝 ≥ 𝑐 𝑗𝑝 ∀𝑗∀𝑝 (4)
𝑥 𝑖𝑗 𝑛 𝑗 =1 = 1 ∀𝑖 (5)
𝑥 𝑖𝑗 𝑛 𝑖=1 = 1 ∀𝑗 (6)
𝑥 𝑖𝑗 ∈ 0,1 ∀𝑖∀𝑗 (7)
𝑐 𝑗𝑝 , 𝑟 𝑗𝑝 ∈ 𝑅 ∀𝑗∀𝑝

Where, Objective function (1) minimizes the total delay. Constraint (2) establishes that the delay or idle time of the product positioned in j is its own delay (or idle time) in addition to overload of the position𝑗 -1. Constraints ( 3) and ( 4) indicate that 𝑟 𝑗𝑝 = 𝑚𝑎𝑥 0, 𝑐 𝑗𝑝 since the objective function takes into account only the delays (and not idle times). Constraint [START_REF] Chaves | Hybrid metaheuristic for the assembly line worker assignment and balancing problem[END_REF] guarantees that only one position can be assigned to each product; constraint [START_REF] Hyun | A genetic algorithm for multiple objective sequencing problems in mixed model assembly lines[END_REF] indicates that only one product can be assigned in each position of the sequence; and, finally, constraint [START_REF] Ponnambalam | Genetic algorithms for sequencing problems in mixed model assembly lines[END_REF] requires that the assigned variables are binary.

𝑛

META-HEURISTIC APPROACH

The objective is to find solutions having a similar quality reported by Aroui et al. [START_REF] Aroui | Minimisation des retards dans le séquencement des véhicules sur une ligne d'assemblage multi modèles[END_REF], but faster in computation time. Meta-heuristics are chosen because they are well adapted for difficult problems in industrial context. A well-designed meta-heuristic method can usually provide a solution that is nearly optimal in short execution time.In this research, we propose three meta-heuristics: GA, SA and combination of both (GA-SA). We use two cases for numerical tests: small instances (academic instances) and big instances (Renault Trucks Case). The resultsare compared to assess the performance of each meta-heuristic. The analyses are done based on the quality of the solution and the computational time.

Genetic Algorithm

Genetic Algorithm (GA) is a meta heuristic founded by John Holland in the 1960s [START_REF] Yang | Nature-Inspired Metaheuristic Algorithms[END_REF]. This algorithm is based on the Darwin's theory of natural selection. There is an initial population consisting of different chromosomes. These chromosomes represent the solution: in our case it is a vehicle sequence. This population will evolve because of genetic operations applied to improve the solution. There are three basic elements (operations) in Genetic Algorithm: reproduction, crossover, and mutation. Another element of the GA is called the fitness function. The fitness function shows the ability of individuals to survive in the population [START_REF] Yang | Nature-Inspired Metaheuristic Algorithms[END_REF], [START_REF] Kim | Sequencing in mixed model assembly lines: a genetic algorithm approach[END_REF] and is measured by the objective function to be optimized [START_REF] Kim | Sequencing in mixed model assembly lines: a genetic algorithm approach[END_REF]. For the maximization problem, the fitness is equalor linear to the objective function (f(x)). But for the minimization problem, the fitness is expressed as the inverse of the objective function (1/f(x)) [START_REF] Ahmed | Genetic algorithm for the traveling salesman problem using sequential constructive crossover operator[END_REF]. InSantosa & Willy [START_REF] Santosa | Metoda Metaheuristik Konsep dan Implementasi[END_REF]the fitness function is given as in equation ( 9). Constant 1 ensures that the fitness does not tend to ∞, when the f(x) = 0.

𝐹 𝑥 = 1 1+𝑓 𝑥 (9)
Where,

𝐹 𝑥

: Fitness function, 𝑓 𝑥 : Objective function.

In GA, there are several parameters, such as: number of chromosomes in the population, number of elite chromosomes, crossover probability, and mutation probability. Crossover probability (pc) is usually very high, typically in the range of 0.7 ~1.0. The mutation probability (pm) is usually small (usually 0.001 ~ 0.05). If pc is too small, the crossover occurs sparsely, which is not efficient for evolution. If the mutation probability is too high, the solutions could still 'jump around' even if the optimal solution is approaching [START_REF] Yang | Nature-Inspired Metaheuristic Algorithms[END_REF]. Elite chromosomes are the chromosomes, which will be carried over to the new generation without being modified. The purpose of the elitism is to keep a good solution in the population. The pseudo code for the proposed genetic algorithm is given in Figure 3.

We note that some improvements are brought to the classical genetic algorithm that doesn't perform well using the parameters mentioned by Yang [START_REF] Yang | Nature-Inspired Metaheuristic Algorithms[END_REF]. So, we improve the native GA.

In improved GA, we do several things:

-Crossover and Mutation are always done in each iteration (p c = p m = 1).

-The offspringis only accepted if it has better fitness function than their parent.

-Mutation is doneon the three best chromosomes with different types of mutation (Swap, Flip, and Slide).

Figure 3Pseudo code of Improved Genetic Algorithm for MMSP-W

The details of how to choose the parents, crossover and mutation will be described in the next section.

Chromosome & Genes

The essence of genetic algorithms isthe encoding of an optimization function as arrays of bits or character strings to represent the chromosomes,the manipulation operations of strings by genetic operators, and the selection according to their fitness with the aim to find a solution to the problemconcerned [START_REF] Yang | Nature-Inspired Metaheuristic Algorithms[END_REF]. These arrays of bits or strings are known as genes.

In our case, a chromosome is a sequence of products. And genes are product i in position jin the vehicle sequence.

Figure 3.4 Representation of Chromosome and Genes

Elitism

The best individuals with higher fitness should be preserved and passed onto the next generation. Elitism is a process to select the most fit individual (in each generation) which will be carried over to the new generation without being modified by genetic operators [START_REF] Yang | Nature-Inspired Metaheuristic Algorithms[END_REF].

In this research, we keep 4 individuals in each iteration as elite chromosomes.

Parent Selection &Crossover

Parent is the selected chromosome to do the crossover process. And parent selection is done based on fitness.Parent with high fitness ismatched with another high fitness parent to do reproduction (crossover). Then low fitness parent is matched with another low fitness one. Parent that has been chosen cannot be chosen another time.

Crossover is the exchange of genetic material between parent chromosomes that results in newchromosomes (child). Crossover is done for all chromosomes, except the elite chromosomes. Parentsare matched based on probability. More detailed explanation can be seen in Figure 5. In this example, there are 6 parents. Parent 1 is matched with parent 3 in the first draw. Then, in the second draw among the remaining parents, parent 2 is paired with parent 4 (according to probabilities). And finally, parent 5 is paired with parent 6. With this scheme, all parents have a pair and do the crossover process. No parents are matched with themselves or do the crossover twice.

Figure 5 Parent selection and crossover process

In crossover itself, only one point crossover is used. It is when a single crossover point on both parents' chromosome is selected. All data beyond that point in either chromosomestring is swapped.The crossover point is chosen arbitrarily. In general, the crossover process is done as illustrated in Figure 6.

Figure 6 Crossover process

We have two parents, A and B, and one point crossover after gene number 2 (position 2 is randomly picked). After the crossover operation, we obtain 2 children, X and Y. Child X inherits his first two genes from parent A and the remaining 3 genes from parent B. (the inverse is true for child Y). None of the children is a valid offspring since a gene is repeated twice and there is one product that hasn't been listed. Here we have to revise part before the crossover point (gene no. 2). For the child X, we have to revise the first two genes. And the product that hasn't been listed is product no. 5. Alternatively, we have two options to replace the first two genes. It can be 5-2 or 2-5. Because this part is coming from parent A and parent A has a sequence of 1-2-3-4-5 (see Figure 6), we put the sequence of 2-5 instead of 5-2 for the child X. The same process is done for the second child. When the both childrenare valid, we check their fitness. If the child's fitness is better than the parent's fitness, we replace the worst parent with the child, otherwise the parent stays in the population. This process is done to keep the better individual in the population like in the elitism procedure. So, if the fitness is not good enough, we can say that the child cannot survive in the population or die in the next iteration.

Mutation

Mutation is a genetic operator used to maintain genetic diversity from one generation to the next. In our case, Mutation is done to the three out of four elite chromosomes in each iteration. Mutation is done using three procedures: swap, flip, and slide.All of this process (which gene will replace which gene) is done arbitrarily. We generate two random numbers that correspond to the j th position in the sequence (j th gene in the chromosome).

Figure 7Mutation process

Swap is done by swapping two genes. For example, we generate two random numbers corresponding to positions 1 and 5. Then we swap between product 2 and product 1 (figure 7.a). Flip or inverse is done by inverting the chosen gene sequence.In figure 7.b, we chose randomly position 4 (which is product 3 here) and position 2 (which is product 5), then we invert the sequence of 5-4-3 to 3-4-5.And slide is done by sliding the genes, for example, based on a random pick, we slide position 4 to position 2. Here we get product 3 (which is in position 4) and then we slide it to position 2 of the sequence. So, we change the sequence from 5-4-3 to become 3-5-4 (fig 7.c).

Simulated Annealing

Simulated annealing is founded by Kirkpatrick et al. [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. This algorithm is based on the annealing process of steel. Simulated Annealing (SA) has been widely used for solving combinatorial problems such as VRP, TSP, scheduling and many more. Simulated Annealing is also known as a meta-heuristic that has the capability to get out of local optima: It has a schematic process to accept a worse solution to avoid being trapped in a local optimum.

In the literature review, there are not many researchers that have used the Simulated Annealing to solve the Mixed Model Sequencing Problem. The present work illustrates the application of SA on MMSP-W and provides some numerical tests on the performance of this type of problem. The pseudo code of simulated annealing that we use to solve the MMSP-W is shown in Figure 8.

Figure 8Pseudo code of SA

One advantage of the SA method is that it can accept a worse solution based on the acceptance probability, to avoid local optima. The acceptance probability is formulated as shown in eq. [START_REF] Aroui | Minimizing work overload in mixed model assembly lines: A case study from truck industry[END_REF]. Based on the pseudo code (Figure 8), we have two loops: inner loop (iteration) and outer loop (cycle). Outer loop (cycle) affect the probability of acceptance and inner loop (iteration) affect temperature degradation.The temperature update or temperature degradation happens after several iterations (inner loop), meaning that temperature isn't always decreasing at each iteration (inner loop). One disadvantage of the method is that the performance of the algorithm depends highly on the good choice of the parameters.

𝑃 𝐸 = 𝑒 -𝐸 𝑘𝑇

In this SA we do the neighbourhood search using flip. This procedure are done by flipping or inversing the sequence of selected position. It is the same procedure as explained in the mutation of GA (figure 7.b).

GASA

The idea is to benefit from the strengths of each meta-heuristic to improve the method. The initial population of the GA is randomly generated. If this initial population gives already good solutions, we may hope to reach a good solution faster. Using this assumption, we propose to inject the results of SA as initial population to GA. Concretely, we run the simulated annealing several times, and then we take the results of simulated annealing as the initial population for the Genetic Algorithm. We refer to it as GASA.

NUMERICAL TEST

The testsare done to evaluate the performance of Genetic Algorithm. They are divided into two sets of instances: Small Instances (Academic Instances) and Big Instances (Real Case from Industrial Data). Small instances are used to check whether the algorithm works correctly or not. Big instances are used to check the performance of the algorithm on real case. An instance is characterized by the number of workstations, the number of products to sequence and the process times of products on a given workstation. The big instances come from an industrial case study(Renault Truck's Bourg-en-Bress Plant). Tests for both instances are performed, using MATLAB, in Intel Core i5 -2.4 GHz -4 GB RAM. The results are compared to the results obtained by Mixed Integer Linear Program (MILP) developed in Aroui et al [START_REF] Aroui | Minimisation des retards dans le séquencement des véhicules sur une ligne d'assemblage multi modèles[END_REF], based on the computational time and quality of the result (efficiency). Note that, MILP is solved using CPLEX on the same computer.

The small instances are composed of three sets of data for multi workstation as shown in Table 1 to Table 3. In these small instance sets, there are for 4 workstations and 20 products. Efficiency is calculated as in eq. [START_REF] Aroui | Minimisation des retards dans le séquencement des véhicules sur une ligne d'assemblage multi modèles[END_REF].

𝐸𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑐𝑦 = 1 - 𝐴𝑙𝑔𝑜𝑟𝑖𝑡 𝑚 𝑆𝑙𝑜𝑢𝑡𝑖𝑜𝑛 -𝑂𝑝𝑡𝑖𝑚𝑎𝑙 𝑆𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑂𝑝𝑡𝑖𝑚𝑎𝑙 𝑆𝑜𝑙𝑢𝑡𝑖𝑜𝑛 × 100% (11)
10 trials are run for each of the instance regardless the size of the instances. Tables 1 to 3 show the best run and the worst run among the 10 runs. We note that initial test runs are performed for the SA, to find the most appropriate parameter setting. Based on these tests, the c= 0.8 and Initial temperature = 10.000. For small instances, Genetic Algorithm (GA), Simulated Annealing (SA) and combined Genetic Algorithm -Simulated Annealing (GASA) don't have much difference. In almost all cases, they can reach the optimal solution. But in terms of calculation speed, SA provides the fastest solutions (optimal solutions found in 3 seconds, with very low variability around this mean). The execution time performance of the MILP model, however, depends on the instance set.

For the big instances,we have 9 instances with 56 ~ 61 products and 77 workstations, each corresponding to a day's production.In those instances, GA, SA and GASA behave differently. Table 4 to Table 6 showresults of each algorithm. For these instances, GA gives very good results. Note that the MILP is stopped at the end of 1 hour and the best solution found is recorded on Table 4. The Genetic Algorithm can achieve solutions -3% ~ 7% better than the MILP solutions. The calculation time is also much faster than the MILP. Here we can see that SA provides low execution times. We only need 17 seconds to execute SA. We get good results, even if this result isn't as good as MILP. On the contrary, GA needs more execution time to achieve the same quality result as MILP. For the big instance sets,GASA's calculation time is faster (7-8 minutes) than the MILP model (1 hour of execution) and the quality of the solution issimilar to MILP. For the best runs, GASA improves solutions within the range of 1% to 7%.

We can conclude that, in terms of calculation speed, SA is the best method as can be seen in Table 4. SA only needs 17 seconds of execution.And, in terms of quality of solution, GA and GASA are the best method among the MILP and meta-heuristic methods.

CONCLUSION

We successfully applied several meta-heuristics for the resolution of the Mixed Model Assembly Line Sequencing Problem, in the case of delay minimization. The meta-heuristics achieved good quality solutions with faster computational time compared to the MILP model. The Genetic Algorithm gets the optimal solution in all small instances. For the big instance, GA gives very good results; in an average of 5 ~ 6 minutes it achieves solutions 0% ~ 7% better than the MILP solutions of 1 hour's calculation. The GA's calculation time is also much faster than the MILP. GA is more efficient than MILP for these instances. Despite good performance in small instances, SA is less performing for big instances than GA and MILP, in terms of quality of the obtained results but in less calculation times. GASA, on the big instances, has a calculation time between 7-8 minutes and the quality of the solution is similar to MILP. GASA even improves solutions within the range of 1% to 7% for the best run. GA or GASA are good procedures for industrial users. Both are sufficiently rapid and give better solution than manual solution This research can still be improved on several aspects such as: type of operators, the methods, and the objectives. As reported by Aroui et al [START_REF] Aroui | Minimizing work overload in mixed model assembly lines: A case study from truck industry[END_REF] in their further study, they reveal that there are other types of operators in the industrial case, for instance, operators working only on specific vehicles or operators that are working in tandem with other operators. The meta-heuristics can be extended to model such operators as well. In terms of method, SA can further be improved by modifying the neighbourhood search. Indeed, in this article, we only considered neighbours generated by a flip. We can also look at the Dynamic Programing as an alternative to meta-heuristic method. In terms of objectives, the single objective can be extended to multi objective case to take into account the industrial reality. For instance, minimizing delay and minimizing the part rate usage at the same time can be interesting to model.

  required by product i for operator p 𝑑 𝑖𝑝 : Delay or tardiness for operator type 1, 𝑑 𝑖𝑝 = 𝑡 𝑖𝑝 -𝛾 𝑥 𝑖𝑗 : Equal to 1 if product i is assigned to position j, otherwise 0 𝑐 𝑗𝑝 : Delay or idle time for operator p for product on position j 𝑟 𝑗𝑝 : Total Delay operator pfor product on position j

  <= Old Solution; Solution = New Solution elseif New Solution > Old Solution; Delta = (New Solution -Old Solution); Probability = random Number; Acceptance = exp(-Delta/(k*T)); if Probability <= Acceptance; Solution = New Solution; end end Iteration = Iteration + 1; end Temperature Update = T * c; cycle = cycle + 1; end

Table 1Genetic Algorithm results for 4 workstations 20 products Problems Linear Programming Genetic Algorithm Optimal Solution Calculation Time (Second) Objectives Function Calculation Time (Second) Efficiency

 1Genetic

				Worst	Best	Longest	Shortest	Worst	Best
	1	23.6	1373.00	23.8	23.6	20.25	6.40	99%	100%
	2	16.0	136.00	17.5	16.0	23.62	5.69	91%	100%
	3	5.5	311.00	5.6	5.5	39.58	5.44	98%	100%

Table 2 Simulated Annealing results for 4 workstations 20 products

 2

		Linear Programming			Simulated Annealing
	Problems	Optimal Solution	Calculation Time (Second)	Objectives Function Calculation Time (Second) Worst Best Longest Shortest	Efficiency Worst Best
	1	23.6	1373.00	23.6	23.6	2.48	2.46	100%	100%

Objectives Function Calculation Time (Second) Efficiency

				Worst	Best	Longest	Shortest	Worst	Best
	2	16.0	136.00	16.7	16.2	2.48	2.47	96%	99%
	3	5.5	311.00	5.5	5.5	2.46	2.45	100%	100%

Table 3 GASA results for 4 workstations 20 products

 3

		Linear Programming		Genetic Algorithm with Simulated Annealing (GASA)
	Problems	Optimal Solution	Calculation Time (Second)	Objectives Function Calculation Time (Second) Worst Best Longest Shortest	Efficiency Worst Best
	1	23.6	1373.00	23.6	23.6	54.09	41.72	100%	100%
	2	16.0	136.00	16.2	16.0	49.04	39.79	99%	100%
	3	5.5	311.00	5.5	5.5	43.80	33.20	100%	100%

Table 4 GAresult on industrial data Problem Type Linear Programming Genetic Algorithm Last Found Solution Calculation Time (Second) Objectives Function Calculation Time (Second) Efficiency

 4

				Worst	Best	Longest	Shortest	Worst	Best
	2013-07-09	673.0	3600.00	662.4	634.4	273.88	217.71	102%	106%
	2013-07-10	445.1	3600.00	454.2	433.6	306.04	232.79	98%	103%
	2013-07-11	597.2	3600.00	602.7	593.3	318.89	224.41	99%	101%
	2013-08-13	720.0	3600.00	725.7	707.1	285.77	195.46	99%	102%
	2013-08-14	1566.2	3600.00	1571.6	1451.1	329.46	249.31	100%	107%
	2013-08-20	630.5	3600.00	627.7	617.7	351.68	245.90	100%	102%
	2013-08-21	640.3	3600.00	661.7	623.6	336.14	249.52	97%	103%
	2013-08-22	1248.2	3600.00	1268.9	1237.0	370.20	267.34	98%	101%
	2013-08-23	632.4	3600.00	609.1	600.3	349.69	284.12	104%	105%

Table 5 SAresult big instances

 5

		Linear Programming			Simulated Annealing
	Problem Type	Last	Calculation	Objectives Function	Calculation Time (Second)	Efficiency
		Found	Time
		Solution	(Second)	Worst	Best	Longest	Shortest	Worst	Best
	2013-07-09	673.0	3600.00	760.5	677.9	17.22	17.07	87%	99%
	2013-07-10	445.1	3600.00	505.0	447.7	17.14	17.07	87%	99%
	2013-07-11	597.2	3600.00	647.7	606.8	17.09	17.03	92%	98%
	2013-08-13	720.0	3600.00	849.2	762.5	16.29	16.23	82%	94%
	2013-08-14	1566.2	3600.00	1744.1	1567.5	17.35	17.26	89%	100%
	2013-08-20	630.5	3600.00	704.8	634.6	17.45	17.40	88%	99%
	2013-08-21	640.3	3600.00	749.7	640.9	17.48	17.38	83%	100%
	2013-08-22	1248.2	3600.00	1446.4	1312.6	17.48	17.42	84%	95%
	2013-08-23	632.4	3600.00	672.5	636.2	17.46	17.39	94%	99%

Table 6 GASA result big instances Problem Type Linear Programming Genetic Algorithm with Simulated Annealing (GASA) Last Found Solution Calculation Time (Second) Objectives Function Calculation Time (Second) Efficiency

 6

				Worst	Best	Longest	Shortest	Worst	Best
	2013-07-09	673.0	3600.00	656.8	635.8	430.22	376.50	102%	106%
	2013-07-10	445.1	3600.00	446.2	437.1	450.65	387.72	100%	102%
	2013-07-11	597.2	3600.00	600.1	594.1	453.52	365.68	100%	101%
	2013-08-13	720.0	3600.00	717.7	710.9	421.68	371.96	100%	101%
	2013-08-14	1566.2	3600.00	1526.6	1460.7	471.57	381.66	103%	107%
	2013-08-20	630.5	3600.00	625.6	618.7	443.82	376.97	101%	102%
	2013-08-21	640.3	3600.00	629.4	623.0	486.59	401.36	102%	103%
	2013-08-22	1248.2	3600.00	1263.9	1239.0	468.89	404.55	99%	101%
	2013-08-23	632.4	3600.00	607.1	593.8	482.94	403.11	104%	106%

CIE44 & IMSS'14 Proceedings, 14-16October 2014, Istanbul / Turkey


S Tewari

P Roulleau

C Grenier

F Portier

A Cavanna

U Gennser

D Mailly

P Roche

Robust quantum coherence above the Fermi sea

HAL is

In this paper we present an experiment where we measured the quantum coherence of a quasiparticle injected at a well-defined energy above the Fermi sea into the edge states of the integer quantum Hall regime. Electrons are introduced in an electronic Mach-Zehnder interferometer after passing through a quantum dot that plays the role of an energy filter. Measurements show that above a threshold injection energy, the visibility of the quantum interferences is almost independent of the energy. This is true even for high energies, up to 130 μeV, well above the thermal energy of the measured sample. This result is in contradiction with our theoretical predictions, which instead predict a continuous decrease of the interference visibility with increasing energy. DOI: 10.1103/PhysRevB.93.035420

A new type of quantum device, relying on the onedimensional edge states of the quantum Hall regime, where electrons mimic the photon trajectory of a laser beam, has opened a route towards electron quantum optics and the manipulation of single electron excitations [1][2][3]. Pauli statistics and interactions provide new ingredients for the physics of the electrons which are not relevant for photons. For example, when electrons are injected above the Fermi sea, it is fundamental to understand how their phase coherence will be affected by the injection energy. We explore this issue by first using a quantum dot to inject the carriers at a controllable energy into an edge state. Then an electronic Mach-Zehnder interferometer is used to monitor the quantum coherence of the electronic quasiparticle. We find that above a certain threshold, the coherence is energy independent; it is even preserved at energies 50 times larger than the electronic temperature. This is remarkable, since from simple considerations based on Fermi's golden rule, one would expect that the relaxation rate increases with the injection energy, thus reducing quantum coherence. Indeed, our simulations using recent theories [4] predict a continuous trend of increasing relaxation. While the origin of this coherence robustness remains unidentified, it has a significant bearing for the implementation of quantum information encoded in electron trajectories [5][6][7][8].

The edge states used in this new type of device for electronic quantum optics are obtained by applying a high magnetic field perpendicular to a high mobility two-dimensional (2D) electron gas. When the number of electrons per quantum of flux (the filling factor) is an integer, the transport occurs through one-dimensional channels located at the edge of the sample. The electron motion in these wires is chiral: The electrons drift in one direction with a speed of the order of 10 4 -10 6 ms -1 [9,10], thus compensating the confining electric field with the Lorentz force. Much recent progress has been made in the understanding of decoherence mechanisms and the energy exchanges in the integer quantum Hall regime at a filling factor of 2. In this regime, there are two parallel, adjacent channels at each edge in the sample. The Coulomb interaction has been shown to play a key role in these systems: * patrice.roche@cea.fr (i) The quantum coherence at finite temperature is limited by the thermal charge noise of the environment [11,12]; (ii) there is an energy relaxation in out-of-equilibrium edge states and an energy transfer between the two edge states [13]; and (iii) the visibility of quantum interferences in electronic Mach-Zehnder interferometers exhibits a side lobe structure at finite bias, which is explained by a beating effect between a neutral and a charged excitation shared by the two edge states [4,[14][15][16][17]. Here, we inject a quasiparticle at a well-defined energy above the Fermi sea into an edge channel. We then explore to which extent such a single charge behaves as a free noninteracting particle in this interacting quantum system, thus directly probing the validity of the Landau's Fermi liquid picture [18,19]. While extensively studied in diffusive quantum conductors, this question has yet to be addressed experimentally for the case of one-dimensional chiral conductors. We note, however, that this is a key point for quantum information experiments using electrons transported through edge channels. Experimentally, we use a quantum dot (QD) as an energy filter, and a Mach-Zehnder interferometer to probe the injected quasiparticle's loss of phase coherence.

The QD consists of two quantum point contacts (QPCs) and a plunger gate (see Fig. 1). Figure 2(a) shows a 2D plot of the transmission probability through the dot as a function of the drain-source bias and the plunger gate voltage V P . The transmission probability is defined as dI T /dI 0 , where I T is the transmitted current through the QD and I 0 is the impinging current. From this measurement we deduce the lever arm of the plunger gate α = d /d(eV P ) = 1.46 × 10 -4 , where is the QD energy. As an example, Fig. 3 (top graph) displays the current I T as a function of V P . It is clear that I T is almost constant for energies varying between 0 and eV DS . This signals an absence of excited states for this particular tuning of the quantum dot. We deduced an electronic temperature of 31 ± 1 mK by fitting the resonance of the transmission probability [see Fig. 2(b)] as a function of V P with dI T /dI 0 ∝ cosh -2 (δ/2k B T ) [20], where δ = α(eV P -eV P 0 ) is the energy difference between the QD energy level and the Fermi level, and V P 0 is the plunger gate voltage which maximizes the conductance. Indeed, the intrinsic width of the energy filter is much smaller than the electronic temperature. In practice, the presence of excited states in the dot limited the energy range we explored to 120 μeV, close to the addition energy estimated to be of the order of 150 μeV [see the Supplemental Material (SM) [21]].

Once the QD is characterized, we measure the current through the whole system, constituted by the QD in series with the MZI, and operated in its optimal parameter regime for maximum visibility (see SM [21]). Probing the decoherence at finite energy is simply realized by measuring the relative amplitude of the current oscillations revealed by sweeping the central gate voltage V C of the Mach-Zehnder interferometer. A color plot of the current as a function of V C and V P is displayed in Fig. 3. This plot corresponds to experiments with five different drain-source voltages V DS . One would expect a current injection only for a range of V P values between 0 and V DS (see the inset of Fig. 1), and this is indeed what can be seen in the top graph of Fig. 3, where the window of nonzero current (around 25 pA) increases linearly with increasing source-drain voltage. Sweeping V C while V P is within this window reveals the quantum interferences. The interference pattern is clearly observable in Fig. 3. It is easily noticed that, within the window for nonzero current, the amplitude of the current oscillations appears to change very little, a clear indication that the quantum interferences, and hence the coherence, are more or less independent of the injection energy of the electrons.

To be more quantitative, we plot in Fig. 4 the visibility quantum coherence, or e.g., to the finite distance, 2.7 μm, between the energy filter and the MZI. Indeed, a short relaxation length ( 2.7 μm) would mean that the electrons injected above the Fermi sea are fully relaxed before reaching the entrance of the interferometer. However, this scenario is in contradiction both with recent experimental results and theories, as we will now show.

On the experimental side, experiments at a filling factor of 2 with a hot electron distribution in one channel have demonstrated that the inelastic length is of the order of 3 μm [13]. This length is of the order of the separation between the QD and the MZI, meaning that electrons are only partially relaxed. A second experimental evidence against fast relaxation is the Hong-Ou-Mandel (HOM) experiment with electrons injected one by one into edge states via quantum dots, which are excited with a periodic radio-frequency gate voltage matched to their addition energy [3]. The observed HOM dip, resulting from the collision of two electrons emerging from two separate sources, implies a coherent propagation of the injected electrons. This observation, for which electrons have been injected at energies higher than 1 meV, shows qualitatively that an important fraction of the injected electrons remains coherent, as observed in our experiment.

Furthermore, our quantitative model calculations show that the electrons reaching the MZI are only partially relaxed. We describe our system through an approach similar to Ref. [22], which considers the relaxation of an energy resolved single electron wave packet φ k 0 (x) ∝ e ik 0 x injected in the outer edge channel, and takes into account the coupling of the ν = 2 edge channels [4]. This approach has proved useful for the understanding of unexpected phenomena observed on out-of-equilibrium MZIs [17], and of the observed energy transfer between edge channels [13,23]. The coupling to the electromagnetic environment, responsible for the electronic decoherence, is described by an electrostatic coupling to a metallic electrode, itself connected to the ground by an arbitrary impedance at an arbitrary temperature. This idealized picture is expected to be valid for electrons having a typical injection energy ε 0 = v F k 0 , with their energy width. It is a criterion that is valid in our experiment, and as long as the dephasing environment can be considered as linear and in thermal equilibrium. Based on our previous experiments on phase coherence in the integer quantum Hall effect (IQHE) at ν = 2 [11,12], we assume that the electromagnetic environment of the interfering edge state consists of the co-propagating edge state, which replaces the metallic gate. Furthermore, as we probe the phase coherence of electrons injected at energies above the Fermi sea much higher than the thermal energy, we assume a zero temperature of the neighboring edge mode. The effect of its finite temperature is then simply taken into account in the finite visibility of Aharonov-Bohm (AB) interferences at zero energy. The two main outcomes of this approach to quasiparticle relaxation are the electron distribution δn(k,k 0 ,L) and the elastic scattering probability Z(k 0 ,L) after a finite propagation length L. The visibility is then obtained from δn(k,k 0 L)Z(2k 0 L)dk.

To include the relaxation on the length L QD separating the QD from the MZI, we computed δn(k,k 0 L QD ) and the full theoretical V Theo. ∝ δn(k,k 0 L QD )Z(2kL MZI )dk, where L MZI is the MZI arm length. The inclusion of this length shows that the electronic wave packet is partially (and not totally) relaxed before entering the MZI. The theoretical results are plotted in Fig. 4 for two different values of the coupling parameter θ measured in recent experiments [24,25]. In order to fit our data and the decrease of the visibility at low energy, we chose a drift velocity of 5×10 4 ms -1 . However, the general appearance of the curves does not qualitatively change for reasonable values of the coupling strength and of the drift velocity (see SM [21]).

As expected, but in contradiction with our observations, the theory does not lead to a saturation of the visibility at high energy, but to a continuous decrease, albeit with a slope less pronounced than at low energy. This continuous decrease can be easily understood: For a short range interaction, the environment can absorb any energy; hence, the higher the injection energy, the larger is the probability for energy losses. This argument is in fact quite general and should apply regardless of the microscopic details of the model. Hence, although the other existing models [4,[26][27][28] have not treated explicitly the decoherence of electrons injected in a narrow energy window above the Fermi level, we do not expect they would yield predictions qualitatively different from ours, as they also start from a short range Coulomb interaction between co-propagating edge channels. Last, a recent theoretical paper [29] based on short range interaction reaches the same conclusions. Our experimental observation of a saturation of the high energy visibility suggests the existence of another not yet identified energy scale in the system, or suggests the presence of an additional relaxation mechanism which remains to be identified. A somewhat similar effect has been probed at lower filling factors and much higher energies [30].

The robustness of the quantum coherence, which we have demonstrated here, increases the attraction of electron guns for quantum optics type experiments. At the same time, our results raise questions regarding our understanding of excitations in the integer quantum Hall regime. To further elucidate the relaxation in edge states, the next step in this experiment will be to implement a second quantum dot after the MZI, to make the spectroscopy of the injected electron after its propagation through the interferometer. P.R. would like to thank Eugene Sukhorukov for stimulating discussions. This work has been supported by the French ANR Contract No. 11-BS04-022-01 IQHAR.

FIG. 1 .

 1 FIG. 1. Measured device. Colorized scanning electron microscope view of the sample before the final fabrication step in which the gates and ohmic contacts are connected to larger connecting pads. The sample consists of a quantum dot in series with an electronic Mach-Zehnder interferometer (MZI). The quantum point contacts (QPCs) and gates controlling the dot are in red while those controlling the MZI are in yellow. QPC1 and QPC2 serve as the two beam splitters of the MZI. Note the small ohmic contacts in blue connected to the ground, which prevents spurious quantum interferences inside the interferometer. The chiral trajectory of the excitations in the outer edge state is schematically represented in blue, and the quantum interference takes place between the trajectories (a) and (b). The mesa is 1.2 μm wide, each arm of the Mach-Zehnder interferometer is 7.2 μm long, and the distance between the dot (the energy filter) and the MZI is 2.8 μm. V P is the potential applied on the plunger gate of the dot, allowing a control of the relative position of the energy levels of the dot compared to the Fermi level. The quantum interferences in the MZI are revealed by sweeping the gate voltage V C , which modifies the trajectory length difference between (a) and (b).

FIG. 2 .

 2 FIG. 2. (a) Color plot of the transmission probability through the quantum dot as a function of the plunger voltage V P and the drain source voltage V DS . (b) Transmission probability as a function of V P at zero bias. The continuous line is a fit of this resonance.

I 4 FIG. 3 .

 43 FIG. 3. Current as a function of the plunger gate V P and the central gate voltage V C of the MZI for different values of the drain source bias V DS . The top graph, representing the current as a function of V P , corresponding to a 2D plot cut through line B. The right graph shows the interference pattern as a function of V C corresponding to a cut through line A.

  FIG.4. Visibility of quantum interferences as a function of the energy for five different bias voltages. The vertical lines represent the maximum energy achievable for each bias. The solid and the dashed lines show the calculated visibility for interaction parameters θ = π/2 and θ = π/3, respectively.