A Taky

G Belaud

Using a furrow system for surface drainage under unsteady rain

Keywords: Runoff, furrow system, unsteady rainfall, overland flow, kinematic wave model

Water excess during winter limits crop development on heavy clay soil conditions of the Gharb valley (Morocco). The furrow system to eliminate these negative effects is the adopted solution. This article focuses on the development of a water transfer model through a furrow system during unsteady rainfall event to evaluate the runoff volume resulting from a reference rainy event. This model contains a production function associated to a transfer function. The production function is based on the Green-Ampt infiltration equation. The latter has been adapted to account for unsteady rain conditions and rainfall intermittence. The transfer function is based on the kinematic wave model, the explicit solution of which is coupled with the water excess generated by the production function. Simulated runoff in the furrows is collected by a drainage ditch evacuating the flow outside a plot of 1.3 ha. The similarity between parameters of a furrow irrigation model and those of the production function is advantageously used for model calibration.

The proposed modelling approach shows capabilities to predict water amount and peak discharges evacuated from a plot of around 1 ha by a furrow system under unsteady rainfall events. As an application, it is used to evaluate the ability of the surface drainage system to evacuate the excessive volumes of water under typical rainfalls.

I f : cumulative infiltration during furrow irrigation (m) I rp : cumulated rainfall (m) J: the friction slope (m/m) k : rainfall event number K s : saturated conductivity (m/s) L f : length of the infiltration front (m) m : wall slope of the furrows (-) n : roughness coefficient (s/m 1/3 ) N: number of furrows P: wetted perimeter (m) P 1 : cumulative rainfall for the one year return period (mm) q : infiltration rate (m/s) q c : infiltration capacity (m/s) q r : rainfall intensity (m/s) Q: flow rate in the furrow during runoff (m 3 /s) Q': derivative Q with respect to Y (Q'= dQ/dY) Q p : peak discharge (m 3 /s) R : runoff ( m/s) r = R fs : lateral contribution due to runoff on a width fs (m 2 /s), R h : hydraulic radius (m) S 0 : field slope (m/m) t: current time (s) T a (x) : advance time of the water front during an irrigation event (s) t B : time of rain event beginning (s) t B k : time when the k-th rainfall event begins (s) t E k : time when the k-th rainfall event stops (s) t p : ponding time (s) x: furrow abscissa (m) Y: water depth (m): Z r : a reference soil depth (1 m)  : parameter of the K(h) equation (m -1 )  : parameter of the soil water depletion function  c : capillary length (m)  : parameter ruling the macro pore effect in the irrigation model : saturation deficit  t: calculation time step (s) x: calculation space step (m)  : water content (cm 3 /cm 3 ) θ B k : water content (cm 3 /cm 3 ) at the beginning of the rainfall event k θ E k : water content (cm 3 /cm 3 ) at the end of the rainfall event k

 s : saturated water content (cm 3 /cm 3 )  fc : water content at field capacity (cm 3 /cm 3 )  i : initial water content (cm 3 /cm 3 ) a shape parameter of the furrow= (1+m²) 1/2  : current time between -D and + D in the one year return period rainfall formulation (s)

Note: Numerical values of variables may be given in different units in the text and figures.

Introduction

Rainfall excess on the heavy clayey soils of the Gharb valley (Morocco) strongly limits the production of crops sown during winter. The subsurface drainage system, the main role of which is to maintain the level of the saline groundwater deep enough, does not allow a complete and fast evacuation of this water excess [START_REF] Bouarfa | Water table shapes and drain flow rates in shallow drainage systems[END_REF], due to the absence or degradation of land levelling. Considering energy costs and maintenance problems, sprinkler irrigation has often been replaced by traditional surface irrigation [START_REF] Mailhol | A furrow irrigation model to improve irrigation practices in the Gharb valley of Moroco[END_REF]. Impacts of water excess during winter on crop growing in unlevelled plots are clearly shown in [START_REF] Taky | Diagnostic des pratiques d'irrigation gravitaire et possibilités d'amélioration dans le Gharb au Maroc[END_REF]. This work also highlights the beneficial role of a furrow system on crop production (sugar beet) for both experimental and farmer plots. To encourage farmers adopting modern surface irrigation systems, furrow irrigation tests were conducted on 180 to 240 m long laser levelled plots, and water was supplied using siphons, floppy pipes or rigid gated pipes The furrow system, designed for water supply during the irrigation season, also allows the evacuation of excessive rainfall amounts during winter.

The objective of this article is to present a numerical model allowing a runoff prediction within a furrow system submitted to complex rainfall events. Overland flow modelling under constant rain intensities has been addressed by several researchers using a two dimensional approach [START_REF] Chow | Hydrodynamic modelling of two dimensional water flow[END_REF][START_REF] Zhang | Modelling of two dimensional overland flow[END_REF][START_REF] Taylor | Applicability of StVenant equations for two dimensional overland flow over rough infiltration surfaces[END_REF]. More recently [START_REF] Esteves | Overland flow and infiltration for small plots during unsteady rain: numerical results versus observed values[END_REF] proposed a comparable approach for complex rainfall events. Their overland flow model results from a coupling of the Green-Ampt infiltration solution with the Saint-Venant solution used for the surface transport.

The previous studies all refer to small-scale plots of a few square-meters. At this scale, runoff is subject to micro topography and a 2D-modelling approach is required for the surface water transport. [START_REF] Muños-Carpena | Modelling hydrology and sediment transport in vegetative filter strips[END_REF] also used the Green-Ampt model coupled with the kinematic wave model for sediment transport simulation in vegetative filter strips. [START_REF] Paige | A strip model approach to parameterize a coupled Green -Ampt kinematic wave model[END_REF] used a similar modelling approach with a rainfall simulator and analyzed the sensitivity of the runoff response to the rainfall rate. These two latter researches refer to a plane surface and the kinematic wave solution is not described.

Although soil properties of the experimental plot are assumed to be spatially constant, the originality of the present work is to propose a runoff transfer model through a furrow system model applicable at a plot scale of around 1 hectare. In addition, the production function is based on an infiltration model that uses the same parameters as the ones used to simulate the furrow irrigation. This runoff model is a useful tool to design the ditches for the runoff evacuation of a whole system of furrow irrigated plots. It could be considered as a component of a future integrated modelling approach for simulating furrow irrigation practices at the plot and season scale in the Gharb valley.

Model development

The hydrological model consists of a production function generating a water excess at the scale of an elementary sub-plot (Fig. 1) coupled with an overland flow model transferring the water excess through a furrow system. The simulated plot is composed of N furrows of length FL. These furrows outflow in a ditch at the extremity of which was installed an automatic gauge station (Fig. 2).

Production function

The production function calculates the water excess at the scale of the elementary sub-plot of width fs = 1.5 m (furrow spacing) and of length dx, the transfer space step. This water excess

Author-produced version of the article published in Agricultural Water Management, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014 is assumed to instantaneously supply a furrow element. Since water heights within the furrows are generally low, infiltration along their walls can be neglected. Consequently, infiltration will be globally considered as being governed by a 1D process. This production function needs to manage intermittent rainfall. A rainfall episode, the duration of which can be of one or two days, constitutes continuous rainfall events at different rates separated by no rain periods. For an unsteady rain event, there may exist different periods during which rainfall rate is greater than infiltration capacity. Under such circumstances, several ponding times can be generated while there is at most one ponding time under steady rain conditions. This production function uses the Green-Ampt infiltration model [START_REF] Green | Studies in soil physics I. The flow of air and water through soil[END_REF] as a simplified approach of the infiltration process. It assumes a homogenous soil profile and an initial water content  i uniformly distributed with depth. According to this model, the advancing front in the soil separates a saturated zone and a zone still at initial soil moisture. [START_REF] Mein | Modeling infiltration during a steady rain[END_REF] and [START_REF] Swartzendruber | Infiltration of constant flux-rainfall into soil as analyzed by the approach of Green and Ampt[END_REF] used this model to describe ponded infiltration with constant intensity rainfall conditions. This model was also applied under unsteady conditions [START_REF] Chu | Infiltration during unsteady rain[END_REF][START_REF] Esteves | Overland flow and infiltration for small plots during unsteady rain: numerical results versus observed values[END_REF]. [START_REF] Muños-Carpena | Modelling hydrology and sediment transport in vegetative filter strips[END_REF] also used the Green-Ampt approach with the method of [START_REF] Chu | Infiltration during unsteady rain[END_REF] to account for the water amount infiltrated before ponding. [START_REF] Paige | A strip model approach to parameterize a coupled Green -Ampt kinematic wave model[END_REF] used a same modelling approach with a rainfall simulator and performed a sensitivity analysis of runoff response to the rainfall rate. They showed that runoff is sensitive to rainfall intensity and initial moisture content. They also showed that a unique hydraulic conductivity parameter is not sufficient to correctly simulate the runoff peak on a surface with bare soil and vegetation especially in case of high rainfall intensities. A significant improvement of the simulations was obtained on the basis of a strip model configuration after tuning a second effective conductivity.

Runoff production during a rainfall event

Infiltration is composed of two stages. As described by [START_REF] Mein | Modeling infiltration during a steady rain[END_REF], when the rain intensity q r is lower than infiltration capacity q c , all the rain infiltrates into the soil and actual infiltration rate q is equal to q r . Cumulative infiltration I is equal to cumulative rain I rp . This is valid as long as the rain intensity is lower than the infiltration capacity q c . To calculate the infiltration capacity, the Green-Ampt approach is used, where the ponding time t p is the time at which q = q r = q c . According to Green-Ampt's model, the depth of the saturated zone (in which  =  s ) is given by

L f = I rp / (1)
where L f is the position of the wetting front,  s - i is saturation deficit (saturation moisture minus initial moisture) and I rp is the cumulated rainfall from the rain event beginning (t B ) until the ponding time. The application of the Darcy's equation gives the infiltration capacity for ponded conditions:

f f s s c L h K K q   (2)
where K s is saturated conductivity and h f is the capillary suction at the wetting front. After ponding (t > t p ), infiltration q equals infiltration capacity q c and the cumulative infiltration is calculated as follows:

   dt t q I I c t t rp p ) ( . ( 3
)
The infiltration capacity is calculated according to Eq. ( 2) in which the wetting front position is given by L f = I/ (4) A similar approach was used by [START_REF] Chu | Infiltration during unsteady rain[END_REF] who integrates Eq. ( 3) and gets a relation between I and t. A numerical method is required to determine I then its derivative q. Focussing on the infiltration rate q (rather than the cumulative infiltration), it is more convenient to solve Eqs (2) and (4) at each calculation time step, which is done using an explicit numerical scheme. Finally, runoff R(t) is obtained by the difference between rain intensity and infiltration capacity: R(t) = q r -q c (5) The water excess at the soil surface of the elementary sub-plot is assumed to be quasi instantaneously evacuated by the furrow system and consequently the resulting positive pressure head at the soil surface is not considered. Saturated conductivity K s is assumed constant and not affected by the rainfall rate, the air entrapment phenomenon not being considered here.

Soil parameters

A physical signification of h f has been proposed by some authors. [START_REF] Newman | Wetting front pressure head in the infiltration model of Green and Ampt Water[END_REF], for instance, has shown that

  i h s f dh h K K h 0 ) ( 1 (6)
,where h i is initial capillary suction.

In the frame of the validity of the Green-Ampt model [the delta Soil as defined by [START_REF] Philip | The theory of infiltration I. Sorptivity and algebraic infiltration equation[END_REF][START_REF] Philip | On solving the unsaturated flow equation. 1. The flux-concentration dependent relation[END_REF]] h f ≈ - c (capillary length) using the Gardner function, where saturated conductivity is assumed to increase exponentially with h [START_REF] Gardner | Some steady state solutions of the unsaturated moisture flow equation[END_REF]:

K(h) = K s exp(h) (7)
The value of  -1 =  c is the macroscopic capillary length scale [START_REF] Philip | Steady infiltration from circular cylindrical cavities[END_REF][START_REF] White | Macroscopic and microscopic capillarity length and timescale from field infiltration[END_REF][START_REF] Thony | Field measurement of the hydraulic properties of soil[END_REF] which is a hydraulic characteristic of the soil. From a physical point of view  c represents the average height of the capillary fringe. It is numerically equivalent to the soil water pressure head |h f | at the wetting front [START_REF] Kutilek | Agricultural Water Management[END_REF]. It varies from 3 to 5 cm for sand to 25 cm [START_REF] White | Macroscopic and microscopic capillarity length and timescale from field infiltration[END_REF] and more for clay [START_REF] Kutilek | Agricultural Water Management[END_REF][START_REF] Mailhol | A predictive form of Horton's equation for simulating furrow irrigation[END_REF].

A new procedure to account for rainfall intermittence

When the rainfall stops, the soil water depletion phase starts. Its duration depends on the occurrence date of the next rainfall event. The soil moisture conditions for predicting the infiltration capacity of the subsequent rainfall event need to be updated. Since the rain duration is quite short, generally less than one day, evapotranspiration is neglected. Denoting t B k the time when the k-th rainfall event begins, t E k when it stops; the corresponding soil moistures are denoted θ B k and θ E k . At the end of the rainfall event k, averaged soil water content on a reference depth Z r is given by:

 E k =  B k + I k /Z r ( 8
) where I k = cumulative infiltration (m) during the rainfall event k and Z r is the reference soil depth (m). A reference soil depth Z r = 1 m is assumed to be a valuable depth to account for Hortonian infiltration processes under high clayey soil conditions. At the beginning of the depletion stage (t

≥ t E k ), soil water content decreases rapidly from t E k ( E k ≤  s ) until it
reaches the field capacity, then the depletion is much slower [START_REF] Marcesse | Détermination de la capacité de rétention d'un sol au moyen de l'humidimètre à neutron[END_REF]. The soil water content at field capacity is denoted  fc . Because runoff is specific to the winter period, where

Author-produced version of the article published in Agricultural Water Management, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014

 is usually greater than  fc , the first depletion phase will be only considered. The depletion process is simulated as follows:

θ(t) =  fc + ( E k - fc ) exp[-t-t E k )]
, (9) in which β is an empirical coefficient depending on soil properties and soil conditions. Its value reflects the velocity at which the soil drains. At the beginning of the next rainfall event

k + 1, the soil water content is equal to  B k+1 =  fc + ( E k - fc ) exp [-t B k+1 -t E k )].
From this moisture, the infiltration capacity at the beginning of rainfall event (k+1) is calculated according to Eqs (1) and (2).

Transfer model

Model equations

The slope of the furrows is generally greater than 0.1%. In this range, the kinematic wave model is a good approximation of the Saint-Venant's transfer equations [START_REF] Woodhiser | Unsteady, one-dimensional flow over a plane: the rising hydrograph[END_REF] for the runoff generated by each elementary sub-plot. This model is often used to simulate the advancing front in furrow irrigation [START_REF] Walker | Kinematic wave furrow irrigation model[END_REF] but with different boundary conditions.

The kinematic wave solution combines the following set of equations:

r x Q t A       (10a) 3 / 2 2 2 2 0 h R A Q n J S   (10b)
where A is the wetted area in the furrow (m 2 ), t the current time (s), Q the discharge in the furrow (m 3 /s), r = R fs the lateral contribution due to runoff on a width fs (m 2 /s), S 0 the furrow slope (m/m), J the friction slope, n the Manning roughness coefficient (s/m 1/3 ) and R h is the hydraulic radius (m) defined as A/P (P = wetted perimeter). The boundary conditions are:

At t = 0, A(x,t) = 0 for all x (11a) At the upstream end (x = 0), A = 0 for all t.

(11b) Each furrow supplies the ditch which conveys the whole flow out of the plot. Due to the depth of the ditch, a free drainage boundary condition is set at the furrow downstream end. The kinematic wave solution of Eqs (10a, b) is also used for simulating the flow in the ditch, in which lateral contribution is the outlet discharge of each furrow. This set of equations is solved numerically according to an explicit scheme described below.

Numerical solution

