O Gaiffe

S Euphrasie

B Cretin

P Vairac

Generation of local magnetic field; application to VCSEL

Keywords: . 41.20.Gz Magnetostatics, magnetic shielding, magnetic induction, boundary-value problems -42.55.Px Semiconductor lasers, laser diodes -42.25.Ja Polarization

In this article, we report on the simulation and realization of a magnetic circuit allowing to apply an axial magnetic field higher than 1 Tesla. The proposed electromagnet is used to investigate the influence of a magnetic field on the polarization behavior of a vertical-cavity-surface-emitting laser (VCSEL). A nonlinear polarization behavior can be experimentally observed for currents as low as the threshold current, provided that the magnetic field is strong enough.

Introduction

In the field of microsystems, the use of electric or magnetic fields for excitation or detection is common, and in most cases the required levels are relatively low. More particularly, for the generation of magnetic fields at this scale, microfabrication techniques are used to make micro-coils [1][2][3][4] with ferromagnetic cores. Yet, when it comes to generate a local magnetic field of about 1 Tesla, this type of techniques reach their limits and it is therefore necessary to develop a specific system in this context. We propose in this article the study and implementation of an electromagnet allowing for the generation of local magnetic fields up to 1 Tesla. We use this specific device as part of a study on the effect of an axial magnetic field on the polarization properties of a vertical-cavity surface-emitting laser (VC-SEL). Several experiments actually show that the light generated by VCSELs exhibits a large variety of polarization behaviors [START_REF] Li | Vertical-Cavity Surface-Emitting Laser Devices[END_REF][START_REF] Jansen Van Doorn | [END_REF][7]. Some papers discuss the effect of an axial magnetic field on the polarization behavior of VC-SELs from a theoretical point of view, but there is only a single report on related experimental results [8]. The motivation of this paper is rooted in the progress concerning the understanding of this effect, a future step being to use these results in metrology, as the very long life time and the high coherence length (around 10 m) of VCSELs are very promising in this area.

Electromagnet

Presentation of the magnetic system

We aim at an electromagnet producing a local magnetic field about 1 Tesla in the beam propagation direction. To a e-mail: pascal.vairac@femto-st.fr localize the magnetic field on the VCSEL, we designed a magnetic circuit based on 2500-turn coil with two conic poles. One of the poles was pierced to allow the VCSEL laser beam to go through. An important criterion to take into account to achieve such an electromagnet is that the laser beam has a divergence of about 20 degrees. Therefore, we have to adapt the dimensions of the hole in order not to block off the beam. The VCSEL was then mounted on a circuit board, along with the two conic poles, the distance between the two cones being set to 2 mm to be able to insert easily the VCSEL. The complete electromagnet is thus contained in a volume of 60×60×20 mm 3 . A schematic of the complete electromagnet design is displayed in Figure 1. The exact dimensions of the electromagnet air gap are detailed in Figure 2. h and s are respectively the cone height and the tip radius. The radius of the base B is set at 10 mm. These three parameters condition the half-angle α of the cone. The second cone is pierced by a cylinder (radius r and length of 1 mm), by a cone (height (h-1 mm), cone tip radius r and fixed base radius R equal to 5 mm) and finally by another cylinder (height H and radius R). This geometry allows us to handle the issue of the VC-SEL divergence. In the next paragraph, we will focus on the influence of the h, s and r parameters on the spatial distribution and intensity of the magnetic field B z in the air-gap. Fig. 2. Zoom of the magnetic circuit. s is the radius of the two cone tips, h is their height and r is the hole radius.

Magnetic simulation

Simulations of this magnetostatic problem were carried out with CST ElectroMagnetic Studio. This commercial software is part of a suite based on the Finite Integration Technique (FIT) and the Perfect Boundary approximation PBA (T M ) [START_REF] Krietenstein | Proceedings of LINAC'98[END_REF][START_REF] Spachmann | [END_REF][START_REF] Lau | Proceedings of Asia-Pacific Conference on Applied Electromagnetics[END_REF]. The FIT discretizes the integral form of Maxwell's equations, rather than the differential one, on a pair of dual interlaced discretization grids. The PBA (T M ) allows us to avoid the disadvantage of the staircase approximation of complex boundaries (for instance cone-shaped boundaries).

In our analysis, the current of the coil is noted I b . We choose a maximum current I bmax of 500 mA to avoid temperature problems in the coil. The experiments reported here were performed at a value of 250 mA. Figure 3 shows the intensity of the field B z along the x axis (at z = 1 mm) for two different shapes of magnetic circuit. The curve with open circles corresponds to a geometry composed of two cylinders (h = 5 mm, s = 10 mm and r = 1 mm), while the curve with solid circles corresponds to a geometry composed of two cones (h = 5 mm, s = 1.5 mm and r = 0.5 mm).

As expected, the magnetic field along the x axis is more localized for the geometry based on the cones than when using cylinders. There is no change in the magnetic field intensity B z at the center of the air gap. This can come from two opposing effects. On the one hand, the magnetic field lines are focused by the geometry of the cone, but on the other hand, there is more magnetic energy that escapes from the loop.

At the center of the air gap, the field value changes with the geometry of the cones as shown in Figure 4 shows the evolution of the field intensity B z on the z axis for different values of s. The other parameters used are h = 5 mm and r = 0.4 mm. The hole drilled to allow the beam output has a significant influence on B z . It result in an asymmetric behavior (along the z axis): the field B z is higher towards z = 0 than towards z = 2 mm. Conversely, when the variable s is between 1 and 10 mm (s = 10 mm corresponding to the case of cylinders), changes are negligible. For s ≥ 1 mm, the magnetic field remains almost constant in the first half of the air-gap (0≥ z ≥ 1 mm). However, Figure 4 shows that the intensity of B z decreases when z becomes larger than 1 mm and is even reduced to approximately 50% of its value at z = 2 mm.

If we want to localise the field by reducing the value of s (s = 0.5 mm), the behavior of the curve is different and the field is no longer constant in the first region of the air gap. At the center of the gap, there is also a decrease of B z by approximately 10%. Using different values of h (from 3 to 7 mm) little change is observable. Figure 5 shows the intensity of the field B z as a function of the hole radius r for s = 0.5 mm. We can see that there is little influence of r on B z for z < 1.5 mm.

For z > 1.5 mm, the field B z decreases with increasing r. As expected, when r decreases, the curve tends to be symmetrical. However, when r = 0 (no cylinder but still a cone-shaped hole) the magnetic field greatly increases near the second cone, showing that the hole concentrates the magnetic field B z . Without the hole the curve is logically symmetrical. In the case of s = 1.5 mm (not represented here), a similar trend can be observed, but the decrease in B z for r ≥ 0.15 mm is not as pronounced. To fabricate our electromagnet, we chose the following parameter values: h = 5 mm, s = 1.5 mm and r = 0.5 mm. As the effect of h is indeed negligible, its value was set to 5 mm in order to relieve the fabrication constraints. Setting s = 1.5 mm allows localizing the field while keeping it constant in the first half of the air gap. r = 0.5 mm has been chosen as it can be easily pierced. The VCSEL was then be placed at the centre of the gap.

With a maximum current I bmax = 500 mA, we get a magnetic field around 800 mT. This value is below the desired one, which is 1 Tesla. For that reason, we also used another design composed this time of two coils. The design is the same with the same parameters for the gap except for the magnetic circuit that is longer, hence allowing to place two coils. This has enabled us to exceed 1 Tesla. However, we observe that the magnetic field begins to saturate, as indicated in the following section.

Comparison between simulations and experiments

The measurement of the magnetic field B z is performed with a gaussmeter. The probe is 1.5 mm thick and the sensor is located in the middle of the probe. Figure 6 which remains under the maximum desired dimension. Figure 7 compares the evolution of B z (with two coils) as a function of current I b for simulation results and experimental measurements. Once again, it shows a good agreement between the simulations and the measured field. 3 Application to VCSEL

Background theory

A VCSEL is composed of two distributed Bragg reflectors (DBR) separated by an active layer. The active region consists of one or several quantum wells allowing for laser light generation. A four-level model has been developed by San Miguel, Feng and Moloney [16] which takes into account the polarization of the laser field by including the spin sublevels of the conduction and valence bands of a semiconductor. Travagnin et al. [19] have extended the model to include the effect of a magnetic field. This model treats the quantum-well gain medium as a four-level medium: right circularly polarized light E + corresponds to the transition from J z = -1/2 to J z = -3/2 and left circularly polarized light E -correspond to the transition from J z = 1/2 to J z = 3/2. The total population between the conduction and the valence band has a decay rate γ p due to spontaneous emission and γ j is the spin-mixing rate and account for the mixing of the population associated with right and left circularly polarized emission.

The global laser dynamics is described by the time evolution of the mean fields E ± and of the carrier densities

D ± : Ė± = -κ[(1 -iα)(1 -D ± )E ± + (ǫ l + ıσ l )E ∓ ±(ǫ c + ıσ c )E ± ] (1) Ḋ± = -γ p [(1 + |E ± | 2 )D ± + γ j γ p (D ± -D ∓ ) -(1 + β)]
where κ is the cavity decay rate; ǫ l and σ l are the linear amplitude and phase anisotropy; ǫ c and σ c are the circular amplitude and phase anisotropy and β is a normalized pumping parameter equal to zero at the lasing threshold.

To study the polarization properties of VCSEL generated light, it is useful to consider separately the amplitudes and phases of the fields, writing

E + = |E + | exp(-ıφ + ) (2)
E -= |E -| exp(-ıφ -) (3)
and to introduce the ellipticity angle

χ = arctan |E + | -|E -| |E + | + |E -| , - π 4 < χ < π 4 (4)

Experimental setup

We used a single-mode VCSEL operating at 850 nm, distributed by Avalon Photonics. The total dimension of the VCSEL chip (electrical contacts included) is 230 µm × 280 µm. The electrical contact has a diameter of 80 µm.

The VCSEL chips are glued with silver paste on a thin printed circuit board serving. The advantage of using these VCSELs lies in their large linewidth equal to 30 MHz which allows their use in metrology.

In our experiment, the injection current I is varied between 1.1 I thr and 1.8 I thr with I thr = 2.8 mA. The resulting optical power hence ranges from 0.3 to 0.6 mW. Under these conditions, the TEM 00 mode is present. At zero magnetic field (B z = 0 T), two peaks separated by about ∆v = 12 GHz are visible in the optical spectrum. This may be caused by a birefringence effect due to the cavity anisotropy of the VCSEL resulting from the manufacturing process [START_REF] Choquette | [END_REF]13]. The higher-frequency peak corresponds to the lasing mode and the lower-frequency peak corresponds to spontaneous emission in the orthogonal polarization mode.

The VCSEL is then put between the two poles of the electromagnet (see Fig. 8). A lens is positioned after the electromagnet to collimate the laser beam. To analyse the ellipticity of the lasing mode, we use a combination of a quarter-wave plate and a polarizer. The polarizer is used to select a privileged axis to determine the behavior of the ellipticity angle as a function of the magnetic field.

Experimental results: influence of the field on the VCSEL

