Lamia Belguerras
email: lamia.belguerras@univ-lorraine.fr

Larbi Hadjout
email: lhadjout@yahoo.fr

Smail Mezani
email: smail.mezani@univ-lorraine.fr

Thierry Lubin
email: thierry.lubin@univ-lorraine.fr

A Rezzoug

Study of HTS Magnetic Coupler Using Analytical and Numerical Computations

Keywords: Analytical solution, BSCCO tape, flux concentration, high temperature superconductors, magnetic coupler, permanent magnet

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

agnetic Couplings (or couplers) (MCs) are widely used in industry to transmit a torque from a prime mover to its load without contact through a separation wall which can be air, vacuum, fluid or other media. Such contactless torque transmission devices have found wide applications in nuclear environments, vibration isolation, compressors, blowers, vacuum pumps [START_REF] Cong | Design and Application of Magnetic Coupling Used for Ultra-high Vacuum Robot[END_REF]. There are three types of MCs: Eddy current, Hysteresis and Synchronous, each of which may be of Concentric or Axial topology.

The use of High Temperature Superconductors (HTS) in electromagnetic devices is currently an important subject of research worldwide. Indeed, these materials can be efficiently used for developing high magnetic flux density in machines with large air gap as couplers and gears [START_REF] Jian | A Novel Coaxial Magnetic Gear Using Bulk HTS for Industrial Applications[END_REF][START_REF] Okano | Characteristics of the magnetic gear using a bulk high-Tc superconductor[END_REF].

In the same context, a flux concentration (or flux focusing) concept is an effective solution since the machine produces more torque compared to the ones with surface mounted PMs [START_REF] Enokizono | Development of high density and high efficiency machines[END_REF].

The aim of this paper is to use HTS and PMs to transmit high torque from a cold environment to a warm environment using a superconducting magnetic coupler. For example, the torque produced by a HTS motor can be transmitted to its load via the HTS magnetic coupler of Fig. 1.a. This may constitutes an alternative solution to torque tubes. A 3D view of the studied flux concentration coupler is shown in Fig. 1.b.

Obviously, this doesn't constitute the sole application of a magnetic coupler. Many industrial processes (i.e. nuclear, chemical, pharmaceutical and food industries) need hermetic isolation between the prime mover and a follower connected to the process.

Study of HTS Magnetic Coupler Using Analytical and Numerical Computations

Lamia Belguerras, Larbi Hadjout, Smail Mezani, Thierry Lubin, and Abderrezak Rezzoug The proposed HTS coupler (or a wholly HTS coupler whose both inner and outer rotors are superconducting) may be the unique solution if large distance between the two rotors is required (several centimeters). In fact, superconductors are the only materials which can create enough magnetic field to transmit a torque without contact, and within a reasonable coupler volume, throw a thick wall.

Usually, the magnets of flux concentration permanent magnets rotor are rectangular (as in Fig. 1.b). The distribution of magnetic field can be calculated by dividing the rectangular magnets into a finite number of radial slices. However, this will drastically increase the number of sub-domains and penalizes the computation time. Another solution consists of using conformal mapping to handle the non radial areas of slot domains [START_REF] Wu | Comparison of Analytical Models for Predicting Cogging Torque in Surface-Mounted PM Machines[END_REF].

In this study, the rectangular magnets are replaced by sector ones having the same volume so the magnetic field distribution can be calculated by using analytical method. In this technique, polar coordinates are used and Maxwell's equations are solved by separation of variables in low permeability sub-domains (air regions, magnets, windings) [START_REF] Liu | Analytical solution of air-gap field in permanent magnet motors taking into account the effect of pole transition over slots[END_REF][START_REF] Zhu | An accurate subdomain model for magnetic field computation in slotted surface-mounted permanent magnet machines[END_REF][START_REF] Gu | Effect of slotting in PM electrical machines[END_REF][START_REF] Lubin | Exact analytical method for magnetic field computation in the air-gap of cylindrical electrical machines considering slotting effects[END_REF]. In particular, flux concentration switched PMs machines have been studied in [START_REF] Gysen | Modeling of Flux Switching Permanent Magnet Machines with Fourier Analysis[END_REF].

The first part of this paper introduces the sector PMs FCSMC and the analytical calculation of the magnetic field. The results are validated by finite element (FE) computation.

In the second part, we propose a method to calculate the critical current in the superconducting coil by taking into account the HTS tape characteristics and the flux density on the HTS wire. This flux density is calculated by the analytical and the numerical techniques.

In the last part, we will analyze the capability of the developed analytical model in predicting critical current of HTS coil and rectangular magnets FCSMC performances. A discussion on the obtained results ends the paper.

II. STUDIED SUPERCONDUCTING MAGNETIC COUPLER

The considered coupler is coaxial (or concentric); it has one rotor inside another rotor. The inner rotor consists of different pieces of iron and permanent magnets that are fixed jointly on a non-ferromagnetic shaft; each magnet is magnetized in a reverse direction of its neighboring ones. The outer rotor is composed of a series of HTS coils separated by pieces of glass fiber. The coils are cooled inside a non magnetic cryostat (stainless steel) at 30K (Fig. 2).

The geometrical parameters of prototype FCSMC to be investigated (Fig. 3) are -The inner radius of the inner rotor R 1 -The outer radius of the inner rotor R 2 -The inner radius of the HTS winding R 3 -The outer radius of the HTS winding R 4 -The inner radius of outer rotor yoke R 5

The inner rotor consists of 2p (p is the number of polepairs) sector PMs. The angular position of the i-th PM is defined (Fig. 2) by

p i i       2 with p i 2 1   (1)
Where  is the PM opening angle.

III. PROBLEM FORMULATION AND ASSUMPTIONS

A 2D analytical model is developed where a magnetic vector potential formulation is used in polar coordinates. The relative permeability is considered infinite for the ferromagnetic parts and is equal to 1 for the PMs and air regions. The inner rotor is in a fixed position for all the calculations.

A. Source Terms

