 Abstract -This paper presents an approach for quick calculation of steady-state and transient performances of an axial-field eddy-current coupling. Based on a two-dimensional (2D) approximation of the magnetic field distribution, a simple analytical expression for the transmitted torque is first developed. This expression is valid for low slip values which correspond to the normal working area of such devices (high efficiency). The proposed torque formula is then used to study the steady-state (constant speed operations) and the transient performances of the coupling (small variations of the slip speed). The results are compared with those obtained from 3-D finite elements simulations and tests. It is shown that good agreements are obtained.

Index Terms--Analytical model, magnetic transmission, eddy-current, torque, transient performance.

NOMENCLATURE R 1

Inner radius of the magnets R 2

Outer radius of the magnets R m

Mean radius of the magnets (R m =(R 1 +R 2 )/2) R 0

Inner radius of the copper R [START_REF] Yonnet | Analytical calculation of permanent magnet couplings[END_REF] Outer radius of the copper L Radial length of the magnets (L=R AGNETIC couplings can transmit a torque from a primary driver to a follower without mechanical contact. Therefore, they are well suited for use in isolated systems such as vacuums, high-pressure or cryogenic vessels. Compared to mechanical couplings, they present great advantages such as self protection against overload condition and great tolerance to shaft misalignment [START_REF] Hornreich | Optimal design of synchronous torque couplers[END_REF]- [START_REF] Lubin | Simple analytical expressions for the force and torque of axial magnetic couplings[END_REF].

There are two main types of magnetic coupling: synchronous [START_REF] Hornreich | Optimal design of synchronous torque couplers[END_REF]- [START_REF] Krotsch | Radial forces in external rotor permanent magnet synchronous motors with non-overlapping windings[END_REF] and induction type (eddy-current) [START_REF] Davies | An experimental and theoretical study of eddy-current couplings and brakes[END_REF]- [START_REF] Mohammadi | Design Optimization of Double-Sided Permanent-Magnet Radial-Flux Eddy-Current Couplers[END_REF], both with radial or axial flux configuration. In this paper, we focus on single-sided axial field eddy-current coupling as shown in Fig. 1. It consists of two disks: one is composed of rare earth permanent magnets (axially magnetized), and the other one is equipped with a conducting plate (generally copper). To improve the performances, the flux is closed by two back iron yokes.

The working principle of eddy-current couplings is well known [START_REF] Davies | An experimental and theoretical study of eddy-current couplings and brakes[END_REF]. The eddy currents, which are induced in the conducting plate by the slip speed (= 1 - 2 ), interact with the magnetic field and generate a torque. The normal working range area of such devices corresponds to low slip values where the temperature rise, due to induced-current, is limited. In this condition, the torque is proportional to the slip [START_REF] Mohammadi | Non-linear modeling of eddy-current couplers[END_REF]. As shown in [START_REF] Davies | An experimental and theoretical study of eddy-current couplings and brakes[END_REF], the efficiency of eddy-current couplings is given by  =1-s where s is the slip (s = / 1 ). Fig. 2 shows the prototype of the studied eddy-current coupling placed on the test bench. The physical and geometrical parameters are given in Table I. These parameters have been optimized but the optimization procedure is not developed in this paper. Design optimization of eddy-current couplings can be found in [START_REF] Mouton | Modeling and Optimal Design of an Eddy Current Coupling for Slip-Synchronous Permanent Magnet Wind Generators[END_REF], [START_REF] Potgieter | Optimum design and technology evaluation of slip permanent magnet generators for wind energy applications[END_REF] and [START_REF] Mohammadi | Design Optimization of Double-Sided Permanent-Magnet Radial-Flux Eddy-Current Couplers[END_REF]. With these parameters and considering c = 3 mm for the air-gap, we obtain a torque of about 10 Nm for a slip speed of 150 rpm [START_REF] Niguchi | cogging torque analysis of magnetic gear[END_REF]. A slip speed of 150 rpm corresponds to an efficiency of 95% (=1-s) when the eddy-current coupling is used with a motor running at a nominal speed of 3000 rpm [START_REF] Davies | An experimental and theoretical study of eddy-current couplings and brakes[END_REF]. Eddy-current couplings have been studied for a long time and many papers can be found in the literature [START_REF] Davies | An experimental and theoretical study of eddy-current couplings and brakes[END_REF]- [START_REF] Mohammadi | Analytical modeling and analysis of axial-flux interior permanentmagnet couplers[END_REF].

The study of such devices is usually based on a twodimensional approximation for the magnetic field distribution. The induced currents in the conducting plate are computed by solving a 2D diffusion equation leading to a complex expression for the transmitted torque [START_REF] Canova | Design of axial eddy-current couplers[END_REF]- [START_REF] Wang | A general analytical model of permanent magnet eddy current couplings[END_REF]. More recently [START_REF] Mohammadi | Double-Sided permanent-magnet radial-flux eddy-current couplers: three-dimensional analytical modeling, static and transient study, and sensitivity analysis[END_REF], [START_REF] Mohammadi | Non-linear modeling of eddy-current couplers[END_REF], the magnetic equivalent circuit (MEC) method has been used to predict the performance of eddy-current couplings by taking into account the magnetic saturation effects and all material properties. With the MEC method, it is possible to study eddy-current couplings with complex geometries [START_REF] Mohammadi | Analytical modeling and analysis of axial-flux interior permanentmagnet couplers[END_REF]. In such cases, analytical models have significant issues. The main drawback of the MEC method is that the flux paths must be a priori known in order to define the reluctance expressions which appear in the torque expression [START_REF] Mohammadi | Double-Sided permanent-magnet radial-flux eddy-current couplers: three-dimensional analytical modeling, static and transient study, and sensitivity analysis[END_REF], [START_REF] Mohammadi | Non-linear modeling of eddy-current couplers[END_REF], and [START_REF] Mohammadi | Analytical modeling and analysis of axial-flux interior permanentmagnet couplers[END_REF].

For an engineering purpose, it is important to have a simple and accurate torque formula in order to quickly study the effects of geometrical parameters on the coupling performances. As indicated previously, the working area of eddy-current couplings corresponds to low slip values. Therefore, analysis of such device can be limited to low slip values, which simplify greatly the analytical development.

The goal of this paper is to develop the simplest model which gives acceptable results in the design area (low slips). To achieve this outcome, we have to do some assumptions that will be justified in the next developments. The 3D effects (actual paths of the induced currents) are taken into account by using the Russell correction factor [START_REF] Russell | Eddy currents and wall losses in screened-rotor induction motors[END_REF]. The proposed torque formula, which depends directly on the physical and geometrical parameters, is then used for the analysis of both steady-state and transient operations. It should be noted that transient performance analysis of eddy-current couplings is rarely treated in the literature [START_REF] Davis | Transient performance of eddy-current couplings[END_REF]- [START_REF] Tonoli | Dynamic characteristics of eddy current dampers and couplers[END_REF], and very few experimental results are given [START_REF] Tonoli | Dynamic modeling and experimental validation of eddy current dampers and couplers[END_REF], [START_REF] Bouwer | Modelling and dynamic performance of a direct-drive direct-grid slip permanent magnet wind generator[END_REF].

II. TORQUE EXPRESSION FOR LOW SLIP VALUES

Axial flux eddy-current coupling is an inherent 3D geometry from the modeling point of view. However, in order to simplify the analysis, the 3D problem of Fig. 1 can be reduced to a 2D problem by introducing a cutting plane at the average radius of the magnets R m =(R 1 +R 2 )/2. Fig. 3 shows the resulting 2D model in Cartesian coordinate (y, z) with an infinite length along the x-direction. The 3D effects (actual path of the induced currents) cannot be overlooked for such device and will be taken into account by using the classical Russell correction factor [START_REF] Russell | Eddy currents and wall losses in screened-rotor induction motors[END_REF]. Moreover, for reasons of simplicity, we adopt the following assumptions:

1) The reaction field due to induced currents is neglected because of the low slips assumption,

2) The iron-yokes are considered with infinite magnetic permeability and zero conductivity,

3) Only the first harmonic is considered for the magnetic field distribution, 4) Magnets are axially magnetized with relative recoil permeability 1  r  . As the iron-yokes present an infinite permeability, the tangential component of the magnetic field is zero at z=0 and z=b+c+d (boundary conditions). This limits the number of region where the magnetic field has to be computed. This hypothesis can be regarded as accurate because the thicknesses a and e of the back-iron (see Fig. 3) must be determined to avoid magnetic saturation.

As it can be seen in Fig. 3, the whole domain of the field problem is divided into two regions: the PMs region (region I) and the air-gap with copper region (region II). As the reaction field is neglected, the air-gap and copper can be connected because we have the same governing equation for this region, i.e. Laplace's equation.

A magnetic vector potential formulation (A i =A i (y,z)e x ) is used to solve the problem. From the Maxwell equations and considering the adopted assumptions, we have to solve the following equations in each region where  0 is the vacuum permeability and M z is the axial magnetization of the magnets. The distribution of the axial magnetization M z (y) is plotted in Fig. 4. The axial magnetization can be expressed in Fourier's series. It has been shown in [START_REF] Wang | A general analytical model of permanent magnet eddy current couplings[END_REF] that the first harmonic component plays a dominant role in the torque transmission for eddy-current couplings (around 98% of the total torque). So, in the following development, we will consider only the first harmonic of the magnetization.

  0 4 ( ) sin sin 2 r z B M y y        with m pR   (2)
According to the symmetry of the magnetic field distribution and considering the classical boundary and interface conditions for such problem, solutions of ( 1) is

          sinh ( ) cosh ( , ) 1 cos sinh I c d z A y z K y b c d           (3)             sinh cosh ( ) ( , ) cos sinh II b z b c d A y z K y b c d         (4) with   4 sin 2 rm BR K p    ( 5
)
The induced current density in the conducting plate can be deduced from Lorentz's equation

II II m A R y         x J v B e ( 6
)
where  =  1 - 2 is the slip speed (rd/s) between the two discs. The transmitted torque can be evaluate by the power losses dissipated in the conducting plate From ( 4), ( 6) and ( 7), we can derive a closed-form expression for the electromagnetic torque which depends directly on the physical and geometrical parameters

  12 TT T K K       ( 8
)
where K T is given by Fig. 4. Magnetization distribution along y-direction (region I)

2 3 2 2 2 8 sin 2 sinh ( ) sinh(2 ) 1 2 sinh ( ( )) T r m K B R Ld bd d b c d                (9)
As expected for low slip values, the torque is proportional to the slip speed  (rd/s). The expression has been arranged in order to isolate its dimensional part. It has been shown in [START_REF] Wang | A general analytical model of permanent magnet eddy current couplings[END_REF], [START_REF] Mohammadi | Double-Sided permanent-magnet radial-flux eddy-current couplers: three-dimensional analytical modeling, static and transient study, and sensitivity analysis[END_REF] and [START_REF] Mouton | Modeling and Optimal Design of an Eddy Current Coupling for Slip-Synchronous Permanent Magnet Wind Generators[END_REF] that the 3D endeffects cannot be neglected for such devices. Fig. 5 shows the actual paths of the induced currents in the conducting plate (3D effects).The induced currents flow along a closed path and not only in the x-direction as it was obtained with the 2D model [START_REF] Ravaud | Analytical design of permanent magnet radial couplings[END_REF]. To take into account the end-effects, an efficient correction factor [START_REF] Charpentier | Optimal design of cylindrical air-gap synchronous permanent magnet couplings[END_REF] has been given by Russel and Norsworthy [START_REF] Russell | Eddy currents and wall losses in screened-rotor induction motors[END_REF] (2 /

) tanh( / 2)

1 1 tanh( / 2) tanh( / 2) Russel LL k LL      ( 10
)
The torque expression with the 3D effects can be finally expressed as

' T TK   with ' T Russel T K k K  ( 11
)
where K' T is the torque coefficient. Fig. 6 shows the variation of K' T versus the air-gap length (the other parameters are those given in table I). As expected, K' T decreases with the air-gap length. I).

III. STEADY-STATE ANALYSIS

In this section, we use the torque formula [START_REF] Niguchi | cogging torque analysis of magnetic gear[END_REF] to study the steady-state performances of the coupling (constant speed operation). In order to show the effectiveness of the proposed model, the results obtained with the analytical formula (11) are compared with those obtained from 3D nonlinear finite element analysis (FEA). The finite element simulations (COMSOL Multiphysics®) are conducted by considering the reaction field due to the induced currents and the real geometry of the coupling (i.e. cylindrical structure as shown in Fig. 1). Moreover, the nonlinear magnetic property for the back-irons (B-H curve) has been considered in the finite-element simulations and is shown in Fig. 7.

A. Torque-Speed Characteristic

The torque-speed characteristic of the studied coupling is shown in Fig. 8. As indicated previously, only the low slip values are considered here (high efficiency area). The geometrical parameters are the ones given in Table I with c=3mm. Fig. 8 shows that the torque is well predicted by the analytical formula if we consider the Russel correction factor. In this case, the deviations between the proposed formula and the 3D nonlinear FEA is always below 10%.

