H Arbess

M Bafleur

D Trémouilles

M Zerarka

Combined MOS-IGBT-SCR structure for a compact high-robustness ESD power clamp in smart power SOI technology

Smart power technologies are required to withstand high ESD robustness, both under powered and unpowered conditions, in particular for automotive and aeronautics applications among many others. They are concurrently confronted to the challenges of high-temperature operation in order to reduce heat sink related costs. In this context, very compact high-robustness ESD protections with low sensitivity to temperature are required. To fulfill this need, we studied a new ESD protection structure that combines in the same component MOS, IGBT and thyristor effects. This is achieved by inserting in the same LDMOS device P + diffusions in the drain. We studied the impact of N + /P + ratios on R ON and holding current at high temperatures. Structures optimization have been realized with 3D TCAD simulation and experimentally validated. The proposed structures provide a high ESD robustness with small footprint and reduced temperature sensitivity compared to classical solutions. Original design solutions to improve their immunity to latch-up are also presented.

I. INTRODUCTION

Within the context of a sustainable development and to face the challenges of fossil energy shortage, transports are becoming more electrically driven in order to optimize overall energy consumption efficiency. This requires new generations of power devices and electronic circuits. With the advent of wide band gap semiconductors such as GaN and SiC, the driving circuitry that is still realized using silicon technologies, has to be placed as close as possible to the power devices and should be able to operate at high temperature (≥ 200°C).

To overcome isolation issues induced by high temperatures, implementing a full dielectric isolation (Silicon-On-Insulator, SOI) in silicon Smart Power technologies is becoming more and more interesting. It allows managing very high blocking voltages, reducing parasitic bipolar effects and increasing circuits speed.

In addition to high temperature constraint, ESD requirements from the integrated circuit to the system are particularly severe: for example, in automotive applications, 2-8kV in contact according to the IEC 61000-4-2 standard [START_REF]Electromagnetic compatibility (EMC) -Part 4-2: Testing and measurement techniques -Electrostatic discharge immunity test[END_REF] is often required. It has to be reminded that this ESD standard requires tests both for unpowered and powered systems. For the latter testing, the effect of temperature has then to be taken into account in the ESD design window. Protection should provide a high failure current but also work at high temperature, its performance being as insensitive as possible to this parameter.

In this paper, we propose a study of an ESD protection solution based on combining MOS and IGBT structures in the same device and allowing SCR triggering to reach a high ESD robustness and a low sensitivity to temperature.

In section II the proposed device and its technology are described. A first analysis based on simulation is carried out in section III. Finally, in section IV, the experimental results are thoroughly analyzed including triggering overshoot and latch up risk.

II. SMART POWER SOI TECHNOLOGY AND PROPOSED ESD PROTECTION

The chosen technology is a smart power SOI (TFSMART 1) provided by Telefunken Semiconductors [START_REF]Release: SIS1_0.32_071005[END_REF]. It is a 0.8 µm Bipolar CMOS DMOS merged technology on SOI. The buried oxide and active silicon thicknesses are 500 nm and 2µm, respectively. The process offers 0.8µm five-volts CMOS devices, lateral npn-and pnp-transistors, P-type and N-type lateral DMOS devices with complete deep trench isolation between wells and vertical oxide isolation to handle wafer/substrate for different supply voltages

Regarding the N-LDMOS devices, TFSMART 1 library contains two types of structures. The first one is named Source Body Short (SBC), since the source N + diffusion is shorted to the P-type body diffusion. This short is implemented by inserting within the source N + diffusion a P + body contact. The second one, called Source Body Open (SBO), has a Shallow Trench Isolation (STI) between the source and the P + body contact (Fig. 1). This latter structure allows biasing the body at a different voltage from the source.

The TFSMART 1 library proposes a high-voltage active MOS-based ESD-power-clamp. This is actually a 25V N-LDMOS with increased gate-drain ballasting distance. It is composed of 11 identical cells that are formed of two 150 µmfingers with a central drain diffusion, the body being connected to the source. The corresponding surface is 302 µm x 183 µm and it provides a 2 kV HBM robustness according to TFSMART1 foundry data [START_REF]Release: SIS1_0.32_071005[END_REF] that correlates with our measurement. We observed a robustness equal to 1.74 A at 25°C and 1.34 A at 200°C for 100 ns current pulses. This is a dramatic robustness reduction of 25% with temperature. The performance of this MOS-based protection structure, namely its on-resistance, is very sensitive to the temperature [START_REF] Maloney | Novel clamp circuits for IC power supply protection[END_REF]. According to our measurements, the on-resistance increases from 4.4 Ω at 25°C to 7.14 Ω at 200°C. If this temperature behavior is not taken into account during a circuit's ESDprotection-design, it could induce detrimental effects such as a lower failure current and even could fail providing the expected protection due to the lack of compliance with the ESD design window. To compensate this effect, the size of the power clamp protection would have to be increased. To cope with these issues, we contemplated developing a new ESD protection. The basic concept presented in this paper consists in combining MOS and bipolar effects in the same structure in order to compensate the detrimental effects of the temperature. To improve the ESD robustness, we also take advantage of the turn-on of the parasitic thyristor of the structure that will provide a very low on-resistance.

We implemented within the structure of the N-LDMOS a lateral IGBT by inserting in the drain region both N + and P + diffusions (Fig. 2).

The equivalent electrical schematic is presented in Fig. 3. The protection is triggered via its gate to allow protecting low voltage circuitry and a nearly triggering of the structure (MOS then IGBT) to compensate the intrinsic SCR slow response.

Such MOS-thyristor combination was already proposed in our laboratory (LAAS) in 1985 by H. Tranduc to improve the performance of power device structures [START_REF] Tranduc | Le transistor-thyristor métal-oxyde-semiconducteur (T2 MOS)[END_REF]. More recently a pure LIGBT was proposed to implement an efficient ESD power clamp in an SOI smart power technology [START_REF] Gevinti | Novel 190V LIGBTbased ESD protection for 0.35µm Smart Power technology realized on SOI substrate[END_REF]. The authors demonstrated the improved performance compared to an NDRIFTMOS device. However, combining the three devices'action (MOS, IGBT, SCR) in a single device was not yet studied to the best of our knowledge.

Fig. 2 :Structure of the "elementary cell" of the proposed LDMOS-LIGBT.Here 50%-IGBT, i.e. the P + /N + ratio at the drain side is equal to 1. Two silicon runs (RUN01 and RUN02) have been made in order to optimize the behavior of these new devices.

We use a naming convention where the N-LDMOS is considered as the reference structure and the IGBT percentage represents the proportion of P + diffusion among the full device drain width. For example, a 20%-IGBT device has 20% of P + and 80% of N + in its drain. Thus, 0%-IGBT is the reference N-LDMOS (no IGBT action) and 100% -IGBT is the full LIGBT (no NMOS action). An illustration of a 50%-IGBT is shown in Fig. 2 where the N + /P + diffusions ratio is thus equal to 1.

In a previous work presented in [START_REF] Arbess | High-temperature operation MOS-IGBT power clamp for improved ESD protection in smart power SOI technology[END_REF], we designed a test chip with various structures going from a N-LDMOS to a LIGBT and mixed structures with an increased IGBT percentage of 33, 50 and 66%. The proposed structures provide a high ESD robustness (> 5.5 A) with a ten times smaller silicon area compared to the LDMOS-based ESD power clamp from the library. Another beneficial effect of this design is that the performance of the ESD protection becomes much less sensitive to temperature.

The main problem we encountered in this initial study was the small value of the SCR holding current. In the next section, we propose solutions to overcome this issue, firstly optimized using 3D TCAD simulations and then experimentally validated on silicon.

III. INVESTIGATION ON HOLDING CURRENT OPTIMIZATION BASED SOURCE SIDE TCAD ANALYSIS

The holding current of the parasitic SCR is defined by the following parameters and equation [START_REF] Troutman | Latchup in CMOS technology: the problem and its cure[END_REF]:

! = ! ! ! ! !! .! !" !! ! ! ! !! .! !" ! ! ! ! !! (1)
Where β n and β p are the current gains of the parasitic NPN and PNP bipolar transistors and I NW and I PW the currents flowing into the N-well and P-well, respectively.

To increase I H , as we can hardly modify the current gain of the two bipolar transistors, the classical way is to reduce as much as possible the resistances of the two wells [START_REF] Mergens | High holding current SCRs (HHI-SCR) for ESD protection and latch-up immune IC operation[END_REF]. In the proposed structure, this means on the one hand optimizing the N + /P + ratio in the drain and, on the second hand, introducing P + diffusions into the N + source diffusion to reduce the P-well resistance;

To find out the best trade-off between high I H and low onresistance, we performed 3D simulations. They show that varying the N + /P + ratio in both source and drain regions allows controlling and increasing I H whereas preserving the ESD performance. However measurements showed that I H is almost independent of the N + /P + ratio at the drain side, which means that the SCR behavioris mainly controlled by the parasitic NPN bipolar transistor (N+ source, P-well, N-well). Therefore, the most efficient way to increase I H consists in engineering the source side of the device that defines this parasitic bipolar transistor.

All the structures used in the first test run were of SBO type (Fig. 1) with external short circuit between source and body contact since the original LDMOS was implemented that way.

A first solution to optimize the I H value is to use an SBC-type structure that allows significantly reducing the P-Well resistance (base of the NPN bipolar transistor). Using Sentaurus TCAD, we simulated a 3D structure to check the impact on the holding current. The structure under study is based on an NLDMOS 25 V. We could not have access to accurate doping profiles for this technology. To perform the 3D simulation on the NLDMOS 25V, we performed an extrapolation from another SOI technology and did some fitting according to TLP measurement. As a consequence, the simulation results can only be considered as qualitative to compare the relative efficiency of different design solutions.

To perform the three-dimensional (3D) simulations, which are known to be time consuming, we used a quasi-stationary simulation and an elementary 3D cell of 3µm-width to minimize computation time. In addition, since the simulations are carried out in the low current regime, around the holding current point, self-heating does not have a very large impact on the device behaviour; as a consequence the heat equations were not included in the simulations.