The permanent magnets buried into the rotor of the machine are tangentially magnetized. The magnetization direction is alternated from one magnet to another. Hence, the magnetization vector of the i th permanent magnet regions is

  e B M r i 0 ) 1 - (  (2)
where B r is the remanent flux density of the magnet.

The HTS coils in the outer rotor are supplied by direct current (Fig. 4). The Fourier series expansion of the current density J of the HTS winding is given by

  odd n , cos ) ( 1      n n n J J    (3) 1 0 , 2 n sin 4 0          w w n J n J     (4)
Where J 0 is the dc current density in the superconductor and  the load angle (angle between the magnetic fields created by the HTS rotor and the PM rotor).  w is the ratio between the HTS coil cross section width (R 3  w π/2p) and half of the pole pitch (R 3 π/2p).

B. Analytical expressions of the magnetic field in the different subdomains

The studied structure is divided into a shaft (domain I), airgap (domain II), HTS winding (domain III), outer cryostat wall (domain IV) and permanent magnet regions (domains i), Fig. 3. Notice that the air gap includes the inner wall of the cryostat (non magnetic) and a mechanical clearance to allow the rotation.

In order to simplify the analytical expressions in the different sub-domains, we adopt the following notations

z z z a b b a ) b , a ( P               (5) z z z a b b a ) b , a ( E               (6)
In polar coordinates, for the i-th permanent magnet subdomain region and in the HTS winding domain, Poisson's equations apply

r B A r i i ) 1 (    (7) J A III 0     (8)
For the sub-domains I, II and IV, Laplace's equations have to be solved

0   I A (9) 0   II A (10) 0   IV A (11)

C. Boundary conditions

The boundary conditions for the i-th PM domain (Fig. 5) are

0     i i A    and 0          i i A (12)
The continuity conditions of the magnetic vector potential between the i-th PM and sub-domains I, II are

) , R ( A ) , R ( A I i   1 1  ( 13
)
) , R ( A ) , R ( A II i   2 2  ( 14
)
The boundary condition for the region III is given by

3 3 R r III R r II r A r A        (15)
The boundary conditions for the region IV are given by The continuity conditions of the magnetic vector potential between the sub-domain III and sub-domains II, IV are

) , R ( A ) , R ( A II III   3 3  ( 18
)
) , R ( A ) , R ( A IV III   4 4  (19)
0 θ A i θ i    ) g( ) f(  IV R 4 R 5 R 1 R 2 R 3 II 𝜃 I PM Domain i r PM (i) HTS rotor III ) , R ( A ) , R ( A III II   3 3  0 θ A i θ i      ) , R ( A ) , R ( A III IV  

IV

The continuity of the tangential magnetic field at the interface between the permanent magnet and the shaft, can be written as

) ( f r A R r I      1 (20) with                elsewhere if B r A f i i i 0 , ) 1 ( ) ( i r R r 1      (21)
The continuity of the tangential magnetic field at the interface between the permanent magnet and the air gap can be written as

) ( g r A R r II      2 (22) with                elsewhere if B r A g i i i 0 , ) 1 ( ) ( i r R r 2      (23)

IV. ANALYTICAL SOLUTIONS

To solve the partial differential equations, the method of separation of variables is used in cylindrical coordinates.

A. General Solution of Posson's equation in the i-th PMs subdomain (region i)

In the i th PM sub-domain, the problem to solve is

r B A r r A r r A r i i i i ) 1 ( 1 1 2 2 2 2 2                         i i R r R 2 1 (24)
The general solution of (24) can be written as

           1 2 1 / 1 / 2 1 / 2 / 0 0 ) ( cos( ). ) , ( ) , ( ) , ( ) , ( (
ln ) 1 ( ) , ( k i k k i k k k i k i i r i i k R R E R r E B R R E R r E A r B A B r r A              (25)
where k is a positive integer, ( ) is given by ( 6).

The constants , , and are determined using a Fourier series expansion of shaft and air-gap vector potentials ( 𝜃), ( 𝜃) over the PMs region interval [ i , 

i +]. They are given by

            i i d R A B R R B A I r i i i )) , ( ) 1 ( ( 1 ln 1 1 1 0 0 (26)             i i d R A B R R B A II r i i i )) , ( ) 1 ( ( 1 ln 2 2 2 0 0 (27)                   i i d k R A B R A i I r i i k )) ( cos( ) , ( ) 1 ( 2 1 1 (28)                   i i d k R A B R B i II r i i k )) ( cos( ) , ( ) 1 ( 2 2 2
(29)

B. General solution of Laplace's equation in the shaft subdomain (region I)

The shaft sub-domain and its associated boundary conditions are shown in Fig. 5. The equation to solve is

0 1 1 2 2 2 2 2           I I I A r r A r r A for          2 0 0 1 R r ( 30
)
By taking into account the boundary condition (20), the general solution of the magnetic vector potential in the shaft can be written as

)) n sin( B ) n cos( A ( R r n R ) , r ( A I n n I n n I                1 1 1 (31)
Where n is a positive integer, the coefficients and are determined by using a Fourier series expansions of (𝜃)

over the interval [0, 2].        2 0 2 2 d ) n cos( ) ( f A I n (32)        2 0 2 2 d ) n sin( ) ( f B I n (33)

C. General solution of Laplace's equation in the air gap (region II)

In the air-gap sub-domain, the problem to solve is

0 1 1 2 2 2 2 2           II II II A r r A r r A for          2 0 3 2 R r R ( 34
)
According to the boundary conditions ( 15) and ( 22), the general solution of the magnetic vector potential in the air-gap sub-domain can be written as

                       1 2 3 2 2 3 2 3 2 1 2 3 2 3 3 2 3 2 n n n II n n n II n n n n II n n n II n II n sin ) ) R , R ( ) R , r ( n R D ) R , R ( ) R , r ( n R C ( n cos ) ) R , R ( ) R , r ( n R B ) R , R ( ) R , r ( n R A ( ) , r ( A    (35)
Where n is a positive integer, ( ) and (

) are defined by ( 5) and ( 6). The coefficients , , and are determined using a Fourier series expansions of

       2 0 2 2 d ) n cos( ) ( g A II n (36)          2 0 3 2 2 d ) n cos( r A B R r III II n (37)      d ) n sin( ) ( g C II n   2 0 2 2 (38)          2 0 3 2 2 d ) n sin( r A D R r III II n (39)

D. General Solution of Poisson's equation in the HTS winding sub-domain (region III)

The problem to solve in the HTS winding region is given by the following equation

     J A r r A r r A III III III 0 2 2 2 2 2 1 1           for          2 0 4 3 R r R (40)
By taking into account the boundary conditions ( 18) and [START_REF] Barth | Influence of Shear Stress on Current Carrying Capabilities of High Temperature Superconductor Tapes[END_REF], the general solution for the magnetic vector potential is

                                            
n r R R R r E D R R R r E C n r R R R r E B R R R r E A r A        (41) With               3 3 4 4 4 3 3 4 3 4 4 3 4 4 3 3 4 4 2 R X R r R X R r r X R X R X R R ) R , R ( ) R , r ( E R X R X R R ) R , R ( ) R , r ( E ) r ( n n n n n n n n n n n n n n n n                                                           (42) and                     otherwhere 0 1,3,5... i , ip n if 16 2 n if 1 r ln 4 16 ) ( 2 0 2 0 r J r J r X n n n   (43)
Where n is a positive integer, p is the number of pole-pairs and (

) is defined by [START_REF] Liu | Analytical solution of air-gap field in permanent magnet motors taking into account the effect of pole transition over slots[END_REF].

The coefficients , , and are determined using a Fourier series expansions of A II (R 3 ,θ) and A IV (R 4 ,θ) over the air-gap interval [0, 2].

       2 0 3 2 2 d ) n cos( ) , R ( A A II III n (44)        2 0 4 2 2 d ) n cos( ) , R ( A B IV III n (45)        2 0 3 2 2 d ) n sin( ) , R ( A C II III n (46)        2 0 4 2 2 d ) n sin( ) , R ( A D IV III n (47)

E. General Solution of Laplace's equation in the outer cryostat wall sub-domain (region IV)

The equation to solve is

0 1 1 2 2 2 2 2           IV IV IV A r r A r r A for          2 0 5 4 R r R (48)
According to the boundary conditions ( 16) and ( 17), the general solution of the magnetic vector potential in the outer cryostat wall sub-domain can be written as

                    1 5 4 5 4 5 4 5 4 n n n IV n n n IV n IV n sin ) R , R ( ) R , r ( P n R B n cos ) R , R ( ) R , r ( P n R A ) , r ( A    (49)
Where n is a positive integer, (

) and ( ) are defined by ( 5) and [START_REF] Liu | Analytical solution of air-gap field in permanent magnet motors taking into account the effect of pole transition over slots[END_REF]. The coefficients and are determined using a Fourier series expansions of

         2 0 4 2 2 d ) n cos( r A A R r III IV n (50)          2 0 4 2 2 d ) n sin( r A B R r III IV n (51)

F. Electromagnetic torque

The electromagnetic torque is obtained using the Maxwell stress tensor. A circle of radius R e in the air-gap is taken as the integration path so the electromagnetic torque is expressed as follows Using the coefficients of the magnetic field solution (35) in the air gap (region II), the electromagnetic torque expression (52) can be rewritten as

        2 0 0 2 d ) , R ( B ) ,
) Z Y X W ( LR T n n n n n e e     1 0 2   (53)
where In order to validate the proposed model and for the geometrical parameters given in Table I, the analytical results have been compared with 2-D finite-element simulations obtained using FEM software [START_REF] Meeker | Finite Element Method Magnetics Version 4[END_REF]. The chosen parameters lead to non-saturated iron parts of the coupling so the linearity assumptions adopted for the analytical model are fulfilled. Fig. 6 shows the magnetic flux lines under no-load (φ=0°) and load (φ=57°) conditions. The flux lines at no-load are symmetrical but it's not anymore the case under load conditions.

) R , R ( ) R , R ( R R B ) R , R ( ) R , R ( R R A W
       ( 54
) ) R , R ( ) R , R ( R R D ) R , R ( ) R , R ( R R C X n e n e II n n e n e II n n 2 3 2 3 3 2 3 2        ( 55
) ) R , R ( ) R , R ( R R D ) R , R ( ) R , R ( R R C Y n e n e II n n e n e II n n 2 3 2 3 3 2 3 2       ( 56
) ) R , R ( ) R , R ( R R B ) R , R ( ) R , R ( R R A Z
The waveforms of the flux density distribution (radial and tangential components) in the middle of the air-gap are presented in Fig. 7 and Fig. 8. It can be seen that the radial and tangential flux density waveforms under load conditions are deformed when compared to the no-load situation.

The static torque as well as the reluctance torque vs. the load angle  are presented in Fig. 9. Notice that the reluctance torque has been computed by turning off the PMs, viz. B r =0.

For the considered geometrical parameters, the pull-out torque is obtained for a load angle nearly equals to =57°.

The reluctance torque represents approximately 10% of the total torque (Fig. 9

.b).

A very good agreement is noticed between the finite element predictions and the analytical results.

VI. THE HTS WINDING

This section is devoted to the analysis of the HTS winding properties. It is well known that HTS materials exhibit a non linear relation between the electric field E and the current density J. This non-linearity is well described by the power law [START_REF] Bruzzone | The index n of the voltage-current curve, in the characterization and specification of technical superconductors[END_REF]

n c c J J E E          ( 58
)
Where J c is the critical current density and n the index value of the experimental U(I) curve [START_REF] Zhang | Investigations on N Values of Long Bi2223/Ag Tapes[END_REF].

E c is the critical electric field corresponding to the critical current with the arbitrary criteria E c =10 -4 V/m. J is the current density in the coil.

J c and n depend on the cooling temperature, on the magnetic flux density components and on some mechanical constraints. So, the electric loadability (maximal allowable current avoiding thermal limits) also depends on the same physical properties. At a given operating temperature of the HTS material, the critical current density only depends on the electromagnetic and mechanical constraints. Nevertheless, one has to care about local hot-spot temperatures in the HTS winding. In fact, important local loss densities (due to a non-homogeneous cooling or a failure of the cooling system) can cause the quench of the HTS winding [START_REF] Lee | Thermal Quench in HTS Double Pancake Race Track Coil[END_REF], [START_REF] Shoykhet | Quench in High-Temperature Superconducting Motor Field Coils: Computer Simulations and Comparison With Experiments[END_REF].

Obviously, these thermal issues are of great importance but they are beyond the scope of the present paper so they will not be discussed further.

A. Influence of the flux density distribution on the electric loadability of the HTS winding

For the considered device, the HTS winding consists of 2p race-track coils made from BSCCO tapes (Fig. 10). Each coil consists of N turns connected in series. The critical current of the coil is defined as the current flowing in the coil which causes a coil's terminal voltage equals to its critical value.

Since a 2D electromagnetic model is considered here, we neglect in this study the influence of the end-winding on the coil's voltage. So we only take into account the active length of the coupler. To describe the dependence of Jc and n values on the flux density components, the Kim model is used [START_REF] Kim | Critical persistent currents in hard superconductors[END_REF]  

               0 2 2 2 0 1 B B B k J B J // c c (59)   0 0 1 n B B n B n    ( 60
)
Where J c0 , n 0 and k are respectively the critical current density in zero field, the n value in zero field of the SC tape and the anisotropy factor (<1). B 0 , B n0 and β are positive constants characterizing the material at the given cooling temperature.

and are the values of the component of the flux density perpendicular and parallel to the largest surface of the tape (Fig. 10).

Equations ( 59) and ( 60) show that the J c and n values are mainly affected by the perpendicular flux density.

To rigorously determine the voltage over the coil's terminals, the actual stranded topology of the winding cross section has to be considered. The total voltage of the coil is then the sum of the elementary voltages corresponding to each turn of the HTS winding. This strict computation is not feasible using the developed analytical model and becomes very complicated by FE analysis since the generated mesh growth rapidly when the number of turns increases.

To overcome the problem, we use a bulk topology of the coil for which we define an average current density (as in (3)) given by

coil coil turns S I N J  ( 61
)
where I coil is the current flowing in each conductor and S coil is the cross section area of the coil.

This approach has been successfully used in our lab [START_REF] Bendali | Dimensionnement d'un moteur supraconducteur HTc[END_REF] since measured U coil (I coil ) characteristics are consistent with those computed using the bulk approximation.

In polar coordinates, the voltage across the coil is the sum of the voltages of the "in" and "out" of the plane bundles of conductors.

                out in S out S in coil turns coil d d r r E d d r r E S N L U     r , r , (62)