B. Influence of Geometrical Parameters

The torque formula is now used to study the influence of several geometrical parameters (copper thickness, polepairs number, and the radial length of the copper). Fig. 10 shows the torque versus copper thickness. We have fixed the air-gap length c=3mm and the slip speed at 150 rpm (the other parameters are those given in Table I). We can observe that an optimal value exists for the copper thickness (T ≈ 10 Nm for d ≈ 5mm). This optimal value is well predicted by the analytical formula (11) when compared to 3D FEA.

The torque versus the pole-pairs number is shown in Fig. 11. For the analysis, we have considered c=3mm, the other geometrical parameters are kept constant and are given in Table I. Once again, we can observe that an optimum value exists for the pole-pairs number (p≃4). This optimum value is well predicted by using the torque formula. The error is never greater than 15% compared to 3D FEA.

As indicated in the nomenclature,  is the ratio of the radial copper end-lengths to the radial magnet length (active length): =0 corresponds to the same radial dimensions for the copper and the magnets (H=L) whereas =1 corresponds to H=2L. Fig. 12 shows the torque versus the ratio  obtained with the analytical formula and with 3D FEA. We can observe that the torque increases with the copper end-lengths and then becomes constant for  ≥ 1 (an improvement of about 75% is obtained with =1 compared to =0). For the studied coupling (Table I), we have therefore chosen =1, that corresponds to 30 mm additional length for the copper compared to the active length of the magnets (15 mm on each side as indicated in Table I). Fig. 12 clearly shows the ability of the torque formula to take into account the ratio λ. The error between analytical model and 3D FEA is not greater than 17% (for λ=0), and we can observe that the trend is good.

C. Axial Force Expression

The main drawback of axial-field couplers is the large axial force between the two discs [START_REF] Lubin | Simple analytical expressions for the force and torque of axial magnetic couplings[END_REF]. This force must be known because it affects the bearing losses of the drive. For eddy-current couplings, the attraction force is maximal when the reaction field is null (i.e. slip speed equal to zero). The proposed analytical model can then be used to derive a closed-form expression for the axial force. By using the Maxwell stress tensor, the axial force expression is

  2 2 2 2 2 21 2 0 4 sinh ( ) sin 2 sinh ( ( )) r b F B R R b c d          (12)
By using [START_REF] Lubin | Experimental and theoretical analyses of axial magnetic coupling under steady-state and transient operation[END_REF], we obtain F = 581N with the geometrical parameters given in Table I and considering c=3mm. The 3D FEA gives F = 515N, which corresponds to an error of about 13% on the axial force prediction.

D. Experimental validation

For the experimental validation, we have manufactured an axial-field eddy-current coupling prototype as shown in Fig. 2. The geometrical parameters of the prototype are given in Table I. The magnetic coupling is inserted between a DC motor and a load as shown in Fig. 14 (DC motor of 3kW, 3000 rpm). For the test, one shaft is locked (copper side i.e. Ω 2 = 0 in Fig. 14), and the other shaft can rotate (magnets side). The speed is controlled by changing the voltage applied to the armature of the DC motor. The slip speed is measured with an encoder as shown in Fig. 14. The torque is obtained by the armature current measurement (the torque constant of the DC motor is 1,35Nm/A).

Fig. 13 shows a comparison between the analytically predicted torque-slip speed characteristics and the measured data. As shown in Fig. 13, three values for the air-gap length were considered (c=3mm, c=5mm, and c=7mm). It can be noticed that the experimental measurements are in good agreement with the one obtained with the proposed torque formula. The deviations between the analytical formula and the test results are not greater than 10%.

However, we can observe in Fig. 13 that the error between the analytical predictions and tests becomes larger when the speed increases. This difference is due to the wellknown reaction field which can no longer be overlooked for a speed greater than 200 rpm (for the studied coupling). In order to predict accurately the torque-speed characteristic for larger slip values, a more complex model which takes into account the reaction field effect must be used [START_REF] Canova | Design of axial eddy-current couplers[END_REF]- [START_REF] Wang | A general analytical model of permanent magnet eddy current couplings[END_REF], but this is out of the scope of the present paper.

IV. TRANSIENT ANALYSIS

The transient analysis of the magnetic coupling is obtained from the equation of motion for rotating rigid bodies. Fig. 14 shows the scheme of the test bench. The magnetic coupling is placed between two DC motors. One motor is used as the prime mover (attached to the magnets side of the coupling) and the other DC motor is used to apply a load torque (DC generator). The prime mover rotates at Ω 1 and the load runs at Ω 2 . In steady state condition, the difference between the two speeds is a function of the load torque. Two encoders (4096 pulses/revolution) have been placed on the test bench to measure the speed on both sides of the coupling during the transient. From Fig. 14, the motion equations are given by

1 1 1 1 2 2 2 2 DC load d J B T T dt d J B T T dt           (13)
where J 1 and B 1 denote the moment of inertia and the viscous damping coefficient of the DC motor, J 2 and B 2 denote the moment of inertia and the viscous damping coefficient of the load. T DC is the DC motor torque and T load is the external load torque.

Here, we consider only small variations of the slip speed. As shown in [START_REF] Amati | Dynamic behaviour of torsional eddy-current dampers: sensitivity of the design parameters[END_REF]- [START_REF] Tonoli | Dynamic modeling and experimental validation of eddy current dampers and couplers[END_REF], eddy-current couplings can be modeled as a first order transfer function (small signal model), where K' T is the DC gain and τ c is the time constant (s is the Laplace domain variable)

' () ( ) 1 T c K Ts ss    with 1 c TMAX p    ( 14
)
Ω TMAX is the slip speed (rad/s) corresponding to the maximum torque of the torque to speed characteristic (Ω TMAX can be obtained by using a model which take into account the reaction field [START_REF] Canova | Analytical modeling of rotating eddycurrent couplers[END_REF]). For the studied coupler (Table I), we obtain Ω TMAX = 84 rad/s which gives τ c = 2.3ms. This time constant should be compared to the one given in ( 17) which takes into account the mechanical part of the system (τ=38.5ms for c=3mm). For the studied coupler, the ratio τ/ τ c ≃ 16.7. Therefore, τ is the dominant time constant of the system and τ c can be neglected in the transient analysis. In the next developments, we assume that the torque transmitted by the magnetic coupling is proportional to the slip speed (the coupler is seen as a purely viscous torque). By considering [START_REF] Niguchi | cogging torque analysis of magnetic gear[END_REF], the equation of motion ( 13) can be re-written as

    '' 1 1 1 1 2 '' 2 2 2 2 1 T DC T T T load d J B K T K dt d J B K K T dt             (15)
In the following sub-sections, we study three types of transient response. The first test concerns the speed response to a step input on the torque of the DC motor when the copper face of the coupling is locked (brake operation). In the second test, we study the transient performance during start-up (under no-load condition) for different values of the air-gap length. The last test concerns the speed responses on both sides of the coupling when an abrupt load variation is applied.

A. Speed Response to a Step Input on T DC

In order to study the transient behavior of the coupling, a first test consists to block one part of the magnetic coupling (the load part is locked i.e. Ω 2 = 0 in Fig. 14), and to apply a step input on T DC . In this case, the motion equation ( 15) can be simplified to

  ' 1 1 1 1 T DC d J B K T dt      (16)
For this test, the eddy-current coupling acts as a brake and is equivalent to a supplementary viscous damping coefficient which is equal to the torque coefficient K' T given by [START_REF] Niguchi | cogging torque analysis of magnetic gear[END_REF]. As indicated in ( 9) and [START_REF] Charpentier | Optimal design of cylindrical air-gap synchronous permanent magnet couplings[END_REF], the torque coefficient depends directly on the coupling geometrical parameters and more particularly on the air-gap value (Fig. 6). For an air-gap c = 3mm, we obtain K' T = 0.65Nm.s/rad as shown in Fig. 6.

Equation ( 16) is a well-known first order linear differential equation. The speed response to a step input on T DC is given by [START_REF] Canova | Analytical modeling of rotating eddycurrent couplers[END_REF]. We can define the time constant ( ) of the system and the final speed value ( 1 ). As the viscous damping coefficient of the DC motor B 1 is much lower than K'

T , it can be neglected (B 1 = 0.002Nm.s/rad whereas K' T is greater than 0.1 Nm.s/rad as shown in Fig. 6). The viscous damping coefficient B 1 has been measured by the no-load test. The moment of inertia J 1 in ( 17) is equal to 0.025kg.m 2 and has been measured by the classical freewheeling operations. 1 for simulation and test. In all cases, the error is less than 9%. As expected, the final speed values given in Table II correspond to the ones obtained during the steady-state analysis (Fig. 13) for a torque of 5Nm.

B. Transient Performance During Start-Up

For this test, the two sides of the coupling can rotate. The principle is to start the DC motor by connecting its armature circuit to a DC source of low voltage. The DC motor starts from standstill to 400 rpm in steady-state condition. The test is done under no-load condition (T load =0 in (15)).

Figs. 16-17 compare the speed response obtained with the analytical model and with experimental measurements for two values of the air-gap length (c = 3mm, and c = 7mm). The simulation results have been obtained by using [START_REF] Nehl | Nonlinear two-dimensional finite element modeling of permanent magnet eddy current couplings and brakes[END_REF]. We can observe in Fig. 16 and Fig. 17 that a delay time exists for the speed between the two rotors of the magnetic coupling. This delay time increases when the airgap length increases. This is due to the decrease of the torque coefficient K' T with the air-gap length, as shown in Fig. 6. This delay time has to be taken into account for servomechanism application. Because of the no-load condition, the speed value in steady state is almost identical on both sides of the coupling (a little speed difference can be observed in Fig. 17 for a large air-gap i.e. 7mm).

It is worth noting that no oscillation appears on the speed responses during the transient, as it was for synchronous coupling with permanent magnets on both sides [START_REF] Lubin | Experimental and theoretical analyses of axial magnetic coupling under steady-state and transient operation[END_REF]. An advantage of eddy-current couplings compared to synchronous couplings is that it is free from torque ripple as shown in [START_REF] Potgieter | Optimum design and technology evaluation of slip permanent magnet generators for wind energy applications[END_REF], and can be used as vibration isolation system.

Fig. 18 shows the torque transmitted by the coupler during the transient corresponding to the speed variation given in Fig. 17. We can observe that the torque reaches a maximal value of around 2Nm (this point corresponds to the maximum slip speed during the transient) and then is almost equal to zero in steady-state (no-load condition).

C. Transient Performance with Load variation

The transient performance of the eddy-current coupling for a sudden application of load torque is now studied. Before the load torque is applied, the DC motor is in steady-state condition and runs with a speed of 750 rpm under no-load condition. For the abrupt application of the load torque, we use a DC generator by connecting a resistance to its terminal armature winding (the mechanical parameters of the DC generator are B 2 =0.0015Nm.s/rad and J 2 =0.01kg.m 2 ). Once again, two values of the air-gap length are considered for this test. Figs. [START_REF] Choi | Analytical magnetic torque calculations and experimental testing of radial flux permanent magnet-type eddy current brakes[END_REF][START_REF] Mohammadi | Double-Sided permanent-magnet radial-flux eddy-current couplers: three-dimensional analytical modeling, static and transient study, and sensitivity analysis[END_REF] show the speed responses (simulation and experimental results) for a sudden load torque of 3.2Nm. After a transient of around 0.4s, the two rotors are once again in steady state but with a slip speed which depends on the air-gap value. As expected, the relative speed between the two rotors increases with the air-gap value.

Overall, we can observe closed agreements between the analytical and experimental results, which confirm the validity of the proposed model. However, it is important to keep in mind that the torque formula is only valid for low slip values and the proposed model can only predict the transient performance for small variations of the slip speed and under certain conditions indicated previously.

V. CONCLUSION

In this paper, we have developed a simple torque formula that can be easily used for a first evaluation of eddy-current couplings performances, for both steady-state and transient operations. The proposed model is valid for low slip values but it corresponds to the normal working area of eddy-current couplings (limited temperature and high efficiency). Comparisons with 3-D finite-element simulations have shown that the torque formula is able to predict the influence of the geometrical parameters with a good precision. Therefore, it can be used as an efficient tool in the design process of eddy-current couplings.

A prototype has been built and the comparison between analytical results and tests has shown the precision of the model. An error of less than 10% has been obtained on the torque-speed characteristic. Through detailed transient analysis, we have shown that the magnetic coupling causes time delays between the two rotors during start-up. This time delay depends on the air-gap value and must be taken into account for servomechanism applications and every transient behavior.

  2 -R 1 ) H Radial length of the copper (H=R 3 -R 0 )  Copper end-lengths to magnets length  = (H-L)/L a Thickness of the back-iron (magnets side) b Magnets thickness c Air-gap thickness d Copper thickness e Thickness of the back-iron (copper side) α PMs pole-arc to pole-pitch ratio p Pole-pairs number B r Remanence of the permanent magnets  Conductivity of the conducting plate (copper)   = p/R m I. INTRODUCTION

