U Zaghloul

G J Papaioannou

F Coccetti

P Pons

R Plana

Effect of Humidity on Dielectric Charging Process in Electrostatic Capacitive RF MEMS Switches Based on Kelvin Probe Force Microscopy Surface Potential Measurements

Keywords: Dielectric charging, Humidity, Silicon nitride, Kelvin Probe Force Microscopy

In this article we investigate the effect of relative humidity on dielectric charging/discharging processes in electrostatically actuated MEMS devices. The assessment procedure is based on surface potential measurements using Kelvin Probe Force Microscopy (KPFM) and it targets in this specific work PECVD silicon nitride films in view of application in electrostatic capacitive RF MEMS switches. Charges have been injected through the AFM tip and the induced surface potential has been measured under different relative humidity levels. The impact of the charge injection duration and the bias level as well as bias polarity applied during the charge injection step, have been explored. Finally, the effect of the dielectric film thickness under different relative humidity levels has been assessed through depositing SiN films with different thicknesses.

INTRODUCTION

A very few work has been performed in order to study the effect of humidity on the dielectric charging process in RF-MEMS switches and therefore a full understanding of this effect is still missing. Surface charging has been reported to be more detrimental to switch lifetime than bulk charging [1]. Besides, based on measuring the shift in the switch actuation voltage, bulk charging was found to dominate in dry air, while surface charging was reported to increase linearly with increasing humidity [START_REF] Peng | Proc. Asia-Pacific Microwave Conference[END_REF]. Finally, it is reported that under comparable humidity levels and electric fields, switches employing silicon dioxide dielectric are less susceptible to surface charging than switches employing silicon nitride [1].

Recently, Kelvin Probe Force Microscopy (KPFM) has proven to be a very efficient method in assessing the dielectric charging process [START_REF] Lamhamdi | the 14 th International Conference on Solid-State Sensors[END_REF][START_REF] Herfst | the 46 th Annual International Reliability Physics Symposium[END_REF][START_REF] Belarni | [END_REF][START_REF] Zaghloul | MRS Fall Meeting[END_REF][START_REF] Zaghloul | the 15 th International Conference on Solid-State Sensors, Actuators, and Microsystems[END_REF]. For example, the surface potential distribution was found to be more confined in thinner dielectric films when the same bias is applied during charge injection, independently on the substrate nature [START_REF] Zaghloul | MRS Fall Meeting[END_REF]. Yet, all previous work done up to now which made use of the KPFM for dielectric charging assessment has been performed in ambient air and hence the effect of humidity and the surface charge have not been addressed. On the other sides, the removal of the water-related layer from both surfaces on the sample and the AFM tip is reported to be very important to improve the reliability of the KPFM measurements [START_REF] Ono | [END_REF]. At low humidity ≈10% RH trapped surface charges could be detected in the KPFM measurements while such charging effects become less prominent at high humidity ≈80% RH [9]. Finally, the adsorbed water layer on the sample surface is considered to shield the surface potential [10] This work therefore presents the effect of humidity and in sequence the surface charge on the charging/discharging processes in PECVD silicon nitride films based on KPFM methodology. The KPFM measurements have been performed under different relative humidity levels. Charges have been injected using different bias levels and various injection durations for both positive and negative bias polarities. Finally, the effect of the dielectric film thickness under different relative humidity levels has been investigated.

KPFM EXPERIMENTAL SETUP

In our KPFM experiments, charges are first injected in single points over the SiNx surface by applying voltage pulses of defined amplitude, Up, and duration, Tp, to the AFM tip (Fig 1 .a). Then the surface potential decay is measured with time. The analysis of the results is performed by measuring the full width at half maximum (FWHM) as well as the surface potential peak directly after charge injection for each injected point charge (Fig. 1.b). The effect of the charge injection duration, Tp, has been investigated through injecting charges for different durations in the range from 0.01sec to 100sec. Besides, the impact of the applied bias amplitude, Up, has been explored through applying Up pulses of 10V and 20V amplitude. Finally, the experiments have been performed for both positive and negative Up.

The surface potential measurements have been performed under two different relative humidity levels; in ambient air (RH ≈ 40%) and under dry air flow (RH ≈ 6%). The samples which have been measured under dry air flow have been first dehydrated at 120ºC for two cycles during three hours under vacuum just before the KPFM measurements. The SiN films were deposited by high frequency (HF) power supply (13.56 MHz) PECVD method over a Ti/Au layers which have been evaporated over low resistivity silicon substrate. In order to investigate the effect of dielectric thickness, two different film thicknesses have been investigated; 200nm and 300nm.

The KPFM measurements have been performed using the Digital Instrument Dimension-3100 equipped with Nanoscope IV controller. SCM-PIT conductive AFM tips from Veeco, which are Antimony n-doped silicon and coated with platinum/iridium, have been used. Here it should be mentioned that the measured surface potential varies from one tip to another due to the differences in geometry and intrinsic characteristics of each tip. Therefore, the same tip has been used for all measurements in order to be able to compare the results for different samples precisely. The employed tip-sample separation in our measurements was 20 nm. Moreover, the parameters used during the KPFM measurements were optimized using the scheme presented in [9]. Finally, the KPFM experiments have been performed several times for each sample under different charge injection conditions and the error percentage in the results presented in this work is kept less than 5%.

RESULTS AND DISCUSSION

Charges have been injected through the AFM tip in both contact and tapping modes. The impact of the relative humidity levels on the surface potential decay as well as on the lateral charge spreading for charges which have been injected in contact mode will be addressed in the following section. Then, the effect of the humidity on the surface potential distribution will be investigated for charges that have been injected in tapping mode for different duration and under different bias amplitudes.

Effect of humidity on surface potential decay and redistribution

For this set of experiments, charges have been injected in contact mode using a pulse amplitude, Up, of 10V and a pulse duration, Tp, of 60sec and the measurements have been performed for both the 200nm and 300nm SiN films. Figure 2(a,b) presents a comparison between the time dependence of the surface potential measured under different relative humidity levels, for both investigated SiN samples. For both investigated SiN thicknesses and under both 6% and 40% relative humidity levels, the surface potential is found to decay exponentially with time following the stretched exponential law,

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - β τ t exp
, where τ is the process time constant and β (0≤β≤1) is the stretch factor [START_REF] Zaghloul | the 15 th International Conference on Solid-State Sensors, Actuators, and Microsystems[END_REF]. As can be observed, the decay time constant in case of measurements under dry air is much longer than in case of measuring the surface potential in ambient air. In case of 6% RH , a KPFM signal of around 73% and 58% of the initial potential at t=0 for the 200nm and 300nm samples respectively, could be detected even after more than 11 hours of measurements. On the contrary, in case of measurements performed under 40% RH, the measured KPFM potential is around 17% and 13% of its initial potential at t=0 for the 200nm and 300nm samples respectively, after 35 minutes only. There are several reasons of possible charge decay, which are usually classified into two groups: external and internal charge decay. The external decay is due to the neutralisation of surface charges, which can be caused by the adsorption of counterions from the surrounding medium or by charge loss due to its electrical conductivity. In the case of air as environment the amount of available counterions and the electrical conductivity are both influenced by the humidity. On the other side, the internal decay is caused by the small but non-zero electrical conductivity of the dielectric material and by diffusion of charge carriers. Diffusion plays mostly a minor role in charge decay [START_REF] Sessler | Electrets[END_REF]. As surface charges are predominantly subjected to external charge decay and volume charges are predominantly subjected to internal decay which is not affected by a change of the environmental conditions, the huge observed difference in the decay time constant between measurements under 6% and 40% RH is attributed to the excess surface charge decay in case of higher relative humidity levels due to a possible thin water film over the dielectric surface. This water film has been considerably removed during the dehydration step under vacuum performed for the SiN samples measured under dry air flow.

Comparing the values of the decay time constants for the 200nm and 300nm SiN films for KPFM measurements performed under dry air (Fig. 2a) and in ambient air (Fig. 2b), it can be observed that the time constant is larger in the thinner dielectric films than in the thicker ones [START_REF] Zaghloul | the 15 th International Conference on Solid-State Sensors, Actuators, and Microsystems[END_REF], independently of the employed relative humidity. The longer time constants in case of thin dielectric films can be attributed to the fact that the thinner material will be affected more by the defective first deposited SiN layers. Such layers will contain more complex defects with large capture cross section located deep in the band gap. In the case of the thick dielectric film on the contrary the top deposited layers are expected not to be affected by these disordered bottom layers. In other words, the defects in these top layers are expected to exhibit simpler structure and have smaller capture cross section.

On the other side, the lateral charge redistribution of the injected charges has been investigated through monitoring the evolution of the potential distribution with time, Fig. 2 (c,d). No lateral charge spreading has been observed when the KPFM experiments have been performed under dry air with 6% RH even after 11 hours (Fig. 2c). On the contrary, lateral charge redistribution is clearly observed after 35 minutes only in case of surface potential measurements performed in ambient air with 40% RH (Fig. 2d), in spite of the low reported mobility of electrons in SiN [START_REF] Herfst | the 46 th Annual International Reliability Physics Symposium[END_REF]. This can be attributed to two main reasons. First, the dehydration step which has been performed for the SiN samples measured under dry air flow. Such annealing and cooling steps under vacuum are expected to remove a considerable part of the water film which already exists over the dielectric surface. Second, the KPFM measurements under dry air will prevent the water condensation over the dielectric surface. These two reasons will thereby inhibit lateral charge migration due to the electrical conductivity of the mentioned possible water film [START_REF] Knapp | [END_REF].

