Claire Olivier

Matthieu Nolwenn

Corentine Lara

Cette thèse n'aurait pas été possible sans les discussions, le soutient et l'aide d'un grand nombre de personnes que je me dois de remercier ! Tout d'abord, il y a Frank Fournel qui m'a bien vendu le sujet de la thèse, alors que je l'avais au départ contacté pour un autre sujet. Merci également Frank d'avoir suivi cette thèse de près tout en me laissant suffisamment d'autonomie pour que cette thèse soit bien la mienne. Ensuite bien sûr vient Élisabeth Charlaix qui a dirigé cette thèse et surtout qui a apporté une expertise impressionnante au niveau de la mécanique des fluides. Merci encore pour toutes ces journées passées à réfléchir ensemble sur le sujet et la simulation.

Je tiens aussi à remercier tous les membres du jury de ma thèse, pour commencer, son président, François Rieutord, dont les travaux ont été un excellent point de départ et avec qui j'ai eu plaisir à échanger. Ensuite mes deux rapporteurs, Étienne Barthel et Matteo Ciccotti, qui ont lu mon manuscrit avec attention et qui m'en ont fait des retours constructifs et mélioratifs ! Merci également à Roy Knechtel d'avoir fait le déplacement pour la soutenance, mais également pour m'avoir invité à participer à la conférence WaferBond 2013 à Stockholm. Et bien sûr encore une fois merci à Frank et Élisabeth qui étaient aussi membres de mon jury.

Ces remerciements continus avec toutes les personnes des laboratoires du DTSI et du DCOS. En premier lieu, Fabrice Geiger, Chrystel Deguet et Véronique Carron pour m'avoir accueilli au sein du DTSI/SSURF/LSJ et avoir permis à cette thèse de se dérouler dans d'excellentes conditions, tant du cotév du matériel et des infrastructures que du coté de l'ambiance et de la vie au laboratoire. Puis tous les gens avec qui j'ai pu travailler et échanger, merci à Hubert, Vincent, Loïc, Momo, Christophe, Laurent, Maman, Floch, Del, Caro, ... Sans oublier tous les thésards, alternants et stagiaires et autres avec qui j'ai également pu travailler, mais surtout partager de bon moments de détente, tant lors des pauses au bureau que lors de sorties en soirées et week-end. Je pense notamment à Xavier, Raphaël, Ludo, Marie, Toshi, Lyvia, Sylvain, Marwann, mais je ne peux pas citer tout le monde car la liste serait bien trop longue ! Il y a bien sûr aussi toutes les personnes extérieures au CEA qui m'ont soutenu, parfois même sans le savoir, durant ces trois ans. D'abord les membres de l'association Excalibur Dauphiné, mais aussi tous les scouts avec qui j'ai pu partager de bons moments de vie et ma famille et mes amis, bien entendu. Merci donc à Jeff, Jean d'arc et Val, Emrys, le Caro et tous les autres membre d'XK, à Delphine, Nicolas,

Introduction

Unlike a lot of PhD thesis, the title of this thesis, Direct bonding of patterned surfaces, is quite short. However like all of them it needs to be explained. Direct bonding is a process by which two sufficiently flat and clean surfaces can bond to each other without any added adhesive layer. In theory any materials can be bonded by direct bonding, however the flatness and cleanness requirement are such that in practice this process is often limited to materials used for micro-electronic devices in a clean room environment. Thus the main applications and studies of direct bonding are related to the micro-electronic industry.

The patterned surfaces in the title are also related to the micro-electronic industry, especially for the fabrication of Micro-Electro-Mechanical Systems (MEMS), and simply are silicon wafers with small cavities etched on their surface. Those silicon wafers with cavities are then bonded to each other or to a plain wafer to make MEMS devices (Fig. 1). The fabrication of these devices is expensive and it would be quite useful to have guideline when designing new devices to know in advance if direct bonding will be possible. This thesis aims to study the direct bonding dynamics over patterned surfaces. By building a simulation model of the bonding wave propagation we expect to be able to give guideline regarding the design of cavities.

Inlet Outlet Chamber

Valves

Piezoelectric actuator

Pump Diaphragm The bonding of two substrates can be described as follow: The two wafers are put on top of one another in a standard clean room environment. A thin air layer (≈ 30 µm) is trapped between the two wafers forming an air cushion. To initiate the bonding a small pressure is applied on the top wafer to bring it into contact with the bottom wafer. This creates a first bonded area that will propagate over the entire wafer area by expelling the trapped air. The frontier between the bonded and unbonded area will be called the bonding front or the bonding wave. The propagation of the bonding wave is observed by looking through the silicon wafers with an infrared (IR) camera (see Fig. 2), indeed silicon is transparent to infrared light with wavelength above 1 µm. Other bonding processes exist in the micro-electronic industry, such as anodic bonding, eutectic bonding, thermocompression bonding, etc. The main advantages of direct bonding over those are mainly due to its bonding interface properties. First of all the position of the interface is well defined because of the absence of a bonding layer. This interface can also easily have good optical or electrical transparency. It can also withstand high temperature treatment and the bonding process in itself allows quite a good throughput. But those advantages have an important counterweight: the required surface quality. The use of complex (and thus expensive) processes and dedicated (clean) tools is required to obtain or keep the appropriate surface quality.

Following this introduction the manuscript is divided in five chapters. Chapter I quickly describes the state of the art of direct bonding. Chapter II is a description of the simulation model build for this study. Chapter III presents all the experimental tools and processes used for the bonding experiments. Then starts the study of the dynamics of the direct bonding wave in Chapter IV with an experimental validation of the simulation results. Finally an in depth study of the energy of adhesion E a responsible of the bonding wave propagation is lead in Chapter V. At the end of the manuscript the reader will find a collection of appendix, most of them describing the different software used or made during this thesis, which will be a most useful resource for any potential future user.

Nomenclature Variables and Constants

E

Young's modulus Chapter I

E e = E

State of the Art Introduction

This thesis main topic is direct bonding which is a really well studied topic, especially since the mid-1980s. Firstly, the state of the art of direct bonding will be presented through a look at the history, main uses and experimental parameters of direct bonding.

To quantify the quality of a direct bonding two energies can be measured. First the bonding energy (E c ), which is the energy required to separate two bonded wafer, and the adhesion energy (E a ), which is the energy that drives the bonding wave. The section on the measurements of those energies will describe them in more depth. Then a look at the existing simulations and studies of the bonding wave will show the need of our own simulation model. Finally we will see that the existing studies of cavities direct bonding focus on the quality of the fluid seal that they offer.

I.1 Direct bonding I.1.1 History

As described in the general introduction, direct bonding is the bonding of two sufficiently flat and clean surfaces without any added adhesive layer. The first mention of this phenomenon is from Desagulier in 1792 [START_REF] Dowson | History of tribology[END_REF] who bonded two sphere of lead to each other, with the help of an applied external pressure. Then Lord Rayleigh published in 1936 [START_REF] Rayleigh | A Study of Glass Surfaces in Optical Contact[END_REF] a study on glass surfaces bonding to each other at room temperature, a phenomenon known by the optician and used to bond optical elements with good transparency. In the 1980s two different groups, Shimbo and co-workers from Toshiba [START_REF] Shimbo | Silicon-to-silicon direct bonding method[END_REF] and Lasky from IBM [START_REF] Lasky | Wafer bonding for silicon-on-insulator technologies[END_REF] investigated the direct bonding of silicon wafers. At this point the micro-electronic industry knew a rapid growth and many groups investigated direct bonding, also called wafer bonding by some research groups. Stengl [START_REF] Stengl | A model for the silicon wafer bonding process[END_REF] is the first to present a mechanism of the direct bonding process, using hydrogen bonds of water adsorbed on the surfaces to be bonded. In 1999 Tong and Gösele [START_REF] Tong | SemiConductor Wafer Bonding: Science and Technology[END_REF] published a book that gives a really good view of what was known about direct bonding. A few years later an important review article by Christiansen et al. [START_REF] Christiansen | Wafer direct bonding: From advanced substrate engineering to future applications in micro/nanoelectronics[END_REF] gives an updated view of the wafer direct bonding industry, research and future applications.

I.1.2 Main uses

Direct bonding has many usage in the micro-electronic field. It can be used to make advanced substrate like Silicon on Insulator (SOI), first by back etching [START_REF] Lasky | Wafer bonding for silicon-on-insulator technologies[END_REF] and then using the Smart Cut™1 process [START_REF] Bruel | Silicon on insulator material technology[END_REF] which saves a lot of material. The SOI fabrication process is actually no more than a layer transfer process, which is another important application of the direct bonding. Pushing further this layer transfer capability another important application is the stacking of active layers to form 3D devices driven by the More than Moore technologies. This layer stacking is even more interesting with direct bonding due to the possibility to bond heterostructures at low temperature [START_REF] Moriceau | Low temperature direct bonding: An attractive technique for heterostructures build-up[END_REF] and to provide electrical interconnection [START_REF] Beilliard | Chip to wafer copper direct bonding electrical characterization and thermal cycling[END_REF]. By bonding a wafer with cavities to another one or to a flat wafer, direct bonding can be used to make MEMS such as pressure sensor [START_REF] Petersen | Silicon fusion bonding for pressure sensors[END_REF] or to make hermetic package for any 3D MEMS [START_REF] Tan | Low temperature wafer-level bonding for hermetic packaging of 3D microsystems[END_REF].

I.1.3 Experimental parameters

When bonding two wafers many experimental parameters must be considered. They were extensively studied over the years and while the effects of some are now quite defined the influence of some others is just observed or suspected. We can separate them in two groups, the parameters of the wafers and their surfaces and the environmental parameters.

Wafers parameters:

Bow, waviness and roughness: As already mentioned, to perform the direct bonding of two surfaces those must be flat enough. This requirement can be described by three distinct parameters, each concerning a different scale (Fig.I.1): at the largest scale, the bow (Fig. I.1a) of the wafer is the height difference between the center and the edge of the wafer. It can cause the bonding to be impossible or to stop before completion. Turner's work [START_REF] Turner | Modeling of direct wafer bonding: Effect of wafer bow and etch patterns[END_REF] shows that to bond a 100 mm diameter silicon wafer with a bow of 40 µm to a flat identical silicon wafer, both 750 µm thick, a minimal adhesion energy (E a ) of 12 mJ m -2 is required. The waviness (Fig. I.1b) occurs at an intermediate scale and is not well studied, thus no numerical limits has been found. Finally the roughness (Fig. I.1c), at the smallest scale, is more studied. It is usually measured using an atomic force microscope (AFM) and can be monitored by fine polishing or chemical etching. The work of Rayssac [START_REF] Rayssac | Etude du collage par adhesion moleculaire hydrophile : application au controle de l'énergie de collage[END_REF] gives a limit of acceptable roughness for the bonding of two 100 mm diameter oxidized silicon wafers of 5 Å of root mean square roughness (R rms ) as measured by AFM on a 5 µm square. Miki [START_REF] Miki | Effect of nanoscale surface roughness on the bonding energy of direct-bonded silicon wafers[END_REF] prefers to use the bearing ratio to characterize the roughness of the surfaces and find a linear correlation between the bearing ratio and the bonding energy (E c ) of the samples.

Mechanical properties:

If bonding an easily deformable wafer, the energy required to deform the wafer into close contact will be less. To quantify the rigidity of a wafers the parameter D, called the flexural rigidity is commonly used [START_REF] Landau | Theory of Elasticity[END_REF][START_REF] Rieutord | Dynamics of a Bonding Front[END_REF][START_REF] Navarro | Direct silicon bonding dynamics: A coupled fluid/structure analysis[END_REF][START_REF] Timoshenko | Theory of Plates and Shells[END_REF]: Where E e is the reduced Young's modulus, E the Young's modulus, ν the Poisson's ratio and t w the thickness of the wafer.

D = E e • t w
D has an influence when deforming the wafers which can be because of their bow or waviness as seen above. But even with perfectly flat wafers, when the bonding propagates while expelling the trapped air cushion, the wafers are deformed and thus the parameter D is important even when the bow or waviness are not.

Surface contaminations: Two kind of contaminations are usually present on the wafer surfaces, the organic contamination and the particles contamination [START_REF] Tong | SemiConductor Wafer Bonding: Science and Technology[END_REF]. To remove those, wet chemical cleaning, chemical-mechanical polishing (CMP) or dry etching processes can be used. Both type of contamination can create defects at the bonding interface. For example the presence of a particle of diameter H can create a unbonded area of diameter 2R d (Fig. I.2) as follow (adapted from [START_REF] Tong | SemiConductor Wafer Bonding: Science and Technology[END_REF]):

R d = 1 3 • E e • t w 3 E a 1 /4 • H 1 /2 (I.2)
Where E e is the reduced Young's modulus of the wafers, t w their thickness and E a the adhesion energy. For a standard pair of wafers (See table II.2 for the parameters value) bonded with E a = 40 mJ m -2 around a particle of diameter H = 1 µm the predicted unbonded area diameter is 2R d = 0.97 cm. Eq. I.2 is true for R d > 2 • t w , for particles inducing smaller unbonded area another equation should be used, which gives an unbonded radius R d much smaller, see section 3.2 at pages 42-43 in [START_REF] Tong | SemiConductor Wafer Bonding: Science and Technology[END_REF].

Surface treatment:

In direct bonding the bonding is driven by the links formed between the two surfaces. The chemical nature of the surface is thus important. It has often been observed that the presence of adsorbed water molecule at the surface of the wafers [START_REF] Stengl | A model for the silicon wafer bonding process[END_REF] helps to bond the wafers, so that a hydrophilic surface will bond better than a hydrophobic surface, at equivalent surfaces roughness. With an appropriate surface cleaning, Si-OH groups can be formed at the surface and gives a good hydrophilic surface [START_REF] Tong | SemiConductor Wafer Bonding: Science and Technology[END_REF].

Environmental parameters:

Temperature: Because of the hight thermal conductivity (149 W m -1 K -1 ) of silicon and the relatively thin thickness of the wafers, the temperature of the bonding atmosphere or of a heated chuck is also the temperature of the bonded wafers. The elevation of temperature can influence the viscosity of the air, the surface chemistry and the mechanical properties of the wafers. So its effect on the bonding is quite difficult to predict.

Pressure, viscosity and relative humidity:

The air properties influence mainly the propagation of the bonding front. Indeed the air is expelled by the bonding front propagation and thus acts as a brake. It is quite straightforward to understand that an increased viscosity will increase the brake effect and slow the bonding front propagation. The influence of those parameters is more discussed in section I.3.

Direct bonding is a process used to bond wafers without any added adhesive layer. It requires surfaces with good flatness, waviness and roughness and is influenced by many other parameters such as the flexural rigidity D of the wafers, the cleanness of the surfaces and the properties of the atmosphere in which it takes place.

I.2 Measurements of bonding and adhesion energies

To characterize a direct bonding two different energies can be used. To measure them, mechanical models are used to extract the energy value (by unit area), often from a direct measurement of an unbonded length. Because it was first used, the bonding energy (E c ) is still more widely used than the adhesion energy (E a ) which is more relevant to this work.

The bonding energy:

The bonding energy (E c ) is the energy required to separate two bonded wafers. This energy is widely used to characterize the quality of a direct bonding. To measure E c the double cantilever beam (DCB) under prescribed displacement test is most often used. It is also well known as the Maszara method as he first used the Gillis and Gilman work for direct silicon wafer bonding energy measurement [START_REF] Maszara | Bonding of silicon wafers for silicon-on-insulator[END_REF]. For this measure a thin blade (usually of a few hundreds micrometers) is inserted between the bonded wafers. This insertion leads to the debonding of the wafers over a length L. Measuring this length (usually by infrared observation through the wafers, see Fig. I.3) the E c of two bonded wafers can be calculated with [START_REF] Gösele | What determines the lateral bonding speed in silicon wafer bonding?[END_REF]:

E c = 3 8
E e1 t w1 3 • E e2 t w2 3 E e1 t w1 3 + E e2 t w2 3 • t b 2 L 4 (I.3)
With E e = E/(1 -ν 2 ) the reduced Young's modulus of the wafers, E the Young's modulus, ν the Poisson's ratio, t w their thickness, with index 1 and 2 each referring to a different wafer and t b the thickness of the blade. This method is so often used thanks to its much more practical and easier setup compared to the other possible methods, such as the tensile test method, the chevron method, etc. [START_REF] Vallin | Adhesion quantification methods for wafer bonding[END_REF]. However the Maszara method yield a high error in the results, mainly due to the wafer geometry [START_REF] Turner | Accurate characterization of wafer bond toughness with the double cantilever specimen[END_REF] and the debonding length error which is elevated at the power four in Eq. I.3. The use of beams cut in a wafer helps to reduce the geometry error and the debonding length needs to be correctly calculated, see Subsec. III.3.2, to further decrease the final error. Moreover the silicon being an anisotropic material the work from El-Zein et al. [START_REF] El-Zein | Evaluation of GIC of a DCB specimen using an anisotropic solution[END_REF] can be used to further improve the calculation accuracy. Another thing to take into account when measuring E c is the humidity of the environment. Indeed Fournel et al. [START_REF] Fournel | Measurement of bonding energy in an anhydrous nitrogen atmosphere and its application to silicon direct bonding technology[END_REF] remind that the water molecule usually present in the debonding atmosphere lead to a stress corrosion of the bonding link and thus influence the bonding energy measurement interpretation. In this work an experimental setup is also described to avoid any stress corrosion by doing the energy measurements in an anhydrous nitrogen atmosphere (See Subsec. III.1.5).

The adhesion energy:

The second energy that one can measure to characterize a direct bonding is the adhesion energy (E a ), which is the energy released by the bonding of the two wafers. In other words E a is the energy that actually drives the bonding front forward. The asymmetry between E c and E a was hinted at by Miki in 2003 [START_REF] Miki | Effect of nanoscale surface roughness on the bonding energy of direct-bonded silicon wafers[END_REF], and Grierson [START_REF] Grierson | Characterization of Hysteresis of Surface Energy in Room-Temperature Direct Bonding Processes[END_REF] studied the hysteresis between the two energy at different relative humidity (RH). His setup consists of a micro-sized beam whose distance to a flat silicon surface is monitored. When bringing the beam closer to the support, the bonding front advances and the measure of the unbonded length lead to the value of E a . Then the beam is taken away from the support to measure the E c which is higher than E a . Navarro [START_REF] Navarro | Direct Wafer Bonding Dynamics[END_REF] also has an experiment to measure this E a , involving a bonding over an obstacle of known height and the use of a mechanical deformation model to deduce the E a from the deformation profile over the obstacle. In 2006 while studying the bonding front propagation of bowed wafers, Turner [START_REF] Turner | Mechanics of direct wafer bonding[END_REF] reported that E a increases over time, observing a short time fast increase and then a slow increase over a long time period.

Apart from those studies of the E a of direct bonding, more general studies of the adhesion of solid exist, such as the work from Charlaix [START_REF] Charlaix | Adhesion forces between wetted solid surfaces[END_REF] which study the adhesion forces between wetted solid surfaces with a force apparatus. Another interesting paper from Bocquet [START_REF] Bocquet | Moisture-induced ageing in granular media and the kinetics of capillary condensation[END_REF] the kinetics of capillary condensation and its effect on granular media cohesion.

Two different energies exist in wafer bonding characterization, E c and E a . E a is the most important energy for this work as it is the energy that drives the bonding front forward. So far the methods used to measure it are not so simple. A long-time evolution of E a has been reported.

I.3 Simulations and studies of the direct bonding wave

The observation and study of the direct bonding wave, i.e. the propagation of the bonding front, is important to understand the direct bonding mechanism and its limitations. The relevant energy here is the E a , however the distinction between E c and E a was not made in the earlier work presented below. For coherence we will replace any mention of E c by E a when applicable.

Experimental observations:

In the mid-1990s two groups studied the dependence of the bonding wave velocity with some experimental parameters. Gösele et al. [START_REF] Gösele | What determines the lateral bonding speed in silicon wafer bonding?[END_REF] mainly studied the influence of the pressure and nature of the gas between the wafers. They observed that the bonding speed decreases when the pressure increases or when the viscosity increases. They also tested two different wafer thickness (380 µm and 905 µm), without any apparent effect on the bonding velocity. However Bengtsson et al. [START_REF] Bengtsson | The influence of wafer dimensions on the contact wave velocity in silicon wafer bonding[END_REF] specifically studied wafer thickness influence, with 100 mm diameter wafers ranging from 270 µm to 5000 µm in thickness. They observed a decreasing bonding speed with increasing thickness and the thickest wafer did not bond.

1D models: A 1D model of the propagation of the bonding wave was proposed by Rieutord [START_REF] Rieutord | Dynamics of a Bonding Front[END_REF]. This model assume a steady state system of two symmetrical wafers being brought together by an adhesion energy and pushing away a fluid along a line (See Fig. I.4). It gives the following dependances for the bonding speed:

v = A 3 /4 9 • E a 5 /4 η • t w 3 /4 • Λ 1 /2 E e 1 /4 (I.4)
With v the velocity of the bonding wave, A a constant close to 1, E a the adhesion energy, η the fluid dynamic viscosity, t w the thickness of the wafers, Λ the mean free path of the fluid particles (cut-off distance) and E e the reduced Young's modulus. Very recently the work from Navarro [START_REF] Navarro | Direct silicon bonding dynamics: A coupled fluid/structure analysis[END_REF] also proposed a 1D model, with a similar system of an air layer trapped between two wafers but with a transient study showing different step of the bonding propagation and the evolution of the wafer profile over time. 2D models: Two 2D models exist in the litterature, the first by Turner [START_REF] Turner | Modeling of direct wafer bonding: Effect of wafer bow and etch patterns[END_REF][START_REF] Turner | Mechanics of direct wafer bonding[END_REF] and the second by Kubair [32,[START_REF] Kubair | Cohesive zone model for direct silicon wafer bonding[END_REF][START_REF] Kubair | Multiscale mechanics modeling of direct silicon wafer bonding[END_REF]. Both models are static model which means that they provide only an energy limitation for the bonding or a final state and that they do not take into account the fluid flow between the wafers. Turner works on wafer with bow and succeeds to accurately determine, prior to bonding, where the bonding wave will stop due to the lack of energy to further deform the bowed wafers into contact. His study also uses a pattern of etched cavities to decrease the adhesion energy when bonding bowed wafers. The model from Kubair use a Cohesive Zone Model to tackle the problems of local effects of patterning, surface roughness and waviness and also the asymmetry between E c and E a . For that a system of two wafers bonded by a water interlayer is used which can be seen as a direct bonding with an especially thick layer of adsorbed water.

1D dynamic models exist, and 2D static models too. But the absence of 2D dynamic model of the bonding wave propagation is to be noted.

I.4 Direct bonding and cavities