Fig. 1 .Fig. 2 .

 12 Fig. 1. Geometry of the studied axial-field eddy-current coupling (p = 5) with its geometrical parameters.

Fig. 3 .

 3 Fig. 3. 2D model of the eddy-current coupling at the mean radius of the magnets Rm =(R1+R2)/2.

Fig. 5 .

 5 Fig. 5. Actual paths for the induced-currents (3D effects)

Fig. 7 .

 7 Fig. 7. B(H) curve of the back-irons ( M-27 steel) used in 3D FEA

Fig. 8 .

 8 Fig. 8. Torque-speed characteristic for the studied coupling (with c=3mm): comparison between analytical results and 3D FEA.The flux density (z-component) and the eddy-current distributions at the surface of the copper obtained with 3D FEA are respectively shown in Fig.9aand Fig.9b. These results correspond to a slip speed of 150 rpm. The maximum value of the flux density on the surface of the copper is around 0.55T (with alternate polarity).

Fig. 9 .

 9 Fig. 9. 3D FEM analysis: (a) flux density distribution (z-component) on the surface of the copper, (b) eddy-current distribution at 150 rpm.

Fig. 10 .

 10 Fig. 10. Torque versus copper thickness at 150 rpm for c = 3mm and p=5.

Fig. 11 .

 11 Fig. 11. Torque versus number of pole pairs at 150 rpm for c = 3mm.

Fig. 12 .

 12 Fig. 12. Torque versus copper end-length at 150 rpm for c=3mm.

Fig. 13 .

 13 Fig. 13. Measured torque-slip speed characteristics for three values of the air-gap length (c=3mm, c=5mm, and c=7mm).

Figs. 15

 15 Figs. 15 show the speed responses  1 (simulations and tests) when we apply a step in load torque DC =5 Nm at t = 0s. The electromagnetic torque of the DC motor is controlled by the armature current regulation (closed loop control). As shown in Fig. 15(a), a step input T DC =5 Nm corresponds to a step input on the DC armature current of 3.7A (torque constant of 1.35Nm/A). Three values of the air-gap length are considered for this test (3mm, 5mm and 7mm). Fig. 15(b) shows that the transient behavior is well predicted by the analytical model compared to the experimental results.Table II compares the time constant  and the final speed value 

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Scheme of the test bench for the transient analysis

Fig. 16 .Fig. 17 .

 1617 Fig. 16. Speed response during start-up for c = 3mm: (a) experimental result, (b) simulation result.

Fig. 18 .Fig. 19 .

 1819 Fig. 19. Speed response to a sudden load torque of 3.2Nm for c = 3mm: (a) experimental result, (b) simulation result.

Fig. 20 .

 20 Fig. 20. Speed response to a sudden load torque of 3.2Nm for c = 7mm: (a) experimental result, (b) simulation result.

  Outer radius of the magnets Inner radius of the conducting plate Outer radius of the conducting plate Thickness of the back-iron (magnets side)

		TABLE I
	PARAMETERS OF THE STUDIED EDDY-CURRENT COUPLING
	Symbol	Quantity	value
	R1	Inner radius of the magnets	30 mm
	R2		60 mm
	R0		15 mm
	R3		75 mm
	a		10 mm
	b	Magnets thickness	10 mm
	c	Air-gap length	variable
	d	Conducting plate thickness	5 mm
	e	Thickness of the back-iron (copper side)	8 mm
	α	PMs pole-arc to pole-pitch ratio	0.9

p Pole-pairs number 5 Br Remanence of the permanent magnets (NdFeB) 1.25 T  Conductivity of the conducting plate (copper) 57 MS/m

TABLE II COMPARISON

 II BETWEEN TEST AND SIMULATION (TIME CONSTANT AND FINAL SPEED VALUE)


Pierre Bousquet
email: pierre.bousquet@univ-amu.fr

AND Lorenzo Brasco
email: lorenzo.brasco@univ-amu.fr

Vesa Julin
email: vesa.julin@jyu.fi

LIPSCHITZ REGULARITY FOR LOCAL MINIMIZERS OF SOME WIDELY DEGENERATE PROBLEMS

Keywords: Mathematics Subject Classification. 35J70, 35B65, 49K20. Key words and phrases. Degenerate elliptic equations, Anisotropic problems, Lipschitz regularity

We consider local minimizers of the functional N i=1 ˆ(|ux i |δi) p + dx + ˆf u dx, where δ1, . . . , δN ≥ 0 and ( • )+ stands for the positive part. Under suitable assumptions on f , we prove that local minimizers are Lipschitz continuous functions if N = 2 and p ≥ 2, or if N ≥ 2 and p ≥ 4.

F(u; Ω ′ ) = N i=1 ˆΩ′ (|u x i | -δ i ) p + p dx + ˆΩ′ f u dx, u ∈ W 1,p loc (Ω), Ω ′ ⋐ Ω.
Here Ω ⊂ R N is an open set, 2 ≤ p < ∞, δ i ≥ 0, ( • ) + stands for the positive part and f ∈ L p ′ loc (Ω) where p ′ = p/(p -1). This functional F stands for a model case of a more general class of problems, with specific growth and monotonicity assumptions. For the sake of clarity, the results in this paper are only stated for F. However, their proofs can be easily adapted to embrace general functionals having a similar structure.

The functional F naturally arises in problems of Optimal Transport with congestion and anisotropic effects, see for example [START_REF] Brasco | Congested traffic equilibria and degenerate anisotropic PDEs[END_REF][START_REF] Brasco | On certain anisotropic elliptic equations arising in congested optimal transport: local gradient bounds[END_REF] for some motivations. These two papers contained among others some regularity results for local minimizers of (1.1). For instance [START_REF] Brasco | On certain anisotropic elliptic equations arising in congested optimal transport: local gradient bounds[END_REF]Main Theorem] proves that if f ∈ L ∞ loc (Ω), then u is "almost Lipschitz", i.e. u ∈ W 1,r loc (Ω) for every r ≥ 1. On the other hand, in [START_REF] Brasco | Congested traffic equilibria and degenerate anisotropic PDEs[END_REF] it is proved that if f ∈ W 1,p ′ loc (Ω), then

(1.2) (|u x i | -δ i ) p 2 + u x i |u x i | ∈ W 1,2 loc (Ω), i = 1, . . . , N.
However, it must be mentioned that to the best of our knowledge, Lipschitz regularity of local minimizers is still unknown. More surprisingly, even the case δ 1 = • • • = δ N = 0 does not seem to be fully understood. Only very recently some results have been obtained in this case, see [START_REF] Birindelli | Lipschitz regularity for solutions of the pseudo p-Laplace Poisson equation[END_REF][START_REF] Demengel | Lipschitz interior regularity for the viscosity and weak solutions of the pseudo p-Laplacian[END_REF].

Observe that local minimizers of (1.1) are local weak solutions of the anisotropic degenerate equation

(1.3) N i=1 (|u x i | -δ i ) p-1 + u x i |u x i | x i = f,
which reduces to the Poisson equation for the so-called pseudo p-Laplacian when δ 1 = • • • = δ N = 0, i.e.

(1.4)

N i=1 |u x i | p-2 u x i x i = f.
The terminology "pseudo p-Laplacian" appears in [START_REF] Belloni | The pseudo p-Laplace eigenvalue problem and viscosity solution as p → ∞[END_REF]. We just point out that such an operator already appeared in J.-L. Lions's monograph [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF], where existence issues for solutions to evolutions equations are tackled.

In order to neatly explain the difficulty of the problem, we now recall some class of functionals for which the Lipschitz property for local minimizers is known to be true. The first one is given by (1.5) ˆG(∇u) dx, with G enjoying a p-Laplacian-type structure at infinity. This means that there exist c, C > 0 and m ≥ 0 such that G verifies the ellipticity condition We refer the reader to [START_REF] Brasco | Global L ∞ gradient estimates for solutions to a certain degenerate elliptic equation[END_REF][START_REF] Celada | Existence and regularity of minimizers of nonconvex integrals with pq growth[END_REF][START_REF] Chipot | Linearization at infinity and Lipschitz estimates for certain problems in the calculus of variations[END_REF][START_REF] Esposito | On the Lipschitz regularity for certain elliptic problems[END_REF] and [START_REF] Fonseca | An existence result for a nonconvex variational problem via regularity[END_REF] for example. For completeness, we also point out the papers [START_REF] Colombo | An excess-decay result for a class of degenerate elliptic equations[END_REF][START_REF] Colombo | Regularity results for very degenerate elliptic equations[END_REF] and [START_REF] Santambrogio | Continuity in two dimensions for a very degenerate elliptic equation[END_REF] for related regularity results on the term ∇G(∇u), when m > 0.

Another type of well-studied functionals having some similarities with F is given by (see for example [START_REF] Bildhauer | A regularity theory for scalar local minimizers of splitting-type variational integrals[END_REF][START_REF] Bildhauer | Variational integrals with a wide range of anisotropy[END_REF] and [START_REF] Fusco | Some remarks on the regularity of minima of anisotropic integrals[END_REF]Section 4])

(1.8) ˆ G(∇u) dx, with G(z) = N i=1 (µ + |z i | 2 ) p i 2 .
Here µ > 0 and 1 < p 1 ≤ p 2 ≤ • • • ≤ p N are possibly different exponents. When the p i are not equal, such a functional belongs to the class of problems with non standard growth conditions, whose systematic study started with the paper [START_REF] Marcellini | Regularity of minimizers of integrals of the Calculus of Variations under non standard growth conditions[END_REF] by Marcellini. In this case we can infer local Lipschitz continuity if the exponents p i are not "too far apart" (see the above mentioned references for more details). However, our functional F does not fall neither in the class of the functional (1.5) nor in that of (1.8). Indeed, observe that in our case

F (z) = N i=1 (|z i | -δ i ) p + p ,
verifies (1.7), but (1.6) crucially fails to hold, since for every m > 0, there always exists z such that |z| = m and the least eigenvalue of D 2 F (z) is 0. Observe that this phenomenon already occurs for the pseudo p-Laplacian, i.e. when δ 1 = • • • = δ N = 0. Indeed, the main difficulty of the problem is that the region where ellipticity fails is unbounded.

For the same reason, F is not of the type (1.8), since already in the standard growth case

2 ≤ p 1 = p 2 = • • • = p N we have 0 < min |ξ|=1 D 2 G(z) ξ, ξ , z ∈ R N .
When one allows µ = 0 in (1.8), the corresponding functional becomes degenerate along the axes z i = 0, like in the case of the pseudo p-Laplacian. This case has been considered in the pioneering paper [START_REF] Uralt'seva | The boundedness of the gradients of generalized solutions of degenerate quasilinear nonuniformly elliptic equations[END_REF] by Uralt'seva and Urdaletova. There the Lipschitz character of minimizers has been shown under some restrictions on the exponents p 1 , . . . , p N , by using the so-called Bernstein method.