Here it should be pointed out that the evolution of the surface potential from KPFM measurements in ambient air showed that the charges practically do not diffuse across the surface during charge dissipation and the decay of the charge was attributed to transport and trapping into the bulk of the dielectric [START_REF] Lamhamdi | the 14 th International Conference on Solid-State Sensors[END_REF][START_REF] Herfst | the 46 th Annual International Reliability Physics Symposium[END_REF]. Besides, it should be emphasized that although lateral charge redistribution is observed in our KPFM measurements performed in ambient air (40% RH), this cannot justify the previously observed surface potential fast decay (Fig. 2b). Therefore, the relation between the lateral charge spreading and the surface potential decay is still under further investigation.

Effect of humidity on surface potential distribution

The results presented in this section are for charges which have been injected in tapping mode for different duration, Tp, different bias amplitude, Up, and with both positive and negative bias polarities (see Fig. 1). The investigated time range for Tp is from 0.01sec to 100sec while two different bias amplitudes have been used; 10V and 20V. The FWHM and hence the surface potential distribution resulting from charges which have been injected in different points, are measured directly after the charge injection step. Figure 3 presents a sample of these results where a row of single point charges separated by 4 µm each has been injected and measured under dry air flow (Fig. 3a) and in ambient air (Fig. 3b) using different pulse duration, Tp, and under pulse amplitude of 20V. The charge injection in different points has been performed using the NanoScript software available from Veeco, which controls the tip movements on the sample surface during the charge injection step. It is clearly observed from the figure that both the FWHM as well as the surface potential amplitude, Us, increase remarkably with increasing the pulse duration.

Figure 4 shows the FWHM versus the charge injection duration ,Tp, for the 200nm SiN sample measured under dry air (6% RH) and in ambient air (40%) for charges which have been injected for different pulse amplitude levels (10V and 20V) with both positive and negative bias polarities. It is clear from the figure that the FWHM increases with increasing the charge injection duration, Tp. This is found to apply for both 10V and 20V bias amplitudes, positive and negative bias polarities, and finally for measurements under both 6% and 40% RH levels. Besides, from the comparison between the FWHM for Up of 20V and 10V, at either 6% or 40% RH levels (Fig. 4a), it is noticed that as the applied Up increases, the FWHM gets wider and hence the surface potential distributions becomes less confined.

the effect of relative humidity on the surface potential distribution, Fig. 4a highlights that the FWHM is smaller in case of KPFM measurements performed under 6% RH than in experiments performed in 40% RH for both Up of 20V and 10V. In other words the measured surface potential becomes more confined as the relative humidity decreases. Moreover, the rate of increase of FWHM with the injection time is found to be more pronounced with higher bias at higher relative humidity levels. These observations have been found to apply as well when using negative charge injection, Fig. 4b. The impact of relative humidity on the FWHM and hence on the surface potential distribution can be directly related to the difference in the conductivity of the possible water film and hence the amount of the surface charge which exists on the dielectric surface as the relative humidity changes. Basically, the measured surface potential distribution follows the electric field distribution which exits during the charge injection step between the AFM tip and the SiN sample ground plane. As reported in [13], the distribution of the electric field just under the probe becomes less concentrated as the adsorbed water increases. Hence, as the conductivity as well as the amount of surface charge are expected to be smaller for the KPFM measurements performed at the 6% RH due to the reasons explained earlier, the electric field distribution and consequently the induced surface potential are expected to be more confined comparing to the KPFM measurements occurred in ambient air. Another contribution could be attributed to the expansion of the apparent tip curvature radius by the adsorbed water layer at higher relative humidity levels [10]. Comparing the FWHM values when positive or negative bias Up are applied to the AFM tip reveals that the surface potential distribution resulting from the negative charge injection is more confined than in case of applying positive bias, Fig. 4(c,d). Moreover, this has been observed for both the investigated bias levels, 20V (Fid. 5c) and 10V (Fig. 4d). This can be attributed to an asymmetric charging behavior which have been observed in both symmetric [14] and asymmetric [15] MIM capacitors. Moreover, comparing Fig. 4c and Fig. 4d reveals that the difference between the FWHM measured in 6% RH and 40%, is more pronounced as the applied bias gets smaller. As discussed earlier the FWHM increases as the applied bias, Up, increases. Therefore, for the 20V applied bias the FWHM is expected to be already large and then the charge spreading can be fairly enhanced by the relative humidity variation. While in case of the 10V bias, the FWHM is already small enough to be affected sharply by the variation in the relative humidity levels. In other words, the effect of changing the relative humidity on the FWHM is not very obvious being overwhelmed by the impact of the high bias applied.

CONCLUSION

The effect of the relative humidity on the dielectric charging process in PECVD silicon nitride films has been investigated based on the surface potential measurements from Kelvin Probe Force Microscopy. Two different humidity levels have been investigated; 6% under dry air flow and 40% in ambient air. Various parameters have been investigated which are the dielectric film thickness, charge injection duration and charge injection bias amplitude and polarity. The investigation reveals that the potential distribution was found to be more confined as the relative humidity decreases especially at smaller charge injection bias. The possible water film as well as the excess surface charge which exists at high relative humidity levels is found to contribute to the fast charge decay as well as the lateral charge redistribution. The surface potential was found to decay following the stretched exponential low for KPFM experiments performed under the two investigated humidity levels. Moreover, thinner dielectric films were found to have longer time constants, independently of the employed relative humidity level. Finally, asymmetry has been observed when positive and negative bias is applied during the charge injection.

Figure 1 .

 1 Figure 1. The procedure of the KPFM experiment and result analysis (a) charge injection in single point, (b) resulting surface potential distribution.

Figure 2 .

 2 Figure 2. Impact of relative humidity on surface potential decay and redistribution with time. Surface potential decay for SiN films of 200nm and 300nm thickness at (a) 6% RH and (b) 40% RH. Redistribution of the surface potential with time for the 200nm SiN film at (c) 6% RH and (d) 40% RH.

Figure 3 .

 3 The measured surface potential at (a) 6% RH and (b) 40% RH for charges injected in three different positions each for different charge injection time, Tp, (100sec, 10sec and 1 sec from left to right respectively), pulse amplitude, Up, is 20V and the SiN film thickness is 200nm.

Figure 4 .

 4 The measured FWHM under 6% RH and 40% for the 200nm SiN sample, versus the charge injection time, Tp, for (a) different positive pulse amplitude, 10V, 20V (b) different negative pulse amplitude, -10V, -20V, (c) different bias polarities 20V, -20V, (d) different bias polarities,10V, -10V.

ACKNOWLEDGEMENTS

This work has been partially supported by the following projects: the French ANR project FAME (PNANO-059), POLYNOE project funded by European Defense Agency (B-0035-IAP1-ERG), and the SYMIAE project funded by the Fondation STAE.


Adrien Broué
email: adrien.broue@novamems.com

Jérémie Dhennin

Cédric Seguineau

Xavier Lafontan

Christel Dieppedale

Jean-Michel Desmarres

Patrick Pons

Robert Plana

Methodology to Analyze Failure Mechanisms of Ohmic Contacts on MEMS Switches IRPS_ABroue_2009_final_NON corrumpu.doc Methodology to Analyze Failure Mechanisms of Ohmic Contacts on MEMS Switches

Keywords: electrical contacts, contact heating, failure mechanisms, switches, microelectromechanical systems (MEMS)

This paper demonstrates the efficiency of a new methodology using a commercial nanoindenter coupling with electrical measurement on test vehicles specially designed to investigate the micro contact reliability. This study examines the response of gold contacts with 5 µm² square bumps under various levels of current flowing through contact asperities. Contact temperature rising is observed leading to shifts of the mechanical properties of contact material, modifications of the contact topology and a diminution of the time dependence creep effect. The data provides a better understanding of micro-scale contact physics especially failure mechanisms due to the heating of the contact on MEMS switches.

INTRODUCTION

Research on contact degradation for microelectromechanical system (MEMS) switches has been driven by the necessity to reach a high-reliability level for micro-switch applications. The two main failure mecanisms observed during aging of the devices are the stiction of the free-standing element and the increase of the electrical contact resistance (Rc). The mechanisms leading to failure are plentiful and localized at the metal contact interface. The root causes of those issues stem from interactions between the deformation of the contact asperities, the current flow and the heating at the contact [START_REF] Patton | Fundamental studies of Au contacts in MEMS RF switches[END_REF] and lead to melting, material transport, adhesion, friction, wear, arcing, fracture, built-up of insulating films… The challenge is to set apart one event from each other.

The behavior of the current flow at contacting surfaces depends on several mechanical and operating parameters such as thermal (σ t ) and electrical (σ e ) conductivity of contact materials, mechanical properties of the contact spots (E and H), roughness (Ra), adhesion energy of contact surfaces (γ) , mechanical stiffness (k), humidity (Rh), contact voltage (V C ), contact force (F C ) (cf. Fig. 1).

Contact parameters leading to modifications of the topology of a micro contact surface

This study suggests a way to identify the dominant failure mechanisms by coupling specific tools and dedicated test vehicles under variant operational conditions. The purpose of this paper is to point out the usability of this experimental set-up to improve the understanding of the micro-contact physic to lead the way toward the high-reliable contact with a low electrical contact resistance.

CONTACT THEORY AND FAILURE MECHANISMS