The higher-frequency peak, corresponding to the laser mode, has a polarization with a small ellipticity angle 2χ ≤ 0.06 rad. Figure 9 shows the ellipticity angle for I = 1.1 I thr as a function of the magnetic field B z (I b = 100, 200, 300, 400 and 500 mA). The values of B z used here were obtained from the measurements reported in Figure 6. The experimental curve can be divided into two distinct regions (above and below B z = 500 mT). Both exhibit a linear behavior, though with different slope angles. The overall behavior is hence nonlinear. Under the same experimental conditions, Jansen Van Doorn et al. [8] also observe a linear variation of the ellipticity angle for B z < 500 mT, i.e. in region 1, but with a lower slope, meaning that the influence of the field B z on χ is less significant. The presented results can be explained by the value of the birefringence (12 GHz) which is about twice as large. For stronger magnetic fields (B z > 500 mT), region 2 also shows another linear behavior with a higher slope.

A parameter usually employed in the literature to quantify the influence of a magnetic field on the polarization in VCSELs is the Verdet constant defined by:

V = n c κσ l ∂(2χ) ∂B (5)
Taking n= 3.5 for the average value of the refractive index and κ = 1.10 12 [15,14], we obtain 0.808 rad.cm -1 .T -1 for the first region and 1.479 rad.cm -1 .T -1 for the second. As a comparison, Jansen et al. measured a Verdet constant of 0.97 rad.cm -1 .T -1 for B z < 500 mT, i.e. in region 1. Similar experiments were carried out taking I = 1.8 I thr .

The "nonlinear" effect observed in Figure 9 is accentuated with increasing injection current I. Figure 10 shows the ellipticity angle for I = 1.8 I thr as a function of the magnetic field B z (I b ranging from 0 to 1000 mA with 100 mA steps). B z was obtained from the measurements reported in Figure 7. We observe that the bifurcation point, located between the two regions, moves towards the origin when the injected current increases. The results displayed in Figures 9 and 10 can allow us to retrieve the circular phase anisotropy. From the simplified model presented in [8]:

Ė± = κ ± ıσ c (B)E ± -ıσ l E ∓ (6)
with the eigenvalue tan(2χ) = ±σ c /σ l . It is therefore possible to determine σ c (B). The obtained values of the circular anisotropy can be used to evaluate the influence of the injection current on the ellipticity. σ c can be reinjected in the theoretical model (1), in order to determine the theoretical values of tan(2χ). To do this, we use the following parameters: σ l = πκ∆v = 0.037, γ p = 10 9 s -1 , γ j = 10 11 s -1 [17,18], α = 3, β = 0.1 and ǫ l = ǫ c = 0. For I = 1.1 I thr , the difference between the theoretical and experimental results is about 6%. For I = 1.8 I thr (ie β = 0.8) the difference is greater and is about 30%. These experimental investigations show the limits of the theoretical model in the case of high injection currents and demonstrate that an extended study is necessary to fully account for the observed nonlinear behavior.

Conclusion

We have described in this article the study, the simulation and the development of a specific electromagnet that can generate a local magnetic field higher than 1 Tesla. This magnetic circuit allowed to investigate the polarization behavior of VCSELs in the presence of an axial magnetic field. The experimental results showed that the dependence of this behavior on the magnetic field is non-linear, which can be observed for currents a low as the threshold current if the applied field is higher than 0.5 Tesla. This effect becomes more significant with increasing current. These experimental results are potentially interesting for the scientific community working in this area, since most of the reported works in this field remain limited to theoretical investigations. More generally, we believe that this study of a specific electromagnet for local scale generation of magnetic field can be useful for the characterization or excitation of MEMS.

Fig. 1 .

 1 Fig. 1. Representation of the electromagnet

Fig. 3 .

 3 Fig.3. Intensity of the magnetic field Bz along the x axis (at z = 1 mm) for two different shapes of magnetic circuit. Open circles : two cylinders; solid circles : two cones.

Fig. 4 .

 4 Fig. 4. Intensity of the magnetic field Bz along the z axis for increasing values of s with r = 0.4 mm and h = 5 mm.

Fig. 5 .

 5 Fig. 5. Intensity of the magnetic field Bz along the z axis for increasing values of r with s = 1.5 mm.

Fig. 6 .

 6 Fig. 6. Intensity of the magnetic field Bz at z = 1 mm as a function of the coil current I b for a single coil. Solid curve: simulation results; Solid circles: experimental measurements.

Fig. 7 .

 7 Fig. 7. Intensity of the magnetic field Bz at z = 1 mm as a function of the coil current I b for two coils. Solid curve: simulation results; Solid circles: experimental measurements.

Fig. 8 .

 8 Fig. 8. Experimental set-up

2 Fig. 9 .

 29 Fig. 9. Ellipticity angle 2χ as a function of magnetic field Bz for injection current I = 1.1 I thr .

Fig. 10 .

 10 Fig. 10. Normalized ellipticity angle 2χ as a function of magnetic field B for injection current I = 1.8 I thr .

  The reason is attributable to the saturation of E + (right polarization) and E -(left polarization). As in the previous case, for stronger magnetic fields, i.e. for B z > 500 mT, the behavior remains linear, with a Verdet constant equal to 1.1 rad.cm -1 .T -1 . However, the difference between the slopes at an injection current of I (1.1 I thr and I = 1.8 I thr ) respectively is approximately 25% for the highest values of B z .

		0.35
		0.3
	2χ	0.25
	Ellipticity	0.15 0.2	1			2
		0.1
		0.05
		0	0.2	0.4	0.6	0.8	1
				Magnetic field Bz (T)

The authors wish to express their thanks to S. Benchabane and D. Vernier for many stimulating discussions. This work is supported by La région de Franche-Comté of France.


Carlos A F Fernandes
email: ffernandes@ist.utl.pt

José A P Morgado
email: japmorgado@gmail.com

José B M Boavida
email: jose.boavida@gmail.com

Optimisation of an Asymmetric Three Phase-Shift Distributed Feedback Semiconductor Laser Structure Concerning the Above-Threshold Stability

Keywords: Distributed feedback laser, flatness, matrix-methods, mode selectivity, phase-shift, side-mode suppression ratio, threshold gain

This paper shows that accurately optimised asymmetric three phase-shift (3PS)-distributed feedback (DFB) laser structures can strongly improve the stability of the single-longitudinal mode (SLM) operation, described by the mode selectivity and the flatness of the photon density profile, through an extended range of current injection, when compared to optimised symmetric 3PS-DFB structures reported elsewhere.

This study reveals its importance in modern high bit-rate optical communication systems, by enhancing the possibility of attaining high performance DFB lasers, in easily fabricated structures.

The procedure, based on matrix techniques, aims at the description of the optimal design of the laser structure and it is described step-by-step.

Above-threshold calculations have been accomplished to evaluate the performance of the optimised asymmetric 3PS-DFB structure, namely: the mode selectivity (S), the flatness, the lasing wavelength, the optical power, and the side-mode suppression ratio (SMSR) evolutions with the current injection. For a current injection five times bigger than the threshold current, substantially improvements in S (five times bigger) and in the SMSR (about 9 dB higher) are achieved when compared to similar, but symmetric, DFB structures.

Optimisation of an Asymmetric Three Phase-Shift Distributed Feedback Semiconductor Laser Structure Concerning the Above-Threshold Stability I. INTRODUCTION Nowadays, distributed-feedback (DFB) lasers play a decisive role in high bit-rate optical communication systems using direct modulation [START_REF] Bornholdt | 40 Gbit/s directly modulated passive feedback DFB laser for transmission over 320 km single mode fibre[END_REF]- [START_REF] Wedding | 43 Gbit/s transmission over 210 km SMF with a directly modulated laser diode[END_REF], where they should present single longitudinal mode (SLM) operation over the largest range of current injection. Therefore, laser structures with a very large normalised mode selectivity (S), combined with sufficiently flat intra-cavity fields -usually quantified by a flatness parameter (F) -are demanded at threshold, so that modal instabilities induced above threshold by spatial hole burning (SHB) can be avoided [5, p. 159], [6, p. 249], [7, p. 123], [START_REF] Fessant | Gaussian-like tapered grating quarter wave-shifted DFB semiconductor lasers for high-power single-mode operation[END_REF]. These constraints put increasing demands on the design of directly modulated lasers with adequate laser performance, in easily fabricated structures.

To avoid the degradation of SLM operation with the current injection, various non-conventional DFB laser structures have been proposed. Some of them point to the use of gain coupling structures with high selectivity and reduced SHB [5, p. 173], [9, p. 185]. However, they are generally associated with technological processes that are hard to implement [5, p. 173], [START_REF] Kapon | The effect of complex coupling coefficient on distributed feedback lasers[END_REF], [START_REF] David | Basic analysis of AR-coated, partly gain-coupled DFB lasers: the standing wave effect[END_REF]. Therefore, index-coupled structures have been used more often, namely: with multiple phase-shifts (MPS-DFB) placed along the cavity [START_REF] Agrawal | Distributed feedback lasers with multiple-phase-shift regions[END_REF], [START_REF] Ogita | Long cavity multi-phase-shift Distributed Feedback Laser for linewidth narrowing[END_REF], [START_REF] Ogita | Optimum design for multiple phase shift Distributed Feedback Laser[END_REF], [7, pp. 124-135]; with distributed coupling coefficients (DCC-DFB) in standard quarter wavelength-shifted (QWS) structures [7, p. 136-141]; and with chirped gratings in corrugation-pitch-modulated (CPM-DFB) structures [START_REF] Fessant | Threshold and above-threshold analysis of corrugationpitch-modulated DFB lasers with inhomogeneous coupling coefficient[END_REF]. Moreover, improvements can be achieved in structures with adequate combinations of distributed phase-shift (PS) -MPS or CPM structures -and distributed coupling strength -DCC structures -along the cavity [7, p. 141-145], [5, p. 169], [START_REF]Gaussian-like tapered grating quarter wave-shifted DFB semiconductor lasers for high-power single-mode operation[END_REF]. However, multi-section structures associated with axial variations in the pitch or/and in the coupling coefficient induce additional technological problems, making more difficult their monolithic integration [5, p. 169]. For this reason, MPS-DFB structures have been considered up to now as the most feasible solution [START_REF] Morthier | A new DFB laser diode with reduced spatial hole burning[END_REF]- [START_REF] Agrawal | Modelling of distributed feedback semiconductor lasers with axially varying parameters[END_REF]. Besides, recent studies have reported that the design of the PSs in those kind of lasers can be accurately set [START_REF] Bastard | Realization of a distributed phase-shifted glass DFB laser[END_REF]. Unfortunately, the flatness improvement of the symmetric MPS-DFB structures is achieved at the expense of a drastic reduction in the mode selectivity [7, p. 145]. On the other hand, and as far as the authors know, the asymmetric MPS-DFB structures have not been, up to now, subjected to an accurate optimisation [START_REF] Ogita | Long cavity multi-phase-shift Distributed Feedback Laser for linewidth narrowing[END_REF], [START_REF] Ogita | Optimum design for multiple phase shift Distributed Feedback Laser[END_REF].

This paper aims at to show that optimised asymmetric 3PS-DFB structures may strongly improve the laser performance above threshold, by reducing the modal instabilities induced by SHB, when compared with currently reported symmetric DFB structures [7, p. 145].