Where E in (r, θ) and E out (r, θ) are respectively the electric fields over the cross sections of the coil (S in ) and (S out ). The double integral in (62) is calculated numerically. By using the power law (58) associated to the Kim model ( 59) and (60), equation (62) becomes

) FE Analytical                                              out in S B n c S B n c w turns coil rdrd B B J J rdrd B B J J R R N p U     ) ( // ) ( // 2 3 2 4 c , , E L 2 (63)
For the studied device and using the analytical model, the parallel and perpendicular flux densities are respectively and (Fig. 10).

The critical voltage U c over the coil's terminals is (64)

B. Mechanical stress due to bending HTS materials properties are very sensitive to mechanical stresses like bending, torsion or shear [START_REF] Oomen | Transposed-Cable Coil & Saddle Coils of HTS for Rotating Machines: Test Results at 30 K[END_REF][START_REF] Barth | Influence of Shear Stress on Current Carrying Capabilities of High Temperature Superconductor Tapes[END_REF]. In addition to their harmful effects on the mechanical integrity of HTS materials, these stresses irreversibly decrease the critical current density.

The heads of race track coils should not be too short as they are subjected to bending and torsion during their manufacture. A minimal bending diameter D b is then required which constitutes a geometrical constraint on the coil opening (Fig. 10). Hence, the HTS coil width to half pole pitch ratio ( w ) has to fulfill the following condition

3 1 R D p b w     (65)
The race track coil being made from rectangular tapes of thickness e w , the number of layers N ct in the tangential direction  is limited to a maximal value of

w w f ct e R p K N 3 2    (66)
where is the filling factor (in the tangential direction) of the HTS coil.

C. Critical current of the HTS coupling

The studied coil is made of BSCCO tapes cooled at 30K. The electrical and geometrical characteristics of the BSCCO tape are given respectively on Tables II and III. These characteristics correspond to the last generation BSCCO tape (of type H for dc applications) manufactured by Sumitomo [20].

For K f =0.8 at no-load (φ=0°), Fig. 11 shows the coil's terminals voltage vs. the current. The critical voltage of one BSCCO coil is U c =21 mV. The critical current obtained for the studied HTS coupling is about 277A, which represents 52% of the self-field critical current of the used BSCCO tape (I c0 =531A). The analytical predictions are in good agreement compared to the numerical computations. Fig. 12 presents the critical current of the BSCCO coil for different load angles φ of the superconducting rotor. It can be seen that the maximum value of the critical current is about 462A (0.87 I c0 ) obtained for a load angle nearly equals to φ=33° (electrical angle). For the pull-out torque position corresponding to φ=57°, the critical current is about 441 A (0.83 I c0 ).

The minimum value of the critical current, obtained for φ=0°, is about 277 A (0.52 I c0 ). This limit allows fixing the operating current of the coupler at a value lower than this minimal critical current in order to prevent thermal failures or to allow some mechanical over-load. So from now, the rotor of the studied coupler will be composed of rectangular magnets embedded in a ferromagnetic yoke. Indeed, we propose an analysis of the capability of the developed analytical model in predicting the performances of FCSMC with rectangular PMs.

Notice that superconductors create a strong magnetic field which saturates the iron parts, in particular in the inner rotor. The FE computations are then performed both in the linear and in the non-linear cases. The ferromagnetic material used in this study (US Steel Type 2-S) has a B(H) curve shown in Fig. 13.

To allow the analytical computation, the rectangular magnets are replaced by sector ones having the same volume (Fig. 14). The rectangular magnet opening at r=R 1 , noted β', couldn't exceed p R / 1  . We define the ratio between R 1 and R 2 as

2 1 R R r   (67)
The sector PM, having the same volume as the rectangular one, has to fulfill the following condition

          p p r r m 2 sin 1 4      (68)
For all what follows and in order to predict the performances of FCSMC with rectangular PMs, linear FE, non-linear FE and analytical computations are performed and compared.

The following values have been adopted: Fig. 15 presents the voltage-current curve of the HTS coil at no-load (φ=0°). The critical current obtained using the analytical model is about 277 A; this value is very close to the one computed using the linear FE model. The non-linear FE computation gives a critical current of about 272 A which is 1.8 % lower than the one given by the analytical model. These good predictions suggest that the critical current is mainly limited by the self field created by the HTS winding current. Indeed, for reasonable saturation levels, this current gives the same flux density distribution across the coil area whatever the geometrical configuration (rectangular or sector PMs).

Fig. 16 shows the flux lines under load conditions (φ=57°,

I coil =80%I c =222.4 A).
We notice some leakage flux in the shaft which mainly depends on the values of α r and α m .

As shown in Fig. 17, the air gap radial flux density harmonics are globally higher with the rectangular PMs configuration. For the most important harmonic orders [START_REF] Cong | Design and Application of Magnetic Coupling Used for Ultra-high Vacuum Robot[END_REF][START_REF] Okano | Characteristics of the magnetic gear using a bulk high-Tc superconductor[END_REF][START_REF] Wu | Comparison of Analytical Models for Predicting Cogging Torque in Surface-Mounted PM Machines[END_REF][START_REF] Zhu | An accurate subdomain model for magnetic field computation in slotted surface-mounted permanent magnet machines[END_REF][START_REF] Lubin | Exact analytical method for magnetic field computation in the air-gap of cylindrical electrical machines considering slotting effects[END_REF], the relative difference between the non-linear FE computations and the analytical predictions doesn't exceed 7%.

Position of pull-out torque

For a current value I coil =222.4 A, Fig. 18 presents the static torque vs. load angle computed by the different methods in use. The linear FE prediction gives a pull-out torque of about 198 kNm while the non-linear FE computation gives 195 kNm. The analytical model underestimates the pull out torque since it predicts 171 kNm. Furthermore, the load angle corresponding to the pull out torques is φ=55° using the FE models whereas this angle is 13° higher (viz. φ=68°) when using the analytical model. One can also note that the torque computed analytically is lower than the FE one whatever the load angle. Nevertheless, the analytical model allows a correct estimation of the pull-out torque.

A. Influence of p and α r on the critical current and the torque capabilities of FCSMC

By fixing R 2 =461 mm as in Table I, the critical current computed for different values of the number of pole pairs p and φ=57° is presented in Fig. 19. The other parameters remain unchanged; in particular  w is constant so the overall volume of the HTS material is also constant whatever the value of p. It can be seen from Fig. 20 that the FE models as well as the analytical model predict identical critical current values for p=2÷5. Beyond, the analytical model overestimates the numerical one by about 9 %. Furthermore, one can see that the critical current increases with p. This is due to the fact that the local values of the magnetic flux density over the HTS coil decrease as p increases.

Regarding the influence of the ratio  r , the calculations showed that it has no influence on the critical current values.

For a current corresponding to 80% of I c in the HTS windings, Fig. 20 presents the variation of the pull-out torque vs. the number of pole pairs p for 2 values of  r . It can be seen that the torque decreases monotonically as p increases, and the analytical computations follow the same tendency as the nonlinear FE ones.

In order to analyze the influence of the parameter  r , Fig. 21 shows the evolution of the pull-out torque when  r varies between 0.2 and 0.8. It is clear that an optimal value of  r witch maximizes the torque exists. This value is about  r =0.55 (170 kNm) for p=4 and  r =0.7 (130 kNm) for p=8.

B. Optimized torque transmission capabilities of the FCSMC

The parametric studies showed that optimal configurations in terms of torque transmission and electric loadabilty of the HTS winding exist. These studies have been conducted with a moderate active length of HTS tape. To increase the torque produced by the coupler, we need to use more HTS material as to increase the airgap magnetic field. By doing so, one can face some problems related to the magnets demagnetization. Hence, we have to ensure that the used NdFeB magnets are not demagnetized by systematically checking the values of the tangential magnetic field (which corresponds to the direction of the magnetization).

Several simulation studies (analytical and numerical) allowed us to obtain an optimized FCSMC whose main dimensions are presented in Table IV. The minimal value of the critical current is obtained at φ=0° (I c =305A).

The static torque curve, obtained for a HTS coil current I coil =0.8I c =243 A, is presented in Fig. 22. The pull-out torque reaches 425 kNm. A 5cm back-iron thickness is enough to close the flux which leads to an active torque density of about 450 kNm/m 3 .

As an example, this coupler is able to transmit the torque of a 5 MW-120 rpm ship propulsion motor. It has clearly the potential to be used in high power application systems like wind turbines, marine propulsion…

VIII. CONCLUSION

The paper presents a new type of torque transmitter which can replace favorably, by avoiding any material contact, a mechanical torque tube usually used in superconducting machines. The use of high critical current BSCCO tapes at 30K and flux concentration PMs rotor allows to obtain a high torque density magnetic coupler.

An exact 2-dimensional analytical method to calculate the magnetic field distribution in FCSMC with sector magnets has been proposed and validated through FE computations. The study takes into account the electrical and mechanical constraints of BSCCO tape at 30K. The results presented in this paper have shown that the developed model is able to predict the critical current of the BSCCO coil and the performances of the superconducting coupler with a good accuracy when compared to the FE predictions.

In the flux concentration machine the permanent magnets are usually rectangular, so we have been led to introduce an equivalent sector shape magnets in order to obtain analytically the results presented here. The different results showed that, by considering the HTS constraints, the analytical model developed in this paper predicts the value of the pull-out torque within a maximum error of about 20% when compared to the finite element calculations of the rectangular PM topology. Furthermore, different parametric studies have been performed. It has been shown that the analytical predictions follow the same tendency as the numerical ones in evaluating the electromagnetic torque and the critical current.

An optimized FCSMC has been also presented leading to a very important torque density coupler. Hence, it can be used in high power electromechanical systems.

Fig. 1 .

 1 Fig. 1. (a) Schematic diagram showing the use of a superconducting magnetic coupler, (b) 3D view of the studied flux focusing coupler.

Fig. 2 .

 2 Fig. 2. 2D geometry of the studied superconducting magnetic coupler (p=2)

Fig. 5 .

 5 Fig. 5. Boundary conditions

4 4 Fig. 4 .

 444 Fig. 4. Current density distribution in the HTS windings along the θ direction (φ = 0)

  g(θ) over the air-gap interval [0, 2].

  over the air-gap interval [0, 2].

  radial and tangential flux densities in the air gap and L is the active axial length of the FCSMC.

Fig. 6 .

 6 Fig. 6. Magnetic flux lines under (a) no-load (φ=0) and (b) load conditions (φ=60°)

Fig. 8 .Fig. 7 .

 87 Fig. 8. Radial (a) and tangential (b) flux densities in the middle of the air gap (Re=345 mm) under load conditions (φ=60 0 )

Fig. 10 .Fig. 9 .

 109 Fig. 10. 3D view part of a race-track BSCCO coil showing the tapes geometrical arrangement

Fig. 11 .

 11 Fig. 11. Voltage vs. current of HTS coil of FCSMC (φ=0°)

.

  These choices allow a reasonable saturation level in the iron parts.

Fig. 14 .Fig. 12 .

 1412 Fig.14. Replacing rectangular PMs by sector ones R 1

Fig. 17 .

 17 Fig. 17. Radial flux density waveforms (a) and their harmonic spectra (b) at load conditions (φ=57° elect, Icoil=222.4 A).

16 .

 16 PMs) FE linear (Rectangular PMs) FE non-linear (Rectangular PMs) Magnetic flux lines of FCSMC with (a) sector and (b) rectangular PMs under load conditions (φ=57°élec, Icoil=222.4 A).