First, the micro-contact physically differs from the macro-contact due to the roughness of the contact surface. Only high points on each surface come in contact thereby the effective contact area, named asperities or a-spots, is largely smaller than the apparent one. The constriction resistance mainly governs the contact resistance as the current flow is constricted through the small asperities leading to the electrical contact. For a single circular contact spot of radius a, Rc can be written as [START_REF] Holm | Electrical Contacts-Theory and Applications[END_REF]

R c = ρ 2a ( 1
)
where ρ is the resistivity of the contact material.

During the first contact establishment between the two surfaces, the load applied is generally higher than the yield stress of the contact material. Thus the deformation of the contact asperities is predominantly plastic. The contact area and the contact load can be linked to the radius of the contact spot a using Abbott and Firestone's plastic contact model [START_REF] Holm | Electrical Contacts-Theory and Applications[END_REF] π π

H F A a C C = = ( 2
)
where Ac is the contact area and H the Meyer hardness of the softer material. Due to the surface roughness, the effective contact resistance can be obtained as a first approximation by summing all the contact radii of all the individual asperities as the effective contact radius a eff of the effective contact area A eff . This effective contact radius can be substituted into (1) and ( 2) to take into account the multiple asperity feature of the micro contact.

Moreover the small contact force available in microswitches yields generally to contact spots with a size comparable to the mean free path in the contact materials. Because of that, mechanical, electrical and thermal properties shift from bulk to a specific physics of thin films lead by the geometry of the asperities. For example, the plasticity of the contact asperities (a-spots) in micro-contacts is rather ruled by the nucleation of dislocations than by their motion unlike macro-switches [START_REF] Fortini | Asperity contacts at the nanoscale: comparaison of Ru and Au[END_REF]. In the same manner, the heating of the contact spot is extremely localized when the current flows through the contact whereas the device level remains at room temperature [START_REF] Jensen | Effect of nanoscale heating on electrical transport in RF MEMS switch contacts[END_REF]. The contact spot temperature T c has already been expressed by Holm as a function of the contact voltage V c [START_REF] Holm | Electrical Contacts-Theory and Applications[END_REF]

0 4 2 T L V T c c + = ( 3
)
Where L=2.45 x 10 -8 W. /K² is the Lorentz constant and T 0 is the ambient temperature.

The second important point to be addressed is the reproducibility of the contact resistance over the actuations, and thus along its whole lifetime (cf. Fig. 2). The reliability of the switch depends on its ability to withstand some degradations occurring at the contact interface. Three types of phenomena can be studied: the mechanical (cold welding, strain hardening), electrical (arcing, hot welding, annealing) and chemical ones (formation of insulating films at the extreme surface), all inducing modifications of the topological, mechanical and electrical properties of the contact. As the experimental set-up developed here can be used to vary several parameters, the understanding of potential failure mechanisms can be made easier if the effort is put upon the identification of the key parameters for each mechanism. Moreover, these phenomena should be well described to explain any discrepancy between experimental measurements and the predicted values by the available models.

In addition, the contact temperature (Tc) of contact a-spots could be linked with the contact force depending on the contact voltage [START_REF] Kwon | Investigation of the electrical contact behaviors in Au-to-Au thin-film contacts for RF MEMS switches[END_REF]. Indeed the plastic deformation of the asperities during the contact formation proceeds more rapidly when the softening temperature is reached [START_REF] Fortini | Asperity contacts at the nanoscale: comparaison of Ru and Au[END_REF]. Thus the effective contact area increases inducing a drop of the contact resistance. However the softening of the metal at the asperities of contact reduces the strain hardening of the a-spots and could accelerate the aging of the contact by the activation of thermal failure mechanisms. Therefore the softening temperature of contact metals is probably one of the most important critical parameter for a stable and low contact resistance.

EXPERIMENTAL METHODS

Description of the experimental set-up

The test described in this paper is based on methods combining a nanoindenter and a high resolution source meter for the determination of the electrical contact resistance versus the contact force and the displacement of the free standing part of the contact. The actuation load is thus reproducible and known with a good accuracy (10 -9 N and 10 -10 m), two points which remain hard to achieve with the actual actuations (for example, electrostatic actuation forces often drift due to dielectric charging). The nanoindenter's spherical tip is used as a mechanical actuator and can be applied on the sample with a 20nm accuracy positioning thanks to a piezoactuated XY table. Thus a pseudo-AFM scan is performed by probing the surface with a fixed stiffness of contact between the tip and the surface (cf. Fig. 3).