The furrow is discretized with a constant length step dx. The simulation period is discretized with a time step dt. The discretized variables A and Q at the i-th space step and j-th time step are denoted A i j and Q i j . The terms of Eq (10a) are approximated at the first order as follows:

)] ( ) [( 2 1 1 1 1 1 n i n i n i n i A A A A dt t A           (12a) )] ( ) [( 2 1 1 1 1 1           n i n i n i n i Q Q Q Q dx x Q (12b) ) ( 2 1 1 n n r r r    (12c)
Author-produced version of the article published in Agricultural Water Management, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014

At the (n+1)-th time step, the values at the (i+1)-th abscissa are searched as functions of the values in i and those of the n-th time step. Two unknown terms need to be determined,

1 1   n i A and 1 1   n i
Q . The finite differences in the system of Eqs (12a, b, c) are replaced using the following first-order approximations, in which the unknown is the water depth Y i+1 n+1 :

n j n j n j n j n j A Y Y A A ' ) ( 1 1      , j = i, i+1, (13a) k j k i k i k i k i Q Y Y Q Q ' ) ( 1 1      , k = n, n+1, (13b)
(Delete the dots as multiplication sign in both equations) in which A' and Q' denote the derivatives of A and Q with respect to the water depth Y. The calculation of Q' uses the Manning roughness formula (10b). The details of the calculation are given in appendix. Finally, replacing the partial derivatives in Eq. ( 10a) with the approximated finite differences leads to the following equation:

  dt r r Q Y Q Y Y dx dt A Y Y A Y Q dx dt A Y n n n i n i n i n i n i n i n i n i n i n i n i n i n i ) ( ' ' ) ( ' ) ( ' ' ' 1 1 1 1 1 1 1 1 1 1 1                           (14) which gives 1 1   n i Y explicitly. To avoid numerical instabilities, the Courant number, x A t Q   , of
this explicit scheme should be close to 1 [START_REF] Cunge | Practical aspects of computational river hydraulics[END_REF].

Field Experiments

The studied plot belongs to an experimental station in the Gharb valley managed by ORMVAG (Kenitra, Morocco). The soil is mainly clayey (65-70%) with sand (20-30%).

During summer it develops large cracks under high water deficit conditions [START_REF] Mailhol | A furrow irrigation model to improve irrigation practices in the Gharb valley of Moroco[END_REF]. The plot was 230 m long and 50 m wide and was laser levelled with a slope S 0 = 0.2%. Sugar beet was cultivated for two years: 2004-2005 and 2005-2006. For both, it was planted at mid November after a furrow ridging, then furrow irrigated. The furrow spacing fs was 1.5 m. Average furrow shape was trapezoidal with a side slope m ≈1/2, and a bottom width and a maximal depth of 0.15 m. At the extremity of the plot, runoff from furrows was collected by a trapezoidal ditch with a bottom width of 0.5 m and a slope of 0.25%.

A complete meteorological station located on the experimental site recorded rainfall events thanks to an automatic rain-gauge system with a recording time step of 15 min. To measure the outflow, a triangular flume equipped with an ultrasonic sensor was installed in November 2004 at the downstream end of the ditch. Instantaneous outlet discharges were automatically recorded at a time step of 15 min by a SAB600 (LMU) system supplied by an electrical solar panel.

Five significant runoff events (peak discharge Q p ≥ 1 L/s) occurred and were monitored during the two seasons. Two events had a peak discharge close to 8 L/s (March 2005, January 2006), a third one had a peak close to 4 L/s (March 2006) while the 4 th (November 2004) and the 5 th (February 2006) had a peak of 2 L/s and 1 L/s respectively. For soil water content monitoring, two sites (upstream and downstream part of the plot) were equipped with two 30-cm CS6115 Cambell TDR probes inserted into the soil in the vicinity of the crop lines.

During the irrigation campaign of 2005, two irrigation events were monitored, on April 2 nd (3 rd irrigation) and June 6 th (6 th irrigation). Just before these two irrigation events, the soil water content profile was established by the gravimetric method on 1 m depth. Inlet discharge supplied by the gated pipe system was measured by means of trapezoidal flumes (RBC type) well adapted for low discharge (Walh et al., 2005). Total input discharge was controlled by a volumetric valve too. The advance trajectories of 18 furrows were monitored

Author-produced version of the article published in Agricultural Water Management, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014 by determining the advance times T a at nine abscissas. This experimental setup allowed the calibration of the hydraulic characteristics of the soil.

Model Calibration and Validation

The model consists of a production function, which yields the volume of water drained out of the plot, and a transfer function that represents the propagation of the flow from the plot to the outlet. For the production function, five parameters need to be calibrated: θ s , θ i , ,  c , and K s .

These parameters are first calibrated, so that the model correctly represents the runoff volume. Some of these parameters are hydraulic characteristics of the soil and may be adjusted independently of the rainfall events. To this end, these parameters were calibrated on the irrigation events of April and June 2005. In a second step, the transfer function was calibrated so that it correctly represented the delay process within the drainage system. The transfer function introduces one more parameter, the Manning roughness coefficient. The surface flow conditions are rather different during furrow irrigation and during rainfall, and this Manning coefficient should be calibrated separately in each condition.

Calibration during irrigation events

The capillarity fringe  c and the hydraulic conductivity K s are hydraulic characteristics of the soil that can be determined during irrigation events. To this end, we use the infiltration equation developed by [START_REF] Mailhol | A predictive form of Horton's equation for simulating furrow irrigation[END_REF] to simulate the furrow irrigation process:

I f (t) = 0.9  c  [ 1 -exp (K s (  c  ) -1 t )] + K s t ( 15
)
where  = θ s -θ i and  is an empirical factor set at values greater than or equal to 80 for heavy cracking soils, and close to 10 without or under weak cracking soil conditions. Unlike the Green-Ampt equation, Eq. ( 15) enables to account for the role played by the macroporosity that can be quasi instantaneously filled by an important water amount [START_REF] Mailhol | A furrow irrigation model for real time applications on cracking soils[END_REF]. The  parameter indeed allows us to reduce or to increase the role of macroporosity in the infiltration process, giving thus more flexibility to the model. Parameters  c and K s are involved in the exact analytical advance solution of the water balance equation.

They are derived from the advance monitoring of an irrigation event by an inverse method [START_REF] Mailhol | A predictive form of Horton's equation for simulating furrow irrigation[END_REF][START_REF] Mailhol | SOFIP: A field-scale modelling of water management and crop yield for furrow irrigation[END_REF]. Saturated soil water content θ s = 0.45 was obtained from laboratory experiments during which the TDR calibration was performed by the method of [START_REF] Quinones | Comparison of three calibration procedures for TDR soil moisture sensors[END_REF]. The initial soil moisture θ i is measured (average value on 1 m depth) just before irrigation. The irrigation event of April 2005 is used for parameter calibration while that of June is used for parameter validation. The average discharge at the furrow head is 1.15 L/s for both irrigations. For the downstream discharge (RBC flume measurements), the coefficient of variation is about 10%. For the advance times T a (x), the coefficient of variation is between 7 and 17%, which is consistent with the variation on the head discharges. Figure 3a shows the calibration result on the advance trajectory (Ta) of the median furrow for a value of 0.18, while Fig. 3b presents the validation step for  For both irrigation events, the roughness coefficient n is set to 0.04 [START_REF] Mailhol | A furrow irrigation model for real time applications on cracking soils[END_REF][START_REF] Mailhol | A predictive form of Horton's equation for simulating furrow irrigation[END_REF]Mailhol et al., 2005). The derived parameters from the furrow irrigation model are c  = 40 cm and s K = 0.3 cm/h.

Author-produced version of the article published in Agricultural Water Management, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014

Model calibration during rainfalls

The runoff event of March 2005 is used for model calibration, the other ones being used for model verification. The runoff of March 2005 event represents 7 mm out of a total rain of 21 mm fallen the 4 th of March in 10 hours, from 3 am to 1 pm. The events were simulated with a space step dx = 1.5 m. The time step was set to 1 s, which insured Courant number to be close to 1. Since the rain depth is measured every 15 min, the rain intensity is assumed constant during this time interval.

The calibration is composed of two steps. The first refers to the parameters of the production function while the second refers to n, the transfer model parameter. For the first step, the production function was disconnected from the transfer function. The parameters of the production function are set to the values obtained from irrigation monitoring. For the soil parameters, a value of 40 cm for  c (assimilated to -h f in the Green-Ampt model) and 0.3 cm/h for K s are adopted. Initial soil water content θ i is derived from TDR measurements. The value which is taken into consideration corresponds to the beginning of the rainfall event initiating runoff. In March 2005, the value θ i = 0.43 is obtained (Fig. 4), which is close to saturation.

Parameter governing the fast drainage function, is obtained from the soil water depletion curve (Fig. 5) from the end of the considered rainy event. The following set of parameters is obtained: = 0.00007, θ i = 0.43, c = 40 cm, K s = 0.3 cm/h. The corresponding simulations resulted in slightly under-estimated values of the runoff (6 mm instead of 7 mm). Indeed, infiltration process under furrow irrigation is different from that existing under a rainy event. Under furrow irrigation, water infiltrates the soil on the furrow surface, resulting in a 2D process, while infiltration occurs on the inter-furrow area during rainfalls. After tuning parameter K s , a slight improvement was obtained with K s = 0.2 cm/h instead of 0.3 cm/h, the value derived from the furrow model. The difference in these K s values can result from a difference of the soil temperatures between winter and summer (or end of spring) according to [START_REF] Heasom | Hyrologic modelling of a bioinfiltration best management practice[END_REF], causing differences in water viscosity. For instance, for 2005, the daily average temperature of January is 8°C while that of May is 19°C, causing the water viscosity to be 1.4 times lower in May than in June. This is consistent with variations of K s .

In a second step, the transfer function was calibrated by tuning the Manning roughness coefficient. This coefficient should be different between irrigation conditions and rainfall ones. Indeed, the water depth is much lower during rainfalls than under furrow irrigation, the lumps in the furrows causing the relative roughness to be much higher in the case of rainfall. Also, the transfer model integrates the runoff on the elementary sub-plot partially covered by the crop as well as the runoff in the furrows. The calibration phase yields a value of n = 0.2 during rainfall for all events. This value of n is higher than the one usually observed in furrows, which is consistent with the remarks above. Comparable values but in the case of runoff on a plane surface with vegetation were found by Muñoz-Carpena et al. (1998), [START_REF] Paige | A strip model approach to parameterize a coupled Green -Ampt kinematic wave model[END_REF], [START_REF] Mailhol | SPFC : a tool to improve water management and hay production in the Crau region[END_REF]. In fact, surface conditions (roughness of the soil surface, crop height) can cause some variation of the roughness coefficient with time, affecting the transfer dynamics. Such a refinement could be addressed with a large number of rainfall events, but it is not significant in the present study. The sensitivity analysis presented further shows the influence of the Manning roughness parameter on the outflow hydrograph. The quality of the simulation is attested by the high values of the coefficient of efficiency of Nash-Sutcliffe (ASCE, 1993), denoted CE, obtained after model calibration (Fig. 6). A CE value close to 1 indicates that predicted discharge is close to the measured one.

The other events are used for the validation step. For the rainy event of January 2006, the initial soil moisture is close to saturation ( i = 0.43) at the beginning of the rain, causing the runoff process. Maximal discharge is well simulated, whereas the depletion phase is partially under-predicted (Fig. 7). For November 2004, the initial moisture is lower than  s , Author-produced version of the article published in Agricultural Water Management, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014 unlike for the event of January 2006, and the rainfall intensity is too weak to initiate a significant runoff, as can be seen in Fig. 8. A comparable behaviour is observed for the event of March 2006. At the beginning of the rainfall event, no significant runoff occurs, due to a rather low value of the initial moisture ( i = 0.38), compared to that of January 2005 ( i = 0.43) and low rainfall intensity conditions. Runoff is generated when the soil moisture reaches  i ≈ 0.44 (Fig. 9). However, the peak-discharge is largely under-predicted. It seems that the furrows also act as a shallow drain system collecting flows from the superficial soil layer. The contribution of these interflows is all the more significant as the intensity of the rain generating direct runoff is low. This phenomenon is not taken into account in the model.

For the event of November 2004, the total rain height is low, but the rain intensity is rather large (15 mm/h) and, in this case, the response of the furrow system is correctly simulated. This example also points out the importance of the rainfall sampling time, which is 15 min. The rainfall distribution is unknown within this time step, and, for instance, the maximum measured value of the rainfall height (2.4 mm) for the event of February 2006 could have fallen in 5 min, resulting in a precipitation rate close to 30 mm/h, which is much larger than K s . Despite some discrepancy between the measured and simulated outflows for the events of February and March 2006 (Figs. 9 and 10), the runoff volumes are correctly simulated as shown in Fig. 11.

Sensitivity Analysis

A sensibility analysis is conducted on the first event (March 2005) used for model calibration (Fig. 12). One can see that Fig. 12a, parameter governing soil water depletion, is much less sensitive under than above its optimal value. This highlights the role of the soil water depletion model. Peak discharge decreases by 25% when  is multiplied by 10. A comparable variation (Fig. 12b) is obtained when K s is increased by 50%. Parameter | h f | (i.e.  c ) appears to be much less sensitive in this case (Fig. 12c), since the value  i is close to saturation at the beginning of the simulation. One can note a strong dissymmetry regarding n, the roughness coefficient. Increasing n by 25% results in a peak discharge decreased by 18%, whereas decreasing n by 25% results in a discharge peak increased by 40% (Fig. 12d).

Application to drainage issues

Water excess during winter is a constraint for the production system in the Gharb Valley irrigation system. The furrow irrigation system may facilitate the evacuation of the excessive water volumes during the winter rainfalls. The efficiency of such a system can be analyzed under typical rainfalls, such as events of one-year return period. These events are characterized by significant rainfall intensities and should be correctly simulated with the present model. As an illustration, a typical rainfall is generated thanks to the method proposed in the French national flood synthesis in 1981 for an event of a one-year return period. The method uses the coefficients of the Montana formula. The latter gives the rainfall intensity as a function of time (q r = a t -b ) for a given return period characterized by specific values of a and b. The values a = 35 mm/h and b = 0.6 for a one-year return period were obtained by [START_REF] Hammani | Modélisation couplée du drainage souterrain et du drainage de surface: Application à la conception du drainage agricole dans le périmètre irrigué du Gharb (Maroc)[END_REF] for the Gharb region. A characteristic time D is defined (the lag-time in the case of a watershed), generally set to the time between the maximum rainfall intensity and the peak discharge of a representative event. For a time τ between -D and D, the rain intensity is defined as follows:

Author-produced version of the article published in Agricultural Water Management, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014

7 3 10 2 2 2 2 14 20 14 ) 1 ( ) ( b b r D D D D b a q                                  . ( 16
)
The cumulative rainfall between - and + is given by:

b D D a P                  1 7 / 3 2 1 ] 2 2 ) [( ) , (      . ( 17
)
The value of D corresponding to the highest observed discharges (March 2005, January 2006) is close to 0.75 h. In order to generate a complex rainfall, the synthetic hyetograph is assumed to be preceded by moderate rainfalls (Fig. 13) filling progressively the soil water reserve as for the event of March 2005. Two different initial moisture values were used ( i = 0.35 and  i = 0.40). This example shows that, even when the simulation starts with  i far from saturation, an intense rainy event generate a substantial runoff. Depending on initial soil conditions, a peak discharge in the range of 22 to 35 L s -1 ha -1 can be reached. In the experimental plot, the evacuation system consists of trapezoidal ditches with a bottom width of 50 cm, a maximal depth of 25 cm, a side slope of 0.5 and a bottom slope of 0.2%. A Manning roughness coefficient of 0.05 leads to a conveyance capacity of 39 L/s. This conveyance capacity allows the one year return simulated event to be evacuated whatever the value of the initial moisture be. However, in the Gharb Valley, the previous ditch characteristics are designed to collect the runoff of plot groups of 6 hectares, and then most ditches have a largely insufficient conveyance capacity. Even if the furrow system is efficient to evacuate water from the plots, this excessive water should be conveyed out of the cropped lands; the present modelling approach can be used to design such an evacuation system.

Conclusion

This paper analyzed the behaviour of a furrow system under complex rainfall events. A numerical model was developed, combining the kinematic wave solution for the water transfer through the furrow system and a production function based on the Green-Ampt equation. The model was used to evaluate peak discharge and water amounts that can be evacuated from plots of approximately 1 ha. A calculation procedure was developed to account for complex rainfall events. The model was calibrated using flow measurements obtained during irrigation events and five rainfall events in 2005 and 2006. The soil parameters obtained during the irrigation events give satisfactory results when applied to the simulation of rainfall events. The parameter used in the flow propagation (the Manning roughness coefficient) within the furrows needs to be calibrated separately under irrigation and rainfall, since the flow characteristics are largely different. In the validation process, the model showed some limitations in the case of small rainfall intensities. These limitations may be addressed by adopting shorter sampling times for the rainfall monitoring, and by representing the interflows which may become significant when rain intensities are low.

As an application, the model was used to simulate the behaviour of the furrow system under a one-year return period rainfall. The results provide the peak discharge that should be evacuated out of a parcel or a group of parcels, and therefore can be used to design ditches able to drain the excessive volumes of water out of the cropped area during a rainy period.

Appendix: calculation of the derivatives dA/dY and dQ/dY

In a cross section, assumed trapezoidal, water depth is denoted Y, bottom width is denoted B and side slope is denoted m. The wetted area is given by: A = (B + mY)Y.

(A1) The derivation A'=dA/dY is straightforward:

A' = B + 2mY . (A2) The wetted perimeter P is given by: P = B+2σY (A3) where σ = (1+m²) 1/2 and the hydraulic radius R h =A/P. From Eq. ( 10b), one gets Q as a function of Y:

0 3 / 2 3 / 5 0 3 / 2 ) 2 ( ²) ( 1 S Y B mY BY n S R n A Q h       (A4)
and then its derivative Q'=dQ/dY:

                           3 / 5 3 / 2 0 2 ² 4 2 ² ) 2 ( 5 3 ' Y B mY BY Y B mY BY mY B n S Q    . ( A5
)
After a few algebraic manipulations, one gets: Author-produced version of the article published in Agricultural Water Management, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10. 1016/j.agwat.2009.02.014

                   Y B Y m BmY Y B B Y B mY BY n S Q     2 ² 16 10 6 ² 5 2 ² 3 ' 3 / 2 0 . ( A6

Fig. 3 .

 3 Fig. 3. Calibration of  c and K s on the median furrow for the 3 rd irrigation event (a) and validation on the 6 th irrigation event (b) of 2005.

Fig. 6 .

 6 Fig. 6. Model calibration on the event of March 2005 (Ks = 0.2 cm/h, n = 0.2,   = 0.43 ; CE = 0.965). Measured rainfall on a 15 min time step is converted here in mm/h.

Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .Fig. 11 .Fig. 12 .Fig. 13 .

 78910111213 Fig. 7. Model verification on the event of January 2006: (Ks = 0.2 cm/h ; n = 0.2,   ; CE = 0.889)

  Author-produced version of the article published in Agricultural WaterManagement, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014 Graphical calibration of  involved in Eq. (9) from the runoff event of March the 4 th in 2005. -produced version of the article published in Agricultural Water Management, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014

	0,3 0,35 0,4 0,45 Fig. 4. Evolution of volumetric water content measured by TDR (averaged values of the two 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Time (min)  (cm 3 /cm 3 ) Runoff event sites) from 03/03 at 0 h AM to 03/06 at 3.00 AM in 2005. Author-produced version of the article published in Agricultural Water Management, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014 y = -0,0001x -0,8409 -0,9 -0,89 -0,88 -0,87 -0,86 -0,85 -0,84 -0,83 -0,82 0 50 100 150 200 250 300 350 Time (min) Ln() Fig. 5. Author

  Author-produced version of the article published in Agricultural WaterManagement, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014

		10																0
	Q (L/s)	2 4 6 8		Simulated Measured											5 10 15 20	Rainfall (mm/h)
		0																25
		0	85	170	255	340	425	510	595	680	765	850	935	1020	1105	1190	1276	1360
									Time (min)

Author-produced version of the article published in Agricultural WaterManagement, 2009, 96 (7),1128-1136. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0378377409000596 doi:10.1016/j.agwat.2009.02.014


Lorenzo Brasco
email: brasco@cmi.univ-mrs.fr

Giovanni Franzina
email: franzina@science.unitn.it

ON THE HONG-KRAHN-SZEGO INEQUALITY FOR THE p-LAPLACE OPERATOR

Keywords: 2010 Mathematics Subject Classification. 35P15, 35P30, 49R05 Nonlinear eigenvalue problems, Hong-Krahn-Szego inequality, Stability for eigenvalues

Given an open set Ω, we consider the problem of providing sharp lower bounds for λ2(Ω), i.e. its second Dirichlet eigenvalue of the p-Laplace operator. After presenting the nonlinear analogue of the Hong-Krahn-Szego inequality, asserting that the disjoint unions of two equal balls minimize λ2 among open sets of given measure, we improve this spectral inequality by means of a quantitative stability estimate. The extremal cases p = 1 and p = ∞ are considered as well.

Introduction

In this paper, we are concerned with Dirichlet eigenvalues of the p-Laplace operator -∆ p u := -div(|∇u| p-2 ∇u), where 1 < p < ∞. For every open set Ω ⊂ R N having finite measure, these are defined as the real numbers λ such that the boundary value problem -∆ p u = λ |u| p-2 u, in Ω, u = 0, on ∂Ω has non trivial (weak) solutions. In particular, we are mainly focused on the following spectral optimization problem

(1.1) min{λ 2 (Ω) : |Ω| = c},
where c > 0 is a given number, λ 2 (•) is the second Dirichlet eigenvalue of the p-Laplacian and | • | stands for the N -dimensional Lebesgue measure. We will go back on the question of the well-posedness of this problem in a while, for the moment let us focus on the particular case p = 2: in this case we are facing the eigenvalue problem for the usual Laplace operator and, as it is well known, Dirichlet eigenvalues form a discrete nondecreasing sequence of positive real numbers 0 < λ 1 (Ω) ≤ λ 2 (Ω) ≤ λ 3 (Ω) ≤ . . . , going to ∞. In particular, it is meaningful to speak of a second eigenvalue so that problem (1.1) is well-posed and we know that its solution is given by any disjoint union of two balls having measure c/2. Moreover, these are the only sets which minimize λ 2 under volume constraint.

Using the scaling properties both of the eigenvalues of -∆ and of the Lebesgue measure, we can reformulate the previous result in scaling invariant form as follows This "isoperimetric" property of balls has been discovered (at least) three times: first by Edgar Krahn ([25]) in the '20s, but then the result has been probably neglected, since in 1955 George Pólya attributes this observation to Peter Szego (see the final remark of [START_REF] Pólya | On the characteristic frequencies of a symmetric membrane[END_REF]). However, almost in the same years as Pólya's paper, there appeared the paper [START_REF] Hong | On an inequality concerning the eigenvalue problem of membrane[END_REF] by Imsik Hong, giving once again a proof of this result: it has to be noticed that Hong's paper appeared in 1954, just one year before Pólya's one. For this reason, in what follows we will refer to (1.2) as the Hong-Krahn-Szego inequality (HKS inequality for short).

We briefly recall that for successive Dirichlet eigenvalues of the Laplacian, much less is known: apart for a result by Bucur and Henrot giving existence for k = 3 in the class of quasi-open sets (see [START_REF] Bucur | Minimization of the third eigenvalue of the Dirichlet Laplacian[END_REF]), in general existence, regularity and characterization of optimal shapes for a problem like (1.1) are still open issues. As for existence, a general (positive) answer has been given by Buttazzo and Dal Maso (see [START_REF] Buttazzo | An existence result for a class of shape optimization problems[END_REF]), for "box-constrained" spectral optimization problems, i.e. problems of the type min{λ k (Ω) :

|Ω| ≤ c, Ω ⊂ D},
where D is a given open bounded subset of R N .

For the case of the p-Laplace operator, this is clearly a completely different story: the nonlinearity of the operator and the lack of an underlying Hilbertian structure complicate things a lot. For example, though there exists a variational procedure to produce an infinite sequence of eigenvalues of -∆ p (the so called eigenvalues of Ljusternik-Schnirelmann type, see [START_REF] Fučik | Spectral Analysis of Nonlinear Operators[END_REF][START_REF] García Azorero | Existence and nonuniqueness for the p-Laplacian: nonlinear eigenvalues[END_REF] for example), up to now it is not clear whether the resulting variational spectrum coincides with the whole spectrum of -∆ p or there exist some other eigenvalues; negative answers were given in [START_REF] Binding | Variational and non-variational eigenvalues of the p-Laplacian[END_REF][START_REF] Drábek | On variational eigenvalues of the p-Laplacian which are not of Ljusternik-Schnirelmann type[END_REF] for slightly different nonlinear eigenvalue problems. Moreover, it is not even known whether or not the collection of the eigenvalues of -∆ p forms a discrete set.

This said, while it is easy to define the first eigenvalue λ 1 , in principle it becomes quite difficult even to start speaking of the second eigenvalue, the third one and so on, since discreteness of the spectrum is not guaranteed. However, as it is well known, it turns out that also in the case of -∆ p one can speak of a second eingevalue λ 2 : this means that there is a gap between λ 1 and λ 2 , as for p = 2. Moreover, this second eigenvalue is a variational one (see Section 2 for more details).

The main aim of the present paper is the study of the spectral optimization problem (1.1) for a general 1 < p < ∞: as we will see, the Hong-Krahn-Szego inequality still holds in the case of the p-Laplace operator (Theorem 3.2). Namely, any disjoint union of two equal balls minimizes the second eigenvalue of -∆ p among sets of given measure, that is

(1.3) |Ω| p/N λ 2 (Ω) ≥ 2 p/N ω p/N N λ 1 (B)
. The proof runs very similarly to the case p = 2 and it is based exactly on the same two ingredients, which still hold in the nonlinear setting:

• the Faber-Krahn inequality (see next section) for the first eigenvalue of -∆ p ;

• the fact that for a connected open set the only eigenfunction of constant sign is the first one.

We will then turn our attention to the stability issue: indeed, when dealing with shape optimization problems having unique solution (possibly up to some suitable group of rigid transformations, like rotations or translations, for example), a very interesting and natural question is to know whether this optimal shape is stable or not. For example, specializing this question to our problem (1.1), we are interested in addressing the following issue:

|Ω 0 | = c and λ 2 (Ω 0 ) ≃ min{λ 2 (Ω) : |Ω| = c} ?
=⇒ Ω 0 "near" to the optimal shape

In this paper, we give a positive answer to this question, by proving a quantitative version of (1.3). By "quantitative" we mean the following: actually, inequality (1.3) can be improved by adding a reminder term, which measures (in a suitable sense) the distance of the generic set Ω from the "manifold" of optimizers O, i.e. the collection of all disjoint unions of two equal balls. Then the result we provide (Theorem 4.2) is an improvement of (1.3) of the type

|Ω| p/N λ 2 (Ω) -2 p/N ω p/N N λ 1 (B) ≥ Φ(d(Ω, O)), where d(•, O) is a suitable "distance" from O and Φ : [0, ∞) → [0, ∞)
is a continuous strictly increasing function, with Φ(0) = 0: more precisely, in Theorem 4.2 we prove a result like this with d equal to the L 1 distance of characteristic functions (a variant of the so-called Fraenkel asymmetry, see Section 4) and Φ given by a power function. This quantitative estimate guarantees that if the second eigenvalue λ 2 of a set Ω is almost equal to the λ 1 of a ball having measure |Ω|/2, then Ω is almost the disjoint union of two equal balls, i.e. we have stability of optimal shapes for our spectral optimization problem. Our analysis will cover the whole range of p: indeed, we will show that the same proof can be adapted to cover the cases p = 1 and p = ∞ as well, where λ 2 becomes the second Cheeger constant and the second eigenvalue of the ∞-Laplacian, respectively (see Section 5 for the precise definitions). In the case of the first eigenvalue λ 1 , we recall that quantitative results of this type have been derived in [START_REF] Hansen | Isoperimetric inequalities in potential theory[END_REF][START_REF] Melas | The stability of some eigenvalue estimates[END_REF][START_REF] Povel | Confinement of Brownian motion among Poissonian obstacles in R d , d ≥ 3[END_REF] (linear case, p = 2) and [START_REF] Bhattacharya | Some observations on the first eigenvalue of the p-Laplacian and its connections with asymmetry[END_REF][START_REF] Fusco | Stability estimates for certain Faber-Krahn, Isocapacitary and Cheeger inequalities[END_REF] (general case, 1 < p < ∞).

We point out that though problem (1.1) is a very natural one also for -∆ p , we have not been able to find in literature any paper recording a proof of (1.3): only after the completion of this work, we found a related recent paper by Kennedy ([24]), dealing with problem (1.1), but for the second eigenvalue of -∆ p with Robin boundary conditions. For this reason, we decided to write properly the complete proof of (1.3); on the contrary, the quantitative stability results of Theorem 4.2 and Theorem 5.2 in this paper are certainly new, though probably not sharp, except for the case p = ∞ (see the discussion in Section 6).

We conclude this introduction with the plan of the paper: in order to make the work as selfcontained as possible, Section 2 recalls the basic facts about the first two eigenvalues of -∆ p that we will need in the following; in Section 3 we prove the Hong-Krahn-Szego inequality for λ 2 , while Section 4 provides a quantitative version of the latter, thus extending to the nonlinear case a result recently proven in [START_REF] Brasco | Sharp stability of some spectral inequalities[END_REF]. In Section 5, for the sake of completeness, we consider the shape optimization problem (1.1) for the "extremal" cases, i.e. for p = 1 and p = ∞: in this case, the first two eigenvalues λ 1 and λ 2 become two purely geometrical objects and we study stability of optimal shapes for them. Finally, Section 6 concludes the paper with some examples, remarks and conjectures concerning the sharpness of the quantitative estimates derived in this work.

Tools: the first two eigenvalues of -∆ p

Given an open set Ω ⊂ R N having finite measure and p ∈ (1, ∞), we define the L p unitary sphere

B p (Ω) = {u ∈ L p (Ω) : u L p (Ω) = 1},
and we indicate with W 1,p 0 (Ω) the usual Sobolev space, given by the closure of C ∞ 0 (Ω) with respect to the norm

u = Ω |∇u(x)| p dx 1/p .
If for a certain λ we have that there exists a non trivial u ∈ W 1,p 0 (Ω) satisfying (2.1) -∆ p u = λ|u| p-2 u, in Ω, in a weak sense, i.e.

Ω |∇u(x)| p-2 ∇u(x), ∇ϕ(x) dx = λ Ω |u(x)| p-2 u(x) ϕ(x) dx, for every ϕ ∈ W 1,p 0 (Ω),
then we call λ a Dirichlet eigenvalue of -∆ p in Ω: correspondingly u will be a Dirichlet eigenfunction of -∆ p . In particular, observe that for every such a pair (λ, u) there results

Ω |∇u(x)| p dx = λ Ω |u(x)| p dx.
Though we will not need it in the sequel, we recall that it is possible to show the existence of a diverging sequence of eigenvalues of -∆ p , see [START_REF] Fučik | Spectral Analysis of Nonlinear Operators[END_REF][START_REF] García Azorero | Existence and nonuniqueness for the p-Laplacian: nonlinear eigenvalues[END_REF].

Remark 2.1. Observe that in general solutions of (2.1) are just in C 1,α (see [START_REF] Dibenedetto | C 1+α local regularity of weak solutions of degenerate elliptic equations[END_REF]): in particular they may not posses second derivatives even in a weak sense (this is the case only when p ≤ 2, see [START_REF] Acerbi | Regularity of minimizers of non-quadratic functionals: the case 1 < p < 2[END_REF]). Then eigenfunctions in general are not classical solutions of the equation (2.1).

The first Dirichlet eigenvalue of the p-Laplacian of a set has the following variational definition

(2.2) λ 1 (Ω) = min u∈Bp(Ω)∩W
λ 1 (Ω) = min u∈Bp(Ω)∩W 1,p 0 (Ω) Ω |∇u(x)| p dx ≤ Ω |∇v(x)| p dx = λ.
The first important result that we need concerns the simplicity of λ 1 on a connected open set.

Theorem 2.3. Let Ω ⊂ R N be an open connected set, having |Ω| < +∞. Then λ 1 (Ω) is simple, i.
e. the corresponding eigenfuctions form a 1-dimensional linear space.

Proof. A very short and elegant proof of this fact can be found in the paper [START_REF] Belloni | A direct uniqueness proof for equations involving the p-Laplace operator[END_REF] by Belloni and Kawohl: actually, the result in [START_REF] Belloni | A direct uniqueness proof for equations involving the p-Laplace operator[END_REF] is stated for the case of Ω being a bounded set, but it can be easily seen that this hypothesis plays no role and the same proof still works for Ω having finite measure.

Throughout the paper, we will use the following convention: when Ω is a disconnected open set, the set of its Dirichlet eigenvalues is made of the collection of the eigenvalues of its connected components. The eigenvalues are obtained by gathering and ordering increasingly the eigenvalues on the single pieces; correspondingly, each eigenfunction is solution of (2.1) on a certain connected component and vanishes on the others.

The following result plays a crucial role: it asserts that any eigenfunction having constant sign is the first one on some connected component of the open set Ω.

Theorem 2.4. Let Ω ⊂ R N be an open set, having finite measure. Let u ∈ W 1,p 0 (Ω) be a Dirichlet eigenfunction relative to some eigenvalue λ. If u has constant sign in Ω, then λ = λ 1 (Ω 0 ) for some connected component Ω 0 of Ω, i.e. u is a first eigenfunction of Ω 0 . In particular λ = λ 1 (Ω) if Ω is connected.
Proof. For connected sets, this has been first proven in [START_REF] Ôtani | On the first eigenvalue of some quasilinear elliptic equations[END_REF] for the case of smooth domains, while for a simple proof in the general case we refer to [START_REF] Kawohl | Positive eigenfunctions for the p-Laplace operator revisited[END_REF]Lemma 3.1]. Again, it has to be noticed that in [START_REF] Kawohl | Positive eigenfunctions for the p-Laplace operator revisited[END_REF], the set Ω is required to be bounded for this result to hold: this guarantees that for such a set Ω there exists a sequence of smooth sets