Though the growth conditions considered are more general than ours, the type of degeneracy is again weaker than that admitted in F (see the next subsection for more comments on the result of [START_REF] Uralt'seva | The boundedness of the gradients of generalized solutions of degenerate quasilinear nonuniformly elliptic equations[END_REF]).

About the restriction p ≥ 2 considered in this paper, it is noteworthy to observe that for 1 < p < 2 our functional has a p-Laplacian-type structure. Indeed, in this case p -2 < 0 and thus (1.6) is satisfied with m = 0, i.e.

D 2 F (z) ξ, ξ = (p -1) N i=1 (|z i | -δ i ) p-2 + |ξ i | 2 ≥ (p -1) |z| p-2 |ξ| 2 , ξ ∈ R N , z ∈ R N , while of course |DF (z)| ≤ |z| p-1 , z ∈ R N .
Then in this case local minimizers are locally Lipschitz continuous by 1 [16, Theorem 2.7].

1.2. Main results. In this paper, we prove the following results.

Theorem A (Two dimensional case). Let N = 2 and p ≥ 2. Let f ∈ W 1,p ′ loc (Ω), where p ′ = p/(p-1). Then every local minimizer U ∈ W 1,p loc (Ω) of the functional F is a locally Lipschitz continuous function.

Theorem B (Higher dimensional case). Let N ≥ 2 and p ≥ 4. Let f ∈ W 1,∞ loc (Ω). Then every local minimizer U ∈ W 1,p loc (Ω) of the functional F is a locally Lipschitz continuous function. Let us now spend some words about the methods of proofs. The preliminary step in both cases is an approximantion argument. Namely, the functional F is replaced by a regularized version F ε , for a small parameter ε > 0. This permits to infer the necessary regularity on the solutions u ε of the regularized problem, in order to justify the manipulations needed to obtain a priori Lipschitz estimates uniform in ε. Then one aims at taking these estimates to the limit as ε goes to 0. However, one should pay attention to the fact that F is not strictly convex when at least one δ i = 0. Thus a sequence of solutions u ε may not necessarily converge to the selected local minimizer. In [START_REF] Brasco | On certain anisotropic elliptic equations arising in congested optimal transport: local gradient bounds[END_REF] a penalization argument was used to fix this issue. Here on the contrary, we use a simpler argument, based on the fact that the lack of strict convexity of t → (|t| -δ i ) p + is "confined" (see Lemma 2.3). The core of the proof of Theorem A is the a priori Lipschitz estimate of Proposition 4.1. Such an estimate is achieved by means of a Moser's iteration technique applied to the equation solved by the partial derivatives u x j of the local minimizer. More precisely, we look at power-type subsolutions of this equation, i.e. quantities like |u x j | s for s ≥ 1. This is a standard strategy for equations having a p-Laplacian-type structure, but as already said our operator does not have such a structure and this entails several additional difficulties.

As explained in the introduction of [START_REF] Brasco | On certain anisotropic elliptic equations arising in congested optimal transport: local gradient bounds[END_REF], the main difficulty of this method is that the Caccioppoli inequality we get for |u x j | s is quite involved. Indeed, due to the particular structure of D 2 F , in principle we have a control only on a "weighted" norm of ∇|u x j | s , the weights being dependent on all the other components u x i of the gradient (see Lemma 3.6 below). Roughly speaking, what we control in the Caccioppoli inequality is a quantity like

N i=1 ˆ|u x i | p-2 |u x j | s+1 x i 2 .
For the diagonal term, i.e. when i = j, we can combine the x j -derivative of u x j with the weigth |u x j | p-2 and simply recognize the x j -derivative of yet another power of u x j . Since we would like to have a control on the full gradient of such a power of u x j , we still miss all the x i -derivatives (i = j) of this function. To overcome this difficulty, we use in a crucial way the Sobolev property (1.2) together with Hölder inequality, in order to "cook-up" suitable Caccioppoli inequalities for all these missing terms. Surprisingly enough, even if the functional F has p-growth in every direction, we rely on the anisotropic Sobolev inequality due to Troisi (see [START_REF] Troisi | Teoremi di inclusione per spazi di Sobolev non isotropi[END_REF]) in order to produce an iterative scheme of reverse Hölder inequalities. This procedure works for N = 2, but it seems to be limited just to the two dimensional case (see Remark 4.2 below). 1 To be more precise, for 1 < p < 2 the function F is not C 2 . However, this is not an issue, since the result of [START_REF] Fonseca | An existence result for a nonconvex variational problem via regularity[END_REF]Theorem 2.7] holds for convex functions satisfying a qualified form of uniform convexity for |z| ≥ m. This coincides with (1.6) if the function is C 2 , but it is otherwise more general.

In contrast Theorem B is valid in every dimension, but we need the restriction p ≥ 4. This second result partially superposes with the already mentioned [25, Theorem 1] by Uralt'seva and Urdaletova. However, it should be noticed that the monotonicity assumptions on the operator 2 made in [START_REF] Uralt'seva | The boundedness of the gradients of generalized solutions of degenerate quasilinear nonuniformly elliptic equations[END_REF] does not allow for δ i > 0. Moreover, the result in [START_REF] Uralt'seva | The boundedness of the gradients of generalized solutions of degenerate quasilinear nonuniformly elliptic equations[END_REF] is stated for p > 3, but a careful inspection of the proof reveals that the same condition p ≥ 4 is needed there as well 3 .

Both the proofs of Theorem B and that of [START_REF] Uralt'seva | The boundedness of the gradients of generalized solutions of degenerate quasilinear nonuniformly elliptic equations[END_REF]Theorem 1] are based on a priori Lipschitz bounds, obtained by means of pointwise estimates in the vein of Bernstein method. However, computations are not the same and we believe ours to be slightly simpler. In [START_REF] Uralt'seva | The boundedness of the gradients of generalized solutions of degenerate quasilinear nonuniformly elliptic equations[END_REF] the first step is to look at the equation solved by a concave power of u, given by the function

w = (u + u L ∞ + 1) γ , 0 < γ < 1.
Then they consider the equation solved by (some function of) ∇w. There is an extra term in this new equation coming from the concave power which crucially leads to the result.

Here on the contrary we obtain the Lipschitz estimate by directly attacking equation (1.3). The main point is to consider the equation satisfied by the quantity

|∇u| 2 + λ u 2 ,
for a suitably large paramater λ. We notice that this is exactly the same test function used to prove classical gradient estimates for linear uniformly elliptic equations (see for example [START_REF] Han | Elliptic partial differential equations[END_REF]Proposition 2.19]).

One of the drawbacks of these two strategies is the assumption on f , which does not seem to be optimal. Indeed, we expect the result to be true under the natural hypothesis f ∈ L q loc (Ω) with q > N .

1.3. Plan of the paper. In Section 2 we set notations and preliminary results needed throughout the whole paper. In particular, we introduce there a regularized version of the problem which will be useful in order to get the desired Lipschitz estimate. Then Section 3 is devoted to prove some Caccioppoli-type inequalities for the gradient of the solution of the regularized problem. The proof of Theorem A is contained in 4, while Section 5 contains the proof of Theorem B. Two appendices containing some technical results complement the paper.

Acknowledgements. A quick but stimulating discussion with Nina Uralt'seva in June 2012 led to a better understanding of the paper [START_REF] Uralt'seva | The boundedness of the gradients of generalized solutions of degenerate quasilinear nonuniformly elliptic equations[END_REF], we thank her. Guillaume Carlier is warmly thanked for his interest in this work. Part of this paper has been written during the conferences "Journées d'Analyse Appliquée Nice-Toulon-Marseille" held in Porquerolles in May 2014, "Nonlinear partial differential equations and stochastic methods" held in Jyväskylä in June 2014 and "Existence and Regularity for Nonlinear Systems of Partial Differential Equations" held in Pisa in July 2014. Organizers and hosting institutions are gratefully acknowledged.

Preliminaries

g i (t) = 1 p (|t| -δ i ) p + , t ∈ R, i = 1, . . . , N,
2 See equation (8) of the paper [START_REF] Uralt'seva | The boundedness of the gradients of generalized solutions of degenerate quasilinear nonuniformly elliptic equations[END_REF]. 3 This comes from hypothesis (5) in [START_REF] Uralt'seva | The boundedness of the gradients of generalized solutions of degenerate quasilinear nonuniformly elliptic equations[END_REF]. Also observe that this condition contains a small typo, mi-2 should be replaced by mi -2.

where 0 ≤ δ 1 , . . . , δ N are given real numbers. We will also define (2.1) δ = 1 + max{δ i : i = 1, . . . , N }.

Remark 2.1 (Smoothness of g i ). When p is an integer and

δ i > 0, g i is of class C p-1,1 . When p ∈ N, then g i ∈ C [p],p-[p] (R) where [ • ] denotes the integer part.
Remark 2.2 (The limit case p = 2). Observe that for p = 2 and δ i > 0, we have

g i ∈ C 1,1 (R) ∩ C ∞ (R \ {1, -1}), but g i ∈ C 2 (R).
In this case, like in [START_REF] Brasco | On certain anisotropic elliptic equations arising in congested optimal transport: local gradient bounds[END_REF] a smoothing around |t| = δ i would be necessary, notably for the result of Lemma 2.8 below. However, in order not to overburden the presentation, for the sequel we will assume for simplicity p > 2 (see [9, Section 2] for more details).

We are interested in local minimizers of the following variational integral

(2.2) F(u; Ω ′ ) = N i=1 ˆΩ′ g i (u x i ) dx + ˆΩ′ f u dx, u ∈ W 1,p loc (Ω),
where f ∈ L p ′ loc (Ω) and Ω ′ ⋐ Ω. We recall that u ∈ W 1,p loc (Ω) is said to be a local minimizer of F if for every Ω ′ ⋐ Ω we have

F(u; Ω ′ ) ≤ F(u + ϕ; Ω ′ ),
for every ϕ ∈ W 1,p 0 (Ω ′ ). We first observe that F is not strictly convex, unless δ = 1. Thus minimizers are not unique in general. The following result guarantees that it will be sufficient to prove the desired result for one minimizer.

Lemma 2.3 (Propagation of regularity). Let B ⋐ Ω be a ball and ϕ ∈ W 1,p (B). Let u 1 , u 2 ∈ W 1,p (Ω) be two solutions of (2.3) min F(v; B) : v -ϕ ∈ W 1,p 0 (B) . Then it holds (2.4) |(u 1 ) x i | -|(u 2 ) x i | ≤ 2 δ i , a. e. in B, i = 1, . . . , N.
In particular, if a minimizer of (2.3) is (locally) Lipschitz, then this remains true for all the other minimizers.

Proof. Let us suppose that (2.4) is not true. Then there exists i 0 ∈ {1, . . . , N } such that

E i 0 := x ∈ B : |(u 1 ) x i 0 | -|(u 2 ) x i 0 | > 2 δ i ,
has strictly positive measure. We then set u s = (1 -s) u 1 + s u 2 for some s ∈ (0, 1) and observe that this is admissible in (2.3). In view of Lemma A.1 in Appendix A,

g i 0 (1 -s) (u 1 ) x i 0 + s (u 2 ) x i 0 < (1 -s) g i 0 ((u 1 ) x i 0 ) + s g i 0 ((u 2 ) x i 0 ), a. e. inE i 0 .
Thus we get

F(u s ) < (1 -s) F(u 1 ) + s F(u 2 ) = F(u 1 ) = F(u 2
), which gives the desired contradiction.

We will also need the following regularity result, which is reminiscent of [START_REF] Stampacchia | On some regular multiple integral problems in the calculus of variations[END_REF]. A more general result of this type can be found in [START_REF] Bousquet | Global Lipschitz continuity for minima of degenerate problems[END_REF].

Theorem 2.4 ([6]). Let B ⊂ R N be a ball, ϕ ∈ C 2 (B) and f ∈ L ∞ (B). Let u ∈ W 1,1 (B) ∩ L ∞ (B) be a solution of min ˆB H(∇v) dx + ˆB f v dx : v -ϕ ∈ W 1,1 0 (B) ,
where

H : R N → [0, ∞) is a C 2 convex function such that for some µ > 0 (2.5) D 2 H(z) ξ, ξ ≥ µ |ξ| 2 , ξ, z ∈ R N .
Then u ∈ W 1,∞ loc (B). 2.2. Approximation scheme. We now introduce a regularized version of the original problem. We set

(2.6) g i.ε (t) = g i (t) + ε 2 t 2 = 1 p (|t| -δ i ) p + + ε 2 |t| 2 , t ∈ R.
Let U be a local minimizer of F. We also fix a ball

B ⋐ Ω such that 2 B ⋐ Ω as well.
Here 2 B denotes the ball having the same center as B scaled by a factor 2. For every 0 < ε ≪ 1 and every x ∈ B, we set

U ε (x) = U * ̺ ε (x)
, where ̺ ε is a smooth convolution kernel, supported in a ball of radius ε centered at the origin.

Then by definition of U ε there exists 0 < ε 0 < 1 such that for every 0

< ε ≤ ε 0 (2.7) U ε W 1,p (B) = ∇U ε L p (B) + U ε L p (B) ≤ ∇U L p (2 B) + U L p (2 B) =: C 1 .
Finally, we define

F ε (v; B) = N i=1 ˆB g i,ε (v x i ) dx + ˆB f ε v dx,
where

f ε = f * ̺ ε .
The following preliminary result is standard.

Lemma 2.5 (Basic energy estimate).

There exists a unique solution u ε to the problem

(2.8) min F ε (v; B) : v -U ε ∈ W 1,p 0 (B) .
The following uniform energy estimate holds

(2.9) ˆB |∇u ε | p dx ≤ C 2 ,
for some constant

C 2 = C 2 (N, δ, p, |B|, C 1 , f L p ′ (2 B) ) > 0.
Proof. We start by observing that a solution u ε exists, by a standard application of the Direct Methods. Uniqueness then follows from strict convexity of the integrand

(2.10) L ε (x, u, z) = N i=1 g i (z i ) + ε 2 |z| 2 + f ε (x) u,
in the gradient variable. In order to prove (2.9), we use the minimality of u ε , which implies

F ε (u ε ; B) ≤ F ε (U ε ; B). This gives N i=1 ˆB g i,ε ((u ε ) x i ) dx ≤ N i=1 ˆB g i,ε ((U ε ) x i ) dx + ˆB |f ε | |u ε -U ε | dx.
We now use the fact that

(2.11) 1 p 1 2 p-1 |t| p -δ p ≤ g i,ε (t) ≤ 2 p |t| p + p -2 2p .
The lower bound in (2.11) follows from

|t| p ≤ 2 p-1 ((|t| -δ i ) p + + δ p i
), while the upper bound is a consequence of Young inequality. This implies