In section II, the transfer matrix method (TMM) is reviewed, both for threshold and above-threshold regimes. Special attention is given to the details concerning the adequate strategies searching for efficient convergence procedures above threshold. The optimisation of the 3PS-DFB structure is carried out in section III at threshold regime, aiming at the calculation of the structural parameters concerning the optimised laser structure. Still in this section, the most usual laser figures of merit, like the mode selectivity, the flatness and the side mode suppression ratio (SMSR) are investigated for the optimised 3PS-DFB laser structure in the above threshold regime. A comparative analysis with other structures referred elsewhere is also performed. Finally, the main conclusions are summarised in section IV.

II. THE TRANSFER-MATRIX-METHOD

A. Threshold Analysis

To perform the TMM-based model for laser threshold analysis, the laser cavity with length L is divided into M concatenated sections, which are identified by the constancy of its structural parameters. For the m-th section with length L m , those structural parameters are: the corrugation period Λ m , the amount of feedback per unit length κ m and the phase of the section grating with respect to the left side of the cavity Ω m .

Each section is described by two counter propagating electrical field waves described by their complex amplitudes E R (z) and E S (z), which allow the internal electrical field intensity E(t, z), to be determined according to

E(t, z) ∝ E R (z) + E S (z) • exp (j ω t) , ( 1
)
where t is the time, z is the z-axis coordinate, { • } is the real part operator, j = √ -1 and ω is the field angular frequency. Equation (1) assumes that the longitudinal laser axis coincides with the z-axis. In the TMM the column matrices related to E R (z) and E S (z) components are considered for the same spatial position (see Fig. 1).

On the basis of the coupled wave equations, the transfer matrix for the m-th section of the one-dimensional DFB laser structure indicated in Fig. 1 is given by [7, p.105], [START_REF] Fessant | Threshold and above-threshold analysis of corrugationpitch-modulated DFB lasers with inhomogeneous coupling coefficient[END_REF] and it links the column matrices related to the complex electric fields of the wave solutions at z m and z m+1

T(z m+1 /z m ) = ⎡ ⎣ t m) 11 t m) 12 t m) 21 t m) 22 ⎤ ⎦ , ( 2
) E R (zm) E R (z m+1 ) E S (zm) E S (z m+1 ) z z m+1 zm Λm κm Ωm
E R (z m+1 ) E S (z m+1 ) = T(z m+1 /z m ) • E R (z m ) E S (z m ) , ( 3
)
where t m)

11 , t m) 12 , t m)
21 and t m)

22 are given, respectively, by

t m) 11 = ξ m -ρ 2 ξ -1 m 1 -ρ 2 ζ m ; t m) 12 = - ρ ξ m -ξ -1 m e -j Ωm 1 -ρ 2 ζ m t m) 21 = ρ ξ m -ξ -1 m e j Ωm 1 -ρ 2 ζ -1 m ; t m) 22 = - ρ 2 ξ m -ξ -1 m 1 -ρ 2 ζ -1 m , (4)
with ξ m = e γ m • (zm+1-zm) and ζ m = e j βm • (zm+1-zm) .

The propagation constant, β m , and the complex propagation constant, γ m , are given, respectively, by

β m = π Λ m ; γ m = (α -j δ m ) 2 + κ 2 m , (5)
where ρ m and δ m are given by

ρ m = j κ m α -j δ m + γ m ; δ m = δ + π 1 Λ 1 - 1 Λ m , (6)
with α and δ being, respectively, the gain and detuning for the propagation modes taking the left section as a reference. Ω m is given by [START_REF] Fessant | Threshold and above-threshold analysis of corrugationpitch-modulated DFB lasers with inhomogeneous coupling coefficient[END_REF]]

Ω m = Ω 1 + 2 • m-1 k=1 π Λ k • L k ; 2 ≤ m ≤ M . ( 7
)
The fields at both ends of the cavity are connected by the elementary matrix product

E R (L) E S (L) = T Total • E R (0) E S (0) , ( 8
)
where

T Total = 1 m=M T(z m+1 /z m ) . ( 9
)
The formulation of the transfer matrices for other modified DFB laser structures is straightforward, as far as the changes are correctly translated to the matricial formalism. Namely, the inclusion of the PS discontinuity ϕ, is given by the following matrix [START_REF] Ghafouri-Shiraz | Distributed Feedback Laser Diodes and Optical Tunable Filters[END_REF] e j ϕ 0 0 e -j ϕ , [START_REF] Kapon | The effect of complex coupling coefficient on distributed feedback lasers[END_REF] assuming that the field discontinuity is usually small along the plane of the PS. In this case the matrix (10) should be included in the matrix product T Total indicated in the equation ( 9), at the correspondent z position. The oscillation condition corresponds to the vanishing of the incoming waves and it is determined by the following requirement

t Total 22 (α, δ) = 0 , ( 11
)
where t Total

22

is the fourth element of the matrix T Total , given by [START_REF]Fiber-Optic Communication Systems[END_REF]. The solutions are the mode gain, α, and the detuning, δ, and are related to the modes that are allowed to propagate inside the cavity. For the main mode their values are, respectively, the threshold gain, α th , and the threshold detuning, δ th . For a grating with a first-order Bragg diffraction, the mode gain and the detuning can be expressed, respectively, as [7, p. 151]

α(z) = Γ g(z) -α loss 2 (12)
and

δ(z) = 2 π λ n(z) - 2 π n g λ λ B • (λ -λ B ) - π Λ(z) , ( 13
)
where Γ is the optical confinement factor, α loss is the total loss, n is the effective index, λ B is the Bragg wavelength, λ is the lasing mode wavelength, n g is the group effective index and g is the material gain, given by [7, p. 151]

g(z) = A 0 • (N (z) -N 0 ) - -A 1 • λ -λ 0 -A 2 N (z) -N 0 2 . ( 14
)
In [START_REF] Ogita | Optimum design for multiple phase shift Distributed Feedback Laser[END_REF], N is the carrier concentration, A 0 is the differential gain, N 0 is the carrier concentration at transparency (g = 0), λ 0 is the peak wavelength at transparency and A 1 and A 2 are parameters used in the parabolic model assumed for the material gain. Using the first-order approximation for the effective index n, one obtains [7, p. 151]

n(z) = n 0 + Γ ∂n ∂N N (z) , ( 15
)
where n 0 is the effective index at zero carrier injection and ∂n/∂N is the differential index. The photon concentration (S) and N are coupled together through the steady-state carrier rate equation [7, p. 152]

I q V act = A N(z) + B N 2 (z) + C N 3 (z)+ + v g g(z) S(z) 1 + ε S(z) , ( 16
)
where I is the injection current, q is the modulus of the electron charge, V act is the volume of the active layer, A is the spontaneous emission rate, B is the radiative spontaneous emission coefficient, C is the Auger recombination coefficient, ε is a non-linear coefficient to take into account saturation effects and v g = c/n g is the group velocity, with c being the free space velocity.

In an purely index-coupled DFB laser cavity, which happens to be the case considered along this paper, the mutual interaction between the coupled waves E R (z) and E S (z) can be neglected in the rate of total power change [7, p. 59], [START_REF] Kapon | The effect of complex coupling coefficient on distributed feedback lasers[END_REF]. Therefore, the local photon density inside the cavity can be expressed as [7, p. 152]

S(z) ≈ 2 ε 0 n(z) n g λ h c • c 2 0 E R (z) 2 + E S (z) 2 , ( 17
)
where ε 0 is the free space permittivity, h is the Planck's constant and c 0 a dimensionless coefficient that allows the determination of the total electric field at the above threshold regime, taking into account that the normalization

E R (0) 2 + E S (0) 2 = 1 (18)
has been imposed. The equation ( 18) and the boundary conditions imposed at the left facet allow the calculation of the two counter running waves, E R (z) and E S (z), at z = 0. The use of the TMM allows the calculation of the longitudinal electric field profile. The output power at the right facet can be determined as [7, p. 152]

P = d w Γ • v g • h c λ • S(L) , (19)
where d and w are the thickness and width of the active layer, respectively.

From the solutions of the oscillation condition [START_REF] David | Basic analysis of AR-coated, partly gain-coupled DFB lasers: the standing wave effect[END_REF], α th and δ th are determined. Using equations ( 12)-( 15), the carrier concentration at threshold (N th ), the effective index at threshold (n th ), the threshold wavelength (λ th ) and λ 0 are successively evaluated. Threshold current (I th ) is obtained from [START_REF]Gaussian-like tapered grating quarter wave-shifted DFB semiconductor lasers for high-power single-mode operation[END_REF], assuming that S is negligible at threshold. Within this assumption, the z dependence described in eqs. ( 12), ( 14) and ( 15) is also neglected. The same assumption is valid for eq. ( 13), except for CPM structures where a z dependence is included in Λ(z).

B. Above-Threshold Analysis

In the above-threshold regime, S(z) is high enough to induce important non-uniformities in N (z) and n(z). Despite the SHB effect can be minimised by an adequate design of the DFB structure, the interdependence of S(z), N (z) and n(z) induces strong longitudinal inhomogeneities that forces the division of each section into several sub-sections, in order to ensure a correct evaluation of the above-threshold characteristics. According to Ref. [7, p. 153], for a cavity with L = 500 µm, at least a total of M = 5000 subsections are needed.

In this paper, the numerical procedure for the abovethreshold calculations follows closely the method developed in Refs. [7, p. 149], [START_REF] Fessant | Threshold and above-threshold analysis of corrugationpitch-modulated DFB lasers with inhomogeneous coupling coefficient[END_REF]. However, in order to ensure a quick convergence in the evaluations of the laser characteristics, an adequate strategy is now proposed, which is fully described below.

1) Lasing mode analysis: For each bias current I, the numerical above-threshold calculations concerning the lasing mode are summarised as follows:

(i) Successive (G × G) grids are created in the (c 0 , λ) plane. The i-th grid is centred at c i) 0c , λ i) c
and it is enclosed in the region defined by the limits c i)

0min , c i) 0max , λ i) min and λ i) max . For the initial grid (i = 1) 1 λ 1) c = λ th (20) c 1) 0c = h c (I -I th ) 2 q V act v g g th ε 0 n th n g λ th E R (0) 2 + E S (0) 2 (21) c 1) 0min = c 1) 0c -Δc 1) 0 ; c 1) 0max = c 1) 0c + Δc 1) 0 (22) λ 1) min = λ 1) c -Δλ 1) ; λ 1) max = λ 1) c + Δλ 1) , ( 23
)
where Δc

1) 0 = c 1)
0c /10 and Δλ 1) = 0.1 nm seem adequate for the majority of the DFB structures, when G ≈ 10. However, a readjustment of Δc 1) 0 and Δλ 1) may, occasionally, be necessary in order to prevent the mode hopping. This is a critical aspect of the proposed analysis, since an inadequate choice would mislead the results, leading to solutions with no physical meaning;

(ii) For each one of the (G × G) pairs of the i-th grid c i)

0 k , λ i) l
with k ; l = 1 . . . G, the equations ( 14) -( 17) are self-consistently solved, in order to determine the material gain, carrier density, photon density and effective index for each one of the j sub-section, respectively, g j , N j , S j and n j , with 1 ≤ j ≤ M ;