{Ω k } k≥1 such that Ω = k≥1 Ω k and lim k→∞ λ 1 (Ω k ) = λ 1 (Ω).
For the sake of completeness, we briefly sketch how this property remains true also under the assumption |Ω| < ∞: then the rest of the proof can be carried over exactly as in Lemma 3.1 of [START_REF] Kawohl | Positive eigenfunctions for the p-Laplace operator revisited[END_REF]. Observe that setting Ω k := Ω ∩ {x : |x| < k}, we have

Ω = ∞ k=1 Ω k = ∞ n=1 Ω k,n , since each Ω k is a bounded open set, with Ω k,1 ⊂ Ω k,2 ⊂ • • • ⊂ Ω k and Ω k,n smooth. Collecting all of them {Ω n,k } n,k≥1
, we obtain an exhaustion of Ω by means of smooth sets: moreover relabeling these sets by {Ω j } j≥1 and extracting (if necessary) a suitable subsequence, it is not difficult to see that we also have lim

j→∞ λ 1 (Ω j ) = λ 1 (Ω).
Then the proof can be carried on exactly as in Lemma 3.1 of [START_REF] Kawohl | Positive eigenfunctions for the p-Laplace operator revisited[END_REF].

Finally, if Ω is disconnected, then λ has to be a Dirichlet eigenvalue of a certain connected component Ω 0 ; correspondigly u is an eigenfunction of Ω 0 , having constant sign. Then it sufficies to apply the first part of the proof to conclude.