The computed I H should be scaled according to the real structure width from the results obtained with the 3 µm-width elementary cell. According to equation ( 1) I H is inversely proportional to the well resistances via I PW and I NW currents, and a scaling coefficient can be computed. According to the way we compute it, its value can be in the range of 6 to 10.

Fig. 4 shows that I H increases from 5 mA for a SBO structure to 8 mA for a SBC one for simulated 3µm-wide 50%-IGBT devices at room temperature and V GS = 3V.

As the increase of I H from 5 to 8 mA is not sufficient, another interesting solution consists in locally inserting a P + diffusion into the source either as close as possible to the gate (Optimized mixed structure SBC Fig. 5 middle), or even totally replacing the N + source diffusion (Fig. 5 right). In that latter case, the MOS channel will not be able to form then resulting in a structure with a reduced channel and the local elimination of the parasitic NPN bipolar transistor. We named this latter structure "SBC structure with channel reduction".

Fig. 4: TLP characteristics difference between SBO mixed structure, SBC mixed structure, optimized mixed structure that illustrated in Fig. 5 and mixed structure with reduced channel by 30% at room temperature and at V GS = 0V. The inset in the right shows the zoom of the three first curves. should be in the range of 100 mA therefore meeting the latch up specification. The experimental result will be detailed and discussed in section VI-b.

IV. OPTIMIZED STRUCTURES ELECTRICAL CHARACTERIZATION RESULTS

Starting from the simulation results we have designed a number of optimized mixed structures. We designed structures with a N + /P + ratio in the drain varying between 0 and 4 (0, 1/4, 1/3, 1/2, 1, 2, 3 and 4) which corresponds respectively to 100, 80, 75, 66.6, 50, 33.3, 25 and 20%-IGBT device according to the naming convention defined in section II.

For each N + /P + ratio, the IGBT P + diffusion can be implemented as a single segment or partitioned into several ones. For example, we designed two configurations for the 80%-IGBT structure, and four configurations for the 75%-IGBT. However, this does not lead to significant changes of the devices behavior and thus will not be detailed here.

For the design of the mixed structures, we started from 25V ESD NLDMOS of SBO type from TFSMART1 library. To reduce the PWell resistance, we also designed different configurations for the source side of the device. We introduced P + islands of 3.2 µm abutting the N+ source and locally replacing the STI.

Both drain side and source side configuration will be detailed later on in the next sections.

To characterize the optimized structures of the second test chip called RUN02, we used the same techniques as for the first test chip (RUN01): an analogic curve tracer to extract holding current and a100 ns pulse generator (TLP bench) to characterize the ESD structures robustness, trigger voltage and serial resistance, both at low and high temperatures.

a. Robustness results and temperature effects

Measurements are done until failure defined when the leakage current becomes significant (>1µA) at room temperature. For a given N + /P + ratio, we did not notice any significant impact of drain or source configuration on the robustness (Table 1). We noted that the robustness of the structures increases with the percentage of P + in the drain with one exception, which is the LIGBT. Indeed the presence of P + in the drain allows the triggering of the parasitic thyristor associated with the IGBT, and then a current distribution within the volume of silicon. Thus in a mixed structure, most of the current flows through the thyristor region of the component (lower resistance), so if this region is larger, the structure can hold more power and the opposite is true. Between the 80%-IGBT structure and the LIGBT one, the robustness decreases from 10 A to 7.5 A. A linear extrapolation of the robustness with the percentage of P + into the drain would result in 13 A for the LIGBT robustness.

Fig. 6(a) presents the variation of the normalized failure current and resistance as a function of P + drain percentage. In this Figure, we see that although its on-state resistance is smaller than all mixed structures, the robustness of the LIGBT is lower. From Fig. 6(b), we can see that the LIGBT power failure is almost equal to that of the 25%-IGBT structure. Note that the power to failure is calculated by multiplying the failure current by the failure voltage.

This exception case (LIGBT) will be analyzed in detail in the next sub-section. Moreover, measurements up to failure show that the robustness does not change by changing the gate bias. It is the same for the on-state resistance. Indeed, increasing the voltage V GS , only changes the low-current part of the curve (linear and saturated) but right after the triggering of the thyristor, the gate bias has no more impact. Fig. 7 shows the TLP characteristics as a function of temperature for grounded gate 20%-IGBT (1P4N), 80%-IGBT (4P1N) and LIGBT structures. The results show a decreased robustness when increasing the temperature for all components. This decrease is accompanied by an increase of the on-resistance which in turn results in increased self- IGBT heating and then to early failure. Fig. 7: TLP characteristics of 1P4N structure (20% IGBT) (up) and 4P1N structure (80% IGBT) (center) and LIGBT (bottom) with grounded gate and as a function of temperature. We can remark the three part of structures name in the legend: the first is the drain configuration, the second is the segments number of the mentioned diffusion in the drain and the third number in the second part is the source configuration.

b. Investigation on the LIGBT unexpected behavior

Fig. 8 shows a photoemission view of the first configuration of 1P2N (33% IGBT) structure where the IGBT is divided into three segments. This Figure shows a photo-emission captured at a current level of 4 A. It clearly shows that the power is distributed over the different regions of the IGBT structure (3 light-emitting areas) that correspond to the three P + zones of this structure. At this current level, we observe a small current crowding phenomenon. It may be noted also that in this structure, which has two fingers, a larger current flows in the upper finger (light emitting more intense) than in the lower one. The upper finger is the first to be triggered due to a lower source access resistance. This difference between upper and lower fingers make the robustness smaller than be. A first explanation for this low experimental I t2 value of the LIGBT could be the occurrence of a current focalization in the structure that does not occur in the case of mixed structures. However, an observation in photoemission of the LIGBT at 6A shows that the current is homogeneous, even at the periphery of the drain.

We analyzed the TLP voltage and current waveforms across the LIGBT and mixed structure. It shows that in the LIGBT, the SCR needs more time to propagate throughout the structure and have therefore the lowest snapback voltage. Indeed, Fig. 9 presents the TLP characteristic of a LIGBT for V GS = 5V on the left figure and current waveforms (bottom right) and voltage (upper right) with a value of I TLP = 5A. Fig. 10 shows the same curves, for the 80%-IGBT structure with one P + segment in the drain. In the first case, the snapback voltage firstly reaches 75 V at t = 25 ns, and then continues to decrease linearly until t = 50 ns. In the second case, the snapback voltage is steep. So we can conclude that this delay in the LIGBT generates a phenomenon of self-heating that is more important than in a mixed structure. The results in the previous section showed a high voltage across the DUT. To evaluate the actual overvoltage seen by the DUT before the triggering of the protection structure, we carried out transient TLP measurement using a Very Fast TLP system with a pulse duration of 2.5 ns and a rise time equal to 100 ps. The results show the decrease of the overvoltage with increasing gate bias (Fig. 11). This decrease becomes more significant with a higher percentage of IGBT structure.

The high value of the overvoltage can be a cause of failure of the component or circuit to be protected if greater than the failure voltage of the component. Several methods can be considered to solve this problem and provide a perspective of this work [START_REF] Bourgeat | Local ESD protection structure based on silicon controlled rectifier achieving very low overshoot voltage[END_REF] .

The variation of overvoltage with P + (IGBT) percentage in the drain is presented in Fig. 12 at V GS = 4V and at room temperature. We can remark that the overvoltage decreases with the increase of P + percentage. This overvoltage decrease maybe related to the reduction of the on-resistance of the structure with the P + percentage.

d. SCR holding current and latchup risk

To improve latch up immunity, we focused our efforts on the source engineering by implementing P + islands to reduce the Pwell body resistance (Fig. 13). These islands are implemented using two configurations. In the first one (config 1), their width constitutes 1/8 of the total width of the structures and in the second (Config 2) 1/4. A third configuration consists in removing the STI between the P + body contact and the N + source diffusion that corresponds to a Source Body Shorted (SBC) type structure. Table 2 shows the variation of holding current for three different source configurations for a voltage V GS of 3 V. The differences in I H values between Config 1, Config 2 and SBC are not significant enough to draw a conclusion about the impact of the source configurations.

P+ percentage in the drain

Then we analyzed the impact of the drain configuration, i.e., the P + percentage variation in the drain (Fig. 14). In this figure, there are several points for a single percentage of P + , these points correspond to different segmentations of the P + for the same percentage of this latter in the drain. Although there is a slight increase of I H for structures having a P + percentage higher than the one in the 2P1N one (66% IGBT), the higher holding current is obtained with the lowest percentage of P + in the drain. The holding current exhibits a minimum value for 2P1N (67% IGBT). This can be explained by the fact that the P + percentage in the drain has two opposite effects. Increasing the P + diffusion ratio into the drain increases the quantity of holes injected in the low-doped region thus decreasing in the Nwell resistance by conductivity modulation. As a consequence a higher current is required in the NWell for the PNP bipolar transistor to turn-on. Concurrently, the increase of P + in the drain results in a greater hole current in the Pwell (source side), favoring an easier NPN bipolar transistor turn-on and then the thyristor triggering. According to the realized measurements (Fig. 14), 2P1N structure configuration corresponds to the threshold where there is no predominance of either of these effects (Nwell-resistance modulation or NPN base current increase).

e. Impact of channel reduction

For the study of the structure with reduced channel, we used non-ballasted structures that were originally designed as improved power devices. As a result, the ESD robustness is almost halved by the absence of the ballast resistor: for example, for the 4P1N structure it decreases from 10 A with the ballast resistor down to 6.5A without it.

Fig. 15 shows these specific configurations (source side) for 1P2N structure without ballast drain resistance to study the effect of that area in this technology. They are respectively from top to bottom: 1P2N-20PC, 33%-IGBT structure with20% channel reduction, 1P2N-10PC,33%-IGBT structure with 10% reduction channel reduction and 1P1N-1P1N, 50%-IGBT structure with SBC source configuration.

To implement this channel reduction, we replaced 10% (20%) of the N + diffusion of the source by P + diffusion. To study the channel reduction effect on the holding current value, we made three different measurements for each structure using the curve tracer to get more accurate values.

We studied the effect of the channel reduction on the robustness and the on-resistance of the structures. The TLP measurement showed that there is no change in resistance or robustness.