N i=1 ˆB |(u ε ) x i | p dx ≤ C N i=1 ˆB |(U ε ) x i | p + ˆB |f ε ||u ε -U ε | dx + C
where C = C(N, p, δ, |B|) > 0 depends on N, p, δ and |B| only. By using

f ε L p ′ (B) ≤ f L p ′ (2 B)
and (2.7), standard computations lead to the desired conclusion.

Lemma 2.6 (Regularity of the minimizer I). Let u ε still denote the unique minimizer of (2.8).

Then we have u ε ∈ L ∞ (B). Moreover, if f ∈ L ∞ loc (Ω), then there exists a constant M independent of ε such that (2.12) u ε L ∞ (B) ≤ M.
Proof. We use again (2.11). This implies that the integrand (2.10) satisfies

(2.13) c |z| p -f ε L ∞ (B) |u| -C ′ ≤ L ε (x, u, z) ≤ 1 c |z| p + f ε L ∞ (B) |u| + C ′ , with c = c(N, p) > 0, C = C(N, p) > 0 and C ′ = C ′ (p, δ) > 0. Thus u ε ∈ L ∞ (B) by [18, Theorem 7.5 & Remark 7.6]. If f ∈ L ∞ loc (Ω), we have f ε L ∞ (B) ≤ f L ∞ (2 B)
. By (2.13) we thus get that L ε satisfies growth conditions independent of ε. Then by using again the a priori estimate of [18, Theorem 7.5] and (2.9) we get the desired conclusion.

Remark 2.7. The previous L ∞ estimate uniform in ε will be needed in the proof of Theorem B.

The following result is not optimal, but it is suitable to our needs. Lemma 2.8 (Regularity of the minimizer II). Let u ε still denote the unique minimizer of (2.8). We have

u ε ∈ C k loc (B), where k = 2, if 2 < p ≤ 3, 3, if p > 3.
Proof. We divide the proof in two parts.

Local Lipschitz regularity. By Lemma 2.6, we know that u ε is bounded. Then the local Lipschitz continuity is a plain consequence of Theorem 2.4, applied with

(2.14) F ε (z) = N i=1 g i (z i ) + ε 2 |z| 2 , z ∈ R N ,
which verifies (2.5) with µ = ε > 0.

Local higher regularity. Let B ⋐ B be a ball and set ℓ = ∇u ε L ∞ ( B) , which is finite thanks to the previous step. By optimality, we have that u ε solves the elliptic equation

(2.15) div(∇F ε (∇u ε )) = f ε , in B,
where F ε is as in (2.14). Since we have

ε |ξ| 2 ≤ D 2 F ε (∇u ε ) ξ, ξ ≤ ε + (p -1) ℓ p-2 |ξ| 2 , on B,
we can infer u ε ∈ W 2,2 loc ( B) by a standard differential quotients argument (see for example [START_REF] Giusti | Metodi diretti nel calcolo delle variazioni. (Italian) [Direct methods in the calculus of variations[END_REF]Theorem 8.1]). This in turn permits to find the equation locally solved by ∇u ε , by differentiating (2.15). Thus ∇u ε ∈ C 0,σ loc ( B) by the celebrated De Giorgi-Moser-Nash Theorem, for some σ > 0. It remains to observe that F ε ∈ C k,α , where k is as in the statement and

α = min{p -2, 1}, if 2 < p ≤ 3, min{p -3, 1}, if p > 3.
Then [START_REF] Giusti | Metodi diretti nel calcolo delle variazioni. (Italian) [Direct methods in the calculus of variations[END_REF]Theorem 10.18] implies that u ε has the claimed regularity properties.

Lemma 2.9 (Convergence to a minimizer). With the same notation as before, we have

lim ε→0 u ε -u L p (B) = 0,
where u is a solution of

(2.16) min F(ϕ; B) : ϕ -U ∈ W 1,p 0 (B) .
Proof. By (2.9), there exists a sequence {ε k } k∈N converging to 0 as k goes to ∞ and a function u ∈ W 1,p (B) such that {u ε k } k∈N converges weakly to u in W 1,p (B) and strongly in L p (B). The function U ε = U * ̺ ε is of course admissible for the approximated problem and thus

lim inf k→∞ F ε k (U ε k ; B) ≥ lim inf k→∞ F ε k (u ε k ; B) ≥ lim inf k→∞ F(u ε k ; B) ≥ F( u; B),
where we used the weak lower semicontinuity of F. We then observe that by using the strong convergence of U ε to U and inequality (A.2) in Appendix A, we get

lim k→∞ |F ε k (U ε k ; B) -F(U ; B)| ≤ lim k→∞ N i=1 ˆB |g i ((U ε k ) x i ) -g i (U x i )| dx + lim k→∞ ε k 2 ˆB |∇U ε k | 2 + lim k→∞ ˆB |f ε k U ε k -f U | dx = 0,
and thus

F(U ; B) = lim k→∞ F ε k (U ε k ; B) ≥ F( u; B).
By definition of local minimizer, the function U itself is a solution of (2.16), then the previous inequality implies that u is a minimizer.

Local energy estimates for the approximating problem

For the ball B ⋐ Ω we consider the regularized problem (2.8). We still denote by u ε its unique solution, which verifies the Euler-Lagrange equation

(3.1) N i=1 ˆg′ i,ε ((u ε ) x i ) ϕ x i dx + ˆfε ϕ dx = 0, ϕ ∈ W 1,p 0 (B).
From now on, in order to simplify the notation, we will systematically forget the subscript ε on u ε and simply write u.

We now insert a test function of the form ϕ = ψ x j ∈ W 1,p 0 (B) in (3.1), compactly supported in B. Then an integration by parts lead us to

(3.2) N i=1 ˆg′′ i,ε (u x i ) u x i x j ψ x i dx -ˆfε ψ x j dx = 0,
for j = 1, . . . , N . This is the equation solved by u x j .

Caccioppoli-type inequalities.

In what follows we use the parameter δ defined in (2.1). The general Caccioppoli inequality for an important class of subsolutions is given by the following result.

Lemma 3.1. Let Φ : R → R + be a C 2 convex function such that

(3.3) Φ ′ (t) ≡ 0 for |t| ≤ δ.
Then there exists a constant C 3 = C 3 (p) > 0 such that for every Lipschitz function η with compact support in B, we have

N i=1 ˆAj g ′′ i,ε (u x i ) Φ(u x j ) x i 2 η 2 dx ≤ C 3 N i=1 ˆAj g ′′ i,ε (u x i ) |Φ(u x j )| 2 |η x i | 2 dx + C 3 ˆAj |f ε | 2 Φ ′ (u x j ) 2 + Φ ′′ (u x j ) Φ(u x j ) η 2 dx + C 3 ˆAj Φ(u x j ) 2 |η x j | 2 dx, (3.4)
where we set A j = {x ∈ B : |u x j | ≥ δ}.

Proof. In (3.2) we take the test function 4 ψ = ζ Φ ′ (u x j ), with Φ : R → R + as in the statement and ζ nonnegative Lipschitz function with support in B. We thus obtain N i=1 ˆAj

g ′′ i,ε (u x i ) Φ(u x j ) x i ζ x i dx + N i=1 ˆAj g ′′ i,ε (u x i ) u 2 x i x j Φ ′′ (u x j ) ζ dx = ˆAj f ε ζ Φ ′ (u x j ) x j dx.
Finally, we test the previous equation against ζ = η 2 Φ(u x j ), where η is again a Lipschitz function with support in B. Then we get N i=1 ˆAj

g ′′ i,ε (u x i ) Φ(u x j ) x i 2 η 2 dx + S(η) ≤ 2 N i=1 ˆAj g ′′ i,ε (u x i ) Φ(u x j ) x i Φ(u x j ) η |η x i |dx + ˆAj |f ε | η 2 Φ(u x j ) Φ ′ (u x j ) x j dx.
where we have introduced the sponge term

S(η) = N i=1 ˆAj g ′′ i,ε (u x i ) u 2 x i x j Φ ′′ (u x j ) Φ(u x j ) η 2 dx.
From the previous inequality, by Young inequality in the first term on the right-hand side

N i=1 ˆAj g ′′ i,ε (u x i ) Φ(u x j ) x i 2 η 2 dx + 2 S(η) ≤ 4 N i=1 ˆAj g ′′ i,ε (u x i ) Φ(u x j ) 2 |η x i | 2 dx + 2 ˆAj |f ε | η 2 Φ(u x j ) Φ ′ (u x j ) x j dx. (3.5)
We now estimate the term containing f ε . We first observe

ˆAj |f ε | η 2 Φ(u x j ) Φ ′ (u x j ) x j dx ≤ ˆAj |f ε | |Φ ′ (u x j )| Φ(u x j ) x j η 2 dx + 2 ˆAj |f ε | |Φ ′ (u x j )| Φ(u x j ) η |η x j | dx + ˆAj |f ε | Φ ′ (u x j ) x j Φ(u x j ) η 2 dx.
On the set A j we have

(3.6) g ′′ j,ε (u x j ) ≥ (p -1)
. Let us consider the first term above containing f ε :

ˆAj |f ε | |Φ ′ (u x j )| Φ(u x j ) x j η 2 dx ≤ 1 2 τ ˆAj |f ε | 2 |Φ ′ (u x j )| 2 η 2 dx + τ 2 ˆAj Φ(u x j ) x j 2 η 2 dx ≤ 1 2 τ ˆAj |f ε | 2 |Φ ′ (u x j )| 2 η 2 dx + τ 2 (p -1) ˆAj g ′′ j,ε (u x j ) Φ(u x j ) x j 2 η 2 dx.
The last term can be absorbed in the left-hand side of (3.5), by taking τ = (p -1)/2. The second term containing f ε is simply estimated by Young inequality

ˆAj |f ε | |Φ ′ (u x j )| Φ(u x j ) η |η x j | dx ≤ 1 2 ˆAj |f ε | 2 |Φ ′ (u x j )| 2 η 2 dx + 1 2 ˆAj Φ(u x j ) 2 |η x j | 2 dx,
while for the last one we use the sponge term S(η) to absorb the Hessian of u. Namely, we have

ˆAj |f ε | Φ ′ (u x j ) x j Φ(u x j ) η 2 dx = ˆAj |f ε | |u x j x j | Φ ′′ (u x j ) Φ(u x j ) η 2 dx ≤ τ ˆAj u 2 x j x j Φ ′′ (u x j ) Φ(u x j ) η 2 dx + 1 τ ˆAj |f ε | 2 Φ ′′ (u x j ) Φ(u x j ) η 2 dx ≤ τ p -1 S(η) + 1 τ ˆAj |f ε | 2 Φ ′′ (u x j ) Φ(u x j ) η 2 dx.
In the last estimate we used again (3.6). The term τ /(p -1) S(η) can then be absorbed in the left-hand side. This concludes the proof.

If we allow for derivatives of f ε on the right-hand side of (3.4), the previous estimate is simpler to get. In this case we can allow for more general subsolutions. Lemma 3.2 (Right-hand side in a Sobolev space). Let Φ : R → R be a C 1 convex function. Then there exists a constant C 3 = C 3 (p) > 0 such that for every Lipschitz function η with compact support in B, we have

N i=1 ˆg′′ i,ε (u x i ) Φ(u x j ) x i 2 η 2 dx ≤ C 3 N i=1 ˆg′′ i,ε (u x i ) |Φ(u x j )| 2 |η x i | 2 dx + C 3 ˆ|(f ε ) x j | |Φ ′ (u x j )| |Φ(u x j )| η 2 dx. (3.7)
Proof. Let us suppose for simplicity that Φ ∈ C 2 . If this were not the case, a standard smoothing argument will be needed, we leave the details to the reader. We start observing that equation (3.2) can also be written as

N i=1 ˆg′′ i,ε (u x i ) u x i x j ψ x i dx + ˆ(f ε ) x j ψ dx = 0, j = 1, . . . , N. (3.8)
Then we take in (3.8) the test function ψ = ζ Φ ′ (u x j ) as before, with Φ as in the statement and ζ a nonnegative Lipschitz function supported in B. We obtain

N i=1 ˆg′′ i,ε (u x i ) Φ(u x j ) x i ζ x i dx ≤ -ˆ(f ε ) x j Φ ′ (u x j ) ζ dx.
Finally, we take again ζ = η 2 Φ(u x j ), to get

N i=1 ˆg′′ i,ε (u x i ) Φ(u x j ) x i 2 η 2 dx ≤ 2 N i=1 ˆg′′ i,ε (u x i ) Φ(u x j ) x i Φ(u x j ) η |η x i | dx + ˆ (f ε ) x j |Φ ′ (u x j )| |Φ(u x j )| η 2 dx.
By using Young inequality as before, we get