As explained previously many applications exist with the direct bonding of patterned wafers [START_REF] Christiansen | Wafer direct bonding: From advanced substrate engineering to future applications in micro/nanoelectronics[END_REF][START_REF] Schmidt | Wafer-to-wafer bonding for microstructure formation[END_REF]. However few study the direct bonding process over cavities. Turner uses cavities to reduce the bonding area and Mack [START_REF] Mack | Gas tightness of cavities sealed by silicon wafer bonding[END_REF][START_REF] Mack | Analysis of Bonding-Related Gas Enclosure in Micromachined Cavities Sealed by Silicon Wafer Bonding[END_REF] investigates the sealing properties of such a direct bonding. But nothing was found that study the dynamics of the bonding wave over cavities, nor about limitations to the design of cavities that could arise.

Conclusion

Through the presentation of the direct bonding state of the art it was show that seemingly simple process of bonding two wafers is complex to fully monitor. Indeed the high surface quality requirement is only a part of all the parameters that influence the direct bonding. Then the characterization methods of the direct bonding quality were presented. And again a seemingly simple method, the insertion of a blade between two bonded wafers, appeared more complicated, with influences from the geometry and properties of the materials, the atmosphere RH, ... The important difference between E c and E a was also explained, and it was noted that E a is more relevant to this thesis. Finally the presentation of the existing 1D dynamic models and 2D static models of direct bonding showed the need to develop a dedicated 2D dynamic model for the study of the direct bonding of patterned surfaces.

Introduction

In this work, simulation was first looked at as a tool to foresee the results of the direct bonding dynamic of patterned surfaces. Indeed designing a pattern for a specific purpose, then making it into a lithographic mask to test it out is quite expensive. It would be quite an edge to have a tool allowing us to predict if a particular patterned surface will correctly bond to another one. But as seen on section I.3, there is no existing tool available for our specific purpose. For building such a tool, a model of the direct bonding of surfaces is needed. The making of this model lead to many questions, and choices had to be made to keep a practical working model. This chapter will first describe the simple model we have decided to use before explaining how we used the commercially available Comsol ® software to turn this model into a simulation. Then variations to the simple model will be presented.

II.1 Description of the simple model

When making a model for simulation purposes a number of steps are to be followed. First, after clearly describing the physical system to model, the number of space dimensions and the geometry of the model are set. Then the physics of the model need to be defined, along with the equations that represent them. Finally the model needs an initial state and boundary conditions to be complete and usable in a simulation.

II.1.1 The physical system

As mentioned previously, the bonding study of patterned surfaces is the study of the direct bonding of a patterned Si wafer to a flat Si wafer. Our base model will just be of the bonding of two flat Si wafers cut into beams, the cavities of the pattern will be added later on (See II.3.4). The bottom beam is perfectly flat and fixed. The upper beam will be deformed during the bonding and will force the air between the two beams to flow out (Fig. II.1). The bonding wave will be moving from the left side to the right side of the system.

Top beam

Bottom beam

Air flow

II.1.2 Geometry

For our model a two-dimension space (2D) was required to be able to work on edge effects, which are due to the edges of the pattern or just the edges of the Silicon wafers to be bonded. Thus we consider a plane (x,y) intercepting the top flat surface of the bottom beam and a variable z is used to describe the position of the top beam. The geometry of the model is a simple rectangle of width w and length 2R (Fig. II.2). The base values used are w = 2 cm and R = 100 mm, which describes a beam of length equal to the diameter of a 200 mm wafer. This simple geometry does not consider the moving edge of the bonding wave, it's position will be known by post-treating the results of the following simulation (See II.2.3). It can also be noted that this simple geometry as several symmetry axis. The one cutting the beam in two and passing by the middle of the edges 1 and 4 will be of particular interest later on.

II.1.3 Physics

This simple model is animated by two physics: the mechanical deformation of the upper beam and the compression of the fluid between the two beams. Each of these physics must be written as an equation to turn the model into a simulation. We will also need equations describing the main forces applied to the beams: the action of gravity on the floating upper beam and the adhesion force between the beams.

The mechanical deformation:

The top beam is considered as a thin plate as its thickness is small compared to its width and length (t w w and t w 2R). The deformation applied on this thin plate are considered to be small deflection, i.e. deflection smaller than the thickness of the plate (t b t w ). Thus we can use the equation of equilibrium of a thin plate from [START_REF] Landau | Theory of Elasticity[END_REF] to describe its small deformation:

D∆ 2 z = F (II.1)
Here z is the altitude of the plate, D = Etw 3 12(1-ν 2 ) its flexural rigidity, with E the Young's modulus of the plate, ν it's Poisson's coefficient and t w its thickness. The Laplace operator ∆ is the two dimensional operator ∆ = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 , applied twice here as it is squared. The F on the right side of the equation is the sum of the forces applied on the plate by unit area. Considering the aforementioned physics we can write:

F = P visc + G + M att (II.2)
Where P visc is the force due to the fluid relative pressure between the plates, G the gravity force by unit area and M att the attraction force between the plates.

The pressure of the fluid: When the two plates come closer to each other they must expel the fluid trapped between them to do so. Thus the fluid applies a resistive action against the downward movement of the upper plate. This resistance can be expressed as the relative pressure of the fluid between the plates. This relative pressure can be expressed as [START_REF] Hocking | The effect of slip on the motion of a sphere close to a wall and of two adjacent spheres[END_REF]:

div(z 3 ∇P visc ) = 12η ∂z ∂t (II.3)
Where ∇ is the two dimensional gradient operator and η is the fluid dynamic viscosity. This equation diverge when z goes to 0, so a cut-off distance will be used to make our problem solvable. This cut off distance will be called z 0 and will be of the same order of magnitude as the molecular mean free path of the considered fluid. Indeed the equation II.3 is a fluid equation and it is reasonable to stop taking into account the fluid at a scale at which it is better described as free atoms moving around independently.

The gravity: The gravity simply pulls the upper wafer down as follow:

G = -ρgh (II.4)
With ρ the density of the plate, h it's thickness and g the gravitational constant.

Tha attraction force: As explained earlier in I.1, the direct bonding attraction force is viewed as Van der Waals force. But it has also been shown that the roughness of the bonded surfaces plays a role in this bonding. That and the fact that a numerical expression of the Van der Waals forces is highly repulsive (which is a problem for the convergence of the simulation solutions) lead us to consider our own model for the attraction force. A simple well of potential energy E p was considered and derived to obtain a conservative force M att . A Gaussian-like well was considered:

E p = -E a exp -(z -z 0 ) 2 β 2 (II.5)
Where the coefficient β is proportional to the well width (the width at half-height is about 1.68 × β) and is thus related to the force range. The cut-off distance z 0 , introduced above, is also the equilibrium distance of the force, where neither repulsion nor attraction occurs. Finally E a is the adhesion energy. Indeed the work (W ) of the conservative force, derived from this well, needed to bring two separate plates from infinity to the equilibrium distance is

W separated→bonded = E p(separated) -E p(bonded) = E a .
The force M att is thus written as:

M att = -2 • E a (z -z 0 ) β 2 exp -(z -z 0 ) 2 β 2 (II.6)
The Symmetry: One can note that the symmetry axis previously mentioned is also an axis of symmetry for every physics when the bonding wave propagate from one of the short edge as draw on Fig. II.2. This means that during the simulation we will be able to solve half our model and then get the full solution by simply applying a symmetry transformation to the partial solution.

Equations system: Those four equations give us a system of two equations with two variables, z and P visc , to solve. This system is of the fourth order, due to the squared Laplace operator and is non linear:

       D∆ 2 z = P visc + G + M att div(z 3 ∇P visc ) = 12η ∂z ∂t (II.7)
From now on, the following notation will be used for the variables and their derivative:

z → u ; P visc → v ; ∂ 2 u ∂x 2 → P ; ∂ 2 u ∂y 2 → Q
Moreover the time and space derivative form of every variable will be noted using a subscript letter as for example:

∂u ∂t → u t ; ∂ 2 u ∂x 2 → u xx → P ; ∂ 2 u ∂y∂x → u yx

II.1.4 Initial state

To run a simulation with our model an initial state is needed. This initial state needs to describe the position and speed of the upper plate and the value and time derivative of the relative pressure of the fluid between the plates. In other words we must describe z, ∂z ∂t , P visc and ∂P visc ∂t at t = 0 for the whole geometry. An unmoving initial state is considered, so both time derivatives are nulls and the system is at an equilibrium so the relative pressure is also null. Thus only the initial position of the upper plate is left to be described, and will be taken as a profile along the x dimension, while being flat along the y dimensions. Three initial states are considered as presented on Fig. II.4. The free fall and stylus pressure initial state have high gradient at the first contact point initiation, so are not really practical to use in a simulation. So for the simple model the retracting blade initial state will be used, with a simple second order profile:

u = t b 4R 2 x 2 + z 0 (II.8)
Where t b is the thickness of the retracting blade. This 2 nd order profile was used instead of the equilibrium 3 rd order profile given by Timoshenko [START_REF] Timoshenko | Theory of Plates and Shells[END_REF]:

u = t b 2 2 - 3 * (2R -x) L + 2R -x L 3 + z 0 (II.9)
with L the unbonded distance due to the inserted blade. This unbonded distance can be calculated with Eq. I.3, using E a to obtain the closing L:

E a = 3 8
E e1 t w1 3 • E e2 t w2 3 E e1 t w1 3 + E e2 t w2 3 • t b 2 L 4 (II.10)
this equation is then simplified to our case of a single moving wafer, thus taking the fixed wafer of infinite thickness (t w2 >> t w1 ) giving:

E a = 3 8
E e1 t w1 3 • t b 2 L 4 = 9 2 D • t b 2 L 4 (II.11)
Which gives the expression of L:

L = 9 2 1 /4 t b 1 /2 • D 1 /4 E a 1 /4
(II.12)

To have L = 2R a blade of thickness t b = 1.6 mm is needed but a blade this thick is out of our hypothesis range of small deflection. So the simple 2 nd order profile was used to have a full unbonded length whatever the thickness of the blade. This initial profile can be seen as a step just before the start of bonding. The use of this simpler profile will be justified in the result Subsec. IV.2.2.

Stylus pressure Free fall

Retracting blade

II.1.5 Boundary conditions

The last things that we need are boundary conditions. Our geometry has 4 boundaries, numbered from 1 to 4 as seen on Fig. II.2. Using the symmetry to reduce our geometry the edge number 3 become the axis of symmetry with the edges 1 and 4 with a length of w/2. The system calls for mechanical and fluid boundary conditions for each edge. Different kind of conditions will be considered, four mechanical conditions and three fluid conditions:

Mechanical boundary conditions:

Free edge: This mechanical conditions is given to any edges free to move. A general system of equations for this conditions is found in [START_REF] Landau | Theory of Elasticity[END_REF] (Equations 12.6 and 12.7):

       -∂∆u ∂n + (1 -ν) ∂ ∂l [cos θ sin θ(P -Q) + (sin 2 θ -cos 2 θ)u xy ] = 0 ∆u + (1 -ν) [2u xy sin θ cos θ -P sin 2 θ -Q cos 2 θ] = 0
Where n is the direction of the vector normal to the edge, l the direction of the tangent vector and θ the angle between the direction x and n. Theses equations are simplified when considering edges parallel to x or y. Thus for the edges parallel to x:

       ∂(P +Q) ∂y + (1 -ν) ∂uxy ∂x = 0 P + Q + (1 -ν)(-P ) = 0        P y + Q y + (1 -ν)P y = 0 Q + νP = 0        Q y + (2 -ν)P y = 0 Q + νP = 0 (II.13)
And for the edges parallel to y:

       P x + (2 -ν)Q x = 0 P + νQ = 0 (II.14)
Fixed edge: The fixed edge condition is to set a certain value for u while keeping its first space derivative along its normal vector null. For example u = z 0 and u x = 0 for an edge parallel to x.

Pivoting edge:

The pivoting edge conditions is also to set a certain value for u but this time it is the second space derivative along the normal vector which must be null. Which gives for example u = z 0 and P = 0 for an edge normal parallel to x.

Symmetry: A symmetry condition says that the same things happen on each side of the edge. This special condition is expressed as a null flux along the normal to the edge.

Fluid boundary conditions:

Reservoir: The fluid is free to go out at an edge with a reservoir condition and the volume of fluid outside the system is such that the external pressure of fluid is stable no matter how much fluid transit at this edge. Which leads to the expression v = 0.

Watertight:

The fluid cannot escape by this edge, which can be expressed as null flux.

Symmetry: A symmetry condition says that the same things happen on each side of the edge. This special condition is expressed as a null flux along the normal to the edge.

The edge 1 will always be our starting point for the bonding wave propagation. It will always stay in contact with the fixed bottom beam and so will have pivot mechanical boundary condition. A watertight condition could be used for edge 1 but a reservoir condition will be used most often as the small gap of z 0 imposed here will not let much fluid escape, and the condition is far better in term of convergence. Edges 2 and 4 are free edges with a reservoir condition and edge 3 is a symmetry edge.

The simple model consists of the bonding of a single moving rectangular plate on top of a unmoving plate. The displacement of the top plate is described with a 4 th order differential equation with 3 applied forces, the attraction forces, the gravity and the fluid pressure between the plates. The initial state is an equilibrium unmoving state with a simple 2 nd order profile along x for the position of the top wafer. Mechanical and fluid related boundaries conditions are used, such as clamp or free edges and reservoir or sealed edges. A numerical simulation based on this model will then be run.

II.2 Working with Comsol ®

To run our simulation the commercially available Comsol ® software was chosen. This choice was made because it is an all-in-one software simulating tool, including the drawing, the meshing, the solving and the post processing features as well as being oriented to solve multi-physics problems. A base license was acquired for our work and we worked with the math module, defining our problem with the aforementioned equations. First our model needed to be described in Comsol ® which will be succinctly explained, then convergence was obtained by choosing the right simulation parameters and we could start to use the results of the simulation. In the appendix A one can find a full example of creating this simple model in Comsol ® .

II.2.1 Translating the model

The first easy steps where to choose a two-dimensional space to work in, define our base functions (i.e. M att ) and parameters (See table II.2) and draw the geometry.

It is not possible to solve the two linked differential equations II.1 and II.3 as such in Comsol ® because Eq. II.1 is of 4 th order. We have to define two additional function and there associated differential equations to be able to write Eq. II.1 as a 2 nd order equation. Those two new functions are P and Q with their corresponding equations P = u xx and Q = u yy . Thus we have to implement in Comsol ® 4 differential equations of 2 nd order maximum, with their initial and boundary condition. This system of four equations is reminded below in Eq. II.15:

             D∆ 2 u = v+G+Matt D (1) P = u xx (2) Q = u yy (3) div(u 3 ∇v) = 12η ∂u ∂t (4) (II.15)
Translating the equations: Four equations have to be translated to the general form of differential equation of Comsol ® , which is for a general variable g:

e a ∂ 2 g ∂t 2 + d a ∂g ∂t + ∇ • Γ = f (II.16)
Noting that the scalar product with the operator ∇ is a divergence operator. Four inputs, e a , d a , Γ and f are needed for each equation. The indices 1,2,3 and 4 will respectively be associated with the variable u, P , Q, and v. Let it also be reminded that ∇ • ( ∇U ) = ∆U which applied to Eq. II.1 gives:

∇ • ( ∇(∆(u))) = v + G + M att D
To identify it to Eq. II.16, one defines the following input:

Γ 1 = ∇(∆(u)), f 1 = v + G + M att D
, e a = 0, and d a = 0 and to reduce the 4 th order of Γ 1 to a second order the function P and Q are used:

∆(u) = u xx + u yy = P + Q ⇒ ∇(∆(u)) =   P x + Q x P y + Q y  
And so equation (1) of the system II.15 is written in the general form of Eq. II.16 with:

Γ 1 =   P x + Q x P y + Q y   , f 1 = v + G + M att D
, e a = 0, and d a = 0

The same procedure is used to write the equations (2), (3) and (4) of the system II.15 in the general form of Eq. II.16 with:

Γ 2 =   u x 0   , f 2 =
P, e a = 0, and d a = 0

Γ 3 =   0 u y   , f 3 = Q, e a = 0, and d a = 0 Γ 4 =   u 3 v x u 3 v y   , f 4 = 12ηu
t , e a = 0, and d a = 0

Translating the initial state: For the initial state Comsol ® requires the value of each variables and there first time derivative. As stated in II.1.4 the derivatives will be taken null and the equation of the initial profile (Eq. II.8) quickly lead to the following initial state expressions:

at t = 0 : v = 0, u = t b 4R 2 x 2 + z 0 , P = t b 2R 2 , Q = 0 at t = 0 : v t = 0, u t = 0, P t = 0 and Q t = 0
Translating the boundary conditions: The different kind of boundary conditions that we can use in Comsol ® are:

• 0-flux conditions, which means that -n • Γ = 0, n being the vector normal to the boundary.

• Dirichlet conditions, which calls for a constant value for the variable g of the current equation.

• Constraint conditions, which calls for a constraint equation equal to 0.

Moreover the Individual dependant variable will be used for the Dirichlet and constraint conditions. Four boundary conditions are needed for each edge. Going back to our previous boundary conditions we translate each of these conditions in the form of the three type of conditions used by Comsol ® (see table II.1 for a summary). For the mechanical boundary conditions, one must give three boundary conditions in Comsol ® , one for each mechanical equation, i.e. the equation ( 1), ( 2) and (3) of the system II.15. The fluid boundary condition only requires a unique condition for equation (4) of the system II.15.

Free edge: Eq.II.13 relative to the mechanical free edge gives the constraint conditions Q y + (2 -ν)P y = 0 and Q + νP = 0 for an edge parallel to x. The third condition will be a 0-flux conditions, indeed we always have -n • Γ 2 = 0 on those same edges. See table II.1 for the edges parallel to y.

Fixed edge:

In the case of an edge of normal x, fixed at a height z 0 , we get the constraint conditions u = z 0 and u x = 0 and a 0-flux condition as -n • Γ 3 = 0 is always true on the edges of normal x.

Pivoting edge:

The pivoting edge condition of an edge of normal x, pivoting at a height z 0 , is expressed as two constraint conditions u = z 0 and u x = P and a 0-flux condition as -n • Γ 3 = 0 is always true on the edges of normal x. Symmetry: By definitions all symmetry conditions are written as 0-flux conditions.

Reservoir: This simple fluid condition is written as a Dirichlet condition v = 0.

Watertight: For the watertight fluid condition, the derivatives along the normal to the edge must be null, so a constraint condition v x = 0 is used for an edge of normal x.

II.2.2 Simulation parameters

With our model and all its equations ready we now need to mesh our geometry and to choose and configure a solver. As we are working with dynamics phenomenon, a transient study is chosen, with its time dependent solver. When solving the equation we need a convergence to have a solution and for that we have to monitor the mesh and the time step to obtain a realistic solution while keeping a reasonable computing time. The most common convergence problem arise from a great variation in a small space or in a short time, both relatively to the mesh size and the time step.

Condition Variable Edge of normal x [START_REF] Dowson | History of tribology[END_REF][START_REF] Lasky | Wafer bonding for silicon-on-insulator technologies[END_REF] Edge of normal y (2,3)

u P x + (2 -ν)Q x = 0 Q y + (2 -ν)P y = 0 Free edge P P = -νQ (P + νQ = 0) 0-flux Q 0-flux Q = -νP (Q + νP = 0) u u = z 0 u = z 0 Fixed edge P u x = 0 0-flux Q 0-flux u y = 0 u u = z 0 u = z 0 Pivot P P = 0 0-flux Q 0-flux Q = 0 u P Symmetry Q 0-flux 0-flux v Reservoir v v = 0 v = 0 Watertight v v x = 0 v y = 0
Table II.1: Summary of the boundary conditions. The pivot and fixed edge are considered at a height z 0 .

Mesh In a simulation the mesh is the space division or sampling of the geometry. If the mesh elements size is too big it is possible to miss small effects or not to converge due to the relative size of the gradient. But if it is too small it takes forever to compute. Indeed the smaller the mesh the more point you have to calculate. The total number of points is limited by the hardware used to compute the simulation. Two type of elements can be used to mesh a geometry, triangle and square. With our rectangle geometry the square shape is much better suited. After a few try it was observed that an element size of 0.5 cm works best for our simple model.

Time step

The time step we refer to is the time division or sampling of the total time of the experiment. This total time is usually of a few seconds. If we take a time step too big we risk convergence problems as the difference between the preceding step and the next one can be too important. But if we divide the time too much we have too much time step to calculate and the compilation takes forever. A time step of a tenth of a second was found to work fine with the base model.

II.2.3 Post processing the results

After all those steps, we finally have a running simulation which gives us results. But what are they and how do we use them? For each time step we have the value of each of our four variables. Most of the time the data will be extracted from the simulation and handled in IGOR Pro (see appendix D). The only things that will be used directly from Comsol ® are the 2D-surface graph and the animation. We will mostly be interested in the variable u with which the position of the bonding front can be determined. For that we simply look for the points at which z 0 -tol < u < z 0 + tol, with tol being a small tolerance. Then the position of the point farthest to the right gives the position of the bonding wave which will be plotted along time. The velocity of the bonding wave is given by the derivative of the position curve (Fig. II.5).

II.3 Variations to the simple model

The simple model presented above is just a base for the simulation. Below are presented modifications that can be brought to this model.

II.3.1 Wafer geometry

One of the modification that comes to mind is to go for the simulation of a full wafer.

The geometry stays simple but the boundary condition becomes quite complicated. Nevertheless it can be found in [START_REF] Landau | Theory of Elasticity[END_REF] and it was tried. But we couldn't get it to converge and we did not push this study farther as the use of such a geometry would have brought even more problems of scale differences and thus convergence.

II.3.2 Two moving wafers

Another possible modification is to consider two moving wafers. Rieutord's model used a symmetric system with two moving plates, and we showed that there is a movement of the bottom plate (see Subsec. IV.2.1). For that, a second set of parameters were added for the position of the second plate and its two 2 nd order space derivative. Some problems were raised when in some cases the bottom plate got to negative values of altitude. The second problem with this addition was also the extra time that it took to run it. So this modification to the simulation was not developed more, especially because the simulation results with only one moving wafer can be numerically adjusted, as it will be discussed in Chapter IV.

II.3.3 Initial curvature

In theory as our model consider a perfectly flat plate bonded at one edge, it will always bond to the end. But in reality each wafer has its own bow. To put this bow in our model we bond the perfectly flat top plate to a fixed curved bottom plate. Different profiles are used as a curvature, from a single bow to a sinus wave. Each time the equation of the profile is subtracted from the altitude in every existing equation. This addition works well and its only problem is that we actually fix the final shape of the assembly whereas in reality if one of the wafer as the profile described by the equation and the other is perfectly flat, the assembly final shape will stand somewhere in the middle. The following equations describe respectively the simple bow shape, (Eq. II.17) with b the bow at the center of the beam, and the waviness of the beam, (Eq. II.18) expressed as a sinusoid of amplitude a and with n period over the length of the beam. Both initial curvature profiles are plotted in Fig.