Indeed, the simulation showed that the holding voltage for reduced channel structures is much larger than the others. Therefore, if we assume that the maximum power failure is constant for these structures, their robustness should be reduced by a percentage proportional to the increase of the holding voltage. Fig. 16 shows the I-V characteristics for 1P2N (percentage of P + into the drain = 33.3%) and for LIGBT measured at room temperature and with a voltage V GS of 5 V. It can be noted from this figure that the difference between a normal structure and a structure with a 10% channel reduction is hardly noticeable. However, the 20% channel reduction has a significant effect. This was observed for all the structures. In Fig. 16, for the configuration 20PC, and after thyristor triggering, it can be noticed that the characteristic of the LIGBT snaps back several times. This behavior may be due to the different phases of current distribution in the structure. Each snap back shows that there is a new current filament that opens. In addition to the increase in holding current, reducing the channel also allows increasing the holding voltage.

V. CONCLUSION

We proposed a new mixed ESD component MOS-IGBT improved to operate at high temperature and have a high robustness. The basic concept is to trigger the structure as an IGBT but to also allow its parasitic SCR triggering. To do so, we combined within a same N-LDMOS a lateral IGBT by implementing P + diffusions in the N + drain region with different N + /P + ratios. The proposed structures provide a high ESD robustness (> 5.5 A) with a ten times smaller silicon area compared to the initial LDMOS-based power clamp. The best measured robustness is equal to 10 A or 1,7x10 5 A/cm². Another beneficial effect of this design is that the performance of the ESD protection becomes much less sensitive to temperature. We proposed original design solutions to improve the immunity to latch-up of these structures by engineering the source side, on the one hand, and providing an additional control viathe gate, on the other hand. These solutions were validated through comparative 3D simulations and ona silicon test vehicle. The most efficient solution to control the SCR holding current is the channel reduction with IGBT percentage equal to 33 %. The highest measured robustness is equal to 10 A or 1.7x10 5 A/cm² with IGBT percentage equal to 80 %. The proposed optimized elementary cell can already be used either in series or in parallel to provide the required latch-up immunity.

Fig. 1 :

 1 Fig. 1: SBC and SBO types structures.

Fig. 3 :

 3 Fig. 3: Equivalent electrical schematic of the mixed structure LDMOS-LIGBT.

Fig. 5 :

 5 Fig. 5: Source side optimization of an50%-IGBT:SBC mixed structure (a), Optimized SBC mixed structure (b) and SBC mixed structure with reduced channel (c).The result of the simulation for the optimized mixed structure (Fig.5middle) shows that I H increases from 8 mA for a SBC structure to 12 mA. For the reduced channel structure I H raises up to 40 mA. However, in this latter case, the holding voltage increases to 7.5 V to be compared to 4 V for the three first ones. With such holding voltage increase there is a risk to lower the ESD robustness of the structure. Finally, taking into account the simulation-scaling coefficient, the resulting I H

Fig. 6 :

 6 Fig. 6: Failure current and on-resistance variation as a function of percentage of P + diffusion in the drain (a) and, failure power variation as a function of the percentage of P + in the drain (b) at room temperature and grounded gate

Fig. 8 :

 8 Fig. 8: Photoemission of the 1P2N structure with 4 A of current and at room temperature.

Fig. 9 :-Fig. 10 :

 910 Fig. 9: TLP Characteristics of the LIGBT structure (left) and voltage and current waveforms (right) for ITLP = 5A and VGS = 5V

Fig. 11 :

 11 Fig. 11: Waveform of the voltage across the structure 3P1N (75% IGBT) for many gate polarizations at constant temperature and initial signal amplitude 100V.

Fig. 13 :

 13 Fig. 13: Various source configurations for a structure with ballast drain resistance.

  Fig. 14: Holding current variation with P+ percentage in the drain at room temperature, for several values of VGS and for the three configurations of source.

Fig. 15 :

 15 Fig. 15: Diagram of the different configurations (source side) of the 1P2N structure without ballast resistor

Fig. 16 :

 16 Fig. 16: I-V characteristics for the 1P2N structure (up) and IGBT (bottom) at room temperature, with VGS = 5V and for several source configuration.

Failure current (A) Drain config. 1 Drain config. 2 Drain config. 3 Drain config. 4

	Source Config 1	9.33	9.73	9.54	9.93
	Source Config 2	9.52	9.73	9.14	9.38
	Source Config 3	N.A.	N.A.	9.33	N.A.

Table 1 :

 1 Failure current for the different source and drain configurations of the structure 3P1N (75% IGBT). The third source configuration is implemented with only one drain configuration.

Table 2 :

 2 Holding current variation for the different configurations

			Source configuration
		Config 1	Config 2	SBC
	Drain configuration		Holding current (mA)
	1P4N (20% IGBT)	28	27,5	27,5
	1P3N (25% IGBT)	29	28,5	28
	1P2N (33% IGBT)	27,5	27	26,5
	1P1N (50% IGBT)	24	24	22
	2P1N (66% IGBT)	22	23	22
	3P1N (75% IGBT)	22,5	23,5	21,5
	4P1N (80% IGBT)	23	23	22
	LIGBT	26	26	24

ACKNOWLEDGEMENTS