We give now a precise definition of what we mean by the second eigenvalue of -∆ p : this definition keeps into account the multiplicity of the first eigenvalue. As we will see, this is necessary in order to properly deal with our spectral optimization problem (1.1).

Definition.

Let Ω be an open set having finite measure. Then its second eigenvalue is given by

(2.3) λ 2 (Ω) =    min{λ > λ 1 (Ω) : λ is an eigenvalue } if λ 1 (Ω) is simple λ 1 (Ω) otherwise.
Some words about the consistency of this definition are in order, when λ 1 is simple. Indeed, using Theorem 2.4 it can be proven that if Ω is connected, then λ 1 is isolated in the spectrum, the latter being a closed set ([26, Theorem 3]): this shows that the minimum in (2.3) is well-defined. On the other hand, if Ω consists of infinitely many connected components, one only has to check that the collection of the first eigenvalues on the single components cannot accumulate at any value λ ≥ λ 1 (Ω), but this follows by combining the assumption |Ω| < ∞ and the Faber-Krahn inequality (see below). This shows again that the minimum in (2.3) is meaningful.

Remark 2.5. By definition, the nodal domains of an eigenfunction u are the connected components of the sets {x : u(x) > 0} and {x : u(x) < 0}. If Ω is connected, we recall that every eigenfunction corresponding to λ 2 has exactly two nodal domains (see [START_REF] Cuesta | A nodal domain property for the p-Laplacian[END_REF]), in which case by Theorem 2.4 we can infer λ 2 (Ω) = min{λ > λ 1 (Ω) : λ admits a sign-changing eigenfunction}.

For the sake of completeness, we recall that one can give a variational characterization also for λ 2 : in order to introduce it, we need some further notations. Given a pair of functions u, v ∈ B p (Ω) ∩ W 1,p 0 (Ω), let us denote by Γ Ω (u, v) the set of continuous (in the W 1,p topology) paths in

B p (Ω) ∩ W 1,p 0 (Ω) connecting u to v, i.e. Γ Ω (u, v) = γ : [0, 1] → B p (Ω) ∩ W 1,p 0 (Ω) : γ is continuous and γ(0) = u, γ(1) = v .
Theorem 2.6. Let Ω ⊂ R N be an open set having finite measure, not necessarily connected. Let u 1 ∈ B p (Ω) ∩ W 1,p 0 (Ω) be one of its first eigenfunctions. Then λ 2 (Ω) has the following mountain pass characterization

(2.4) λ 2 (Ω) = inf γ∈Γ Ω (u 1 ,-u 1 ) max u∈γ([0,1]) Ω |∇u(x)| p dx.
If λ 1 (Ω) is not simple, this characterization is independent of the particular u 1 we choose.

Proof. If Ω is connected, this has been proven in the paper [START_REF] Cuesta | The beginning of the Fuĉik spectrum for the p-Laplacian[END_REF], to which we refer for the proof.

Here, we just show how (2.4) can be extended to the case of general open sets: anyway, since we will not need this result in the sequel, the uninterested reader may skip the proof at a first reading. Let us take Ω not connected, then the following alternative holds: either λ 1 (Ω) is simple or not.

Case λ 1 simple: in this case, a first eigenfunction

u 1 ∈ B p (Ω) ∩ W 1,p 0 (Ω) is unique and let us consider a second eigenfunction u 2 , still normalized by u 2 L p (Ω) = 1.
If u 2 is sign-changing, then it is supported on some connected component Ω 0 of Ω such that λ i (Ω) = λ i (Ω 0 ), i = 1, 2: in particular the mountain pass characterization of λ 2 (Ω 0 ) holds, with the maximum performed on the restricted class of curves Γ Ω 0 (u 1 , -u 1 ) ⊂ Γ Ω (u 1 , -u 1 ). Thus setting

λ := inf γ∈Γ Ω (u 1 ,-u 1 ) max u∈γ([0,1]) Ω |∇u(x)| p dx
on the one hand we have λ ≤ λ 2 (Ω 0 ) = λ 2 (Ω), while on the other hand we get λ 1 (Ω) < λ, since λ gives a Dirichlet eigenvalue of -∆ p in any case (see [START_REF] Cuesta | The beginning of the Fuĉik spectrum for the p-Laplacian[END_REF]Section 2]). Summarizing, we obtain λ 1 (Ω) < λ ≤ λ 2 (Ω) which gives the thesis in this case, thanks to (2.3). On the contrary, if λ 1 (Ω) is simple but u 2 has constant sign, then we have λ 1 (Ω) = λ 1 (Ω 0 ) and λ 2 (Ω) = λ 1 (Ω 1 ), with Ω 0 and Ω 1 distinct connected components. We construct a special element of Γ Ω (u 1 , -u 1 ), a continuous path γ defined as follows