II.6. y = b - b R 2 (x -R) 2 (II.17) y = a sin(x n • π R -π) (II.18)
It should be remembered that to satisfy the small deflection hypothesis of the model bot a and b need to be kept small compared to the wafer thickness t w = 725 µm

II.3.4 Cavities

The addition of cavities is essentials to this thesis. To add cavities their shape was drawn inside the geometry and in those area the force M att was removed while a depth p was added to the variable u in the Eq. II.3. The main problem with those is that they have a relatively small scale and force us to be wary of the mesh. Generally cavities with a geometry similar to Fig. II.7 will be used. For this kind of cavity geometry a set of parameter has been set: the depth p of the cavities, their number nb c , their width l c , and the distance between the cavities w i . Other cavity geometry will be used and will be described in Chapter III. Whatever their design most of those have a symmetry axis at the center of the beam. A few additions to the simple model are considered. The wafer geometry did not worked, while the two moving wafers addition was not used a lot, mainly because of numerical problems. However the initial curvature and cavities addition will be used later.

Conclusion

A simple model was described along with some additions. Using this model and the Comsol ® software a simulation was built. This time-dependant simulation gives the position of the moving plate being bonded and from this position the velocity of the bonded wave is extracted.

Chapter III

Experimental setups Introduction

The experimental part of this thesis is aimed at validating the simulation and calibrating it. Those experiments were conducted in a clean room of the CEA-LETI in Grenoble. The base workflow of our experiments starts with a surface preparation of the silicon wafers to be bonded. Indeed, as we want to evaluate only the bonding wave, bare silicon wafers can be used without taking care of the bonding defects which could appear after thermal annealing [START_REF] Tong | SemiConductor Wafer Bonding: Science and Technology[END_REF]. Moreover, MEMC <001> 200 mm bright new silicon wafer completely fulfil the direct bonding physical specification (bow, waviness, roughness). To bond them correctly it is just needed to remove the residual particle contamination and to optimize the surface chemistry. After this slight surface cleaning, the silicon wafers are bonded while recording the advance of the bonding wave. Those bonded wafers then need to be cut into beams to suit the simulation geometry and for bonding and adhesion energy measurement accuracy. Once the beams are cut, various measurements are conducted. A description of the tools used for the experiments will follow, and then a more detailed workflow with variation will be presented. The various measurements made on the beams will then be described and a last section will focus on the implementation of the cavities.

III.1 Tools

Many of the following tools are used daily by the clean room technicians for operations similar to the ones we will use. Nevertheless some operations needed some tweaking or were completely new. Unless otherwise noted the author was trained to use the tools and was the one using them. The tools are sorted by uses in their description below.

Their name is the one used to identify them in the cleanroom and often contains brand or manufacturer name.

III.1.1 Surface cleaning and chemical treatment FSI Magellan

The FSI Magellan (Fig. III.1a) is an automated wet bench. This tool is mainly dedicated to direct bonding surface preparation. It allows us to have very clean surfaces with less than 10 particles measured by a surfscan between 90 nm and 500 nm and no particle bigger than 500 nm. It is then completely within the direct bonding specification. This tool also allows us to adapt the surface chemistry to the EVG850LT The EVG850LT is an automated bonding tool with some precursor steps like scrubbing and plasma treatments. All our bondings were done manually so we did not use it, except once for a plasma treatment.

Manual wet bench A few manual wet benches (Fig. III.1b) are available for a great range of study. In our case they were mostly use to etch some beams. When the etching was done by HF it was made by someone else as non-permanent LETI employees do not generally have the proper authorization. A set of two recipes will be used:

Si anisotropic etching: TMAH (Tetramethylammonium hydroxide) 12.5 % at 80 • C De-oxidising: HF 50 %, fast oxide etching speed, superior to 1 µm min -1

III.1.2 Oxidation

Tempress The Tempress is a high temperature furnace mostly used in this work for oxidation. A recipe containing the duration and temperature of each step is chosen and the wafers loaded before running the thermal treatment. In our case wet oxidation were performed at 1150 • C for 11 h when aiming for a 3 µm oxide layer. As a side effect of the growth at high temperature and because of the difference in coefficient of thermal dilatation between Si and SiO 2 the thin oxide films are stressed. Those stressed films can be used to induce a curvature to the wafer, by removing only one of them.

III.1.3 Thickness and depth measurements

Woolam The Woolam is an ellipsometer measurement tools. It is used to measure the thickness of thin film without contact. We use it to measure and check the oxide thickness obtained after thermal oxidation.

WYKO The Wyko is an optical interferometer measurement tools. It has two modes, one to measure steps height and another one focused on the roughness of the surface.

It was mainly used to measure the depth of the cavities.

AFM Atomic Force Microscopes (AFM) are used on a daily basis to image the surface of silicon wafer. The few AFM measurement performed on our sample were done by lab technicians. The images were used to measure the roughness of the bonded surface.

FRT The FRT is a tool using two confocal microscopes to measure the position of the upper and bottom surface of a sample. Its usual function is to measure the bow and the waviness of wafers or bonded wafers. But working with the tool owner, we have developed a procedure to record the movement of the wafers while they were bonded. However one of the limitations we faced was that the sensor had to move along a line, which is not such a problem as we know its speed. Fig. III.2 illustrate our measurement procedure, with the two sensors moving at a slow velocity against the bonding wave to be sure to encounter it.

Free fall Initiation of the bonding

Bonding wave propagation

III.1.4 Wafer cutting Disco4

The Disco4 is a dicing tool. It is used to cut wafers and pair of bonded wafers. It uses a rotating blade, 0.2 mm thick in our case. Most cuts were of 6 to 7 beams 2 cm wide per pair of bonded wafers. Both the distance between the cuts and their number (and also their direction if needed) can be changed at will. The bonded wafers are first taped to a frame and placed in the tools. Then the edges of the wafer (or any other reference point) are set-up and the recipe, once checked, is launched. When the cutting is over the beams are cleaned and removed from the tape to be stored for future use.

III.1.5 Bonding stations

Although automatic bonding station exist and were available (i.e. EVG850LT), manual bonding station were used. Both because unusual wafer shape were bonded (beams) and because it is easier to set-up a camera to record the bonding wave propagation.

Manual bonding station

This manual bonding station is a simple support with a small slope and guiding pins to keep the wafers aligned. A window is cut under the sample to allow the lights from under the support to shine through the sample and to use a camera to record an IR image. One of the very first task of this thesis was to set-up and configure a new camera. It is a grayscale camera recording 15 frames per second (fps) when scanning its full resolution of 2448 × 2058 pixels. To work with the beams, a special support was designed and the mechanical lab made it. The special support has a rectangular transparent plastic window to let the light pass and a clamping apparatus on one side. A plastic stylus is available to initiate the bonding, and a plastic blade to debond the samples when needed. The clamp is used to keep one side of the beam pair bonded when debonding most of the area by the insertion of the thick debonded blade. Indeed it was observed that if the beam pair is completely debonded lots of edge defect appear when rebonding the beams. Those defects make it impossible to correctly study the bonding wave. The cutting step is responsible and it is supposed that despite the water cooling of the rotating blade a local annealing occurs at the cut side. So when the beams are debonded the opening does not specially occurs at the bonding interface at the sides, but form a pattern. When rebonding, the slightest misalignment produce bonding edge defect as the pattern does not fit perfectly.

BAG

The BAG is a glove box (Boite À Gant in French) with an anhydrous nitrogen atmosphere. An airlock with vacuum and a desiccator are used to maintain the atmosphere anhydrous. A sensor measures the residual water concentration. Usually we work at a concentration below 0.1 ppm of water. The BAG wasn't designed to bond wafer but only to measure their bonding energy [START_REF] Fournel | Measurement of bonding energy in an anhydrous nitrogen atmosphere and its application to silicon direct bonding technology[END_REF] so the support has no slope but the light and camera are present.

MIL01

This MIL01 was acquired by the lab late in my thesis. Like the BAG it is designed to measure bonding energy but in a humid atmosphere. With some small tweaking it can be used to record bonding waves. It is equipped with the same camera than the simple bonding station. Its particularity is that it has a vaporized water and an nitrogen supply with a controller and sensors allowing the user to monitor the relative humidity (RH) inside its chamber.

Many tools are used for the experiments, and in most cases the author was trained to use them to be autonomous for most experiments.

III.2 Workflow

The tools described in the previous section are used to bring fresh out-of-the-box 200 mm silicon wafer to a state suited for direct bonding and to analyze this bonding. The simplest flow is made of 4 steps: First the surface preparation is done, followed by the bonding of the full wafers, then those bonded wafers are cut into beams and finally the beams are studied.

III.2.1 Step 1: Surface preparation

As seen previously, to be able to bond two substrates by direct bonding some surface requirement have to be met. The surfaces have to be clean and flat. Even thought it is possible to bond two bare Si wafers out of the box it is always better to have them cleaned, at least to remove the residual particle contamination and to setup the hydrophobic or hydrophilic character of the surfaces. For that a simple run into the FSI Magellan with either a hydrophilic or hydrophobic recipe is enough. Moreover in some cases it can be interesting to lower the initially quite good roughness of the bonding surface (usually around 1.3 nm RMS, measured on a 5 × 5 µm 2 AFM scan).

For that the easiest way is to grow an SiO 2 oxide of a few hundred nm which will give a bonding surface with a slightly lower roughness (i.e. a higher RMS value). Then to further decrease the roughness a few hundred nm of oxide are removed with a chemical etch, leaving an oxide surface with an increased RMS value. The oxidation is handled in the Tempress while the slow de-oxidation is done in the FSI Magellan.

III.2.2 Step 2: Full wafer bonding

After the surface preparation the wafers have to be bonded. The direct bonding is done on the manual bonding station. The two wafers are aligned by eyes, using the notch as a reference point and the initiation point is given by a light pressure with a stylus at the edge opposing the notch. The propagation of the bonding wave is recorded and the eventual bonding defects are localized. The pair of bonded wafers is then ready to be cut into beams.

III.2.3 Step 3: Cutting into beams

After bonding, the bonded wafers are cut with the Disco4. Usually beams 2 cm wide are cut, but some studies were conducted with beams from 1 cm to 6 cm wide.

III.2.4 Step 4: Beam study

When the beams are ready their study can begin. The bonded beams are first open and then re-bonded. The two main studies are the measure of the velocity of the bonding wave and the measure of the adhesion energy (E a ). Both are done using the bonding stations and the camera, the following sections (Subsec. III.3.1 and Subsec. III.3.2) will describe those measures in depth. On top of that some bonding energy (E c ) were measured, a few profiles of the deformation due to the bonding wave propagation were recorded with the FRT and in rare cases where we decreased the roughness, some beams were measured by AFM after all the other studies.

III.2.5 Additional steps: Patterning and Curvature

To expand this simple workflow, some additional steps can be added for specific purposes. The first is the addition of a pattern of cavities on one surface. In that case a wafer as to be patterned before the first cleaning step and when cutting the beams it is necessary to know precisely were to cut in order not to damage the pattern. The full patterning process with its specifics step is described in Sec. III.4. The second additional steps are needed when an initial curvature wants to be given to a wafer. For that an SiO 2 oxide layer is grown on the wafer prior to bonding, this growth results in two stressed oxide films. The wafer is then bonded and chemically de-oxidized, but thanks to the bonding only the backside external oxide is removed. As briefly explained in Subsec. III.1.2 this removal of only one of the stressed oxide layers results in the apparition of a bow. The de-oxidation step is done after cutting the wafers into beams, the lateral etching is small enough to be neglected.

The basic worflow is as follow:

A surface preparation step is done prior to the bonding of the full wafers, the pairs of bonded wafers are then cut into beams for subsequent study (Sec. III.3). On top of that basic workflow two additional step are considered, the etching of cavities to create patterned surfaces (Sec. III.4) and the curvature of a beam using a stressed oxide film on one side of a beam.

III.3 The main measurements

Three main measurements will be used experimentally. The first two, the measurement of the bonding wave velocity and of E a , will be used a lot and will be of great help to calibrate our simulation. The third, the measurement of the deformation profile won't be used as much, as the simulation model of only one moving wafer is quite far from the actual two moving wafers system measured here. The three measurements protocol will now be described.

III.3.1 The bonding wave velocity

The bonding wave velocity is the most used measurement in this thesis and will be extensively used and discussed in Chapter IV. Before the start of this work the velocity of the bonding wave was simply measured by starting a hand chronometer at the initiation and stopping it when the wafers were fully bonded, which was done by following the live stream from the camera of the bonding station (See Fig. 2). Although this method is fast and can give a rapid idea of the speed of the wave, it is not enough accurate for the purpose of our work. So at the start of this work a new camera was acquired and installed at the manual bonding station. When we developed the graphical user interface for the camera we took care of adjusting the recording feature to our purpose, see appendix B. Thus for each bonding we can record the advance of the bonding wave as a set of individual picture in a folder, with the time stamp of each picture written in a text file. With this output we then made a second piece of software with Matlab ® to exploit those image files. This home-made software uses the files from the camera and a calibration image for scale and gives the position of the bonding wave on each image along with the speed between images by reading the text file with the time stamps. It also generate a picture showing the position of the wave at each time, which will be called a chronographie, see Fig. III.5. The software user interface is presented in appendix C. The data file exported from the images treatment software is then post-processed in IGOR Pro, see appendix D. This post-processing calculate the average bonding velocity along three lines parallel to the bonding wave propagation direction. When bonding wafers, the recording is started after the two wafers are on top of one another, but before the initiation begins, whatever the initiation method used (gravity, stylus). When bonding beams, the recording sometimes includes the partial debonding, but usually begins at the initiation of the bonding by removal of the debonding blade.

Measure error: When experimentally measuring v, the value obtained has a measure error. Two kind of velocity measures must be differentiated. The first is the measure on full beams, by successive debonding and bonding cycles, detailed in Subsec. IV.2.3, and which experimentally gives an average velocity with an error of ±10 %. The second is the measure on beams with cavities on which the batch image treatment does not work because of the cavities. In this case the measure is usually done by considering only a small portion of the beam, and dividing manually the distance between the bonding fronts by the time between the images. Those measures are usually not plotted with their error which should be calculated individually because of the huge range of velocity considered.

III.3.2 Bonding and adhesion energies

As described in section I.2 there are two different energy measurements: the bonding (E c ) and adhesion (E a ) energy. In this work we will focus on E a measurements as it is the one relevant to the bonding dynamics and it will be extensively used and discussed in Chapter V. The setups used by Grierson [START_REF] Grierson | Characterization of Hysteresis of Surface Energy in Room-Temperature Direct Bonding Processes[END_REF] and Navarro [START_REF] Navarro | Direct Wafer Bonding Dynamics[END_REF] (presented in section I.2) are both quite complicated to setup. For our E a measurements a really simple method is used. A razor blade is put between the two wafers before they are bonded, so when the bonding occurs it is stopped by the inserted blade. The stopping distance is measured (See Fig. III.6) and the E a calculated from it. The advantage of using this setup is that it is similar to the DCB method used for E c measurement. The measure of the stopping distance is done by taking an image of the stopped wave with the same setup than when recording a bonding wave propagation and process the image with our home-made Matlab ® software (see appendix C). The distance L c measured is not the full unbonded length: as seen on Fig. III.6, two additions must be done. First the length b, which correspond to the blade edge length and is measured by optical microscopy, is added. Then a small distance a must be considered because the first interference fringe correspond to a gap of λ /4, λ is the wavelength of the infrared light going through the silicon (λ = 1 µm), but complete bonding is when there is no gap left. The distance a is deduced from an equation of the deformed profile from El-Zein [START_REF] El-Zein | Evaluation of GIC of a DCB specimen using an anisotropic solution[END_REF]. The full calculation of E a is done in an excel sheet based on Fournel [START_REF] Fournel | Measurement of bonding energy in an anhydrous nitrogen atmosphere and its application to silicon direct bonding technology[END_REF] work using El-Zein [START_REF] El-Zein | Evaluation of GIC of a DCB specimen using an anisotropic solution[END_REF] equations. Concerning the bonding energy measurement a few were made, usually after an adhesion energy measurement, by pushing a blade between to bonded wafers. The system obtained is the same as above and is treated the same way with a total unbonded length of L = L c + a + b. In both energy measurements, beams are more accurate than full wafers because of the equation hypothesis used to obtain the energy from the unbonded distance.

Measure error: When experimentally measuring E a , the value obtained has a measure error. To quantify this error let us consider the measure of E a for a pair of similar wafers. At the first order E a is written as:

E a = 3 16
E e t w 3 t b 2 L 4
This expression is only composed of product and quotient so the relative error ∆Ea Ea is easy to write as a sum of relative errors:

∆E a E a = ∆E e E e + 3 ∆t w t w + 2 ∆t b t b + 4 ∆L L
The relative error on the reduced young modulus E e is considered to be really low thanks to the use of the anisotropic El-Zein equation, and so gives a 1% error. The relative error from the thickness of the wafers is better than the specification of the manufacturer and reads:

∆t w t w = 5 725
The measure of the thickness of the blade used as a spacer is measured with a mechanical apparatus and gives:

∆t b t b = 5 310
Finally the error on the measure of L is the sum of the error on L c , a and b. The error on L c is due to the image treatment and can be as much as one interference fringes off. For a typical measure this correspond to an error ∆L c = 0.5 mm for a measured L c = 60 mm. The measure of the blade edge length, by optical microscopy, gives b = 1 mm and ∆b = 50 µm. The error on a is more difficult to obtain, but taking a worst case scenario it is possible to estimate it at ∆a a = 10 %, with a = 1.5 mm. So the full relative error on L is written:

∆L L = δL c + ∆a + ∆b L c + a + b = 0.

III.3.3 The deformation profile

The deformation profile measurements are done using the FRT (see Subsec. III.1.3). Those measurements were done both on plain wafers and beams. The restriction of having the sensors moving during the measure is not so important as the moving speed of the sensor was set really slow (0.2 mm s -1 ). Thus the sensors can be considered as almost immobile as the bonding wave velocity is at least 20 mm s -1 for the sample used. The measurement are made of 10 000 points spaced by 1 µm for a total length of 10 mm with a point measured every 5 µs.

The measurement protocols of the bonding wave velocity, the E a and the deformation profile were carefully setup to achieve quality measurements. The bonding wave velocity measurement will be extensively used and discussed in Chapter IV. It uses a home-made software to record images of the bonding wave propagation, then a Matlab ® script was developed to treat this images to locate the position of the bonding wave and finally the exported data are postprocessed using an IGOR Pro script. The E a measurement will be extensively used and discussed in Chapter V. It uses the same home-made image recording software, then our Matlab ® script to get the measured unbonded distance, and finally an excel sheet is used to calculate the final E a value.

III.4 Patterned wafers

To study the bonding of patterned surfaces we indeed need patterned surfaces. To obtain them it was first necessary to design the pattern then transfer it to silicon wafers. The transfer of the design on the silicon is done by photo-lithography. The patterned need to be drawn on a photo-resistive mask. After a very first batch with a simple array pattern to try out all the steps we were ready to design three different patterns for specific studies.

III.4.1 Photo-lithography principle

It was chosen to use a photo-lithography etching process with supple mask for practicality and cost management. Fig. III.7 shows the principle of this technique. First the silicon wafer is covered by a photo-sensible resin and a mask with the design to be transferred. The stack is then exposed to an ultraviolet light, which is blocked by the mask except where the pattern is. The UV light weakens the resin which is then easily removed. After removing the mask only the silicon and the patterned resin remains. The silicon can then be dry etched, the depth of the etching is time monitored. Finally the resin is stripped away from the sample, which is then ready for the classical steps, beginning with a surface cleaning.

III.4.2 Patterns design

The three patterns are visible on Fig. III.8. The tournament mask is designed to study the influence of the number and width of cavities on the bonding velocity, with a fixed bonding surface ratio. The triangle mask is designed to study the effect of a reduced bonding surface. The trench mask is designed to study the crossing capability of the bonding wave. The free and open source software vector graphics editor Inkscape is used to draw the patterns. The photo-resistive mask are to be printed on a supple mask, which is simply a transparent plastic film. The black ink will act as the photo-resistive mask. Each mask has seven 2 cm wide areas which will be the future beams, the center area is free of pattern for reference measurements (mainly E a ). The alignment of the mask and the Si wafer will be done with the notch. The seven beams are referred to as B 1 to B 7 from left to right, with the notch at the bottom. Here are the precise descriptions of each pattern parameters:

Tournament pattern: This patterned is covered with cavities that run parallel to the bonding wave propagation direction, with a fixed surface but an increasing number of cavities. On each beam the cavities number is doubled twice. The pattern of B 1 is as follow: five 1 mm wide and 40 mm long cavities which become ten 0.5 mm wide and 40 mm long cavities and then twenty 0.25 mm wide and 40 mm long cavities. This gives a beam surface with 25 % of cavities. For B 2 the starting point is ten 1 mm wide cavities, so 50 % of cavities. For B 3 the starting point is fourteen 1 mm wide cavities, so 70 % of cavities. On the right side the surface ratio are the same, but with fewer cavities, so two and a half 2 mm wide cavities for B 5 , five 2 mm wide cavities for B 6 and seven 2 mm wide cavities for B 7 .

Triangle pattern: This pattern is covered by triangles, with 18 mm bases and 140 mm height. The only difference between them is the position of the vertices. Three designs are identical, B 1 and B 7 , B 2 and B 6 and B 3 and B 5 . When bonding such a pattern from the vertex toward the base the bonding area decrease progressively and reach 10 % at the base of the triangles (18 mm of cavity for a beam of 20 mm).

Trench pattern: This pattern is an array of trenches that the bonding wave will have to cross. The trenches cover the full width of the beams, and are identical for the first three beams and the last three beams. The trenches are eavenly spaced: 4.99 mm from the middle of a trench to the middle of the next one on the left side and 4.98 mm on the right side. On the left side the first trench is 100 µm wide and each subsequent trench 20 µm wider, so the 26 th and last trench is 600 µm wide. On the right side the first trench is 200 µm wide and each subsequent trench 40 µm wider, so the 26 th and last trench is 1200 µm wide.

III.4.3 Patterned wafers workflow

Only one set of patterned wafers was made. For that a box of 25 Si wafers was run through the standard clean room wafer treatment. A specific workflow was designed to go from the box of bare wafers to the box of patterned wafers. Once the flow is translated into the LETI production software, the clean room technicians take care of the whole process. The 25 wafers will be referenced by their slot position in the box and named S 1 to S 25 . Each mask will be used on 8 wafers and the 25 th wafer will be a reference wafer. The simplified, i.e. without all the steps needed because of the internal working needs, workflow is as follow:

1. Oxidation: [S 1 to S 25 ] A thin SiO 2 oxide layer is grown on the wafers to protect the underlying clean Si surface.

Photo-lithography:

[S 1 to S 24 ] A resin layer is deposited on the wafers and a pattern from a mask is transferred to this resin, in three batches: The tournament mask is used for S 1 to S 8 , the triangle mask for S 9 to S 16 and the trench mask for S 17 to S 24 .

Etching: [S

Resin stripping:

[S 1 to S 24 ] The resin is stripped from the wafers.

Cavity cleaning: [S even

]At this point we got he box back, but before removing the protective oxide layer we tried to clean the cavities. Indeed little needle like defects (see Fig. III.9) are present in the cavities which is due to the use of a supple mask. To remove those, the sample were bathed into a TMAH solution for 15 min. Only the wafer with an even number were treated in this manner as previous tests showed that those defects did not prevent the bonding and because the TMAH treatment is adding depth to the cavities.

De-oxidation:

The protective oxide layer is removed in the FSI Magellan.

Depth measurement:

Finally the depth of each wafer cavities is measured with the WYKO.

After that the patterned wafers are ready to follow the previous workflow steps and be bonded to any other wafer.

Conclusion

All the experimental tools have now been described. Most of them are used by the author after a standard formation, but some were tweaked or completely made to suit our purpose. Along with the simulation we now have all that is needed to start the study of the bonding wave, which will start by the study of the bonding wave of plain wafers.

Chapter IV

The bonding wave velocity Introduction

The study of the bonding front velocity is important to better understand the mechanisms of the bonding front propagation. Firstly the coherence of our simulation model with the existing 1D law of velocity is checked, to fully justify our simulation choice. Then a study of the initiation and homogeneity of the bonding wave propagation is conducted. With the results of those preliminary studies, a comparison of experimental data and simulation results lead us to propose a 2D law of the direct bonding velocity of flat beams. This 2D law will finally be compared with the bonding of beams with cavities or bow. By default, the experimental conditions and simulation parameters used are those described as the standard ones in the previous chapters.

After checking the coherence of our simulation and presenting the initiation and homogeneity of the bonding wave propagation, a 2D law of the bonding velocity of flat beams is obtained and compared to the bonding of beam with cavities and bow.

Unless specified otherwise the bonding experiment are done with beams of width w = 2 cm cut in a pair of bonded wafer of diameter 2R = 20 cm and of individual thickness t w = 725 µm cleaned and left with an hydrophilic surface after their passage in the FSI Magellan. Likewise, unless specified, the simulations use the parameters from Tab. II.2 as a base, with a clamped edge, one symmetry edge and two free edges.

IV.1 Dynamic auto-coherence of the model

IV.1.1 Hypothesis and simulation choices

As it was explained in Chapter II a few choices and hypothesis were made when building our simulation model. The mechanical hypothesis of the small deformation (t b t w ) of a thin plate (t w w and t w 2R) are easy to justify by making sure that our simulation and experiment work with the required dimension. The choices we made are to use a simple adhesion force, to simulate only one moving plate and to find the bonding wave position by post-processing the simulation data. So before using our model in 2D, a comparison to an existing equivalent 1D analytical model is made.

IV.1.2 Comparison to the steady state 1D law

The steady state 1D law of the bonding wave velocity from Rieutord [START_REF] Rieutord | Dynamics of a Bonding Front[END_REF] (See chapter I) is used for the comparison. This law is obtained from the bonding of a symmetrical system of two moving infinitely large plates. The adhesion is a simple energy (the surface energy in the original paper, but the E a will be used here), but an experimental measure of the bonding profile gives a condition at the bonding front. The bonding velocity is supposed to be in a steady state. This all lead to the following 1D law of the bonding velocity (v 1D ), rewritten with our notations:

v 1D = A 1D E a 5 /4 z 0 1 /2 ηD 1 /4 ; A 1D = 0.95 3 /4 9 * 12 1 /4 (IV.1)
With A 1D a numerical factor, E a the adhesion energy, z 0 the cut-off distance (taken as the mean free path of the air particles), D the flexural rigidity of the plates and η the dynamic viscosity of the air.

To compare this law with our simulation a simulation is used, with a clamped edge at the initiation side, a free edge at the opposite side and symmetry conditions for the two lateral sides, to model an infinitely wide plate. With those conditions the results of our model should be the same that those from the 1D analytical model, with probably a pre-factor to account for some difference between the models. Simulations are run and as expected a steady state is obtained after a short initiation step (more details in the next section, Sec. IV.2) which was post-processed to get the bonding velocity, noted v ∞ . For each parameters of Eq. IV.1 (E a , z 0 , D and η) at least three values of v ∞ are calculated, by changing one parameter value while keeping the others constant at their base value, reminded on Tab. IV.1. The velocities v ∞ and v 1D are plotted along each parameter at the convenient power predicted by Eq. IV.1 to get a linear behaviour, for example E a Each graph shows that, as predicted, both the velocity v 1D and v ∞ are linear when plotted against the expected parameters. The slope of each linear function is indicated in the legend of each graph. By comparing each of those slopes the following expression is found:

v 1D = v ∞ √ 2 (IV.2)
So our simulation is coherent with the 1D law with a single factor of √ 2 as a difference. This factor can come from the difference between the two moving plates model and single moving plate model. It could be interesting to try to make our simulation work with two moving plate, but it was not achieved so far. However this result confirm our model choice of the adhesion force and post-processing of the bonding front position. The last thing to check is that the β parameter introduced by the adhesion force does not influence the velocity.

Our simulation model coherence is checked by the comparison with an equivalent 1D analytical law from the literature. The relation v

∞ = v 1D
√ 2 is found. Thus our simulations choices and hypothesis are justified, with an added limit of β ≤ 50 nm due to our adhesion force model.

IV.2 Initiation and homogeneity of the bonding wave

Before studying the velocity of the bonding wave, a few results regarding the initiation and homogeneity of the bonding wave are presented.

IV.2.1 Thickness of the air cushion

As previously explained different methods exist to initiate a direct bonding wave propagation. For two of them, initiation by stylus pressure and by gravity fall, a thin air cushion exists between the wafers. Using the FRT, described in Subsec. III.1.3, the evolution of the thickness of pairs of full wafers being bonded is recorded. Both initiation methods are used, but to correctly capture the bonding wave passage we need to know from where the bonding starts.

For the gravity fall a small bump is placed under the bottom wafer, to locally raise it and create a preferential initiation point for the bonding wave, see Fig. IV.3. As explained previously, the slow lateral translation of the FRT sensor is mandatory but is neglected in regards to the bonding velocity. The measure takes place on the side opposite the initiation point, at about 12 cm from it.

Bonding FRT sensor

Initiation area

Side view

Top view A t = 0 s, the measure is started, then around t = 20 s, the top wafers is released and starts its fall on the bottom wafer. The thickness of the assembly is first above the measurement limit of the FRT, then decreases rapidly when the top wafer is released. About 5 s after the release, the fall becomes really slow. Then the bonding wave initiates and propagates across the wafers. The rapid increase and decrease of the thickness around t = 116 s is the signature of the bonding wave macroscopic influence on the position of the wafers. The thickness stabilises just after t = 120 s, which means the wafers are bonded. The 4 s interval between the bonded state and the start of the deformation could very well correspond to the 12 cm between the initiation point and the sensor being travelled by the bonding wave at 30 mm s -1 . So the start of the deformation is the influence of the initiation of the bonding wave at the other side of the wafer. So the thickness of the air cushion before bonding is the 20 µm thickness difference between unbonded and bonded measured thickness. In case of a bonding initiated by a stylus pressure, the same setup is used with the initiation point on the side opposite the sensors side. The measure of the thickness, and the individual beam position are plotted on Fig. IV.5, all shifted to a common final value of 0 for easier observation. Two time intervals are defined on the graph, the longer, of 7 s, is the time between the first influence of the stylus on the system and the wafers being bonded at the sensor position. The shorter 6 s time interval is the time between the actual initiation of the bonding wave and the wafers being bonded at the sensor position. Observing the individual positions of the wafers at the beginning of the longer time interval, it is seen that at first both wafers go down while the thickness remain constant. This is because when first applying pressure with the stylus, the whole suspended support is pushed down. Then at the beginning of the shorter time interval, while the system keep on going down the thickness start to increase, which means that the initiation as started and the two wafers start to move away from each other at the sensor position because of a lever effect. The stylus is moved away from the system after the passage of the bonding wave under the sensor, when both wafers are going up while the thickness stays the same, at the bonded thickness. In this case the thickness of the air cushion, which is again the thickness difference before and after the bonding, is of 40 µm.

A simple simulation of the free fall of a Si 20 cm side square plate is run and gives a maximal air cushion thickness at its center as plotted on Fig. IV.4. The simulation gives a maximal air cushion thickness of 37 µm 10 s after the release of the top wafer. This simple simulation gives a correct order of magnitude for the air cushion which strengthens our results.

Those FRT measure also give some information on the vertical movement of the two wafers during the bonding propagation, which are going to be briefly studied below.

Vertical movement of the bottom wafer

The FRT measures used to determine the thickness of the air cushion can also be used to show that contrary to our simulation the two wafers do move during bonding. The position of the wafers of the measure with a stylus pressure initiation are already plotted on [START_REF] Tong | SemiConductor Wafer Bonding: Science and Technology[END_REF] measure is due to a derive of the sensor and can be safely ignored. In the same fashion the small bump occurring on both wafers after t = 121 s is a movement of the whole system, without any interest to us.

The observation of both measures shows that when the bonding starts at the other side of the wafer, the two wafers start to move away from each other at the sensor position. This movement can be explained by a double lever effect, with the adhesion force bringing the wafers closer on one side and a local air overpressure acting as a pivot, thus driving the other side of the wafer away from each other. Then the wafers start to get closer as the bonding wave, and the pivot point, travel toward the sensor, until they finally bond to one another. In the case of the stylus initiation the downward movement of the bottom wafer is limited by the support, while it is much more free to move when resting on a bump on one side as seen on The FRT measures are stopped here as we did not focus on the bonding profile shape, as it was not found to be important for the study of the bonding wave velocity. Moreover those studies are made on full wafers, while all our other studies consider beam bondings. However a further study would be of great interest to compare to and to improve the study of Navarro [START_REF] Navarro | Direct Wafer Bonding Dynamics[END_REF] on the bonding profile shape dynamic.

IV.2.2 Beam bonding: from initiation to steady state

Experimental study: When studying the bonding wave of beam bonding, the experimental initiation method used is always the removal of a blade, in order to prevent the formation of many edge defects which makes all measurements impossible, as seen in Subsec. III.1.5. The main parameter of such a bonding initiation is the thickness of the retracting blade. Fig. IV.7 shows the recorded velocity of the bonding of a beam, with an initiation with a thick blade (in red, about 5 mm thick) and with a thin blade (in blue, about 300 µm thick). The recording of the velocity spans only about 2 s for the thin blade because a thin blade induce a shorter unbonded length which is bonded in less time than when bonding almost the full beam length with the thick blade. For both bonding a stable velocity is reach after an initiation step. When giving the bonding velocity of a sample in the rest of the thesis, an average value of this stable velocity post initial step will be used, see appendix D for details.

The thick blade bonding as an initial velocity visibly higher than the stable average velocity, and last a little less than 2 s. With the thin blade the initial step is really short. Past the initiation step the stable velocity is of the same order whatever the blade thickness, with a velocity of 27 mm s -1 for the thick blade and 26.5 mm s -1 for the thin blade. So the variation of blade thickness in the initiation step does not influence the bonding wave velocity. The following experiments will use the thick blade for experimental reasons.

The higher initial bonding velocity when using a thick blade comes from an initial state with some extra energy. Indeed with the thick blade and the clamp at the other extremity of the beams, the top beam is deformed more than the equilibrium profile suggest. When retracting the blade the extra mechanical energy of deformation is released and contributes to the initial bonding wave velocity. Simulation study: As discussed in Subsec. II.1.4, the initial profile used for the simulation is a simple 2 nd order profile. To justify its use instead of the 3 rd order profile of the plate theory, simulations were run with different initial states. Blade thickness must be kept low to achieve convergence with the simulation, but the experiment showed that the steady state velocity did not depend on the blade thickness. Fig. IV.8 show the position of the bonding front along time of the bonding simulation. The bonding velocity is obtained by looking at the slope of the linear fit of the front position points. Blade thickness of 50, 100, 200 and 300 µm were simulated for the 2 nd order profiles and a blade of 300 µm was used for the 3 rd order profile.

The initiation step seems quite short for all initial bonding condition. Indeed the position front is quite linear really fast. To be sure to calculate the bonding velocity of the stable state the linear fit is done between t = 2 s and t = 4.5 s. The exact velocities obtained are displayed on Fig. IV.8, but they are all really close to 13.80 mm s -1 . The important difference of initial position between the 3 rd and 2 nd order profiles is due to the fact that our expression of the 2 nd order profile gives a full debonding whatever the blade thickness while the realistic 3 rd order profile has a debonding length linked to the blade thickness, see Eq. II.12. However this starting position does not influence the steady state velocity. So neither the profile order nor the blade thickness in the initial step influence the steady state bonding velocity, which is consistent with the experimental observations, and confirm that we can use the 2 nd order profile for our simulations.

IV.2.3 Homogeneity of repeated bonding cycles

The workflow described in Sec. III.2 present a first bonding, with full wafers, before the study of the bonding of the cut beams. This study of the bonding of beams is done by partially debonding the cut sample then let it rebond, and study this rebonding wave. As explained previously the samples are not completely debonded, thanks to the use of our special beam support with a clamp at one side, to prevent the apparition of many edge defects which makes the study of the bonding impossible. This partial debonding and rebonding is called a bonding cycle and can be cycled a few times. The more cycles are done the more the chance for a defect to appear is important.

To study the evolution of the bonding wave velocity for such bonding cycles, a pair of bonded beams is debonded and bonded 8 times, while also alternating the thickness of the debonding blade between the thick and thin blade of the previous discussion. For each cycle the velocity is measured at w /4, w /2 and 3w /4 and respectively called the left, middle and right velocity. Those velocities are average values of the steady state velocities measured by the batch treatment of images. Fig. IV.9 shows the result of these measurements. On top of the already discussed small variation due to the blade thickness, a velocity variation due to the position of the measure and another one due to the number of repeating cycles can be seen. Those variations are both small enough (under the 4% error bar of the average value of all the measurements) to let us use the average value as a unique value of the velocity of the bonding wave. This result justify the use of the bonding cycle's velocity measurements: when giving a value of the bonding velocity from an experiment, an average value of 3 to 5 successive bonding will be used. However the quite good dispersion of ±4 % showed here is the best we obtained, usually the error to consider is ±10 %. However as already mentioned not all experiments cannot be batch treated, especially the experiment with cavities. In those cases the velocity is measured simply by manually dividing the distance between the bonding front of two images, as far apart as possible from each other, by the time difference between those images.

The experimental measurements of the evolution of wafer pairs thickness being bonded evaluate the thickness of the air cushion present between the wafers prior to bonding between 20 to 40 µm, depending on the time after which the bonding initiate after the release of the top wafer. The vertical movements of both the top and bottom wafers are also observed.

An experimental study of the initiation of the bonding of beams showed that after a relatively short initial step, the bonding velocity reach a steady state. In the following discussion, except when specified, the value of this steady state velocity will be used. The same experimental study showed that the thickness of the blade used to initiate the bonding by its removal, does not influence the bonding wave steady state velocity. So in the next experiments the thick blade will be used by default to have an almost complete debonding and thus be able to measure the velocity on a longer distance. A simulation study of the same kind justifies the use of the simple 2 nd order profile given by Eq. II.8 and of a thin blade as initiation condition. The bonding velocity variations, measured by batch image treatment, both across a beam width and of successive bonding cycle, are small enough to let us use a unique average value for the velocity, with a 10% error bar.

IV.3 A 2D law of the direct bonding velocity of flat beams

After checking the coherence of our model in 1D and studied the initiation and homogeneity of the bonding wave propagation, it is time to study 2D bonding. In fact the comparison of the shape of a bonding wave of two flat wafers (Fig show the need of a 2D model. The second bonding wave is in advance in the area with cavities, which means that the bonding velocity is higher in this area. To begin the study of this 2D effect, the bonding of beams of varying width is studied.

IV.3.1 Velocity versus beam width variations

To study the bonding of beams of varying width, both simulation and experiments will be used.

Experimental study: Beams of width from 1 to 6 cm are cut from wafers as usual.

Then a value of the velocity of the bonding wave is obtained by averaging the value of six bonding cycle measurements. The resulting velocity is plotted along the beam width on Fig. IV.11. As expected the result show that the velocity increases when the width of the beams decreases. Moreover the adhesion energy of those beams was also measured (See Chapter V) and a unique value of E a = 25 mJ m -2 was found for all of them. So the velocity evolution seems to be only due to the width difference between the samples.

Simulation study: Similarly, simulations are run with w changing from 0.1 to 6 cm. Three different values of E a (25, 40 and 60 mJ m -2 ) are also used for each w. The other parameters are the base parameters, found in Tab. II.2. As already observed with the experiments the velocity increases when the width decreases, so the velocity is plotted against 1 /w, see Fig. IV.12. A good 2 nd order fit is found for the three different values of E a which leads to consider a 2 nd order law, to describe the velocity of the 2D model, of the following form:

v = v 0 + v 1 w + v 2 w 2 (IV.3)
Taking a beam of infinite width, Eq. IV.3 is reduced to the parameter v 0 which will be compared with v 1D and v ∞ from the 1D study, when the full expression of v 0 will be known in the next Section.

Comparison: A rapid comparison of the velocity evolution with beam width between experiment and simulation is plotted on Fig. IV.13. It shows a similar evolution which is promising but further study and calibration are needed which will now be developed.

IV.3.2 Velocity versus the other parameters

To continue the study of the velocity done with w, the influence of the other system parameters variation is now studied.

The other parameters:

In light of the previous study, the other parameters are the adhesion energy E a , the flexural rigidity D of the beam, the cut-off distance z 0 and the dynamic viscosity η. As with the 1D study, the parameter β does not influence the velocity of the bonding wave as shown by the results in Tab. IV.2. As previously with E a , for each other parameters three values are taken (the base value and one value above and under it) and simulations are run for every value of w used previously. So plots similar to Fig. IV.13 are obtained and 2 nd order fits gives a value for v 0 , v 1 and v 2 for each case. To obtain their exact expression a dimensional analysis is conducted.

Dimensional analysis:

The dimension of a variable will be noted between square brackets. For example the dimension of E a is noted [E a ] = M T 2 , expressed in the Mass, Length, Time fundamental system (MLT system). So all the parameters can be expressed in the MLT system:

[v] = LT -1 ; [E a ] = M T 2 ; [D] = M L 2 T -2 ; [η] = M L -1 T -1 ; [z 0 ] = L
Eq. IV.3 gives the dimension for v 0 , v 1 and v 2 :

[v 0 ] = LT -1 ; [v 1 ] = LT -2 ; [v 2 ] = LT -3
From what is already known from the 1D law we write:

[v 0 ] = [E a ] α 0 [z 0 ] δ 0 [η] γ 0 [D] ω 0 ; [v 1 ] = [E a ] α 1 [z 0 ] δ 1 [η] γ 1 [D] ω 1 ; [v 2 ] = [E a ] α 2 [z 0 ] δ 2 [η] γ 2 [D] ω 2
Which after a few steps gives the three following systems of 3 equations with 4 variables :

       α 0 = 1 + ω 0 δ 0 = 2ω 0 γ 0 = 1 ;        α 1 = 1 + ω 1 δ 1 = 1 + 2ω 1 γ 1 = 1 ;        α 2 = 1 + ω 2 δ 2 = 2 + 2ω 2 γ 2 = 1
More information are needed to fully get the expression of v 0 , v 1 and v 2 . For v 0 the 1D study gives the solution α 0 = 5 /4, δ 0 = 1 /2 and ω 0 = 1 /4. For v 1 and v 2 the value of the coefficient ω was hinted at by the result of the simulations as seen on Fig. IV.17 and gives ω 1 = 0 and ω 2 = -1 /4. So we finally get:

             α 0 = 5 /4 δ 0 = 1 /2 ω 0 = 1 /4 γ 0 = 1 ;              α 1 = 1 δ 1 = 0 ω 1 = 1 γ 1 = 1 ;              α 2 = 3 /4 δ 2 = 3 /2 ω 2 = -1 /4 γ 2 = 1
And the following expressions for v 0 , v 1 and v 2 :

v 0 = A 0 E a 5 /4 z 0 1 /2 ηD 1 /4 ; v 1 = A 1 E a z 0 η ; v 2 = A 2 E a 3 /4 z 0 3 /2 ηD -1 /4 (IV.4)
Full expressions of v 0 , v 1 and v 2 : To check the coherence of those expressions with the previous simulation results, the value of v 0 , v 1 and v 2 obtained from the fits are plotted against the expected parameters. For example v 0 , v 1 and v 2 obtained from the fits of Fig a , E a and E 3 /4 a respectively. It can be seen there that the linear dependence expected from Eq. IV.4 is correct. To finally get a full expressions of v 0 , v 1 and v 2 the only thing left is to find the value of the dimensionless factors A 0 , A 1 and A 2 . Those value are extracted from the slope of the linear fits used to check the coherence of Eq. IV.4 to get:

v 0 = 0.08 • E a 5 /4 z 0 1 /2 ηD 1 /4 ; v 1 = 0.4 • E a z 0 η ; v 2 = 2 • E a 3 /4 z 0 3 /2 ηD -1 /4 (IV.5)
Finally the full 2D law of the direct bonding velocity of flat beams is written:

IV.3. A 2D
v(w, E a , z 0 , η, D) = 0.08 • E a 5 /4 z 0 1 /2 ηD 1 /4 + 0.4 • E a z 0 ηw + 2 • E a 3 /4 z 0 3 /2 ηD -1 /4 w 2 (IV.6)
As mentioned earlier, if considering w = ∞ in the 2D law only the velocity v 0 remain. So v 0 is the 1D component of the 2D law. Moreover a comparison of v 0 with the other 1D velocity gives:

v 0 = v ∞ = √ 2 • v 1D .
As expected the 1D velocitiy from different sources are consistent with one another.

Reduction to a dimensionless equation:

To help understand this 2D law a reduction to a dimensionless equation is performed on Eq. IV.3 by dividing it by v 0 . Firstly, dividing v 1 and v 2 by v 0 gives:

v 1 v 0 = 5 • z 0 1 /2 • D 1 /4 E a 1 /4 ; v 2 v 0 = 25 • z 0 • D 1 /2 E a 1 /2 (IV.7)
Those two equations reveal that v 1 /v 0 = v 2 /v 0 and a distance d can then be defined as:

d = 5 • z 0 1 /2 • D 1 /4 E a 1 /4 (IV.8)
So we can reduce Eq. IV.3 to a a dimensionless equation:

v v 0 = 1 + d w + d w 2 (IV.9)
and rewrite the 2D law from Eq.IV.6:

v(w, E a , z 0 , η, D) = 0.08 • E a 5 /4 z 0 1 /2 ηD 1 /4   1 + d w + d w 2   (IV.10)
The reduction distance d: To help understand the physical meaning of the reduction distance d it is first compared to the radius of curvature of the dynamic profile presented by Rieutord [START_REF] Rieutord | Dynamics of a Bonding Front[END_REF]:

u dynamic (x) = l -2 /3 Ax 5 /3 ; l = D 3ηv 1D
But nothing there could be find that looked like d. So the radius of curvature of a static profile is then considered because with the fluid cut-off distance z 0 a static profile should occur first. Then when the distance between the two wafers is greater than z 0 the dynamic profile occur. So at a really small scale a static profile can be considered with a final thickness of z 0 as show on The static profile used comes from the thin plate theory [START_REF] Timoshenko | Theory of Plates and Shells[END_REF]:

u static (x) = z 0 2 2 - 3(L -x) L + L -x L 3
where L is the unbonded distance(from Eq. II.12):

L = 9 2 1 /4 z 0 1 /2 • D 1 /4 E a 1 /4
And thus: So d appears to be about three times the static unbonded distance L, more precisely d ≈ 3.4L. So when w is smaller than d the lateral air flow is greater than the air flow in front of the bonding wave, thus increasing the bonding velocity.

IV.3.3 Comparison of the 2D simulation law with experimental data

Now that we have the complete expression of the 2D law (Eq. IV.6) we can plot our experiment data against it. The parameters η and D are fairly well known and defined and z 0 is not an experimental parameter, so only the parameters w and E a are available for easy study. The experimental data from the study of the variation of the beam width were already presented and another set of data of the evolution of the velocity with the variation of E a is used. This last set of data will be more described in Chapter V. Remembering that the simulation bonds a single moving wafer to a fixed one while in the experiments the two wafers are free to move, another correction can be tried. Considering a symmetrical system the distance of the moving wafers to the final state should be divided by two which can be roughly translated by dividing the viscosity by 4 and the adhesion energy by 2. But once again, the law is a bit closer to the experimental point but the correlation to the trend worsen, as can be seen on Finally it was found that a correction giving good results is to add an offset value of 24 mJ m -2 to the E a as can be seen on Fig. IV.22. A calibration of the simulation with a systematic offset seems strange. Indeed the only difference known is due to the fact that the two wafer move during the experiment, as shown by the FRT measures in Subsec. IV.2.1, and this difference would be proportional to the value of E a .

However the systematic addition of 24 mJ m -2 to the simulation input adhesion energy E a to fit the experiments can also be seen as a systematic error of the adhesion energy experimental measure. The systematic 24 mJ m -2 missing in the experimental measure could very likely come from the interaction with the thin metal blade used to obtain the E a measure. Indeed the same blade is used to make all the E a measures, so this would explain a systematic offset. For example if the blade is slightly twisted, this twist would change the constraint along the y axis (which is the axis parallel to the short edges of the beam), which are not taken into account in the model used to measure E a . So this systematic offset which works well for all comparison for our experiments, should be calibrated if using another experimental setup.

Using the dimensionless equation (Eq. IV.9), the dimensionless law is plotted along with the experimental points on Fig. IV.23. As expected the dimensionless law fit well the experiment points. This dimensionless law will now be used to try to predict the bonding velocity of beam with cavities or bow. A 2D law of the direct bonding velocity of flat beams is obtained from simulation and dimensional analysis, see Eq. IV.6. It is then compared to experiment results and adjusted to best fit those by adding an offset value to E a . This systematic offset value could come from the metal blade used to measure E a , and thus not included in the value measured. A dimensionless equation (Eq. IV.9) is also written with a reduction distance d linked to the unbonded distance L of a static profile existing below z 0 .

IV.4 Application of the 2D law to beams with cavities or bow

The 2D law previously built is a 2D law of the bonding velocity of full flat beams. However the bonding of beam with cavities can be studied by comparing this 2D law with results of beams bonding with cavities and also with bow. As a general rule, bonding experiment with cavities are much more likely to have default. This is probably mainly due to the poor quality of the cavities we obtained from the cheap supple masks. So, many samples were not usable and the tendency to have default appear after a few cycles of debonding and rebonding limited the number of cycles used. Also the oxide layer grown on the wafer without cavity to potentially obtain a bow, hinder the batch image treatment capability of our Matlab ® script. Finally, the bonding wave velocity of beams with cavities can be really fast compared to what was previously observed which make it hard to know if or when a steady state is reached.

The following discussions will discuss the bonding of beams of width w = 2 cm, cut from three main cavity designs presented in Subsec. III.4.2. The samples are named by a slot number (S 1 to S 2 5) followed by a beam number (B 1 to B 7 ). After the study of triangle, tournament and trench cavity designs, the bonding velocity of bowed beam will be studied.

IV.4.1 Bonding of a triangle cavity

The experimental bonding of a triangle cavity is done from the side of the summit of the triangle, toward its base as shown on This variation of E a change the value of both v 0 and d accordingly, and for the law the previous offset is added to E a (w c ). A simulation of the bonding of the sample S 13 B 6 is also run (with the offset of 24 mJ m -2 added to E a ) and the resulting velocity is plotted on Fig. IV.25, with the experimental points and the corrected dimensionless law. It can be seen that the law predict quite well the experimental velocity with the first point (i.e. the smaller w c ) a bit higher than the law prediction. The simulation velocity is even closer to the experiment results, but also fail to predict the higher experimental velocity at the smaller w c .

So the 2D law is not too far off to predict the simulation velocity of a triangle cavity bonding, but a direct simulation is better still. Indeed the 2D law does not take into account the part of the beam that does not bond. But this part of the beam does exist and must be bent and does not correspond to the reservoir edge condition of the 2D law. However the simulation of the exact shape of the experiment, while much more time-consuming than a simple calculation using the 2D law, better represent the actual experimental measure. The value of w observed here are above 2d, but the study of the tournament cavities that will be done right after, nicely complete the study with value of w from about 0.1d to 1.6d.

IV.4. APPLICATION OF THE 2D LAW TO BEAMS WITH CAVITIES OR BOW65

IV.4.2 Bonding of tournament cavity

The experimental bonding of a tournament cavity is done from the side with less cavities toward the side with more cavities, so from zone 1 (Z 1 ) to zone 3 (Z 3 ) as illustrated on Fig. IV.26. It is important to remember that the tournament cavities where designed to have w c constant across the 3 zones. However the actual width to consider for comparison with the 2D law depends on the number of cavities, this width is noted w i and is the width of the individual bonding area between the cavities. For example for the zone 1 drawn on Fig. IV.26 where there is only one cavities we should take w i = w c /2. For the zones where the actual individual bonding widths are not equal, the higher width value is used. Although 6 different sets of tournament cavities were designed (see Subsec. III.4.2) only 3 will be de discussed here, B 1 , B 5 and B 6 . The exact designs are reminded in Tab. IV.3 along with the value of w i and the total bonding width ratio w c /w for B 1 , B 5 and B 6 . The other beam designs did not work well experimentally, their small value of w i made it sometimes impossible to correctly see the interferences fringes or their bonding velocity was too fast to be captured. The bonding velocity of the three zones of 5 samples (S 2 B 5 , S 3 B 1 , S 3 B 6 , S 8 B 1 , S 8 B 6 ,) will be presented and discussed.

Number of cavities Cavity width (mm)

w i (mm) Simulations of the 5 experimental samples were designed, using the exact experimental designs, and run. The simulation of each zone was done separately for a total number of 15 simulations to get the full results. Unlike our usual simulations, the geometry of B 5 does not have a symmetry line y = w /2, so for once the full width was simulated. The input simulation E a use the experimental measured values, around 35 mJ m -2 for all the samples, plus the offset of 24 mJ m -2 . For all the other parameters the base value of Tab. II.2 are used.

w c /w (%) Z1 Z2 Z3 Z1 Z2 Z3 Z1 Z2 Z3 w = 2 cm
Both the experimental and simulation bonding velocity of the 5 samples are plotted with the dimensionless law on The first general observations are that the experimental data does not follow the dimensionless law for every case, but that the simulation predict the experiments quite well. Indeed like for the triangle case, the simulation does take into account the full plate, with bonding and unbonding area alike, with the deformation and air flow modifi-cation expected. So for a prediction purpose a simulation, indeed more time consuming, is best suited. Once again the systematic offset of 24 mJ m -2 added to the simulation E a gives great results, further confirming our interpretation hypothesis.

As the measured E a of all the samples are really close (around 35 mJ m -2 ), the main difference between the two B 1 samples is their cavities depth. So the higher velocity of S 8 B 1 compared to S 3 B 1 is supposed to come from the deeper cavities. The increase in velocity with deeper cavities is also observable when comparing the samples S 3 B 6 and S 8 B 6 . Those comparisons of B 1 samples and B 6 samples also point out that when w i decreases the velocity difference increases. Indeed as w i decreases, it is normal that the velocity increases and so does the importance of the available cavity volume for a better air flow.

The experimental data point of sample S 8 B 6 on zone 3 is missing because we did not have two images showing a bonding front, the bonding wave was too fast for the camera sampling rate. The prediction of the simulation confirm this result as it predict a velocity of 462 mm s -1 . With such a speed, the wave would cross the 40 mm of zone 3 in a time t = 8.66 × 10 -2 s. The camera used should record 15 images per second in theory, but the actual count is closer to 12 images per second, so an image every 8.33 × 10 -2 s. So it is normal that the velocity could not be measured, and even with 15 images per second the chance of getting an image at the beginning of zone 3 and thus get a second image in zone 3 is really small.

The cavities depth influence on the bonding velocity observed with the previous results can also be specifically studied. However the simulations run for this purpose were run before the offset on E a was found, so the comparison between simulation and experiment is not really possible. The simulation were run at E a = 40 mJ m The results confirm that the velocity increases when the cavities depth, and so their volume, increases. Both simulations have the same total volume of cavities, but the distance between the cavities w i is smaller for the design B 1 . The higher velocity of the simulation B 1 Z 1 is due to the smaller w i .

IV.4.3 Bonding of trench cavity

The experimental bonding of a beam with trench cavities is done from the side with the shallowest trench toward the side with the widest trench. So the bonding wave must cross wider and wider trenches. The bonding front propagation stops after crossing a few trenches, even after a wait of about 270 s on S 17 B 1 , with cavities depth 0.75 µm. However if a stylus is used to bring the top beam in contact to the bottom beam, beyond the stopping trench, the bonding occurs without defect, so the stopping of the bonding is not due to the presence of defect ahead of the bonding front preventing its propagation, but really linked to the crossing of the trench. Sometimes, a second bonding wave starts at the other side after the first wave stopped (S 17 B 2 ) or before it stopped (S 22 B 3 ).

On the part where the bonding wave crosses the trenches, the bonding wave slow almost to a stop while crossing a trench, then a normal velocity is observed until the wave reaches the next trench. This gives a step evolution of the bonding wave position as seen on However, contrary to the expected behavior, no linear correlation between the time needed to cross a trench and its width is experimentally observed. Indeed the 180 µm wide trench is crossed faster than both the 160 and 200 µm wide ones. Both this seemingly incoherent crossing time of trenches and the stopping of the wave on some trenches is due to the experimental beams. Indeed to propagate the bonding after a trench, the top beam needs to come close enough to the edge of the trench for the adhesion force to bond the first point and propagate. But as the real beams are not perfectly flat, with a possible waviness, the rebonding point might take longer to come into the adhesion energy range, or not be able at all. This phenomenon is quite sensible as the adhesion force range is really short.

The average value of the velocity given between the trenches is lower than the reality because of the frequency of the image recorded. Indeed for the velocity between the 180 and 200 µm wide trenches only four images are recorded. On the 1 st image the wave is at the left trench and on the 4 th it is stopped at the right trench. Thus by dividing the distance between the trenches by the time difference from the 4 th image to the 1 st because the time difference is overestimated, the calculated average velocity is lower than the real one. Thus the difference between the 2D law velocity and the real velocity is even lower than the visible difference seen on Fig. IV.29. This means that the velocity between the trenches is well predicted by the 2D law.

The simulation of a bonding wave crossing trenches needs some modifications from the classical cavity simulation because of the really small width of the trenches. Indeed because of this small size of cavity the meshing size must be very small and lead to an very high number of nodes. So a simulation is designed with a beam of w = 2 cm as usual, but with a length 2R = 4 cm. Five trenches were drawn every 5 mm with width value of T 1 = 100 µm, T 2 = 160 µm, T 3 = 200 µm, T 4 = 300 µm and T 5 = 500 µm. An adhesion energy E a = 30 + 24 mJ m -2 while the other parameters use their usual value, found in Tab II. The simulation results show a similar trend than the experimental measure, which is that the bonding front slows almost to a stop when crossing the trenches. However because this simulation changed much from our usual simulation the value of the velocity is not relevant. A study of the mesh for this particular simulation would be required to get the correct velocity value. But the purpose here is just to show a general trend.

The simulations predict that the bonding wave should cross trenches wider than observed experimentally. This is not so surprising because in the simulation the beam equilibrium state is a flat state, so basically as the bonding started it will bond to the end. To stop a bonding wave in the simulation one needs to add some bow. The influence of the addition of bow on the beam will now be studied. But not on a simulation of trenches cavities because our model is not really suited to study it due to the simple adhesion force model used.

IV.4.4 Bow and velocity

Experiment

As already explained, experimentally a bow is added to some beams by thermally growing an oxide layer around the whole wafer then removing it only on the external backside of the beam. The resulting system is one with a thin compressed thin film of oxide on the interface side which induces a small bow to the beam. However with our thickest oxide, bow of only 120 µm maximum were measured for 200 mm long beam. With such a small bow no influence is really expected on the bonding wave velocity. So an experimental setup with a spacer placed under the bottom beam was designed to induce bow of greater amplitude as shown on Fig. IV.31. By changing the height H of the spacer the bow can be monitored, but it was observed that the plastic support under the bottom beam was also deformed by the spacer so the exact value of the bow of the bottom beam is not known accurately. Experimentally, a pair of beams is debonded and rebonded successively with different spacers of varying thickness H. An average velocity is calculated for the zone before, above and after the spacer and the results are written in Tab. IV.4, along with the 2D law value of the velocity with the measured E a = 56 mJ m -2 . For those measures both the plastic support and the bottom beam are clamped on both side, so it is possible to evaluate the bow of each. Taking a Young's modulus of the plastic of E p = 1 GPa with a thickness of t p = 3 mm it can be calculated that the bow b of the silicon beam is about a third of the spacer height, with E = 160 GPa the silicon beam Young's modulus and t w = 725 µm its thickness.

b = H E p t p 3 E p t p 3 + Et w 3 ≈ H 30.7 100
The velocity values of this table are fairly close and only a small decrease of the velocity at the higher bow could be noted. The approximated bow value of the bottom beam show that contrary to our expectation the bow of this setup is not much higher that the bow induced by a thin oxide layer. An improved setup which gives a bow b much closer to the spacer height H is used in Chapter V for the study of the adhesion energy measure of bowed beam.

The quasi nonexistent evolution of the velocity is unexpected and not understood as even those small bows should induce a relatively important energy cost. Indeed when bonding a flat beam to a bowed beam the flat beam must be deformed and this deformation costs some additional elastic energy. This energy cost should decrease the local adhesion energy and thus the local velocity. The energy cost of the bow obtained here can be calculated. Timoshenko [START_REF] Timoshenko | Theory of Plates and Shells[END_REF] writes the energy cost of the elastic deformation of a plate as:

E l = M 2 2D (IV.13) With M = -D ∂ 2 u
∂x 2 the moment of curvature and D the flexural rigidity of the plate. The profile u of the bowed beam depends on the way the beam is obtained. For a bow obtained by clamping a beam over a spacer one must consider the Eq. II.1 of the mechanical equilibrium of the plate, with no forces applied on the plate. Calling R b the distance between the clamp and the spacer and taking the clamp at x = 0 the following limit conditions are used:

at x = 0 : u(x) = 0 and u x (x) = 0 at x = R b : u(x) = b and u x (x) = 0
The resolution of the system gives a profile u(x) between x = 0 and x = R b of:

u(x) = b R b 2 3x 2 - 2x 3 R b x ∈ [0, R b ] (IV.14) So the moment of curvature M (x) = -D ∂ 2 u ∂x 2 is: M (x) = - 6bD R b 2 1 - 2x R b x ∈ [0, R b ]
and finally the elastic energy E l (x) is written:

E l (x) = 18b 2 D R b 4 1 - 2x R b 2 x ∈ [0, R b ] (IV.15)
The maximal value of E l is reached at x = 0 and x = R b . This maximal value is calculated for our case, with b = 200 µm, D = 5.54 J and R b = 9 cm and is E l(max) = 60.8 mJ m -2 . This maximal elastic energy is higher than the adhesion energy of this sample, E a = 56 mJ m -2 , so the wave should stop around the spacer position. Even if it does not stop it should really slow when crossing the spacer region. This incoherence of the experimental observation is not yet understood.

Simulation

As explained in Chapter II, in the simulations the bow of the bottom beam used is a simple 2 nd order profile (see Eq. II.17). This profile corresponds to a bow obtained with the oxide layer setup. Simulations of the bonding of a flat beam over a fixed bowed bottom beam are run. The value of the parameters used are the base value from Tab. II.2. With the simple 2 nd order profile, the curvature of the bow is constant and so is the elastic energy induced by the bow and a steady state velocity is observed. The steady state velocity of simulations with different bow are gathered in Tab. IV.5.

The elastic energy E l cost when bonding over a bow, is a part of E a that is not used to drive the bonding front. So the bonding velocity of a beam with a bow b with a plane beam with an adhesion energy E a should be similar to the bonding velocity of beams without bow but with an adhesion energy E a(eq) reduced by the elastic energy E l due to the bow. As previously E l is calculated with Eq. IV.13 but this time the profile u of the plate is given by the simple 2 nd order Eq. II.17, which gives us a constant value of M :

M = 2Db R 2
With b the bow value and R the half length of the plate and D its flexural rigidity. And so the elastic energy is:

E l = 2Db 2 R 4 (IV.16)
The elastic energy, calculated for the value of bow used in the simulations, and the 2D law prediction of the bonding two flat beam with E a(eq) = E a -E l are added to Tab. IV.5. The first observation is that the simulations show a decrease in velocity when the bow increases. Moreover as expected the bonding velocity of a pair of flat beams with an adhesion energy E a(eq) = E a -E l is similar to the bonding energy of the simulations with a bow.

Bow Velocity with bow Elastic energy stored E

l Velocity with E a -E l (µm) (mm s -1 ) (mJ m -2 ) (mm
It is important to remember that for all cavities bonding many bonding defects were observed during the experiments, limiting the available experimental data, and are imputed to the poor quality of our cavities. The 2D law obtained for the bonding front velocity of full flat beams is compared with the bonding velocity of beams with cavities or bow. The bonding velocity of triangle cavities is fairly well predicted by a simulation of the exact design. It is normal that the direct simulation predict better than the 2D law because the law does not take into account the unbonding area influence. For the bonding of tournament cavities this unbonding area influence is more important the 2D law does not fit well. However the 2D law can give an idea of the bonding velocity by a rapid numerical calculation, while a simulation takes a lot more time to design and run. The results also show the importance of the influence of the cavities depth p (and so their volume), especially at high velocity, the bonding velocity increasing when p increases. The bonding across trench cavities show a good fit with the law when bonding between the cavities. Moreover both the simulations and experiments show that the velocity decrease almost to a stop when crossing the trenches. Finally the bonding of flat beam over a bowed beam is studied. The experiments results are unexpected and no yet understood. However the simulation shows, as expected, a decrease in velocity when the bow increase. The decrease of the velocity is linked to the elastic energy cost of the bonding over a bowed beam.

Conclusion

This chapter studied the bonding wave velocity of beams. After checking the autocoherence of the simulation model described in Chapter II, the thickness of the air cushion present between the wafers prior bonding was determined. Then it was checked that the initiation step does not influence the steady state bonding wave velocity, for both experiments and simulations and that a single samples can be bonded and debonded a few times without significant impact.

Through the study of the influence of the beam width w, its flexural rigidity D, the adhesion energy E a , the fluid dynamic viscosity η and the cut-off distance z 0 on the bonding wave velocity v, a 2D law was obtained, see Eq. IV.6. This law was compared and calibrated with the experiment by adding an offset of 24 mJ m -2 . A dimensionless expression of the law was obtained by introducing a reduction distance d, which was found to be proportional to the unbonded distance L of a small scale static profile resting on the small cut-off length z 0 .

The 2D law obtained for the bonding front velocity of full flat beams is compared with the bonding velocity of beams with cavities or bow. The bonding velocity of triangle cavities is fairly well predicted by a simulation of the exact design. It is normal that the direct simulation predict better than the 2D law because the law does not take into account the unbonding area influence. For the bonding of tournament cavities this unbonding area influence is more important the 2D law does not fit well. However the 2D law can give an idea of the bonding velocity by a rapid numerical calculation, while a simulation takes a lot more time to design and run. The results also show the importance of the influence of the cavities depth p (and so their volume),

Introduction

The adhesion energy E a is really important when studying the direct bonding. In Chapter I it was seen that it is a fairly recent notion, indeed before only the bonding energy (E c ) was considered. But a hysteresis exist between the bonding energy E c measured when opening a bonding and the adhesion energy E a measured when closing a direct bonding. The measurement method used is described in Subsec. III.3.2 with a relative measure error of ±10 %. This chapter will first check the static auto-coherence of the simulation model, then confirm the homogeneity of the E a measurement on a full wafer and on beams with cavities. Then a study of the time dependence of the E a will be conducted and an explanation of the phenomena proposed. Follow a study of the influence of the relative atmosphere humidity (RH) with a first approach for a mechanism. Finally the link between bow and E a will be examined.

V.1 Static auto-coherence of the model

The simulation model presented in Chapter II and extensively used in Chapter IV to study the velocity should also be capable of predicting the static equilibrium profile of a pair of beam being bonded on a stopping inserted blade (see Fig. III.6). The simulations should give profiles similar to the analytical mechanical profile from Eq. II.9 given by Timoshenko [START_REF] Timoshenko | Theory of Plates and Shells[END_REF] for an infinitely wide plate. This analytical profile predict an unbonded distance L that is written in Eq. II.12.

Simulations are run, with different blade thickness to compare their profiles to the analytical one. The simulations use the base parameters of Tab. II.2, with the usual pivot condition on the left edge, a free edge condition for the lateral edge and a symmetry condition at the center of the beam. The edge with the blade uses a pivot condition at a total height z 0 +t b . A dynamic simulation is run with an initial state close to the predicted final state to limit the number of time step needed before reaching an unmoving equilibrium state. Three different blades thickness are used and the resulting profiles are plotted on simulation profile at the beam edge (so on the free edge). An additional simulation is run with two symmetry axis to get the profile of a plate of infinite width, for each case, but plotted only for t b = 300 µm. The zoomed static profiles show that a small difference exists between the analytical profile and the simulations. The simulation profiles of the free edges show that the edge of the plate goes a little under the equilibrium distance z 0 . This is mainly because the adhesion model repulsion is much weaker than the real solid interaction repulsion but it is required for convergence purpose. This influence of the beam width is also linked to the Poisson's ratio. However the simulation of infinitely wide wafer shows an almost perfect fit to the analytical profile, as shown on Fig.

V.3.
The simulated unbonded distance L for each case allow to calculate the E a corresponding. All the values are gathered in Tab. V.1. It shows that the small difference between finite and infinite beam width is not so critical, indeed by taking the value at the center of the beam the relative error are around 5 %. This result increases our trust in our simulation. The static auto-coherence of the model is checked. The profiles from the analytical 1D equation and the simulation of an infinite width beam fit almost perfectly. The simulation of 2 cm wide beam show a relative error between the input E a and the calculated E a at the center of the beam around 5 %.

t b (µm) w = ∞ w = 2 cm w = ∞ w = 2 cm E a(input) =

V.2 Homogeneity of the adhesion energy

The experimental measure of E a is not a real destructive measure method; however the inserted blade has a good probability of damaging the surfaces, thus creating bonding defect when attempting to bond it later. So it would be great not to have to measure the E a on every samples. Moreover in some cases, like the bonding of bowed beam, no direct model allowing to measure E a exist. So it would also be great to be able to measure E a on a reference sample and know the value of E a on a batch of samples. By having a look at the measure of E a across a single wafer and the measure of E a on beams with cavities we will show how it is possible to reduce the number of E a measures needed.

V.2.1 E a across a single wafer

Our samples come from cutting a pair of bonded wafer into beams. Because E a is defined as an energy by unit area, and the area of a wafer is quite homogeneous, it is expected that the E a should be similar across a wafer. Five beams of different widths are cut in a single pair of bonded wafers, and the energy adhesion measured at the center of each beams. The results are plotted on Fig. V.5. The variation of E a across the different sample is really small as expected. So it is possible to measure E a on a single sample and take the value as the value for all the samples cut from the same wafer pair. Usually it is also observed that wafers treated in the same batch have E a value close to each other, but sometimes a pair showed a very different value, so for security a sample of each wafer is measured.

V.2.2 E a with cavities

When describing the designs of the cavity patterns in Subsec. III.4.2 it was mentioned that a reference beam, numbered B 4 was kept clear of cavities on each pattern. Thus it is possible to measure the E a of the whole wafer with this sample B 4 , as shown just previously.

For the beam with cavities, the actual bonding area is just a portion of the total bonding area. So trying to measure an unbonded distance and calculate the associated adhesion energy on a beam with cavities is a bit different. For example, measuring an adhesion energy on a triangle cavity would give an unbonded distance L c linked to the corresponding width w c of the bonding front as drawn on Fig. V.6. The resulting apparent adhesion energy E ac should be linked to the full E a as:

E ac × w c = E a × w ⇔ E ac E a =

V.3 Time dependence of the adhesion energy

As reported by Turner [START_REF] Turner | Mechanics of direct wafer bonding[END_REF] a time dependence of the adhesion energy can be observed. With our setup the long time observation was easy to make by leaving the sample on the bonding station and taking images at regular time interval over a long time period.

Then the evolution of the unbonded distance L versus time can be measured and the corresponding E a calculated. The first general observation is that this time evolution of E a is observed only on hydrophilic samples, the couple of observations made with hydrophobic samples showed no evolution at all over a whole week-end. Moreover observations of the evolution of E a of hydrophilic samples in the dry atmosphere of the BAG (an anhydrous atmosphere glove box, see Subsec. III.1.5) showed a really small amplitude of variation compared to the experiment in standard clean room atmosphere.

Those first observations strongly suggest the importance of water on the phenomena. A standard hydrophilic sample of width w = 2 cm is bonded on an inserted blade of thickness t b = 310 µm and the subsequent unbonded distance L is recorded and plotted on The results show that E a increases with time until it reaches a maximal value, by describing a S curve. However the comparison of the maximal value of E a and minimal value of E c still show an hysteresis between the two variables. This increase in E a is really important as the value goes from about 40 mJ m -2 to a little more than 80 mJ m -2 . However this time evolution does not has an impact on our previous result as the energy that actually drive the bonding wave forward is the E a measured at the very moment on the stopping of the bonding wave.

Considering the importance of the water for this time evolution shown by the first observations, two hypotheses can be proposed for the phenomenon: a diffusion phenomenon or a thermally activated phenomenon. Both hypotheses are going to be tested.

It is important to note that for shorter data treatment time, the following results will only be presented with the value of the unbonded distance L.

V.3.1 Diffusion hypothesis

The diffusion hypothesis is the hypothesis which link the long time evolution of E a to a diffusion phenomenon toward the closet edge of the system which are the side of the beam. The time constant of a diffusion is linked to the diffusion coefficient and the diffusion length. To test the diffusion hypothesis two parameters are tested. First by changing the thickness of the inserted blade, the profile at the bonding front is changed, i.e. more open with a thicker blade, and it could influence a little the time to reach the final state as the water has a bigger height to diffuse. Secondly the beam width w is decreased, and in case of a diffusion phenomenon the time to reach the final state should be shorter as the diffusion distance to the outside reservoir is divided by two. For easier comparison each set of data is plotted with its own axis range, so that both the first and last values of each experimental L is plotted at the same level.

The first observation is that the change in blade thickness has no influence on the time evolution. This is not so surprising as the diffusion distance does not change when the blade thickness is changed. The second observation is that the experiment with w = 1 cm shows a slower time evolution. However this slower time evolution is not the expected behavior for a shorter diffusion length. As the temperature in these experiments is not monitored, we think that this unexpected result could be due to a different temperature in this specific experiment as compared to the other ones. The general conclusion is that we have not been able to highlight any effect of a diffusion phenomenon in this long time evolution of the adhesion energy.

V.3.2 Thermal activation hypothesis

To test the thermal activation hypothesis we needed to heat the sample. A glass pie plate put head down is used to create a small heating compartment, while being transparent to both visible and infrared lights. A thermocouple is put inside to read the atmosphere temperature under the plate. The heating is simply provided by the lamp used to see by infrared through the sample. At maximum power those lamps get quite hot. As the silicon thermalisation is relatively fast, it is assumed that the temperature read with the thermocouple device is the temperature of the sample. The air under the pie plate is heated by the lamp before introducing the sample. Then a usual long time measure of E a is conducted, with the lamps still heating the sample. So the temperature increases with time. The evolution of the measured L is plotted, with a few temperature measures on The first observation is that the sample reaches an equilibrium state after only around t = 3.5 ks, so about three times faster than for the non heated cases. Thus the hypothesis of a thermally activated phenomena, such as the thermally activation of capillary bridge, seems quite possible [START_REF] Noel | Sliding Velocity Dependence of Adhesion in a Nanometer-Sized Contact[END_REF]. The second observation is that at t = 328 s a sudden decrease of L is observed. This sudden change is not understood. As the importance of the water was indeed observed, the next study will be about the importance of the debonding atmosphere relative humidity RH on the measured value of E a .

The evolution of E a with time was observed and studied. This evolution is a long time evolution with a S curve when plotted against the log of time. It reaches an equilibrium state after about 10 ks in a standard clean room atmosphere. But the equilibrium state is reached faster when the experiment is done in a heated atmosphere which gives weight to the hypothesis of a phenomenon driven by the thermal activation of capillary bridges.

V.4 Adhesion energy and relative humidity

The previous study shows the importance of the water on the long time study of E a , indeed with a hydrophobic sample no evolution of E a is observed. The following study aims to understand how the relative humidity (RH) of the measuring atmosphere affects the adhesion energy E a . To measure E a at different RH, the MIL01, described in Subsec. IV.2.2, is used. However the setup does not monitor very accurately the value of RH in the MIL01 chamber, and the value of RH given should be considered with a error RH ± 5 %.

V.4.1 Observation of the phenomena

The first study was made by measuring E a at high humidity, over a long time to reach equilibrium, and then decrease RH. The results of the measure of the unbonded distance L is plotted on Fig. The very important observation is that decreasing the humidity RH after the equilibrium has been reached does not change the unbounded length. This tends to show that the equilibrium value of E a reached at RH = 80 % is also an equilibrium value for the lower value of RH = 50 % and RH = 10 %.

A second observation is the rapid evolution occurring 30 s after closing the cover of the chamber. We attribute this rapid evolution to the small increase in humidity (indeed measured by our sensor) due to the closing. The diffusion coefficient of water vapor in the atmosphere is D dif f = 24 × 10 -6 m s -1 . Considering the beam width w = 2 cm, this gives a diffusion time t dif f = w 2 /D dif f = 16 s. The time scale is comparable to the experimental feature.

A final observation is that after t = 40 s the long time evolution of the unbounded length is comparable to the one previously discussed, although slightly faster.

The second study is the reverse experiment: E a is measured over a long time at a low RH, then RH is increased to 50 % then to 80 %. The data is shown in Fig. V.12. The important result is that to the contrary of the previous experiment, increasing the humidity changes the value of the steady state E a reached for each value of RH. Therefore the value of E a obtained at RH = 10 % (25 mJ m -2 ) and at RH = 50 % (62 mJ m -2 ) cannot be considered as absolute equilibrium value for those humidity, as a value of 85 mJ m -2 could be obtained by decreasing the humidity from RH = 80 %.

We develop in the next paragraph a possible mechanism based on the capillary condensation of liquid bridges between rough surfaces for explaining this asymmetrical behavior of E a with respect to the change of RH.

V.4.2 Capillary condensation between rough surfaces

We attribute the time evolution of the adhesion energy E a in a humid atmosphere to the thermally activated capillary condensation of liquid water bridging the two surfaces in the contact line region. The phenomenon of capillary condensation can be characterized by the Kelvin radius, which is the radius of mean curvature of a liquid-vapor interface at equilibrium with an under-saturated vapor.

Kelvin radius

The Kelvin radius (r K ) is often used in the fields of porous materials and granular materials study. This Kelvin radius is used in the Kelvin relation and gives the critical radius under which a small cylindrical hole of a porous material will be filled by a liquid at equilibrium with the vapor. Similarly, for the study of granular materials, this Kelvin radius can be seen as the maximal curvature radius of the liquid bridge forming between the individual particles of the granular materials [START_REF] Bocquet | Moisture-induced ageing in granular media and the kinetics of capillary condensation[END_REF].

The expression of r K , adapted from [START_REF] Bocquet | Moisture-induced ageing in granular media and the kinetics of capillary condensation[END_REF], is used to calculate its value:

r K = - γ nk B T ln RH 100 (V.1)
where γ is the surface tension of the liquid-vapour interface, n is the number density of molecules in the liquid phase, k B the Boltzmann constant and T the temperature. With this expression it is possible to calculate r K in function of RH, which is plotted in Fig. V.13. The Kelvin radius does increase with increasing value of RH and is of the order of magnitude of the peak to valley roughness of our silicon surfaces. This result gives weight to the hypothesis that the long time evolution of E a in humid atmosphere is related to the kinetics of capillary condensation between the silicon wafers.

Proposed mechanism for measured adhesion energy in humid atmosphere

We propose a mechanism of a thermally activated capillary condensation between rough surfaces, on the model of Bocquet [START_REF] Bocquet | Moisture-induced ageing in granular media and the kinetics of capillary condensation[END_REF] for granular media and Noel [START_REF] Noel | Sliding Velocity Dependence of Adhesion in a Nanometer-Sized Contact[END_REF] for the adhesion of various surfaces including silicon wafers.

The general idea is that surface roughness creates energy barriers for the condensation of the wetting phase, which is here liquid water. These energy barriers induce a slow kinetics of the condensation, as the time needed to cross a barrier increases exponentially with its amplitude. We propose that the long time evolution of the adhesion energy is related to this slow kinetics of liquid condensation.

When the humidity is increased, E a increases rapidly at first, then more and more slowly as higher energy barriers have to be crossed. It eventually saturates in a metastable state where the remaining liquid phase, which should be stable at the current value of humidity, cannot condensate due to a too high barrier.

The value of this metastable E a reached for each value of the humidity, is a growing function of RH, because the height of the roughness-induced energy barriers is a function of the Kelvin radius and decreases when the latter grows. However when RH is decreased, the liquid bridges formed previously are still stable and thus the adhesion energy does not decrease, as can be seen on At a humidity RH = 80 % the value of the Kelvin radius r K is larger than the RMS value of the surface roughness. At equilibrium, the liquid phase is expected to fill completely the space between the two silicon surfaces in contact. The value of the energy barriers are globally lower, and induces a more rapid kinetics, as can be seen on Fig. V.12.

A schematic example of the formation of the capillary bridges is given on Fig. V.14. The E a measured is linked to the total area linked by capillary bridges. So the bigger the bridges are (due to the increased r K of higher RH), the higher E a is. Moreover the more bridges there is (due to the new condensation happening along time), the higher E a is.

RH = 10%

(a) At 10 % relative humidity (RH) the capillary bridges at the contact area condensate immediately (green lines). Then after some time, the energy barrier to condensate the other meta stable bridges is crossed and those bridges condensate (dotter green lines). The resulting E a is higher when the meta stable bridges are formed. Attempting to separate the two surfaces to measure E c gives a higher energy due to the presence of covalent bond at the contact points (red dots).

RH = 50%

(b) At higher RH = 50 %, due to the RH dependence of the Kelvin radius, the meta stable capillary bridges are bigger. The bridges at the contact points condensate immediately (blue lines). If the RH of 50 % is reached after reaching the equilibrium at 10 %, the existing bridges grow immediately to their equilibrium and the remaining meta stable bridges condensate after a time. If bonding directly at RH = 50 % all the meta stable bridges condensate one by one after a time.

Figure V.14: Example of the formation of capillary bridges between two rough surfaces in regards to humidity and time

The increase of the relative humidity RH of the atmosphere is found to increase the adhesion energy E a value. The proposed mechanism for this behavior involve the Kelvin radius r K . This radius, which gives the maximal radius of curvature of the capillary water bridges able to form between two surfaces, increases with RH. In the experimental condition the order of magnitude of r K is comparable to the bonded surface roughness, both have a nanometer scale. So when RH increase so does the number of capillary bridges possibly formed at the bonding front and thus the total E a measured.

V.5 Adhesion energy and bow limit

As already mentioned in Subsec. IV.4.4 the bonding of bowed wafers cost some extra energy. Indeed when bonding a flat beam to a bowed beam the flat beam must be deformed and this deformation costs some additional elastic energy. But if this elastic energy cost is higher than the available bonding energy, the bonding should not occur. So the bow of a beam is a possible limit for the bonding of a pair of beam.

An experiment is designed to measure this effect. We bond on a bowed beam, a flat beam with a triangle cavity (see Subsec. V.2.2) in order to have an effective adhesion energy which depends on the local position of the bonding front. We expect the bonding front to stop when the local effective adhesion energy (E ac ) is equal to the elastic energy due to the bow.

To have a bowed beam with enough amplitude the setup with a clamped beam resting on a spacer, as already shown on Fig. IV.31, is used. The top beam has a triangle cavity. The apparent adhesion energy E ac is proportional to the length w c of the local bonding front (w -w c is the local width of the cavity). To get the value of the elastic energy due to the bow, one should just need to measure the value of w c at which the bonding front stops and calculate E ac = E a w c /w with E a the nominal adhesion energy. The E a is measured on the reference beam B 4 of this sample, as explained in Subsec. V.2.1.

This experiment is done with the sample S 17 B 3 and for four different values of spacer height, H = 200, 350, 400 and 500 µm. Its adhesion energy E a is measured on sample S 17 B 4 and gives E a = 30 mJ m -2 . The measured values of w c and calculated values of E ac are gathered in Tab. V.2. For this experiment the spacer is put in a way that should reduce the deformation of the plastic support under the bottom wafer. The bow induced in the bottom beam by the spacer should thus be higher than the H /3 described in Subsec. IV.4.4, but still lower than the full spacer height H.

The first observation is that with the smaller spacer the bonding front does not stop. The triangle cavity is designed to give a minimal value of w c(min) = 2 mm at the base of the triangle cavity. The minimal adhesion energy associated is E ac(min) = 3 mJ m -2 . So for the spacer 200 µm as the wave does not stop, E l should be smaller than 3 mJ m -2 .

The second observation is linked to the relative position of the bonding wave compared to the spacer. For the spacer H = 500 µm the wave stops before the spacer, while for H = 350 µm and H = 400 µm it stops after the spacer. The stopping of the wave after the spacer is strange if the clamp conditions are strict. Indeed as previously calculated, the elastic energy E l of beam clamped on both side and resting on a spacer given by Eq. IV. [START_REF] Miki | Effect of nanoscale surface roughness on the bonding energy of direct-bonded silicon wafers[END_REF] The third observation is that even with the more favorable case of taking b = H /3 only the case with the spacer H = 200 µm is coherent with the experimental observation.

Despite our best effort, the study of the bonding of a flat beam on a bowed beam, does not yet achieve to match a model with the experimental observations. So our observations and measures of the stopping of the bonding wave are not able to give us a general prediction.

The bonding of a flat beam over a bowed beam is studied. A top flat wafer with a triangle cavity is used to have an effective adhesion energy which depends on the local position of the bonding front. The stopping of the bonding wave due to the elastic energy cost is measured. However no model could be find that is able to predict this behavior accurately.

Conclusion

This chapter studied the adhesion energy of bonding wafers. The static auto-coherence of our model was first check and a small influence of the beam width w on the unbonded distance L due to the Poisson's ratio effect is observed.

It was experimentally checked that the adhesion energy E a is homogeneous across a single wafer pair, within the 10 % error measure. Moreover when bonding beams with cavities the local apparent adhesion energy E ac can be obtained from the full E a , by using the corresponding bonding front width ratio w c /w.

Then the effect of water on E a were studied. First the long time evolution of E a observed was attributed to the thermally activated formation of capillary bridges. This phenomena of the capillary bridge formation was also used to explain the observation of the evolution of E a with the relative humidity RH of the atmosphere.

Finally the bonding of a flat beam over a bowed beam is studied. A top flat wafer with a triangle cavity is used to have an effective adhesion energy which depends on the local position of the bonding front. The stopping of the bonding wave due to the elastic energy cost is measured. However no model could be find that is able to predict this behavior accurately.

Conclusion

This thesis studied the direct bonding of patterned surface using a set of three patterns of cavities on silicon wafer beams. The aim was to be able to predict the bonding of a patterned surface to reduce the cost of testing each new pattern design experimentally.

For the prediction purpose a model of the bonding was built. This model is a 2D model of a top mobile beam bonding on a fixed bottom beam, with a fluid pressure due to the air between the two beams. Numerical simulations of this model were performed with the Comsol ® Multiphysic simulation software.

The simulation was first used to predict the bonding of two flat beams of width w but without any cavities. This studied resulted in a very good prediction tool in the form of a 2D law of the velocity depending on the beams width.

Then the velocity of beams with cavities was studied. The study of triangle and tournament cavities showed that the velocity could not be simply predicted by the previous 2D law of flat beams without cavities. However the direct simulation of the bonding with the cavities gave great predictions of the experimental velocity. An important result here was that by adding cavities, deep enough to have a good air flow, the bonding velocity of perfectly flat beams always increases.

The study of the bonding across trench cavities showed that between the trenches the velocity is well predicted by the 2D law. However the important influence of small variations of the waviness of the samples on the rebonding criterion after a trench is not predictable by our model.

The simulations predict that for perfectly flat beams, the bonding should always be possible. But perfectly flat wafers do not exist in reality, so the bonding of bowed wafer was studied and simulated. Observations of the stopping of the bonding wave were made but no predictive tool could be found. Further experimental studies should be made.

The main experimental limit of the bonding of beams with cavities observed were the high presence of local defect around the edges of the beams or cavities and the difficulty to cross trenches. In the first case, the local defects are probably due to the low quality of our cavities and beam edge, which would be easy to improve. For the crossing of trenches, the results tend to show that the presence of a small bonding bridge across the trench would greatly improve the capability of the wave to cross the trench.

While studying the bonding velocity the importance of the adhesion energy E a was highlighted and so a study of E a was also conducted. A method to measure E a was developed. The adhesion energy was found to be homogeneous across a pair of silicon wafers of 200 mm of diameter. Then the influence long time evolution of E a and its dependence with the relative humidity are both explained by a proposed mechanism of a thermally activated formation of capillary bridges between the rough surfaces.

As a final conclusion we propose a few guidelines to design patterns of cavities. The CHAPTER V. CONCLUSION bonding of perfectly flat cavities should always be possible, whatever the cavities used.

Limitations to the bonding of real wafers are due to the elastic energy cost of deforming the non-perfectly flat wafers. This limit is reached easily when a wave must cross a trench, so a design with a small bonding guide to help cross the cavity will work best. The width of this wave guide should be chosen by considering the bow of the wafer. Indeed the second important design rule is to keep a bonding area big enough to have more adhesion energy than the elastic energy cost due to non-flat wafers deformation.

To go beyond this thesis two main perspective can be suggested. The first, at a relative short term, would be to study further the bonding of bow wafer, with a better and more controlled experimental setup. Then to be able to give more precise pattern design rules, a model focusing on the mechanical deformation and the adhesion energy, but without the time dependent effect of the fluid, should be looked into.

A.2 Global definition

Once the model is created the first thing to do is to set a few Global definition, that will be used through the whole simulation. For the base model those definitions are a set of parameters, an analytic function and a piecewise function, each define in a node. To create the three nodes right click on the Global definition node and left click on the wanted sub nodes.

A.2.1 Parameters

The parameters needed are those defined in Tab II. Under the Definition drop down menu, the Expression field defines the actual mathematical expression of the function, taken from Eq. II.6. The Arguments field is the list of all the variable of the defined Expression, so only one variable, named u in the Expression above.

A.2.2 Analytic function

Under the Units drop down menu, the Arguments field is used to define the expected unit of the input variable. The square brackets are not necessary here as Comsol ® know that this field expects a unit. The Function field is used to define the unit of the output result of the function.

A.3 Geometry

After setting the Global definition, the Geometry must be defined. For the base model the geometry is really simple as it is just a rectangle. Because of the symmetry of the system explained in Chapter II the rectangle is of the dimension of half a beam. The Rectangle node is added by right clicking on the Geometry node.

A.4 Equations

The next step is to define all the equations for the four physics, i.e. the four General form PDE g, g1, g2 and g3. Each time the general equation must be entered along with its initial condition and boundary limits.

A.4.1 General form

As explained in Chapter II, the four equations of the system of Eq. II.15 must be put into the general form used by Comsol. In each General form PDE nodes under the four physics nodes the value of the vector Γ, the source term f , and the first and second order time derivative coefficient d a and e a must be entered. Those values are defined in Subsec. II.

A.5 Mesh

The mesh of this base model is defined by two nodes. The first node is the Size node, visible on Fig. A.10a, which is set to the Custom mode and defines the maximum and minimum element size, the maximum element growth rate, the resolutions of curvature and narrow regions. The second node is the Free Quad node, visible on Fig. A.10b, which mesh the Remaining domain with free quad mesh elements two times smaller than the size node definition in the direction x and y. This mesh was found to work properly with our base model but others (such as a Mapped mesh node) could also be used. To build the mesh two buttons are used, the first one is the blue building button with a red box, which build the mesh up to the selected node, and the build all button, with only the blue building on the icon.

Introduction

At the beginning of this thesis a new camera was acquired by the lab to observe and record the bonding waves. This appendix describe the user interface and possibility of the home-made image recording software, developed to work with the new camera. This appendix should prove especially useful for future user at the lab.

B.1 Initialisation

The software is started by using the short-cut CamCapture found on the desktop on the computer linked to the camera. If the camera was powered a few minutes before the start of the software, the main window that appears at when the software is run already display the camera images flux, as seen on Fig. B.1. If the image is black, wait a few second and check that the lamps below the sample are lit. If a message No GigE Vision Cameras Found appears it means that the camera and the software did not automatically linked. The easier solution is then to close the software, wait a few more second for the camera to start, and re launch the software. Another way would be to manually start the link in the camera settings, described in Sec. B.3.

B.2 The user interface

The main window of Fig. B.1 is divided in two main part. On the left side is the display area, were the image flux from the camera is displayed, with the live TimeStamp of the current image written in the top left corner. On the right side are all the buttons and control used by the user to monitor the images displayed and recorded. At the top of the control zone two groups called Settings and Zoom can be found. Under those group a third group, Stats, displays some live information about the current frame rates and number of image recorded. Then the larger and last group, the Recording group, is used to set the recording all the recording options. The groups briefly described here will be detailed below when the need arise.

value should be used all the time and so will not be detailed here.

The sliding cursor allow a direct change of the Exposure Time of the camera, i.e. the time during which the camera collect lights for a single image. Raising the value value of the Exposure Time give a brighter image and can be well suited for single image recording. However when recording a sequence of images, one should make sure that the Exposure Time stays below the time between two consecutive image capture. To adjust the brightness of the recording image, our setup also offers the possibility to directly change the lamp brightness.

The Camera button open a new window, with several tabs, which allow to set many option for the camera. The two tabs shown on

C.6 Results

The results tab (Fig. C.6) is used to plot and save the results of the batch analysis. The Calculer button use the time step in the text file to calculate the speed of each point between to images. The Corriger button correct some know errors that can happen when the camera lost some images and the time stamp are incorrect. Clicking on Afficher will plot the results in a 3D graph. The x axis (labelled Temps en s) is the time in second, the y axis (labelled Numéro du point) is the number of the point (so related to its position) and the z axis (labelled Vitesse en mm/s) the velocity of the point in mm s -1 , the v max field limit the maximum velocity value to be displayed (in mm s -1 ). Finally a click on Enregistrer will save the results, with the raw data in a text file on top of the simulation parameter and an Matlab ® .fig file of the plot. The raw data will be post-processed using IGOR Pro scripts, see appendix D.

Conclusion

The home-made Matlab ® software presented here was developed to analyze recordings of direct bonding wave. This presentation of the software shows how to use it and that the initial goal is met. However it is a specific piece of code and one has to be careful when using it as it is not aimed at a large audience. Of course it is best used with images from the CamCapture home-made recording software (Appendix B) and the results best processed with the specific IGOR Pro script (Appendix D).

Conclusion

Igor Pro is used to plot all the graphs of this thesis and to get the bonding velocity of both experimental and simulation data using two scripts written by the author and briefly described in this appendix.

Résumé

Le collage direct est un procédé par lequel deux surfaces suffisamment planes et propres peuvent se coller sans ajout d'un adhésif. Le collage direct de surfaces structurées est souvent utilisé pour la fabrication de système mécanique microélectronique (MEMS), où une plaque de silicium avec des cavités est collée à une autre plaque de silicium. La fabrication de ces dispositifs est chère et il serait utile d'avoir une ligne directrice lors du dessin de structures afin de savoir à l'avance si le collage direct sera possible.

Un modèle de simulation 2D pour le collage direct de deux substrats est développé et utilisé pour étudier l'influence des cavités sur la vitesse de propagation de l'onde de collage. Les prédications données par des simulations avec Comsol ® sont en bonne cohérence avec les mesures expérimentales et une loi en 2 dimensions de la vitesse de collage est obtenue. Le collage de plaques parfaitement planes avec des cavités serait toujours possible. Les limitations lors du collage de vraies plaques sont dues au coût de l'énergie élastique pour déformer les plaques non parfaitement planes. Cette limite est atteinte facilement quand l'onde de collage doit traverser une tranchée, dans ce cas un dessin avec un petit guide de collage pour aider à traverser la cavité fonctionnera mieux. La taille de ce guide d'onde doit être choisis en considèrent la flèche de la plaque. En effet la seconde règle importante du dessin est de garder une surface de collage suffisante pour avoir plus d'énergie d'adhésion que le coût en énergie élastique dû à la déformation des plaques non parfaitement planes.

L'énergie d'adhésion est un important paramètre du collage direct, car c'est l'énergie qui permet l'adhésion. Cette énergie d'adhésion est différente de l'énergie de collage la plus répandues qui est l'énergie requise pour séparer deux plaques précédemment collées. Dans cet ouvrage une méthode simple de mesure d'adhésion est proposée. Une mesure de l'évolution de l'énergie d'adhésion sur un temps long nous mène à proposer un mécanisme d'évolution lié à la formation de ponts capillaires entre des surfaces rugueuses.

  Don't be content with the what but get to know the why and the how. Robert Baden-Powell Rovering to success (1922) iii 1 Example of a MEMS device made using direct bonding: A piezoelectric micro pump . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Infrared (IR) images of a direct bonding wave propagation from initiation to completion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.1 Three parameters of the flatness of a wafer: bow, waviness and roughness, with range of classical values. . . . . . . . . . . . . . . . . . . . . I.2 Particle induced defects. . . . . . . . . . . . . . . . . . . . . . . . . . . I.3 Setup used to measure E c . . . . . . . . . . . . . . . . . . . . . . . . . . I.4 Sketch showing the velocity profile of the air escaping the closing gap between the two deformed wafers, from [17] . . . . . . . . . . . . . . . II.1 The base physical system to model. . . . . . . . . . . . . . . . . . . . . II.2 Geometry of the model. . . . . . . . . . . . . . . . . . . . . . . . . . . II.3 Plot of the potential energy well and its derived conservative attraction force. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II.4 The three possible initial states. . . . . . . . . . . . . . . . . . . . . . . II.5 Position of the bonding front versus bonding time. . . . . . . . . . . . . II.6 Profile of the bow and sinusoid initial curvature. . . . . . . . . . . . . . II.7 Geometry with four cavities of width l c and evenly spaced. . . . . . . . III.1 Surface cleaning and chemical treatment tools. . . . . . . . . . . . . . . III.2 Illustration of our special FRT measurement. . . . . . . . . . . . . . . . III.3 Manual bonding station and special beam support. . . . . . . . . . . . III.4 Special bonding stations. . . . . . . . . . . . . . . . . . . . . . . . . . . III.5 Chronographie of the direct bonding of two Si wafers. . . . . . . . . . . III.6 Adhesion energy measurement setup, with L c the unbonded length as seen by the camera, a the unbonded length before the first interference fringe and b the unbonded length of the blade edge. . . . . . . . . . . . III.7 Photo-lithography principle . . . . . . . . . . . . . . . . . . . . . . . . III.8 The three masks design . . . . . . . . . . . . . . . . . . . . . . . . . . . III.9 Optical microscope image of the defects after etching . . . . . . . . . . IV.1 Comparison of the evolution of the velocity from the 1D law and our simulation with an infinite plate. . . . . . . . . . . . . . . . . . . . . . IV.2 Velocity versus β, the parameter linked to the range of the attraction force model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV.3 Schematic views of the gravity fall on a bump FRT measure. . . . . . . IV.4 FRT measure of a gravity fall initiation. . . . . . . . . . . . . . . . . . IV.5 FRT measure of a stylus initiation. . . . . . . . . . . . . . . . . . . . . xi xii LIST OF FIGURES IV.6 FRT measurement of the position of the top and bottom wafers of a wafer pair bonded by gravity fall on a bump. . . . . . . . . . . . . . . . IV.7 Experimental initiation of beam bonding: effect of blade thickness. . . . IV.8 Simulation of the initiation of beam bonding: effect of blade thickness and initial profile law. . . . . . . . . . . . . . . . . . . . . . . . . . . . IV.9 Experimental measure of the velocity evolution of repeated bonding cycles. IV.10Comparison of the bonding wave shape for wafer bonding without and with cavities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV.11Experimental velocity versus beam width. . . . . . . . . . . . . . . . . IV.12Simulation velocity versus inverse of beam width. . . . . . . . . . . . . IV.13Simulation and experimental velocity versus beam width. . . . . . . . . IV.14Plot of computed value of v 0 , v 1 and v 2 versus E a 5 /4 , E a and E a 3 /4 respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV.15Plot of computed value of v 0 , v 1 and v 2 versus z 0 1 /2 , z 0 and z 0 3 /2 respectively. IV.16Plot of computed value of v 0 , v 1 and v 2 versus 1 /η. . . . . . . . . . . . . IV.17Plot of computed value of v 0 , v 1 and v 2 versus D -1 /4 , D and D 1 /4 respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV.18Physical meaning of the reduction dimension d. . . . . . . . . . . . . . IV.19Comparison of the 2D law with experimental data. . . . . . . . . . . . IV.20Comparison of the 2D law with doubled z 0 with experimental data. . . IV.21Comparison of the 2D law, with half adhesion energy and a fourth of the viscosity, with experimental data. . . . . . . . . . . . . . . . . . . . IV.22Comparison of the 2D law, with an offset adhesion energy, with experimental data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV.23Comparison of the dimensionless law with experimental data. . . . . . . IV.24Schematic view of the bonding of beams with triangle cavities. . . . . . IV.25Velocity of the bonding of a triangle cavity. . . . . . . . . . . . . . . . . IV.26Schematic view of the bonding of tournament cavities. . . . . . . . . . IV.27Velocity of the bonding of tournament cavities. . . . . . . . . . . . . . IV.28Velocity of the bonding of tournament cavities with cavities depth variation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV.29Experimental velocity of the bonding wave crossing trenches. . . . . . . IV.30Simulation of a bonding wave crossing trenches. . . . . . . . . . . . . . IV.31Bottom beam bow obtained by clamping it over a spacer. . . . . . . . . V.1 Static profile of a bonding beam stopped on a blade of thickness t b . . . V.2 Zoomed static profile of a bonding beam stopped on a blade of thickness t b = 100 µm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V.3 Zoomed static profile of a bonding beam stopped on a blade of thickness t b = 300 µm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V.4 Zoomed static profile of a bonding beam stopped on a blade of thickness t b = 500 µm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V.5 Measure of E a on beams of varying width, but cut from the same wafer pair. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V.6 Schematic view of adhesion energy measure on beam with cavities. . . . V.7 Adhesion energy measure on beams with triangle cavities (E ac ). . . . . V.8 Long time span E a and E c measure. . . . . . . . . . . . . . . . . . . . . V.9 Long time span E a measure, with w and t b variation. . . . . . . . . . .

Figure 1 :

 1 Figure 1: Example of a MEMS device made using direct bonding: A piezoelectric micro pump

Figure 2 :

 2 Figure 2: Infrared (IR) images of a direct bonding wave propagation from initiation to completion

( 1 -ν) 2 ReducedD = Ee•h 0 3 12

 123 Flexural rigidity of the beams R Radius of the wafers = half length of the beams w Width of the beams R rms Roughness RMS of the wafers E c Bonding energy (opening)

3 12 = E • t w 3 12( 1 -ν 2 )

 312312

Figure I. 1 :

 1 Figure I.1: Three parameters of the flatness of a wafer: bow, waviness and roughness, with range of classical values.

Figure I. 2 :

 2 Figure I.2: Particle induced defects.

Figure I. 3 :

 3 Figure I.3: Setup used to measure E c .

Figure I. 4 :

 4 Figure I.4: Sketch showing the velocity profile of the air escaping the closing gap between the two deformed wafers, from[START_REF] Rieutord | Dynamics of a Bonding Front[END_REF]

Figure II. 1 :

 1 Figure II.1: The base physical system to model.

Figure II. 2 :

 2 Figure II.2: Geometry of the model, with 2D space axis, boundary numbers (1-4), rectangle dimension (w, 2R) and bonding wave propagation direction.

  Fig. II.3 is a plot of both the energy well and the force with E a = 40 mJ, z 0 = 50 nm and β = 10 nm.

Figure II. 3 :

 3 Figure II.3: Plot of the potential energy well and its derived conservative attraction force.

Figure II. 4 :

 4 Figure II.4: The three possible initial states, with a flat profile along y for the top plate and a fixed bottom support.

Figure II. 5 : 1 .

 51 Figure II.5: Position of the bonding front versus bonding time. The slope b of the data gives the bonding velocity 13.8 mm s -1 .

Figure II. 6 :

 6 Figure II.6: Profile of the bow and sinusoid initial curvature.

Figure II. 7 :

 7 Figure II.7: Geometry with four cavities of width l c and evenly spaced.

  (a) FSI Magellan. (b) Manual chemical wet bench.

Figure III. 1 :

 1 Figure III.1: Surface cleaning and chemical treatment tools.

Figure III. 2 :

 2 Figure III.2: Illustration of our special FRT measurement.

  (a) Manual Bonding station. (b) Special beam support.

Figure III. 3 :

 3 Figure III.3: Manual bonding station and special beam support.

  (a) Anhydrous glove box (BAG). (b) Chamber with controlled RH (MIL01).

Figure III. 4 :

 4 Figure III.4: Special bonding stations.

Figure III. 5 :

 5 Figure III.5: Chronographie of the direct bonding of two Si wafers. The inconsistency of the gap between lines is due to the loss of images from the camera, not an actual local change in bonding velocity.

Figure III. 6 :

 6 Figure III.6: Adhesion energy measurement setup, with L c the unbonded length as seen by the camera, a the unbonded length before the first interference fringe and b the unbonded length of the blade edge.

Figure III. 7 :

 7 Figure III.7: Photo-lithography principle

  1 to S 24 ] The patterns are etched on the wafers, first through the thin protective oxide layer and then in the Si. Four different depths target are used: 1. [S 1 , S 2 , S 9 , S 10 , S 17 and S 18 ] Target of 1 µm deep cavities. 2. [S 3 , S 4 , S 11 , S 12 , S 19 and S 20 ] Target of 5 µm deep cavities.

  Figure III.8: The three masks design

Figure III. 9 :

 9 Figure III.9: Optical microscope image of the defects after etching

5 / 4 ,

 54 see Fig. IV.1. E a 40 mJ m -2 D 5.54 J η 1.85 × 10 -5 Pa s z 0 50 nmTable IV.1: Base parameters value used for 1D simulation.

2 .

 2 Once again the velocity is obtained by simulation for different values of β and plotted on Fig. IV.2. It can be seen that Velocity (mm/s) v ∞ Linear fit, slope = 0.63718 1D law, v 1D , slope = 0.45242 (a) Velocity versus E a 5 /4 . v ∞ Linear fit, slope = 51.426 1D law, v 1D , slope = 36.194 (b) Velocity versus z 0 1 /v ∞ Linear fit, slope = 0.017529 1D law, v 1D , slope = 0.012421 (c) Velocity versus D -1 /4 .

  Simulation points, v ∞ Linear fit, slope = 2.1198e-7 1D law, v 1D , slope = 1.4972e-7(d) Velocity versus 1 /η.

Figure IV. 1 :

 1 Figure IV.1: Comparison of the evolution of the velocity from the 1D law and our simulation with an infinite plate.

Figure IV. 2 :

 2 Figure IV.2: Velocity versus β, the parameter linked to the range of the attraction force model.

Figure IV. 3 :

 3 Figure IV.3: Schematic views of the gravity fall on a bump FRT measure.

  Figure IV.4: FRT measure of the thickness of a wafer pair bonded by gravity fall on a bump initiation.

Figure IV. 5 :

 5 Figure IV.5: FRT measure of the thickness of a wafer pair bonded by stylus pressure initiation.

  Fig. IV.5, and those from the measure with gravity fall initiation are plotted on Fig. IV.6. The small slope present on Fig. IV.

Figure IV. 6 :

 6 Figure IV.6: FRT measurement of the position of the top and bottom wafers of a wafer pair bonded by gravity fall on a bump.

  Figure IV.7: Experimental initiation of beam bonding: effect of blade thickness (thin ≈ 300 µm, thick ≈ 5 mm).

Figure IV. 8 :

 8 Figure IV.8: Simulation of the initiation of beam bonding: effect of blade thickness (from 50 to 300 µm) and initial profile law (2 nd or 3 rd order).

  Figure IV.9: Experimental measure of the velocity evolution of repeated bonding cycles.

  . IV.10a) and that of a wafer with cavities parallel to the bonding front propagation direction (Fig.IV.10b)

  (a) The bonding wave shape on wafer without cavities. The dotted blue line is the bonding front shape without cavities.(b) The bonding wave shape on wafer with cavities parrallel to the bonding wave propagation direction.

Figure IV. 10 :

 10 Figure IV.10: Comparison of the bonding wave shape for wafer bonding without and with cavities.

Figure IV. 11 :

 11 Figure IV.11: Experimental velocity versus beam width.

Figure IV. 13 :

 13 Figure IV.13: Simulation and experimental velocity versus beam width.

  . IV.12 are plotted on Fig. IV.14, against E 5 /4

Figure IV. 14 :

 14 Figure IV.14: Plot of computed value of v 0 , v 1 and v 2 versus E a 5 /4 , E a and E a 3 /4

Figure IV. 15 : 1 / 2 , z 0 and z 0 3 / 2

 151232 Figure IV.15: Plot of computed value of v 0 , v 1 and v 2 versus z 0 1 /2 , z 0 and z 0 3 /2 respectively.

Figure IV. 16 :

 16 Figure IV.16: Plot of computed value of v 0 , v 1 and v 2 versus 1 /η.

Figure IV. 17 :

 17 Figure IV.17: Plot of computed value of v 0 , v 1 and v 2 versus D -1 /4 , D and D 1 /4 respectively.

FigFigure IV. 18 :

 18 Figure IV.18: Physical meaning of the reduction dimension d.

IV. 3 .

 3 A 2D LAW OF THE DIRECT BONDING VELOCITY OF FLAT BEAMS 61

  Fig. IV.19 show the plot of those two sets of data with the 2D law plotted using the base parameters, reminded in the included table.

  Velocity versus adhesion energy.

E

  a E a(measured) = 25 mJ m -2 D 5.54 J w w (measured) = 2 cm z 0 50 nm η 1.85 × 10 -5 Pa s (c) Base parameters used.

Figure IV. 19 :

 19 Figure IV.19: Comparison of the 2D law with experimental data.

  Fig. IV.21.

  Velocity versus adhesion energy.

Figure IV. 20 :

 20 Figure IV.20: Comparison of the 2D law with doubled z 0 with experimental data.

  law with E a /2 and n/4 2D law base (a) Velocity versus width.

  law with E a /2 and n/4 2D law base (b) Velocity versus adhesion energy.

Figure IV. 21 :

 21 Figure IV.21: Comparison of the 2D law, with half adhesion energy and a fourth of the viscosity, with experimental data.

  Velocity versus adhesion energy.

Figure IV. 22 :

 22 Figure IV.22: Comparison of the 2D law, with an offset adhesion energy, with experimental data.

Figure IV. 23 :

 23 Figure IV.23: Comparison of the dimensionless law with experimental data.

Fig. IV. 24 .Figure IV. 24 :

 2424 Figure IV.24: Schematic view of the bonding of beams with triangle cavities.

  Figure IV.25: Velocity of the bonding of a triangle cavity.

Figure IV. 26 :

 26 Figure IV.26: Schematic view of the bonding of tournament cavities.

3 :

 3 Tournament cavities design parameters of beam B 1 , B 5 and B 6 .

  Fig. IV.27. As previously the value of v 0 and d use weighted value of E a (Eq. IV.12), but this time w c does not change with the position of the bonding front.

Figure IV. 27 :

 27 Figure IV.27: Velocity of the bonding of tournament cavities. The cross linked by the dotted line represent the experimental data, the different symbols represent the simulation data corresponding to each experimental parameters.

2 Figure IV. 28 :

 228 Figure IV.28: Velocity of the bonding of tournament cavities with cavities depth variation. With small cavities depths the velocity is slowed.

  Fig. IV.29.

Figure IV. 29 :

 29 Figure IV.29: Experimental velocity of the bonding wave crossing trenches. Data from sample S 17 B 1 with a measured E a = 30 mJ m -2 . The 2D law predict the velocity for a sample without trenches. The experimental value of the velocity is underestimated in comparison to the real velocity, see text.

2 .

 2 The position of the bonding front is plotted versus time on Fig. IV.30. The bonding front velocity, calculated by a simple ratio of the distance crossed by the bonding wave every 5 simulation time steps of 5 × 10 -3 s, is also plotted on Fig. IV.30.

5 Figure IV. 30 :

 530 Figure IV.30: Simulation of a bonding wave crossing trenches. The position of the bonding front (in red) is almost constant while the trenches are crossed. The trenches width are T 1 = 100 µm, T 2 = 160 µm, T 3 = 200 µm, T 4 = 300 µm and T 5 = 500 µm.

Figure IV. 31 :

 31 Figure IV.31: Bottom beam bow obtained by clamping it over a spacer.

Figure V. 1 :

 1 Figure V.1: Static profile of a bonding beam stopped on a blade of thickness t b . With base parameters value of E a = 40 mJ m -2 , D = 5.5476 J, z 0 = 50 nm and a simulation beam width w = 20 mm.

  Figure V.2: Zoomed static profile of a bonding beam stopped on a blade of thickness t b = 100 µm.

Figure V. 3 :Figure V. 4 :

 34 Figure V.3: Zoomed static profile of a bonding beam stopped on a blade of thickness t b = 300 µm.

Figure V. 5 :

 5 Figure V.5: Measure of E a on beams of varying width, but cut from the same wafer pair.

Figure V. 6 :

 6 Figure V.6: Schematic view of adhesion energy measure on beam with cavities.

Figure V. 8 :

 8 Figure V.8: Long time span E a and E c measure.

Figure V. 9 :

 9 Figure V.9: Long time span E a measure, with w and t b variation.

FigFigure V. 10 :

 10 Figure V.10: Long time span E a measure, with a sample in a heated atmosphere.

Figure V. 12 :

 12 Figure V.12: Long time span E a measure, with RH increasing.

Figure V. 13 :

 13 Figure V.13: Value of the Kelvin radius (r K ) versus RH. Calculated from Eq. V.1 with T = 22 • C, γ = 72 mJ m -2 and n = 33.36 × 10 27 m -3 . The typical value of the RMS roughness of our silicon surfaces is 0.2 nm, with a peak to valley of 2 nm.

  Fig. V.11. Moreover the energy hysteresis observed in Fig. V.8 when trying to reopen the bonding front, also prevent any decreased in the observed E a .

8 10 18 Figure V. 15 :

 81815 Figure V.15: Stopping adhesion energy E ac and elastic energy E l induced by a bow from clamp. The three highlighted point of E ac are the three experimental stopping points. The elastic energy are symmetric one the left and right side of the spacer positioned at x = 0.085 m, and only the left side is plotted. The more favorable case is considered, with b = H /3.

  (a) Add physics full window. (b) Units tab.

Figure A. 2 :

 2 Figure A.2: Add physics windows.

2 .

 2 For each parameter, enter the Name, Expression and Description in the corresponding field on the centre zone, after selecting the parameter node, as shown on Fig. A.4. In the Expression field, units must be entered between square brackets, and Comsol displays in the Value column the value of the parameters in the SI units.

Figure A. 4 :

 4 Figure A.4: Parameters window.

  An analytic function is used to define the adhesion, see Fig. A.5a. The Function name defined is used to call the function later on.

  Fig. A.6 shows the full window with the plotted rectangle in the right graphical zone. To plot the geometry click on the blue building button on the right top corner of the center zone.

Figure A. 6 :

 6 Figure A.6: Geometry drawing window.

  2.1 and Fig. A.7 shows the example for the equation g of variable u. A.4.2 Initial value Then the initial values are set, again those values are defined in Subsec. II.2.1 and Fig. A.8 shows the example for the equation g of variable u. appear in order to set the Individual dependent variables option one must click on the eye button on the right top of the Model Builder zone and check the Advanced Physics Options line, as shown on Fig. A.9a.

  (a) Mesh size definition window. (b) Mesh free quad elements definition window.

Figure A. 10 :

 10 Figure A.10: Mesh definitions windows.

Figure A. 14 :

 14 Figure A.14: 2D plot group window.

  (a) Surface window. (b) Height expression window.

Figure A. 15 :

 15 Figure A.15: Top beam display definition windows.

Figure A. 17 :

 17 Figure A.17: Profiles plot: 1D line graph window.

Figure A. 18 :

 18 Figure A.18: Bonding front position plot: 1D line graph window.

  Fig. B.2 shows a zoomed view of the three group with options.

  Fig. B.3 contain the most important camera option. The Source tab, see Fig. B.3a, is used to select the camera from which the image flux is displayed. If, as described in Sec. B.1, the camera is not started long enough before the software, it won't be recognised automatically. But selecting the camera in the drop down menu and checking the Acquire video checkbox should display the camera flux on the main windows, after clicking on OK. Apart from that the rest should be set as seen on Fig. B.3a. The Format tab, see Fig. B.3b, offer many options; however only the scan area options should be used. Changing the scan area allow to scan only a part of the camera full resolution, which allow to boost a little the maximum recording frame rate and save some disk space by not recording useless pixels. The first time one change the scan area, it important to first stop the acquisition by unchecking the Acquire video checkbox on the Source tab.With the camera used here only the Height of the scan area and the OffsetY can be changed. After changing the scan area the image displayed on the main window is the reduced image that will actually be recorded. However the zoom options might need to be changed to view a non-deformed image on the main window.

  (a) Source options (b) Format options

Figure B. 3 :

 3 Figure B.3: Details of Camera button options

Figure C. 4 :

 4 Figure C.4: Measurements tab, before batch analysis

Figure C. 5 :

 5 Figure C.5: Measurements tab, during batch analysis

Figure C. 6 :

 6 Figure C.6: Results tab

Figure D. 2 :

 2 Figure D.2: Screen-shot of the result of the script calculating the velocity of simulation bonding

  To calibrate this model and check its accuracy experiments are needed. The next chapter will describe those experiments.

	Name Expression	Value	Description
	Beam properties
	R	100 [mm]	0.1 m	Radius of the wafer / Half beam length
	w	20 [mm]	0.02 m	Width of the beam
	t w	725 [µm]	7.25E-4 m	Thickness of the beam
	E	160 [GPa]	1.6E11 Pa	Young's modulus
	ν	0.29	0.29	Poisson's ratio
	ρ	2.33 [g m -3 ]	2330 kg m -3 Density
	E e	E/(1 -ν 2 ) 1.7469E11 Pa Reduced Young's modulus
	D	E e × t 3 w /12	5.5476 J	Flexural rigidity
	Fluid
	η	1.85E-5[Pa s] 1.85E-5 Pa s Dynamic viscosity
	Mass
	g	9.81[m s -2 ]	9.81 m s -2	Gravitational constant
	G	-ρ × g × t w	16.572 Pa	Gravitational pressure
	Attraction force
	z 0	50 [nm]	5E-8 m	Equilibrium / cut-off distance
	E a	40 [mJ m -2 ]	0.04 J m -2	Adhesion energy
	β	10 [nm]	1E-8 m	Coefficient of the reach of adhesion
	Simulation
	t b	150 [µm]	1.5E-4 m	Thickness of the blade
	tol	15 [nm]	1.5E-8 m	Tolerance of the bonded distance
	Cavities
	p	10 [µm]	1E-5 m	Depth of the cavities
	l c	1 [mm]	1E-3 m	Width of the cavities
	nb c	5	5	Number of cavities
	w i	2.5 [mm]	2.5E-3 m	Distance between the cavities
	Initial curvature
	b 1	30 [µm]	3E-5 m	Bow of the plate
	a	5 [µm]	5E-6 m	Amplitude of the sinusoid
	n	2	2	Number of period
		Table II.2: Summary of the parameters.

Table IV .

 IV 5: Simulation velocity with different bow value. With E a = 40 mJ m -2 . E l is calculated with Eq. IV.16. The velocity with E a -E l should be similar to the velocity with bow, see text.

				s -1 )
	0	13.8	0	13.8
	200	13.3	4.44	12.03
	400	7.3	17.75	6.88