This work has been sponsored by Fondation de Recherche pour l'Aéronautique et l'Espace (http://www.fnrae.org/) within the Framework of the collaborative project COTECH.


Pietro Belotti
email: pbelott@clemson.edu

Sonia Cafieri
email: sonia.cafieri@enac.fr

Jon Lee
email: jonlee@us.ibm.com

Leo Liberti
email: liberti@lix.polytechnique.fr

Andrew Miller
email: andrew.miller@math.u-bordeaux1.fr

On the composition of convex envelopes for quadrilinear terms *

 we prove that having fewer groupings of longer terms yields tighter convex relaxations. In this paper we give an alternative proof of the same fact and perform a computational study to assess the impact of the tightened convex relaxation in a spatial Branch-and-Bound setting.

Introduction

One of the most crucial steps of the spatial Branch-and-Bound algorithm for solving Mixed-Integer Nonlinear Programming (MINLP) problems is the lower bound computation. When the MINLP is factorable, it is possible to construct a convex relaxation automatically by means of a particular type of lifting reformulation (called MINLP standard form [START_REF] Smith | A symbolic reformulation/spatial branch-and-bound algorithm for the global optimisation of nonconvex MINLPs[END_REF][START_REF] Liberti | Reformulations in mathematical programming: a computational approach[END_REF]) first proposed in [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF] and then exploited in most existing sBB algorithms [START_REF] Ryoo | A branch-and-reduce approach to global optimization[END_REF][START_REF] Adjiman | A global optimization method, αBB, for general twice-differentiable constrained NLPs: I. Theoretical advances[END_REF][START_REF] Smith | A symbolic reformulation/spatial branch-and-bound algorithm for the global optimisation of nonconvex MINLPs[END_REF][START_REF] Liberti | Writing global optimization software[END_REF][START_REF] Tawarmalani | Global optimization of mixed integer nonlinear programs: A theoretical and computational study[END_REF][START_REF] Belotti | Branching and bounds tightening techniques for non-convex MINLP[END_REF]. If we consider polynomial problems, higher order monomials are recursively rewritten as products of monomials of sufficiently low order for which a tight convex relaxation (possibly the convex envelope) is known. Each lower order monomial is replaced by an added variable, and an equality constraint defining the added variable in terms of the monomial it replaces is adjoined to the MINLP. This operation is carried out recursively until the MINLP consists of a linear objective, some linear constraints, and several defining constraints of the form w j = h j (x, w) for all j in some appropriate set J, where the functions h j represent monomials. To obtain a convex relaxation, each defining constraint is replaced by a set of constraints defining the convex relaxation of its feasible set, thus yielding a convex relaxation for the whole problem.

Let B = [x L , x U ]. The quadrilinear feasible set S 4 = {(w 1 , x 1 , x 2 , x 3 , x 4 ) | w 1 = x 1 x 2 x 3 x 4 } ∩ B over a box can be lifted in many different ways according to the way associativity is applied: the grouping ((x 1 x 2 )x 3 )x 4 , for example, yields the set S 2,2,2 = {(w 1 , w 2 , w 3 , x 1 , x 2 , x 3 , x 4 ) | w 2 = x 1 x 2 ∧ w 3 = w 2 x 3 ∧ w 1 = w 3 x 4 } ∩ B, whereas the grouping (x 1 , x 2 , x 3 )x 4 yields S 3,2 = {(w 1 , w 2 , x 1 , x 2 , x 3 , x 4 ) | w 2 = x 1 x 2 x 3 ∧ w 1 = w 2 x 4 } ∩ B. Since convex/concave envelopes exist in explicit form for both bilinear [START_REF] Al-Khayyal | Jointly constrained biconvex programming[END_REF][START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF] and trilinear terms [START_REF] Meyer | Trilinear monomials with mixed sign domains: Facets of the convex and concave envelopes[END_REF][START_REF] Meyer | Trilinear monomials with positive or negative domains: Facets of the convex and concave envelopes[END_REF], we can derive two different convex relaxations of S 4 . The first, S2,2,2 , consists in replacing the bilinear constraints w i = x j x k appearing in S 2,2,2 by the corresponding bilinear envelopes. The second, S3,2 , consists in replacing the trilinear terms with the trilinear envelope and the bilinear term with the bilinear envelope. A question then arises naturally: which one is tighter?

In [START_REF] Cafieri | On convex relaxations of quadrilinear terms[END_REF] we proved that S3,2 ⊆ S2,2,2 and performed a computational study of the containment of the convex relaxations when different parameters were varied. In this paper we provide an alternative proof (based on formal grammars) of the same result, and then test the impact of the tightened convex relaxation S3,2 using sBB.

The rest of this paper is organized as follows. In Sect. 2 we present the main motivations of this work and a literature review on convex relaxations for multilinear monomials and their impact on a sBB algorithm. In Sect. 3 we propose a theoretical framework, based on concepts from the formal languages theory, to compare convex relaxations of multilinear monomials obtained as a composition of convex envelopes of lower-degree monomials. In Sect. 4 we discuss some computational experiments aimed at comparing different convex relaxations of quadrilinear terms in a spatial Branch-and-Bound setting. Concluding remarks are given in Sect. 5.

Motivation and literature

The above discussion implies that deriving convex relaxations that are as strong as possible (i.e., that approximate the convex hull as closely as possible) for multilinear monomials can be critically important for the performance of a spatial branchand-bound algorithm designed to globally solve nonconvex polynomial optimization problems. Because of this, numerous efforts have studied the convex hulls of sets defined by lower order product terms and the use of these convex hulls in recursively factorized formulations (such as the MINLP standard form defined above).

Four valid inequalities for the three-dimensional set

S 2 = {w, x 1 , x 2 : w = x 1 x 2 , x ∈ [x L , x U ]
} were proposed by [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF], and later [START_REF] Al-Khayyal | Jointly constrained biconvex programming[END_REF] showed that these four inequalities suffice to describe the convex hull. At present most global MINLP solvers that use general sBB methods (among recent examples see [START_REF] Belotti | Branching and bounds tightening techniques for non-convex MINLP[END_REF], [START_REF] Lin | The global solver in the LINDO API[END_REF], [START_REF] Sahinidis | Baron 8.1.1: Global optimization of mixed-integer nonlinear programs[END_REF]) use the convex hull for recursively defined instances of S 2 to define the polyhedral relaxations that are solved at each node of the branch-and-bound tree.

However, it may be thought that limiting solvers to the use of envelopes defined by simple bilinear terms may result in convex approximations for the original problem that are less strong (perhaps much less so) than those that exploit envelopes for more complex expressions. For problems involving multilinear multinomials defined by products of more than two variables, this consideration has motivated research into the envelopes of trilinear functions [START_REF] Meyer | Trilinear monomials with mixed sign domains: Facets of the convex and concave envelopes[END_REF][START_REF] Meyer | Trilinear monomials with positive or negative domains: Facets of the convex and concave envelopes[END_REF]. Comparing the use of convex envelopes for bilinear and trilinear forms in building convex approximations for MINLPs motivated the study in [START_REF] Cafieri | On convex relaxations of quadrilinear terms[END_REF], and comparisons involving more general functional forms motivate the present article.

A natural generalization of bi-and tri-linear functions are functions that are known to have vertex polyhedral convex envelopes. (The convex envelope of an n-dimensional function f (x) is said to be vertex polyhedral if its domain X is a polyhedron, and if every extreme point of the convex hull of {(x, f (x) : x ∈ X} is defined by an extreme point of X itself.) In [START_REF] Meyer | Convex envelopes for edge-concave functions[END_REF] Meyer and Floudas generalized the approach developed for trilinear functions to functions with vertex polyhedral convex envelopes. Essentially, their approaches can be thought of as enumerative methods that consider all possible combinations of n + 1 extreme points of X (equivalently, extreme points of conv({(x, f (x)) : x ∈ X})), and then establish conditions under which the hyperplane defined by such a set of points defines a linear inequality satisfied by all the other extreme points of conv({(x, f (x)) : x ∈ X}). Such an inequality is then valid for {(x, f (x)) : x ∈ X} and facet-defining for the convex hull of this set.

General multilinear functions (i.e., any function composed of a sum of products of variables, in which the degree of each variable in each product is 0 or 1) were shown to have vertex polyhedral convex envelopes by [START_REF] Rikun | A convex envelope formula for multilinear functions[END_REF]. An implication of this result is that many of the concepts mentioned in the preceding paragraph can be used for general mutilinear functions; their use is not limited to monomial products (for example). The extension of such results to define convex envelopes for multilinear functions (and generalizations of them) has been discussed in [START_REF] Sherali | Convex envelopes of multilinear functions over a unit hypercube and over special discrete sets[END_REF][START_REF] Tardella | Existence and sum decomposition of vertex polyhedral convex envelopes[END_REF][START_REF] Tardella | On the existence of polyhedral convex envelopes[END_REF][START_REF] Tawarmalani | Convex extensions and convex envelopes of l.s.c. functions[END_REF], among other references.

Empirical testing of the approaches mentioned above (beyond the use of bilinear envelopes defined by [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF]) has been limited, but recently authors have begun exploiting some of these concepts to solve quadratically constrained quadratic programs, in which sums of bilinear products often figure prominently. In particular, the authors of [START_REF] Bao | Multiterm polyhedral relaxations for nonconvex, quadratically constrained quadratic programs[END_REF] discuss how to dynamically generate facets of the convex hull of the sum of bilinear products in order to define a stronger relaxation of the original MINLP, and they report that strengthening the formulation with such inequalities can significantly improve the performance of BARON [START_REF] Sahinidis | Baron 8.1.1: Global optimization of mixed-integer nonlinear programs[END_REF], which by default uses only McCormick envelopes to exploit multilinear terms in defining convex relaxations. Even more recently, [START_REF] Luedtke | Some results on the strength of relaxations of multilinear functions[END_REF] provides rigorous bounds for how much the approach of [START_REF] Sahinidis | Baron 8.1.1: Global optimization of mixed-integer nonlinear programs[END_REF] (and, implicitly, of [START_REF] Sherali | Convex envelopes of multilinear functions over a unit hypercube and over special discrete sets[END_REF]) can strengthen the relaxations defined by the use of McCormick envelopes, and also provides numerical results illustrating that these bounds are tight.

It is important to note that the bounds defined by [START_REF] Luedtke | Some results on the strength of relaxations of multilinear functions[END_REF] apply only to problems that have sums of bilinear products, but not quadratic terms (i.e., if the quadratic function in a given constraint is represented by f Q (x) = x T Qx, the bounds defined in [START_REF] Luedtke | Some results on the strength of relaxations of multilinear functions[END_REF] are valid for problems in which the diagonal elements of Q are all 0). Moreover, computational experience seems to confirm that the smaller the absolute values of elements on the diagonal of Q are in comparison to the off-diagonal elements, the more important the role played by strong convex relaxations for bilinear functions becomes in defining strong relaxations for the MINLP. (Defining effective relaxations for nonconvex quadratically constrained problems in which the diagonal elements of Q are large requires, in addition to the techniques described in this section, other methods that are fundamentally different. References that discuss solving nonconvex quadratically constrained problems with large diagonal absolute values include [START_REF] Anstreicher | Semidefinite programming versus the reformulation-linearization technique for nonconvex quadratically constrained quadratic programming[END_REF][START_REF] Bao | Multiterm polyhedral relaxations for nonconvex, quadratically constrained quadratic programs[END_REF][START_REF] Kim | Second order cone programming relaxation of nonconvex quadratic optimization problems[END_REF][START_REF] Saxena | Convex relaxations of non-convex mixed integer quadratically constrained programs: Extended formulations[END_REF][START_REF] Saxena | Convex relaxations of non-convex mixed integer quadratically constrained programs: Projected formulations[END_REF], and the references contained therein.)

An unresolved issue that is directly related to much of the research on multilinear functions described above is the question of whether or not it is possible to define a description of the convex envelope of multilinear functions that does not require the explicit a priori enumeration of all of the extreme points of the domain. More formally, given an n-dimensional function f (x) = ∏ n i=1 x i over a domain B = [x L , x U ], is it possible to define a set of criteria that 1) each facet of the convex envelope must satisfy, and 2) can be checked in time polynomial in n? Most of the approaches described above, as well as the motivation of this article, are based on the implicit assumption that the answer to this question is no. However, only a comparatively small number of research efforts (e.g., [START_REF] Sherali | Convex envelopes of multilinear functions over a unit hypercube and over special discrete sets[END_REF][START_REF] Luedtke | Some results on the strength of relaxations of multilinear functions[END_REF]) have addressed this question directly. Moreover, their consideration of this question has been limited to establishing criteria for x L and x U that are sufficient to guarantee that the answer is yes.

Computational complexity theory, and in particular results of [START_REF] Karp | On linear characterizations of combinatorial optimization problems[END_REF] suggests that a short (i.e., polynomial in n) description of the convex envelopes of multilinear functions can be defined if and only if the following optimization problem is polynomial solvable:

min n ∏ i=1 x i - n ∑ i=1 c i x i (1) s.t. x L i ≤ x i ≤ x U i , i = 1, ..., n, (2) (3)
where c ∈ R n is some rational vector. It seems that this problem is likely to be NP-complete unless fairly restrictive assumptions on x L and x U are satisfied. For example, generalizing some of the results of [START_REF] Luedtke | Some results on the strength of relaxations of multilinear functions[END_REF], in [START_REF] Namazifar | Valid inequalities, separation, and convex hulls for bounded multilinear functions[END_REF] the authors show that it is possible to solve the above optimization problem in polynomial time if there exists a constant a < 1 such that ax L i = x U i for i = 1, ..., n. It is also clear that slightly more general conditions can be established. However, the authors of [START_REF] Namazifar | Valid inequalities, separation, and convex hulls for bounded multilinear functions[END_REF] conjecture that the above optimization problem in NP-complete in general, and the complexity of this problem remains an important open question in the area of how best to approximate the convex envelopes of functions involving multilinear terms.

We will next turn to the general question of when, and how, one approach to defining convex relaxations of factorable functions can be shown to yield relaxations that are stronger than those generated by another approach. The primary contribution of this article is to establish a general result concerning this issue. We should perhaps first note, however, that this contribution does not tell us how much stronger the dominant formulation will be; this is necessarily an empirical question. Moreover, the comparative ease with which different relaxations can be solved is also a necessarily empirical criterion, and in general both of these considerations must be weighed in considering relaxation to use in a given situation.

The composition of convex envelopes

In this section we prove that a stronger relaxation is obtained when one replaces "large terms" with tight convex relaxations instead of breaking up such terms in sums/products of smaller terms before replacing each small term with its respective convex relaxation. Although we find that this is quite an intuitive result, because of the inherently recursive nature of factorable functions and of the fact that we deal with a recursive symbolic procedure for constructive the convex relaxation, we did not find it easy to prove this result formally. For this purpose, we use theoretical tools that are well known to the formal languages community but perhaps not so commonly found in the optimization literature: this is why we detail every step and attempt to be somewhat didactical in presentation, alternating formal statements to informal explanations and examples. To the well-versed in such matters, a brief glimpse to the section might suffice to understand our strategy: assign a special semantic value (the corresponding convex relaxation) to each operator node of an expression tree, define the semantics of the composition operator, and finally compare the resulting relaxation with the tight convex relaxation given for the composite operator at "atomic" level.