The 4-wires electrical measurements are coupled with the nanoindenter software in order to control the level of load and the displacement of the tip by taking into account variations of the electrical resistance of the contact (cf. This experiment is also useful for studying hot switching or cold switching modes. In this way we are able to discriminate the failure mechanisms induced by the current or by the mechanical cycling of the contact. Moreover, the level of deformation sustained by the contact asperities (elastic, elasto-plastic or plastic) can be pointed out, depending on the applied load and actual contact area.

Test vehicle description

Specific test vehicles have been designed for contact analysis, enhancing the extraction of characteristic curves, and making possible the comparison between different contact shapes or materials. As illustrated in Fig. 5, the testing device is composed of a bridge suspended over the contact line. A 3µm square bump is processed underneath. The stiffness of the bridge is about 480 N/m. In addition, the test structures are designed to measure only the contact resistance as the resistance of the line is canceled by the 4-wire measurement process. The material contact is Au-to-Au because gold is the most popular material for electrical contact due to its high conductivity and its high resistance to oxidation.

TEST DESCRIPTION AND RESULTS

A campaign of tests has been made on these test vehicles to study the self heating of the contact with increasing level of current flowing through the contact. These performance tests enable us to compare the temperature of the contact with the analysis of the contact resistance versus the contact force up to 150 µN by sweeping the contact current from 1 mA to 100 mA and measuring the contact voltage in a cold-switching way. Moreover, the holding plateau at load max is fixed at 20s to study the creep of the contact. The load rate is about 1µN/s to allow the temperature to reach steadystate. The following results are based on the same 11 successive tests and are presented to emphasize the various observations allowing by this experimental setup. Fig. 6 shows typical contact force-contact resistance curves for a test structure with a 3µm square dimple at 11 levels of contact current. This curves point out the relationship between the contact force and normal load. The contact resistances decrease with the increasing contact force as a function of the contact current and are described by the equation

Sensitivity to current increase

x C C AF R - =
where A and x are constant parameters. The parameter x decreases from 0.46 to 0.07 with contact current rising and does not fit with (2) because this equation does not take into account the roughness of the contact surface. The second observation is that the higher the current is, the lower is the dependence of the contact resistance to the load applied. The initial contact resistance remains The origin of this contact resistance reduction are thermally induced with the increase of the contact temperature with the current until reaching the softening temperature. Thus the annealing of the contact leading to a reduction of the contact hardness and the roughness of the contact is lower after softening of the surfaces. Thereby, the deformation of the contact asperities is higher with increasing current making easier the breakdown of insulating films, and enlarging the total contact area.

Contact temperature

Fig. 7 shows the contact temperature as a function of the contact current. The published softening temperature for gold contact is ~100°C, corresponding to a contact voltage of 70-80 mV for contact near room temperature [START_REF] Jensen | Effect of nanoscale heating on electrical transport in RF MEMS switch contacts[END_REF]. We observe that when the contact current is increased after the first contact softening from 40 mA, the contact resistance continues to decrease keeping the contact temperature roughly constant. The contact temperature increases with constant slope from 1mA to 40mA. At this point the contact temperature stops to increase and fluctuates between 80°C and 120°C (65 mV to 75 mV). Thus when the softening temperature is reached, the contact temperature does not depends strongly on contact current at high current level.

Topological modification of the contact surface

As shown in Fig. 8, the number of asperities n, can be evaluated by the way of the mechanical measurements provided by the methodology. The stiffness of the structure can be measured by evaluating the load versus displacement slope and the occurrence of a new mechanical contact can be observed: A slight increase of this value is due to the deformation of the asperity of contact. When the contact force is high enough, new asperities come into contact, inducing a slight sudden modification of the loading configuration, that is to say a discontinuity in the stiffness recorded. We can observe that the effective contact area spreads with the number of the a-spots in contact while contact loads increase, while the contact resistance decreases. The total number of a-spots coming into contact also increases with higher contact current but the cross sectional area of a-spot also increases as they are pressed down. These two phenomena are in competition and it is hard to set apart one from each other. However we can conclude that the softening of the contact material leads to an increase of the number of contact a-spots.

A representation of the roughness of the contact bump in the test vehicle is obtained from AFM measurements (Fig. 9). The first picture shows the contact bump of a structure untested coming from the same run as the test used for the test campaign previously described. The second picture is the AFM imaging of the bump after the 11 successive tests (Fig. 6). We can observe the flattening of the asperities due to plastic deformation of the a-spots during the loading of the contact. The measurements point out that the rms roughness of the top of the contact bump decrease from 18.1 nm to 13.8 nm. The ratio of the real contact area for each current level to the initial contact area for Ic = 1mA is evaluated from (1) and ( 2 We observed that the contact area increases with a slope which is steeper from 50 mA. Indeed the inflexion point of 50 mA matchs with the rise of the contact temperature to the softening temperature. The plastic deformation is higher when the softening temperature is reached leading to a larger contact area in relation with a larger topographical modification of the contact surfaces. The variation of the contact resistance is recorded during the holding plateau for the 11 tests at load max (~143 µN) for each current level. The results are shown in the Fig. 11. The initial deformation is enough to support the load max applied on the contact and time dependant creep effects lead to additional asperity deformation during the holding plateau.

Contact creep with increased contact current

The variation of the contact resistance is directly linked to the contact area variation. We observed that the creep of the contact asperities is lower at high current although the contact material is softer. It could be explained by the fact that high current density in the asperities leads to a quick establishment of large contact area. Thus the larger contact area leads to a reduction of applied pressure on each a-spot and prevents to creep due to less stress at asperity contacts [START_REF] Patton | Fundamental studies of Au contacts in MEMS RF switches[END_REF].

CONCLUSIONS

This new test facility had enabled to perform tests under realistic conditions in a static and dynamic way.

An emphasis was placed on the role of the current increase leading to contact heating in the understanding of the contact physics. Various phenomena were observed with asperity heating effects on contact resistance versus contact force dependence, stabilization of the contact temperature from the softening temperature of the contact material, modification of the surface topography, and diminution of the creep effect. All this mechanisms are commonly report as failure vectors in the contact switch lifetime.

And finally, this tests bring to light the usability of this methodology to investigate the contact reliability and performance under several test conditions. The test vehicles used in this study could be manufactured with several contact materials with various bump designs. The investigation on these two conception parameters governing the robustness of the contact and the contact geometry, which has a strong influence on the heat distribution across the contact, is underway. The knowledge gained will be the identification of the failure mechanism for each contact design in order to bring recommendations for further improvements of the contact reliability.

Figure 1 .

 1 Figure 1. Typical lifetime of micro switches leading to contact failure

Figure 2 .

 2 Figure 2. Picture of the scan performed with the nanoindenter on the test structure

  Fig 4).

Figure 3 .

 3 Figure 3. Principle of the electrical test performed with a nanoindentation tip

Figure 4 .

 4 Figure 4. Scheme of the test structures

Figure 5 .

 5 Figure 5. Contact resistance versus contact force as a function of the current flowing through the contact

  cycle at each increasing step of current.

Figure 6 .

 6 Figure 6. Contact temperature versus the current flowing through the contact asperities

Figure 7 .

 7 Figure 7. Evaluation of the number of asperities: number of steps on the stiffness versus contact force as function of the contact current

Figure 8 .

 8 Figure 8. Roughness measurements performing by AFM with nontested bump (Ra = 18.1 nm and RMAX = 114 nm) (a) and tested bump (Ra = 13.8 nm and RMAX = 90.0 nm) (b)

  ) at load max (Fig 10).

  Figure 9. Contact area ratio (

Figure 10 .

 10 Figure 10. Variation of the contact resistance during versus contact current during the holding plateau at FC ~ 143 µN for t = 20S

µm 2 µm

ACKNOWLEDGMENT

This work was supported by the ANR (Agence Nationale de la Recherche) in the frame of the FAME project. The authors would also like to thanks the CEA-LETI for the manufacturing of the test vehicles.


H Achkar

D Peyrou

F Pennec

K Yacine

A Ferrand

P Pons

M Sartor
email: marc.sartor@insa-toulouse.fr

R Plana
email: plana@laas.fr

Use Matlab Code to Link and Pilot Softwares Used for MEMS Simulation

This paper presents the advances done in the field of simulation and design applied to MEMS. It will be discussed some techniques to link different finite element softwares, especially COMSOL and ANSYS, in order to profit from all the advantages offered by each finite element code. Then, an explication is done on how we piloted ABAQUS using a Matlab routine. In a third and very important section, the use of reverse engineering method is described and applied to RF-MEMS (Radio Frequency-Micro Electro Mechanical System) switches. The applications for each of these methods will be explained.

Introduction

The MEMS are micro electromechanical systems which can do different jobs such as sensing, actuating or doing an electrical function. More specifically, RF MEMS characteristics have demonstrated during the last ten years a very attractive potential to allow the introduction of intelligence in the RF front end architecture [START_REF] Rebeiz | RF MEMS Switches and switch circuits[END_REF] through analogue signal processing techniques. Nevertheless, those micro devices with moveable mechanical structures still have some issues to be solved before they can be successfully introduced in the industrial level and available in the market.

First the reliability issue that is related to the type of degradation and to the actuation methods. Until nowadays, cantilevers and membranes are actuated using electrostatic, thermal, magnetic and/or piezoelectric energy with all what every type of actuation brings in terms of advantages and inconvenient. It is to be mentioned that today, even that most of the people are turning to piezoelectric actuation since it's actuated at a very low voltages, electrostatic actuators offer the best trade-off [START_REF] Tschulena | Market analysis for MEMS and Microsystems III, 2005-2009. A Nexus task force[END_REF] [3] [4] [START_REF] Van Sprengen | A comprehensive model to predict the charging and reliability of RF-MEMS switches[END_REF] [START_REF] Yuan | Acceleration of dielectric charging in RF-MEMS capacitive switches[END_REF] this is if we exclude hybrid (electrostatic and piezoelectric) actuation.

Second is the MEMS process flow simulation which is until these days not complete and need to be developed more in order to become reliable and precise. This process flow simulation is very important since it drives the development of new devices and circuits. For that issue, a good knowledge of the material properties and specifically the mechanical properties is needed.

Third issue comes from the lake of efficient, complete and easy to use simulator covering the totality of the MEMS design procedure, from individual MEMS component design to complete system simulation. On the component level, finite element analysis (FEA) methods, offers high efficiency and are widely used in modelling and simulating the behaviour of MEMS. Since MEMS are multiphysics coupling thermal, mechanics, piezoelectricity, electrostatics, initial stress, mechanical contact and electromagnetic, it is difficult to find software competitive in all the fields.

As for us, we are interested to ameliorate the reliability of our components and to develop the simulation platform in order to reduce the trial and error loops in the V shape conception flow shown in figure 1. Reducing the number of experimental steps reduces the cost, while developing the simulators and the simulation tools ameliorate the design phase and by consequence reduces the time to market.

In this paper we will explain our contribution in each of the above issues using the developed simulation techniques. At the beginning we will explain the work done on piloting software using Matlab, especially Abaqus (COMSOL is already piloted by Matlab and we are trying to apply it to ANSYS). This work came as a result of the need in material characterisation, specifically in nanoindentation tests where number of simulations on Abaqus is needed and each model depends from the results obtained from the previous simulation and its offset from the experimental results. A Matlab code will compare experimental results to the simulation; modify the parameters in our basic model, then runs a new simulation. This code can be used to optimize structures in Abaqus. Concerning reliability, the reverse engineering method is described. This method enabled us to regenerate the real surface of two metals in contact in order to study the impact of the surface roughness on the failure of the component, it will also permits us to approach from the reality of the contact in capacitive switches to analyse the good quality of the capacitance in the down state and to study the effect of the initial flatness of the bridge on the actuation voltage and the initial deflection. Finally, since there is no software expert in all the fields at a time we have developed, always using Matlab for programming, a method to link the softwares that we find useful for us that are ANSYS and COMSOL. An ongoing work in our lab is studying how to link HFSS to other thermal finite element software.

Piloting ABAQUS

This section is divided into two parts; the first subsection will explain briefly the technique of nanoindentation test, its interest and the problem that alerted a necessity to control the model file in ABAQUS as well as the result file. In a second subsection we will explain the method and the codes we used to solve the problem.

Nanoindentation technique and problem

The basic idea of the nanoindentation test is to bend a micro structure, usually a membrane that is made of the material that we want to characterise and elaborated under known fabrication conditions. [START_REF] Xiaodong | A review of nanoindentation continuous stiffness measurement technique and its applications[END_REF] [START_REF] Staussa | Determining the stress-strain behaviour of small devices by nanoindentation in combination with inverse methods[END_REF] Geometry model 3D

Fabrications & Tests

No Yes Figure 1: V shape design flow, the time taken to finish this cycle is called time to market.

A nanoindenter equipped with a precise system to measure the deflection and the force at the center of the tested membrane will give us the experimental characteristic curves as figure 2.b shows. A characteristic curve (figure 2.a) of a membrane is the graph of the displacement at the centre as a function of the applied force at the same place. This characteristic curve obtained experimentally will be fitted to simulated results in order to bring out the mechanical characteristics such as Young's module, initial stress in the membrane and the ultimate stress. Since the membranes are very delicate and have a thickness of the order of 1µm and the area of the surface not more than an mm 2 , we need precise equipment for measurement in order to do an accurate fitting of the experiments with simulations. The problem of fitting is that its time consuming due to iterations and the long simulation time.

Pilot ABAQUS using Matlab

After experiments, we usually run simulations with different values of the parameters to be characterized, which are Young modulus and the initial stress, in a way to bound the experimental curve with two numerical curves, as in figure 3.a, then interpolate the values of the parameters. The initial stress is calculated from the initial slope of the characteristic curve since for small deflections the stiffness mostly depends on the stress state in the membrane [START_REF] Espinosa | Identification of Residual Stress State in an RF-MEMS Device[END_REF]. The problem is here; a small error on the value of the initial stress will change the shape of the characteristic curve for large displacement and induces an additional error on the measurement of Young modulus which is measured in the large displacement zone. This error due to initial stress comes from a bad fitting of the curves using the classical way of interpolating the Young's modulus from the two envelope curves. A Matlab code was wrote for this issue in order to fit well the two curves using optimization techniques and to reduce the time of convergence of the model to the experiments since additional iteration costs us the time of an additional simulation.

The code used for this service consists of sub routines, called functions, that the fitting loop will call one after the other in order to analyses the data and act on the model to modify it. Figure 4.a shows the flowchart of the code while figure 4.b is the graph of convergence for the fitting algorithm used. The idea is to automate the process of characterization with the minimum intervention of the user. After running the program, we have to enter the ABAQUS model file name then once the simulation is turned on, the user will be asked about the algorithm of fitting. For the fitting, we programmed two algorithms, a precise one that we called complete algorithm but it is time consuming and a fast algorithm that we called fast fitting but less precise. For the complete algorithm consisting of the minimization of the distance between two curves using the probabilistic method issued from monte carlo method called "simulated annealing", we converged to the experimental curve with 5 iterations and a final relative error of 0.86%. The fast fitting consisting of only 4 steps method that varies in the first 2 steps the first parameter (initial stress calculated from the initial stiffness of the characteristic curve) and in the next 2 steps the second parameter (Young's modulus), we converged after 4 iterations with a relative error of 0.71%. These results showed that the fast algorithm was better than the complete one since both methods converged with a relative error less than 1% which is excellent while it took more time for the complete method. In terms of results, even though the relative error for both methods was almost the same, the values of the parameters were different, of about 5% on the initial stress and of 13% on Young modulus. We kept the choice of the algorithm for the user who knows better which is nearer to the reality.

Reverse Engineering method

The basic idea is really simple; reverse engineering is a method to create a 2D or 3D virtual model of an existing physical part for use in CAD software. The reverse engineering process involves measuring the real object and then reconstructing it as a finite element model with the real geometry. The object can be C f C r measured using any scanning technology like optical profilometer or Atomic Force Microscopy (AFM) or mechanical profilometer. The measured data alone, usually represented as a point cloud, lacks topological information and is therefore processed and modelled into a more usable format such as a stereolithography format which is an ASCII file used in standard CAD softwares.

Application in capacitive switches

In capacitive switches, the value of the down state capacitance plays an important role in the definition of the passive frequency band and the isolation of the switch. Usually, the measured value of the down state capacitance is lower than that calculated theoretically by considering a zero air gap between the electrodes and a perfect contact between the two surfaces. The value of this capacitance is affected by the surface roughness and the dielectric charging [START_REF] Palasantzas | Surface roughness influence on the pull-in voltage of microswitches in presence of thermal and quantum vacuum fluctuations[END_REF]. We will model analytically the real capacitance by adding a capacitance with air in series to the theoretical capacitance in order to consider the air gap effect in our calculation as shown in figure 5.

Beam profile under electrostatic actuation up to pull in voltage

Contact between the beam and the dielectric (pull in time)

The down state capacitance is modelled by two capacitances in series :

Application to calculate capacitance

This analytical model is validated using numerical technique for simulation. The measurements done on the surface roughness of the dielectric and the bridge in the contact zone were regenerated into a 3D model as shown in figure 6. We used an optical profilometer (VEECO) to capture three-dimensional (3D) data points of the dielectric surface. Then, using Matlab functions we insert the generated surface as the top surface of a block. For CAD modelling, the segmented data are further transformed into individual surfaces. Several mathematical schemes for representing geometrical shapes are available. We choose to convert the closed surface into stereolithography format because this is an ASCII file used in most standard CAD softwares. Moreover, generated solid model are described by a list of the triangular surfaces which can be easily implemented on Matlab. The last step is to import the generated solid into CAD software (COMSOL) to make the down-state capacitance analysis.

Figure 6 illustrates the full method and data flow : dielectric roughness is scanned by the optical profilometer (A) to generate an ASCII file (B) representative of 3D coordinates points of the surface. Then solid model (C) and object file (E) as STL, VRML or Comsol geometry can be obtained using Matlab. To enhance the shape, one can use (D) some filters (filter2 for MatLab or a CAD software, CATIA for example). After this, the model is ready to be imported to finite element software (F). We used multiphysics software provided by COMSOL, to set up the electrostatic model (materials and boundaries conditions (F) -meshing (G), results and post processing (H)). The value of the capacitance obtained by simulation and the one obtained analytically are about 4% which is very good while it showed far values when compared with the measurements. This difference can be explained with the charging effect of the dielectric and neglecting of the deformation of the asperities during the contact. The table below summarizes the results obtained.

Simulated

Calculated Measured Theoretical 0.999335 pF 0.9571255 pF 1.41 pF 3.03 pF

Application to calculate resistance

The reverse engineering method is also used to simulate the contact between two rough surfaces in order to calculate the contact area and by that deduce the contact resistance. The methodology used to regenerate the model is the same as for the capacitance simulation but this time we will model the roughness in contact with a flat plate and using ANSYS. The contact resistance is a very important issue to solve in DC contact RF MEMS. It has been identified that most of the limitations are related to the quality and the repeatability of the contact that drive the RF performance (insertion losses, isolation, power handling) and the reliability. This numerical method is used to predict the real contact area as a function of the contact force between rough surfaces starting from real shape of the surfaces that come in contact, and then use analytical expressions to extract the electrical contact resistance. The novel approach relies on a reverse engineering method to generate the real shape of the surface. The post-processing generates the distribution of the contact pressure on the contact surface as shown in figure 7.b. From the pressure distribution and size of each contact spot it is possible then to use an analytical expression to extract the electrical contact resistance.

Linking Ansys to COMSOL

To have a complete simulation platform in order to predict in a good way the behaviour of a design, we are interested to link different software in order to profit from the advantages offered by each one of them. We will explain in what follows the method flow to link COMSOL to ANSYS using Matlab.

Explication of the linking method

As already stated, we need to link ANSYS and COMSOL in a way to exchange deformed geometry and treat it in the convenient software. This section describes the full method in order to regenerate geometries starting from the deformation obtained from a mechanical simulation. Finally meshing either the air between the membrane and the plane or meshing the membrane itself depending on the application.

COMSOL 3.3 appeared as software that is capable to do piezoelectric simulations combined to mechanical structure and contact problems, without forgetting thermal and initial stress effects. The problem appears here, when we need to simulate complicated contact problems and where COMSOL solver fails most of the time.

When talking about contact in COMSOL, it is somehow critical and difficult to converge the simulation specially that we have to vary the contact parameters and it's still not sure to converge. The convergence of the contact model using COMSOL failed; one have then to find a way to link both ANSYS and COMSOL in order to use ANSYS contact models. By exchanging the deformed geometry between them and then running the convenient simulation using the appropriate tool, some simulation problems can be solved.

A very practical way to run simulations on two different softwares, is to pilot it using Matlab. To achieve our purpose, we created in Matlab some functions which make the piloting and the link between the softwares easier. In what follows, the flow chart shows the key steps of the method by using a simple membrane of 200x200µm 2 subjected to constant pressure and blocked at its corner. We can begin the cycle from where we want, either from ANSYS or from COMSOL. In the flow chart, we began from COMSOL, with a simple membrane under pressure as explained in figure 8. The geometry is drawn, and then the loads and boundary conditions are applied. A mechanical simulation is run in order to obtain the deformation of the membrane. Considering that a contact analysis is needed after the deformation of the membrane, so we export the deformation of the membrane into a text file. In the text file there will be the displacement of each node of the chosen surface of interest as well as some data that we don't need. Matlab will treat this file to clean unneeded information and to reorganize data in the file. The treated file will be used then to create a cloud of keypoints, by creating it one by one, and then generating surfaces from the keypoints. In ANSYS now, we will mesh the regenerated geometry, we will define the physics that we need to study (in our case it's the contact) and then run a new simulation. The deformation is once again output in a text file, giving the displacement value for each node of the surface of interest. Another treatment and organization of the output data file, before being injected in the function that we created in Matlab which regenerates the surface and export it to COMSOL. The geometry is ready to be used in COMSOL, new objects can be added to it, for example electrode, before going into meshing and then boundary condition, to finish into a new simulation. The loop can turn as much times as we need, passing from software to another, until we attain our goal.

Conclusions and perspective

This paper showed the modelling techniques that we used to advance the MEMS simulation and to ameliorate the reliability of the components by characterizing the material with higher precision and modelling the contact state of the switch using the reverse engineering method.

Piloting ABAQUS permitted us to save time in nano-indentation tests since the curves fitting method is done automatically and it increased the precision of the characterization. As perspectives, the pilot of ABAQUS can be integrated on a production line of a clean room in order to automatically characterize the material. A profilometer can detect the dimensions of the membrane in order to create the first model for the fitting loop, then the nano-indentor will be automatically positioned to perform the test, the experimental results will be used then by the program in order to give the mechanical characteristics.

The reverse engineering method was described and it was used to simulate the capacitance in the down state of a capacitance in order to validate the analytical formula for the air gap effect on the capacitance loss. It was also applied to contact simulation in order to extract the electrical contact resistance. Simulating the contact resistance replaced the expensive tests that needed to be done for different materials in contact and for different surface roughness. This method will be used later to simulate the effect of the real initial shape of the bridge on the deformation due to initial stress and on the actuation voltage as well as the capacitance quality in the down state.

Finally we showed a method to link COMSOL to ANSYS in order to profit of the advantages offered by different simulators. In a previous paper we did a study on the numerical platform validation and we

Figure 2 Figure 3 :

 23 Figure 2: a) Experimental characteristic curve in a nanoindentation test. b) Experiment apparatus.