E a = 46 mJ/m 2 on S 13 B 4

  Figure V.7: Adhesion energy measure on beams with triangle cavities (E ac ). Plotted as a ratio of the E a measured on beams without cavities (B 4 ) versus the ratio of bonding front width w c /w.It was experimentally checked that the adhesion energy E a is homogeneous across a single wafer pair, within the 10 % error measure. Moreover when bonding beams with cavities the local apparent adhesion energy E ac can be obtained from the full E a , by using the corresponding bonding front width ratio w c /w.

	S 13 B 5
	S 13 B 6
	S 13 B 7
	E a = 52 mJ/m 2 on S 14 B 4
	S 14 B 5
	S 14 B 6
	S 14 B 7

  V.11, where RH started at 80 % and is decreased twice, first at 50 % then at 10 %, after waiting a long time before each decrease. The unbonded distance L is plotted instead of E a to save some time, indeed each E a must be calculated individually. The discussion can be made with either L or E a , remembering thatL 4 ∝ 1 /E a .Figure V.11: Measure of E a over a long time, with RH decreased twice. The raw measure of the unbonded distance L is plotted instead of E a , see text. The time axis is a log axis.

		80
	Unbonded length L (mm)	75 70 65 60	Cover closed at t = 10 s		RH 80 %	RH 50 %	RH 10 %
		55
			10 s	100 s	1 ks	10 ks
					Time (s)