Alphabets, languages and grammars

An alphabet A is a set of symbols. We let A * be the set of all finite sequences of elements of A . A formal language L is a subset of A * . A language L is decidable if, given a string s ∈ A * , there exists a finite algorithmic procedure that decides whether s ∈ L or not.

Informally, decidability of a language is concerned with its syntax: is a string a valid element of the language or not? Having decided what a language is, we have to decide what it says: to every string there corresponds a semantic value, which, in the theory and language of Zermelo-Fraenkel, is usually a set. In this setting, our formal language is the set of all valid functions f (x) that can be written as finite strings of simbols in infix notation. The semantic values assigned to f (x) are sets such as

{(w, x) ∈ R n+1 | w = f (x) ∧ x L ≤ x ≤ x U } (exact semantics) and {(w, x) ∈ R n+1 | w ∈ R( f , x L , x U ) ∧ x L ≤ x ≤ x U } (relaxed semantics) where R( f , x L , x U
) is a convex relaxation of the exact semantics. Since the cardinality of our language is countably infinite, we cannot explicitly assign exact/relaxed semantics to each function in the language. Instead, we recall that a decidable language has finite procedure for recognizing strings in the language: for each of the (finitely many) operations specified by this procedure we define a corresponding operation on the semantic values involved, thus obtaining a semantic definition for the whole language.

To this effect, we make use of possibly the best known device for specifying the syntax of a formal language L , i.e. a formal grammar. This is a quadruplet Γ = (Σ , N, P, S) such that:

• Σ ⊆ A is the set of terminal symbols • N is a set of nonterminal symbols (N ∩ Σ = / 0) • P is a set of rewriting, or production rules (P ⊆ (Σ ∪ N) * N(Σ ∪ N) * → (Σ ∪ N) * ) • S ∈ N is the start symbol.
In practice, one recursively applies the production rules to the start symbol as many times as possible, generating strings in (Σ ∪ N) * . Those generated strings that are in Σ * are strings of the language L . If a string in A * is not in the set of all strings in Σ * that the grammar generates, then it is not in L . From this, we conclude that aa, ab, ba, bb are in the language specified by the grammar. It must be remarked that formal grammars can also be given for languages which are not decidable (for example if the recursion does not terminate); this is one such grammar: the repeated application of p 3 yields longer and longer strings all involving the nonterminal symbol S. ⊓ ⊔

Mathematical expression language: syntax

We now formally define our function language through the use of a formal grammar. We use an alphabet A = X ∪ K ∪ B ∪ O where X = {x 1 , . . . , x n } is the set of symbols denoting original variables, K is the set of all computable numbers, B = {'(', ')'} and O is a finite set of operators {+, -, ×, ÷, ↑, √ , log, exp, sin, cos, tan}, where +, × are binary operators,can be unary or binary and ↑ is the (binary) power operator. The grammar Γ is defined as follows. The start symbol is F , N = {F }, Σ = A , and P is:

F -→ x i ∈ X (4) F -→ k ∈ K (5) F -→ (F ) (6) F -→ (-F ) (7) F -→ log(F ) (8) F -→ exp(F ) (9) F -→ sin(F ) (10)
F -→ cos(F ) (11) F -→ tan(F ) (12) F -→ (F -F ) (13) F -→ (F ÷ F ) (14) F -→ (F ↑ F ) (15) F -→ (F + F ) (16) F -→ (F × F ) (17)
Notice that rules (4)-( 5) are given in schematic form: i.e. the string on the left of the arrow is not in (Σ ∪ N) * , but it is possible to define "sub-languages" that decide whether a string is in X or in K.

Example 2. In order to recognize that the string

F ≡ x 1 + ((x 2 ↑ 2) + (x 3 × (x 4 × log(x 1 )))
) is in L we can apply the production rules as follows (there are other possible orders in which the rules can be applied yielding the same result):

F -[1] → (F + F ) by (16) -[2] → (F + (F + F ))
by ( 16)

-[3] → (F + ((F ↑ F ) + F ))
by ( 15)

-[4] → (F + ((F ↑ F ) + (F × F )))
by ( 17)

-[5] → (F + ((F ↑ F ) + (F × (F × F ))))
by ( 17)

-[6] → (F + ((F ↑ F ) + (F × (F × log(F ))))) by (8) -[7] → (x 1 + ((x 2 ↑ F ) + (F × (F × log(F ))))) by (4) -[8] → (x 1 + ((x 2 ↑ F ) + (F × (F × log(F ))))) by (4) -[9] → (x 1 + ((x 2 ↑ F ) + (x 3 × (F × log(F ))))) by (4) -[10] → (x 1 + ((x 2 ↑ F ) + (x 3 × (x 4 × log(F ))))) by (4) -[11] → (x 1 + ((x 2 ↑ F ) + (x 3 × (x 4 × log(x 1 ))))) by (4) -[12] → (x 1 + ((x 2 ↑ 2) + (x 3 × (x 4 × log(x 1 )))))
by ( 5)

-[13] → x 1 + ((x 2 ↑ 2) + (x 3 × (x 4 × log(x 1 ))))
by [START_REF] Cafieri | On convex relaxations of quadrilinear terms[END_REF].

We need to apply 13 rewriting rules in order to recognize that F ∈ L . ⊓ ⊔

Mathematical expression language: semantics

We are now going to use the formal grammar Γ to assign semantic values to strings. Informally, we assign different sets to the different occurrences of the symbol F in each production rule, in such a way that the set assigned to F appearing in the left hand side of each rule is defined in terms of the sets assigned to the symbols F appearing in the right hand side. More precisely, for a production rule ρ in ( 4)-( 17) of the form F → T , where T ∈ (Σ ∪ N) * , let ν(ρ) be the number of occurrences of the symbol F in the string T . Let X 0 (ρ) be the set assigned to the symbol F appearing on the left hand side of ρ, and for all i ∈ {1, . . . , ν(ρ)} let X i (ρ) be the set assigned to the i-th occurrence of the symbol F in T .

Exact semantics

The exact semantics of L is defined according to the following rules.

F -→ x i ∈ X : X 0 = [x L i , x U i ] F -→ k ∈ K : X 0 = {k} F -→ (F ) : X 0 = X 1 F -→ (-F ) : X 0 = {(w, x) | w = -x ∧ x ∈ X 1 } F -→ log(F ) : X 0 = {(w, x) | w = log(x) ∧ x ∈ X 1 } F -→ exp(F ) : X 0 = {(w, x) | w = exp(x) ∧ x ∈ X 1 } F -→ sin(F ) : X 0 = {(w, x) | w = sin(x) ∧ x ∈ X 1 } F -→ cos(F ) : X 0 = {(w, x) | w = cos(x) ∧ x ∈ X 1 } F -→ tan(F ) : X 0 = {(w, x) | w = tan(x) ∧ x ∈ X 1 } F -→ (F -F ) : X 0 = {(w, x 1 , x 2 ) | w = x 1 -x 2 ∧ ∀i ∈ {1, 2} x i ∈ X i } F -→ (F ÷ F ) : X 0 = {(w, x 1 , x 2 ) | w = x 1 /x 2 ∧ ∀i ∈ {1, 2} x i ∈ X i } F -→ (F ↑ F ) : X 0 = {(w, x 1 , x 2 ) | w = x x 2 1 ∧ ∀i ∈ {1, 2} x i ∈ X i } F -→ (F + F ) : X 0 = {(w, x 1 , x 2 ) | w = x 1 + x 2 ∧ ∀i ∈ {1, 2} x i ∈ X i } F -→ (F × F ) : X 0 = {(w, x 1 , x 2 ) | w = x 1 x 2 ∧ ∀i ∈ {1, 2} x i ∈ X i }
A meta-linguistic note: the naming of the semantic values X 0 , X 1 , X 2 must be local to each rule. Otherwise, if the same rule ρ is applied twice, we might get two different definitions assigned to the same name X 0 (ρ). In order to obtain a consistent naming, we observe that the recursive nature of string recognition in L is finite, so the different strings of (Σ ∪N) * generated during the recognition procedure can be listed in the order of rewriting, as in Example 2. For a string f ∈ L let r( f ) be the length of this list. For all k ≤ r( f ), we can now let X k 0 be the semantic value assigned to F appearing in the left hand side of the production rule ρ being applied at the k-th rewriting step, and let X k 1 , . . . , X k ν(rho) be the sets assigned to the various occurrences of F in the right hand side of ρ.

As will appear clear in Example 3, some of the semantic sets will be projections of other semantic sets on some of their coordinates. For every semantic set X we shall therefore let V (X) be the sequence of variable symbols in terms of which X is defined (so that X ⊆ R |V (X)| ), and for all W ⊆ V (X) let π(X,W ) be the projection of X on the w coordinate (if W = {w}, we write π(X, w)).

which means that convex relaxations should get tighter when the definition intervals get smaller.

We remark that R is a symbol in the metalanguage, in the sense that it should be replaced by an actual description of the convex sets assigned to each operator (in other words, it stands for the sentence "for all possibly ways of defining convex relaxations of operators. . . "). A typical definition of R Γ used by most sBB solver codes (e.g. ooOPS [13] and Couenne [START_REF] Belotti | Branching and bounds tightening techniques for non-convex MINLP[END_REF], both based on a grammar very similar to Γ ) is as follows: for all linear operators, R Γ applied to that operator is the same as the exact semantics (because, as an affine space defined over a cartesian product of intervals, it is convex). The log, exp operators are concave/convex univariate, and hence R Γ is defined as a convex subset of R 2 delimited by the function itself and the secant at the interval endpoints [START_REF] Liberti | Writing global optimization software[END_REF]; for piecewise convex/concave functions we employ the convex envelope defined in [START_REF] Liberti | Convex envelopes of monomials of odd degree[END_REF]; for trigonometric functions it is easy to work out convex relaxations/envelopes using secants and convex/concave portions of the functions themselves. We remark that providing convex/concave relaxations/envelopes of convex/concave functions and piecewise convex/concave functions suffices to define R Γ over all univariate monomials of the form x k where x ∈ I ∈ I. For bilinear products, we employ the well-known McCormick envelopes:

R(×, [w L 1 , w U 1 ], [w L 2 , w U 2 ]) = {(w 0 , w 1 , w 2 ) | w 0 ≥ w L 1 w 2 + w L 2 w 1 -w L 1 w L 2 ∧ w 0 ≥ w U 1 w 2 + w U 2 w 1 -w U 1 w U 2 ∧ w 0 ≤ w L 1 w 2 + w U 2 w 1 -w L 1 w U 2 ∧ w 0 ≤ w U 1 w 2 + w L 2 w 1 -w U 1 w L 2 ) ∧ w 1 ∈ [w L 1 , w U 1 ] ∧ w 2 ∈ [w L 2 , w U 2 ]}.
It is easy to check that the above definition of R satisfies [START_REF] Meyer | Trilinear monomials with mixed sign domains: Facets of the convex and concave envelopes[END_REF]. The relaxed semantics of L is defined according to the rules:

F -→ ⊕(F , . . . , F ) : X 0 = R Γ (⊕, I 1 , . . . , I α(⊕) ).
Relaxed semantics can be combined following grammatic production rule composition in much the same way as exact semantics can, by noticing that when X is a convex subset of R n , the projection of X on one coordinate axis is always an interval (because projection preserves convexity). Now let F be a valid string of L : then F is a mathematical expression with, say, x = (x 1 , . . . , x n ) as variable symbol arguments corresponding to a certain mathematical function f : R n → R. Then we can certainly add the following rule to Γ :

ρ ′ ≡ F -→ F(F , . . . , F n ), (19)
yielding an extended grammar Γ ′ , and still obtain L as generated language. The advantage is that Γ ′ allows recognition of the string F in one step, and assignment of a special relaxed semantics to F (instead of relying on the composition of relaxed semantics of substrings of F through the production rules). This is useful for those operators which do not appear in the list of production rules, but for which we have a tight convex relaxation (or a convex envelope).

Comparison of relaxed semantics

Let F ∈ L represent an n-ary function such that ρ ′ , defined as in [START_REF] Meyer | Convex envelopes for edge-concave functions[END_REF], is not a production rule of Γ . Define A ′ as A ∪ {F} and Γ ′ as Γ with ρ ′ as an added production rule. Assume that the given relaxed semantics for F in Γ ′ is included in the computed relaxed semantics for F in Γ (which is usually the case in practice, for otherwise we would not add the "useless" rule ρ ′ to Γ ′ ), i.e. that, for all I 1 , . . . , I ℓ ∈ I,

R Γ ′ (F, I 1 , . . . , I ℓ ) ⊆ R Γ (F, I 1 , . . . , I ℓ ). ( 20
)
Theorem 1. For all strings T ∈ L that are functions of p variable symbol arguments and for all I 1 , . . . , I p ∈ I, we have R Γ ′ (T, I 1 , . . . , I p ) ⊆ R Γ (T, I 1 , . . . , I p ).

Proof. If recognition of T through Γ ′ never involves rule ρ ′ , both grammars yield the same relaxed semantics. Otherwise, consider the last time that ρ ′ is used on T : then Γ ′ matches a string F which is an operator F of n arguments. Let J 1 , . . . , J ℓ be the relaxed semantics assigned to each of the n arguments. Since this the last time ρ ′ is used, each of the J i (i ≤ n) are the same whether we use Γ or Γ ′ , which means that, by [START_REF] Namazifar | Valid inequalities, separation, and convex hulls for bounded multilinear functions[END_REF], [START_REF] Meyer | Trilinear monomials with mixed sign domains: Facets of the convex and concave envelopes[END_REF], any relaxed semantics involving J Γ ′ will be contained in the same relaxed semantics with J Γ ′ replaced by J Γ . Thus, if the statement holds from the (k + 1)-st to the last time rule ρ ′ is used, the k-th time ρ ′ is used the argument intervals of the relaxed semantics in Γ ′ must be contained in the argument intervals of the corresponding relaxed semantics in Γ .

J Γ ′ = R Γ ′ (F, J 1 , . . . , J ℓ ) ⊆ R Γ (F, J 1 , . . . , J ℓ ) = J Γ . By
In particular, we have the following.

Corollary 1. If F(x 1 , x 2 , x 3 ) = x 1 x 2 x 3 and we assign to F the relaxed semantics given by the trilinear envelopes given in [START_REF] Meyer | Trilinear monomials with mixed sign domains: Facets of the convex and concave envelopes[END_REF][START_REF] Meyer | Trilinear monomials with positive or negative domains: Facets of the convex and concave envelopes[END_REF], the convex relaxation obtained through Γ ′ is at least as tight as that obtained through Γ for any mathematical function in L .

Proof. Assumption [START_REF] Namazifar | Valid inequalities, separation, and convex hulls for bounded multilinear functions[END_REF] holds by definition of convex envelope.

Computational results

In this section, we computationally evaluate the tightness of convex relaxations for quadrilinear monomials obtained combining bilinear and trilinear convex envelopes in different ways. Specifically, we consider relaxations of the following four sets:

S 222 = {(x, w) ∈ R 4 × R 3 | x i ∈ [x L i , x U i ] ∧ w 1 = x 1 x 2 , w 2 = w 1 x 3 , w 3 = w 2 x 4 }, S222 = {(x, w) ∈ R 4 × R 3 | x i ∈ [x L i , x U i ] ∧ w 1 = x 1 x 2 , w 2 = x 3 x 4 , w 3 = w 1 w 2 }, S 32 = {(x, w) ∈ R 4 × R 2 | x i ∈ [x L i , x U i ] ∧ w 1 = x 1 x 2 x 3 , w 2 = w 1 x 4 }, S 23 = {(x, w) ∈ R 4 × R 2 | x i ∈ [x L i , x U i ] ∧ w 1 = x 1 x 2 , w 2 = w 1 x 3 x 4 }.
In [START_REF] Cafieri | On convex relaxations of quadrilinear terms[END_REF] numerical experiments were carried out in order to evaluate the relative tighteness of the four considered relaxations. The comparison was mainly made in terms of volume of the corresponding enveloping polytopes (projected onto R 5 to have comparable results) on a set of randomly generated instances. It showed that the smallest values of volumes correspond to relaxations involving the composition of trilinear and bilinar envelopes, and in particular the best results for more than 80% of the considered instances were obtained using relaxation S 23 . Numerical experiments on some real-life problems were carried out using a bound evaluation algorithm, whose purpose is to assess the quality of the proposed alternative bounds for quadrilinear terms. This "partial sBB" algorithm at each branching step only records the most promising node and discards the other, thus exploring a single branch up to a leaf. The best bounds were obtained using a relaxation involving a trilinear envelope.

In the present paper, we further investigate the strenght of the proposed relaxations in a sBB algorithm. To that effect, we implemented the computation of the four relaxations for quadrilinear monomials in COUENNE [START_REF] Belotti | Branching and bounds tightening techniques for non-convex MINLP[END_REF]. Computational experiments were carried out running COUENNE on 7 instances of the Molecular Distance Geometry Problem (MDGP) [START_REF] Liberti | Molecular distance geometry methods: from continuous to discrete[END_REF], the problem of finding an embedding x : V → R 3 of the vertices V of a weighted graph G = (V, E) such that all the edge weights d uv (for {u, v} ∈ E) are equal to the Euclidean distances x ux v . The MDGP mathematical programming formulation is:

min x ∑ {u,v}∈E x u -x v 2 -d 2 uv 2 , (21)
a nonconvex NLP involving polynomials of fourth degree. In our experiments we impose a time limit equal to 4 hours. Results were obtained on a 2.4 GHz Intel Xeon CPU of a computer with 8GB RAM shared by three other similar CPU running Linux. For the smallest MDGP instance, the optimal solution is computed within the time limit using all the considered relaxations. A comparison of CPU time is reported in Table 1 and shows that the time needed to solve the problem when relaxation S 23 is used is 81% smaller than the time needed using S 222 , that is the second best time to solve the problem. For the other instances, for which the optimal solution is not reached within the time limit, we compare the (lower) bounds obtained with the four relaxations. Results are shown in Table 2. These results confirm the results obtained in [START_REF] Cafieri | On convex relaxations of quadrilinear terms[END_REF]. It appears that the best bounds are always obtained using a relaxation involving a trilinear envelope and, in 5 cases out of 6, correspond to relaxation S 23 . The sBB based on this relaxation gives bounds which are significantly better than the ones obtained using a relaxation based on the composition of bilinear envelopes, in particular on the largest instances. For the first instance in Table 2 the optimal solution is found with relaxations S 222 and S 23 23 on MDGP instances. Bounds were obtained within a 4h time limit. The best values are reported in bold face. The symbol (*) denotes optimal solutions found. Solutions were obtained on a 2.4 GHz Intel Xeon CPU of a computer with 8GB RAM shared by three other similar CPU running Linux.

Conclusion

We analyzed four different convex relaxations for quadrilinear monomials, obtained by the composition of the known convex envelopes for bilinear and trilinear monomials. Starting from theoretical as well as computational results given in [START_REF] Cafieri | On convex relaxations of quadrilinear terms[END_REF], we further investigated these relaxations. We provided an alternative proof of the fact that a relaxation of k-linear terms that employs a successive use of relaxing bilinear terms (via the bilinear convex envelope) can be improved by employing instead a relaxation of a trilinear term (via the trilinear convex envelope). We computationally evaluated the impact of the tightened convex relaxations in a spatial Branch-and-Bound algorithm on a set of instances of a real-life problem.

Example 1 .

 1 Consider the alphabet {a, b} and the grammar given by N = {S} where S is the start symbol, Σ = {a, b} and the production rules p 1 = aS → bS, p 2 = Sb → Sa, p 3 = S → aSSb, p 4 = SS → / 0 . We repeatedly apply p 1 , . . . , p 4 to the start symbol, obtaining the situation below:

Table 1

 1 within the time limit. It took 8311.97 seconds in the first case and 7063.73 seconds in the second one. Comparison of CPU time (seconds) obtained by running couenne with relaxations S 222 , S222 , S 32 , S 23 on the smallest MDGP instance. The best value is reported in bold face. Solutions were obtained on a 2.4 GHz Intel Xeon CPU of a computer with 8GB RAM shared by three other similar CPU running Linux. ( * ) 199.864 200.45 228.574 ( * ) lavor6 93.4905 135.899 84.9467 144.399 lavor7 2.75184 90.3962 70.9786 207.255 lavor8 24.5401 95.0223 36.421 334.968 lavor10 -266.843 -105.584 -91.4539 93.6579 lavor20 -1571.58 -1215.7 -589.636 -1146.5

	Instance S 222	S222	S 32	S 23
	lavor3 311.934 372.306 475.835 58.0872
	Instance	S 222	S222	S 32	S 23
	lavor5 228.574

Table 2

 2 Comparison of lower bounds obtained by running couenne with relaxations S 222 , S222 , S 32 , S

The second and the fourth authors gratefully acknowledge financial support under ANR grant 07-JCJC-0151.

Example 3. The exact semantics of F, as defined in Example 2, is derived as follows.

Replacing symbols where possible, we obtain a definition of the exact semantics of our string in function of only six sets and ten variables (4 original variables and 6 added variables):

Suppose now we consider an enriched alphabet A ′ with one more 4-ary operator ⊗ such that ⊗(x 1 , . . . , x 4 ) = x 1 + x 2 2 + x 3 x 4 log(x 1 ), and an extended grammar with one more production rule ρ ′ ≡ F -→ F + F ↑ 2 + F × F log(F ). The generated language L ' is identical to L because we showed previously that L contains strings as that appearing in the right hand side of ρ ′ even without the production rule ρ ′ . However, using the extended grammar, the string F can be recognized in only one step. By replacement of the appropriate variable symbols w ℓ , the exact semantics {(w,

with the extended grammar is precisely the projection of X 1 0 on the subspace of R 10 spanned by (w 1 , x 1 , . . . , x 4 ). ⊓ ⊔

Relaxed semantics

We now define the relaxed semantics of L . Whereas in the exact semantics we assigned to each string the set of values taken by the corresponding function as its arguments range in the appropriate (recursively defined) sets, the relaxed semantics assigns to strings convex relaxations of such sets. To this end, we shall describe an operator R Γ that computes the convex relaxation of a set using the composition of production rules in Γ . For each operator ⊕ ∈ O, let α(⊕) be its arity (the number of its arguments). Denote α(⊕) by ℓ, I the class of all closed and bounded intervals in R, and let I 1 , . . . , I ℓ ∈ I; then we use the notation R Γ (⊕, I 1 , . . . , I ℓ ) to indicate a convex relaxation in R n+1 of the exact semantic value of ⊕, i.e. the set {(w 0 , w 1 , . . . , w ℓ ) | w 0 = ⊕(w 1 , . . . , w ℓ ) ∧ ∀i ≤ n (w i ∈ I i )}. We impose a consistency (monotonicity) requirement: ∀⊕ ∈ O, I 1 , . . . , I ℓ , J ∈ I s.t. ∃i ≤ ℓ J ⊆ I i R Γ (⊕, I 1 , . . . , I ℓ ) ⊇ R Γ (⊕, I 1 , . . . , I i-1 , J, I i+1 , . . . , I ℓ )), [START_REF] Meyer | Trilinear monomials with mixed sign domains: Facets of the convex and concave envelopes[END_REF]


Andreas Savin

Towards a systematic way to correct density functional approximations

In order to simulate the exact universal density functional, approximations are nowadays constructed by permitting more flexibility in its ansatz. In view of the difficulty of defining a systematically improvable form for it, this paper argues that an alternative way could be considered. It falls within the class of hybrid functionals with multi-determinant wave functions. The parameter controlling the hybridization is considered as variable. The invariance of the exact result with respect to changes in this variable is used to introduce information about the system under consideration, and to correct the density functional result. The construction considered in this paper accelerates convergence from the model system to the physical one, in the vicinity of the latter. The method, at the present level of implementation, should be seen as a starting point for further development, and not necessarily as a computationally advantageous tool.

INTRODUCTION

The celebrated Hohenberg and Kohn paper [1] suggested that there was a systematic way of going beyond the local density approximation by using a gradient expansion. However, it turned out quickly that the first correction, the gradient expansion approximation, did not produce an improvement [2]. It took more than a decade to make progress, using generalized gradient expansions [3], and even more to arrive at hybrid methods [4]. Presently, the wide-spread strategy to improve density functionals is to use parameters (as a rule fitted on experimental data), and validate the newly produced functional by comparing the result with large sets of experimental reference data. The initial goal of a method providing a systematic improvement has thus been transformed into that of finding a more flexible ansatz that hopefully approaches the exact functional.

The aim of the present paper is less ambitious, in the sense that no approximation for the universal density functional is searched. However, it retains a density functional approximation as a starting point, and tries to improve on it systematically. Whether such an approach can lead to a practical method needs yet to be shown.

The approach considered in this paper consists in the following steps:

1. construct a model Hamiltonian, using a hybrid density functional approximation, 2. calculate the energy of this model system, and its first derivative with respect to a parameter defining the hybridization, 3. make an estimate of the correction to the density functional approximation using the first derivative.

This sequence can be continued by either using higher derivatives, or by considering new values of the parame- * Electronic address: andreas.savin@lct.jussieu.fr ter. Thus, an improvement is achieved by adding more information of the specific system under study, that can be provided by the family produced by the variation of the parameter of the model. The structure of the paper is the following. First, the method is described. Next, some simple numerical examples are given. Finally, other similar approaches are mentioned.

The central result of the paper is given by the equations 24 and 25.

METHOD Model Hamiltonian from density functional calculation

In hybrid density functional calculations, the ground state energy is computed via:

E 0 = min Ψ ( Ψ|T + V ne + W (µ)|Ψ + Ēhxc [n Ψ , µ]) (1)
T is the operator for the kinetic energy, V ne is the local one-particle operator for external potential (for the interaction between electrons and nuclei), and W is a local two-particle operator. Ēhxc is an approximation to the universal Hartree-exchange-correlation functional, calculated from the density n. The subscript to n indicates that is obtained from the N -electron wave function Ψ. Notice the presence of the parameter µ. As W depends on it, Ēhxc depends not only on n, but also on µ.

In the Kohn-Sham method, W = 0. In the present formalism, as one considers W = 0, the minimizing wave function Ψ gets a multi-determinant character. For practical applications it will be important to use such forms of W that do not introduce large errors when only a few determinants are used. For example, already a weak interaction can switch on the coupling between nearly degenerate states of the non-interacting system. In this paper, Ψ has, for W = 0, multi-determinant character, although in literature the term "hybrid" is often used for the approximation in which Ψ is described by a single Slater determinant.

The stationarity condition with normalization constraint for Ψ gives a Schrödinger equation:

H(µ)Ψ(µ) = E(µ)Ψ(µ) (2)
The eigenfunction Ψ(µ) is the solution to the minimization problem eq 1. The Hamilton operator

H(µ) = T + V (µ) + W (µ) (3)
contains a local one-particle operator,

V (µ) = V ne + N i v hxc (r i , µ) (4)
the local potential v hxc being obtained from a functional derivative,

v hxc (r, µ) = δ Ēhxc [n, µ]/δn(r) (5)
In the following, it will be assumed that W → 0 for µ → 0, and that W → V ee , the electron-electron interaction, when µ → ∞. A correctly constructed approximation Ēhxc vanishes when µ → ∞, and in this limit, E 0 becomes the exact ground state energy. Please notice that the model system described by the above Hamiltonian also produces excited states, and that in the limit µ → ∞ these become the exact excited states.

For exemplification,

W (µ) = i<j erf(µr ij )/r ij (6)
is used in this paper; r ij is the distance between electrons. A local density approximation (µ-dependent LDA [5][START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF][START_REF] Paziani | [END_REF]) is used for the exchange-correlation functional.

Improving the estimate of the exact energy

As the potential V (µ) is known once an approximation was chosen for Ēhxc (eqs 4, 5), we can assume that E(µ), eq 2, can be obtained accurately. For example, for W = 0 (in a Kohn-Sham calculation using some approximate exchange-correlation potential) it is simply a sum over occupied orbital eigenvalues. As E(µ) is, in general, not equal to the exact energy, E(∞), a correction

Ē(µ) = E(∞) -E(µ) (7)
is needed. The correction given by the density functional approximation,

Ē(µ) = E 0 -E(µ) (8) = Ēhxc [n Ψ(µ) , µ] -d 3 r n(r)v hxc (r, µ) (9)
Ē(µ) is only an estimate of Ē(µ) because Ēhxc used in eq 1 is not exact.

We discuss now how to obtain some information about the correction, Ē(µ), without using the density functional approximation Ēhxc for it. We choose to make the correction accurate for large, but finite µ by using the derivative of Ē(µ) with respect to µ, Ē (µ). We first notice that the Hellmann-Feynman theorem permits a relatively cheap accurate determination of E (µ) without re-computing Ψ(µ),

E (µ) = Ψ(µ)|∂ µ V (µ) + ∂ µ W (µ)|Ψ(µ) (10)
Second, we notice that E(µ) + Ē(µ) = E(∞) does not depend on µ, and thus, although Ē(µ) is not known, accurate information about Ē (µ) is available through the derivative of E(µ) with respect to µ, E (µ),

E (µ) + Ē (µ) = 0 (11)
We can gain in accuracy if we posses information about the behavior of Ē(µ). For example, assume that around a given µ, we know how the exact correction behaves,

Ē(µ) = cf (µ) + . . . ( 12
)
where f is a known function, and the constant c is system-dependent. In particular, we not only know that for W of eq 6, at large µ, Ē(µ) vanishes, but also know the decay rate with µ (see, e.g., [5,8,9]),

Ē(µ) = cµ -2 + O(µ -3 ) (13)
We now use the derivative of Ē(µ) with respect to µ,

Ē (µ) = cf (µ) + . . . (14)
to extract the coefficient c = Ē (µ)/f (µ), and obtain [10] Ē