N i=1 ˆg′′ i,ε (u x i ) Φ(u x j ) x i 2 η 2 dx ≤ 4 N i=1 ˆg′′ i,ε (u x i ) Φ(u x j ) 2 |η x i | 2 dx + 2 ˆ|(f ε ) x j | |Φ ′ (u x j )| |Φ(u x j )| η 2 dx.
This concludes the proof.

A Sobolev estimate.

In what follows we set

W j = δ 2 + (|u x j | -δ) 2 + . Lemma 3.3.
There exists a constant C 4 = C 4 (p, δ) > 0 such that for every Lipschitz function η with compact support in B, we have

N i=1 ˆ ∇W p 4 i 2 η 2 dx ≤ C 4 N i,j=1 ˆW p-2 2 i W j |η x i | 2 dx + C 4 N j=1 ˆ|(f ε ) x j | W j η 2 dx. (3.9)
Proof. We choose the function Φ(t) = t in (3.7). First observe that Φ(u x j ) x i 2 = u 2

x i x j , and that by Lemma A.3 for |t| ≥ δ

g ′′ i,ε (t) ≥ (p -1) δ -δ i δ p-2 δ 2 + (|t| -δ) 2 + p-2 2 ≥ 1 δ p-2 δ 2 + (|t| -δ) 2 + p-4 2 (|t| -δ) 2 + .
(3.10)

In the second inequality above we also used that p ≥ 2 and δ -δ i ≥ 1. Then, by (3.10) we have 5

g ′′ i,ε (u x i ) Φ(u x j ) x i 2 ≥ 1 δ p-2 W p-4 2 i (|u x i | -δ) 2 + u 2 x i x j = c ∂ x j W p 4 i 2 ,
where c = c(δ, p) > 0. We further observe that

|Φ(u x j )| = |u x j | ≤ √ 2 W j , and 6
(3.11)

g ′′ i,ε (u x i ) = (p -1) (|u x i | -δ i ) p-2 + + ε ≤ c W p-2 2 i
, where c = c(p) > 0. Then we get the desired result by summing (3.7) over j = 1, . . . , N .

In what follows, we will use for simplicity the notation

E ϕ dx := 1 |E| ˆE ϕ dx.
Corollary 3.4. There exists a constant C ′ 4 = C ′ 4 (p, δ, N ) > 0 such that for every pair of concentric balls B R 0 ⋐ B ρ 0 ⋐ B, we have

N j=1 1 R N -2 0 ˆBR 0 ∇W p 4 j 2 dx ≤ C ρ 0 R 0 N -2 ρ 0 ρ 0 -R 0 2 N j=1 Bρ 0 W p 2 j dx + C ρ 2 p-1 -N +2 0 N j=1 ˆBρ 0 |(f ε ) x j | p ′ dx.
(3.12)

Proof. It is sufficient to insert the test function

η(x) = min 1, (ρ 0 -|x|) + ρ 0 -R 0 ,
in (3.9) and then use Hölder and Young inequalities in the right-hand side. These give

N i,j=1 ˆW p-2 2 i W j |η x i | 2 dx ≤ 1 (ρ 0 -R 0 ) 2 N i,j=1 ˆBρ 0 W p 2 i dx p-2 p ˆBρ 0 W p 2 j dx 2 p
, 5 Observe that the inequality holds true everywhere, not only on Ai, since Wi is constant outside Ai. 6 We use that tδi = (tδ) + (δδi) ≤ (tδ)+ + δ, which implies

(t -δi)+ ≤ (t -δ)+ + δ. and N j=1 ˆ|(f ε ) x j | W j η 2 dx ≤ C N j=1 ρ 2 p-1 0 ˆBρ 0 |(f ε ) x j | p ′ dx + C ρ 2 0 N j=1 ˆBρ 0 W p 2 j dx,
which concludes the proof.

Remark 3.5 (Uniform Sobolev estimate). From the previous result, we obtain that if f ∈ W 1,p ′ loc (Ω), then for every i = 1, . . . , N the function W p/4 i enjoys a W 1,2 loc (B) estimate independent of ε, thanks to (2.9) and

f ε W 1,p ′ (Bρ 0 ) ≤ f W 1,p ′ (2 B) .
3.3. Power-type subsolutions. We still use the notation

W j = δ 2 + (|u x j | -δ) 2 + .
Then we have the following result. Lemma 3.6. There exists a constant C 5 = C 5 (p) > 0 such that for every s ≥ 0 and every Lipschitz function η with compact support in B, we have

N i=1 ˆg′′ i,ε (u x i ) W s+1 2 j x i 2 η 2 dx ≤ C 5 N i=1 ˆW p-2 2 i W s+1 j |∇η| 2 dx + C 5 (s + 1) 2 ˆ|f ε | 2 W s j η 2 dx, j = 1, . . . , N. (3.13)
Proof. In equation (3.4) we make the choice 7

Φ(t) = δ 2 + (|t| -δ) 2 + s+1 2 ,
for s ≥ 0 which satisfies hypothesis (3.3). Observe that by definition we have

Φ(u x j ) = W s+1 2 j , so that Φ(u x j ) x i 2 = W s+1 2 j x i 2 .
Thus the left-hand side of (3.4) coincides with

N i=1 ˆg′′ i,ε (u x i ) W s+1 2 j x i 2 η 2 dx.
We now come to the right-hand side: 7 Observe that this function is not C 2 , but only C 1,1 near t = δ ot t = -δ. This is not a big issue, since in any case Φ ′′ stays bounded as |t| → δ, thus we can use (3.4) for a regularization of Φ and then pass to the limit at the end.

N i=1 ˆg′′ i,ε (u x i ) |Φ(u x j )| 2 |η x i | 2 dx = N i=1 ˆg′′ i,ε (u x i ) W s+1 j |η x i | 2 dx ≤ C N i=1 ˆW p-2 2 i W s+1 j |η x i | 2 dx,
thanks to (3.11). For the other two terms, by using the definition of Φ we simply have

ˆAj |f ε | 2 Φ ′ (u x j ) 2 + Φ ′′ (u x j ) Φ(u x j ) η 2 dx + ˆAj Φ(u x j ) 2 |η x j | 2 dx ≤ C (s + 1) 2 ˆ|f ε | 2 W s j η 2 dx + N i=1 ˆW p-2 2 i W s+1 j |η x j | 2 dx,
where we used that

Φ ′ (t) 2 + Φ ′′ (t) Φ(t) ≤ C (s + 1) 2 δ 2 + (|t| -δ) 2 + s = C (1 + s) 2 Φ(t) 2 s s+1 ,
and

W s+1 j ≤ N i=1 W p-2 2 i W s+1 j , which follows from W i ≥ 1.
In particular, we get an estimate for the diagonal terms, corresponding to i = j.

Corollary 3.7. There exists a constant C 6 = C 6 (p, δ) > 0 such that for every s ≥ 0 and every Lipschitz function η with compact support in Ω, we have

ˆ W p 4 + s 2 j x j 2 η 2 dx ≤ C 6 N i=1 ˆW p-2 2 i W s+1 j |∇η| 2 dx + C 6 (s + 1) 2 ˆ|f ε | 2 W s j η 2 dx, j = 1, . . . , N. (3.14)
Proof. We fix j, by keeping only the term i = j and dropping all the others in the left-hand side of (3.13), we get

ˆAj g ′′ j,ε (u x j ) W s+1 2 j x j 2 η 2 dx ≤ C 5 N i=1 ˆW p-2 2 i W s+1 j |∇η| 2 dx + C 5 (s + 1) 2 ˆ|f ε | 2 W s j η 2 dx.
We now observe that again by Lemma A.3 on A j we have

g ′′ j,ε (u x j ) ≥ p -1 δ p-2 δ 2 + (|u x j | -δ) 2 + p-2 2 = p -1 δ p-2 W p-2 2 j
, and that

W p-2 2 j W s+1 2 j x j 2 = 2 + 2 s p + 2 s 2 W p 4 + s 2 j x j 2 ≥ 2 p 2 W p 4 + s 2 j x j 2
, so that the conclusion follows.

Proof of Theorem A

The core of the proof of Theorem A is the a priori estimate of Proposition 4.1 below. We postpone its proof and proceed with that of Theorem A.

Proof. Let Ω ′ ⋐ Ω and set d = dist(Ω ′ , ∂Ω). We take r 0 ≤ d/10, then Ω ′ can be covered by a finite number of balls centered at points in Ω ′ and having radius r 0 . Let B r 0 := B r 0 (x 0 ) ⋐ Ω be one of these balls, it is clearly sufficient to show that

∇U L ∞ (Br 0 ) < +∞.
To this aim we take the solution u ε of the regularized problem (2.8) in the ball B := B 4 r 0 (x 0 ). Observe that by construction we have 2 B = B 8 r 0 (x 0 ) ⋐ Ω. Then there exists ε 0 = ε 0 (d) > 0 such that for every 0

< ε ≤ ε 0 f ε W 1,p ′ (B) ≤ f W 1,p ′ (2 B) .
By using estimate (4.2) below with R 0 = 2 r 0 and ρ 0 = 3 r 0 we get

(4.1) ∇u ε L ∞ (Br 0 ) ≤ C, for every 0 < ε ≤ ε 0 ,
where C > 0 depends only on p, δ, r 0 , f W 1,p ′ (2 B) and the constant C 2 in (2.9). We then observe that by Lemma 2.9, we can find a sequence {ε k } k∈N converging to 0 and such that {u ε k } converges strongly in L p (B) and weakly in W

(4.2) (u ε ) x i L ∞ (Br 0 ) ≤ C 7 R 0 R 0 -r 0 4 J (u ε , f ε ; R 0 , ρ 0 ) 2   B R 0 |(u ε ) x i | p dx 1 p + δ   ,
where C 7 = C 7 (p, δ) > 0 is a constant that only depends on p and δ and

J (u ε , f ε ; R 0 , ρ 0 ) = ρ 0 ρ 0 -R 0 2 Bρ 0 |∇u ε | p dx + δ p + ρ 2 p-1 0 ˆBρ 0 |∇f ε | p ′ dx + ρ 2 p 0 ˆBρ 0 |f ε | 2 p ′ dx 1 p ′ . (4.3)
Proof. For notational simplicity, we write again u in place of u ε . We still use the notation

W j = δ 2 + |u x j | -δ 2 + , j = 1, 2.
We give the proof for u x 1 , the one for u x 2 being exactly the same. By (3.14) we already know that

ˆ W p 4 + s 2 1 x 1 2 η 2 dx ≤ C 6 2 i=1 ˆW p-2 2 i W s+1 1 |∇η| 2 dx + C 6 (s + 1) 2 ˆ|f ε | 2 W s 1 η 2 dx,
where η is any Lipschitz function supported on B and such that 0 ≤ η ≤ 1. We add the term

ˆ|η x 1 | 2 W p 2 +s 1 dx,
on both sides of the previous inequality and observe that

ˆ W p 4 + s 2 1 x 1 2 η 2 dx + ˆW p 2 +s 1 |η x 1 | 2 dx ≥ 1 2 ˆ W p 4 + s 2 1 η x 1 2 dx.
We thus obtain

ˆ W p 4 + s 2 1 η x 1 2 dx ≤ C 2 i=1 ˆW p-2 2 i W s+1 1 |∇η| 2 dx + C (s + 1) 2 ˆ|f ε | 2 W s 1 η 2 dx, (4.4) with C = C(p, δ) > 0.
The main problem of the Caccioppoli inequality (3.13) is that apparently we can not use it to control the missing term

W p 4 + s 2 1 x 2 .
Thus there is an obstruction to derive estimates for ∇W

p 4 + s 2 1
which could lead to an interative scheme of reverse Hölder inequalities. In order to overcome this problem, we observe that

W p 4 + s 2 1 x 2 = p + 2 s p W p 4 1 x 2 W s 2
1 .

Then if we fix 1 < q < 2, by Hölder inequality with exponents 2/q and 2/(2 -q), we have

ˆ W p 4 + s 2 1 x 2 q η q dx 2 q ≤ p + 2 s p 2 ˆ W p 4 1 x 2 2 η 2 dx ˆspt(η) W q 2-q s 1 dx 2-q q .
The precise value of q will be specified later. We now add the term

ˆW p q 4 + s q 2 1 |η x 2 | q dx 2 q
, on both sides of the previous inequality and observe that by triangle inequality

ˆ W p 4 + s 2 1 x 2 q η q dx 2 q + ˆW p q 4 + s q 2 1 |η x 2 | q dx 2 q ≥ 1 2 ˆ W p 4 + s 2 1 η x 2 q dx 2 q
.

Thus we get

ˆ W p 4 + s 2 1 η x 2 q dx 2 q ≤ C (1 + s) 2 ˆ W p 4 1 x 2 2 η 2 dx ˆspt(η) W q 2-q s 1 dx 2-q q + C ˆW p q 4 + s q 2 1 |η x 2 | q dx 2 q , (4.5)
with C = C(p) > 0. For 0 < r < R < R 0 , we now take η ∈ W 1,∞ 0 (B R ) to be the standard cut-off function