(iii) The equations ( 12) and ( 13) are solved in order to determine the lasing mode gain and detuning for the j-th sub-section, respectively, α j and δ j . The transfer matrix of the j-th sub-section, T (z j+1 /z j ), is than calculated;

(iv) Using the TMM, the two counter-running waves at the output of the j-th sub-section, E R j and E S j , are obtained. For the M -th sub-section, the discrepancy found between those values and the laser right facet boundary condition is represented by ε i) kl . This value is evaluated and stored for each pair c i)

0 k , λ i) l
of the i-th grid. The error associated to the i-th grid is given by ε

i) = min ε i) kl ; (v) Whenever ε i) = ε i-1) , the central pair remains the same c i+1) 0c = c i) 0c , λ i+1) c = λ i)
c , but new limits are required for the next grid. c 0 and λ discretizations should be reduced, for instance: 1) , the pair associated with the min ε i) kl is chosen as next central 1 Notice that, according to (18), c

Δc i+1) 0 = Δc i) 0 /10 and Δλ i+1) = Δλ i) /10. Whenever ε i) < ε i-
1) 0c is numerically equal to h c (I -I th ) / (2 q Vact vg g th ε 0 n th ng λ th ). pair c i+1) 0c , λ i+1) c
, while c 0 and λ discretisations remains unchangeable. For i = 1, ε i-1) is taken as the error associated with the central pair c

1) 0c , λ 1) c .
For each one of the (G × G) pairs c i+1)

0 k , λ i+1 l
, the steps (i-v) are repeated until ε i+1) ≤ ε min , where ε min is a preset error value (less than 10 -14 , as indicated in Ref. [7, p. 156]).

Since the gain α j and the detuning δ j are z-dependent, the lasing characteristics for each biasing current are associated with their mean values along the cavity, given by

α(I) = 1 M M j=1 α j (I) (24)
δ(I) = 1 M M j=1 δ j (I) . ( 25
)
Notice that the sequential analysis (i)-(v) assumes a onemode propagation laser behaviour. This approach is itself a good assumption, since the present analysis focuses on DFB structures that must ensure SLM operation. Otherwise, different strategies must be adopted.

When studying the I influence on the laser characteristics, a considerable CPU time reduction can be achieved if, instead of using ( 20), λ 1) c = λ is assumed, where λ is the solution found in a previous bias current.

2) Side mode analysis: Now, the procedure starts by assuming the threshold solution (α, δ), associated with the side mode under analysis. Steps (i)-(v) are carried out, with the exception of the calculations concerning n(z). Since the effective index profile is settled by the main-mode distribution, n(z) is the distribution obtained in subsection B.1. Therefore, while the steady-state carrier rate equation ( 16) is still valid in the side mode analysis, the equation ( 15) is no longer used. This means that the calculated S(z) and N (z) distributions have no physical meaning in this context, since the main mode has been suppressed in this analysis. Only the solutions δ(I), α(I) and λ(I) shall be kept.

To summarise, a method of determining the laser main mode characteristics such as N (z), S(z), n(z), λ, P , α and δ versus the biasing current has been described. Besides, the assessment of the side modes λ, α and δ is also possible. Therefore, crucial information about the mode stability, the mode selectivity and the SMSR with the current injection can be made available for the DFB structure under analysis.

III. THE ASYMMETRIC 3PS-DFB LASER

A. Structure Description and Parameters Definition

The laser structure under analysis is represented in Fig. 2. This is a multi-section anti-reflection (AR)-coated with uniform grating period Λ m = Λ and uniform coupling coefficient κ m = κ . The PS discontinuities, ϕ i , are located along the DFB laser structure. Their positions are represented by a normalised parameter given by where z i is the ϕ i position. A purely index-coupled laser structure, which assures that κ is real [7, pp. 56-57], is assumed, with L = 500 µm, Λ = 227.039 nm, Ω 1 = 0 rad and PSP 2 = 0.5.

PSP i = z i L i = 1 ; 2 ; 3 , ( 26
)
The normalised mode selectivity and the flatness of the electric field are given, respectively, by [7, p. where I(z) is the normalised electric field intensity at an arbitrary position, z, given by2

I(z) = E R (z) 2 + E S (z) 2 E R (0) 2 + E S (0) 2 , ( 28
)
and I is its average value along the cavity. Besides, (α th L) is the normalised threshold gain and (α L) is the normalised gain of the main side mode. The laser structural and material parameters used along the paper are summarised in Table I Ref. [7, p. 157].

B. Structure Optimisation

The objective of the 3PS-DFB laser structure optimisation is both to maximise S and to minimise F at threshold, by varying simultaneously and independently the following set of variables (decision variables): (κ L), ϕ 2 , PSP 1 , ϕ 1 , PSP 3 and ϕ 3 .

Initially S ≥ 0.25 S min = 0.25 and F ≤ 0.05 F max = 0.050 is required, since these limits have been generally accepted as good boundaries in order to reach a

z L 0 PSP 1 PSP 2 = 0.5 PSP 3 1 ϕ 1 ϕ 2 ϕ 3 Λ AR-coated AR-coated
Fig. 2 -A Step 1 and Step 2 are sequentially repeated (Steps 1-2) until no improvements on S and F are achieved. The optima values for the structural parameters PSP 1 , ϕ 1 , PSP 3 , and ϕ 3 are found, assuming κ L = 2 and ϕ 2 = 90 • . New optima boundaries (S min , F max ) for the next step are settled.

Step 3 An optimisation of (κ L) and ϕ 2 is performed. Therefore, assuming the laser structure defined at the end of Steps 1-2, the optimisation of both functions S(κ L, ϕ 2 ) and F(κ L, ϕ 2 ) is performed by varying simultaneously and independently κ L and ϕ 2 in their entire ranges:

1 ≤ κ L ≤ 3 and 0 • ≤ ϕ 2 ≤ 180 • . The following selection criterium is used: S κ L, ϕ 2 ≥ S min and F κ L, ϕ 2 ≤ F max ;
Steps 1-2 and Step 3 are repeated until no improvements on S and F are achieved. This means that optima values for S and F, at which correspond the best 3PS-DFB laser structure (PSP 1 , ϕ 1 , ϕ 2 , PSP 3 , ϕ 3 and (κ L)), have been obtained. In all steps, and whenever necessary, an argument based on the smallest threshold gain is used in order to decide optima solutions. Fig. 3 shows the contour maps for S(PSP 3 , ϕ 3 ) and F(PSP 3 , ϕ 3 ), when PSP 3 and ϕ 3 vary along their entire range, Table II summarises the results for S, F and α th L achieved for three different laser structures: the optimised 3PS-DFB, the QWS-DFB and the symmetric 3PS-DFB of Ref. [7, pp. 128-135] indicated as a possible optimum design for 3PS-DFB laser structures. As far as F is concerned, the optimised 3PS-DFB laser structures are clearly advantageous. This is not surprising, since the inclusion of several PS along the laser cavity flattens the field distribution. However, it is worth noticing that the optimised asymmetric 3PS-DFB structure reaches higher S than those obtained for both the QWS-DFB and the symmetric 3PS-DFB structures.

A threshold analysis for the 3PS-DFB structure has been presented. Nevertheless, even if a structure presents an adequate performance at threshold, an above-threshold analysis is essential in order to assess the rate at which the SHB deteriorates the laser features with the increasing current.

C. Above-Threshold Analysis of the Optimised Structure

The S(z) distribution in the asymmetric 3PS-DFB laser cavity can be observed in Fig. 4 QWS-DFB lasers [7, p. 159]. In the symmetric case of Ref. [7, p. 163] these local maxima coincide with the PS positions. In the optimised asymmetric 3PS-DFB laser the local maxima are displaced to the cavity facet vicinity. When compared to the QWS-DFB laser, both 3PS-DFB lasers show smaller differences between the central photon density and the escaping photon densities at the facets. This difference is especially small in the optimised asymmetric structure. This fact benefits the laser performance as far as the emitted power is concerned, as it shall be seen below.

A comparative analysis between the QWS and the 3PS-DFB lasers is available in Figs. 567.

Fig. 5(a) shows the laser mode selectivity versus current injection. Similar mode discriminations at threshold are shown for the QWS-DFB laser and for the asymmetric 3PS-DFB laser. Nevertheless, the mode selectivity has a severe reduction with increasing biasing currents in the QWS case, showing that this laser is strongly affected by the SHB. On the contrary, the asymmetric 3PS-DFB laser is very stable in the high power regime. For the symmetric 3PS-DFB laser, it is apparent that the stability related to flatter photon profiles was obtained at the expense of a great reduction in the mode selectivity at threshold (about half the value of the QWS). Concerning the mode gain discrimination, the asymmetric 3PS-DFB structure is, undoubtedly, the best option. Once more, it shall be emphasised that an adequate profile design for the 3PS-DFB structure revealed crucial for the improvement of the laser performance in the high power regime. Fig. 5(b) focuses on the evolution of the flatness with current injection for the two 3PS-DFB lasers, showing a monotonically decreasing function for both structures. It shall be emphasised that in both cases the flatness lies within the selection criterium commonly referred for lasers with 500 µm length long: F < 0.05 [7, p. 131]. On the contrary, the QWS-DFB laser presents high non-uniformities in the photon profile (F = 0.30 at threshold, and F = 0.079 for I = 5 I th ).

Fig. 6(a) summarises the results obtained for the lasing wavelength versus normalised current injection for the three lasers. The QWS laser shows the larger dynamic change in the lasing wavelength. On the other hand, both 3PS-DFB lasers are associated with more stable outputs, foreseeing their potential use as optical carriers. In the current range 1 ≤ I/I th ≤ 5, relative variations (Δλ/λ th ) in the emitted wavelengths of 9.4 × 10 -4 %, 1.5 × 10 -3 % and 5.5 × 10 -3 % are observed for the asymmetric, the symmetric and the QWS lasers, respectively. Fig. 6(b) compares the light-current stationary characteristic for the three lasers. Under similar normalised current injections, the asymmetric 3PS-DFB laser shows larger values for the optical power output (measured at the right facet), which establishes higher external differential efficiencies for this laser. As previously referred, this is due to the smaller difference between the central photon density and the escaping photon densities at the facets, for the asymmetric laser structure.

The measurement of the laser spectral characteristics is a way of checking its single mode stability. Fig. 7(a) shows the normalised spontaneous emission power for I = 1.5I th and I = 5I th , for the asymmetric 3PS-DFB laser. High values of SMSR are achieved, for both bias currents. Besides, it is worth noticing that the "blue-shift" in wavelengths is negligible.

The inset of Fig. 7(a) shows the δ-α plot for the modes that are allowed to propagate inside the cavity at threshold. This figure points out two positive side-modes that are very close in frequency (encircled by a dashed line), which originates the broadening of the spectrum around 1546.4 nm. Another relevant aspect lies with the fact that, near 1547.7 nm, the spectral amplitude of the dominant mode remains at a high value with current injection, showing no severe mode competition in the high power regime. This is pin-pointed in Fig. 7(b), where the SMSR of the asymmetric DFB laser is maintained throughout the range of biasing current under analysis. Besides, with the increasing current, the other two lasers show a strong degradation of SMSR, which becomes lower than the required 30 dB for the SLM operation [5, p. 147], over the most part of the current range.

IV. CONCLUSIONS

An designed to single-mode operation has been analysed in the above-threshold regime. A correct optimisation of the PS discontinuities profile at threshold ensures the best mode gain margin with flattened intracavity intensity profile. Based on TMM and on the carrier rate equation, above-threshold numerical simulations have shown that the optimised 3PS-DFB structure achieves very good SLM operation and maintains these characteristics up to high biasing currents. Comparative studies with a similar optimised symmetric 3PS-DFB laser, referred elsewhere in the literature, have shown that the laser proposed in this paper, presents remarkable increases in the mode selectivity (about five times for I/I th = 5), in the SMSR (about 9 dB for I/I th = 5) and in the lasing wavelength stability (about two times in the overall current range).

Fig. 1 -

 1 Fig. 1 -A simplified schematic diagram for a one-dimensional DFB laser structure section, placed between z = zm and z = zm+1.

  131 and p. 128] S = (α L) -(α th L) ;