Table V .

 V is maximal at the clamp and spacer position. So if the wave reach H (µm) w c (mm) E ac (mJ m -2 ) 2: Stopping adhesion energy E ac of a bowed triangle cavity bonding front. E ac is calculated from the experimental measure of the stopped bonding front length w c . E ac = E a w c /w with E a the nominal adhesion energy. the spacer it should be able to bond until the second clamp as it the adhesion energy E a c decrease more slowly than E l as can be seen on Fig. V.15.

	200	-	-
	350	6.66	10
	400	9.87	14.8
	500	12.02	18

TM: trade mark of SOITEC S.A.

Remerciements

especially at high velocity, the bonding velocity increasing when p increases.

The bonding across trench cavities shows a good fit with the law when bonding between the cavities. Moreover both the simulations and experiments show that the velocity decrease almost to a stop when crossing the trenches. Finally the bonding of flat beam over a bowed beam is studied. The experiments results are unexpected and no yet understood. However the simulation shows, as expected, a decrease in velocity when the bow increase. The decrease of the velocity is linked to the elastic energy cost of the bonding over a bowed beam.

Through the whole chapter the importance of the adhesion energy E a was observed and the next chapter will study it in depth.

Chapter V The adhesion energy

Appendix A

Comsol ® Introduction

The simulations of this thesis are run with the software Comsol ® Multiphysics. A base license is used as only general form equations from the base math module are used. The version used is the version 4.3b and the following appendix shows how to build the base simulation, presented in Chapter II. This appendix will be most useful for any following user.