Figure 4 :

 4 Figure 4: a) Flowchart of the fitting routine. b) Graph showing the relative error with respect to the number of iterations for the complete method.

•Figure 5 :

 5 Figure 5: in the top figure we show the switch in contact with the dielectric and below is the model.

Figure 6 :

 6 Figure 6: The description of the data flow for the reverse engineering method.

Figure 7 :

 7 Figure 7: a) The regenerated model in ANSYS for contact resistance simulation. b) Simulation results.

  Figure8: flow chart of linking ANSYS to COMSOL.

  Starting from a drawn geometry in any of the softwares, a mechanical simulation is done.The mechanical deformation is written in an output file in .txt format.Data treatment in Matlab, and rearrangement of the outputs.From a set of point cloud, we create surfaces that will form the deformed geometry.Contact analysis or any other application (ex: CFD, …) is done in ANSYS. Once more, the deformation information is written in a file .txt.Data treatment in Matlab and rearrangement of the outputs. Regenerated geometry with possibility to add parts to the deformed shape.

offered a table that lists the advantages and drawbacks of the softwares, this will help understand why we needed the link. In a future step, we seek to do an automatic link between the softwares and get rid of the user's assistance in order to reduce the lost time. A new need in our group of work is to link the thermal results obtained from HFSS with ANSYS or COMSOL.