Fig. 3 -

 3 Fig. 3 -(a) S(PSP3, ϕ 3 ) and (b) F(PSP3, ϕ 3 ), assuming κ L = 1.7, PSP2 = 0.5, ϕ 2 = 60 • , PSP1 = 0.127 and ϕ 1 = 110 • 7 • . Values for S ≥ Smin or F ≤ Fmax are represented by a solid line. Smin = 0.64 and Fmax = 0.012 are considered. (× solution: PSP3 = 0.64 and ϕ 3 = 100 • , which corresponds to S = 0.78 and F = 0.010)

I = 1 .5 th 4 -Fig. 5 -

 145 Fig. 5 -(a) Selectivity and (b) Flatness vs current injection for the optimised asymmetric 3PS-DFB, the symmetric 3PS-DFB of Ref. [7, p. 128] and the QWS-DFB laser structures, with L = 500 µm. In (b), the QWS-DFB laser flatness falls outside the axis limits.

Fig. 6 -

 6 Fig.[START_REF] Agrawal | Semiconductor Lasers[END_REF] -power vs current injection for the optimised asymmetric 3PS-DFB, the symmetric 3PS-DFB of Ref.[7, p. 128] and the QWS-DFB laser structures, with L = 500 µm.

Fig. 7 -

 7 Fig. 7 -(a) Above-threshold normalised spontaneous emission spectra under two different biasing currents and (b) Side-mode suppression ration vs current injection for the optimised asymmetric 3PS-DFB, the symmetric 3PS-DFB of Ref. [7, p. 128] and the QWS-DFB laser structures, with L = 500 µm.

TABLE I -

 I Summary of laser parameters

	Laser parameter	Value
	Material parameters
	Spontaneous emission rate, A	2.5 × 10 8 s -1
	Bimolecular recombination coefficient, B	1.0 × 10 -16 m 3 • s -1
	Auger recombination coefficient, C	3.0 × 10 -41 m 6 • s -1
	Differential gain, A 0	2.70 × 10 -20 m 2
	Gain curvature, A 1	1.50 × 10 19 m -3
	Differential peak wavelength, A 2	2.70 × 10 -32 m 4
	Internal loss, α loss	4.0 × 10 3 m -1
	Refractive index at zero injection, n 0	3.41351524
	Carrier density at transparency, N 0	1.23 × 10 24 m -3
	Differential index, d n/d N	-1.8 × 10 -26 m 3
	Group velocity, vg	8.33 × 10 7 m • s -1
	Nonlinear gain coefficient, ε	1.5 × 10 -23 m 3
	Structural parameters
	Active layer width, w	1.5 µm
	Active layer thickness, d	0.12 µm
	Cavity length, L	500 µm
	Optical confinement factor, Γ	0.35
	Grating period, Λ	227.039 nm

  simplified schematic diagram for the 3PS-DFB laser structure with non-equal spaced and non-identical PSs.stable SLM operation[7, pp. 128, 131]. After each step, the selection criteria are adjusted by fixing tighter limits, i.e., higher S min and smaller F max . The starting point is a QWS-DFB laser structure ϕ 2 = 90 • with κ L = 2. The restrictions PSP 1 ≤ PSP 2 and PSP 2 ≤ PSP 3 are assumed.The step-by-step procedure can be summarised as follows:Step 1: One PS, ϕ 1 , is added in the 1st half of the cavity. The optimisation of both functions S PSP 1 , ϕ 1 and F PSP 1 , ϕ 1 is performed by varying simultaneously and independently PSP 1 and ϕ 1 in their entire ranges: 0 ≤ PSP 1 ≤ 0.5 and 0

• ≤ ϕ 1 ≤ 180 • . The following selection criterium is used: S PSP 1 , ϕ 1 ≥ S min and F PSP 1 , ϕ 1 ≤ F max . This step leads to the definition of a region in the plane (PSP 1 , ϕ 1 ), from which a solution is chosen and new boundaries (S min , F max ) are settled; Step 2: Assuming PSP 1 and ϕ 1 achieved in Step 1, another PS, ϕ 3 , is placed in the 2nd half of the cavity. The optimisation of both functions S PSP 3 , ϕ 3 and F PSP 3 , ϕ 3 is performed by varying simultaneously and independently PSP 3 and ϕ 3 in their entire ranges: 0.5 ≤ PSP 3 ≤ 1 and 0 • ≤ ϕ 3 ≤ 180 • . The following selection criterium is used: S PSP 3 , ϕ 3 ≥ S min and F PSP 3 , ϕ 3 ≤ F max , where (S min , F max ) are the new boundaries settled in Step 1;

TABLE II -

 II S, F and α th L for several laser structures. κ L = 1.7 is assumed for the 1st structure. κ L = 2 is assumed for the 2nd and 3rd structures. th = 23.4 mA. This threshold current is similar to those reported in Ref.[7, pp. 159 and 162], for the QWS-DFB and the symmetric 3PS-DFB lasers, respectively, I th = 19.8 mA and I th = 21.8 mA.

	Laser structure	S	F	α th L
	Optimised 3PS-DFB	0.78	0.010	1.18
	QWS-DFB	0.73	0.30	0.70
	3PS-DFB PSP1 = 0.25 ; ϕ 1 = 60 • PSP3 = 0.75 ; ϕ 3 = 60 • PSP2 = 0.50 ; ϕ 2 = 60 •	0.34	0.012	0.98
	corresponding to I

Notice that, according to[START_REF] Kimura | Linewidth reduction by coupled-phaseshift distributed feedback lasers[END_REF],I(z) is numerically equal to E R (z) 2 +


Ari Sihvola

Six-dimensional view of dielectric mixtures as metamaterials

This paper discusses homogenization formulas and mixing rules that connect effective dielectric properties of mixtures with their microstucture. The special point of view in this paper is that all parameters appearing in mixing formulas are treated with equal status. The parameters that can be varied are the permittivities (which are allowed to be negative as well as positive) of the inclusion, environment, and composite, the volume fraction, spatial dimension, and an abstract parameter which chooses the given mixing formula type. The analysis helps in the interpretation of predictions from classical mixing rules when they are extended into the domain of negative-permittivity metamaterials.

Introduction

For understanding natural materials and designing engineered composites, electromagnetic research has provided homogenization theories and mixing formulas that give estimates for the macroscopic electric and magnetic properties of these materials [START_REF] Landauer | Electrical conductivity in inhomogeneous media[END_REF][START_REF] Neelakanta | Handbook of electromagnetic materials[END_REF][START_REF] Milton | The theory of composites[END_REF]. The classical use of mixing formulas has been to calculate the effective permittivity of a sample of heterogeneous medium, given the permittivies of the components and their fractional volumes in the mixture.

Methods to compute or even to estimate effective material parameters are needed in the remote sensing of natural environment where the dielectric response of geophysical media needs to be understood [START_REF] Ulaby | Microwave remote sensing, Active and passive[END_REF], in composite materials design [START_REF]Selected papers on linear optical composite materials[END_REF], and especially within the emerging field of metamaterials research [START_REF]Metamaterials: physics and engineering explorations[END_REF].

This presentation takes a new look at the homogenization problem, generalized in two respects. First, the restriction to the class of ordinary, conventional materials is relaxed, and the permittivities are allowed to reach also negative values as well as positive. Secondly, the usual way of looking at mixing rules as explicit formulas to calculate the effective permittivity when the properties of the mixture are given, is replaced by a multidimensional point of view. In other words, mixing rules are seen as connected relations where the parametrizations in all dimensions are placed on equal footing. The rank of this multidimensional world is six: these parameters are the three permittivities (environment phase, inclusion phase, effective medium), fractional volume, dimensionality of the space, and an abstract parameter ν connected to the type of mixing rule. On the other hand, in some other respects the degreees of freedom in the mixtures under consideration are restricted: there are only two component phases making the mixture, these materials are isotropic (their permittivity is equivalent to a scalar), and the underlying geometry is based on spherical building blocks.

Mixing relations

A great number of homogenization theories exist to explain the relation between the effective permittivity of a mixture ǫ eff as a function of the inclusion permittivity ǫ i , environment permittivity ǫ e , and the inclusion volume fraction f . A generalized mixing formula which attempts to capture a large amount of these homogenization principles, for the case of isotropic spherical inclusions, reads [START_REF] Sihvola | Electromagnetic mixing formulas and applications[END_REF]

ǫ eff -ǫ e ǫ eff + 2ǫ e + ν(ǫ eff -ǫ e ) = f ǫ i -ǫ e ǫ i + 2ǫ e + ν(ǫ eff -ǫ e ) , in 3 dimensions (1)
In this formula, the role of the unitless parameter ν is to pick various mixing rules which are known in the literature. Different choices of ν recover the following results: ν = 0 gives the Maxwell Garnett rule [START_REF] Maxwell | Colours in metal glasses and in metal films[END_REF], ν = 2 gives the Bruggeman (symmetric) formula [START_REF] Bruggeman | Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen, I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen[END_REF], and ν = 3 gives the Coherent potential mixing rule [START_REF] Kohler | Some applications of the coherent potential approximation[END_REF]. It is obvious that since the geometry of the mixture (spherical inclusions, random arrangement) and the component material responses are isotropic, also the effective continuum is isotropic and ǫ eff is scalar.

This relation (1) applies for three-dimensional (d = 3) spheres. Generalized into any spatial dimension d, the formula looks like

ǫ eff -ǫ e ǫ eff + (d -1)ǫ e + ν(ǫ eff -ǫ e ) = f ǫ i -ǫ e ǫ i + (d -1)ǫ e + ν(ǫ eff -ǫ e ) , in d dimensions (2)
Here the number of parameters has grown to six: the quantities are ǫ eff , ǫ i , ǫ e , f, ν, d.

The appearance of Equation ( 2) is certainly different from the way the well-known mixing formulas like Maxwell Garnett or Bruggeman are normally given. 1 To calculate the effective permittivity from this equation, a second-order equation needs be solved. The physically acceptable solution for this second-order equation has to be chosen with the correct sign of the imaginary part of ǫ ′ eff -jǫ ′′ eff : if the constituents are passive (ǫ ′′ e ≥ 0, ǫ ′′ i ≥ 0), the effective medium cannot be active (ǫ ′′ eff has to be positive). 2The only case when no roots are necessary in the explicit expression for ǫ eff is when ν = 0 which is the Maxwell Garnett formula. The solution with ν = 0 reads