η(x) = min 1, (R -|x|) + R -r ,
then by multiplying (4.4) and (4.5) we get

ˆ W p 4 + s 2 1 η x 1 2 dx ˆ W p 4 + s 2 1 η x 2 q dx 2 q ≤ C 1 (R -r) 2 2 i=1 ˆBR W p-2 2 i W s+1 1 dx + (s + 1) 2 ˆBR |f ε | 2 W s 1 dx × (s + 1) 2 ˆBR W p 4 1 x 2 2 dx ˆBR W q 2-q s 1 dx 2-q q + 1 (R -r) 2 ˆBR W p q 4 + s q 2 1 dx 2 q . (4.6)
We now estimate the terms appearing in (4.6). To this aim, it will be useful to introduce the quantity

I(W 1 , W 2 , f ε ; R 0 ) = 2 i=1 B R 0 W p 2 i dx + ˆBR 0 ∇W p 4 i 2 dx + R 2 p 0 ˆBR 0 |f ε | 2 p ′ dx 1 p ′ . (4.7)
Then we start with the first term on the right-hand side of (4.6). Observe that

2 i=1 ˆBR W p-2 2 i W s+1 1 dx = ˆBR W p 2 1 W s 1 dx + ˆBR W p-2 2 2 W 1 W s 1 dx.
We use Hölder inequality in conjunction with Sobolev-Poincaré inequality 8 , to get

ˆBR W p 2 1 W s 1 dx ≤ C B R 0 W p 2 1 dx + ˆBR 0 ∇W p 4 1 2 dx R 2 p ′ 0 ˆBR W s p 1 dx 1 p ≤ C I(W 1 , W 2 , f ε ; R 0 ) R 2 p ′ 0 ˆBR W s p 1 dx 1 p , (4.8)
8 Since we are in dimension N = 2, we have W 1,2 (BR 0 ) ֒→ L 2 p ′ (BR 0 ). Then we have

ˆBR W p 4 2 p ′ dx 1 p ′ dx ≤ C R 2 p ′ B R W p 4 2 dx + ˆBR ∇W p 4 2 dx , with a constant C = C(p) > 0.
and

ˆBR W p-2 2 2 W 1 W s 1 dx ≤ C   2 i=1 ˆBR 0 W p 4 i 2 p ′ dx 1 p ′   ˆBR W s p 1 dx 1 p ≤ C 2 i=1 B R 0 W p 2 i dx + ˆBR 0 ∇W p 4 i 2 dx R 2 p ′ 0 ˆBR W s p 1 dx 1 p ≤ C I(W 1 , W 2 , f ε ; R 0 ) R 2 p ′ 0 ˆBR W s p 1 dx 1 p , (4.9)
for some constant C = C(p) > 0 depening only on p.

The term containing f ε in (4.6) is estimated as follows. Observe that 9

W 1,p ′ (B R 0 ) ֒→ L 2 p ′ (B R 0 ), since 2 p ′ < (p ′ ) * ,
then by Hölder's inequality and the definition of

I(W 1 , W 2 , f ε ; R 0 ) ˆBR |f ε | 2 W s 1 dx ≤ ˆBR 0 |f ε | 2 p ′ dx 1 p ′ ˆBR W s p 1 dx 1 p ≤ I(W 1 , W 2 , f ε ; R 0 ) R -2 p 0 ˆBR W s p 1 dx 1 p (4.10)
For the last term on the right-hand side of (4.6), by Hölder inequality and estimate (4.8) we have

ˆBR W p q 4 + s q 2 1 dx 2 q ≤ C R 2 2 q -1 0 ˆBR W p 2 1 W s 1 dx ≤ C R 2 2 q -1 p 0 I(W 1 , W 2 , f ε ; R 0 ) ˆBR W s p 1 dx 1 p , (4.11)
where C = C(q) > 0.

Finally, for the left-hand side of (4.6), we have

ˆ W p 4 + s 2 1 η x 1 2 dx ˆ W p 4 + s 2 1 η x 2 q dx 2 q ≥ T 2 q ˆ W p 4 + s 2 1 η q * dx 4 q * . (4.12)
Here we used the anisotropic Sobolev-Troisi inequality (see Appendix B) for the compactly supported function W (p+2 s)/4 1 η. The exponent q * is defined by

q * = 2 q 2 -q , where 1 q = 1 2 1 2 + 1 q ,
so that q = 4 q 2 + q and q * = 4 q 2 -q , the constant T q only depends on q and it converges to 0 as q goes to 2. 9 The exponent 2 p ′ is well-defined even in the case p = 2.

By using (4.8), (4.9), (4.10), (4.11) and (4.12) in (4.6), we then arrive at

ˆBr W p 2 +s 1 2 q 2-q dx 2-q q ≤ C R 0 R -r 2 I(W 1 , W 2 , f ε ; R 0 ) R -2 p 0 ˆBR W s p 1 dx 1 p +(s + 1) 2 I(W 1 , W 2 , f ε ; R 0 ) R -2 p 0 ˆBR W s p 1 dx 1 p × (s + 1) 2 I(W 1 , W 2 , f ε ; R 0 ) ˆBR W q 2-q s 1 dx 2-q q + R 0 R -r 2 R 2 2 q -1 p -1 0 I(W 1 , W 2 , f ε ; R 0 ) ˆBR W s p 1 dx 1 p , (4.13)
for a constant C = C(p, q, δ) > 0. We now choose 1 < q < 2 as follows (4.14) q = 2 p p + 1 .

Observe that with such a choice, we have

q 2 -q = p and 2 q - 1 p -1 = 0.
We further observe that

≤ C I(W 1 , W 2 , f ε ; R 0 ) 2 R 0 R -r 2 + (s + 1) 2 2 R -2 p 0 ˆBR W s p 1 dx 2 p
.

By using that R 0 /(R -r) ≥ 1 and (s + 1) ≥ 1 and introducing the notation ϑ = p s, then the previous estimate finally gives

W 1 L 2 ϑ (Br) ≤ C I(W 1 , W 2 , f ε ; R 0 ) R 0 R -r 2 ϑ p + 1 2 p ϑ R -1 ϑ 0 W 1 L ϑ (B R ) , (4.15)
possibly for a different constant C = C(p, δ) > 0. This is the iterative scheme of reverse Hölder inequalities needed to launch a Moser's iteration.

We then fix the two radii R 0 > r 0 > 0 of the statement and consider the sequences

r k = r 0 + R 0 -r 0 2 k and ϑ k = 2 ϑ k-1 = 2 k ϑ 0 = 2 k-1 p.
Then iterating (4.15) infinitely many times with R = r k and r = r k+1 , we get

W 1 L ∞ (Br 0 ) ≤ C R 0 R 0 -r 0 8 I(W 1 , W 2 , f ε ; R 0 ) 4 B R 0 W p 2 1 dx 2 p
, for some constant C = C(p, δ) > 0. We notice that u 2

x 1 ≤ W 1 ≤ u 2 x 1 + δ 2 , by definition of W 1 . Then we obtain with simple manipulations

u x 1 L ∞ (Br 0 ) ≤ C R 0 R 0 -r 0 4 I(W 1 , W 2 , f ε ; R 0 ) 2   B R 0 |u x 1 | p dx 1 p + δ   ,
for a possibly different constant C = C(p, δ) > 0. Finally, by Corollary 3.4 the term I(W 1 , W 2 , f ε ; R 0 ) defined in (4.7) can be estimated as follows

I(W 1 , W 2 , f ε ; R 0 ) ≤ C ρ 0 ρ 0 -R 0 2   2 j=1 Bρ 0 |u x j | p dx + δ p   + C ρ 2 p-1 0 ˆBρ 0 |∇f ε | p ′ dx + ρ 2 p 0 ˆBρ 0 |f ε | 2 p ′ dx 1 p ′ .
This concludes the proof.

Remark 4.2. Observe that the previous strategy does not seem to work for N ≥ 3. Indeed, in this case we would have N -1 missing terms, i.e.

∂ x i W p 4 + s 2 1 , i = 2, . . . , N.
By proceeding as before for each of these terms, i.e. combining (3.9) and Hölder inequality, one would have on the left-hand side the term

ˆ W p 4 + s 2 1 η x 1 2 dx N i=2 ˆ W p 4 + s 2 1 η x i q dx 2 q
, which in turn can be estimated from below by Sobolev-Troisi inequality by

ˆBr W p 4 + s 2 1 q * dx 2 q * .
The right-hand side would still contain the term

ˆBR W q 2-q s 1 dx 2-q q .
The exponent q * is now defined by

q * = N q N -q , where 1 q = 1 N 1 2 + N -1 q , so that q = 2 N q 2 N + q -2 and q * = 2 N q 2 N -q -2 . Then Moser's iteration would work if q * 2 s > q 2 -q s ⇐⇒ q < 2 N -1 .
Of course, when N ≥ 3 the last condition does not fit with the requirement q > 1.

Proof of Theorem B

Proof. The proof is the same as that of Theorem A. The essential point is the uniform Lipschitz estimate of Proposition 5.1 below, which replaces that of Proposition 4.1.

Proposition 5.1 (Uniform Lipschitz estimate, p ≥ 4). Let N ≥ 3 and p ≥ 4. For every pair of concentric balls B r 0 ⋐ B R 0 ⋐ B, we have

(5.1) ∇u ε L ∞ (Br 0 ) ≤ C 8 ,
where

C 8 = C 8 (N, p, δ, R 0 -r 0 , M, f W 1,∞ (2 B) ) > 0 does not depend on ε.
Here the constant M is the same appearing in (2.12).

Proof. As usual, for notational simplicity we simply write u in place of u ε . By Lemma 2.8, we get that u is indeed a local C 3 solution of the equation (3.1) in B, i.e. it verifies

(5.2) N i=1 g ′ i,ε (u x i ) x i = f ε , in B ′ ,
for every B ′ ⋐ B. This means that pointwise we have

(5.3) N i=1 g ′′ i,ε (u x i ) u x i x i = f ε in B ′ .
We now derive the previous equation with respect to x j and obtain (5.4)

N i=1 g ′′′ i,ε (u x i ) u x i x i u x i x j + g ′′ i,ε (u x i ) u x i x i x j = (f ε ) x j , in B ′ .
We introduce the following linear differential operator (5.5)

L[ψ] = N i=1 g ′′′ i,ε (u x i ) u x i x i ψ x i + g ′′ i,ε (u x i ) ψ x i x i ,
then (5.4) can be simply written as L[u x j ] = (f ε ) x j . Also observe that

L[ϕ ψ] = ϕ L[ψ] + ψ L[ϕ] + 2 N i=1 g ′′ i,ε (u x i ) ϕ x i ψ x i .
Thus for u 2 x j we obtain

L[u 2 x j ] = 2 u x j L[u x j ] + 2 N i=1 g ′′ i,ε (u x i ) u x j x i 2 = 2 u x j (f ε ) x j + 2 N i=1 g ′′ i,ε (u x i ) u 2 x j x i .
By linearity of L we thus get

L |∇u| 2 = 2 N j=1 u x j (f ε ) x j + 2 N i,j=1 g ′′ i,ε (u x i ) u 2 x j x i .
We now fix a pair of concentric balls B r 0 ⋐ B R 0 ⋐ B as in the statement of Proposition 5.

1. Let ζ ∈ C 2 0 (B R 0 ) be a function such that 0 ≤ ζ ≤ 1 and (5.6) ζ = 1 on B r 0 , |∇ζ| 2 ≤ C (R 0 -r 0 ) 2 ζ and |D 2 ζ| ≤ C (R 0 -r 0 ) 2 ,
and consider in B R 0 the equation for the function ζ |∇u| 2 + λ u 2 . The crucial parameter λ will be chosen later. By using the product rule for L and its linearity, we get

L ζ |∇u| 2 + λ u 2 = ζ L |∇u| 2 + |∇u| 2 L[ζ] + 2 N i,j=1 g ′′ i,ε (u x i ) (u 2 x j ) x i ζ x i + λ L[u 2 ] = 2 N j=1 u x j (f ε ) x j ζ + 2 N i,j=1 g ′′ i,ε (u x i ) u 2 x j x i ζ + |∇u| 2 L[ζ] + 2 N i,j=1 g ′′ i,ε (u x i ) (u 2 x j ) x i ζ x i + 2 u λ L[u] + 2 λ N i=1 g ′′ i,ε (u x i ) u 2 x i .
By using the expression (5.5) of L and the equation ( 5.3), we can rewrite the previous identity as follows

L ζ |∇u| 2 + λ u 2 = 2 I + 2 ζ G 1 + 2 G 2 + 2 λ G 3 + G 4 , (5.7)
where we used the notation

I = λ u f ε + ζ N j=1 u x j (f ε ) x j , G 1 = N i,j=1 g ′′ i,ε (u x i ) u x j x i 2 G 2 = N i,j=1 g ′′ i,ε (u x i ) (|u x j | 2 ) x i ζ x i + |∇u| 2 2 N i=1 g ′′ i,ε (u x i ) ζ x i x i , G 3 = N i=1 g ′′ i,ε (u x i ) |u x i | 2 ,
and