The default main windows of Comsol ® is divided in three zone, the left zone contain the Model builder tree, the center zone display the option of the current selected node, and the right zone is the graphic zone where plots and geometry are drawn, see Fig. A.4. The Model builder tree shows the simulation structure and each line is called a node. Each node can have parent nodes or sub nodes. Left clicking on a node displays the option for this node on the center zone. Right clicking on a node offers various option such as renaming it, creating sub nodes or executing a node. This appendix will guide the reader from the creation of a new model file to the exportation of the simulation results.

A.1 New model

To start a new model, simply open the Comsol ® software or use the New option of the File menu of an already opened Comsol ® instance. The creation of a new model starts by the choice of the space dimension, the physics and the study type.

A.1.1 Select space dimension

When a new model is created the centre zone displays the space dimension window (Fig. A.1. For our model the 2D option is selected and clicking on the blue right arrow validates the choice and brings the next window.

A.1.2 Add physics

After selecting the 2D space dimension, the physics have to be selected. As explained on Chapter II, our model will solve the system of Eq. II.15. This system of four equations of the 2 nd order will use four mathematical General form PDE.

A.1.3 Select study type

The

A.2.3 Piecewise function

The piecewise function is created to be able to detect the bonding front position when post-processing the results. It is defined in the same way than an analytic function except that this time the Expression field is replaced by an intervals field where the user must define the function for each interval. Comsol ® offers extrapolation and smoothing options to help with the function definition process. The function defined on Fig . A.5b is a test function which returns 1 when the argument x is near the value of z 0 , with a small tolerance, and 0 otherwise. . This field requires a list of the times at which the system must be solved. Comsol ® provides an easy way to write this times list by using the range function, which has three argument. The first argument is the initial time, which is 0 in our case, the second argument is the time step, chosen to be 0.1 s and the last is the final time step, equal to 5 s for this case.

Once the time step and range are defined the user is strongly advised to save the current project. Indeed the computing time takes some time and if the computer crashes during the computation the whole model could be lost. Then the simulation can be launched using the Compute button, which is the green equal sign icon button. This base model computation time is around 2 h on the computer used. However it is possible to run around 3 or 4 simulations in parallel without significantly changing the computing time, so different parameters variation can be tested at the same time. Also it is important to remember that once the computation has started the only possible action is to switch between the different graphical windows on the right side area. The Model builder tree and the center option zone are frozen. It is however possible to stop the computation by clicking on the small red cross next to the progression bar in the right bottom corner of the Comsol ® window. During the computation the Progress tab next to the Messages tab display some information on the advancement of the computation. Two convergence graphs are also plotted by default and should be checked if the computing time seems too important or if the progress bar does not advance.

A.7 Results

Once the computation is over the user is again strongly advised to save the current work, which now contain the solution of the simulation, to prevent any unexpected loss. The exact computation time is written in the Messages tab under the graphical area. The data of the simulation can now be post-processed. A non-exhaustive example of the post-processing done in Comsol ® are presented below, followed in the next section by the exportation of the raw data for post-processing in another software.

A.7.1 Data set

The first possible step is to define a new data set using a Mirror 2D node to be able to plot the whole geometry of the system, using the half that was computed and building the other half by symmetry. As shown on

A.7.2 Derived value

The second step is to define a value derived from the data set to locate the position of the bonding wave at any given time. For that a Line maximum node (under the Derived Values parent node) can be used. Fig. A. [START_REF] Turner | Modeling of direct wafer bonding: Effect of wafer bow and etch patterns[END_REF] shows the option window of the node. The line 3 selected is the symmetry line so the choice of the data set is irrelevant here, but the Time selection must be set to All to have the position of the front at all the computed time. The expression of this value use the previously defined piecewise function nearZ0(x) that returns 1 when the top beam is close enough to z 0 , i.e. when it is bonded, and 0 otherwise. By multiplying this function by the position x and taking the maximum value, the returned value is the position of the last bonding point, which is the position of the bonded front. The yellow equal button is then used to evaluate the function and store the data in a table.

A.7.3 2D plot

A nice possibility with Comsol ® is to be able to plot nice 2D and 3D graph, for example the plot of the position of the top beam. A 2D plot node is used to define a new graph, on which different surface or contour plot can be added using Surface or Contour sub nodes. For example on the 2D Plot u graph, seen on

A.7.4 1D plot

A few 1D plot are also interesting to plot, and uses the same structure of a parent graph node followed by sub node each describing a plot on the graph.

The first 1D plot is the plot of the profile of the top beam position, taken on the symmetry axis, for all the time step considered (Fig. A.17). To plot the profile at every time step the Time is set to All on the data set selection of the parent plot 1D Plot profiles option.

The second 1D plot is the plot of the position of the bonding front along time (Fig. A.18). It is a 1D plot that uses the value from the table created when evaluating the Line maximum derived value, so it displays the position of the bonding front along time on the symmetry line. The slope of this graph would give the velocity of the bonding wave, but we don't know if there is a way to get its value inside Comsol ® , so the data will be exported to find this velocity value.

A.8 Export

Conclusion

This appendix is a guide for the user who wants to reproduce the base simulation described in this thesis. The user is then welcome to implement any of the additions to the model proposed in Chapter II or of its own making.

Appendix B

Home-made image recording software . By default the checkbox Zoom fit is checked, and so the camera flux is displayed to fit in the window, which by default has the same image size ratio than the camera. But if changing the size of the main window, or more often the scan area of the camera, the image size ratio are not the same and the displayed image is deformed. Thus one need to uncheck the Zoom fit checkbox and use the now enabled Zoom + and Zoom -button the adjust the display.

All the zoom options only affect the displayed image in the main window; the recorded images will always be of the same size as the camera scan area.

B.5 Recording options Conclusion

This simple interface was developed to have a few interesting option for the recording of single or multiple images from the camera. The live display is not stopped during the recording so the user can check the general position of the bonding wave being recorded. The data will then be batch treated (if a sufficient contrast is present) with the Matlab based software, presented next in appendix. C

Appendix C

Velocity and adhesion energy measurements with Matlab ® Introduction

To analyze all the images obtained experimentally a home-made software was made with Comsol ® code. This appendix will quickly describe its working and its possibilities. It will also serve as a guide for future user at the lab. For those users, two executable are available, one for Windows 32 bit, the other for Windows 64 bit: Mesure vitesse32.exe and Mesure vitesse64.exe. They were compiled using the version R2013a of Matlab ® and require the version 7.17 of the MCR (Matlab Compiler Runtime) to be installed if Matlab ® is not installed. The current version of this software is version 2.2 and is indicated in the bottom left corner of the window.

C.1 Working principle and pre-requisite

This program was designed to be used with images coming from the manual bonding station camera, recorded using the CamCapture software (See appendix B). Those images are grayscale images and need to be numbered for batch analysis, with a base name and a six digit number, so starting at name000001 for example. A text file named name.txt containing the time stamp of each image is also required in the folder containing all the images. When analysing an image it is the color gradient that is looked at. The image then gives a collection of numbered objects corresponding to the higher gradient area.

C.2 Calibration

The loading screen allows the user to reload a previously saved study. To start a new study go to the calibration tab (Fig. C.1) by clicking on the Suivant button. The calibration is intended to work with an image of a disk of known diameter, if you do not have such an image it is possible to manually enter a calibration on the Paramètres tab (Fig. C.3). To calibrate the tool, first load the calibration image, using the button Parcourir to get the path to the image and then click on Afficher to load and display it. Then select a rectangle zone including the calibrating disk by clicking on Dessiner un rectangle then draw it by clicking and dragging on the image. The rectangle can then manually be modified using the four boxes and redrawn using the Actualiser le tracer button. The frame Traitement de l'image allows the user to adjust the border detection parameters. Sigma is usually taken between 1 to 12 and the higher it is the more the detected border must have a sharp edge. Seuil de détection can be set to auto or usually taken between 0.01 to 0.4 and the higher it is the more the detected border must be contrasted. Finally the Petit objet set the size limit (in pixel) under which objects are ignored. Both Actualiser button refresh the display with the current parameters values. For this calibration step one wants to detect only the edge of the disc, or a part of it, indeed the Calibrer button will try to fit all the objects to an ellipse and calibrate both axis to the real size (in mm) indicated in the box Diamètre réelle. The axis of the ellipse will be displayed in blue and the ellipse itself in red on top of the original image for quick checking of the calibration. If the calibration seems good proceed to next step by clicking Suivant, else start over and try to adjust the parameters to detect only the edge of the calibrating disc.

C.3 Selections

The selection tab (Fig. C.2) is used to select the images to be analysed and how. The first step on this tab is to select an image file using the Parcourir button. The other field of the top frame will be automatically filled with the information from this first image. The base name is under Racine du nom des images, the first image number under Numéro de la première image, the last image number under Numéro de la dernière image, the format of the image file under Format and the full name of the text file under Nom du fichier texte. Those field can be manually adjusted, especially the starting and last image number which probably won't be correct as it is. Then click Afficher to display the image. A selection rectangle needs to be defined following the same steps than in the calibration section (Sec. C.2). Here the aim is to detect the edge of the bonding wave which will move, hence the starting selection must include both the starting and last position of the bonding wave (usually the full wafer or beam). Finally the last frame (Sens de l'onde) indicates the propagation direction of the bonding wave. It is a great data analysis and graphing tool, with lots of programming option. It was first aimed at working with wave type data, i.e. 1D evenly spaced data, so its main data format is a wave object. An important feature is also its really good graphing options and the possibility to export the graphs in a vector format for high-quality scientific graphs. Along the thesis quite a few function and script were developed but only the two batch treatment scripts will be briefly described here. The version 6.36 of Igor Pro was the last update installed during this thesis.