C Deledalle

S Vaiter

G Peyré

J Fadili

C Dossal

Proximal Splitting Derivatives for Risk Estimation

This paper develops a novel framework to compute a projected Generalized Stein Unbiased Risk Estimator (GSURE) for a wide class of sparsely regularized solutions of inverse problems. This class includes arbitrary convex data fidelities with both analysis and synthesis mixed ℓ 1 -ℓ 2 norms. The GSURE necessitates to compute the (weak) derivative of a solution w.r.t. the observations. However, as the solution is not available in analytical form but rather through iterative schemes such as proximal splitting, we propose to iteratively compute the GSURE by differentiating the sequence of iterates. This provides us with a sequence of differential mappings, which, hopefully, converge to the desired derivative and allows to compute the GSURE. We illustrate this approach on total variation regularization with Gaussian noise and to sparse regularization with poisson noise,o automatically select the regularization parameter.

Introduction

This paper focuses on unbiased estimation of the ℓ 2 -risk of recovering an image f 0 ∈ R N from low-dimensional noisy observations y = Φf 0 + w, where w ∼ N (0, σ 2 Id P ). The linear bounded imaging operator Φ : R N → Y = R P entails loss of information so that P < N or is rank-deficient for P = N , and the recovery problem is typically ill-posed.

In the following we denote f (y) ∈ R N the estimator of f 0 from the observations y ∈ Y. More specifically, we consider an estimator f (y) defined as a function of coefficients x(y) ∈ X (where X is a suitable finite-dimensional Hilbert space) that solves x(y) ∈ argmin x∈X E(x, y) where the set of minimizers is nonempty. Here E(x, y) is an energy functional parameterized by the observations y ∈ Y. In some cases (e.g. total variation regularization), one directly has f (y) = x(y), but for sparse regularization in a redundant synthesis dictionary, the latter maps coefficients x(y) to images f (y). We then make a distinction between f (y) and x(y) in the following.

This work proposes a versatile approach for unbiased risk estimation in the case where x(y) is computed by proximal splitting algorithms. These methods have become extremely popular to solve inverse problems with convex non-smooth regularizations, e.g. those encountered in the sparsity field.

Previous Works

Unbiased Risk Estimation. The SURE [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] is an unbiased ℓ 2 -risk estimator. For denoising Φ = Id, it provides an unbiased estimate SURE(y) of the risk E(||f (y)f 0 || 2 ) that depends solely on y, without prior knowledge of f 0 . This can prove very useful for objective choice of parameters that minimize the recovery risk of f 0 . A generalized SURE (GSURE) has been developed for noise models within the multivariate canonical exponential family [START_REF] Eldar | Generalized SURE for exponential families: Applications to regularization[END_REF]. In the context of inverse problems, this allows to compute the projected risk E(||Π(f (y)f 0 )|| 2 ) where Π is the orthogonal projector on ker(Φ) ⊥ . Similar GSURE versions have been proposed for Gaussian noise and special regularizations or/and inverse problems, e.g. [START_REF] Pesquet | A SURE approach for digital signal/image deconvolution problems[END_REF][START_REF] Vonesch | Recursive risk estimation for non-linear image deconvolution with a wavelet-domain sparsity constraint[END_REF]. SURE and GSURE have been applied to various inverse problem, to estimate the reconstruction risk E(||Φ(f (y)f 0 )|| 2 ), for wavelet denoising with linear expansion of thresholds [START_REF] Luisier | The SURE-LET approach to image denoising[END_REF], wavelet-vaguelet non-iterative thresholding [START_REF] Pesquet | A SURE approach for digital signal/image deconvolution problems[END_REF], synthesis sparsity [START_REF] Dupé | A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF][START_REF] Vonesch | Recursive risk estimation for non-linear image deconvolution with a wavelet-domain sparsity constraint[END_REF][START_REF] Giryes | The projected GSURE for automatic parameter tuning in iterative shrinkage methods[END_REF] and analysis sparsity [START_REF] Deledalle | Unbiased risk estimation for sparse analysis regularization[END_REF].

Unbiased Estimation of the Degrees of Freedom. A prerequisite to compute the SURE or GSURE is an unbiased estimate of the degrees of freedom df(y). Roughly speaking, for overdetermined linear models, df(y) is the number of free parameters in modeling f (y) from y. There are situations where df(y) can be estimated in closed-form from f (y). This occurs e.g. in synthesis ℓ 1 regularization, as established in the overdetermined case in [START_REF] Zou | On the "degrees of freedom" of the lasso[END_REF], and extended to the general setting in [START_REF] Kachour | The degrees of freedom of penalized l1 minimization[END_REF]. These results have been extended to analysis sparsity (for instance total variation) [START_REF] Vaiter | Robust sparse analysis regularization[END_REF][START_REF] Tibshirani | The solution path of the generalized Lasso[END_REF].

When no closed-form is available, df(y) can be estimated using data perturbation and Monte-Carlo integration, see e.g. [START_REF] Ye | On measuring and correcting the effects of data mining and model selection[END_REF]. Alternatively, an estimate can be obtained by formally differentiating the sequence of iterates provided by an algorithm that converges to f (y). As proposed initially by [START_REF] Vonesch | Recursive risk estimation for non-linear image deconvolution with a wavelet-domain sparsity constraint[END_REF] and refined in [START_REF] Giryes | The projected GSURE for automatic parameter tuning in iterative shrinkage methods[END_REF], this allows to compute the GSURE of sparse synthesis regularization.

Proximal Splitting Algorithm. Convex optimization problems that appear in imaging applications usually enjoy a lot of structure, that is efficiently captured using proximal splitting methods. These methods are tailored to tackle large non-smooth convex optimization problems and are now mainstream in image processing. See for instance [START_REF] Combettes | Fixed-Point Algorithms for Inverse Problems in Science and Engineering[END_REF] for an overview of these methods. The precise algorithm to be used depends on the specific structure of the problem. The forward-backward (FB) algorithm handles the sum of a smooth and a simple function (for which the proximal mapping can be computed in closed form), see for instance [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. The Douglas-Rachford (DR) algorithm [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] does not make any smoothness assumption, but requires that the functional is split into a sum of simple functions. The generalized foward-backward (GFB) algorithm [START_REF] Reguet | Generalized forward-backward splitting[END_REF] is an hybridization between FB and DR, thus enabling to take into account a smooth function and an arbitrary number of simple functions. Lastly, let us mention primaldual schemes, such as the one recently proposed by Chambolle and Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], that enables to minimize compositions of linear operators and simple functions. In this paper, we illustrate the versatility of our method by applying it to both primal (FB, DR and GFB) and primal-dual (CP) algorithms. The proposed methodology can however be adapted to any other proximal splitting method.

Contributions

The main contribution of this paper is a a new risk estimation for arbitrary sparse regularizations when the solution is computed using a proximal splitting algorithm. This extends the previous iterative computation methods [START_REF] Vonesch | Recursive risk estimation for non-linear image deconvolution with a wavelet-domain sparsity constraint[END_REF][START_REF] Giryes | The projected GSURE for automatic parameter tuning in iterative shrinkage methods[END_REF] to a broader class of regularizations (e.g. total variation) without the need to resort to Monte-Carlo integration (such as [START_REF] Lusier | [END_REF]) which are numerically costly.

Differentiating Iterative Schemes

Iterative minimization scheme. The computation of x(y) as a minimizer of some energy E(x, y) is usually solved using an iterative scheme, initializing x (0) (y) (for instance to 0) and then updating

x (ℓ+1) (y) = ψ(x (ℓ) (y), y), (1)
where ψ : X × Y → X is a mapping, and so that x (ℓ) (y) converge to a solution x(y) that satisfies the (non-expansive) fixed-point equation

x(y) = ψ(x(y), y). (2)
Note that we have made here explicit the dependency of the iterates x (ℓ) = x (ℓ) (y) with the observation y.

A simple, but instructive example, is the case where x → E(x, y) is a C 1 function with 1/L Lipschitz gradient, in which case one can use ψ(x, y) = xτ ∇ 1 E(x, y) where 0 < τ < 2/L and ∇ 1 is the gradient of E with respect to the first variable. However, the functional E we are interested in are typically not C 1 . The remaining sections of this paper develops more advanced proximal splitting schemes to define the iterations.