G 4 = N i=1 g ′′′ i,ε (u x i ) u x i x i 2 λ u u x i + |∇u| 2 ζ x i .
We proceed to estimate separately each term on the right-hand side of (5.7).

The term I.

For this, by Young inequality we get

I ≥ -λ u L ∞ (B R 0 ) f ε L ∞ (B R 0 ) -N ζ |∇u| |∇f ε | ≥ -λ M f L ∞ (2 B) - N p ζ |∇u| p - N p ′ ζ ∇f p ′ L ∞ (2 B) , (5.8)
where M is the constant appearing in (2.12). We also used that

f ε W 1,∞ (B R 0 ) ≤ f W 1,∞ (2 B) .
The term G 1 .

This is a positive term and for the moment we simply keep it. It will act as a sponge term, in order to absorb (negative) terms containing the Hessian of u.

The term G 2 .

This can be estimated by Young inequality and (5.6) as follows

G 2 ≥ - τ 2 N i,j=1 g ′′ i,ε (u x i ) u 2 x j x i ζ 2 x i - 1 2 τ N i,j=1 g ′′ i,ε (u x i ) u 2 x j - |∇u| 2 2 N i=1 g ′′ i,ε (u x i ) |ζ x i x i | ≥ -τ C 2 (R 0 -r 0 ) 2 ζ G 1 - |∇u| 2 2 C (R 0 -r 0 ) 2 + 1 τ N i=1 g ′′ i,ε (u x i ),
where τ < 1 is a small positive parameter. We then observe that the last term can be further estimated by using (5.9)

g ′′ i,ε (u x i ) ≤ (p -1) |∇u| p-2 + 1, so that |∇u| 2 N j=1 g ′′ i,ε (u x i ) ≤ N ((p -1) |∇u| p + |∇u| 2 ) ≤ N p -1 + 2 p |∇u| p + N p -2 p .
In the end we get (5.10)

G 2 ≥ -τ C ′ 1 ζ G 1 - C ′ 2 2 τ |∇u| p - C ′ 2 2 τ , where C ′ 1 = C ′ 1 (C, R 0 -r 0 ) > 0 and C ′ 2 = C ′ 2 (p, N, C, R 0 -r 0 ) > 0. The term G 3 .
By using the form of g i,ε , the convexity of the map m → m p-2 and recalling the definition (2.1) of δ, we have

G 3 = N i=1 g ′′ i,ε (u x i ) |u x i | 2 ≥ (p -1) 1 2 p-3 N i=1 |u x i | p -δ p-2 |∇u| 2 .
By further applying Young inequality to estimate the term |∇u| 2 and using that

N i=1 |u x i | p ≥ N 2-p 2 |∇u| p , we end up with (5.11) G 3 ≥ C ′′ 1 |∇u| p -C ′′ 2 , where C ′′ i = C ′′ (p, N, δ) > 0, i = 1, 2. The term G 4
This is the most delicate term and it is precisely here that the condition p ≥ 4 becomes vital. First we have

G 4 ≥ - N i=1 |g ′′′ i,ε (u x i )| |u x i x i | 2 λ |u| |u x i | + |∇u| 2 |ζ x i | .
Then we observe that by Cauchy-Schwarz inequality (recall the definition (2.6) of g i,ε ) we have

|∇u| 2 N i=1 |g ′′′ i,ε (u x i )| |u x i x i | |ζ x i | ≤ c |∇u| 2 N i=1 (|u x i | -δ i ) p-2 + u 2 x i x i |ζ x i | 2 1 2 × N i=1 (|u x i | -δ i ) p-4 + 1 2 ≤ c N i=1 (|u x i | -δ i ) p-2 + u 2 x i x i |ζ x i | 2 1 2 |∇u| p 2
for some constant c > 0 depending on p and N only. In the last inequality we used that

(|u x i | -δ i ) p-4 + ≤ |u x i | p-4 ,
since p ≥ 4. By further using (5.6), the definition of g i,ε and Young inequality, from the previous inequality we get

|∇u| 2 N i=1 |g ′′′ i,ε (u x i )| |u x i x i | |ζ x i | ≤ c (ζ G 1 ) 1 2 |∇u| p 2 ≤ τ C ′′′ 1 ζ G 1 + C ′′′ 1 τ |∇u| p ,
for some constant C ′′′ 1 = C ′′′ 1 (C, N, p, R 0 -r 0 ) > 0. Similarly, we have Collecting all the estimates.

2 λ N i=1 |g ′′′ i,ε (u x i )| |u x i x i | |u| |u x i | ≤ c λ M N i=1 g ′′ i,ε (u x i ) u 2 x i x i 1 2 N i=1 (|u x i | -δ i ) p-4 + |u x i | 2
We now go back to (5.7) and use (5.8), (5.10), (5.11) and (5.12). Then we get

L[ζ |∇u| 2 + λ u 2 ] ≥ 2 -τ 2 C ′ 1 + C ′′′ 1 ζ G 1 -C ′′′ 2 λ 2 p p+2 G p p+2 1 + 2 λ C ′′ 1 - C ′ 2 τ -2 ζ N p - C ′′′ 1 τ -C ′′′ 2 |∇u| p -2 M λ f L ∞ (2 B) + 2 ζ N p ′ ∇f p ′ L ∞ (2 B) + C ′ 2 τ + 2 λ C ′′ 2 .
We now choose τ small enough, in order to make the coefficient of ζ G 1 strictly positive. Then we also choose λ ≫ 1 large enough, so that |∇u| p as well has a strictly positive coefficient. Observe that the choices of τ and λ only depend on the relevant data of the problem and are in particular independent of ε. By setting for simplicity λ 2 p/(p+2) = Λ 2 , this gives (5.13)

L[ζ |∇u| 2 + λ u 2 ] ≥ c 1 ζ G 1 -c 2 Λ 2 G p p+2 1 + c 1 |∇u| p -c 2 ,
where c 1 > 0 and c 2 > 0 are constants that depend only on f W 1,∞ (2 B) , R 0 -r 0 , N, p, δ and the constant M appearing in (2.12).

Let us now consider the maximum of the function ζ |∇u| 2 + λ u 2 in B R 0 . If this maximum is assumed at x 0 ∈ ∂B R 0 , then we get max Thus at the maximum point x 0 we have

L[ζ |∇u| 2 + λ u 2 ] = N i=1 g ′′ i,ε (u x i ) ζ |∇u| 2 + λ u 2 x i x i ≤ 0.
By combining this with (5.13), we then get where the constant T q is given by T q = (2 -q) 2 4 q 2 -(2 -q) 2 > 0.

c 2 ≥ c 1 ζ G 1 -c 2 Λ 2 G p p+2 1 + c 1 |∇u| p .
Proof. We first observe that for every α, β > 1, by basic calculus we have By Hölder inequality we then have

|u(x 1 , x 2 )| α = α ˆx1 -∞ u x 1 (t,
ˆR2 |u x 1 | |u| α-1 dx ˆR2 |u x 2 | |u| β-1 dx ≤ ˆR2 |u x 1 | 2 dx 1 2 ˆR2 |u x 2 | q dx 1 q × ˆR2 |u| 2 (α-1) dx 1 2
ˆR2 |u| q q-1 (β-1) dx q-1 q

.

We now choose α and β in such a way that 2 (α -1) = α + β and q q -1 (β -1) = α + β, that is α = q + 2 2 -q and β = 3 q -2 2 -q .

Observe that with these choices we have α + β = 4 q/(2 -q). Thus from (B.2) we get (B.1), with

T q = 1 α 1 β = (2 -q) 2 4 q 2 -(2 -q) 2 ,
as desired.

1 . Introduction 1 . 1 .

 111 Overview. This paper is devoted to prove Lipschitz continuity for local minimizers of the anisotropic functional(1.1)

(1. 6 )

 6 D 2 G(z) ξ, ξ ≥ c |z| p-2 |ξ| 2 , |z| ≥ m,and the growth condition (1.7) |DG(z)| ≤ C |z| p-1 , |z| ≥ m.

2. 1 .

 1 Definitions and basic results. Let Ω ⊂ R N be an open set and p ≥ 2. In what follows we set for simplicity

Proposition 4 . 1 (

 41 Uniform Lipschitz estimate, N = 2). Let N = 2 and p ≥ 2. Then for every triple of concentric balls B r 0 ⋐ B R 0 ⋐ B ρ 0 ⋐ B and i = 1, 2 we have

≥ W 2 s p 1 ,since W i ≥ 1 .

 11

2

 2 |∇u| p , for some constant C ′′′ 2 = C ′′′ 2 (N, p, M ) > 0. By keeping everything together, we get (5.12)G 4 ≥ -τ C ′′′ 1 ζ G 1 -

ζ 0 ζ |∇u| 2 + λ u 2 =

 022 |∇u| 2 + λ u 2 ≤ max B R ζ(x 0 ) |∇u(x 0 )| 2 + λ u(x 0 ) 2 ≤ λ u 2 L ∞ (B R 0 ) ≤ λ M,thanks to the fact that ζ = 1 on B r 0 and ζ ≡ 0 on ∂B R 0 . This would prove the local Lipschitz estimate.In order to conclude, let us now assume that x 0 ∈ B R 0 , then we get∇ ζ |∇u| 2 + λ u 2 = 0 at x = x 0 ,andD 2 ζ |∇u| 2 + λ u 2 ≤ 0 at x = x 0 .

2 G 1 p2 G 1 - 2 τ -p 2 p + 2 Λ 2 + 2 , 0 |∇u| 2 ≤ max Br 0 ζ |∇u| 2 + 2 ≤

 211222202022 We multiply the previous by ζ(x 0 ) p/2 > 0, then by Young inequality once again we getc 2 ζ(x 0 ) p 2 ≥ c 1 ζ(x 0 ) p+2 2 G 1 -c 2 Λ 2 ζ(x 0 ) p+2 p+2 + c 1 ζ(x 0 ) p 2 |∇u(x 0 )| p ≥ c 1 -c 2 p p + 2 τ Λ 2 ζ(x 0 ) p+2 c 1 ζ(x 0 ) p 2 |∇u(x 0 )| p .and use that ζ ≤ 1, we finally get(5.14) ζ(x 0 ) |∇u(x 0 )| 2 ≤ C, with C = C(N, p, δ, f W 1,∞ (2 B) , M, R 0 -r 0 ) > 0.By using this bound we get maxBr λ u 2 ≤ ζ(x 0 ) |∇u(x 0 )| 2 + λ u(x 0 ) C + λ u(x 0 ) 2 ≤ C + λ M 2 ,which gives the desired conclusion.

x 2 )

 2 |u(t, x 2 )| α-2 u(t, x 2 ) dt, and|u(x 1 , x 2 )| β = β ˆx2 -∞ u x 2 (x 1 , s) |u(x 1 , s)| β-2 u(x 1 , s) ds. Thus |u(x 1 , x 2 )| α+β ≤ α β ˆR |u x 1 (t, x 2 )| |u(t, x 2 )| α-1 dt ˆR |u x 2 (x 1 , s)| |u(x 1 , s)| β-1 ds .If we now integrate over R 2 and use Fubini Theorem on the right-hand side, we get(B.2) ˆR2 |u| α+β dx ≤ α β ˆR2 |u x 1 | |u| α-1 dx ˆR2 |u x 2 | |u| β-1 dx .

Appendix A. Some properties of the functions g i The functions g i have the following convexity property.

Lemma A.1. For every t 1 , t 2 ∈ R such that |t 1 -t 2 | > 2 δ i we have

Proof. The function g i is convex so that the inequality ≤ holds true for every t 1 , t 2 . If g i ((1 -s) t 1 + s t 2 ) = (1 -s) g i (t 1 ) + s g i (t 2 ), then g i is affine on the segment [t 1 , t 2 ]. This can only happen when

They also satisfy the following Lipschitz-type estimate.

Lemma A.2. Let p ≥ 2. For every t 1 , t 2 ∈ R and i = 1, . . . , N , we have

Proof. By basic calculus we have

Thus we get the conclusion.

The following basic estimate has been used various times.

Lemma A.3. Let p ≥ 2. For every i = 1, . . . , N and every T ≥ δ i , we have

for every |t| ≥ T.

Proof. For T = δ i there is nothing to prove, thus we can suppose that T > δ i . We use the elementary inequality

This implies that for every |t| ≥ T , we have

By multiplying everything by (T -δ i )/T and raising to the power p -2, we get the desired conclusion.

Appendix B. An anisotropic Sobolev inequality in dimension 2

In the proof of Proposition 4.1 we used Sobolev-Troisi inequality. For the reader's convenience, we give a proof of the particular case we needed.

Lemma B.1. Let 1 < q < 2, then for every u ∈ C ∞ 0 (R 2 ) we have

,