γ(t) = ϕ t , with ϕ t (x) = cos(π t)u 1 (x) + t(1 -t)u 2 (x) (| cos(π t)| p + t p (1 -t) p ) 1/p , x ∈ Ω, for all t ∈ [0, 1].
It is easy to see that γ has the following properties

γ(t) ∈ B p (Ω) ∩ W 1,p 0 (Ω), for every t ∈ [0, 1], and
γ(0) = u 1 , γ(1) = -u 1 ,
i.e. the curve γ is admissible for the variational problem (2.4). Hence we get

λ ≤ max t∈[0,1] Ω |∇ϕ t (x)| p dx = max t∈[0,1] | cos(π t)| p λ 1 (Ω 0 ) + t p (1 -t) p λ 1 (Ω 1 ) | cos(π t)| p + t p (1 -t) p ≤ λ 1 (Ω 1 ) = λ 2 (Ω),
where we used λ 1 (Ω 0 ) < λ 1 (Ω 1 ). Thus we get λ 1 (Ω) < λ ≤ λ 2 (Ω) and we can conclude as before.

Case λ 1 multiple: if Ω is not connected and its corresponding first eigenvalue is not simple, we just take two linearly independent first eigenfunctions u 1 , u 2 ∈ B p (Ω) ∩ W 1,p 0 (Ω), which are thus supported on different connected components of Ω. Repeating the construction of the curve γ above, we obtain

λ 1 (Ω) ≤ λ = inf γ∈Γ Ω (u,-u) max u∈γ([0,1]) Ω |∇u(x)| p dx ≤ λ 1 (Ω), which shows that λ = λ 1 (Ω) = λ 2 (Ω).
Observe that if we exchange the role of u 1 and u 2 , we still arrive at the same conclusion, thus proving that in this case formula (2.4) is independent of the choice of the particular first eigenfunction.

Remark 2.7. It is useful to recall at this point that usually the variational eigenvalues {λ k } k≥1 of -∆ p are defined through a minmax problem on B p ∩ W 1,p 0 for the integral Ω |∇u| p , involving the concept of Krasnosel'skii genus: the previous result gives in particular that for k = 2 this characterization coincides with the one given by (2.4).

Finally, since our aim is that of considering a particular class of shape optimization problems involving the spectrum of -∆ p , we conclude this introduction by recalling some further properties of λ 1 and λ 2 that we will need in the sequel. In particular, they are monotone decreasing with respect to set inclusion, while as for their scaling properties we have

λ i (t Ω) = t -p λ i (Ω), t > 0, i = 1, 2,
which in particular implies that the shape functional Ω → |Ω| p/N λ i (Ω) is scaling invariant. Thus the two problems

min{λ i (Ω) : |Ω| = c} and min |Ω| p/N λ i (Ω), i = 1, 2,
are equivalent, in the sense that they both provide the same optimal shapes, up to a scaling. For i = 1, the solution to the previous problem is given by any ball: this is the celebrated Faber-Krahn inequality. The classical proof combines the Schwarz symmetrization with the so called Pólya-Szegő principle (see [START_REF] Henrot | Extremum problems for eigenvalues of elliptic operators[END_REF]Chapter 3], for example).

Faber-Krahn Inequality. Let 1 < p < ∞. For every open set Ω ⊂ R N having finite measure, we have

(2.5) |Ω| p/N λ 1 (Ω) ≥ ω p/N N λ 1 (B),
where B is the N -dimensional ball of radius 1 and ω N := |B|. Moreover, equality sign in (2.5) holds if and only if Ω is a ball. In other words, for every c > 0 the unique solutions of the following spectral optimization problem

min{λ 1 (Ω) : |Ω| = c},
are given by balls having measure c.

The Hong-Krahn-Szego inequality

In this section, we are going to prove that the disjoint unions of equal balls are the only sets minimizing λ 2 under volume constraint, i.e. we will prove the Hong-Krahn-Szego inequality for the p-Laplacian. The key step in the proof is the following technical result: this is an adaptation of a similar result for the linear case p = 2 (see [7, Lemma 3.1], for example).

Lemma 3.1. Let Ω ⊂ R N be an open set with |Ω| < ∞. Then there exists Ω + , Ω -disjoint subsets of Ω such that (3.1) λ 2 (Ω) = max{λ 1 (Ω + ), λ 1 (Ω -)}.
Proof. Let us take u 1 , u 2 ∈ B p (Ω) ∩ W 1,p 0 (Ω) a first and second eigenfunction, respectively: notice that if λ 1 (Ω) is not simple, we mean that u 1 and u 2 are two linearly independents eigenfunctions corresponding to λ 1 (Ω). We can distinguish two alternatives:

(i) u 2 is sign-changing;

(ii) u 2 has constant sign in Ω. Let us start with (i): in this case, by Remark 2.5 u 2 has exactly two nodal domains

Ω + = {x ∈ Ω : u 2 (x) > 0}
and

Ω -= {x ∈ Ω : u 2 (x) < 0},
which by definition are connected sets. The restriction of u 2 to Ω + is an eigenfunction of constant sign for Ω + , then Theorem 2.4 implies that u 2 must be a first eigenfunction for it. Replacing Ω + with Ω -, the previous observation leads to

λ 2 (Ω) = λ 1 (Ω -) = λ 1 (Ω + ).
which implies in particular (3.1) in this case.

In case (ii), let us set

Ω + = {x ∈ Ω : |u 1 (x)| > 0} and Ω -= {x ∈ Ω : |u 2 (x)| > 0}.
Using Theorem 2.4, we have that Ω + and Ω -have to be two distinct connected components of Ω: in addition u 1 , u 2 are eigenfunctions (with constant sign) of Ω + and Ω -, respectively. Then

λ 1 (Ω -) = Ω - |∇u 2 (x)| p dx = Ω |∇u 2 (x)| p = λ 2 (Ω).
Clearly, we also have λ 1 (Ω + ) = λ 1 (Ω) ≤ λ 2 (Ω), which finally gives (3.1) also in this case.

We are now ready for the main result of this section.

Theorem 3.2 (HKS inequality for the p-Laplacian). For every Ω ⊂ R N open set having finite measure, we have

(3.2) |Ω| p/N λ 2 (Ω) ≥ 2 p/N ω p/N N λ 1 (B),
where B is the N -dimensional ball of radius 1 and ω N := |B|. Moreover, equality sign in (3.2) holds if and only if Ω is the disjoint union of two equal balls. In other words, for every c > 0 the unique solutions of the following spectral optimization problem

min{λ 2 (Ω) : |Ω| = c},
are given by disjoint unions of two balls, both having measure c/2.

Proof. With the notation of Lemma 3.1, an application of the Faber-Krahn inequality yields

(3.3) λ 2 (Ω) = max{λ 1 (Ω + ), λ 1 (Ω -)} ≥ max{λ 1 (B + ), λ 1 (B -)},
where B + , B -are balls such that

|B + | = |Ω + | and |B -| = |Ω -|.
Thanks to the scaling properties of λ 1 , we have

λ 1 (B + ) = ω N |Ω + | p/N λ 1 (B) and λ 1 (B -) = ω N |Ω -| p/N λ 1 (B),
so that from (3.3) we obtain

λ 2 (Ω) ≥ ω p/N N λ 1 (B) max{|Ω + | -p/N , |Ω -| -p/N }. Finally, observe that since |Ω + | + |Ω -| ≤ |Ω|, we get (3.4) max{|Ω + | -p/N , |Ω -| -p/N } ≥ |Ω| 2 -p/N
, which concludes the proof of the inequality.

As for the equality cases, we start observing that we just used two inequalities, namely (3.3) and (3.4): on the one hand, equality in (3.3) holds if and only if at least one among the two subsets is a ball, say Ω + = B + , with λ 1 (B + ) ≥ λ 1 (Ω -); on the other hand, if equality holds in (3.4) then we must have

|Ω + | = |Ω -| = |Ω|/2.
Since Ω + and Ω -have the same measure and the one with the greatest λ 1 is a ball, we can conclude that both have to be a ball, thanks to the equality cases in the Faber-Krahn inequality.

The stability issue

We now come to the question of stability for optimal shapes of λ 2 under measure constraint. In particular, we will enforce the lower bound on |Ω| 2/N λ 2 (Ω) provided by the Hong-Krahn-Szego inequality, by adding a remainder terms in the right-hand side of (3.2). At this aim, we need to introduce some further tools: given an open set Ω ⊂ R N having |Ω| < ∞, its Fraenkel asymmetry is defined by

A(Ω) = inf 1 Ω -1 B L 1 |Ω| : B is a ball such that |B| = |Ω| .
This is a scaling invariant quantity such that 0 ≤ A(Ω) < 2, with A(Ω) = 0 if and only if Ω coincides with a ball, up to a set of measure zero. Then we recall the following quantitative improvement of the Faber-Krahn inequality, proven in [START_REF] Bhattacharya | Some observations on the first eigenvalue of the p-Laplacian and its connections with asymmetry[END_REF] (case N = 2) and [START_REF] Fusco | Stability estimates for certain Faber-Krahn, Isocapacitary and Cheeger inequalities[END_REF] (general case): for every Ω ⊂ R N open set with |Ω| < ∞, we have (4.1)

|Ω| p/N λ 1 (Ω) ≥ ω p/N N λ 1 (B) [1 + γ N,p A(Ω) κ 1 ]
, where γ N,p is a constant depending only on N and p and the exponent κ 1 = κ 1 (N, p) is given by

κ 1 (N, p) = 3, if N = 2, 2 + p, if N ≥ 3.
Remark 4.1. One may ask wheter the exponent κ 1 in (4.1) is sharp or not: introducing the deficit

F K(Ω) := |Ω| p/N λ 1 (Ω) ω p/N N λ 1 (B) -1,
one would like to prove the existence of suitable deformations {Ω ε } ε>0 of a ball B, such that

lim ε→0 F K(Ω ε ) = 0 and lim ε→0 A(Ω ε ) κ 1 F K(Ω ε ) = ℓ = {0, +∞}.
i.e. the asymmetry to the power κ 1 and the deficit have the same decay rate to 0. At least in the case p = 2, the answer should be no, since the conjectured sharp exponent is 2 (see [5, pag. 56]), while κ 1 (N, 2) ≥ 3: at present, a proof of this fact still lacks.

In the case of the Hong-Krahn-Szego inequality, the relevant notion of asymmetry is the Fraenkel 2-asymmetry, introduced in [7]

A 2 (Ω) = inf 1 Ω -1 B 1 ∪B 2 L 1 |Ω| : B 1 , B 2 balls such that |B 1 ∩ B 2 | = 0, |B i | = |Ω| 2 , i = 1, 2 .
The main result of this section is the following quantitative version of Theorem 3.2.

Theorem 4.2. Let Ω ⊂ R N be an open set, with |Ω| < ∞ and p ∈ (1, ∞). Then (4.2) |Ω| p/N λ 2 (Ω) ≥ 2 p/N ω p/N N λ 1 (B) [1 + C N,p A 2 (Ω) κ 2 ]
, with C N,p constant depending on N and p only and κ 2 = κ 2 (N, p) given by