ǫ eff = ǫ e + d f ǫ e ǫ i -ǫ e ǫ i + (d -1)ǫ e -f (ǫ i -ǫ e ) (3)
The advantage of the equation ( 2) is that it can be solved for any of the six parameters. As mentioned, the solution for ǫ eff (except in the Maxwell Garnett case) leads to a second-order equation. The situation is the same when ǫ e is solved. But for the other four parameters (ǫ i , f, d, ν), the solution is a first-order equation. This fact helps sometimes in computations for the effective permittivity when instead of directly solving for ǫ eff one may solve for f . The resulting (ǫ eff -f )-connection can be found out by varying ǫ eff between ǫ e and ǫ i and writing down the corresponding f . 3 It is true that with present-day computational methods to solve relatively simple equations like (2), such procedures to facilitate algebraic operations are not as essential as in earlier times. However, the idea of multidimensionality in a single formula and the possibility to view it from each parameter's point of view helps to glean physical interpretations of dielectric mixture behavior.

Application of the six-dimensionality

Indeed, the connections in Equation ( 2) offer many possibilities to study the characteristics of various mixing rules. In this section, the multidimensionality is applied to illustrate some interesting phenomena in the response of dielectric materials.

Percolation

One of such interesting examples is the phenomenon of percolation [START_REF]Percolation processes and structures[END_REF][START_REF] Grimmett | Percolation[END_REF][START_REF] Mclachlan | Electrical resistivity of composites[END_REF][START_REF] Stauffer | Introduction to percolation theory[END_REF]. In particular, percolation appears in different mixing rules in very varying ways. Percolation is a nonlinear phenomenon in the sense that a very abrupt change occurs in the macroscopic behavior. The geometry of the matter where percolation takes place is special: even with small changes in the fractions of the components forming the material, the structure behaves totally differently.

To see how the threshold of percolation comes from the mixing formula, let the dielectric contrast ǫ i /ǫ e in Equation ( 2) become large. In that case, the equation is approximately

ǫ eff -ǫ e ǫ eff + (d -1)ǫ e + ν(ǫ eff -ǫ e ) = f (4)
This gives the following expression for ǫ eff :

ǫ eff = ǫ e 1 + f (d -1 -ν) 1 -f (1 + ν) (5)
However, this relation is only valid as long as ǫ eff /ǫ e is small, i.e. only up to those volume fraction values where percolation appears. Clearly this result breaks down as the denominator reaches the value zero. That very point is interpreted as the percolation threshold point

f = f c = 1 1 + ν (6)
Depending on ν, the threshold varies between f c = 1 (Maxwell Garnett) and f c = 0 (ν = ∞). The well-known result of the percolation threshold of 0.333 for the Bruggeman symmetric formula is a simple corollary. Another noteworthy fact is that the threshold does not depend at all on the dimension of the space (d). 4A fine-tuned approach to percolating mixtures can be followed using the generalized mixing formula by McLachlan and others [START_REF] Mclachlan | Modeling the permittivity of composite materials with a general effective permittivity equation[END_REF]. The result (6), however, gives a straightforward connection between the mixing rule and percolation threshold.

Divergences of the effective permittivity

The recipe (2) is especially powerful in the case of metamaterials when the restriction on the permittivity to be positive is relaxed. Then the predictions for the macroscopic response are very strongly dependent on ν, in other words vary greatly with respect to different mixing rules. Singularities may appear, and complex effective permittivities may arise from real-valued components. Relation (2) can be solved for any of the six parameters, and the conditions for, for example, the singularities, can be calculated.

The solutions of (2) for ǫ eff read

ǫ eff = A ± √ A 2 -4ǫ e νB 2ν (7)
A = ǫ i [f (1 + ν) -1] -ǫ e [d -1 -2ν + f (1 + ν)] (8) B = ǫ i [f (1 -d + ν) -1] + ǫ e (1 -f )(1 -d + ν) (9)
This shows that there always exists a finite solution for ǫ eff for finite ǫ e , ǫ i , d, f, ν if only ν = 0. 5 The problematic case is the Maxwell Garnett formula which corresponds exactly to ν = 0, and is hence the only member of the mixing family which can predict divergences.

Let us take a closer look at the Maxwell Garnett infinities for ǫ eff . The divergence happens exactly when the dimension of the space is

d = (1 -f ) ǫ e -ǫ i ǫ e (10)
To have an infinity for ǫ eff requires of course that ǫ i and ǫ e have different signs, since otherwise condition [START_REF] Kohler | Some applications of the coherent potential approximation[END_REF] is not satisfied for any d ≥ 1 and 0 ≤ f ≤ 1.

Furthermore, the divergence requirement is stricter for higher densities. In other words, the larger f is, the higher the negative contrast ǫ i /ǫ e needs to be to cause an infinity in ǫ eff .

Mixtures of positive and negative permittivities

As a further example of peculiar predictions of various mixing formulas, let us analyze a mixture where the environment has a negative (relative) permittivity of ǫ e = -1, and the inclusions are air (ǫ i = +1). Let d = 3 which corresponds to three-dimensional spheres as inclusions. Figures 1 and2 display the predictions of the various models for ǫ eff of this mixture.

The figures show very interestingly the differing character of the various mixing rules. This variation is not obvious for "ordinary" mixtures with all the components having

Ν 0 Ν 1 Ν 2 Ν 3
Figure 1: The effective permittivity in the complex plane for a mixture where the background permittivity is ǫ e = -1 and the inclusion permittivity ǫ i = +1, as the volume fraction goes from 0 to 1. The predictions of the different mixing rules differ considerably: Maxwell Garnett (ν = 0) does not predict any losses, and for formulas with ν exceeding the values 0.5, the effective permittivity starts to display an imaginary part for a certain mixing ratio. Figure 2: The same as in Figure 1 for the real (solid) and imaginary (dashed) parts of the effective permittivity as functions of the volume fraction f . positive permittivities. The most striking observation is that even if in the calculated example both phases composing the mixture were real, many of mixing rules predict that the composite have complex effective permittivity. It has been pointed out earlier in the literature that the Bruggeman formalism predicts lossy continuum for nondissipative component media [START_REF] Mackay | A limitation of the Bruggeman formalism for homogenization[END_REF], along with a strong conclusion that the Bruggeman formalism is not physically applicable for two-component isotropic mixtures.

Another observation concerning the Bruggeman prediction (ν = 2) for the effective permittivity is the symmetry of the curve. The mixture under consideration is special (ǫ i = -ǫ e ), which leads to the fact that the 50-50 mixture (f = 0.5) gives ǫ ′ eff = 0, and the imaginary part has its maximum. The equal treatment of the host and guest materials in Bruggeman homogenization is a well-appreciated fact in the modeling of mixtures with positive-permittivity components.

Figures 1 and2 cover a wider range of mixing rules with respsct to this appearance of dissipation in homogenization. It is obvious that Bruggeman formula predicts losses for the present mixture under study. However, the figures show that Bruggeman is not the homogenization principle giving the highest amplitude in losses: the case ν = 1 has stronger losses (and numerically seems to be the lossiest of all mixing rules for this mixture). Furthermore, there is a range of mixing rules (in the framework of the classification of Equation ( 1)) which behave like Maxwell Garnett formula and predict ǫ ′′ eff = 0 for all possible volume fractions 0 ≤ f ≤ 1. This is the case for all mixing principles with ν < 1/2. Figures 1 and2 give support to this as far as the present mixture ǫ i = -ǫ e = 1 is concerned but it turns out to be valid also for any other values ǫ i = -ǫ e . 6In Figure 3, the limits are mapped for which parameter values a given mixing rule predicts a non-zero imaginary part for the effective permittivity. This describes the case when the relative permittivity of the inclusions is ǫ i = +1 (like in Figures 1 and2). The shaded area gives the background permittivity values for given volume fraction f for which ǫ eff will become complex.

As can be seen, the increase in ν has the effect that the shaded area increases. The Maxwell Garnett limit of vanishing ν can be seen as a zero-thickness complexity-causing area in the (f, ǫ e ) plane. But the cost of the vanishing area is that the effective permittivity is predicted to be infinitely large. For Maxwell Garnett, the area shrinks to a curve which matches exactly the divergence expression [START_REF] Kohler | Some applications of the coherent potential approximation[END_REF].

The complementary case is shown in Figure 4. There the corresponding map as in Figure 3 is drawn for the mixture where negative permittivity 3D spheres occupy a volume fraction f in air (ǫ e = +1). In other words, the role of the environment and inclusions are interchanged between Figures 3 and4. The area where ǫ eff becomes complex is again shaded. Also here, the increase of ν will make the shaded area larger. It is worth noting the behavior of the complexity regimes of the small-ν predictions: the f = 1 limit goes 0.2 0.4 0.6 0. Figure 3: A map for emergence of complex effective permittivity for a mixture where a negative continuum is pierced by three-dimensional spherical holes (with relative permittivity ǫ i = +1). The shaded area gives the range values of the background permittivity ǫ e and the volume fraction for which a non-zero imaginary part appears in ǫ eff for different mixing rules, characterized with the parameter ν.

to ǫ i = 0, whereas in the complementary case the negative-permittivity component was forced to go towards negative infinity.

Another peculiarity in Figure 4 is the behavior of the Coherent potential mixing formula (ν = 3). Figure 4 shows that for some volume fraction values, the shaded area extends to positive values (although very small values) of the inclusion permittivity. This is a characteristic that is absent from the complementary mixture, as is seen in Figure 3.

Conclusion

The six-dimensional description of the prediction of the effective permittivity of mixtures offers new possibilities to describe and analyze macroscopic properties of heterogeneous media. Especially in the case of metamaterial mixtures where the permittivity is allowed to extend into negative values, the descriptions (1) and (2) can simplify the description of emergent phenomena and even provide new insights into the properties of metamaterials.

It is true that when the possibility of negative permittivities is combined with the com- Figure 4: A similar map as in Figure 3, for the complementary mixture: the background permittivity is air (ǫ e = +1) and the inclusion permittivity varies in the negative regime. plications of mixing effects, sometimes counterintuitive predictions can be encountered, such as singularities, infinities, and complex-valued response parameters. For example Bruggeman mixing rule has been criticized for such behavior. It may be, however, a too strict conclusion to throw away the Bruggeman mixing rule for this type of positivenegative permittivity mixtures, only because of the dissipation that is created in the macroscopic permittivity. Maxwell Garnett rule as an alternative is also problematic with infinities for the continuum permittivity, which hurts everyday expectations. And these infinities are absent in mixing rules with ν = 0, like Bruggeman mixing.

Another reason to tolerate mixing rules with complexity-predicting properties is that, after all, mixtures like the ones considered in Figures 1 and2 are in the domain of metamaterials, having unconventional material parameters. And it is the very ethos of metamaterials that emergent properties appear in the higher-level homogenization. So why not accept emergent dissipation? As a final reminder, it is important to keep in mind the limitation to spherical geometry of the analysis. Mixing formulas have been written also for ellipsoidal inclusion shapes [START_REF] Sihvola | Electromagnetic mixing formulas and applications[END_REF]. However, to account for such geometries would require additional parameters into Equation (2).

