Amory Bisserier

Sylvie Galichet

Fuzzy Piecewise Linear Regression

published or not. The documents may come

I. INTRODUCTION

Model identification is based on a general principle which consists in determining, among candidate functions, the function of the inputs that best explains the output according to a given criterion. Assuming a particular class of models such as linear functions, splines, rule-based systems, neural networks, ..., the best candidate is determined from the available information, usually a set of observations of the input and output variables. Classical identification techniques assume perfect knowledge of input and output values. It means that the observations are supposed to be both precise (point-valued) and certain. However, there are situations in which this assumption is not realistic, especially when the information about the output value is obtained through measuring devices with limited precision. In the framework of fuzzy modeling in which it is possible to handle imprecise representations using fuzzy set theory as proposed by Zadeh [START_REF] Zadeh | Fuzzy Sets[END_REF], the assumption of perfect data becomes even paradoxical. Nevertheless, most fuzzy model identification techniques used in practice, especially in fuzzy control, still consider crisp data. In this context, the main objective of the paper is to revisit some theoretical works about fuzzy regression techniques [START_REF] Diamond | Fuzzy regression analysis, Fuzzy sets in decision analysis, operations research and statistics[END_REF] and to propose some slight improvements for making their practical use in fuzzy model identification easier. According to [START_REF] Diamond | Fuzzy regression analysis, Fuzzy sets in decision analysis, operations research and statistics[END_REF], fuzzy regression techniques can be classified into two distinct areas. The first, proposed by Diamond [START_REF] Diamond | Fuzzy Least Squares[END_REF], is an adaptation of the classical least squares method and the second, called possibilistic regression, reduces the prob-lem of finding fuzzy coefficients of a regression model to one of mathematical programming. The latter, on which we will focus here, was introduced by Tanaka, Uejima and Asai [START_REF] Tanaka | Linear regression analysis with fuzzy model[END_REF] in a linear context. Several improvements on the original method can be found in the literature. For example, Tanaka, Hayashi and Watada [START_REF] Tanaka | Possibilistic linear regression analysis for fuzzy data[END_REF] propose different expressions of the criterion to be optimized and different formulations of the constraints to be satisfied for possibility and necessity estimation models. Still in a linear context, Tanaka and Ishibuchi [START_REF] Tanaka | Identification of possibilistic linear systems by quadratic membership functions of fuzzy parameters[END_REF] extend their approach for dealing with interactive fuzzy parameters. Furthermore, the complete specification of regression problems highly depends on the nature of input-output data [START_REF] Diamond | Fuzzy regression analysis, Fuzzy sets in decision analysis, operations research and statistics[END_REF]. Some works are thus devoted to crisp inputcrisp output data [START_REF] Savic | Evaluation of fuzzy linear regression models[END_REF] while others [START_REF] Sakawa | Mutiobjective fuzzy linear regression analysis for fuzzy input-ouput data[END_REF] consider fuzzy inputfuzzy output data. Most commonly, a mixed approach (crisp input -fuzzy output) is chosen [START_REF] Tanaka | Linear regression analysis with fuzzy model[END_REF]. That is the formalism we adopt here in a linear context with the idea of keeping a simple model, possibly invertible ( [START_REF] Boukezzoula | Nonlinear Internal Model Control: Application of Inverse Model Based Fuzzy Control[END_REF], [START_REF] Boukezzoula | Inverse Controller Design for Interval Fuzzy Systems[END_REF]). The expressiveness of linear regression models with fuzzy triangular parameters is limited. One weakness is the fact that the fuzziness of the model output varies in the same way than the absolute value of the inputs. It follows that it is impossible to have a decreasing (resp. increasing) spread of the model output for positive (resp. negative) inputs. This restriction is acceptable in a measurement context where it is usual to express percentage relative errors. However, when fuzziness is considered as an intrinsic characteristic of the system to be modeled, the assumption that the higher the input, the higher the fuzziness attached to the model output, is open to criticism. Finally, as classical fuzzy regression models are not able to represent any tendency of output spread, they become more imprecise than necessary in some situations. As a consequence, in piecewise fuzzy regression problems, in which collected data can have any kind of spread tendency, actual identification methods are clearly insufficient. For example, Yu, Tzeng and Li [START_REF] Yu | A general piecewise necessity regression analysis based on linear programming[END_REF], [START_REF] Yu | General fuzzy piecewise regression analysis with automatic change-point detection[END_REF] propose a piecewise model in which each fuzzy linear piece is constructed according to its difference with the previous piece but the above-mentioned problem is not tackled. Furthermore, from a computational point of view, the choice of a differential model form is not very convenient. According to above-discussed points, the originality of the presented work is triple:

1) A modified form of fuzzy linear model is proposed so that any kind of spread tendency can be represented.

Fuzzy piecewise linear regression

A. Bisserier, S. Galichet, R. Boukezzoula

2) The formulation of the linear program used to identify the model introduces a modified criterion that assesses the model fuzziness independently of the collected data.

3) These concepts are used in a global identification process in charge of building a piecewise model able to explain every kind of output evolution. This paper is organized as follows. Section II presents the proposed fuzzy linear model and then addresses its identification. Section III presents the global identification process of a piecewise model. Applications on several examples are shown in Section IV. Then, some discussion about the model validity is developed in Section V. Finally, conclusions and perspectives are presented in Section VI.

II. IDENTIFICATION OF A MODEL ON AN INTERVAL

A. Notations and concepts

We consider symmetrical triangular fuzzy numbers whose notation is (m,r) with m the modal value and r the radius. A representation of this kind of fuzzy numbers is shown in Fig. 1. Let us also introduce the notion of α-cut, with α a cut level between 0 and 1. The α-cut of the fuzzy number A is the interval [A] α composed of elements whose membership in A is greater or equal to α, i.e.:

(

The α-cut of A at α = 0 is the support of the fuzzy number, noted (see Fig. 1.).

B. Linear fuzzy model representation

The objective of this part is to propose a convenient form for a fuzzy linear model. For the sake of simplicity but without lack of generality, we consider here a single input model.

A classical fuzzy linear model with a single input is defined on its domain D = [x min ,x max ] by:

(

) 2
where x is the input variable, , Y(x) the corresponding output, and A 0 = and A 1 = the two fuzzy coefficients. The operator represents the sum of two fuzzy numbers, whose result is also a fuzzy number, defined by its modal value and its radius. As A 0 and A 1 are symmetrical triangular fuzzy numbers, and x a crisp input, Y(x) is also a symmetrical triangular fuzzy number. So, its modal value and spread, i.e. m Y(x) and r Y(x) are given by:

(

It is obvious from (3) that the variation of m Y(x) with respect to x only depends on . So, by choosing the sign of , it is possible to represent either increasing or decreasing m Y(x) with respect to x∈D. On the contrary, from ( 4), it appears that the variation of r Y(x) with respect to x only depends on the sign of x since is positive. It follows that when x is positive, r Y(x) is increasing, whereas when x is negative, r Y(x) is decreasing. So, in model ( 2), it is possible to have any kind of variation of the output modal value with an appropriate sign of , but the variation of the output radius is limited by the sign of the input x.

Consequently, if we want to obtain a regression model whose output can have any kind of radius variation for any sign of x, it is necessary to slightly modify the fuzzy model variable to be able to control its sign. With the idea of keeping a linear behavior, we propose to shift the original input x so as to obtain the desired sign for the shifted variable. So, by defining (xshift) as the new input variable of the fuzzy model, it becomes possible to control the radius variation by tuning the value of shift. Thus, we now consider a fuzzy linear model defined on its domain D = [x min ,x max ] by: [START_REF] Diamond | Fuzzy regression analysis, Fuzzy sets in decision analysis, operations research and statistics[END_REF] The output modal value and radius are given by: (6) [START_REF] Keogh | An Online Algorithm for Segmenting Time Series[END_REF] The output modal value variation with respect to x only depends on the sign of . So, increasing and decreasing tendencies can be represented. Concerning the variation of r Y(x) with respect to x, two cases are possible depending on the sign of . When an increasing r Y(x) is desired, it is necessary that , i.e. . As , it is required that . On the contrary, for a decreasing r Y(x) , it is necessary that , i.e. . This constraint has to be satisfied for any , which means that the shift value has to be chosen so that .

In the remaining, a shift convenient value is chosen. More precisely, the bound of the domain D that satisfies the constraint associated to the tendency desired for r Y(x) is selected. In other words, for a model whose output has an increasing radius with respect to x and for a decreasing radius (see TABLE I). It amounts to positioning the zero of the shifted model either at x min or at x max . By doing so, the output Y(x) of the fuzzy model ( 5) is equal to A 0 at one

A [ ] α x μ A x ( ) α ≥ { } = m 1 α - ( )r , m - 1 α - ( )r + [ ] = A [ ] 0 A -A + , [ ] = Membership value x 1 α 0 m A + = m+r A -= m-r m-(1-α )r m +(1-α )r interval [A] α Y x ( ) A 0 A 1 x ⊕ = x D ∈ m A 0 r A 0 , ( ) m A 1 r A 1 , ( ) ⊕ m Y x ( ) m A 0 m A 1 x + = r Y x ( ) r A 0 r A 1 x + = m A 1 m A 1 r A 1 m A 1 Y x ( ) A 0 A 1 x shift - ( ) ⊕ = m Y x ( ) m A 0 m A 1 x shift - ( ) + = r Y x ( ) r A 0 r A 1 x shift - + = m A 1 x shift - x shift - 0 ≥ x shift ≥ x x min x max , [ ] ∈ shift x min ≤ x shift - 0 ≤ x shift ≤ x D ∈ shift x max ≥ shift x min = shift x max =
bound of the domain D.

TABLE I THE TWO MODELS

The shifted model ( 5) can be characterized by its global fuzziness, possibly assessed as the area covered by the fuzzy output when the input varies on its domain [x min , x max ]. By denoting the shifted model input and [w min , w max ] its variation domain, the global fuzziness of the model is given by: .

(

) 8
The model identification procedure described in the next section is based on the following principle: considering two acceptable models, the one with the minimum fuzziness, i.e. the less imprecise is preferred.

C. Fuzzy model identification

The identification of a shifted fuzzy model ( 5) on a given domain D = [x min ,x max ] amounts to determine the shift value and the two fuzzy coefficients A 0 and A 1 . The latter being symmetrical triangular fuzzy numbers, it consists in finding their modal values and and their spreads and . A set of M data is available for the identification. In all input/ output observations , the input x j is a crisp value in D and the corresponding output Y j is a fuzzy number, noted where is the modal value and the radius. The M data are sorted in increasing order of x j .

The first step of the identification concerns the choice of the shift value according to the output radius tendency (see TA-BLE I). The most appropriate tendency is determined from observed data, comparing the initial output radius r init attached to minimal inputs with the final output radius r fin attached to maximal inputs. If r init < r fin , an increasing tendency is chosen, otherwise a decreasing tendency is preferred. The corresponding shift value is determined from TABLE I. The r init and r fin values are estimated by computing mean values from k data, that is r init = mean(e 1 , e 2 , ..., e k ) and r fin = mean(e M-k+1 , ..., e M-1 , e M ).

The next step of the identification concerns the optimization of the fuzzy coefficients A 0 and A 1 . A fuzzy regression approach based on possibilistic concepts as introduced by Tanaka [START_REF] Tanaka | Linear regression analysis with fuzzy model[END_REF] is used. It consists in solving a linear programming problem, with a criterion to be optimized and constraints to be satisfied. Hereinafter, the shifted variable is used in order to simplify the criterion expression and the constraint writing. Note that w j is either positive or negative, depending on the shift value.

If several kinds of constraints can be considered [START_REF] Tanaka | Possibilistic linear regression analysis for fuzzy data[END_REF], we want hereinafter the predicted intervals (resulting of the α-cut of fuzzy numbers) to fully cover the observed ones. Considering the support of observed fuzzy numbers, i.e. α = 0, the whole fuzziness of data is taken into account in the model identification. So, at α = 0, with the j th predicted support, the following constraint is obtained: [START_REF] Savic | Evaluation of fuzzy linear regression models[END_REF] otherwise expressed as:

. ( 10
)
The choice of the criterion to be minimized is also an important issue. The first one introduced by Tanaka, Uejima and Asai [START_REF] Tanaka | Linear regression analysis with fuzzy model[END_REF] is the sum of the spreads of the model coefficients, which leads to the following optimization problem:

Another approach, consisting in minimizing the sum of the spreads of the predicted intervals, also exists:

It is well established in the literature that the criterion ( 12) is more relevant than the criterion [START_REF] Tanaka | Identification of possibilistic linear systems by quadratic membership functions of fuzzy parameters[END_REF]. Indeed, the use of the criterion [START_REF] Tanaka | Linear regression analysis with fuzzy model[END_REF] improves the predictivity of the identified model. However, as the criterion ( 12) is only based on the available data, its minimization does not guarantee that the identified model has the least global fuzziness (8) that could be achieved on the whole domain D. If the identified model is to be used on the whole domain D, it may be more judicious to prefer a model with a lower global fuzziness, i.e. a less imprecise model. In this context, the global fuzziness (8) becomes the criterion to be minimized as proposed in [START_REF] Bisserier | Une vision de la régression linéaire floue au travers de l'arithmétique des intervalles[END_REF]. Actually, the global fuzziness of a model is an intrinsic characteristic of the latter, that should be assessed independently of the identification data, already used to express the constraints of the linear program.

According to this criterion choice, the following optimization problem is obtained:

In order to compare problems ( 12) and ( 13), both criteria are normalized. The normalization of (12) leads to: [START_REF] Yu | General fuzzy piecewise regression analysis with automatic change-point detection[END_REF] where is the average of , while the normalization of (13) results in: output radius variation model used

A 0 A 1 x x min - ( ) ⊕ A 0 A 1 x x max - ( ) ⊕ w x shift - ( ) = 2 w max w min - ( ) r A 0 w max 2 w min 2 - ( ) r A 1 + m A 0 m A 1 r A 0 r A 1 x j Y j , ( ) j , 1 … M , , = y j e j , (
)

y j e j w j x j shift - ( ) = Y ˆj Y ˆj - Y ˆj + , [ ] = y j e j -y j e j + , [ ] Y ˆj - Y ˆj + , [ ] ⊆ m A 0 m A 1 w j r A 0 r A 1 w j + + + y j e j + ≥ m A 0 m A 1 w j r A 0 r A 1 w j + ( ) - + y j e j - ≤ min r A 0 r A 1 + m A 0 m A 1 r A 0 r A 1 , , , min
Mr A 0 r A 1 w j j 1 = M ∑ + m A 0 m A 1 r A 0 r A 1 , , ,
min 2 w max w min - ( ) r A 0 w max 2 w min 2 - ( ) r A 1 + m A 0 m A 1 r A 0 r A 1 , , , min r A 0 r A 1 w + m A 0 m A 1 r A 0 r A 1 , , ,
w w j j , 1 … M , , = (15)
where is the middle point of the interval :

. (

So, if the identification data are uniformly distributed on the domain D, the two criteria are strictly equivalent. If it is not the case, solving problem ( 13) may lead to a model less imprecise than solving problem [START_REF] Tanaka | Linear regression analysis with fuzzy model[END_REF]. This approach can be interesting especially when M is small or when the domain D is bounded and clearly specified.

To sum up, the identification on D = [x min , x max ] of a model using a set of M observed data, is dealt with the following Identification function.

III. PIECEWISE LINEAR MODEL IDENTIFICATION

The fuzzy piecewise identification method is divided into two tasks: structure identification and parameter estimation. The structure identification consists in determining the number of submodels and their intervals of definition according to a segmentation method. In this context, several techniques have been developed to deal with this problem. Of course, any segmentation method can be applied. However, as already pointed out in [START_REF] Keogh | An Online Algorithm for Segmenting Time Series[END_REF], the quality of the model identification depends on the applied segmentation procedure. So, when a segmentation method is available, the identification problem is reduced to a parameters estimation one. That is the approach developed in this paper where the principle proposed by Keogh, Chu, Hart and Pazzani for segmenting time series [START_REF] Keogh | An Online Algorithm for Segmenting Time Series[END_REF] is used. The proposed identification method is used to identify a piecewise fuzzy linear model of the form (17):

(17)
where

, N is the number of submodels which compose the global model, represents the sum of fuzzy numbers. The function is equal to 1 on and to 0 anywhere else. In the identification method implementation some considerations must be taken into account. At the beginning of the process, it is necessary to find a good segmentation of the data set, i.e. to determine in which intervals to identify the different submodels. A changeover between two submodels is characterized either by an abrupt change in output modal values or by an abrupt change in output radius. So, in the segmentation procedure, both changes have thus to be considered. So, in order to get the different , k=1,..., N, the following strategy is applied:

-segmentation according to y j -on each interval got, segmentation according to e j . In this context, the Identification function is applied on each interval given by any segmentation process in order to determine the best model according to the global fuzziness criterion [START_REF] Sakawa | Mutiobjective fuzzy linear regression analysis for fuzzy input-ouput data[END_REF]. The advantage of the proposed method with regard to the one developed by Yu, Tzeng and Li [START_REF] Yu | General fuzzy piecewise regression analysis with automatic change-point detection[END_REF] resides in the two following points:

• from a structural point of view, according to the suitable shifts determination, any kind of output spread tendencies can be represented. • from a practical application point of view, the obtained piecewise model is not differential (the submodels can be identified and used independently). Let's have a few words about the continuity of the global model. Actually, the proposed identification procedure do not require the model continuity because looking for a continuous model may strongly increase the model fuzziness. Indeed, considering two consecutive segments, the more imprecise model induces an additional constraint on the other model when continuity is imposed. Moreover, there is no negative consequence on the use of the model, in case of inversion for example, if an 0-order holder is implemented at the change points.

IV. NUMERICAL EXAMPLES

First, we apply our identification method to the example proposed by Yu, Tzeng and Li [START_REF] Yu | General fuzzy piecewise regression analysis with automatic change-point detection[END_REF]. The model we get is the following (18):

(18) Its global fuzziness (8) is 82.5, whereas the Yu model [START_REF] Yu | General fuzzy piecewise regression analysis with automatic change-point detection[END_REF] has a global fuzziness equal to 99. The model ( 18) is clearly less fuzzy and so can be considered as better. In order to confirm the performance of the obtained models, other criteria have to be considered. In this context, the sum of the spreads of the predicted intervals (see [START_REF] Tanaka | Linear regression analysis with fuzzy model[END_REF]), named Sum, and the Diamond distance [START_REF] Diamond | Fuzzy Least Squares[END_REF] between the observed and the predicted fuzzy numbers, named Dist. are used. For this example, the obtained results are summarized in Table II For j = 1 to M Do w j = x jshift EndFor Solve (13) subject to ( 10) ----------------------------------= Y ˆx ( ) II, it can be stated that the obtained model presents a better fitting of the data, i.e., both Sum and Dist are lower. Then, we apply our identification method to the example proposed by Tanaka and Ishibuchi ([11], example 2). The model we get is the following (19) (see Fig. Finally, we apply our identification technique on synthetic time serie data, with several cases of modal value and spread tendencies. Used data are too numerous to be presented here, so we only explain how the data set is generated (see Table III). On each interval, the input x is increasing with an increment of 1. For each linear submodel, we set a variation for the modal value and another for the spread, then we add a random noise. Note that for the first two intervals, spread is globally constant, and so we only add a noise between 15 and 25 for the first one, and 5 and 10 for the other one, in order to get a change in spread level. The model identified from the data is presented in Table III and shown in Fig. 3. We can see that we get a good identification of a piecewise model on the data set. With the new kind of models with shift, we get a well suited representation of the observed data on each interval, and this, for any output spread tendency.

model = {shift, A 0 , A 1 } min r A 0 r A 1 w m + m A 0 m A 1 r A 0 r A 1 , , ,
w m w min w max , [ ] w m w max w min + ( ) 2 -
A 0 A 1 x shift - ( ) ⊕ = ⊕ Y ˆx ( ) A k0 A k1 x shift k - ( ) ⊕ [ ] 1 x min k x max k , [ ] k 1 = N ∑ = shift k x min k x max k , { } ∈ ⊕ ∑ 1 x min k x max k , [ ] x min k x max k , [ ] x min k x max k , [ ] Y ˆx ( ) 11

V. ABOUT THE VALIDITY OF THE MODEL

The identification method we propose here is based on possibilistic concepts. The main aspect of this kind of method is the use of a linear program. In the literature, several kinds of criteria or constraints are proposed in order to improve the identification of the model. However, all these methods have a similar aspect, i.e. the use of α-cuts in order to finally consider intervals in the optimization procedure. Choosing some α level, a modal value and a radius are determined for each coefficient of the model with the necessity of having inclusion of the α-cut of the observed output in the predicted one (see Fig. 4.). Then, the fuzzy model output is reconstructed, for each value of α between 0 and 1. So, it is possible for any value of the input to determine the corresponding fuzzy output. Two points appears clearly in Fig. 4. Whatever the chosen α level in the optimization problem, for any , the inclusion of the observed intervals in the predicted ones is respected. However, it is not the case anymore for any . Consequently, even if the inclusion is respected at a given α level, it is not the case for the whole fuzzy number after reconstruction.

An idea that can be exploited is to identify the model at α = 1. However, it is well known from the literature that the optimization problem has no solution when the observed modal values are not strictly lined up. Moreover, the higher the α considered for the identification, the wider the support of the predicted fuzzy number. That's why we use in the work presented here a linear program at α = 0, in order to get the most relevant model taking in care the whole fuzziness of the observed data, i.e. the model with the less wide support. However, in our opinion the reconstruction of a fuzzy set from an interval-based identification is not pertinent in a fuzzy point of view. However, as this approach is classically used in the literature, we adopt the same formalism in order to position our technique and results with respect to other existing methods. Knowing this important limit, we work at present on the identification of another form of fuzzy models whose output could contain all fuzzy observed data.

VI. CONCLUSION

With the proposed form of shifted models, it becomes possible to represent output spreads either increasing or decreasing with respect to inputs. Identifying such models leads to models whose fuzziness is possibly lower than usually. Moreover, in the context of piecewise modeling, the double segmentation, one for the modal values of the observed outputs and another for their spreads, contributes towards reducing the global fuzziness of the model. Further works concern the extension of the proposed approach to multi-input problems. Actually, the regression algorithm is easily adaptable, but the way the multidimensional input space should be segmented remains a crucial problem. Another point to be studied is the generalization of the identification procedure to fuzzy inputs, in order to manage several cases of uncertainty on collected data.

Y j [ ] α Y ˆj [ ] α α Y j Y ˆj 1 0 Y α α' α ≤ α'' α >

Fig. 1 .

 1 Fig. 1. Representation of a fuzzy number with an α-cut.

Fig. 2 .

 2 Fig. 2. Identified model on the example of Tanaka and Ishibuchi

Fig. 4 .

 4 Fig. 4. Observed and predicted fuzzy numbers for a given input x j

  .

	model = Identification (data, x min , x max )
	% choice of shift
	Compute r init and r fin
	If r fin < r init Then
	shift = x max	% spread globally decreasing
	Else
	shift = x min	% spread globally increasing
	EndIf
	% identification of A 0 and A 1

TABLE III

 III

					model			Sum	Dist
		Yu model [14]	20	14.5
		Our model (18)	15	7
	Y ˆx ( )	=	(	(	10 1 ,	) 0.5 0.5 , ( +	) x 11 -(	)	) 1 5 11 , [	]
		+	(	(	10 1 ,	) 0.5 1.33 , ( +	)x 11 -(	)	) 1 11 17 , [	]	DATA CONSIDERED AND IDENTIFIED MODEL
											interval	modal tendency	spread	A 0	A 1	shift
											[1,25]	20	random	(24.3,25.6)	(-0.1,0)	1
											[15,25]
											[26,50]	20	random	(23.9,9.2)	(-0.1,0.1)	26
											[5,10]
											[51,75]	1.5x + 1	3/2x	(27.6,4.9)	(1.4,1.5)	51
											[76,100]	1.5x + 1	-3/2x	(102.9,5.6)	(1.4,1.3)	100
											[101,125]	100	2x	(101.6,7.9)	(-0.1,2.1)	101
											[126,150]	100	-x	(103.5,3.9)	(0.1,1.9)	150
											[151,175]	-x	x	(104.3,7.8)	(-0.9,0.9)	151
											[176,200]	-x	-x	(51.4,7.8)	(-1.1,0.8)	200
											[201,225]	-x	1.5x	(51.8,9.4)	(-0.9,1.5)	201

____________________________________________________________________________________________________ Fig. 3. Identified model on a time serie data set


Y Zhang

G Cheng
email: gcheng@opt.ac.cn

G Huo

Y Wang

W Zhao

C Mauclair

R Stoian

R Hui

NOVEL METHODS OF LASER TECHNOLOGIES The Fabrication of Circular Cross Section Waveguide in Two Dimensions with a Dynamical Slit 1

Keywords: numbers: 42.82.Et, 42.81.Qb, 42.65.Tg, 42.55.Rz

A waveguide with a circular transverse profile can be obtained by using a slit beam shaping method. Applications about optical circuit chips in glass require a circular cross section waveguide in two dimensions.

Here we reported to realize a fabrication of circular cross section waveguide in two dimensions by rotating the slit corresponding to the tangent line of the arc of the waveguide. The dependence of cross section on width of the slit was investigated in experiment. The difference of slit width between geometrical optics pre diction and experimental result was explained by considering the focal shift in xy directions due to diffraction effects.

INTRODUCTION

Femtosecond laser optical waveguide writing has attracted a lot of interest as it continuously demon strates its capability to achieve complex three dimen sional photonic structures with very simple equipment in rapid prototyping process [1][2][3]. Longitudinal and transverse writing geometries are possible by adjusting the relative translation between the beam and the sam ple. In the longitudinal geometry, the waveguides are intrinsically symmetric, however the waveguide dimension is limited by the focal length of the focusing objective [4,5], and homogeneous longitudinal waveguide writing requires dynamical aberration com pensation [6]. The transverse geometry provides a much greater flexibility and allows one to write waveguides or photonic circuits of arbitrary length and complexity. However it has the disadvantage of pro ducing a strong asymmetry in the waveguide cross section. The waveguide cross section has a width equal to about twice the beam waist, while it has a dimension in depth of the order of the confocal parameter, which is typically much larger. This prob lem can be overcome by the use of suitable beam shap ing techniques [7,8]. Circle cross section waveguides in silica, polymer and borate glass has been success fully demonstrated in one dimensional writing scheme [9][10][11][12]. Nonlinear effect has been considered and achieved deep writing in fused silica with a mechanical slit [13]. Also a two dimensional deformable mirror to 1 The article is published in the original. shape the spatial profile of an ultrafast laser beam was demonstrated to inscribe circular cross section waveguide in soda lime silica glass slide [14]. In fact, astigmatic writing beam worse the waveguide cross section if the sample is translated parallel to the short axis of the astigmatic beam. So a cylindrical telescope or a fixed slit can't achieve two dimensional circle cross section waveguide which is normal required in planar light wave circuit (PLC) chip, such as direc tional coupler, wavelength division multiplexed (WDM), arrayed waveguide grating (AWG) [15]. In order to obtain circle cross section waveguide in two dimensions, a simple method is to rotate the slit corre sponding to the waveguide configuration. More par ticular, the long axis of the astigmatic writing beam keeps parallel to the tangent line of the waveguide arc. In this article, a mechanical slit fixed on a computer controlled 360° rotation stage is used as a dynamical slit to realize 2D circle cross section waveguide.

THEORY

Since the first report on realization circular cross section waveguide writing with astigmatic laser beam [7], it is the key consideration to understand how to calculate the aspect ratio of truncated Gaussian beam or elliptical Gaussian beam. Geometry optics has been employed based on a condition that the beam waist at one direction is set to equal to the confocal parameter [7,11], and an aspect ratio of about 3 has been obtained. However it seams to be underestimated, low repetition rate laser experimental results (no strong thermal effect) is bigger than the theoretic prediction.

Then diffraction method has been proposed [16]. The writing of circular cross section waveguide has been demonstrated in simulation and experiment with a 12 aspect ratio, which is bigger than most experimental results so far [17,18] if neglecting the material differ ence and nonlinear effects during writing a waveguide with femtosecond laser. In spite of this, diffraction theory is an accurate method to this problem, just need special algorithm to get high precise solution. Here we tried to improve geometrical models by using a trun cated Gaussian beam model. And find that geometri cal optics can't hold any more due to intense diffrac tion. A simple model is proposed to explain how dif fraction affects the focus position and confocal parameter.

In an attempt to improve the model which based on geometry optics, we compare the results under three kinds of beam parameters, plane wave in x and y direc tions, Gaussian distribution in x and y directions, and Gaussian beam in y direction and plane wave in x directions. The position of the slit is shown in Fig. 1. Beam parameters and focusing objective parameters are as following: the diameter of laser beam is 6 mm, the NA of focusing objective is 0.55, the effective focal length of objective is 4 mm, and the wavelength of the writing laser is 800 nm. Here we used 50% maximum as the radius of beam waist considering that the defini tion of Rayleigh length is the distance along the prop agation direction of a beam from the waist to the place where the area of the cross section is doubled, that is, intensity is half of the maximum. For plane wave the beam waist is 0.84λ/NA, and the beam waist of Gaus sian beam goes to ω Y = 1.09λ/NA when truncated ratio goes to infinity, and ω X = 1.03λ/NA [19]. The laser intensity distribution near the laser focus decides the cross section configuration of waveguide, the princi ple to realize a circular cross section waveguide is to keep the bigger beam waist between two beam waists at orthogonal directions (X direction beam waist in this paper) equal to the confocal length (twice Rayleigh length). Here Rayleigh length is derived by the smaller beam waist. Under those conditions, confocal param eter is 2.93 μm, NA x = 0.23, and the aspect ratio is 2.4, for plane wave; and confocal parameter is 4.94 μm, NA x = 0.17, and the aspect ratio is 3.3 for Gaussian beam, and confocal parameter is 4.94, NA x = 0.14, and the aspect ratio is 4.0 under the third condition. Our experimental results are the range of 7.5 to 10. So it shows that plane wave is the poorest and the condition of Gaussian beam in y direction and plane wave in X direction is the most close to experiment.

When considering diffraction effects, we assume there are two independent beams with different beam waists going through a same objective. Gaussian laser beam is focused by a lens and the point of the absolute maximum irradiance of the focused field is located on the axis but at a distance difference Δf, as called focal shift [20][21][22]. An explicit expression of the relative focal shift in focused, apertured Gaussian beam: (1) where

Δf G = -f/(1 + π 2 G 2 ), G = ω 2 /λf, f is focal length, λ is the wavelength of laser. Δf D = -f/{1 + N[1 + (π 2 N/12) 1.51 ] 1/1.51 }, N = α 2 /λf, α = N/G,
a is radius of circular aperture. N is Fresnel number of diffracting aperture, represents the number of Fresnel zones that fill the aperture when the aperture is viewed from the geometrical focus. In our up mentioned parameters, Fresnel number is 704 in y direction when laser beam is 3 mm (equal to the diameter of objective pupil), and focal shift is 8 nm. However in x direction Fresnel number is 12.5 when slit width is 0.4 mm, it contrib utes a focal shift of 30.3 μm which is bigger than con focal length (about 5 μm). Now it is clear that the foci in x direction and in y direction are not same position. Confocal parameter is extended extensively, can't cal culated by Geometry optics. This situation is shown schematically in Fig. 1a. The effective focus of astig matic beam should be located in the between Z X0 and Z Y0 . By use formula (1), we calculated the dependence of focal shift on slits width, shown in Fig. 1b. Focal shift increases exponentially as the width of slit decreases. It is easy to understand that smaller slit pro duces stronger diffraction. However it is difficult to calculate precise results based on Fraunhofer diffrac tion theory [23,24].

Δf Δf G 1 0 . 3α 2 - ( ) exp - [ ] = + Δf D 0.3α 2 /G - ( ) exp [ ] , 0

EXPERIMENT SETUP

The structures presented below were achieved in the transversal configuration, where the sample is translated perpendicular to the laser propagation axis. Polished fused silica parallelepipedic samples are irra diated with 120 fs pulses from an 1 kHz Ti:sapphire ultrafast regenerative amplifier laser systems (Spitfire, Spectra Physics) which delivers an average power of 800 mW at the center wavelength of 800 nm. A com puter driven electromechanical shutter, synchronized with the movements of the positioning system (Physik Instrumente M 405.DG, M 126.DG, M 111.PD) permits the writing of longitudinal structures in the bulk of the sample. A long working distance micro scope objective (Mitutoyo, NA = 0.55, WD = 13 mm) is employed to focus the femtosecond beam in the sil ica glass. A Zernike type positive optical phase con trast microscopy (PCM) system is employed to moni tor the laser irradiated areas. The slits with sizes of 6.0 × 0.5 mm, 6.0 × 0.4 mm, and 6.0 × 0.3 mm are used to generate an astigmatic writing beam. The slit is in the holder of a rotation stage (Physik Instrumente M 060.DG). There are two YZ stages to help optical alignment. One stage makes sure that the writing laser beam goes through the center of the rotation stage, and another one to adjust the position of the slit at the center of the rotation stage. The movement of the rotation stage is decided by the curve of the designed waveguide. Since the effect of spherical aberration is not negligible thus leading to an increase of aspect ratio [16]. The experiment was conducted at a focal depth of 200 μm. The schematic of experimental setup is shown in Fig. 2.

STRAIGHT WAVEGUIDE WRITING

In order to determine the optimum slit width for fabricating waveguides with circular cross sections, groups of channels were written at the same pulse energy but with slit widths ranging from 300 μm up to 500 μm. After irradiation, the sample was side pol ished and checked under an optical microscope in transmission mode. The aspect ratio near 1:1 obtained with 400 and 300 μm slit. Looking back the focal shift shown in Fig. 1, it is easy to understand why there is big error between geometrical optics prediction and experimental results since the focal shift is more bigger than confocal parameter.

In relatively gentle exposure conditions, two inten sity dependent regimes of positive refractive index modifications, hereafter referred to as type I and type II, were mentioned upon irradiation with 800 nm fem tosecond laser pulses in early paper [25,26]. The guid ing properties of type I and type II waveguides have been analyzed in our recent paper [27]. The fantastical nonlinear effects [5,28] such as supercontinuum gen eration in nanoporous doped glass in the structure induced by femtosecond laser is absent in our case. In transversal writing configuration, there provides type I and type II waveguides. And a schematic description of the processing window for each resulting waveguide type in a transversal writing configuration with the slit was given in Fig. 3, under following irradiation condi tions: 100 fs pulse duration, radiation wavelength 800 nm, laser repetition rate 1 kHz. The size of the slit is 6.0 × 0.4 mm, the NA of the focusing objective is 0.55NA, the pupil of the objective is 3 mm. The evalu ation is performed at a working depth centered around 200 μm. The reported pulse length values were mea sured after the focusing objective and were determined by adjusting the compressor position to compensate additional dispersion. Region I produces smooth and low loss type I waveguide, and region 2 is according to type II waveguide which is high optical loss and polar ization sensitive normally.

Figure 4 shows the phase contrast pictures and the cross sections of the waveguides written at different speeds. Black cross sections comes from the scattering or absorption of the waveguides, the resulting phase contrast pictures of the waveguides are white and pos sesses birefringence. When the scanning speeds is beyond 40 μm/s, the cross section is white in trans mission microscopy, which shows a good guiding property. Figure 5 shows the cross section of the waveguides written with 300 μm slit at different scan ning speed. The diameter of the type I waveguides almost keeps constant under different scanning speeds, is about 5.2 μm, this value is very close calcu lated one when taking half maximum full width (HMFW) of a Gaussian beam waist (in our case, it is 4.94). The guiding properties were checked by injec tion HeNe laser, it is a typical type I waveguide, and there is no polarization effect and birefringence.

CURVED WAVEGUIDE WRITING

Based on straight waveguide writing condition, we conducted 90° circular arc waveguide writing. As the rule of thumb, realization of circular cross section arc waveguide fabrication requires a rotating slit. The long axis of the slit keeps parallel to the tangent line of the waveguide arc. Circular arc simplifies this calculation requirement of tangent line of curve. During writing, just keep the angle of slit equal to that of the writing arc. The inset rectangles in Fig. 6a present the direc tion of the slit while writing an arc. Femtosecond laser inscribed waveguides feature very low index changes and always qualify as weakly guiding structures in fused silica. Typically, index changes is between 0.0001 and 0.001 in silica [4,29]. Considering this small refractive index changes leads to huge bend loss at small radius, we wrote an arc with a radius of 500 and 600 μm. A square fused silica sample is employed for convenient mode check of waveguide two cross sections at two ends. The laser power is 1.2 mW. Linear velocity is about 50 μm/s here. Figures 6b and6c are the cross sections of the two ends of the waveguide. When we check out the cross sections of the arc waveguides, the sample is rotated to vertical in compared with the writing condition. So there is a shadow in the cross section pictures. Figure 6d is a portion of the arc waveguide viewing in transmission microscopy.

CONCLUSIONS

We demonstrated how to realize a fabrication of circular cross section waveguide in two dimensions with a dynamical slit. It is helpful to planer light wave circuit (PLC) chip in glass. The dependence of cross section on width of the slit was investigated in experi ment, and processing window for type I waveguide was given under the conditions of our parameters. The dif ference of slit width between geometrical optics pre diction and experimental result was explained by con sidering the focal shift in xy directions due to diffrac tion effects.

Fig. 1 .

 1 Fig. 1. Focal positions of a focused Gaussian beam after a slit (a). Due to Fresnel numbers are different for two orthogral directions, focal shift at x direction is more big ger than y direction. As a results, the real focus of focused Gaussian beam after slit will locate between Z x0 and Z y0 . (b) Shows the dependence of focal shift on radius of aper ture.

Fig. 2 .

 2 Fig. 2. Experimental setup of the femtosecond laser waveguide writing arrangement indicating the irradiation geometry: HWP half wave plate, PS polarized splitter, OBJ microscope objectives, LED light emitting diode. The slit is fixed on a YZ stages which is in the holder of computer controlled rotation stage. Another YZ stages under the rotation stage is used to align laser beam going through the center of the rotation stage.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Processing windows for generating type I transverse waveguiding traces in fused silica with a 400 µm slit.

Fig. 5 .

 5 Fig. 5. The cross sections of the waveguides written by 800 nm femtosecond laser radiation at different scanning speeds. The width of the slit is 300 µm. the writing power is 2 mW after the slit. The scanning speeds in (a)-(d) are 20, 40, 80, and 160 µm/s, respectively.

Fig. 6 .

 6 Fig. 6. A 90° arc waveguide and its two end cross section modes, (a) is schematic view of the photoinscribed arc structure in fused silica using ultrashort pulsed laser radia tion, the inset rectangles shows the directions of the slit. (b, c) Cross sections of the two ends of the waveguide of the 600 µm arc. (d) A portion of the arc waveguide viewing in transmission microscopy. The radius of the arc are 500 and 600 µm in (d), respectively. The writing power is 1.2 mW.

ACKNOWLEDGMENTS

This work was supposed by West Light Foundation of The Chinese Academy of Sciences (no. 0729591213) and Innovative Research International Partnership Project of The Chinese Academy of Sci ences


Mathias Rousset
email: mathias.rousset@inria.fr

On a probabilistic interpretation of shape derivatives of Dirichlet groundstates with application to Fermion nodes

This paper considers Schrödinger operators, and presents a probabilistic interpretation of the variation (or shape derivative) of the Dirichlet groundstate energy when the associated domain is perturbed. This interpretation relies on the distribution on the boundary of a stopped random process with Feynman-Kac weights. Practical computations require in addition the explicit approximation of the normal derivative of the groundstate on the boundary. We then propose to use this formulation in the case of the so-called fixed node approximation of fermion groundstates, defined by the bottom eigenelements of the Schrödinger operator of a fermionic system with Dirichlet conditions

Introduction and results

Throughout this paper, we consider a Schrödinger operator in R d of the form:

H = - ∆ 2 + V, (1)
with a smooth potential V going to infinity at infinity, and acting on real valued functions generically denoted with the letter ψ ('wave functions'). Such functions ψ will be defined up to a real valued multiplicative constant (e.g. in eigenvalue and/or variational problems). We also consider a general family θ → Ω θ

of open smooth domains in R d depending sufficiently smoothly of a parameter θ ∈ R p . The boundary will be denoted ∂Ω θ . Gradients in the space R d will be denoted ∇, and gradients with respect to θ, ∇ θ . The Dirichlet groundstate and its Dirichlet groundstate energy (ψ * θ , E * θ ), are then defined as the unique bottom eigenelements of H, solution to the variational problem:

E * θ := inf     Ω θ ψH (ψ) Ω θ ψ 2 , ψ| ∂Ω θ = 0     = Ω θ ψ * θ H (ψ * θ ) Ω θ (ψ * θ ) 2 . ( 3
)
Calculus of variations detailed in Section 2 then yields the shape derivative of the groundstate energy through the formula:

∇ θ E * θ = - 1 2 ∂Ω θ |∇ψ * θ | 2 r θ dσ, (4)
where in the above σ is the usual surface measure induced by the canonical Euclidean structure R d , and r θ is the shape derivative, i.e. the field

r θ : ∂Ω θ → R p
such that formally the boundary variation writes down:

∂Ω θ+dθ = {x + n(x)r θ (x) • dθ | x ∈ ∂Ω θ } ,
where n(x) is the exterior normal vector at x ∈ ∂Ω θ . If {ψ θ } θ∈R p is a smooth family of smooth functions such that ψ θ (x) = 0 for x ∈ ∂Ω θ , then the shape derivative r θ can also be defined through

r θ (x)∇ψ θ (x) • n(x) = -∇ θ ψ θ (x) ∀x ∈ ∂Ω θ . (5)
Formula ( 5) can be proved as follows: from the Dirichlet conditions, one has for any small h: 4) is a classical result of abstract analytic perturbation of linear operators (see [START_REF] Kato | Perturbation theory for linear operators[END_REF]), but can be proved directly with the variational formulation as in [START_REF] Garcia Melian | On the perurbation of eigenvalues for the p-laplacian[END_REF].

0 = ψ * θ+h x + h • r θ (x)n(x) + O(|h| 2 ) -ψ * θ (x) = h • r θ (x)∇ψ * θ (x) • n(x) + h • ∇ θ ψ * θ (x) + O(h). Differentiability in formula (
In Section 3, we introduce a standard Wiener process (Brownian motion)

t → W t ,
with some given initial distribution in Ω θ . The first exit time of the domain Ω θ is denoted by

τ := inf (t ≥ 0|W t ∈ ∂Ω θ ) . (6)
The long time probability distribution of the latter process with Feynman-Kac weights, and conditioned to remain in the domain

Ω θ is denoted dη * θ Ω θ ϕ dη * θ := lim T →+∞ E ϕ(W T )1 1 T ≤τ e -T 0 V (Ws)ds E 1 1 T ≤τ e -T 0 V (Ws)ds groundstate ψ * θ .
Especially, the pointwise convergence of the normal derivative ∇ψ n θ .n may be hard to achieve. However, for practical situations in high dimension, we do not know any alternative point of view. A clear motivating example of a high dimensional problem is the case of Fermionic systems where the so-called fixed node approximation is used.

In Section 4, we introduce Fermionic groundstates (ψ * F , E * F ) associated to a finite symmetry group S ⊂ O(R d ) of the Hamiltonian H in (1); where O(R d ) denotes the group of isometries. S is simply the permutation goup of identical particles for physical systems. Fermionic groundstates are the solutions to the variational problem:

E * F := inf     R d ψH (ψ) R d ψ 2 , ∀S ∈ S, ψ • S = det(S) ψ     = Ω θ ψ * F H (ψ * F ) Ω θ (ψ * F ) 2 . ( 13
)
Any function ψ verifying the symmetry property ∀S ∈ S, ψ • S = det(S)ψ will be called skew-symmetric, whereas any function ψ verifying ∀S ∈ S, ψ • S = ψ will be called symmetric. Note that existence of (ψ * F , E * F ) follows in our context from the fact that H has a discrete spectrum, and is a classical result of spectral theory for more general potential V (see references in [START_REF] Cancès | Quantum monte-carlo simulations of fermions. a mathematical analysis of the fixed-node approximation[END_REF]); but uniqueness does not hold in general. In practice, (ψ * F , E * F ) is computed using a parametrization of skew-symmetric functions, and a numerical optimization procedure. This is the main problem of computational Quantum Chemistry, and forms a huge scientific field. We refer to [START_REF] Cancès | Méthodes mathématiques en chimie quantique: Une introduction[END_REF] for a mathematical introduction with a consequent bibliography. See also the following two typical papers [START_REF] Umrigar | Energy and Variance Optimization of Many-Body Wave Functions[END_REF][START_REF] Toulouse | Optimization of quantum Monte Carlo wave functions by energy minimization[END_REF] involving wave function optimization using a Monte-Carlo method. Monte Carlo methods in computational Quantum Chemistry are referred to as Quantum Monte Carlo (QMC) methods. For physical systems, the parametrization is given as a finite sum of Slater determinants multiplied by a strictly positive symmetric factor (called the Jastrow fac-Nodal Monte-Carlo. They will require variance reduction techniques exploiting the symmetry structure, in principle such that the variance of ∇ θ E FN θ scales appropriately to 0 in the limit where µ FN θ,λ becomes symmetric (in other words, we seek for an 'asymptotically scaling variance reduction', also called 'zero-variance estimation' in the QMC literature). Ideas are given for future work on this matter.

Let us now position the content of this paper in the context of the QMC literature. First, we recall that efficiently computing in high dimension d ≫ 1 the Fermionic groundstate (13) using Monte-Carlo methods is a fundamental problem with many applications; for instance it amounts to solve the eigenvalue problem for excited eigenstates, where classical power methods fail. A general solution is known to be intractable, and is usually referred to as the sign problem (see Remark 4.2). This explains the necessity of the fixed node approximation. The issue of optimizing the nodes of the trial wave function ψ I θ in the fixed node approximation was pointed out in [START_REF] Ceperley | Fermion nodes[END_REF], where a long discussion on the structure of Fermion nodes and appropriate (from this perspective) trial wave functions is provided. Yet state-of-the art numerical methods optimizing the trial wave function is based on either, (i) the minimization of the VMC energy E I θ in ( 16), as in [START_REF] Umrigar | Energy and Variance Optimization of Many-Body Wave Functions[END_REF][START_REF] Toulouse | Optimization of quantum Monte Carlo wave functions by energy minimization[END_REF]; or (ii) the minimization of the variance of the local energy

E L := V -(ψ I θ ) -1 ∆ 2 (ψ I θ ) (26)
as in [START_REF] Umrigar | Energy and Variance Optimization of Many-Body Wave Functions[END_REF][START_REF] Hongxin | An improved algorithm of fixed-node quantum monte carlo method with self-optimization process[END_REF]. As a consequence, an efficient method optimizing directly the fixed node energy E FN θ in [START_REF] Del Moral | Particle approximations of Lyapounov exponents connected to Schrödinger operators and Feynman-Kac semigroups[END_REF], that is to say the nodes of the trial wave function ψ I θ , remains an unsolved problem and motivates the content of this paper. However, methods to approximately compute the gradient ∇ θ E FN θ were already suggested in the QMC literature in the more general context of the calculation of physical properties (or "forces"). The main goal is to compute the derivative ∇ R E * F,R of the groundstate energy with respect to some parameter R parametrizing the Hamiltonian, typically the potential energy V ≡ V R . A classical example is the following: R is the vector of the nuclei-nuclei distances in molecules. The trial wave function now also depends on R: ψ I θ ≡ ψ I R,θ . Well established methods are available for the exact and variance reduced computation of the variational energy gradient:

∇ R E I θ,R (27)
where E I θ,R is defined by [START_REF] Doucet | Sequential Monte-Carlo Methods in Practice[END_REF]. For instance, methods using coupling (correlated sampling) to reduce variance were tackled in [START_REF] Filippi | Correlated sampling in quantum Monte Carlo: A route to forces[END_REF], and a general construction of variance reduced estimators ('zero variance-zero bias' estimators) was proposed in [START_REF] Assaraf | Zero-variance zero-bias principle for observables in quantum Monte Carlo: Application to forces[END_REF]; see also some applications in [START_REF] Toulouse | Zero-variance zero-bias quantum Monte Carlo estimators of the spherically and system-averaged pair density[END_REF]. The case of the fixed node energy gradient:

∇ R E FN θ,R (28)
where E FN θ,R is defined by ( 15) is more intricate, and requires some approximation since the derivatives ∇ R ψ FN θ,R of the fixed node groundstate remain unknown. Note that computing the gradients ∇ θ E I θ or ∇ θ E FN θ,R can be seen as a particular case of computing respectively [START_REF] Toulouse | Zero-variance zero-bias quantum Monte Carlo estimators of the spherically and system-averaged pair density[END_REF] or [START_REF] Toulouse | Optimization of quantum Monte Carlo wave functions by energy minimization[END_REF], since the nodes of the trial wave function ψ I R,θ depend on R a priori . An approximate formula has already been proposed to compute [START_REF] Toulouse | Optimization of quantum Monte Carlo wave functions by energy minimization[END_REF], see for instance equation (54) in [START_REF] Assaraf | Zero-variance zero-bias principle for observables in quantum Monte Carlo: Application to forces[END_REF], [START_REF] Ceperley | Fermion nodes[END_REF] in [START_REF] Casalegno | Computing accurate forces in quantum Monte Carlo using Pulay's corrections and energy minimization[END_REF], [START_REF] Ceperley | Fermion nodes[END_REF] in [START_REF] Badinski | Total forces in the diffusion Monte Carlo method with nonlocal pseudopotentials[END_REF], and [START_REF] Del Moral | Branching and Interacting Particle Systems approximations of Feynman-Kac formulae with applications to nonlinear filtering[END_REF] in [START_REF] Badinski | Nodal Pulay terms for accurate diffusion quantum Monte Carlo forces[END_REF]. The terms due to the variation of the nodes (which amounts to evaluate ∇ θ E FN θ,R ) are called 'nodal Pulay terms', and were particularly pointed out in [START_REF] Badinski | Total forces in the diffusion Monte Carlo method with nonlocal pseudopotentials[END_REF][START_REF] Badinski | Nodal Pulay terms for accurate diffusion quantum Monte Carlo forces[END_REF]. In the above references, the formula used for the calculation is exactly the formula [START_REF] Reed | Methods of modern mathematical physics. IV. Analysis of operators[END_REF]. However, the interpretation in terms of a stochastic process stopped on the nodes in ( 24)-( 23), the analysis of the symmetry of the associated distribution on the nodes, and the suggestion of an associated Nodal Monte Carlo method are new results.

The following classical textbooks are recommended:

• about spectral theory of operators: [START_REF] Reed | Methods of modern mathematical physics. IV. Analysis of operators[END_REF];

• about elliptic theory of Partial Differential Equations: [START_REF] Gildbarg | Elliptic Partial Differential Equation of Second Order[END_REF];

• about random processes and Feynman-Kac representations: [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF];

• about Monte-Carlo methods in Quantum Chemistry (QMC): [START_REF] Hammond | Monte Carlo Methods in Ab Initio Quantum Chemistry[END_REF][START_REF] Assaraf | A pedagogical introduction to Quantum Monte Carlo[END_REF].

Shape derivatives of Dirichlet groundstates

In this section, some notations and results are recalled concerning Dirichlet groundstates of Schrödinger operators with a smooth potential. Then, formula (4) is proven formally, and references are given for rigorous proofs.

Consider the Schrödinger operator (1) defined on R d with a smooth potential V bounded from below. H defines a self-adjoint operator on the Hilbert space L 2 (R d ). For simplicity, V is assumed to go to infinity at infinity such that H has a compact resolvent, and thus a purely discrete spectrum.

Generalization to operators involving continuous spectrum, although of fundamental importance in Quantum Chemistry, is left as technical extensions of the present work. Then a parametrization of smooth domains (2) is considered for θ ∈ R p such that ∂Ω θ has a smooth boundary for any θ. One assumes then that there exists a set of diffeomorphisms smoothly indexed by θ and such that:

               (θ, x) → R θ (x) is smooth in R p × R d R 0 = Id R θ = Id outside some compact set ∀θ ∈ R p , Ω θ = R θ (Ω 0 ). (29)
In this setting, the shape derivative of θ → Ω θ can be defined as the smooth field:

r θ : ∂Ω θ → R p , verifying: ∀x ∈ ∂Ω θ , r θ (x) = ∇ θ R θ (x) • n(x), (30)
and thus locally for small h ∈ R p :

∂Ω θ+h ∼ {x + h • r θ (x)n(x)| x ∈ ∂Ω θ } .
Now classical results of spectral theory ensures that the Hamiltonian (1) considered in L 2 (Ω θ ) with Dirichlet boundary condition, is self-adjoint with domain D(H) ⊂ H 1 0 (Ω θ ), where H 1 0 (Ω θ ) is the usual Sobolev space of function with Dirichlet conditions and square integrable first order derivatives. Moreover, H has a unique (up to a multiplicative constant) signed groundstate ψ * θ , solution of the variational problem (3), or equivalently solution of the eigenvalue problem:

       H(ψ * θ ) = E * θ ψ * θ ψ * θ | ∂Ω θ = 0 ψ * θ > 0,
where •| ∂Ω θ denotes the usual trace operator on the boundary. The regularity of V then ensures that ψ * θ is smooth on Ω θ and that ∇ψ * θ • n is smooth on ∂Ω θ . The following derivative formula can now be stated.

Lemma 2.1. Consider domains θ → Ω θ verifying [START_REF] Umrigar | Energy and Variance Optimization of Many-Body Wave Functions[END_REF]. Then the Dirichlet energy E * θ solution of the variational problem (3) is differentiable with respect to the parameter θ, and the variation formula (4) holds.

(H -E * θ 0 )(ψ * θ 0 ) = 0, we get (∇ θ E * θ ) θ=θ 0 = 1 2 ∂Ω θ 0 ∇ψ * θ 0 • n ∇ θ 0 ψ * θ 0 dσ Ω θ 0 ψ * θ 0 2 .
Then (5) applied to ψ * θ 0 yields the result. The rigorous proof can be made using a change a variable with the diffeomophism R θ , and then exploiting the variational formulation (see e.g. Theorem 2 in [START_REF] Costantini | Boundary sensitivities for diffusion processes in time dependent domains[END_REF]).

Probabilistic interpretations

In this Section, the probabilistic formulations ( 8)-( 9)-( 11) are proven and detailed. Associated Monte-Carlo methods are recalled with some references.

Consider notations and assumptions of Section 2. Let t → W t be a standard Wiener process with exit time τ from Ω θ defined in [START_REF] Cancès | Quantum monte-carlo simulations of fermions. a mathematical analysis of the fixed-node approximation[END_REF]. The classical probabilistic interpretation of the eigenelements (ψ * θ , E * θ ) is recalled in the following lemma: Lemma 3.1. Assume W 0 is distributed according to ψ init (x)dx ψ init (x)dx where ψ init ∈ L 2 (Ω θ ), and ψ init is non vanishing in each connected component of Ω θ . The groundstate elements can be expressed through the long time behavior of the process with Feynman-Kac weights and conditioned by large exit times. This yields the formulas (7)-( 8)- [START_REF] Ceperley | Monte-Carlo simulation of a many-fermion study[END_REF].

Proof. The result follows from the classical representation of parabolic equations through the Feynman-Kac formula (see [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF]). Let us recall the different steps of the argument. First, consider ϕ

∈ C ∞ c (Ω θ ) a smooth solution (t, x) → u t (x) ∈ C ∞ (R + × Ω θ ) of the parabolic problem:        ∂ t u t = -H(u t ) u t | ∂Ω θ = 0 u 0 = ϕ
Then using Itô calculus, it yields for any φ ∈ C ∞ (R + × Ω θ ):

d φ T -t (W t )e -t 0 V (Ws)ds = ∆ 2 -V -∂ t (φ T -t ) (W t )e -t 0 V (Ws) ds dt + e -t 0 V (Ws) ds ∇φ T -t (W t ) • dW t , (31)
so that applying the latter computation to φ = u between time 0 and the stopping time T ∧ τ = inf(T, τ ) yields:

Ω θ u T (x)ψ init (x)dx Ω θ ψ init (x)dx = E u T -T ∧τ (W T ∧τ )e -T ∧τ 0 V (Ws)ds = E ϕ(W T )1 1 τ ≥T e -T 0 V (Ws)ds .
Now denoting by (ψ * ,n θ , E * ,n θ ) n≥0 the full spectrum of H normalized in L 2 (Ω θ ), the above Feynman-Kac representation reads:

E ϕ(W T )1 1 τ ≥T e -T 0 V (Ws)ds = n≥0 e -E * ,n θ T Ω θ ψ * ,n θ ψ init Ω θ ψ init Ω θ ψ * ,n θ ϕ. ( 32
)
Now using the fact that the groundstate has a spectral gap (E * ,1 θ > E * ,0 θ = E * θ ), the dominant term in (32) when T → +∞ enables to verify that the groundstate has a sign ψ * ,0 θ > 0. Finally, taking the limit ϕ → 1 1 by dominated convergence in the formula above, and then the leading term in the limit T → +∞ yields ( 9). ( 8) follows.

In practice however, a diffusion solution to a stochastic differential equation with a repulsive drift at the boundary ∂Ω is used:

dX t = (ψ I ) -1 ∇ψ I (X t )dt + dW t ψ I | ∂Ω θ = 0, ( 33
)
where ψ I > 0 is a smooth function strictly positive in Ω θ and vanishing on ∂Ω θ . In [START_REF] Cancès | Quantum monte-carlo simulations of fermions. a mathematical analysis of the fixed-node approximation[END_REF], sufficient conditions on the behavior of ψ I near ∂Ω θ and at infinity are given for (33) to be well-posed, and to verify the following variant of ( 8)-( 9):

lim T →+∞ E ϕ(X T )e -T 0 E L (Xs)ds E e -T 0 E L (Xs)ds = Ω θ ϕ ψ * θ ψ I Ω θ ψ * θ ψ I , ( 34
) lim T →+∞ - 1 T ln E e -T 0 E L (Xs)ds = E * θ , (35)
where the so-called local energy is defined by [START_REF] Rousset | On the control of an interacting particle estimation of Schrödinger ground states[END_REF]. The proof of the latter probabilistic interpretation (34)-( 35) is similar to the proof of Lemma 3.1, and is based on the mapping of the Hamiltonian H to a weighted L 2 space through:

H I = (ψ I ) -1 H(ψ I •) = - ∆ 2 -(ψ I ) -1 ∇ψ I • ∇ + E L .
Details and assumptions can be found in [START_REF] Cancès | Quantum monte-carlo simulations of fermions. a mathematical analysis of the fixed-node approximation[END_REF].

Next the probabilistic interpretation of the shape derivative given by formula ( 11) is proven. Proposition 3.2. Assume W 0 is distributed according to ψ init (x)dx ψ init (x)dx where ψ init ∈ L 2 (Ω θ ), and ψ init is non vanishing in each connected component of Ω θ . Assume the boundary ∂Ω θ is smooth and uniformly Lipschitz, and consider the measure defined for any λ < E * θ by [START_REF] Ceperley | Fermion nodes[END_REF]. Then (11) holds.

Proof.

Step 1. Let ϕ ∈ C ∞ (Ω θ ). We claim that the parabolic differential equation with inhomogenous Dirichlet conditions

       ∂ t h t (ϕ) = -(H -λ) (h t (ϕ)) h t (ϕ)| ∂Ω θ = ϕ| ∂Ω θ h 0 (ϕ) = ϕ (36)
has a unique smooth solution for λ < E * converging exponentially fast towards h ∞ (ϕ) unique smooth solution of the elliptic inhomogenous Dirichlet problem:

(H -λ) (h ∞ (ϕ)) = 0 h ∞ (ϕ)| ∂Ω θ = ϕ. (37)
This is a classical consequence of spectral theory, but we recall briefly the basic arguments. Existence of a smooth solution of (37) follows from the fact that (Hλ) -1 can be extended to a bounded operator of L 2 (Ω θ ) so that:

h ∞ (ϕ) = (H -λ) -1 (H -λ) (ϕ) -ϕ,
which indeed is solution of (37). Note that in the above (Hλ) -1 implicitly refer to the operator with homogenous Dirichlet boundary condition, so that (H -λ) -1 •(H -λ) = Id when operating on test function with inhomogenous boundary conditions. Then the homogenous solution of (36) has been solved using spectral decomposition in (32), proving that such homogenous solution is unique and vanishes exponentially fast when λ < E * θ . The latter analysis of the homogenous case proves uniqueness in (36)-(37), as well as exponential long time convergence of the time dependent equation (36).

Step 2. We claim that for any ϕ ∈ C ∞ (Ω θ ),

Ω θ h ∞ (ϕ)ψ init Ω θ ψ init := E ϕ(W τ )1 1 τ <+∞ e -τ 0 (V (Ws)-λ) ds (38)
where h ∞ (ϕ) is the solution to the elliptic partial differential equation with inhomogenous Dirichlet conditions (37). Indeed, using (31) with test function h T -t (ϕ), and stopped at time T ∧ τ = inf(T, τ ) yields:

Ω θ h T (ϕ)ψ init Ω θ ψ init := E h T -τ ∧T (ϕ)(W τ ∧T )e -τ ∧T 0 (V (Ws)-λ) ds .
Now the event 1 1 τ =+∞ has null probability if t → W t is recurrent, and if the latter is transcient then using lim

∞ V = +∞: lim T →+∞ E 1 1 τ =+∞ e -T ∧τ 0 (V (Ws)-λ) ds = 0;
so that taking the limit T → +∞ leads to (38).

Step 3. We claim that

Ω θ h ∞ (ϕ)ψ * θ = - 1 2(E * θ -λ) ∂Ω θ ϕ∇ψ * θ • n dσ.
Since ψ * θ is the groundstate:

ψ * θ = 1 E * θ -λ (H -λ) (ψ * θ ),
so that integration by parts yields:

Ω θ h ∞ (ϕ)ψ * θ = - 1 E * θ -λ ∂Ω θ 1 2 ϕ∇ψ * θ • n dσ + 1 E * θ -λ Ω θ 1 2 ∇ψ * θ • ∇h ∞ (ϕ) + (V -λ)h ∞ (ϕ)ψ * θ = - 1 E * θ -λ ∂Ω θ 1 2 ϕ∇ψ * θ • n dσ + 0.
Step 4. Now using the Markov property of Brownian motion yields:

E ϕ(W τ )1 1 T ≤τ <+∞ e -τ 0 (V (Ws)-λ) ds = E E ϕ( W τ )1 1 τ <+∞ e -τ 0 (V ( Ws)-λ) ds | W 0 = W T 1 1 T ≤τ e -T 0 (V (Ws)-λ) ds ,
where (( W t ) t≥0 , τ ) is an independent copy of (W t ) t≥0 , τ ). Together with Step 2 and Step 3 with the probabilistic interpretation of the groundstate ψ * θ in ( 10)-( 8), it completes the proof.

Fermion groundstates

Consider the Schrödinger operator [START_REF] Assaraf | A pedagogical introduction to Quantum Monte Carlo[END_REF]. Fermionic groundstates (ψ * F , E * F ) associated to (1) are defined with respect to a finite symmetry group S ⊂ O(R d ) of the potential V , where O(R d ) denotes the group of isometries of R d . A symmetry group of V is defined by the property

∀S ∈ S V • S = V.
Fermion systems appears in the special case when S contains symmetries with odd parity: det(S) = -1, so that H can be defined (since the Laplacian operator commutes with isometries) to operate on the Hilbert space of skew-symmetric function:

H = ψ ∈ L 2 (R d ) | ∀S ∈ S, ψ • S = det(S)ψ .
"Fermion" groundstates are then the solutions to the variational problem [START_REF] Moral | Feynman-Kac Formulae[END_REF].

Example 4.1 (Fermions). In the case of N physical quantum particles of Fermionic type in dimension 3, one has R d = R 3N , and a potential of the form:

V (x) := N i=1 V ext (x i ) + 1≤i<j≤N V int (x i -x j ),
where x = (x 1 , . . . , x N ) are the 3 dimensional coordinates of the particles, V ext is an exterior smooth potential that goes to infinity at infinity, and V int a smooth interaction potential that vanishes at infinity. Then the discrete symmetry group is the permutation group S N of physical particles. See [START_REF] Cancès | Quantum monte-carlo simulations of fermions. a mathematical analysis of the fixed-node approximation[END_REF] for the mathematical analyis of Quantum Monte-Carlo (QMC) methods in this context.

Since in this paper, we restrict to smooth operators with compact resolvent, the existence of (ψ * F , E * F ) follows directly from the fact that the spectrum is discrete, and smoothness from the fact ψ * F satisfies an eigenvalue problem of a smooth elliptic operator (see [START_REF] Gildbarg | Elliptic Partial Differential Equation of Second Order[END_REF]). F , E * F ) using Monte-Carlo methods is an untractable problem known as the sign problem. The latter can be summarized as follows. Since H leaves invariant H, eigenfunctions of H in R d are either symmetric or skew-symmetric functions. Thus computing (32) with:

• a skew-symmetric test function ψ I ,

• a non-symmetric positive initial condition ψ init > 0, yields:

E ψ I (W T )e -T 0 V (Ws)ds = ψ * ,n ∈H e -E * ,n T R d ψ * ,n ψ init R d ψ init R d ψ * ,n ψ I ,
and in principle if ϕ is symmetric:

lim T →+∞ E ϕ(W T )ψ I (W T )e -T 0 V (Ws)ds E ψ I (W T )e -T 0 V (Ws)ds = R d ψ * F ψ I ϕ R d ψ * F ψ I .
Unfortunately, this computation relies crucially on the rate of vanishing of the normalisation which is due to the non-symmetry of the initial condition ψ init > 0. Since stochastic processes used to compute such averages quickly forget their initial condition and have a symmetric dynamics in the full state space R d , one is compelled to compute ratios of vanishing averages (of the type 0 0 ), with Monte-Carlo estimators having a non-vanishing statistical variance. This leads to infinite variance when T is large. This forms the sign problem. This problem appears more generally when trying to solve higher eigenvalue problems with Monte-Carlo methods. Although, there is probably no general solution, solving the sign problem for particular situations in high dimension would be considered as a major breakthrough.

In practice, (ψ * F , E * F ) can be approximated using a hybrid methodology in two steps. First, a trial wave function is obtained using an analytical parametrization of H, usually of the form:

ψ I α,θ := J α ψ skew θ . (39)
In the above, J α > 0 is a strictly positive symmetric part called the Jastrow factor and parametrized by α ∈ R n ; and ψ skew θ is the skew-symmetric part parametrized by θ ∈ R p . The most convenient numerical optimization method is then used to solve [START_REF] Moral | Feynman-Kac Formulae[END_REF] in the space of parameters formed by (α, θ), and the solution of the optimization procedure is denoted with the parameters (α 0 , θ 0 ). This yields the set of functions ψ I θ θ∈R p defining the nodal domains through ψ I θ := J α 0 ψ skew θ defined in [START_REF] Gildbarg | Elliptic Partial Differential Equation of Second Order[END_REF]. Use for instance a long time trajectory of a drifted stochastic process of the type (33) with Feynman-Kac weights as in formula (34).

3. Use the latter sample as initial conditions; and sample independent Brownian motions. Stop them when they hit the nodes ∂N θ . Then compute the measure µ FN θ,λ according to (20).

4. Estimate the variations of the fixed node energy ∇ θ E FN θ with [START_REF] Kato | Perturbation theory for linear operators[END_REF].

As a conclusion, we mention two possible strategies of variance reduction that will be necessary for efficient computations. Their development is left for future work.

• First, a drifted diffusion instead of plain Brownian motion may be used in [START_REF] Gildbarg | Elliptic Partial Differential Equation of Second Order[END_REF] to reduce the variance caused by the exponential weights. However, the classical stochastic differential equation (33) where ψ I vanishes at the nodes cannot be used. Indeed, the repulsive drift at the nodes prevent the process from hitting the latter, and ( 20

Remark 4 . 2 (

 42 The sign problem). Computing directly (ψ *

  ) no longer holds. Instead of a skew-symmetric guiding function ψ I , we propose to use as a drift a symmetric function ψ B , for instance an approximation of the Bosonic groundstate ψ * B > 0 of H in R d solution of the eigenvalue problem H(ψ * B ) = E * B ψ * B .

Acknowledgements

The author would like to thank E. Cancès, B. Jourdain and T. Lelièvre for pointing out the problem and fruitfull discussions, as well as the referees for very carefull reading and useful suggestions.

∇ + ψ • n + denotes the exterior normal derivative of ψ in N + θ , and ∇ -ψ • n - the exterior normal derivative of ψ in N - θ . Then the skew-symmetrization of the normal derivative is defined by

and the plain symmetrization with respect to n + is defined by [START_REF] Hongxin | An improved algorithm of fixed-node quantum monte carlo method with self-optimization process[END_REF].

Note that ψ ∈ C 1 (R d ) if and only if ∇ sk ψ • n = 0 on ∂N θ , so that ∇ sk ψ • n can be seen as the gradient discontinuity of ψ at ∂N θ . One can now precise the idea of symmetry breaking on ∂N θ .

, one has the integration by parts formula:

Moreover, ∇ sk ψ • n = 0 on ∂N θ if and only if

In the opposite case, ∇ + ψ • n + and ∇ -ψ • n -are invariant only under the action of the special sub-group S + and a "symmetry breaking" occurs.

Proof. By skew-symmetry of ψ, one has for any S ∈ S in N + θ ∪ N - θ :

S ∇ψ • S = det(S) ∇ψ, and then on ∂N θ , S n + • S = det(S) n + S n -• S = det(S) n -, so that since S T S = Id, it yields on ∂N θ for any S ∈ S :

This yields the symmetry properties of ∇ψ sk •n and ∇ψ sy •n + with respect to S, and of ∇ + ψ • n + and ∇ -ψ • n -with respect to S + .