(µ) = f (µ) f (µ) Ē (µ) + . . . (15)
For our specific example, W of eq 6, at large µ,

Ē (µ) = -2cµ -3 + O(µ -4 ) (16)
and

Ē(µ) = - 1 2 µ Ē (µ) + O(µ -3 ) = 1 2 µE (µ) + O(µ -3 ) (17)
where eq 11 has been used for the last equality.

Let us now consider the behavior at large µ when approximations are used, Ē(µ → ∞), eq 9. For our choice of interaction, eq 6, we know that the correction to the electron-electron interaction term,

Ψ(µ → ∞)|V ee -W (µ → ∞)|Ψ(µ → ∞) is proportional to µ -2 Ψ(∞)| i =j δ(|r i -r j |)|Ψ(∞) (18)
i.e., to µ -2 , and to the system average of the on-top pair density of the exact system (see, e.g., [8,9]). Let us consider µ-dependent approximations for Ēhxc that model pair-densities (like LDA). In them, the error originates from using the system-average of the model pair density, instead of that of the system under consideration. Thus, Ēhxc is proportional to µ -2 ,

Ēhxc ∝ µ -2 + O(µ -3 ) (19)
Eq 9 also contains v hxc , eq 5. The proportionality to µ -2 remains after taking the functional derivative wrt the density, in order to produce v hxc . It follows that Ē also decays as µ -2 , and we can use, for Ē(µ → ∞), the same expression as for the exact correction, eq 17, viz.,

Ē(µ) = - 1 2 µ Ē (µ) + O(µ -3 ) (20)
In this expression, too, the derivative has introduced the information about the proportionality constant.

In order to analyze the error of Ēhxc when µ → ∞ we decompose it into a Hartree, an exchange and a correlation part. In this regime, the approximation probes the on-top pair density of the model, cf. eq 18. The Hartree part is evidently given by the density. At coalescence, the reduced first-order density matrix becomes equal to the density, and thus, for a closed shell, the exchange part of the pair density is also determined by the density. For a closed shell, the only error shows up in the correlation part of the pair density [5,9,11]. For open shells, however, the exchange component is not determined by the density alone, but by the spin-components of the density. This is why density functional approximations rely not on the density alone, but also on spin-density. It is wellknown, however, that the dependence on spin-density, introduces new problems related to size-consistency (see, e.g., [START_REF] Savin | Recent Developments of Modern Density Functional Theory[END_REF][START_REF] Perdew | Density Functional Methods in Physics. Proceedings of a NATO ASI held 1983[END_REF][START_REF] Cohen | [END_REF]). One thus has the dilemma of either seriously losing numerical accuracy, or size-consistency. In the present paper, spin-density is not used, in order to favor size-consistency. For describing the correlation part of Ēhxc errors are present. In LDA, it is modeled using the correlation hole of the uniform electron gas. This is a good approximation, but only an approximation as becomes evident, e.g., for one-electron systems.

We can try to introduce some correct information about the system obtained from the energy derivative, as in eq 17, to replace the false one contained in the model system used in the density functional approximation. First, we consider expanding both Ē and Ē in some (complete) basis,

Ē(µ) = i ci χ i (µ) (21) Ē(µ) = i ci χ i (µ) (22)
We now identify the basis function, χ j , that is responsible for the behavior we are interested in (at large µ).

As we can consider that we can obtain its coefficient cj accurately, we replace the approximate coefficient in the expression of Ē by an accurate estimate of one of the coefficients, cj ,

Ē(µ) ≈ Ē(µ) + (c j -cj )χ j (µ) (23)
To determine the coefficients, we use the derivatives, eq 15, and write

Ē(µ) ≈ Ē(µ) + χ j (µ) χ j (µ) Ē (µ) -Ē (µ) (24)
When choosing W as in eq 6, and making the approximation valid for large µ, by eqs 17 and 20,

Ē(µ) ≈ Ē(µ) - 1 2 µ Ē (µ) -Ē (µ) (25)
The same equation could be produced by requesting that the approximation deviates minimally from the density functional approximation, in a least-squares sense, with the constraint that the derivative of the approximation with respect to µ corresponds to the exact one.

Please notice that in the present construction, density functional approximations were used to construct model systems. The consistency between the correction to the energy and the potential present in density functional theory is not needed. Thus, the present approach can be applied to ground state, and to excited states, as well.

The derivative is an additional information from which we can infer from E in a given point µ the behavior in the vicinity of it. This procedure can be continued with higher derivatives, or in a computationally more convenient way, repeated for different µ points, in order to improve the accuracy in a systematic way. However, in the present paper, only the first correction, using only E (µ) in a single point will be presented.

NUMERICAL RESULTS

In this section, the approximation of eq 25 is tested for some two-electron systems. We know that this approximation becomes exact for large µ. The computations explore the behavior outside this regime. Its importance is related to the ansatz made for W : as long as W is weak, the computational effort is expected to be small.

In order not to bias the results, relatively large atomic basis sets were used (cc-pV5Z [14]). Furthermore, the wave function Ψ(µ) was computed with a full configuration interaction calculation (range-separated configuration interaction program of the Molpro code [START_REF] Werner | Molpro, version 2009.1 a package of ab initio programs[END_REF][START_REF] Schütz | [END_REF][17][18][19][20][21][22][23][24]). All the plots show the error done by approximating Ē(µ) by using eq 25, instead of E(∞) -E(µ), for different values of µ. By construction, the errors vanish at large µ. As the correction is based on an expansion at large µ, one will not expect it to work for µ smaller than ≈ 1. As evident from eq 25, no correction to the density functional expression is present at µ = 0. The first example is that of the hydrogen molecule at the equilibrium distance (fig 1). It shows exactly the features expected. We can see in fig 1 that the approximation has a quality comparable to that short-range LDA, known to work well for larger values of µ.

However, when the hydrogen molecule is stretched, the broken symmetry solution and spin-density has to be introduced in LDA to describe the potential energy curve at LDA level. We see in fig 2 that keeping a correct (zero) spin-density, the LDA error is relatively large when the distance between the hydrogen nuclei increases. In this case, the approximation of eq 25 provides improved results. In the limit of dissociation, the error would be twice that obtained for hydrogen atoms, as sizeconsistency is respected. The next example treats an excited state of the hydrogen molecule, E, F 1 Σ + g . It corresponds to the resonance structure H + . . . H -↔ H -. . . H + . In this case, LDA is not expected to work, because electrons are brought together by correlation, while a transfer of the pair function from the ground state of the uniform electron gas keeps the electrons apart (see fig 3). Furthermore, symmetry breaking cannot help in its standard implementation. However, when using eq 25 the quality is comparable to that seen previously, as long as µ is not too small.

Finally, let us discuss a problem, known to show up for the forms of W (µ), where µ appears as multiplied by the distance between electrons, µr 12 [26]. In such an ansatz one expects a scaling of the error: the error observed for a given µ for the more diffuse system show up at a larger µ for a compact system (because, on average, r 12 is reduced). This can be illustrated by the behavior of the approximation for the He atom (that can be seen as the hydrogen molecule at zero internuclear separation), fig 4. One can see that the shape of the curve in fig 4 is similar to that of the hydrogen molecule at the equilibrium distance. The point where the large-µ approximation becomes critical is marked by oscillations. For He, they show up at a larger value of µ than for H 2 , and and they are more important (also for µ = 0, the Kohn-Sham system, the LDA errors are larger). The curve for the hydrogen molecule and that for the helium atom can be made to be almost superimposed by using a scaling factor of ≈ 2 for one of them.

As for the hydrogen molecule, the method does not bring any correction over the density functional approximation at µ = 0. However, one may notice an improvement over the LDA for small values of µ, also for He.

SUMMARY AND COMPARISONS

It is proposed to consider density functionals as models for producing model energies, E(µ), eq 2. To obtain E(µ), the potential produced by the density functional is used, not the functional itself. To obtain the energy of the physical system, a correction is needed. It is proposed to correct the approximate density functional expression by considering changes in the model obtained by changing µ, introducing this way system-dependent information (see the constant c in eq 12). In particular, one can use the derivative of the model energy with respect to µ, as in eq 15. The same procedure can be applied also for approximations, and the the wrong behavior of the approximation (in the regime considered) can be eliminated (eq 24, or eq 25). It is a first step in correcting approximate density functionals, by introducing more information from model systems.

The approximation introduced in this paper does not require the use of spin-density, and is not restricted to ground states. In order to produce numerical examples, a Hamiltonian containing long-range interaction (erf(µr 12 )/r 12 ) was used. The correction introduced becomes exact as this interaction approaches the physical one, 1/r 12 .

In the regime where the approximation becomes exact, one could as well use eq 17 (as in [10]). The advantage of eq 25 is that, in contrast to eq 17, a correction is active even for small µ. However, for the important limit µ = 0, there is no improvement over the density functional approximation.

Of course, this approach can be extended to density functionals that work better than LDA. One can thus expect improving the quality of the approximation when the interaction is weak.

Finally, it should be mentioned that the form of H(µ) used in this paper is not the only one of interest. For example one could consider model Hamiltonians of the form T (µ) + V (µ) + V ee where T (µ) is a non-local oneparticle operator. It has be shown that it can lead to good approximations [27][START_REF] Krieger | Electron Correlation and Material Properties[END_REF][START_REF] Ayala | [END_REF][30][31].

FIG. 1 :

 1 FIG.1: Error of the approximation, eq 25, as a function of the parameter of the model Hamiltonian (µ) for the hydrogen molecule, in its ground state, at the equilibrium distance (full curve), and error of the LDA functional (dashed curve). The horizontal dotted lines show the limits of "chemical accuracy" (±1 kcal/mol).

FIG. 3 :

 3 FIG. 3: Same as fig 1, but the hydrogen molecule is now in the E, F 1 Σ + g excited state, around the outer minimum of the potential energy curve, 4.2 bohr [25].

FIG. 4 :

 4 FIG. 4: Same as fig 1, but now for the helium atom.

ACKNOWLEDGMENT

J. Toulouse (UPMC Univ Paris 06, France) and the two anonymous referees have contributed by their questions and suggestions to improve the clarity of this paper.