Note, however, that the choice for Bruggeman formula depends on dimension: (for example, Bruggeman corresponds to ν

= 2 for

d = 3, and to ν = 1 for d = 2).[START_REF] Neelakanta | Handbook of electromagnetic materials[END_REF] Here the notation exp(jωt) is taken for the time-harmonic field excitation.[START_REF] Milton | The theory of composites[END_REF] Of course such a procedure returns the whole range of values 0 ≤ f ≤ 1 only in case of real permittivities, and with the assumption that ǫ eff always lies between the environment and inclusion permittivities which unfortunately is not always the case[START_REF] Sihvola | How strict are theoretical bounds for dielectric properties of mixtures[END_REF].

Because for the Bruggeman formula in two dimensions corresponds to ν = 1, the permittivity threshold is then

50% instead of 33%.

Sometimes the solution may be complex for real input parameters, due to the square root.

Also for negative values of ν, ǫ eff seems to remain real-valued. However, a strong conclusion about this claim would require further systematic study.


Michel Bruel Aplinov
email: michel.bruel@aplinov.com

Laser heating of thick layers through the backwards, self-sustained propagation of a steep and steady state thermal front

Keywords: 81, 20, n 81, 40, gh 81, 40, tv 81, 40, wx

We investigate theoretically a mode of heating layers using a laser beam. In this mode thermal energy, the temperature of the layer, propagates in a steady-state self sustained fashion from the bottom of the layer towards the surface and may exhibit a very steep front. The propagation exhibits a constant speed, related to the intensity of the power flux. To achieve this heating mode the absorption coefficient has to remain low in low temperatures and increase rapidly as a function of temperature in higher temperatures. Additionally, a significant temperature increase must be generated to trigger this propagation mode, for example through the presence of a strongly absorbing layer beneath the transparent layer. The mode is well suited to semiconductors, especially silicon. A few examples of applications are proposed.

INTRODUCTION

Laser processing of materials has been studied and used for many years [START_REF] Aplinov | Laser processing of materials[END_REF]. Powerful laser beams have been used with semiconductors as a tool for many different applications, such as annealing of implanted layers, alloying and etching [START_REF] Young | Laser processing of semiconductor materials[END_REF][START_REF] Baeri | Laser annealing of Silicon[END_REF]. The primary effect of applying a laser pulse to a material is heating. Laser energy is provided to the bulk material through the material's surface, leading to a profile of deposited energy peaking at the surface and decreasing with increasing depth. If the material is homogeneous and the absorption coefficient does not depend on temperature, the profile of this decrease is exponential. The thermal energy profile arising from this laser flux is governed by the so called Fick laws. The profile is time-dependant, decreases with increasing depth, and its extension within the bulk material evolves approximately following the square root of time.

In this paper, we present a theoretical evaluation of a new approach using a laser flux to heat up material layers, even very thick layers (several tens of microns). This evaluation is then confirmed by numerical simulation. It is required that the layer material does not absorb the chosen wavelength at the background temperature. Secondly, it is required that the absorption coefficient increases rapidly as a function of temperature in a range of higher temperatures. In these conditions, it will be shown that a self-sustained propagation of a thermal front with a steady-state shape is possible, and that the propagation occurs at a constant velocity (in case of constant power flux) from the depth towards the surface. In addition, the initial conditions required to trigger this phenomenon are shown, so that after a transient period the steady state propagation is reached. After the introduction, the second section of this paper covers the basic thermal aspects of a moving heat source. The third section shows that a heat source whose movement is selfsustained can be created inside a material and the fourth section presents a discussion on the limits of this theoretical work. The fifth section shows numerical simulations based on the conclusions of the first four sections, and we conclude with the sixth section, devoted to possible applications. In this paper, the temperatures will be referenced to the background temperature T 0 , for example, T stands for T-T 0 .

THERMAL ASPECTS IN CASE OF MOVING HEAT SOURCES

The theory of moving heat sources and associated steady thermal regimes have been long studied by different authors such as Rosenthal [START_REF]The theory of moving sources of heat and its applications to metal treatments[END_REF] and Carslaw and Jaeger [START_REF] Carslaw Andj | Conduction of heat in solids[END_REF]. More recent theories have also been developed by Erhard [6] and Modest [START_REF] Modest | Heat Conduction in a Moving Semi-infinite Solid Subjected to Pulsed laser irradiation[END_REF]. These theories generally deal with solving the heat equation in a semi-infinite solid, with a heat flux onto and through its surface that is moved by external means.Some aspects of this theory are presented below. We will consider an infinite homogeneous solid with rectangular co-ordinates XYZ bound to the solid (fig. 1). This solid is subject to a moving heat source which generates heat over a plane perpendicular to the x axis with a rate Φ (unit power per unit surface). The plane moves with velocity v in the direction of the x axis. If we consider the rectangular co-ordinates xyz referring to a fixed origin in the heat source plane and a given point within the solid, the relationship between the co-ordinates in the two systems are: x=X-v.t y=Y and z=Z with time origin taken when the two systems coincide.