κ 2 (N, p) = κ 1 (N, p) • N + 1 2 .
Proof. Thanks to Lemma 3.1, we have existence of two disjoint sets Ω + , Ω -⊂ Ω such that (3.1) holds. We then set

δ + = |Ω + | - |Ω| 2 , δ -= |Ω -| - |Ω| 2 ,
and we observe that it must be δ + + δ -≤ 0, since

|Ω + | + |Ω -| ≤ |Ω|.
To simplify a bit formulas, let us introduce the deficit functional

HKS(Ω) := |Ω| p/N λ 2 (Ω) 2 p/N ω p/N N λ 1 (B) -1.
In order to prove (4.2), we just need to show that

(4.3) HKS(Ω) ≥ C N,p max A(Ω + ) κ 1 + δ + |Ω| , A(Ω -) κ 1 + δ - |Ω| ,
then a simple application of Lemma 4.3 below will conclude the proof. To obtain (4.3), it will be useful to distinguish between the case p ≤ N and the case p > N . For both of them, we will in turn treat separately the case where both δ + and δ -are non positive and the case where they have opposite sign. Finally, since the quantities appearing in the right-hand side of (4.3) are all universally bounded, it is not restrictive to prove (4.3) under the additional assumption (4.4)

HKS(Ω) ≤ 1 4 .
Indeed, it is straightforward to see that if HKS(Ω) > 1/4 then (4.3) trivially holds with constant

C N,p = 1 2 1 2 κ 1 +1 + 1 .
Case A: p ≤ N . In this case the proof runs very similarly to the linear case p = 2 treated in [START_REF] Brasco | Sharp stability of some spectral inequalities[END_REF].

We start applying the quantitative Faber-Krahn inequality (4.1) to Ω + : if we indicate with B the ball of unit radius, recalling (3.1) and using the definition of δ + , we find

γ N,p A(Ω + ) κ 1 ≤ |Ω + | p/N λ 1 (Ω + ) ω p/N N λ 1 (B) -1 ≤ (|Ω| + 2δ + ) p/N λ 2 (Ω) 2 p/N ω p/N N λ 1 (B) - 1
Since p ≤ N , the power function t → (|Ω| + t) p/N is concave, thus we have

(|Ω| + 2δ + ) p/N ≤ |Ω| p/N + 2p N |Ω| p/N δ + |Ω| .
Using this information in the previous inequality, we get

γ N,p A(Ω + ) κ 1 ≤ HKS(Ω) + 2p N δ + |Ω| |Ω| p/N λ 2 (Ω) 2 p/N ω p/N N λ 1 (B)
, that we can rewrite as follows

(4.5) γ N,p A(Ω + ) κ 1 ≤ HKS(Ω) + 2p N δ + |Ω| (HKS(Ω) + 1) .
Replacing Ω + with Ω -, one obtains a similar estimate for Ω -.

Case A.1: δ + and δ -are both non-positive. In this case, it is enough to observe that HKS(Ω) ≥ 0 while δ + ≤ 0, thus from the previous we get

γ N,p A(Ω + ) κ 1 + 2p N |δ + | |Ω| ≤ HKS(Ω).
The same calculations with Ω -in place of Ω + yield (4.3).

Case A.2: δ + and δ -have opposite sign. Let us assume for example that δ + ≥ 0 and δ -≤ 0: the main difference with the previous case is that now the larger piece Ω + could be so large that the information provided by (3.1) is useless. However, estimate (4.5) still holds true for both Ω + and Ω -: using this and the fact that δ + + δ -≤ 0, we can thus infer

HKS(Ω) ≥ - 2p N δ - |Ω| ≥ 2p N δ + |Ω| ,
i.e. the deficit is controlling the error term |δ + |/|Ω|. To finish, we still have to control the asymmetry of the larger piece Ω + in terms of the deficit: it is now sufficient to introduce the previous information into (4.5), thus getting

γ N,p A(Ω + ) κ 1 ≤ HKS(Ω)(2 + HKS(Ω)).
Since we are assuming HKS(Ω) ≤ 1/4, the previous implies that HKS(Ω) controls A(Ω + ) κ 1 , modulo a constant depending only on N and p. These estimates on Ω + , together with the validity of (4.5) for Ω -and with the fact that δ -≤ 0, ensure that (4.3) holds also in this case.

Case B: p > N . Let us start once again with Ω + . Using (3.1) and the quantitative Faber-Krahn (4.1) as before, we get

HKS(Ω) ≥ |Ω| p/N λ 1 (Ω + ) 2 p/N ω p/N N λ 1 (B) -1 ≥ |Ω| 2 |Ω + | p/N (1 + γ N,p A(Ω + ) κ 1 ) -1 .
Then using the definition of δ + and the convexity of the function t → (1 + t) p/N (since p > N ), we have

|Ω| 2 |Ω + | p/N = 1 - δ + |Ω + | p/N ≥ 1 - p N δ + |Ω + | .
Inserted in the previous estimate, this yields

(4.6) HKS(Ω) ≥ γ N,p 1 - p N δ + |Ω + | A(Ω + ) κ 1 - p N δ + |Ω + | .
In the same way, using Ω -in place of Ω + , we obtain a similar estimate for Ω -.

Case B.1: δ + and δ -are both non positive. In this case, in (4.6) we can drop the terms

- p N δ + |Ω + | γ N,p A(Ω + ) κ 1 and - p N δ - |Ω -| γ N,p A(Ω -) κ 1 ,
since these are positive, thus we arrive once again at (4.3), keeping into account that

- δ + |Ω + | ≥ - δ + |Ω| and - δ - |Ω -| ≥ - δ - |Ω| .
Case B.2: δ + and δ -have opposite sign. Let us suppose as before that δ + ≥ 0 and δ -≤ 0. Now the main problem is the term in front of the asymmetry A(Ω + ) in (4.6), which could be negative. Since δ + + δ -≤ 0, applying (4.6) to Ω -we obtain (4.7)

δ + |Ω| ≤ - δ - |Ω| ≤ N p HKS(Ω).
We then observe that if

(4.8) δ + ≤ N p |Ω| 4 ,
we have

1 - p N δ + |Ω + | ≥ 1 - 1 4 |Ω| |Ω + | ≥ 1 2 ,
thanks to the fact that |Ω| ≤ 2 |Ω + |, which easily follows from the assumption that δ + ≥ 0. From (4.6) we can now infer

HKS(Ω) ≥ γ N,p 2 A(Ω + ) κ 1 - p N δ + |Ω + | , then (4.
3) follows as before, using (4.7) and the fact that

- δ + |Ω + | ≥ -2 δ + |Ω| .
This would prove the thesis under the additional hypothesis (4.8): however, if this is not satisfied, then (4.7) would imply HKS(Ω) > 1/4, which is in contrast with our assumption (4.4).

The following technical Lemma of geometrical content completes the proof of Theorem 4.2: this is the same as [7, Lemma 3.3] and we omit the proof.

(4.9) A 2 (Ω) ≤ C N A(Ω + ) + 1 2 - |Ω + | |Ω| + A(Ω -) + 1 2 - |Ω -| |Ω| 2/(N +1)
, for a suitable dimensional constant C N .

Extremal cases: p = 1 and p = ∞

To complete the analysis of our spectral optimization problem in the nonlinear setting, it is natural to give a brief look at what happens for (1.1), when p tends to the extrema of its possible range, i.e. p = 1 and p = ∞: it is known that in these cases, some shape functionals of geometric flavour appear, in place of the eigenvalues of an elliptic operator.

To enter more in this question, we need some definitions: for Ω ⊂ R N open set with |Ω| < ∞, C 1 (Ω) and C 2 (Ω) stand for the first two Cheeger constants, which are defined respectively by

C 1 (Ω) = inf E⊂Ω P (E) |E| and C 2 (Ω) = inf t : there exist E 1 , E 2 ⊂ Ω such that E 1 ∩ E 2 = ∅ and max i=1,2 P (E i ) |E i | ≤ t ,
where P (E) equals the distributional perimeter of a set E if this is a finite perimeter sets and is +∞ otherwise. Also, if |E| = 0 we use the convention P (E)/|E| = +∞. We denote by Λ 1 (Ω) the inverse of the radius r 1 of the largest ball included in Ω, while Λ 2 (Ω) will denote the inverse of the largest positive number r 2 such that there exist two disjoint balls of radius r 2 contained in Ω: it is remarkable to notice that Λ 1 and Λ 2 are indeed two eigenvalues, precisely they coincide with the first two eigenvalues of the ∞-Laplacian (see [START_REF] Juutinen | On the higher eigenvalues for the ∞-eigenvalue problem[END_REF]).

Our interest in these quantities is motivated by the following Theorem, collecting various results about the asymptotic behaviour of λ 1 and λ 2 .

Limiting behaviour of eigenvalues. For every set Ω ⊂ R N , there holds (5.1) lim

p→1 + λ i (Ω) = C i (Ω), i = 1, 2 and lim p→∞ λ i (Ω) 1/p = Λ i (Ω), i = 1, 2.
Proof. The first fact is proven in [START_REF] Fridman | Isoperimetric estimates for the first eigenvalue of the p-Laplace operator and the Cheeger constant[END_REF] and [START_REF] Parini | The second eigenvalue of the p-Laplacian as p goes to 1[END_REF], respectively. For the second, one can consult [START_REF] Juutinen | On the higher eigenvalues for the ∞-eigenvalue problem[END_REF] and the references therein.

Remark 5.1. At this point, one could be tempted to use the previous results for λ 1 , in order to improve inequality (4.1): for example, using the subadditivity of the function t → (1 + t) 1/p , it is not difficult to see that

lim p→∞ F K(Ω) 1/p ≥ |Ω| 1/N Λ 1 (Ω) ω 1/N N -1 ≥ 1 2N A(Ω),
where in the last inequality we used the (sharp) quantitative stability estimate 1 for Λ 1 (see [START_REF] Juutinen | On the higher eigenvalues for the ∞-eigenvalue problem[END_REF], equation (2.6)). Then one could bravely guess that for p "very large", inequality (4.1) has to hold with the exponent κ 1 (N, p) replaced by p, which is strictly small if N ≥ 3: this would prove that 1 The relation between the Fraenkel asymmetry α(Ω) as defined in [START_REF] Juutinen | On the higher eigenvalues for the ∞-eigenvalue problem[END_REF] and our definition is given by A(Ω) = 2 α(Ω).

This explains the discrepancy between our constant 1/(2N ) and the constant 1/N that can be found in [START_REF] Juutinen | On the higher eigenvalues for the ∞-eigenvalue problem[END_REF], equation (2.6).

(4.1) is not sharp, at least for N ≥ 3 and p going to ∞. Needless to say, this argument (and the related one for p → 1) is only a heuristic one, since these limits are not uniform with respect to the sets Ω.

The analogues of problem (1.1) in these extremal cases are the following min{C 2 (Ω) : |Ω| = c} and min{Λ 2 (Ω) : |Ω| = c}.

Once again, they both have (unique) solution given by any disjoint union of two equal balls: for the first one, the reader can see [START_REF] Parini | The second eigenvalue of the p-Laplacian as p goes to 1[END_REF]Proposition 3.14], while the second is easily derived thanks to the geometrical meaning of Λ 2 . In scaling invariant form, these rewrite as