Fig. 15 .

 15 Fig. 15. Voltage vs. current curves of the BSCCO coil for φ=0°

Fig. 21 .

 21 Fig. 21. Pull-out torque versus the ratio r for p=4 and p=8.

Fig. 19 .Fig. 18 .

 1918 Fig. 19. Critical current versus the number of pole pairs (r=0.6)

Fig. 22 .

 22 Fig. 22. Static torque of the optimized FCSMC computed with nonlinear FE (Icoil=243 A)

TABLE I GEOMETRICAL

 I PARAMETERS OF FCSMC WITH SECTOR PMS

		Symbol	Quantity	Value
		p	Number of pole pairs	5
		R5	Outer radius of SC coupling	500 mm
		L	Active axial length	1000 mm
		N	Number of harmonics used
			for magnetic field calculation	20
			in the shaft, air-gap, winding
			and outer cryostat domains
		K	Number of harmonics used
			for magnetic field calculation	30
			in the PM domain
		R4	Outer radius of HTS windings	490 mm
	HTS	R3	Inner radius of HTS windings	481 mm
	outer	αw	HTS coil width to half pole pitch	0.1
	rotor		ratio
		Nturns	Number of series turns	104
		Iccoil	Current of HTS coil	100 A
		LBSCCO	Total length of BSCCO tape	2 km
		R2	Outer radius of PM rotor	461 mm
	PM	R1	Inner radius of PM rotor	276 mm
	inner	β	Opening of PM	0.35 π/p rad
	rotor	Br	Remanence of the PM (NdFeB)	1 T

Lamia Belguerras was born in Algiers, Algeria, in 1979. She received the Dipl.-Ing. and the Magister degrees from the University of Sciences and Technology Houari Boumediene (USTHB), Algiers, Algeria, in 2005 and 2008, respectively. She is currently a Ph.D. student at GREEN Laboratory, University of Lorraine, Nancy, France. Her research interests include the modelling and optimization of electrical machines and superconducting magnetic couplings.

Larbi Hadjout was born in


Antoine Bonnefoy

Valentin Emiya

Liva Ralaivola

Rémi Gribonval

A Dynamic Screening Principle for the Lasso

Keywords: Screening test, Dynamic screening, Lasso, First-order algorithms, ISTA

come   L'archive ouverte pluridisciplinaire

Introduction

The Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] is an optimization problem that aims at finding a sparse solution to a least square problem by minimizing the sum of an ℓ 2 -fitting term and an ℓ 1 -regularization term. Given some observation/signal y ∈ R N and a dictionary matrix D ∈ R N ×K with N ≤ K, this problem writes P(λ, D, y) : x arg min

x 1 2 Dx -y 2 2 + λ x 1 , (1)
where the parameter λ > 0 governs the sparsity of x. We would like to be able to handle [START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF] when both N and K may be large, which occurs in many practical applications resorting to the Lasso: denoising, inpainting or classification. Algorithms relying on first-order information (e.g. gradient) only are particularly suited for these problems, as second-order based methods (e.g. using the Hessian) imply too computationally demanding iterations. These first-order algorithms include primal [START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF][START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF][START_REF] Daubechies | An Iterative Thresholding Algorithm for Linear Inverse Problems with a Sparsity Constraint[END_REF][START_REF] Wright | Sparse reconstruction by separable approximation[END_REF] and primal-dual [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Uzawa | Studies in linear and nonlinear programming, With contributions by Hollis B. Chenery[END_REF] algorithms.

Accelerating these algorithms is yet a key challenge: even though they provably have fast convergence, they remain captive of the dictionary size due to the required multiplications by D and D T over the optimization process. To overcome this limitation, strategies based on screening tests [START_REF] Dai | An ellipsoid based, two-stage screening test for bpdn[END_REF][START_REF] Ghaoui | Safe Feature Elimination in Sparse Supervised Learning[END_REF][START_REF] Tibshirani | Strong rules for discarding predictors in lasso-type problems[END_REF][START_REF] Wang | Lasso Screening Rules via Dual Polytope Projection[END_REF][START_REF] Xiang | Learning sparse representations of high dimensional data on large scale dictionaries[END_REF][START_REF] Xiang | Fast lasso screening tests based on correlations[END_REF] have recently been proposed. They implement two steps: i) locate zeros of x thanks to a screening test and construct the reduced or screened dictionary D 0 which is dictionary D trimmed off of its columns that correspond to the located zeros and ii) solve P(λ, D 0 , y) (see Algorithm 1).

We propose a new screening principle called dynamic screening in order to even more reduce the computational cost of first-order algorithms. We take the aforementioned concept of screening test one step further, and improve existing screening tests by embedding them in the iterations of first-order algorithms. We take advantage of the computation made during the optimization procedure to perform a screening at each iteration with a negligible computational overhead, and we consequently dynamically and iteratively reduce the size of D. To our knowledge, this is the first time such a screening mechanism is envisioned. Opposing perspectives of the proposed dynamic screening and existing static screening are schematized in Algorithms 1 and 2.

Experiments show that the dynamic screening principle significantly reduces the computational cost of the optimization in a large range of λ values. The computational saving reaches up to 90% with respect to the base algorithm, or up to 70% with respect to the algorithm with static screening.

Algorithm 1 Static screening strategy

D0 ← Screen D loop k x k+1 ←Update x k using D0 end loop Algorithm 2 Dynamic screening strategy D0 ← D loop k x k+1 ← Update x k using D k D k+1 ← Screen D k using x k+1 end loop
Section 2 introduces the tools we build our work upon. The dynamic screening principle is then presented and analyzed in Section 3. Section 4 is devoted to numerical experiments. Finally we discuss several extensions that can emerge from this work in Section 5.

In this section, we set the notation, introduce previous works on screening tests for the Lasso and recall state-of-the-art algorithms to solve this problem, pointing out their computational limitations.

Notation

D [d 1 , . . . , d K ] ∈ R N ×K
denotes a dictionary and Ω {1, . . . , K} denotes the set of integers indexing the columns, or atoms, of D. The i-th component of x is denoted x(i). The observation y ∈ R N is assumed to have a sparse representation x ∈ R K in D, i.e. Dxy 2 and x 0 are small. Without loss of generality, observation y and atoms d i are assumed to have unit ℓ 2 norm. The dual problem associated to (1) is [START_REF] Ghaoui | Safe Feature Elimination in Sparse Supervised Learning[END_REF][START_REF] Xiang | Learning sparse representations of high dimensional data on large scale dictionaries[END_REF]:

θ arg max θ 1 2 y 2 2 - λ 2 2 θ - y λ 2 2 (2a) s.t. ∀i ∈ Ω, |θ T d i | ≤ 1. (2b)
A dual point θ is said feasible if it complies with the constraints (2b). The solutions of the primal (1) and dual (2) problems, x and θ respectively, are linked by the relation:

y = Dx + λ θ, ∀i ∈ Ω, | θT d i | < 1 if x(i) = 0 | θT d i | = 1 if x(i) = 0 (3)
We additionally define: d * = arg max d∈{±di} K i=1 d T y, and λ * = d T * y. To avoid the null solution: λ ∈ [0, λ * [.

Screening Tests

The sparsity inducing regularization λ • 1 entails an optimum x that may contain many zeros, and the goal of a screening test is precisely to locate them; an efficient screening test locates many zeros. From the located zeros a screened dictionary D 0 can be defined removing the corresponding atoms, called inactive atoms, from D. Finally the solution of P(λ, D 0 , y) can be readily reconstructed from that of P(λ, D, y). Any optimization procedure using the screened dictionary D 0 therefore computes the solution of P(λ, D, y) at lower computational cost.

Screening tests [START_REF] Ghaoui | Safe Feature Elimination in Sparse Supervised Learning[END_REF][START_REF] Xiang | Fast lasso screening tests based on correlations[END_REF][START_REF] Xiang | Learning sparse representations of high dimensional data on large scale dictionaries[END_REF] are based on an idea emerging from the relation (3) between primal and dual optima, x and θ respectively. According to relation [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]

i veri- fying max θ∈R |θ T d i | < 1.
Screening tests essentially differ from one another on the region R they consider, when R is a sphere [START_REF] Ghaoui | Safe Feature Elimination in Sparse Supervised Learning[END_REF][START_REF] Xiang | Learning sparse representations of high dimensional data on large scale dictionaries[END_REF] max θ∈R |θ T d i | has a closed-form expression and gives the sphere test principle. Spheres that instantiate this principle are described below.

Lemma 1 (Sphere Test Principle [START_REF] Ghaoui | Safe Feature Elimination in Sparse Supervised Learning[END_REF]). If the solution θ of (2)

satisfies ∃{r, c} ∈ R × R N , θ -c 2 ≤ r, then : |c T d i | < 1 -r ⇒ x (i) = 0.
In practice, once a sphere S(c, r) of center c and radius r has been defined, every atom i such that |τ (i)| < 1r are removed, where τ = D T c. The associated screening operator Π c,r (•) is the operator that, given a dictionary D, outputs the corresponding screened dictionary Π c,r (D)

d i s.t. i ∈ Ω, |c T d i | ≥ 1 -r . (4)
The construction of such spheres [START_REF] Ghaoui | Safe Feature Elimination in Sparse Supervised Learning[END_REF][START_REF] Xiang | Learning sparse representations of high dimensional data on large scale dictionaries[END_REF] is based on the following considerations. The dual optimum θ is closer to y/λ than any feasible point in R N . Then from any feasible dual point θ one can construct a sphere, centered on y/λ with radius θ-y/λ 2 , that contains θ. The SAFE method [START_REF] Ghaoui | Safe Feature Elimination in Sparse Supervised Learning[END_REF][START_REF] Xiang | Learning sparse representations of high dimensional data on large scale dictionaries[END_REF], implements the feasible dual point θ = y/λ * .

First-order algorithms for the Lasso

The Lasso problem (1) may be solved with general-purpose first-order algorithms such as ISTA [START_REF] Daubechies | An Iterative Thresholding Algorithm for Linear Inverse Problems with a Sparsity Constraint[END_REF], TwIST [START_REF] Bioucas-Dias | A new twist: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF], FISTA [START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF], SpaRSA [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF], forward-backward splitting [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] or the Chambolle and Pock's primal-dual algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Uzawa | Studies in linear and nonlinear programming, With contributions by Hollis B. Chenery[END_REF]. For the sake of simplicity, ISTA is used as the archetype for first-order algorithms. The extension to all the aforementioned algorithms is described in section 3.

ISTA constructs a sequence of iterates {x k } k≥0 which converges to the optimal x by implementing the update:

x k+1 ← T λ/L k x k - 1 L k D T (Dx k -y) , (5)
where x k is the k-th iterate computed by the procedure, L k is the step size (set according to a backtracking rule see [START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF]), and T t (x) sign(x) max(0, |x|t) is the soft-thresholding operator. In the following update ( 5) is denoted x k+1 ← p k (x k , D) subsequently, in order to ease the notation and emphasize similarities in first-order algorithms. Due to the matrix-vector products involving D and D T , the cost of one update is O(N K), assuming that D has no associated fast transform. In many applications, the dimensions can be large, e.g., K ≥ N ≫ 100. This explains the major interest of reducing the size of the dictionary K without affecting the solution of the Lasso.

Optimizing with Dynamic Screening

Existing screening strategies for the Lasso are static in the sense that they first screen the dictionary and use the screened dictionary to solve the Lasso (see Algorithm 1). We show in this section that calculations made during the optimization procedure can be employed to dynamically and iteratively reduce the dictionary by performing dynamic screening at each iteration. Dynamic construction of better feasible points. Screening tests presented in Section 2.2 build on a feasible dual point. Therefore, producing at each iteration a feasible dual point that is cheap to compute and close to y/λ enable the iterative construction of new SAFE sphere with smaller radius, and thus the iterative construction of more efficient sphere tests.

ISTA directly computes potentially appropriate dual points. Indeed, each update requires the computation of the gradient ∇f

(x) = D T (Dx -y) of the ℓ 2 -fitting term f (x)
Dxy 2 2 . The dual points θ k Dx ky form a sequence {θ k } k≥0 that converges to λ θ. Since θ k is not necessarily feasible, the following dual scaling strategy may be resorted to, in order to give a second sequence of appropriate feasible dual points θk µ k θ k .

Lemma 2 (Dual Scaling [START_REF] Ghaoui | Safe Feature Elimination in Sparse Supervised Learning[END_REF]). Among all feasible scaled versions of θ k , the closest to y/λ is θk = µ k θ k where:

µ k arg min µ∈R µθ k - y λ 2 s.t. D T µθ k ∞ ≤ 1 (6) = min max θ T k y λ θ k 2 2 , 1 D T θ k ∞ , -1 D T θ k ∞ .
Dynamic screening. Embedding this dual scaling strategy within ISTA permits to execute more efficient screening tests almost for free. Indeed, the center of all the spheres is the same, y/λ, hence the test vector D T c is computed only once. Since θ k and D T θ k are computed by the update of ISTA, computing θk requires O(K + N ) operations and computing the sphere radius θky/λ 2 requires O(N ) additional operations. Given that N ≤ K, the computational overhead of the screening test is in O(K). Finally, the total overhead is negligible compared with the O(KN ) operations required for an optimization update.

The resulting ISTA with dynamic screening is presented in Algorithm 3. The input of the algorithm defines: the problem of interest through D, y and λ; the initial state x 0 (set by default to 0); the center c and the function r(•), which compute the radius of the sphere from θ, parameterizing which sphere test is embedded in ISTA-more information on centers and radius functions r(•) are given below in [START_REF] Ghaoui | Safe Feature Elimination in Sparse Supervised Learning[END_REF] and [START_REF] Uzawa | Studies in linear and nonlinear programming, With contributions by Hollis B. Chenery[END_REF]. The algorithm breaks up in two steps: the optimization update (see line 4) where p k (•) returns iterate x k as well as suitable computed vectors θ k and D T θ k ; and the screening step (see line 6-8): at line 6 the feasible dual point θk is computed with the dual scaling strategy, the radius is updated only when decreasing at line 7 and finally the screening operator is applied.

Generalization to other first-order algorithms. As announced in section 2.3, the dynamic screening principle Algorithm 3 ISTA with dynamic screening Require: D, y, λ, x 0 , c, r(.)

1: D 0 ← D,
D k+1 ← Π c,r k+1 (D) ⊲ see (4) 9:
k ← k + 1 10: end while applies to other first-order algorithms as well. Considering that the update does not only modify the iterate x k but some auxiliary variables as well, each first-order algorithm can be describe by its update xk+1 ← p k (x k , D) where xk represents the set of updated variables. Both p k and x are given in Table 1 for a few other representative algorithms.

Table 1 gives two important pieces of information: first all p k (•) have similar computational requirements as ISTA; second the expensive computations required for dynamic screening-computing a new θ k and D T θ k -are provided by these updates. Then dynamic screening applies to these algorithms in the exact same way as it applies to ISTA.

Algorithms

Optimization update xk+1 ← p k (x k , D) TwIST [START_REF] Bioucas-Dias | A new twist: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF] x k+1 ← (1 [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF] same as ISTA except that L k is set with Brazilai-Borwein rule

-α)x k-1 + (α -β)x k xk = {x k , x k-1 } +βT λ x k -D T (Dx k -y) SpaRSA
FISTA [1] x k+1 ← T λ L k z k -1 L k D T (Dz k -y) xk = {x k , z k , t k } t k+1 ← 1 + √ 1 + 4t k 2 z k+1 ← x k+1 + ( t k -1 t k+1 )(x k+1 -x k ) Chambolle- Pock [3] θ k+1 ← 1 1 + σ k (θ k + σ k (Dx k -y)) xk = {x k , xk , x k+1 ← T λτ k x k -τ k D T θ k+1 θ k , τ k , σ k } ϕ k ← 1 √ 1+2γτ k ; τ k+1 ← ϕ k τ k σ k+1 ← σ k ϕ k xk+1 ← x k+1 + ϕ k (x k+1 -x k )
Table 1: Updates for first-order algorithms.

Direct extension to ST3. Section 2.2 presents SAFE spheres, another sphere called ST3 relying on the SAFE sphere have been proposed in [START_REF] Xiang | Learning sparse representations of high dimensional data on large scale dictionaries[END_REF]. Constructed from any feasible dual point θ, ST3 is the sphere centered on c = y/λδd * with radius r(θ) = θy/λ through parameters c and r(•):

SAFE: c = y λ , r(θ) = θ - y λ 2 (7)
ST3: c = y/λ -δd * , r(θ) = θ -y/λ 2 2 -δ 2 (8)
Convergence analysis. First-order algorithms with dynamic screening do not necessarily provide the same iterates as their base version but still converge to the global optimum:

Theorem 3. If a first-order algorithm is proven to converge to the global optimum of the Lasso problem, then its version with dynamic screening converges to the global optimum too.

Proof. As explained in section 2.2, Lasso problems P(λ, D, y) and P(λ, D k , y) for all k ≥ 0 have the same solutions. Since the sequence {r k } k≥0 is non-increasing, Lemma 1 ensures that the set of located inactive atoms is non-decreasing, indeed:

r ≥ r ′ ⇒ (∀i ∈ Ω, |c T d i | < 1 -r ⇒ |c T d i | < 1 -r ′
). This set is upper bounded by the set of zeros in x the solution of P(λ, D, y), so the set of located zeros converges in a finite number of iterations k 0 . Then ∀k ≥ k 0 , D k0 = D k and usual convergence proofs apply.

Numerical Experiments

This section presents experiments used to assess the practical relevance of our approach. The code and data for numerical experiments are released for reproducible research purposes. 1 Runnning Times. We have claimed that, compared with static screening, dynamic screening significantly accelerates the computation of the solution of the Lasso with first-order algorithms. This section evaluates the performance of our method in terms of running times. Note that since each version of the algorithm (no-screening, static screening, dynamic screening) converges to the same optimal x (see Theorem 3), we do not report the value of the objective at convergence .

We measured running times of the algorithm without screening test, with static screening and with dynamic screening. To emphasize the gain, running times are normalized 1 http://pageperso.lif.univ-mrs.fr/~antoine.bonnefoy with respect to running times required by the algorithm without screening. Synthetic data. For experiments on synthetic data, we used two types of dictionaries. The first one is a Gaussian dictionary where observation y as well as all atoms d i are drawn i.i.d. uniformly on the unit sphere by normalizing realizations of N (0, I N ). The second one is the so-called Pnoise introduced in [START_REF] Xiang | Fast lasso screening tests based on correlations[END_REF], for which y and all d i are drawn i.i.d. from the distribution e 1 + 0.1κg and normalized, where g ∼ N (0, I N ), κ ∼ U(0, 1) and e 1 being the first natural basis vector. We set N = 2000 and K = 10000. Audio Data. For experiments on real data we performed the estimation of the sparse representation of audio signals in a redundant Discrete Cosine Transform (DCT) dictionary, which is known to be adapted for audio data. Observations y were taken from 30 music and speech recordings with length N = 1024 and sampling rate 16 kHz. Experiments were run for K = 3N . Experiments. Algorithms were run for several values of λ in order to compute the representation of the observation y in the dictionary D with different sparsity levels. The algorithm stops at iteration k when the ratio between the maximum variation of the functional in (1) and the average of the functional over the M = 10 previous iterations does not exceed the value of ǫ (ǫ = 10 -6 for Gaussian, ǫ = 10 -4 for Pnoise and ǫ = 10 -6 for audio signals).

Figure 1 shows the normalized running times for algorithms with dynamic screening (black squares) and for the corresponding algorithms with static screening (circle) as a function of λ/λ * . Low values account for fast computation. We used Gaussian (left), Pnoise (middle) and Audio (right) data with the ST3 screening test. We plotted the median values over 30 problems. The shaded area contains from 25% to 75% percentiles for plain curves only, in order to illustrate the typical distribution of the values. For all dictionaries, the dynamic screening performs significantly better and is effective in a larger range of λ than the static one. In the audio experiment savings could reach more than 90% over ISTA and up to 70% over ISTA with static ST3 (e.g. for λ/λ * ≈ 0.6). Both static and dynamic screening strategies tend to be more efficient when the dictionary redundancy K/N increases.