Iterative differentiation. The goal is then derive an algorithm to compute ∂x(y), the derivative of y ∈ Y → x(y) ∈ X . A natural way to achieve this goal is to apply an implicit function theorem to [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF]. This is equivalent to applying this implicit function theorem to the first order condition 0 ∈ ∂ 1 E(x(y), y) where ∂ 1 E is the sub-derivative of E with respect to the first variable. This is essentially the methods used in [START_REF] Kachour | The degrees of freedom of penalized l1 minimization[END_REF][START_REF] Vaiter | Robust sparse analysis regularization[END_REF] for the special case of ℓ 1 norm regularization.

While this leads to important theoretical results, the resulting closed form formula are in practice quite unstable and require the computation of x(y) with high precision. Furthermore, the explicit computation of this differential is out of reach for many variational regularization functionals E.

A more practical way to achieve this goal, introduced in [START_REF] Vonesch | Recursive risk estimation for non-linear image deconvolution with a wavelet-domain sparsity constraint[END_REF], and that we pursue here, is to compute iteratively this derivative. This is achieved by deriving formula (1), which allows, for any vector δ ∈ X to compute ξ (ℓ) = ∂x (ℓ) (y)[δ] (the derivative of y → x (ℓ) (y) applied at δ) as

ξ (ℓ+1) = Ψ (ℓ) 1 (ξ (ℓ) ) + Ψ (ℓ) 2 (δ),
where we have defined, for k = 1, 2, the following linear operators

Ψ (ℓ) k = ∂ k ψ(x (ℓ) (y), y).
The following sections are devoted to the extension of this approach to more complicated iterative algorithm that are able to handle non-smooth functionals E.

Differentiating Proximal Splitting Schemes 4.1. Proximal Operator

The proximal operator associated to a proper lower semi-continuous (lsc) and convex function

x → G(x, y) is Prox G (x, y) = argmin z 1 2 ||x -z|| 2 + G(z, y).
A function for which Prox G (x, y) can be computed in closed-form is dubbed simple. A distinctive property of Prox G (•, y) that plays a central role in the sequel is that its is a 1-Lipschitz mapping. When y is fixed, we will denote Prox G (x) instead of Prox G (x, y) to lighten the notation. The Legendre-Fenchel conjugate of G is G * (z, y) = max x x, z -G(x, y). A useful proximal calculus rule is Moreau's identity: x = Prox τ G * (x, y) + τ Prox G/τ (x/τ, y), τ > 0 .

Generalized Forward Backward Splitting

The Generalized Forward Backward (GFB) splitting [START_REF] Reguet | Generalized forward-backward splitting[END_REF] is allows one to solve a variation problem of the form

x(y) ∈ argmin x∈X E(x, y) = F (x, y) + Q i=1 G i (x, y) (3)
under the hypothesis that F is C 1 with 1/L gradient and the G i functions are simple. It reads, for all i = 1, . . . , Q, z

(ℓ+1) i = z (ℓ) i -x (ℓ) + Prox nγG i (X (ℓ) ) where X (ℓ) = 2x (ℓ) -z (ℓ) i -γ∇ 1 F (x (ℓ) )
and

x (ℓ+1) = 1 Q Q i=1 z (ℓ+1) i
. One recovers as a special case the Foward-Backward method [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] when Q = 1 and the Douglas-Rachford [START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] when F = 0.

Please note that in this algorithm, the iterates z

(ℓ)
i , x (ℓ) , x(ℓ) are actually functions of the observations y.

For any vector δ ∈ X , we wish to compute

ξ (ℓ) = ∂x (ℓ) (y)[δ], ζ (ℓ) i = ∂z (ℓ) i (y)[δ] and ξ(ℓ) = ∂ x(ℓ) (y)[δ]
. The iterations on the derivatives reads

ζ (ℓ+1) i = ζ (ℓ) i -ξ (ℓ) + G (ℓ) i,1 (Ξ (ℓ) ) + G (ℓ) i,2 (δ) where Ξ (ℓ) = 2ξ (ℓ) -ζ (ℓ) i -γ(F (ℓ) 1 (ξ (ℓ) ) + F (ℓ) 2 (δ)) and ξ (ℓ+1) = 1 Q Q i=1 ζ (ℓ+1) i
, where we have defined the following linear operators for k = 1, 2:

G (ℓ) i,k = ∂ k Prox nγG i (X (ℓ) , y) and F (ℓ) k = ∂ k ∇ 1 F (x (ℓ) , y) .

Primal-Dual Splitting

Proximal splitting schemes can be used to solve the large class of variational problems

x(y) ∈ argmin x∈X E(x, y) = H(x, y) + G(Kx, y) , (4)
where both x → H(x, y) and u → G(u, y) are proper, lsc, convex and simple functions, and K : X → U is a bounded linear operator. The primal-dual relaxed Arrow-Hurwicz algorithm as proposed in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] to solve (4) reads