The heat source can be represented by a Dirac delta function:

) 1 ( ) ( * ) ( 2 x x S D δ Φ =
It has been shown that in the moving xyz co-ordinates a steady state temperature profile (shape independent of time) is reached [START_REF] Carslaw Andj | Conduction of heat in solids[END_REF]. This is described by the following set of relationships (fig. 1):

) 2 ( / (2b) x/D) exp(-v (x) ........T 0 For x ) 2 ( T (x) .......T 0 For x max max 2D max 2D c Cv T with T a ρ Φ = ⎩ ⎨ ⎧ = > = <
where Φ is the incoming power flux, ρ is the density, C is the heat capacity, D =k/ρC is the thermal diffusion coefficient and k is the thermal conductivity. The profile consists of a plateau temperature and an exponential decrease.

Looking at the general case (steady state), the solid isn't submitted to a Dirac delta function δ(x), but is submitted to a heat source defined by S(x)= Φ * s(x) where s(x) is a physical function (approximately bell type) which is positive or null, continuous and whose integral equals to 1. From the previous set of relationships [START_REF] Young | Laser processing of semiconductor materials[END_REF], further referring to the Τ 2D function, one can derive the relationship between T and x for any kind of heat source.

Let us remember that s(x) can be considered as the convolution product of itself by the Dirac delta function.

) 4 ( )) ( * ( ) ( ) ( s * ) 3 ( ) ( ) ( ) ( s x x s x x x s x δ δ Φ ⊗ = Φ ⊗ =
As the heat equation is linear and as we have chosen constant coefficients, the relationship T(x) between the temperature and x coordinate is the convolution product of s(x) and the Τ 2D (x) function:

) 7 ( * ) ( * ) / ) ( ( * ) ( ) ( ) 6 ( * ) ( * ) ( ) ( (5) (x) ) ( ) ( 2 2D du u s du D u x v Exp u s x T du u x T u s x T T x s x T x x D ∫ ∫ ∫ ∞ ∞ - ∞ ∞ - + - - = ⇔ - = ⇔ ⊗ =
It can be seen that T(x) is the sum of two components .

du u s X * ) ( ∫ ∞
is the adiabatic part, the relationship between T(x) and s(x), when there is no heat diffusion within the solid (D=0).

For example, in the case of a gaussian distribution, s(x) is described by the following relationship:

) 8 ( ) 2 / ( 2 /( 1 ( ) ( 2 0 2 )) * 0 X x Exp X x s gauss - = π Which leads to : (
)

) 10 ( ) 2 / ( ) 2 / 1 ( ) / ) ( ( ) 2 / ( ) 2 / 1 ( / (9) (x) ) /2X (-x * ) 2 / 1 ( ) ( 2 0 2 0 2 0 2 0 max 2D 2 0 2 0 ∫ ∫ ∞ ∞ - - + - - - = ⊗ = x x Gauss Gauss du X u Exp X du D u x v Exp X u Exp X T T T Exp X x T π π π
With the reduced variables:

max 0 0 / ), 2 /( ) / ( ), 2 /( T T X v D X x Gauss Gauss = Θ = = π η π δ one gets ) 11 ( ) 2 / 1 * ( * ) / ( * ) 4 / 1 ( * ) 2 / 1 ( ) * ( * ) 2 / 1 ( 2 π η π δ η δ πη π δ + - - + = Θ Erfc Exp Exp Erfc Gauss π 2 0 X
represents the width of the source profile, defined as the width of a rectangular profile whose height equals the peak value of the actual source profile. D/v is homogeneous to a length and is an important parameter in all moving source issues. The η parameter defines the ratio between the length D/v and the width of the source profile. Figure 2 illustrates the evolution of Θ Gauss (δ,η) with η as parameter . The shape is sigmoidal; the lower parameter η (the smaller D/v compared to the width π 2 0

X

), the closer the temperature curve is to the deposited energy curve and the shorter the distance which is needed to reach the plateau temperature.

HOW TO CREATE A MOVING THERMAL SOURCE WITHIN THE SOLID

The next question is how to create, within the bulk of a solid, a thermal source extending in a steady-state bell shape fashion on both sides of a plane moving with respect to the solid, at a constant speed? The first approach is to imagine a way of creating a power flux from the outside of the solid, which can penetrate the solid and is able to deposit energy in a concentrated fashion at a given depth, changing the outside flux parameters in order to modify this energy deposition depth. A rough approximation could be to use energetic particles with small mass, such as protons, which deposit most of their energy near the end of penetration range of the particles and subsequently generate a concentrated thermal energy around the corresponding depth . Varying the particle velocity varies the depth at which the particles create this thermal source . However, this approach is only a rough approximation and is difficult to implement, especially if one needs to change the particle velocity very rapidly. The second approach, developed in this paper, consists of a situation where the power flux is still created from a source outside the material but the movement of the thermal source, the evolution of the depth at which the power flux is transformed into thermal energy, is no longer governed from outside the solid but is governed by the actual thermal source. In other words, in combination with the incoming flux, the thermal source provides the mechanism which moves the thermal source . The movement of the thermal source inside the solid is selfsustained. Let's consider a photon flux Φ coming from the right side which penetrates a homogeneous solid with an absorption coefficient α ; the power flux is absorbed according to the relationships:

) 13 ( )) ( exp( * ) ( ) ( ) 13 ( )) ( exp( * ) ( ) ( ) 13 ( )) ( exp( * ) ( ) ( ) 12 ( * ) ( / ) ( c x x x x b X X vt X vt X a X X X X dX X X d a a a a a a - Φ = Φ ⇔ - - Φ = - Φ ⇔ - Φ = Φ ⇔ = Φ Φ α α α α
The relationship between Φ and coordinates is the same in both systems of coordinates XYZ and xyz α depends on temperature, especially in semiconductors. As a result, one has to write:

) 14 ( ) ( / ) / ) ( ( )) ( ( ) 14 ( * )) ( ( ) ( / ) ( b x dx x d x T a dx x T x x d Φ Φ = ⇔ = Φ Φ α α ⎪ ⎩ ⎪ ⎨ ⎧ = Φ ⇔ Φ = ∫ ∞ + ∞ - ) 15 ( * ) ( ) 15 ( / ) ( 0 b du u S a dx d x S With ) 17 ( * ) ( / ) ( * ) ( / ) ( )) ( ( ) 16 ( * ) ( * ) ( ) ( 0 ∫ ∫ ∫ ∫ ∞ - ∞ - ∞ - ∞ = = = - Φ = Φ x x x x du u s x s du u S x S x T du u S du u S x α
Associating equation (17) with equation ( 5) provides a relationship between α and Τ.

⎪ ⎩ ⎪ ⎨ ⎧ = ⊗ = ∫ ∞ - x du u s x s x T T x s x T ) 17 ( * ) ( / ) ( )) ( ( (5) (x) ) ( ) ( 2D α
For example: let's consider the heat source whose shape (fig. 3) is described by the following set of equations:

⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ > - - <= <= - - < ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + + + - 1 ( ) ( 0 0 0 0 0 0 2 2 0 c X X for X X x X b X x X for X a X x for X X x X X x X
x s right eq eq lin left eq comp where X eq is set so that the integral of s(x) = 1 The heat source composite profile is a continuous curve with a plateau surrounded by an exponential decrease on the right and a gaussian curve ending as an exponential decrease on the left . It is possible to independently choose X 0 , X left, X right and X lin allowing different shapes. We will name T comp the temperature in the case of this composite profile.

.

) 19 ( ) ( ) ( ) ( T 2 comp X T X s X D Comp ⊗ =
This convolution product can be established either analytically or by numerical calculation. The following reduced variables are now used:

max 0 0 / , / ) / ( / , / , / , / , / T T X v D X X k X X k X X k X X k X x
Comp Comp eq eq lin lin eq eq left left eq right right eq

= Θ = = = = = = η δ
As defined in the former section, X eq represents the width of the source profile and η the ratio between D/v and the width of the source profile. Figure 4 displays the relationship between Comp Θ and the reduced abcisse δ . In this example k right = 0.7, k 0 = 0.15 , k left = 0.15 , k lin = 0.25 , η varying between 0 and 50 Figure 5 shows the relationships between absorption coefficient α and the temperature. The shapes of these curves show an increase of α with the temperature; one can also see that the smaller η is, the more rapid the increase. Coming back to the interpretation of the relationship between α and T, we consider the curve α(T) of fig 5 for η= 0.5, referring to it as α 0.5 (T).

Taking a material with density ρ, heat capacity C and whose absorption for a given wavelength as a function of temperature can be described by the relationship α 0.5 (T), we can irradiate this material through one of its faces with a light flux of the aforementioned wavelength with a power flux Φ. If good initial conditions are chosen, a steady state and self-sustained moving thermal regime will develop with a heat front whose shape is shown on figure 4, corresponding to η = 0.5, moving towards the left side with a velocity v, given by the equation

5 . 0 / ) / ( = = eq X v D η ,
and where the maximum temperature is given by T max = Φ/ρCv. It must be pointed out that as long as the flux Φ remains, the movement of the thermal front is self-sustained.

DISCUSSION

Here, we have shown that a self-sustained, steady-state and steep heat front can exist and propagate at a constant speed in the opposite direction to an incoming power flux.This propagation mode is possible provided that the absorption coefficient is an increasing function of temperature.

In addition, it must be stressed that the absorption coefficient is zero for T = 0. This is obvious when looking at equation ( 17). It becomes

) 20 ( 1 / ) ( * ) ( / ) ( ) 0 ( +∞ = +∞ = = ∫ ∞ ∞ - s du u s s T α must equal 0, otherwise the integral would diverge. ) (∞ s ∫ ∞ ∞ - du u s * ) (
As a result the pure ideal backwards and self-sustained propagation of a steep thermal front cannot exist in an actual physical situation, because it is physically unlikely that an absorption coefficient could strictly equal 0 at T = T0. Nevertheless, there are lots of physical cases which approximate the theoretical conditions so that it is possible to obtain backwards and self-sustained propagation of a steep thermal front with a quasi steady-state shape. This requires the existence of a temperature range at and above the background temperature where the absorption coefficient remains low. For a given experiment, the theoretical evaluation of the temperature increase at the surface at the end of the projected experiment, can provide a good criterion. The lower this temperature increase at the surface, compared with the expected plateau temperature, the closer the thermal profile evolution is to the theoretical case. The transient period,the duration necessary to reach the steady-state propagation, and the creation of conditions leading to a steady-state propagation after a transient period are also key issues from a physical point of view . In an infinite solid whose absorption coefficient is zero at its background temperature, nothing happens when irradiated by the power flux Φneither a steady-state heat front propagation, nor a temperature increase.

A solution is to create a layer called the triggering layer at a given depth, whose initial temperature is well above background temperature so that it is absorbing, or to create a layer where the absorption coefficient is not equal to zero at the background temperature. In both cases, the incoming flux will be absorbed by this triggering layer and will gradually transfer part of its thermal energy to the neighbouring material which in return will become more and more absorbing. The temperature profile will progressively change from the bell-shape (symmetric at t = 0 , and then more and more dissymmetric with an increasing tail) around the triggering layer to a moving steady-state steep thermal front moving at a constant speed towards the surface. The modelling of this transient state is highly complex and lies outside the scope of this paper. However, some general principles can be addressed, namely that the higher the initial temperature or absorption coefficient of the triggering layer, compared to the steady-state plateau temperature, the shorter the transient period and the smaller the distance for the temperature profile to be in a steady state, measured from the triggering layer upwards.

SIMULATION OF SILICON HEATING WITH A CO 2 LASER

The evolution of the thermal profile in bulk silicon irradiated through one of its faces by a CO 2 laser beam was simulated. The wavelength (10.6 µm) of this kind of laser is suitable because the absorption coefficient at room temperature in un-doped silicon is around 2 cm -1 (α L = lattice absorption) and is dependent on the temperature because its major contribution is the free carriers' absorption. The absorption coefficient can be described by [START_REF] Sheik-Bahae | Controlled CO2 laser melting of silicon[END_REF]:

) 22 ( * 10 * 9 . 1 ) ( ) 21 ( )) ( ( * ) ( ) ( 2 2 / 3 20 Cm T T Where T Ni Nd T T L - = + + = σ α σ α
The intrinsic density Ni(T) is given by: The variation of the gap energy with the temperature must be taken into account, it can be approximated by: Fig 6 shows the steep variation of the absorption coefficient as a function of temperature, calculated using equations 21 to 24 . It can be observed that there is a domain in the vicinity of the background temperature where the absorption coefficient remains quite low. Also, there is a domain of higher temperatures with a very strong increase.

The structural example which was chosen is essentially a highly doped (5E19cm

) 23 ( ) / ( * * ) ( 2 / 3 0 kT E Exp T n T N g i - = ) 24 ( ) 636 /( ) * 10 * 73 . 4 ( 17 . 1 2 4 + - = - T T Eg
-3 ) silicon substrate on top of which a very thick layer (150 µm thick), doped to 1E15 cm -3 , was grown. The role of the highly doped substrate is to operate as a trigger, as explained in the earlier section.

Computer software was written to solve the one dimensional heat equation based upon the finite differences method. This was used to simulate the evolution of the thermal profile. In these simulations, variations of heat capacity and thermal conductivity with temperature were not taken into account; constant values (1 J/g°C for heat capacity ,0.5 W/cm°C for thermal conductivity) were chosen. However, this simplification does not change the quality of the results. The incoming flux (it is coming from the left side ) was fixed at 2x10 6 W.cm -2 and the duration to 20 µs. The background temperature being chosen at 300K . These results are presented in figure 7; the temperature profile within the solid is represented at increasing times with step increases of 2 µs. The timing starts when the flux Φ is switched on. By looking at the results, we can see that the steady-state regime is obtained after a transient period of about 4 µs. This is very similar to the steady-state regime described in the previous sections and exhibits a progression speed of 7 µm/µs (700 m/s) with a plateau temperature of 1476 K, which corresponds to a temperature increase of 1176 K. We can also see the gradient of the front. In the linear section, the temperature rises from 600K to 1400K in a distance of approximately 4.5 µm. Inserting numerical values of Tmax, C, ρ and v into equation 2c confirms that the values from the simulation are coherent with the model: Φ/ρcv = 2x10 6 /(2.43x1x700) = 1176 = Tmax

APPLICATIONS

Three main characteristics of this method for heating layers must be emphasized:suitability for heating very thick layers, highly homogeneous heating throughought the thickness of the layer, the very fast process (in the range of ten tenths of a microsecond to some microseconds) and quasi adiabatic heating of the layer (almost no heating of underlying substrate). In the previous section, an example was given using a CO 2 laser beam at a wavelength of around ten microns. Other wavelengths can also work, especially those whose absorption by free carriers is the main mechanism, so that the absorption coefficient depends highly on the temperature. In the case of single crystal-silicon (low doped) for processing layers of medium thickness (for example, 20 microns), YAG lasers or Ytterbium doped fibre lasers at around 1.06-1.08 microns are suitable, provided that they can provide energy pulses of the required power and duration. They are less well adapted to very large thicknesses because of the magnitude of the absorption coefficient at room temperature, lying between 20 cm -1 and 300 cm -1 , depending on the authors. Using low background temperatures such as LN2 could improve the situation by reducing the band gap height, and thus reducing near band edge absorption.Erbium doped fibre lasers at around 1.5 µm and other fibre lasers at around 2 and 3 µm are interesting for single-crystal silicon, but would be particularly interesting for the heat treating of poly crystal materials, as absorption at the YAG wavelength could be too high, especially in the case of polysilicon. It is also possible to melt layers. One can melt silicon either across the entire thickness of the layer or across a portion of the layer thickness.Melting a material buried under a surface could be used for different applications, for example to segregate impurities in a deep region without affecting the surface layer or as a method to improve the crystalline quality of a buried layer by melting it and letting the re-solidification process proceed by epitaxial regrowth using the un-melted top-layer as the template. Melting an entire layer could be useful to re-crystallise non single-crystal layers, either epitaxially or not, in order to get either single-crystal layers or layers with larger grains. This could be useful in the field of photovoltaic materials.

CONCLUSION

It has been demonstrated that it is possible to heat thick ( from micrometers to several tens) layers in a quasi adiabatic fashionin a microsecond range at a very homogeneous plateau temperature. Some experimental work is currently underway to assess this new concept. The existence of high power fibre lasers and their further developments are also of great interest for applications based upon this new concept.

FIGURE CAPTIONS

Fig. 1 :Fig. 2 : 3 :

 123 Fig.1: Schematics of the moving heat source plane

Fig. 4 :

 4 Fig.4: Temperature profile as a function of the reduced abcisse in the moving coordinates system for the composite source profile corresponding to the set of parameters: k right = 0.7, k 0 = 0.15, k left = 0.15 and k lin = 0.25 . η is the parameter.

Fig. 5 :

 5 Fig.5: Relationship between absorption coefficient (normalized to 1) and temperature extracted from equations (17) and (5) for the composite profile with the previous quoted set of parameters. η is the parameter.

Fig. 6a :Fig. 7 .

 6a7 Fig.6a: Absorption coefficient of silicon versus temperature at the 10.6 µm wavelength (linear scale)

  Fig.1

Fig. 5

 5 Fig. 6b

ACKNOWLEDGMENTS

Thanks to P. Mottier for establishing the core of the finite differences software in 2003.Thanks to Philip Ewels for improving the english . .