Note that due to the normalization of running times, Figure 1 cannot be used to draw any conclusion on which of ISTA, FISTA or SpaRSA is the fastest algorithm. Interpretation: Screening progression. To apprehend the effectiveness of the dynamic screening test, we represented how dynamic screening behaves along the iterations. Figure 2 shows on the same scale the evolution of two key values along the iterations: the radius r k (red dashed line); and the relative size of the dictionary K k /K (blue area)-where K k is the size of the screened dictionary at iteration kwhich represents the proportion of atoms remaining in the screened dictionary. Here dynamic ST3 was used in ISTA for a Gaussian dictionary with λ = 0.7λ * . The reduction of the radius induced a nice improvement in the screening. The screening test may be totally inefficient in the first iterations, which shows the advantage of the dynamic screening strategy over the static one.

Discussion and Future Directions

We have shown that the dynamic screening principle is relevant theoretically and practically. Dynamic screening accelerates more first-order algorithms than static screening in the proposed experiments on synthetic and real data, and in a larger range of λ.

Dynamic screening has been shown to work for several algorithms and screening tests, and the question is whether the concept of dynamic screening can be further generalized. The answer is positive: it can be applied to much more algorithms to other screening As far as an optimization process computes the gradient of the ℓ 2 -fitting term and the screening test rests upon the SAFE sphere, e.g. the dome test [START_REF] Xiang | Fast lasso screening tests based on correlations[END_REF], they can combine into a dynamic screening strategy.

The proposed method raises several questions we plan work on, some of them are addressed here as a conclusion. The SAFE test extends to the Group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], but can it be refined dynamically along the iterations of the optimization process in a similar fashion? As in [START_REF] Ghaoui | Safe Feature Elimination in Sparse Supervised Learning[END_REF], we are curious to see how dynamic screening may show up when other than an ℓ 2 fit-to-data is studied: for example, this situation naturally occurs when classification-based losses are considered. As sparsity is often a desired feature for both efficiency (in the prediction phase) and generalization purposes, being able to work out well-founded results allowing dynamic screening is of the utmost importance.

Fig. 1 :

 1 Fig. 1: Normalized running times on synthetic data (left, middle) and real data (right).

Fig. 2 :

 2 Fig. 2: Relative size K k /K and radius r k along the iterations.

  , atoms d i such that | θT d i | < 1 correspond to inactive atoms. If θ were known it would be easy to identify inactive atoms. It is obviously not the case then to locate inactive atoms screening tests build an upper-bound on | θT d i | by constructing a region R ⊂ R N that contains θ and hence satisfies | θT d i | < max θ∈R |θ T d i |. It allows one to remove every atom d

  r 0 ← +∞ 2: while stopping criteria on x k do 3: ........................ Optimization Update ........................ ............................... Screening .................................

	6: 7:	θk+1 ← µ k θ k r k+1 ← min(r( θk+1 ), r k )	⊲ see (6)

4:

{x k+1 , θ k , D T θ k } ← p k (x k , D k ) ⊲ see

(5) 5:

.8:

  2 2δ 2 , where δ = λ * /λ -1. The corresponding screening operator is Π c,r (•) as defined in[START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. Both SAFE and ST3 can be embedded in Algorithm 3

	0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Normalized running times	0.2 Gaussian dictionary 0.4 0.6 0.8 ISTA + dynamic scr. ISTA + static scr. FISTA + dynamic scr. FISTA + static scr. SPARSA + dynamic scr. SPARSA + static scr.	1.0	0.0 0.0 0.2 0.4 0.6 0.8 1.0 Normalized running times	0.2 Pnoise dictionary 0.4 0.6 0.8 FISTA + dynamic scr. FISTA + static scr.	1.0	0.0 0.0 0.2 0.4 0.6 0.8 1.0 Normalized running times	0.2 ISTA + dynamic scr. 0.4 0.6 Audio Data ISTA + static scr.	0.8	1.0
		λ/λ *			λ/λ *			λ/λ *

This work was supported by Agence Nationale de la Recherche (ANR), project GRETA 12-BS02-004-01.

Centre Inria Rennes-Bretagne Atlantique, R.G. acknowledges funding by the European Re-search Council within the PLEASE project under grant ERC-StG-2011-277906


M Marconi

Member, IEEE J Javaloyes
email: julien.javaloyes@uib.es

Member, IEEE S Balle

Member, IEEE. M Giudici

Passive Mode-Locking and Tilted Waves in Broad-Area Vertical-Cavity Surface-Emitting Lasers

Keywords: Mode-Locking, Broad-Area Lasers, VCSELs

We show experimentally and theoretically that an electrically biased 200 µm multi-transverse mode Vertical-Cavity Surface-Emitting Laser can be passively mode-locked using optical feedback from a distant Resonant Saturable Absorber Mirror. This is achieved when one cavity is placed at the Fourier plane of the other. Such non conventional optical feedback leads to the formation of two tilted plane waves traveling in the external cavity with opposite transverse components and alternating in time at every round-trip. Each of these plane waves gives birth to a train of mode-locked pulses separated by twice the external cavity round-trip, while the two trains are time shifted by a round-trip. A large portion of the transverse section of the device contributes to mode-locked emission leading to pulses of approximately 1 W peak power and 10 ps width. We discuss how inhomogeneities in the transverse section of the saturable absorber select the emitted tilted waves, thus leading to tunable emission over 4 nm.

I. INTRODUCTION

L ASER mode-locking (ML) is a fascinating self-organized cooperative effect involving a large number of longitudinal modes that was recently linked to out-of-equilibrium phase transitions [START_REF] Gordon | Phase transition theory of many-mode ordering and pulse formation in lasers[END_REF]. From a practical point of view many applications require sources of short pulses like e.g. medicine, metrology and communications [START_REF] Haus | Mode-locking of lasers[END_REF]. Passive ML (PML) is arguably one of the most elegant method to obtain such pulses. It is achieved by combining two elements, a laser amplifier providing gain and a saturable absorber acting as a pulse shortening element. Under appropriate conditions, the different dynamical properties of the absorption and of the gain favor pulsed emission by creating a limited time window for amplification around an intensity pulse [START_REF]Theory of mode locking with a fast saturable absorber[END_REF], [START_REF]Theory of mode locking with a slow saturable absorber[END_REF]. PML can also be achieved using artificial absorbers like e.g. nonlinear polarization rotation [START_REF] Dorren | Nonlinear polarization rotation in semiconductor optical amplifiers: theory and application to all-optical flip-flop memories[END_REF], Kerr lens mode-locking [START_REF] Ippen | Principles of passive mode locking[END_REF], Crossed-Polarization [START_REF] Javaloyes | Passive mode locking of lasers by crossed-polarization gain modulation[END_REF] or Stark effect modulation [START_REF] Wilcox | Ultrafast optical stark modelocked semiconductor laser[END_REF]. The PML mechanism has led to the shortest and most intense optical pulses ever generated and pulses in the femto-second range are produced by dye [START_REF] Fork | Femtosecond optical pulses[END_REF] and solid-state lasers [START_REF] Keller | Semiconductor saturable absorber mirrors (SESAM's) for femtosecond to nanosecond pulse generation in solid-state lasers[END_REF]. However, the large size of these devices and the need for optical pumping strongly limit their application domain.

transverse section of the VCSEL contributes to mode-locking leading to pulses of 10 ps width and peak power around 1 W. We analyze the mechanism leading to the selection of the tilted waves and we demonstrate a technique for tuning the central wavelength of the mode-locked emission. Our experimental findings are confirmed by a spatially resolved model for the VCSEL and the RSAM taking into account the multiple reflections in the external cavity and our specific imaging conditions.

II. EXPERIMENT

The VCSEL is a 980 nm device manufactured by ULM Photonics [START_REF] Grabherr | Bottom-emitting vcsel's for high-cw optical output power[END_REF]. Its standalone threshold current (J st ) is about We also found that, for the external cavity length considered,

226

the PML regime is bistable with the off solution for J < J st , 227 see [START_REF] Marconi | How lasing localized structures evolve out of passive mode locking[END_REF] for details. In these conditions and after setting the

IV. THEORETICAL RESULTS

326

Our model for the dynamical evolution of the fields E j and 327 the normalized carrier density N j in the quantum well (QW) 328 regions of the VCSEL (j = 1) and the RSAM (j = 2) is deduced in detail in the appendix, and it reads

330 Ė1 = (1 -iα 1 ) N 1 -1 + i∆ ⊥ + c 1 ∆ 2 ⊥ E 1 + h 1 Y 1 , (1) Ė2 = (1 -iα 2 ) N 2 -z + ib∆ ⊥ + c 2 ∆ 2 ⊥ E 2 + h 2 Y 2 ,(2) Ṅ1 = γ 1 J 1 -1 + |E 1 | 2 N 1 + D 1 ∆ ⊥ N 1 , (3) Ṅ2 = γ 2 J 2 -1 + s|E 2 | 2 N 2 + D 2 ∆ ⊥ N 2 , (4)
where

∆ ⊥ = ∂ 2 x + ∂ 2 y is
O j (r, t) = η j E j (r, t) -Y j (r, t) , (5)
where we defined η j = t j 1 / 1 + r j delayed Algebraic Equations [START_REF] Javaloyes | Multimode dynamics in bidirectional laser cavities by folding space into time delay[END_REF] 358

Y 1 (r, t) = t bs O 2 (r, t -τ e ) , Y 2 (r, t) = t bs O 1 (r, t -τ e ) , (6
) with t bs a complex number whose modulus and phase model

Y 1 (r, t) = t bs O 2 (r ′ , t -τ e )e -i ω 0 c r•r' f d 2 r ′ ≡ F τ (O 2 )(7) Y 2 (r, t) = t bs O 1 (r ′ , t -τ e )e -i ω 0 c r•r' f d 2 r ′ ≡ F τ (O 1 )(8)
and where the beam-splitter losses may contain a normal- returns inverted. This can be seen, for instance, using Eq. ( 5)

375

and saying that the RSAM has an effective reflectivity such 376 that O j = r ef f Y j and substituting in Eq. ( 7), which yields

377 Y 1 (r ⊥ , t) = µE 1 (-r ⊥ , t -2τ ) (9)
where we have used that

F τ • F τ [Y 2 ] = αY 2 (-r ⊥ , t -2τ )
= √ 0.942, 1 384 and r

(2)

1 , r (2) 2 =
391 be L d = L z [q 0 (1 -r 1 r 2 )] -1 = 2 µm.
λ (nm) Y 2 =0.0 Y 2 =0.2 Y 2 =0.3 Y 2 =0.4 Y 2 =0.6
Figure 6. Reflectivity as a function of the input field intensity (left) for various values of the wavelength, the resonance being chosen to be λ 0 = 980 nm and (right) reflectivity as a function of wavelength for increasing input power.

Eqs. (1-4) using a semi-implicit split-step method where the spatial operators are integrated using the Fourier Transform.

Such semi-implicit method is particularly appropriate to the stiffness induced by the broad response of the RSAM. The simulation time was 500 single-trips in the external cavity, i.e. 100 ns.

4) RSAM characteristics:

The basic properties of the RSAM can be found assuming an incoming plane wave of amplitude Y at normal incidence and at frequency ω i.e. Y 2 = Y exp (-iωt). The response of the RSAM is given by

z -iω - J 2 (1 -iα 2 ) 1 + s|E 2 | 2 E 2 = h 2 Y. (10)
The effective reflectivities are obtained setting O 2 = rY 2 that is to say r = E 2 /Y 2 -1 yielding the unsaturated and saturated responses as

r u = h 2 -z + iω + J 2 (1 -iα 2 ) z -iω -J 2 (1 -iα 2 ) , r s = h 2 -z + iω z -iω . ( 11
)
By choosing J 2 = a -h 2 we obtain a perfectly absorbing RSAM, i.e. r u = 0 at the frequency ω = δ + α (a -h 2 ). In between these unsaturated and saturated regimes, the effective reflectivity can be obtained by solving Eq. ( 10) numerically.

We represent the results of this procedure for R = |r| 2 in Fig. 6.

5) Defect in the RSAM surface:

We assume that the impurity on the RSAM surface has the effect of locally lowering the reflectivity of the top mirror. As such the parameters a and h 2 are allowed to vary in the transverse dimension. A decrease of r

1 of 0.1 modifies a and h 2 to be a = 12 and h 2 = 24.

We assume a Gaussian profile for such transverse variations wjhose FWHM is ∼ 4.7 µm.

A. Mode-Locking dynamics

The bias current is fixed to J 1 = 0.92, i.e. below the threshold of the solitary VCSEL but above the threshold of the compound device. We depict in Fig. 7, the output time trace for the intensity averaged over the surface of the VCSEL as well as the averaged population inversion in the active region of the VCSEL and in the RSAM. One recognizes the standard PML temporal pattern where the gain experiences depletion during the passing of the pulse followed by an exponential recovery. The same yet inverted pattern is also visible on the RSAM population inversion. We found the pulse width to be

478

The horizontal axis in Fig. 8c) is scaled in normalized spatial 479 frequencies. The transverse value of the wave vector here is We found that by slowly displacing the center of the pinning 483 potential it was possible to tune the normalized transverse 484 wave vector between q x L = 0 and q x L = 36, which

539

The starting point is the scalar Maxwell equation for a 540 monochromatic field as in [START_REF] Mulet | Mode locking dynamics in electrically-driven 689 vertical-external-cavity surface-emitting lasers[END_REF] 541

∂ 2 z + ∆ ⊥ E (ω, r) + ω 2 υ 2 E (ω, r) = -ω 2 c 2 ε 0 P (ω, r) , (12)
where

∆ ⊥ = ∂ 2 x + ∂ 2
∂ 2 z + ω 2 υ 2 -q 2 ⊥ E (ω, q ⊥ , z) = -ω 2 W c 2 ε 0 P (ω, q ⊥ ) δ l .( 13
)
We assumed that the Quantum Well(s) of width W ≪ λ are 547 located at z = l and defined δ l = δ (z -l). In the empty 548 regions where there is no polarization the solution of Eq. ( 13)

549 reads 550 E (ω, q ⊥ , z) = L + e iQz + L -e -iQz if 0 < z < l R + e iQz + R -e -iQz if l < z < L (14)
where the longitudinal wave vector

Q (q ⊥ , ω) = ω 2 υ 2 -q 2 ⊥ .
551

The boundary conditions at the mirrors and at the QW impose

552 that 553 r 1 L -+ t ′ 1 Y = L + (15) r 2 R + e iQL = R -e -iQL (16) L (l) = R (l) = E (ω, q ⊥ , l) (17)
∂ z R (l) -∂ z L (l) = - ω 2 ε 0 c 2 P (18)
where the primed indexes are for transmission and reflection

F 1 (Q) E = - ω 2 2iQε 0 c 2 ΓW P + F 2 (Q) Y (19)
with 559

F 1 = 1 -r 1 r 2 e 2iQL , F 2 = t ′ 1 e iQl 1 + r 2 e 2iQ(L-l) . ( 20
F 1 = F 0 1 + (ω -ω 0 )∂ ω F 1 + q 2 ⊥ ∂ q 2 ⊥ F 1 + • • • (21)
and transform back to space and time using that ω → ω 0 + i∂ t 572 and q 2 ⊥ → -∆ ⊥ . Since F 1 is at a minimum of its modulus,

573

∂ ω F 1 is purely imaginary, and the field evolution can be

P = ε 0 (α -i)g 0 (N -N t )E ( 23
)
where α is Henry's linewidth enhancement factor, g 0 is the 584 material gain coefficient and N t is the transparency carrier 585 density. The carrier density, in turn, obeys the standard de-586 scription as in [START_REF] Mulet | Mode locking dynamics in electrically-driven 689 vertical-external-cavity surface-emitting lasers[END_REF]. Upon normalization, the evolution for the 587 field and carrier density can be written as

588 ∂E ∂t = (1 -iα) N -1 + i∆ ⊥ + c∆ 2 ⊥ E + hY, (24) ∂N ∂t = γ J -1 + |E| 2 N + D∆ ⊥ N, (25)
where γ is the scaled carrier lifetime, J is the current injection

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Experimental Set-up: Temperature-stabilized VCSEL and RSAM. Coll.: Aspheric Lens, BS : Beam Splitter, M: Mirror and D1/D2: Detectors and CCD cameras. The set up is shown in Fig. 1. Both the VCSEL and RSAM are mounted on temperature controlled substrates which allow for tuning the resonance frequency of each cavity; parameters are set for having the emission of the VCSEL resonant with the RSAM. The light emitted by the VCSEL is collected by a large numerical aperture (0.68) aspheric lens and a similar lens is placed in front of the RSAM. A 10 % reflection beam splitter allows for light extraction from the external cavity and to monitor both the VCSEL and the RSAM outputs. Intensity output is monitored by a 33 GHz oscilloscope coupled with fast 10 GHz detector. Part of the light is sent to two CCD cameras; the first one records the near-field profile of the VCSEL, while the second records the VCSEL's far-field profile. The light reflected by the RSAM is also used for monitoring and a third CCD camera records the light on the RSAM surface. The external cavity length is fixed to L = 30 cm. One of the most

  225

228

  system in the off solution, we were able to start PML emission 229 by perturbing optically the RSAM section at the point where 230 one of the two spots appears when the system operates in 231 the PML regime. The local perturbation has been realized by 232 injecting an external coherent beam tuned with the RSAM 233 cavity resonance and having a waist diameter of less than 234 10 µm. 235 III. DISCUSSION OF THE EXPERIMENTAL EVIDENCES 236 These observations disclose a possible explanation about the 237 emergence of the mode-locking regime in our system. From 238 the dynamical point of view, mode-locking is favored when the 239 saturable absorber is more easily saturated than the amplifier 240 [2]. In our scheme, this is achieved when the VCSEL emits a 241 tilted plane wave, which imposes a low power density (hence 242 low saturation) in the gain section and, at the same time, a 243 strong local saturation of the RSAM. The reason is that the 244 VCSEL and the RSAM lie one into the Fourier plane of the 245 other and in these conditions any plane wave emitted by the 246 VCSEL yields a spot on the RSAM and vice-versa. On the 247 other hand, the unsaturated reflectivity of the RSAM is very 248 low (∼ 1 %) at resonance, rising up to ∼ 60 % when fully 249 saturated. Hence, if the VCSEL emits a pulse in the form of a 250 plane wave with a transverse wave vector K ⊥ , all the power 251 will concentrate on a single spot at r s on the RSAM, hereby 252 strongly saturating the RSAM which becomes reflective. The 253 light thus comes back to the VCSEL, where it arrives with a

Figure 4 .

 4 Figure 4. Off-axis position of a single spot in the far-field profile as a function of the spectral emission peak of the VCSEL. Far-field profile is shown for a discrete number of points in the graph. VCSEL emission is in the mode-locking regime for all green points of the graph while the one at the highest wavelength, where a single on-axis spot appears in the farfield profile, corresponds to an irregular dynamics. VCSEL is biased at 700 mA. The transverse wave vector is selected by laterally shifting the RSAM along the back focal plane of its collimating lens. The large size of the RSAM section ( 4 × 4 mm) with respect to the far-field size (0.7 mm at most) renders this operation feasible. The blue line has been obtained plotting rs = λ 0 f 2π λ -2 -λ -2 0 , with λ 0 = 978.65 nm and f = 8 mm. A similar tuning curve is obtained by laterally shifting the collimating lens of the RSAM.

Figure 5 .

 5 Figure 5. VCSEL time traces emission and corresponding far-field emission profile in the mode-locking regime for two values of the VCSEL bias. The upper and the lower panels correspond to J = 703 mA and J = 285 mA.

297

  On the other hand, wave vector selectivity may also arise 298 from any imperfections on the RSAM/VCSEL mirrors that 299 break the spatial invariance. State-of-the-art fabrication pro-300 cess does not fully prevent from formation of small size (a few 301 µm) defects in the transverse plane of semiconductor micro-302 resonators. These are in general local spatial variations of the 303 semiconductor resonator characteristics [19], which in the case 304 of RSAM would be mainly related to a local variation of un-305 saturated reflectivity. Visual inspection of the reflection profile 306 of a transversally homogeneous injected monochromatic field 307 onto several used RSAMs has indeed shown the existence of 308 these local variations. Such inhomogeneities on the RSAM 309 surface may correspond to locations where the RSAM exhibits 310 a lower unsaturated reflectivity, hence leading the system to 311 select the tilted wave which image on the RSAM surface 312 coincide with one of these defect. By laterally shifting the 313 RSAM along the back focal plane of its collimating lens, 314 the position of the defect changes in the far-field plane, thus 315 selecting a tilted wave with another transverse wave vector 316 component. This leads to a continuous scan of the distance 317 between the two spots in the far-field, as shown in Fig. 4318 and a corresponding change of the emission wavelength, as 319 expected for tilted waves. The experimental evidence of Fig. 4 320 seems to indicate that, in our system, the presence of these 321 inhomogeneities is ruling the selection of the transverse wave 322 vector. The role of small imperfections in the RSAM structure 323 in the transverse wave vector selection is confirmed by the 324 theoretical results in the next section.

  325

1 )

 1 the transverse Laplacian and Y j 331 denotes the field injected in device j. In Eqs. (1-4) time 332 and space are normalized to the photon lifetime κ -1 1 and 333 the diffraction length in the VCSEL L d , respectively. The 334 complex parameter z is decomposed as z = a -iδ where 335 a represents the ratio of the photon decay rates between the 336 VCSEL and the RSAM cavities and δ is the scaled detuning 337 between the two cavity resonances. The scaled carrier recovery 338 rates and the biases are denoted γ j and J j respectively. We 339 define the ratio of the saturation intensities of the VCSEL and 340 the RSAM as s. The angular dependence from the reflectivity 341 of the DBR mirrors as well as the one stemming from the 342 Fabry-Perot cavities are contained in the transverse Laplacian 343 and the parameters c 1,2 . The RSAM being a broadband fast 344 absorber, it is characterized by J 2 < 0, γ 2 ≫ γ 1 , a ≫ 1 as 345 well as s ≫ 1. 346 Injected fields and device coupling: In the case of a 347 self-imaging configuration, the link between the two devices is 348 achieved by expressing the emitted fields O j as a combination 349 of the reflection of the injection fields Y j and the self-emission 350 E j , which reads 351