u (ℓ+1) = Prox σH * (U (ℓ) ) where U (ℓ) = u (ℓ) + σK x(ℓ) , x (ℓ+1) = Prox τ G (X (ℓ) ) where X (ℓ) = x (ℓ) -τ K * u (ℓ) , x(ℓ+1) = x (ℓ+1) + θ(x (ℓ+1) -x (ℓ) ) (5)
where u (ℓ) ∈ U, x (ℓ) ∈ X and x(ℓ) ∈ X . The parameters σ > 0, γ > 0 are chosen such that σγ||K|| 2 < 1, and θ ∈ [0, 1] to ensure convergence of x (ℓ) toward a global minimizer of (4). θ = 0 corresponds to the Arrow-Hurwitz algorithm, and for θ = 1 a convergence rate of O(1/ℓ) was established on the restricted duality gap [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

For any vector δ ∈ Y, our goal is to compute the derivatives ξ

(ℓ) = ∂x (ℓ) (y)[δ], υ (ℓ) = ∂u (ℓ) (y)[δ] and ξ(ℓ) = ∂ x(ℓ) (y)[δ].
Using the chain rule, the sequence of derivatives then reads

υ (ℓ+1) = H (ℓ) 1 (Υ (ℓ) ) + H (ℓ) 2 (δ) where Υ (ℓ) = υ (ℓ) + σK ξ(ℓ) , ξ (ℓ+1) = G (ℓ) 1 (Ξ (ℓ) ) + G (ℓ) 2 (δ) where Ξ (ℓ) = ξ (ℓ) -τ K * υ (ℓ) , ξ(ℓ+1) = ξ (ℓ+1) + θ(ξ (ℓ+1) -ξ (ℓ) ) (6)
where we have defined the following linear mappings for k = 1, 2 with ∂ k the derivative w.r.t. the k-th argument

H (ℓ) k (•) = ∂ k Prox σH * (U (ℓ) , y)[•] and G (ℓ) k (•) = ∂ k Prox τ G (X (ℓ) , y)[•].

Discussion on convergence issues

One has to be aware that given that the proximal mappings are not necessarily differentiable everywhere, its differential is actually set-valued. Therefore, one should appeal to involved tools from non-smooth analysis to make the above statements rigorous. We prefer not to delve into these technicalities for the lack of space.

Another major issue is to theoretically ensure the existence of a proper sequence ξ (ℓ) that converges toward ∂x(y) [δ]. Regarding existence, as Prox G (•, y) is a 1-Lipschitz mapping of its first argument. Furthermore, in all the considered application, Prox G (x, •) is also Lipschitz with respect to its second argument. If one starts at an appropriate initialization, by induction, y → x (ℓ) (y) is also Lipschitz, hence differentiable almost everywhere. Note that for sparse synthesis regularization, convergence can be ensured as for ℓ large enough, ∂x (ℓ) (y) is constant equal to ∂x(y). As far as convergence is concerned, this remains an open question in the general case, and we believe this would necessitate intricate arguments from non-smooth and variational analysis. This is left to future research.

Risk Estimator

Projected GSURE. Recall that Π = Φ * (ΦΦ * ) + Φ is the orthogonal projector on ker(Φ) ⊥ = Im(Φ * ), where + stands for the Moore-Penrose pseudo-inverse. Let µ(y) = Πf (y) the projected estimator of Πf 0 . While f (y) is not necessarily uniquely defined, we assume that µ(y) is unambiguously defined as a single-valued mapping of the observation y. This can be ensured under a strict convexity condition on H or G in (4) (see e.g. example ( 9)).

Let µ 0 (y) = Φ * (ΦΦ * ) + y the maximum likelihood estimator. By generalizing the projected GSURE of [START_REF] Pesquet | A SURE approach for digital signal/image deconvolution problems[END_REF] to any linear operator Φ, we have GSURE(y) =||µ 0 (y)µ(y)|| 2σ 2 tr((ΦΦ * ) + ) + 2σ 2 div((ΦΦ * ) + Φf (y)) [START_REF] Eldar | Generalized SURE for exponential families: Applications to regularization[END_REF] where div(g)(y) = tr(∂g(y)) is the divergence of the mapping g : Y → Y. Under weak differentiability of y → µ(y), one can prove that the GSURE is an unbiased estimate of the risk on Im(Φ * ), i.e. E w (GSURE(y)) = E w (||Πf 0µ(y)|| 2 ) (see Appendix A).

Iterative Numerical Computation. One of the bottlenecks in calculating the GSURE(y) is to efficiently compute the divergence term. Using the Jacobian trace formula of the divergence, it can be easily seen that div((ΦΦ

* ) + Φf (y)) = E z ( ∂f (y)[z], µ 0 (z) ) ≈ 1 k k i=1 ∂f (y)[z i ], µ 0 (z i ) (8)
where z ∼ N (0, Id P ) and z i are k realizations of z. Since f (y) and its iterates f (ℓ) (y) are defined as explicit functions of x(y) and x (ℓ) (y) (see Section 6 for a detailed example), the GSURE(y) can in turn be iteratively estimated by plugging ∂x (ℓ) (y)[z i ] provided by ( 6) into (8).

Numerical Results

Total Variation Regularization

Total variation regularization of linear inverse problems amounts to solving

f (y) ∈ argmin f 1 2 ||Φf -y|| 2 + λ||∇f || 1 , (9)
where ∇f ∈ R N ×2 is a discrete gradient. The ℓ 1ℓ 2 norm of a vector field t

= (t i ) N i=1 ∈ R N ×2 , with t i ∈ R 2 , is defined as ||t|| 1 = i ||t i ||.
CP formulation. Problem ( 9) is a special instance of (4) letting x = f , H(x, y) = 0, ∀(x, y) and

K(x) = (Φx, ∇x) and ∀ u = (s, t) ∈ R P × R N ×2 , G(u, y) = 1 2 ||s -y|| 2 + λ||t|| 1 .
Separability of G in s and t entails that

Prox τ G (u, y) = ((1 -τ )s + τ y, T λτ (t)) ,
where T ρ , ρ > 0, is the component-wise ℓ 1ℓ 2 soft-thresholding, defined for i = 1, . . . , N as

T ρ (t) i = max(0, 1 -ρ/||t i ||)t i and ∂T ρ (t)[δ t ] i = 0 if ||t i || ρ δ t,i -ρ ||ti|| P ti (δ t,i ) otherwise , (10)
where P α is the orthogonal projector on α ⊥ for α ∈ R 2 , and T ρ (t) i although not differentiable on the sphere {t i : ||t i || = ρ}, is directionally differentiable there. Therefore

∂ 1 Prox τ G (u, y)[δ s , δ t ] = ((1 -τ )δ s , ∂T λτ (δ t )) and ∂ 2 Prox τ G (u, y)[δ y ] = (τ δ y , 0) .
GFB formulation. It can also be recast as ( 3)

using x = (f, u) ∈ X = R N × R N ×2 , Q = 2
simple functional, and for x = (f, u)

F (x, y) = 1 2 ||Φf -y|| 2 , G 1 (x, y) = λ||u|| 1 , and
G 2 (x, y) = ι C (x) where C = {x = (f, u) \ u = ∇f } .
One has

∇ 1 F (x, y) = (Φ * (Φf -y), 0),
and thus

∂ 1 ∇ 1 F (x, y)[δ f , δ u ] = (Φ * Φδ f , 0) and ∂ 2 ∇ 1 F (x, y)[δ y ] = (-Φ * δ y , 0)
and

Prox τ G1 (x, y) = (f, T λτ (u))
where T ρ is the vectorial soft thresholding defined in eq. ( 10), and thus

∂ 1 Prox τ G1 (x, y)[δ f , δ u ] = (δ f , ∂T λτ (δ u )) and ∂ 2 Prox τ G1 (x, y)[δ y ] = (0, 0) and Prox τ G2 (x, y) = ((Id + ∆) -1 (u + div f ), ∇(Id + δ) -1 (u + div f ))
where ∆ is the Laplace operator and div corresponds to the adjoint operator of ∇, and thus

∂ 1 Prox τ G2 (x, y)[δ f , δ u ] = ((Id + ∆) -1 (δ u + div δ f ), ∇(Id + ∆) -1 (δ u + div δ f )), and
∂ 2 Prox τ G2 (x, y)[δ y ] = (0, 0) .
Figure 1 shows an application of our GSURE to estimate the optimal λ parameter. It is applied to a super-resolution problem, where Φ ∈ R P ×N is a vertical sub-sampling matrix, where (P/N = 0.5), σ = 10 (for an image f 0 with a range [0, 255]). For each value of λ in the tested range, GSURE(y) is computed for a single realization of y using (8) with k = 1 realizations z i . Fig. 2 depicts an application of our GSURE to adjust the value of λ optimizing the recovery for a deblurring problem, where Φ ∈ R N ×N is a convolution matrix (Gaussian kernel of width 2 px), σ = 10 (for an image f 0 with a range [0, 255]). For each value of λ in the tested range, GSURE(y) is computed for a single realization of y using [START_REF] Giryes | The projected GSURE for automatic parameter tuning in iterative shrinkage methods[END_REF] with k = 4 realizations z i . 1

1 Without impacting the optimal choice of λ, the two curves have been vertically shifted for visualization purposes.

Sparse synthesis with block sparsity

Mixed ℓ 1ℓ 2 norm promoting block sparsity of linear inverse problems amounts to solving f (y) = Ψ * x(y) and x(y) ∈ argmin

x 1 2 ||ΦΨx -y|| 2 + λ Q i=1 ||B i x|| 1 , (11)
where Ψ ∈ R N ×N is an orthonormal synthesis dictionary and B i : R N → R N ′ ×B are Q linear operators extracting blocks of size B such as N = N ′ × B. The operators B i are designed such that each of them corresponds to a different partition of R N into N ′ non-overlapping blocks of size B. The

ℓ 1 -ℓ 2 norm of t = (t i ) N i=1 ∈ R N ′ ×B , with t i ∈ R B , is defined as ||t|| 1 = i ||t i ||.
CP formulation. When Q = 1, problem ( 9) is a special instance of (4) letting x = Φ * f , H(x, y) = 0, ∀(x, y) and

K(x) = (Φx, B 1 x) and ∀ u = (s, t) ∈ R P × R N ×2 , G(u, y) = 1 2 ||s -y|| 2 + λ||t|| 1 .
The proximal operator associated to G has already been given in the total variation case. When Q > 1, CP formulation cannot be used and one could instead use the GFB formulation.

GFB formulation. It can also be recast as (3) using x = Φ * f , Q simple functionals, and for all i = 1, . . . , Q

F (x, y) = 1 2 ||ΦΨx -y|| 2 , G i (x, y) = λ||B i x|| 1 .
One has

∇ 1 F (x, y) = Ψ * Φ * (ΦΨx -y),
∂ 1 ∇ 1 F (x, y)[δ x ] = Ψ * Φ * ΦΨδ x and ∂ 2 ∇ 1 F (x, y)[δ y ] = -Ψ * Φ * δ y and Prox τ Gi (x, y) = B * i T λτ (B i x)
where T ρ is the vectorial soft thresholding defined in eq. ( 10), and thus Fig. 3 depicts an application of our GSURE to adjust the size of blocks B optimizing the recovery for a compressed sensing problem, where Φ ∈ R P ×N is a realization of a random matrix distribution (randomized sub-sampling of a random convolution), where (P/N = 0.5), σ = 10 (for an image f 0 with a range [0, 255]). For each size B, Q = B 2 , i.e., all possible partitions have been used, and λ has been set to the value 0.7σ/B. For each value of λ in the tested range, GSURE(y) is computed for a single realization of y using (8) with k = 1 realizations z i .

Conclusion

We obtained proximal splitting derivatives for unbiasedly estimate the projected risk in regularized inverse problems handling both synthesis and analysis sparsity priors as well as mixed norms for block structured sparsity. Its usefulness has been illustrated on automatic choice of the regularization parameter for total variation regularization and automatic choice of the size of blocks for sparse synthesis with block sparsity.

Appendix A. Proof of equation [START_REF] Eldar | Generalized SURE for exponential families: Applications to regularization[END_REF] Assume w → g(w) is weakly differentiable in the sense of [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF], Stein's lemma reads E w w, g(w) = σ 2 div g(w) .

Under weak differentiability of y → µ(y) and using the fact that µ 0 (y) = Πf 0 + Φ * (ΦΦ * ) + w and µ(y) = Πf (y), one has

Figure 1 .

 1 Figure 1. (a) Projected risk and its GSURE estimate 1 (b) y. (c) f (y) at the optimal λ.

Figure 2 .

 2 Figure 2. (a) Projected risk and its GSURE estimate 1 (b) y. (c) f (y) at the optimal λ.6.2. Sparse synthesis with block sparsityMixed ℓ 1ℓ 2 norm promoting block sparsity of linear inverse problems amounts to solving

Figure 3 .

 3 Figure 3. (a) Projected risk and its GSURE estimate 1 (b) µ 0 (y). (c) f (y) at the optimal λ. and thus

∂ 1

 1 Prox τ Gi (x, y)[δ x ] = B * i ∂T λτ (B i δ x ) and ∂ 2 Prox τ Gi (x, y)[δ y ] = 0 .

E||µ 0 2 = 2 = 2 = 2 =

 02222 (y)µ(y)|| E||µ 0 (y)|| 2 -2E µ 0 (y), µ(y) + E||µ(y)|| E||Πf 0 + Φ * (ΦΦ * ) + w|| 2 -2E Πf 0 + Φ * (ΦΦ * ) + w, µ(y) + E||µ(y)|| E||Πf 0 || 2 + σ 2 tr((ΦΦ * ) + ) -2E Πf 0 , µ(y) -2E w, (ΦΦ * ) + Φf (y) + E||µ(y)|| E||Πf 0µ(y)|| 2 + σ 2 tr((ΦΦ * ) + ) -2σ 2 div((ΦΦ * ) + Φf (y)) .