D.1 Velocity of experimental bonding

The first batch data treatment script is a script that use the exported data from the Matlab script, described in Appendix C, that locates the position of the bonding front on the recorded IR images.

The data are imported using the User Procedure → Main function. Fig. D.1 shows the window after the importation of the data and the selection of the steady state area for the average velocity calculation.

At first the Main function plots three velocity profiles, corresponding by default to the position of 1 /4, 1 /2 and 3 /4 of the total width of the imported data. Then the Cursor A and B are dragged and drop to select the steady state region. Finally clicking on the Refresh button on the small control windows display the steady state average velocity of the three profiles.

However sometimes the three default profiles are not well chosen, as the imported data often have some profile where the bonding front position was not correctly followed by the previous batch image treatment, or there was a bonding defect on the original sample. In this case another profile should be selected by using the up/down arrows to change the point number in the control window. To help choose a profile without defect and to have an idea of the homogeneity of the velocity on the sample the second plot of Fig. D.1 shows the average velocity value of all the profile (taken between the two cursor time).

The RtD button stands for Return to Default and the default position of all the windows.

D.2 Post-processing of Comsol ® data

The second batch data treatment script is a script mainly used to get the velocity value of the simulation results. It uses the three data set exported from the Comsol ® simulations (the position of the top beam u and its first and second order time derivative u t and u tt ). The importation of the data is done by clicking on Data → Load Waves → Load Comsol .csv files and selecting one of the three file exported with Comsol ® , the two other files will be automatically opened and read if their names respect the exportation names described in Sec. A.8. After the data importation, the function under User Procedure are used to treat the data.

First the Analyses → All function read from the data files all the information needed such as the time step, the width and number of points of the imported data. Then the Display Control Panel option display the control windows with all the information read previously. Finally the Plots → All function plots a few different plots; however the plot important to the determination of the bonding velocity is only the Position Front plot. On this plot the Cursor A and B are placed as wanted and the Fit button of the Fit Vitesse box of the control window is clicked to calculate the bonding velocity between the cursors, which is displayed in the Vitesse numerical field.

Abstract

Direct bonding is a process by which two sufficiently flat and clean surfaces can bond to each other without any added adhesive layer. Direct bonding of patterned surfaces is often used for the fabrication of Micro-Electro-Mechanical Systems (MEMS), where a silicon wafer with cavities is bonded to a plain wafer. The fabrication of these devices is expensive and it would be useful to have guidelines when designing knew devices to know in advance if direct bonding will be possible.

A 2D simulation model of the direct bonding of two substrates is developed and used to study the influence of the cavities on the bonding wave velocity. The prediction of the simulation run with Comsol ® are in good coherence with the experimental measures and a 2D law of the bonding velocity is obtained. The bonding of perfectly flat wafers with cavities should always be possible. Limitations to the bonding of real wafers are due to the elastic energy cost of deforming the non-perfectly flat wafers. This limit is reached easily when the bonding wave must cross a trench, so a design with a small bonding guide to help cross the cavity will work best. The width of this wave guide should be chosen by considering the bow of the wafer. Indeed the second important design rule is to keep a bonding area big enough to have more adhesion energy than the elastic energy cost due to non-flat wafers deformation.

The adhesion energy is an important parameter of the direct bonding, as it is the energy that drives the adhesion. This adhesion energy is different from the more widely known bonding energy which is the energy needed to separate two previously bonded wafers. In this work a simple method to measure the adhesion is proposed. Long time measurement of the evolution of the adhesion energy lead us to propose a mechanism for its evolution linked to the formation of capillary bridges between rough surfaces.