1 with t j 1 and r j 1

 1 the 352 transmission and reflection coefficients in amplitude of the 353 emitting DBR from the inside to the outside of the cavity. 354 Considering the propagation delay and the losses incurred by 355 the beam-splitter allowing for the light extraction, the link 356 between the two devices in ensured by the following two 357

359

  the losses induced by the beam-splitter and the single trip 360 feedback phase, respectively. 361 However, the image of the VCSEL through its collimator 362 is placed exactly at the object focal plane of the collimator in 363 front of the RSAM. In turn, the RSAM is on the image focal 364 plane of its collimator. In this configuration, and assuming that 365 the collimators have a large diameter, the field injected in one 366 device is the (delayed) Fourier transform in space of the field 367 coming from the other, which is modelized via Kirchhoff's 368 formula for a lens with focal length f as 369

370

  ization constant due to the Kirchhoff's integral. It is worth 371 remarking that Eqs. (7-8) describe the coupling of the two 372 devices to all orders of reflection in the external cavity. They 373 define a delayed map where the field injected into one device 374

378 and µ = t 2 2 )

 22 bs r ef f describes the combined effect of the beam 379 splitter and of the RSAM reflection. Thus, after one round-380 trip, Y 1 overlaps with a spatially reversed copy of itself, and 381 it requires a second round-trip to achieve proper overlap. 382 Parameters: The reflectivities of the DBRs in the VC-383 SEL and the RSAM are taken as r

396 D 1 = 3 )

 39613 We define two 392 numerical domains which are twice the size of the VCSEL 393 and of the RSAM which have both a normalized length 394 L ⊥ = 200. The other parameters are α 1 = 2, α 2 = 0.5, 395 b = 10 -2 ,γ 1 = 10 -3 , γ 2 = 0.1, s = 10, t bs = 0.9, D 2 = 10 -3 c 1 = 10 -2 and c 2 = 10 -6 . 397 Numerical considerations: The simulation of PML 398 lasers is a very demanding problem from the computational 399 point of view: while pulses may form on a relatively short time 400 scale of a few tens of round-trips, the pulse characteristics 401 only settle on a much longer time scale [21]. If anything, 402 the complex transverse dynamics present in our case shall 403 slow down the dynamics even further. Considering a delay 404 of τ e = 2 ns and a time step of δt = 10 -2 implies that 405 one must keep four memory buffers for E j and Y j of size 406 (τ e κ 1 /δt) N x N y ∼ 2×10 3 N x N y with N x and N y the number 407 of mesh points in the two transverse directions. This amounts 408 to 32 Gigabytes with N x = N y = 512. Notice in addition that 409 such values of N x,y are not particularly large if one considers 410 that only half of the mesh discretization in the VCSEL is 411 electrically pumped, the rest being used only for letting the 412 field decay to zero and prevent aliasing, a common problem 413 with spectral methods in broad-area lasers, see for instance 414 [22] for a discussion. The time step of δt = 10 -2 is neither 415 particularly small if one considers the stiffness incurred by the 416 RSAM response. For these reasons we restrict our analysis to 417 a single transverse dimension. We believe that the main spatial 418 feature being lost by such simplification is the description 419 of the non mode-locked two dimensional radial flower mode 420 [17] on the border of the VCSEL, which has already been 421 qualitatively discussed. To conclude, for the sake of simplicity, 422 we used a time delay of τ e = 200 ps (instead of 1 ns) to avoid 423 the temporal complexity of our setup (leading to localized 424 structures, as discussed in [18]), and to concentrate only on 425 the spatial aspects of the problem. We integrated numerically

Figure 7 .

 7 Figure 7. Panels a) and b): Averaged Intensity profile I 1 as a function of time showing stable PML. Panels c) and d) represent the averaged population inversion in the VCSEL and in the RSAM as a function of time.

480∼ 16

 16 which corresponds to 8 transverse oscillations along the 481 active region of the VCSEL, i.e. a wavelength of ∼ 25 µm.

  482

485Figure 8 .Fig. 9 )506Figure 9 .

 899 Figure 8. Panel a): Phase profile of E 1 as a function of space and time. A tilted wave with alternate wave vector from one round-trip to the next is clearly visible. Panel b) Intensity profile of E 1 showing that the transverse intensity profile is uniform. Panel c): Population inversion in the RSAM. The FWHM size of the population inversion spot is ∼ 6 µm and the RSAM is fully saturated locally by the pulse.

516530

  The time traces corresponding to each spot consist of a pulse 517 train with a period 2τ e , and the two trains are delayed by 518 one round-trip one with respect to the other. Accordingly, the 519 two tilted waves are alternatively emitted at every round-trip. 520 We have shown that the mechanism leading to wave vector 521 selection is related to the existence of inhomogeneities in 522 the transverse section of the RSAM. Because the RSAM is 523 in the Fourier plane of the VCSEL near-field, a defect may 524 favor a tilted wave emission with a well-defined transverse 525 component. By shifting the RSAM laterally, the defect moves 526 in the Fourier plane, thus selecting another transverse wave 527 vector and allowing for wavelength tuning of the mode-locked Our approach extends the methodology developed in [23] 531 for single-mode VCSELs under optical injection to the case of 532 broad-area devices. The Quantum-Well active region is consid-533 ered as infinitely thin, and the wave equation is exactly solved 534 inside the cavity for monochromatic plane-waves. Transform-535 ing back to a spatio-temporal representation provides the 536 time-domain evolution of the field inside cavity. Finally, the 537 coupling of the two cavities is included by describing how the 538 injection terms depend on the intra-cavity fields.

  y , P is the polarization of the QW 542 active region and υ = c/n with n the index of refraction. In 543 the longitudinal direction, the cavity is defined by two Bragg 544 mirrors at z = 0 and z = L. 545 Fourier transforming over the transverse coordinates yields 546

554

  processes starting outside of the cavity, r 1 and r 2 the top 555 (emitting) and bottom reflectivities and Y is the amplitude of 556 the external field impinging on the device. After some algebra 557 the equation linking E, P and Y is found to be

F 1 , 2 ω υ l + φ 1 = 2πn 1 , 2 ω υ (L -l) + φ 2 = 2πn 2 where φ 1 , 2 n 2

 1212122 )The modes of the VCSEL correspond to the minima of 560 and the QW is placed in order to maximize the optical 561 confinement factor Γ = 1 + r 1 e 2iQl 1 + r 2 e 2iQ(L-l) for 562 the fundamental mode at q ⊥ = 0. For fixed magnitude of the 563 reflectivities, this is achieved by imposing that 564 are the phases of the reflectivities r 1,2 and n 1 and 565 two integers. Thus the modal frequencies are determined by 566 the effective cavity length L e = L + 2υ (φ 1 + φ 2 ) /ω. Around 567 any modal frequency ω 0 , and in paraxial conditions, F 2 and 568 Γ/Q vary much more slowly than F 1 . In order to have a spatio-569 temporal description of the dynamics of the field, we then fix 570 Γ/Q and F 2 and expand

′ 1 ( 1 +κ ≈ 1 -|r 1 r 2 | 2 | 2υ∂ q 2 ⊥

 111222 P -κτ c E + iL 2 dif f ∆ ⊥ E + t |r 2 |) (-1) m Y + C∆ 2 ⊥ E(22)whereτ c ≈ |r 1 r 2 |2L e /vis the effective cavity transit time, 576 are the cavity losses under normal incidence 577 at resonance, L dif f ≈ |r 1 r Le ω0 is the diffraction length, 578 and C describes the variation of cavity losses with the angle 579 of incidence. 580 The polarization of the QW active region determines the 581 gain and index change induced by the carriers, and we adopt 582 the simple adiabatic approximation 583

589′ 1 Yh 2 = t ′ 3 1 + r 4 1 -r 1 r 2 r 1 r 2 r 3 r 4 , a = 1 -r 3 r 4 1 -r 1 r 2 r 1 r 2 r 3 r 4 , b = r 1 r 2 r 3 r 4

 123121124 above transparency, D is the diffusion coefficient and h 590 describes the coupling of the output field onto the QW and step, we evaluate the field at the laser output O as 593 a combination of the reflection of the injected beam r ′ 1 Y as 594 well as transmission of the left intra-cavity propagating field 595 L -, i.e. O = t 1 L -+ r . Around resonance and using the 596 Stokes relations tt ′ -rr ′ = 1 and r ′ = -r we find, defining 597 η = t 1 / (1 + r 1 ) 598 O = η (-1) m E -Y (27) although the (-1) m is irrelevant in the sense that the QW 599 experiences also a field with a (-1) m which can therefore 600 be removed. In the case of two devices, we scale the photon 601 lifetime, the coupling and the saturation field with respect to 602 the first one leading to the following definitions 603 (28) and s = (g 2 γ 1 ) / (g 1 γ 2 ). 604 ACKNOWLEDGMENT 605 J.J. acknowledges financial support from the Ramon y Cajal 606 fellowship and the CNRS for supporting a visit at the INLN 607 where part of his work was developed. J.J. and S.B. acknowl-608 edge financial support from project RANGER (TEC2012-609 38864-C03-01) and from the Direcció General de Recerca, 610 Desenvolupament Tecnològic i Innovació de la Conselleria

d'Innovació, Interior i Justícia del Govern de les Illes Balears co-funded by the European Union FEDER funds. M.M. and M.G. acknowledge funding of Région PACA with the Projet Volet Général 2011 GEDEPULSE.

waves in the time-averaged near field implies that these waves