|Ω| 1/N C 2 (Ω) ≥ 2 1/N N ω 1/N N and |Ω| 1/N Λ 2 (Ω) ≥ 2 1/N ω 1/N
N , and they both can be improved in a quantitative form, as it is proved in the following.

Theorem 5.2. Let Ω be an open subset of R N having finite measure. Then

(5.2) |Ω| 1/N C 2 (Ω) ≥ 2 1/N N ω 1/N N 1 + h N A 2 (Ω) N +1 ,
where the constant h N only depens on the dimension N . Moreover, for Λ 2 we have

(5.3) |Ω| 1/N Λ 2 (Ω) ≥ 2 1/N ω 1/N N 1 + 1 2N A 2 (Ω) .
Proof. To prove (5.2), we start defining

T Ω = t > 0 : there exist Ω 1 , Ω 2 ⊂ Ω disjoint
) |Ω t + | , P (Ω t -) |Ω t -| ≥ max{C 1 (Ω t + ), C 1 (Ω t -)},
where we used the straightforward estimate C 1 (E) ≤ P (E)/|E|, which is valid for every finite perimeter set E. Now, we introduce the following quantity

D Ω (t) := |Ω| 1/N max{C 1 (Ω t + ), C 1 (Ω t -)} 2 1/N N ω 1/N N -1,
and we proceed exactly as in Case A of the proof of Theorem 4.2: we only have to replace HKS(Ω) by D Ω (t) and the quantitative Faber-Krahn inequality by the following (sharp) quantitative Cheeger inequality (see [START_REF] Figalli | A note on Cheeger sets[END_REF]), (5.5)

|Ω| 1/N C 1 (Ω) ≥ N ω 1/N N 1 + γ N A(Ω) 2 ,
where γ N is a constant depending only on the dimension N . In this way, one arrives at the estimate

D Ω (t) ≥ h N A N +1 2 (Ω), for every t ∈ T Ω , that is |Ω| 1/N t 2 1/N N ω 1/N N -1 ≥ h N A N +1 2
(Ω), for every t ∈ T Ω , thanks to (5.4). Taking the infimum on T Ω on both sides and using the definition of second Cheeger constant, we eventually prove the thesis.

In order to prove (5.3), let us take two disjoint balls B(x 0 , r), B(x 1 , r) ⊂ Ω, whose common radius r is given by Λ 2 (Ω) = r -1 , and set for simplicity O 1 := B(x 0 , r) ∪ B(x 1 , r), then obviously we have

|Ω \ O 1 | = |Ω| -2 ω N r N .
Up to a translation, we can assume that x 0 = (M, 0, . . . , 0) and x 1 = (-M, 0, . . . , 0), for some M ≥ r, then for every t ≥ 1 we define the new centers x 0 (t) := (M + (t -1) r, 0, . . . , 0) and x 1 (t) := ((1t) r -M, 0, . . . , 0): observe that x i (1) = x i , i = 0, 1. Finally, we set

O t := B(x 0 (t), t r) ∪ B(x 1 (t), t
|Ω| -2 ω N r N = |Ω \ O 1 | ≥ |Ω \ O t 0 | ≥ 1 2 A 2 (Ω) |Ω|,
where in the last inequality we used that O t 0 is admissible for the problem defining A 2 (Ω). From the previous, we easily obtain

|Ω| r N ≥ 2 ω N (1 -1/2 A 2 (Ω))
, which finally gives (5.3), just by raising both members to the power 1/N , using the elementary inequality (1t) -1/N ≥ 1 + 1/N t and recalling that r = Λ 2 (Ω) -1 .

Sharpness of the estimates: examples and open problems

In the estimates of Theorem 4.2 and 5.2, we have shown that for every set, the relevant notion of deficit dominates a certain power κ of the asymmetry A 2 . If in addition to this, one could prove that for some sets converging to the optimal shape (i.e. a disjoint union of two equal balls), the deficit and A κ 2 have the same decay rate to 0, then these estimates would turn out to be sharp: we devote the last section to this interesting issue. 6.1. Quantitative Hong-Krahn-Szego inequality. Here, the question of sharpness is quite a delicate issue: first of all, observe that in contrast with the case of the Faber-Krahn inequality, the exponent of the asymmetry κ 2 blows up with N : then, one could automatically guess that κ 2 is not the sharp exponent. However, it has to be noticed that this dependence on N is directly inherited from the geometrical estimate (4.9), which is indeed sharp: let us fix a small parameter ε > 0 and consider the following set

Ω ε = {(x 1 , x ′ ) ∈ R N : (x 1 + 1 -ε) 2 + |x ′ | 2 < 1} ∪ {(x 1 , x ′ ) ∈ R N : (x 1 -1 + ε) 2 + |x ′ | 2 < 1},
which is just the union of two balls of radius 1, with an overlapping part whose area is of order ε (N +1)/2 . We set

Ω ε + = {(x 1 , x ′ ) ∈ Ω ε : x 1 > 0} and Ω ε -= {(x 1 , x ′ ) ∈ Ω ε : x 1 < 0},
and it is not difficult to see that A(Ω ε + ) = O(ε (N +1)/2 ), while on the contrary A 2 (Ω ε ) = O(ε) which means A 2 (Ω ε ) (N +1)/2 ≃ A(Ω ε ). i.e. both sides in (4.9) are asymptotically equivalent, as the area of the overlapping region goes to 0 (see [START_REF] Brasco | Sharp stability of some spectral inequalities[END_REF]Example 3.4], for more details on this example). And in fact one can use these sets Ω ε to show that the sharp exponent in (4.2) has to blow-up with the dimension. Also observe that in the proof of Theorem 4.2, the precise value of κ 1 plays no role, so the same proof actually gives (4.2) with κ 2 = (sharp exponent for (4.1)) × N + 1 2 .

Though we strongly suspect this κ 2 not to provide the right decay rate, currently we are not able to solve this issue, which seems to be quite a changelling one even for p = 2.

6.2. Second Cheeger constant. Also in this case, the exponent N + 1 in (5.2) seems not be sharp in the decay rate of the deficit. In order to shed some light on this fact, we estimate the deficit for C 2 of the same set Ω ε as before: first of all, thanks to the symmetries of Ω ε , it is not difficult to see that C 2 (Ω ε ) = C 1 (Ω ε + ) = C 1 (Ω ε -). Then we have

h N A 2 (Ω ε ) N +1 ≤ |Ω ε | 1/N C 2 (Ω ε ) 2 1/N N ω 1/N N -1 = |Ω ε | 1/N C 1 (Ω ε + ) 2 1/N N ω N -1 ≤ |Ω ε + | 1/N -1 P (Ω ε + ) N ω 1/N N -1,
so that the deficit of this inequality is controlled from above by the isoperimetric deficit of one of the two cut balls. We then estimate the right-hand side in the previous expression: observe that setting ϑ = arccos(1ε), we have for istance

P (Ω ε + ) = N ω N + ω N -1 (sin ϑ) N -1 -(N -1) sin ϑ 0 t N -2 √ 1 -t 2 d̺ ,
and

|Ω ε + | = ω N 1 -N 1 cos ϑ (1 -t 2 )
N -1 2 dt , then

P (Ω ε + ) ≃ N ω N - N -1 N + 1 ω N -1 2 ε N +1 2 and |Ω ε + | 1/N -1 ≃ ω 1-N N N 1 + N -1 N + 1 ε N +1 2
, so that in the end (6.1)

C 1 A 2 (Ω ε ) N +1 ≤ |Ω ε | 1/N C 2 (Ω ε ) -2 1/N N ω 1/N N ≤ C 2 A 2 (Ω ε ) N +1
2 , where we used that A 2 (Ω ε ) = O(ε). Notice that this estimate implies in particular that, also in this case, the sharp exponent is dimension-dependent and it blows up as N goes to ∞.

We point out that the previous computations gives the correct decay rate to 0 of the quantity C 2 (Ω ε ) -C 2 (B), which is O(ε (N +1)/2 ) = O(A 2 (Ω ε ) (N +1)/2 ). Indeed, from the right-hand side of (6.1) we can promptly infer that

C 2 (Ω ε ) = C 1 (Ω ε + ) ≤ N + c ε N +1 2 = C 1 (B) + c ε N +1
2 . Now assume that C 1 (Ω ε + ) ≤ C 1 (B) + ω(ε) for some modulus of continuity ω such that ω(ε) = o(ε (N +1)/2 ) as ε goes to 0, in this case we would obtain

0 ≤ |Ω ε | 1/N C 1 (Ω ε + ) -2 1/N N ω N ≤ -K ε N +1 2 ,
for some constant K > 0 independent of ε. This gives a contradiction, thus proving that

C 2 (Ω ε ) -C 2 (B) ≃ ε N +1 2 .
6.3. Second eigenvalue of -∆ ∞ . On the contrary, it is not difficult to see that the quantitative estimate (5.3) is sharp: still taking the set Ω ε as before, we observe that

Λ 2 (Ω ε ) = Λ 1 (Ω ε + ) = 2 2 -ε ≃ 1 + ε 2 and |Ω ε | 1/N ≃ ω 1/N N 1 - 1 N + 1 ε N +1 2
, while A 2 (Ω ε ) = O(ε) as already observed. Then

|Ω| 1/N Λ 2 (Ω) -ω 1/N N ≃ A 2 (Ω),
proving the sharpness of (5.3). We remark that in this case the sharp exponent does not depend on the dimension, in contrast with the cases p ∈ [1, ∞).

(1. 2 )

 2 |Ω| 2/N λ 2 (Ω) ≥ 2 2/N ω 2/N N λ 1 (B), with equality if and only if Ω is a disjoint union of two equal balls. Here and in what follows, B will always denote a N -dimensional ball of radius 1 and ω N := |B|. Observe that for Ω = B 1 ∪ B 2 , with B 1 and B 2 disjoint balls having |B 1 | = |B 2 |, the first eigenvalue has multiplicity 2, i.e. λ 1 (Ω) = λ 2 (Ω) and these are equal to the first eigenvalue of one of the two balls.

Lemma 4 . 3 .

 43 Let Ω ⊂ R N be an open bounded set. With the notations of Lemma 3.1, we have

  and s.t. max

	i=1,2 It is not difficult to see that if Ω is open, then T Ω = ∅, since Ω contains at least two disjoint small P (Ω i ) ≤ t . |Ω i | balls, which are in particular two sets with positive measure and finite perimeter. Then let us pick
	up a t ∈ T Ω : correspondingly, there exist Ω t + , Ω t -⊂ Ω disjoint and such that (5.4) t ≥ max P (Ω t +

  r), t ≥ 1 i.e. for every t ≥ 1 this is a disjoint union of two balls of radius t r and moreover O t ⊂ O s if t < s. The latter fact implies that the function t → |Ω ∩ O t | is increasing, thus t → |Ω \ O t | is decreasing. We exploit this fact by taking t 0 > 1 such that |O t 0 | = |Ω|: then we have

Acknowledgements. The authors are grateful to Cristina Trombetti for having stimulated the writing of this paper, to Mark S. Ashbaugh for some historical remarks on the Hong-Krahn-Szego inequality and to Peter Takáč for having pointed out to them reference [12]. The first author has been supported by the ERC Advanced Grant n. 226234.