Aakashdeep Gupta

K Nidhin

Member, IEEE Suresh Balanethiram
email: suresh.b@nitpy.ac.in.

Shon Yadav

Sebastien Fregonese
email: sebastien.fregonese@ims-bordeaux.fr

Senior Member, IEEE Thomas Zimmer
email: thomas.zimmer@ims-bordeaux.fr.

Member, IEEE Anjan Chakravorty

S Yadav Is With Globalfoundries

Bangalore S Fregonese

Optimizing Finger Spacing in Multi-Finger Bipolar Transistors for Minimal Electrothermal Coupling

Keywords: SiGe HBT, multi-finger transistor, finger placement, self-heating, thermal coupling, shallow trench isolation, Kirchhoff's transformation

We present a compact modeling framework to optimize finger spacing for improving the thermal stability in multi-finger bipolar transistors with shallow-trench isolation. First, we present an accurate physics-based model for total junction temperature in all the fingers of a transistor. Other than validating the model with 3D TCAD simulations and measured data, we demonstrate its efficacy to achieve finger spacing optimization with the aid of an iterative algorithm. Since the proposed technique is scalable from the viewpoint of the number of fingers within a transistor and their geometries, the proposed framework is found to work seamlessly for various emitter finger numbers.

I. INTRODUCTION

Multi-finger silicon-germanium (SiGe) heterojunction bipolar transistors (HBTs) are preferred in power amplifier applications [START_REF] Erben | Class-a sige hbt power amplifiers at c-band frequencies[END_REF]- [START_REF] Andrews | Compact modeling of mutual thermal coupling for the optimal design of SiGe HBT power amplifiers[END_REF]. A conventional layout of a multi-finger transistor always consists of uniformly spaced emitter fingers with a fixed emitter width (W E ) and emitter length (L E ) [START_REF] Erben | Class-a sige hbt power amplifiers at c-band frequencies[END_REF] (see Fig. 1(a)). High current operations of power amplifiers cause electro-thermal heating inside the device leading to a high junction temperatures [START_REF] Rinaldi | Theory of electrothermal behavior of bipolar transistors: part II-Two-finger devices[END_REF]. Since all the fingers dissipate heat simultaneously, the total junction temperature of a given finger is resulted from the self-heating (due to power dissipation within the finger) and thermal coupling (from power dissipation at adjacent fingers) effects [START_REF] Walkey | Equivalent circuit modeling of static substrate thermal coupling using VCVS representation[END_REF], [START_REF] Zimmermann | Modeling of spatially distributed and sizing effects in high-performance bipolar transistors[END_REF]. Note that all the fingers are electrically isolated but thermally coupled via substrate. The unequal amount of thermal coupling at each finger of the conventional layout (Fig. 1(a)) leads to an uneven total temperature distribution across the fingers and one observes a comparatively higher temperature in the central (or inner) fingers [START_REF] Gao | Thermal design studies of high-power heterojunction bipolar transistors[END_REF]. This can lead to degraded electrical performances and various issues due to electrothermal instability such as snapback and thermal runaway [START_REF] Spina | Thermally induced current bifurcation in bipolar transistors[END_REF] eventually affecting the safe-operating-area and reliability of the device [START_REF] Rinaldi | Theory of electrothermal behavior of bipolar transistors: Part III-Impact ionization[END_REF]. One can overcome these issues with nonuniform finger spacing in a multi-finger transistor structure as shown in Fig. 1(b). A careful optimization of finger spacing can help obtaining nearly uniform temperature at all the fingers. Such an optimization can be carried out either using computationally expensive numerical simulations [START_REF] Metzger | Evaluation of thermal balancing techniques in ingap/gaas hbt power arrays for wireless handset power amplifiers[END_REF], [START_REF] Wang | Boosting up performance of power sige hbts using advanced layout concept[END_REF] or by the use of complex analytical formulation [START_REF] Wan-Rong | Improvement of thermal stability of multifinger power SiGe HBTs using emitter-ballasting-resistor-free designs[END_REF]. A simple, fast and accurate analytical model to solve this optimization problem is not reported so far in the literature.

Existing electrothermal compact models of multi-finger transistors use superposition to compute the rise in the junction temperature of i th finger (∆T i ) above an ambient temperature (T a ) as [START_REF] Walkey | Equivalent circuit modeling of static substrate thermal coupling using VCVS representation[END_REF] [6] [START_REF] Gupta | Static Thermal Coupling Factors in Multi-Finger Bipolar Transistors: Part I-Model Development[END_REF] [14]

∆T i = ∆T ii + n j=1,j =i ∆T ij = P d,i R th,ii + n j=1,j =i c ij P d,j R th,jj (1)
where the first term indicates the self-heating at i th finger and the second term accounts for the thermal coupling on i th finger from all other fingers. P d,i (P d,j ) and R th,ii (R th,jj ) are, respectively, the electrical power dissipation and thermal resistance corresponding to the i th (j th ) finger. Static thermal coupling coefficient (c ij ) quantifies the amount of mutual thermal coupling over finger-i due to the heat source at fingerj and is defined as

c ij = ∆T ij ∆T jj . (2)
However, direct use of (1) on the measured data to extract thermal coupling related model parameters yields erroneous results as reported in [START_REF] Gupta | Extraction of true finger temperature from measured data in multifinger bipolar transistors[END_REF]. The study clearly states that the validity of ( 1) is limited only to a linear case where the thermal conductivity (κ) of a given substrate material is constant. However, the thermal conductivity of substrate Si is temperature-dependent [START_REF] Nghiêm | Electro-thermal simulation based on coupled boltzmann transport equations for electrons and phonons[END_REF]; hence the governing heat diffusion equation becomes non-linear. Therefore, one cannot use (1) as a basis formulation in order to optimize the finger spacing for obtaining a nearly uniform temperature distribution in a multi-finger transistor.

In this paper, we present a simple analytical compact model framework along with an iterative algorithm to obtain a fast and logical solution of this optimization problem. Section II describes modeling framework, followed by model validation with 3D TCAD simulation and measurement data in section III. In section IV we present the results and discussions on spacing optimization using the proposed model and finally conclude in section V.

II. MODELING FRAMEWORK

Figs. 1(a) and (b) show the top view of a shallow trench (ST) isolated five-finger structure with uniform (s) and non-uniform spacing (s i with i =1, 2, 3, 4), respectively. In Fig. 1(b), the spacing between the inner fingers (s 2 and s 3 ) are increased to reduce the total thermal coupling and the total temperature rise at the central finger maintaining a constant the foot print as in Fig. 1(a). Essentially, we kept (n -1) s= n-1 k=1 s k = d, with n as the number of fingers. The problem is to quickly find out the optimized values of [s 1 , s 2 , s 3 , s 4 ] in order to ensure a nearly equal and possibly minimal junction temperature at all the fingers. Note that for every unique spacing combination there exists a unique set of temperature values of all the fingers. We start with an accurate estimation of total temperature rise of a given finger due to self-heating and thermal coupling for any given set of spacing values.

A. Linear self-heating Fig. 2 shows a cross-sectional view of the first three fingers in an n-finger shallow trench (ST) isolated structure with a substrate thickness of H and spacing between the immediate neighbours as s 1 and s 2 . The heat-flow volume of the heating finger (i.e., j = 1 in Fig. 2) is defined with the thermal spreading angles (θ and θ). The z-dependent temperature profile T (z) under the heating finger can be expressed as [START_REF] Rieh | Measurement and modeling of thermal resistance of high speed SiGe heterojunction bipolar transistors[END_REF]

T (z) = T a + P d R th (z) (3)
where R th (z) is z-dependent thermal resistance at an ambient temperature T a and is expressed as

R th (z) = [1/κ(T a )][ z 0 dz/A(x, y, z)] = f G (z)/κ(T a )
. The geometry factor f G (z) is defined in terms of position-dependent heat-flow cross-sectional area, A(x, y, z) = (L E + 2z tan(θ z ))(W E + 2z tan(θ z )), as [START_REF] Balanethiram | Accurate modeling of thermal resistance for on-wafer SiGe HBTs using average thermal conductivity[END_REF]

f G (z) = z 0 dz A(x, y, z) = ln L E (W E +2z tan(θz)) W E (L E +2z tan(θz)) 2 (L E -W E ) tan(θ z ) (4)
where W E and L E are, respectively, the width and length of the heat source. Note that the heat spreading angle θ(z) = θ within ST and θ(z) = θ outside ST. We evaluate (3) to estimate the substrate temperature at any z by considering a constant thermal conductivity of the substrate material (κ(T a )) measured at T a . This yields a linear junction temperature rise at the heat-source (∆T jj,lin = T (z = 0) -T a ) due to the self-heating effect.

B. Linear thermal coupling

It is evident that the junction temperature of the neighboring fingers (i=2, 3, .., n) tends to rise above T a (by an amount ∆T ij,lin ) due to thermal coupling from the heating finger (j=1). In [START_REF] Gupta | Static Thermal Coupling Factors in Multi-Finger Bipolar Transistors: Part I-Model Development[END_REF] we reported a strategy to model the coupling coefficients (c ij ) in device structures without any trench isolations and also for the ST-isolated devices. Although the formulations derived in [START_REF] Gupta | Static Thermal Coupling Factors in Multi-Finger Bipolar Transistors: Part I-Model Development[END_REF] considered κ(T ), the approximations

d Fin-1 s s s s Fin-5 W E L E (a) s 1 s 2 s 4 s 3 ST (b)
A E = W E × L E ). θ DT (a) x z (b) T (zn) T a T (z3) T (z1) T (z2) θ T (zn) T a T (z3) T (z1) T (z2) z = 0 z = H z = 0 T jj T jj T 21 z = H s 1 s 2 W E L E θ ′ T (z4) T (z5) T (z3)

ST ST Si Substrate

Heating finger j = 1

Si Substrate

Sensing finger i = 2, 3, ...n T 31 T 21 T 31 T (z3 = x = s1) = T21 T (zn = x = s1 + s2) = T31
Fig. 2. Schematic cross-section of a five-finger multifinger transistor structure with ST isolation. The heat source is located at z = 0. used in the derivation hold true even for the linear case with constant κ(T a ). The resulting c ij in a no-trench (NT) structure is expressed as [START_REF] Goh | A simple extraction method of thermal resistance in multifinger GaAs HBT[END_REF]

c i j | S T [ % ] ∆T i j , l i n | N T / ∆T i j , l i n | S T T a = 3 0 0 K W E = ( 0 . 2 , 0 . 4 , 0 . 6 , 0 . 8 ) µm x [ µm ] ( a ) ( b ) W i t h -S T T a = 3 0 0 K L E = ( 5 , 1 0 , 1 5 ) µm W E = 0 . 2 µm L E = ( 5 , 1 0 , 1 5 , 2 0 ) µm T C A D f r o m T ( x ) M o d e l ( 3 ) , ( 5 ) a n d ( 6 )
c ij | N T = ∆T ij,lin | N T ∆T jj,lin | N T ≡ T N T (z = x = s ) -T a T N T (z = 0) -T a ( 5
)
where s is the spacing between the heating and sensing finger. Note that the c ij | N T in (5) comes out to be P d -independent due to linear T N T (z). Subsequently, following [START_REF] Gupta | Static Thermal Coupling Factors in Multi-Finger Bipolar Transistors: Part I-Model Development[END_REF] one can express the coupling factor for the structure with-ST as

c ij | ST = c ij | N T ∆T ij,lin | ST ∆T ij,lin | N T ∆T jj,lin | N T ∆T jj,lin | ST ≈ c ij | N T ∆T jj,lin | N T ∆T jj,lin | ST . ( 6
)
Note that all the terms in ( 5) and ( 6) can be directly obtained using (3) along with (4). 4) to obtain the solid lines. From the results of Fig. 3, one can confirm the efficacy of ( 5) and ( 6) in predicting the thermal coupling between a given pair of heating and sensing fingers with any finger spacing.

C. True junction temperature and spacing optimization

The formulations presented in II-A and II-B are linear due to temperature independent thermal conductivity and provide us with the information of temperature rise caused by selfheating and thermal coupling at a given i th finger. The linearity allows one to employ the superposition relation of ( 1) to obtain ∆T i,lin from the sum of ∆T ii,lin and ∆T ij,lin (for all j). However, the substrate thermal conductivity is temperaturedependent in the practical range of operation and follows a relation κ(T ) = β/T α [START_REF] Nghiêm | Electro-thermal simulation based on coupled boltzmann transport equations for electrons and phonons[END_REF]. Following the work in [START_REF] Gupta | Extraction of true finger temperature from measured data in multifinger bipolar transistors[END_REF], we apply Kirchhoff transformation over the already estimated ∆T i,lin to obtain the true total junction temperature as

T i,true = T a 1 + (1 -α) ∆T i,lin T a 1/(1-α) . (7)
Note that α and β are substrate material parameters to model the temperature dependent thermal conductivity. In order to obtain a minimum peak temperature with maximum possible uniformity in the temperature distribution of the multi-finger transistor structure, we embed the above developed model [START_REF] Gao | Thermal design studies of high-power heterojunction bipolar transistors[END_REF] into an optimization algorithm shown in the form of a flow chart in Fig. 4. In the main program we first start by defining the minimum (s min ) and maximum sweeping (s max ) limits for all the spacing variables (s i ). Note that for a practical purpose, s min and s max can be, respectively, set to s uni /2 and 3s uni /2 where s uni = d/n with d as the total spacing between the corner fingers and n as the total number of fingers. Following this, we employ nested loops to generate all the possible combinations of finger spacing variables. In a given iteration, we store a set of spacing variables into an array S. At the if-decision block, we check and select a given S if all its elements add to the fixed value of the parameter d. The main program then calls the T i,true model ( 7) to compute the total junction temperatures at all fingers and stores them in the array T i,true for a selected S. We further store the maximum value of the array T i,true as T i,peak . The algorithm yields optimized set of finger spacing that is associated with the minimum value of T i,peak . All the geometry details and material properties are fed into the program as model parameters.

III. MODEL CALIBRATION USING STRUCTURES WITH UNIFORM FINGER SPACING

A. Calibration with TCAD simulation

We created two types of five-finger transistor structures with silicon substrate in TCAD environment of Synopsys Sentaurus [START_REF] Synopsys | Sentaurus: Sentaurus device user guide, release H-2013[END_REF], one with no-trench and the other with-ST isolation. Uniform finger spacing was kept at d/5 = 2.5 µm. The dimensions of the structures correspond to those of the state-of-the-art SiGe HBTs from STMicroelectronics B55 technology [START_REF] Chevalier | A 55 nm triple gate oxide 9 metal layers SiGe BiCMOS technology featuring 320 GHz f T /370 GHz fmax HBT and high-Q millimeter-wave passives[END_REF]. An emitter area of 0.2x5 µm 2 is maintained for both the structures. We performed 3D thermal simulations of these TCAD structures where all the fingers were heated simultaneously in order to emulate a real operating condition. Note that the substrate material thermal parameters with α = 1.261 and β = 2.076 × 10 5 WK α-1 /m were used in the simulation.

Figs. 5 (a) and (b) compare the P d -dependent T i,true obtained from (7) (solid lines) and the TCAD (symbols) for a five-finger transistor with no-trench and with-ST isolation, respectively. For f G (z) calculation in (4) we used heat spreading angles of 35 o and 48 o , respectively, for within and outside the ST region. The results of our modeling framework show excellent agreement with the TCAD simulation results. Two sets of TCAD data with red circles and black cross (and the respective modeling results) correspond to the 5-finger structures with uniform spacing (s=2.5 µm) and optimized spacing, respectively. Results show that our model predicts the TCAD simulation with high level of accuracy; hence it can be used to optimize the finger spacing for minimizing the overall electrothermal heating instead of rigorous TCAD simulation. In order to compare the maximum achievable accuracy of state-of-the-art existing model, we have also added the total temperature values resulting from (1) (green symbols) after substituting the non-linear TCAD simulated values of ∆T jj and ∆T ij in the case of P d =30 mW. It is observed that the existing model noticeably underestimates the true finger temperature. Optimized spacing related results are discussed in section IV.

B. Calibration with Experimental data

To validate the proposed framework for a fabricated device, we chose the on-wafer measurement data presented in the work of [START_REF] Lehmann | Characterization of the static thermal coupling between emitter fingers of bipolar transistors[END_REF]. Thermal characterizations were carried out for shallow trench-isolated five-finger HBTs from state-of-theart technology of IHP microelectronics. Devices with emitter dimensions of W E × L E = 0.44 × [7.64, 12.68, 27.8] µm 2 and uniform finger spacing of s = 3.23 µm were considered in the study. Following [START_REF] Gupta | Extraction of true finger temperature from measured data in multifinger bipolar transistors[END_REF] we extracted the true junction temperature (T i,true ) values for all the on-wafer devices and used them to validate our model. Note that in a practical device one should account for the upward heat flow through the backend-of-line (BEOL) metallization layers, and also consider the variation in thermal conductivity due to doping and the materials in the heat flow path. The upward heat flow is modeled by a finite BEOL thermal resistance (R th,m ) that appears in parallel with the net substrate thermal resistance R th,s to produce the overall self-heating thermal resistance (R th ) seen from the heating finger. Before the model validation step we used the the extraction strategy detailed in [START_REF] Balanethiram | Extraction of BEOL contributions for thermal resistance in SiGe HBTs[END_REF] that employs the concept of average thermal conductivity of the substrate material to simultaneously extract the substrate κ(T ) related parameters α and β, and the BEOL related R th,m parameter by utilizing multiple measured data points of self-heating junction temperature (T jj ) and thermal resistance (R th,jj ). Since the extracted thermal conductivity parameters are obtained directly from experimental data, they automatically include the effects of substrate material and doping on to the overall thermal behavior. Fig. 6 compares the results of T i,true obtained from our model (solid lines) with those extracted from on-wafer measurements (symbols) for all the three emitter geometries. Excellent model agreement is observed across different emitter geometries. Increased thermal coupling for the inner fingers than the corner fingers results in higher temperature for fingers-3 and 2 than finger-1. R th,m values of 20, 15, 12 (in kK/W) were extracted for the devices with L E of 7.64, 12.68, 27.8 (in µm) respectively. Similarly, the extracted values of α = 1.285 and β = 2.136 × 10 5 WK α-1 /m and they were used for model validation. The inset in the figure demonstrates the highest achievable accuracy in predicting the total junction temperature at finger-3 for the device with L E =7.64 µm by a direct use of superposition of the measured self-heating and coupling temperature data (state-of-the-art model). The prediction becomes significantly erroneous with increasing dissipated power.

IV. FINGER SPACING OPTIMIZATION RESULTS

Excellent model capability and its calibration with TCAD simulation and measured data are demonstrated in section III.

Therefore, we can now apply our model for the purpose of finger spacing optimization in order to ensure minimal over device temperature due to electrothermal heating. We implemented the finger spacing optimization algorithm detailed in section II (see Fig. 4). The algorithm was executed for TCAD multi-finger structures based on STMicroelectronics, as mentioned in section III-A. Appropriate values of model parameters are inserted and optimized spacing values are obtained via iterative calculations. Fig. 5 presents the T i,true values (in black) obtained from TCAD (cross symbols) and model [START_REF] Gao | Thermal design studies of high-power heterojunction bipolar transistors[END_REF] (solid lines) for optimized five-finger structures. The data in black are obtained for optimized set of spacing with values S = [START_REF] Liu | Thermal coupling in 2-finger heterojunction bipolar transistors[END_REF][START_REF] Andrews | Compact modeling of mutual thermal coupling for the optimal design of SiGe HBT power amplifiers[END_REF][START_REF] Andrews | Compact modeling of mutual thermal coupling for the optimal design of SiGe HBT power amplifiers[END_REF][START_REF] Liu | Thermal coupling in 2-finger heterojunction bipolar transistors[END_REF] µm. Note that the corner finger temperatures in optimized structures are more even with lower peak temperature than in uniformly spaced structures (data in red). The optimization routine attempts to obtain a relatively uniform temperature across the fingers. However, for a given footprint where thermal coupling is inevitable, non-uniformity in spatial temperature cannot be completely eliminated. Our algorithm that utilizes analytical formulations, essentially attempts to minimize this non-uniformity. In Fig. 7, for a given emitter geometry and P d we present T i,true values corresponding to the optimized finger spacing obtained from the optimization algorithm for six and seven-finger structures with-ST. Finger temperatures of uniformly spaced structure are also plotted for a direct comparison. Since the thermal coupling at the corner fingers Fin-1 and Fin-6 (or Fin-7) is one sided (either from right or left), their optimized T i,true values are the lowest among all the fingers. The optimized trend of total temperature demonstrated in Fig. 7 is similar to the previously reported ones in [START_REF] Wang | Boosting up performance of power sige hbts using advanced layout concept[END_REF], [START_REF] Wan-Rong | Improvement of thermal stability of multifinger power SiGe HBTs using emitter-ballasting-resistor-free designs[END_REF]. The optimized spacing values obtained here are independent of the power dissipation at the fingers. Fig. 8 presents the P d -dependent T i,true of all four fingers in an optimized seven finger ST-isolated structure with a particular emitter geometry. The data points obtained from the model (solid line) and TCAD simulation (symbols) corresponding to the same finger spacing are also plotted for validation purpose. We have also plotted T i,true values of the central (4 th ) finger from the uniformly spaced structure to highlight the temperature reduction due to optimization.

V. CONCLUSIONS

This work provides a simple finger spacing optimization framework for thermally aware design of multi-finger transistors with shallow trench isolation. The framework comprises of an iterative optimization algorithm and a new physics-based model for accurate prediction of the total junction temperature. The model includes the linear finger temperature rises due to self-heating and thermal coupling effects and employs Kirchhoff's transformation based expression to predict the true total temperature for all the fingers of a transistor. Model is first validated and calibrated against 3D TCAD simulations and subsequently against experimental data after incorporating the effect of back-end-of-line. Excellent model agreement is observed with the values of true finger temperatures from TCAD and the ones extracted from the measurement data of fabricated devices. The model is implemented in an optimization algorithm to achieve the optimized non-uniform finger spacing values in order to minimize the overall nonuniformity in the spatial temperature distribution across the fingers. Obtained values of optimized finger spacing are found to be independent of power dissipation and give rise to an almost uniform temperature distribution across all the fingers of a multi-finger device. The proposed framework is scalable for transistors with different emitter dimensions and any number of fingers.

Fig. 1 .

 1 Fig. 1. Top views of the shallow trench isolated five-finger structures with (a) uniform spacing (of s) and (b) non-uniform spacing. A fixed value of d ensures the same total footprint area for both the structures. Areas shaded in grey and red, respectively, indicates ST isolation and emitter fingers (of area A E = W E × L E ).

Fig. 3 .

 3 Fig. 3. (a) Linear temperature ratio at the sensing fingers of transistors with no-trench to with-ST ( ∆T ij,lin N T / ∆T ij,lin ST ) obtained from TCAD for different W E at a fixed L E = 5 µm and for different L E at a fixed W E = 0.2 µm. (b) A comparison of c ij between TCAD (symbols) and values resulting from the model (solid-lines) of c ij | ST for the case when 1 st finger is heating and rest fingers are sensing. A constant finger spacing of s = 2.5 µm is used in the structures.

  Fig.3(a) demonstrates the validity of the approximation ∆T ij,lin | N T / ∆T ij,lin | ST ≈ 1 (as used in (6)) for various emitter geometries and at a fixed κ(T a ). Since the shallow trench (ST) depth (< 0.4 µm) is usually much less than the finger spacing (around 2.5 µm), the overall coupling temperature in ST-isolated device hardly differs from that of the no-trench (NT) device. One obtains P d -independent c ij | ST from (6) and subsequently ∆T ij,lin | ST using (2). In Fig.3(b) we compare the results from (6) (solid lines) against the TCAD values (with fixed κ(T a )) for five-finger uniformly spaced structures with spacing of 2.5 µm and different emitter geometries. Note that we employ heat spreading angle of 35 o (for volume confined by ST) and 48 o (for volume below ST) in (

Fig. 4 .

 4 Fig. 4. The complete flow chart of the algorithm to obtain the optimized non-uniform finger spacing. The algorithm is valid for any number of fingers and employs true total temperature model as a function. All the geometry details and thermal conductivity parameters are inputs.

Fig. 5 .

 5 Fig. 5. Comparison of total temperature distribution obtained from model (solid lines) and TCAD (symbols) in (a) no-trench and (b) with-ST structures for all the fingers. Data points in red and black result from structures with uniform and optimized spacing, respectively. Green symbols represent existing state-of-the-art model obtained from (1) after substituting non-linear TCAD simulated values of ∆T jj and ∆T ij for uniformly spaced structure.

Fig. 6 .

 6 Fig. 6. True finger temperature of fingers-1, 2 and 3 for a fabricated shallow trench isolated structures with W E = 0.44 µm and L E = [7.64, 12.68, 27.8] µm: comparison of the proposed model against the extracted data from measurements. Inset shows the underestimation in predicting the T i,true values obtained by a direct use of superposition of measured selfheating and coupling temperature data for the central finger.

Fig. 7 .O p t i m i z e d s t r u c t u rFig. 8 .

 78 Fig. 7. Comparison of total temperature distribution obtained from the proposed algorithm in Fig. 4 when applied to ST isolated structures with uniform finger spacing (red circles) and that with optimized finger spacing (black triangles): (a) for a six-finger device with optimized spacing S = [1.73, 3.16, 2.73, 3.16, 1.73] µm and (b) a seven-finger device with optimized spacing S = [1.72, 2.81, 2.97, 2.97, 2.81, 1.72] µm .

This work was supported in part by the EU under Project Taranto under Grant 737454, in part by ISRO project ELE/17-18/176/ISRO/ANJA and in part by DST, India, under Project EMR/2016/004726.


Member, IEEE Suresh Balanethiram
email: sureshbalanethiram@gmail.com.

Member, IEEE Deleep R Nair

Rosario D ' Esposito

Senior Member, IEEE Nihar R Mohapatra

Sebastien Fregonese
email: sebastien.fregonese@ims-bordeaux.fr

Senior Member, IEEE Thomas Zimmer
email: thomas.zimmer@ims-bordeaux.fr.

Member, IEEE Anjan Chakravorty

BEOL Thermal Resistance Extraction in SiGe HBTs

Keywords: SiGe HBTs, self heating, thermal resistance, back-end-of-line (BEOL), parameter extraction, compact models

A prior estimate of the impact of thermal resistance from the back-end-of-line (BEOL) metallization layers is crucial for an accurate circuit design and thermally aware device design. This paper presents a robust technique to extract the thermal resistance component originating from the BEOL metal layers in silicon germanium heterojunction bipolar transistors (SiGe HBTs). The proposed technique is first tested on data generated using analytical equations and later validated with 3D TCAD simulation. The results clearly show that the exact contribution of the BEOL to the overall thermal resistance is captured in the proposed approach. Finally, we verified the method using measured data obtained from fabricated SiGe HBT structures using Infineon B11HFC technology. The extracted parameters show reasonable accuracy and consistency across different emitter dimensions and BEOL configurations.

I. INTRODUCTION

A serious concern about the modern silicon germanium heterojunction bipolar transistors (SiGe HBTs) is the selfheating induced degradation of high-frequency performance and long-term reliability of the device [START_REF] Raghunathan | Bias-and temperature-dependent accumulated stress modeling of mixed-mode damage in SiGe HBTs[END_REF]. This is why a self-heating aware device design has attracted the attention of researchers [START_REF] Rieh | Structural dependence of the thermal resistance of trench-isolated bipolar transistors[END_REF]- [START_REF] Rieh | Scaling of SiGe heterojunction bipolar transistors[END_REF]. On the other hand, to ensure a reliable circuit design, the transistor model is supposed to accurately predict the electrothermal effect using an appropriate thermal sub-circuit. The static thermal network essentially involves two connected elements, a dependent current source equal to the total dissipated electrical power (P d ) and a thermal resistance (R th ). The former is directly obtained from the electrical network, while the latter depends on the transistor structure and material thermal conductivity. Fig. 1 shows a typical transistor structure including the front-end-of-line (FEOL) Si substrate and back-end-of-line (BEOL) metal layers. Since the power dissipation happens mostly at the base-collector junction, the generated heat flows towards both the FEOL substrate and BEOL stack composed of metal layers and inter-layer dielectrics. Therefore, the overall R th has two components in parallel, the FEOL component (R s ) and the BEOL component (R m ) as shown in Fig. 1. Note that the resistor R m additionally includes the thin SiGe base as well as mono and poly-emitter portions (part of upward heat flow). As the thermally insulating inter-layer dielectrics constitute bulk of the BEOL portion along with vias and contacts, a small but non-negligible amount of heat flows upwards [START_REF] Esposito | Dedicated test-structures for investigation of the thermal impact of the BEOL in advanced SiGe HBTs in time and frequency domain[END_REF], [START_REF] Esposito | Thermal penetration depth analysis and impact of the BEOL metals on the thermal impedance of SiGe HBTs[END_REF]. Note that for electrical parameter extraction purposes, special test structures are designed where the emitter is grounded leading to a higher amount of upward heat flow. In case of convective boundary at the chip-ambient interface, thermal ground at the top surface in Fig. 1 assumes that the extracted BEOL thermal resistance additionally includes the effect of the convective boundary condition. Also there can be multiple heat flow path through BEOL as elaborated in [START_REF] Esposito | Innovative SiGe HBT topologies with improved electrothermal behavior[END_REF]. Alternatively, the heat can flow to another device via the BEOL metal line. The fraction of the total heat flow through these BEOL paths depends on the actual structure, material properties and the nearby devices' temperature. This results in a different finite values of R m . From the perspective of device performance, the impact of BEOL design has also been investigated in [START_REF] Esposito | Influence of the BEOL metallization design on the overall performances of SiGe HBTs[END_REF]. While R m is relatively temperature independent, the variation of R s with temperature is significant due to the temperature-dependent thermal conductivity (κ(T )) of Si, which can be expressed as

κ(T ) = κ(T 0 ) T 0 T α . (1)
Here, κ(T 0 ) is the thermal conductivity of Si at the nominal temperature, T 0 . While extracting the thermal resistance parameters for accurate modeling purpose, the parameter α in the Si thermal conductivity model plays a crucial role; and sometimes its value deviates from the physically expected values due to the negligence of a finite R m component [START_REF] Gupta | Extraction of true finger temperature from measured data in multifinger bipolar transistors[END_REF], [START_REF] Huszka | Extraction of compact static thermal model parameters for SiGe HBTs[END_REF]. Also, a proper estimation of the R m component from characterization data can help us obtain a reliable model for R m depending on the number of BEOL metal layers and their design configurations. In specific situations, an accurate estimation of R m can help the modeling engineer to deembed the effects of huge BEOL metal stacks which may not affect the actual device placed in circuit. Indeed, our test structure (used in section IV) additionally has a very short BEOL thermal path to a heat sink (ground plane) allowing us to identify and correct the contribution of R m on devices in actual circuit layout configuration. Different approaches have been reported in the literature to separate the FEOL and BEOL components of thermal resistance [START_REF] Li | Experimental method to thermally deembed pads from r T H measurements[END_REF]- [START_REF] Balanethiram | Extraction of BEOL contributions for thermal resistance in SiGe HBTs[END_REF]. The work in [START_REF] Li | Experimental method to thermally deembed pads from r T H measurements[END_REF] extracted the FEOL

DTI BC junction T j R s T a R m P diss FEOL BEOL Fig. 1
. Schematic representation of a typical transistor cross section including BEOL. The base collector junction is the heat source and is represented as a current source in the equivalent circuit. The total resistance from heat source to ambient in the upward direction is represented by Rm and the same in the downward direction is represented by Rs in the equivalent circuit.

component by de-embedding the BEOL and the pads. The disadvantage of this method is that one requires test structures with multiple M1 bar widths to extract the resistance components. An analytical method to estimate the thermal resistance offered by the BEOL was proposed in [START_REF] Yeats | Inclusion of topside metal heat spreading in the determination of HBT temperatures by electrical and geometrical methods [GaAs devices[END_REF]. However, this method requires the knowledge of exact geometry and material properties of the BEOL stack. Another analytical technique to extract both the FEOL and BEOL components of the thermal resistance was proposed in [START_REF] Balanethiram | Extraction of BEOL contributions for thermal resistance in SiGe HBTs[END_REF]. The method extracts these components by choosing three points in the R th vs T j data. But finding R m and α from characterization data is essentially an inverse problem; hence not only difficult but also sensitive to the data in use. Therefore, it is not advisable to depend on only certain data points. It is found that application of a smoothing function on the overall available data set ensures an extraction of the thermal parameters that are less sensitive to the measurement noise. This paper presents a reliable methodology to simultaneously extract R m and α from characterized data sets. Section II elaborates the extraction methodology, while Section III demonstrates the validation of the proposed methodology using analytically generated synthetic data and TCAD simulation of different transistor structures. In section IV, we apply our technique on measured data obtained from different SiGe HBTs fabricated in state-of-the-art B11HFC process from Infineon and analyze the obtained results. Finally, we conclude in Section V.

II. EXTRACTION TECHNIQUE

Let G th (T a , P d ) be the overall thermal conductivity of the structure as a function of the ambient temperature (T a ) and electrical power dissipation (P d ). This data can be obtained from device characterization [START_REF] Pfost | A practical method to extract the thermal resistance for heterojunction bipolar transistors[END_REF], [START_REF] Balanethiram | Extracting the temperature dependence of thermal resistance from temperature-controlled DC measurements in SiGe HBTs[END_REF] and is often presented in the literature using the conventional R th versus T j plot [START_REF] Rieh | Measurement and modeling of thermal resistance of high speed SiGe heterojunction bipolar transistors[END_REF], [START_REF] Sahoo | A scalable model for temperature dependent thermal resistance of SiGe HBTs[END_REF]. R th = 1/G th is the overall thermal resistance of the structure and T j = T a +P d R th is the operating temperature of the device. However, for near-zero P d the overall thermal conductance of the structure is mostly a function of T a alone and is given as

G th (T a ) = G s (T a ) + G m (2)
where G m = 1/R m is the temperature-independent BEOL thermal conductance. On the other hand, G s = 1/R s is the thermal conductance corresponding to the FEOL substrate component, given as

G s (T a ) = G s0 T 0 T a α . (3)
G s0 = G s (T 0 ) is the FEOL thermal conductance at T 0 . Note that T a and T 0 (usually T 0 =300 K) are known quantities. The unknown quantities are G s0 , α and G m . One can begin by eliminating G m using G th evaluated at T a and T 0 from (2) as

G th (T 0 ) -G th (T a ) = G s0 1 - T 0 T a α . (4)
From the characterized data one can obtain the LHS of ( 4) as a function of T a and use an iterative curve fitting approach to solve for α and G s0 , simultaneously. Using (2) at T a = T 0 , one obtains the temperature-independent G m since G th (T 0 ) and G s (T 0 ) are already known. The robust nature of the extraction approach is due to the fact that the fitting function in the RHS of ( 4) is simple and involves no differentiation of the actual data. Besides, the estimated value of the parameter α is not influenced by R m ; hence, it reflects a physically consistent value corresponding to the FEOL substrate material.

III. VALIDATION OF PROPOSED METHODOLOGY

A. Analytically generated synthetic data

First we test the extraction approach on analytically generated synthetic thermal conductance data based on a temperature-dependent G s and a constant G m . One can use the well-known analytical model reported in [START_REF] Walkey | Equivalent circuit modeling of static substrate thermal coupling using VCVS representation[END_REF] and [START_REF] Balanethiram | Accurate modeling of thermal resistance for on-wafer SiGe HBTs using average thermal conductivity[END_REF] to generate T j -dependent G th data for different values of T a as

G th (T j ) = G s0 • T α 0 (1 -α) • T 1-α j -T 1-α a T j -T a + G m (5)
with, α, G s0 and, G m as model parameters.

The G th is evaluated for T j in the range (T a , 600 K) with α = 1.4, G s0 = 3.357 × 10 -4 W/K, and G m = 0.1 mW/K in (5). The process is repeated for multiple values of T a in steps of 30 K between 300 K and 420 K. Subsequently, the electrical power dissipation in the selected range of T j is calculated, for each T a , as P d (T j ) = G th (T j )•(T j -T a ). The resulting T j (P d ) for different T a values are shown with symbols in Fig. 2a. Similarly, G th is plotted against P d for multiple T a values and extrapolated to P d = 0 point to calculate the G th (T a ). Taking T 0 = 300 K, G th (T 0 ) -G th (T a ) is calculated and non linear regression using ( 4) is done on this data to extract the parameters.

Since there are two parameters in (4), the nonlinear regression is done by fixing one parameter. This makes the equation single-variable and easy to solve. Since the value of G m is less than G s0 , the value of G s0 will lie in the range G th (T 0 )/2 to G th (T 0 ). G s0 is varied in this range, and the corresponding α and mean squared error are recorded. The G s0 and corresponding α for which error is minimum are chosen as the extracted values. Fig. 2b shows the data points used for regression (symbols) and the fit to it (line). While generating the synthetic data, we have used G m = 0.1 mW/K which is exactly returned by the extraction scheme. The extracted values, tabulated in Fig. 2b are highly in agreement with the input values used to generate the data; hence, the model evaluated with these extracted values (lines in Fig. 2a) shows an excellent match with the synthetic T j (P d ) data sets (symbols).

B. TCAD generated data

Using commercially available Synopsys Sentaurus TCAD simulator [START_REF] Synopsys | Sentaurus: Sentaurus device user guide, release H-2013[END_REF], 3D thermal simulations of an HBT structure are carried out. The structure consisted of a substrate with a heat sink maintained at T a at the bottom surface and a BEOL structure with only one metal (Copper) layer. A heat sink is also placed on the top surface of BEOL, and a rectangular heat source of area 5 × 0.2 µm 2 is placed at the location of base-collector junction, 200 nm below the FEOL-BEOL interface. Here 200 nm is the sum of widths of the emitter and the base regions. Note that we approximate the actual nonuniform heat source with a rectangular uniform heat source. Since the proposed extraction method is independent of the uniformity of the heat source, we adopted this approximation, instead of a more rigorous approach followed in [START_REF] Alessandro | Advanced thermal simulation of SiGe: C HBTs including back-end-ofline[END_REF]. The heat source injects power at uniform power density, and all other outer surfaces except the heat sinks are assumed to be thermally insulating. The injected power is varied from 0 to 30 mW, and the temperature at the middle of the heat source is recorded as the junction temperature T j . A constant thermal conductivity is used for BEOL materials to ensure that R m is temperature independent. The temperature dependence of thermal conductivity in Silicon is included using [START_REF] Raghunathan | Bias-and temperature-dependent accumulated stress modeling of mixed-mode damage in SiGe HBTs[END_REF]. The simulations are repeated for three different alpha values: 1.3, 1.4 and 1.5, keeping κ(T 0 ) = 1.702 W/(cm-K) [START_REF] Schröter | Compact Hierarchical Bipolar Transistor Modeling with HICUM[END_REF].

Fig. 3a shows the dissipated power-dependent junction temperature at different ambient temperatures obtained from TCAD simulations with α = 1.4. Thermal conductance G th = P d /(T j -T a ) is calculated from the TCAD data and )) and the structure are the same, the extracted G s0 and R m values remain the same across these devices, as expected. The correlation between symbols (TCAD) and lines (model evaluated with extracted parameters) in Fig. 3a re-confirms that the values extracted are accurate.

C. Generating T a dependent data

The method proposed here requires peak temperature data as a function of dissipated power at many different ambient temperatures. In cases where T j (P d ) for single T a is available, one can generate the required data. This subsection explains how the T j (P d ) data for multiple ambient temperatures can be generated from the data measured at a single T a . The total dissipated power for a given T j and T a can be obtained by multiplying (5) with T j -T a on both sides and following simple algebraic steps one obtains

P d (T j , T a ) = T j -T a R m + G s0 T 1-α j -T 1-α a (1 -α)T -α 0 = T j -T ′ a R m + G s0 T 1-α j -T ′ a 1-α (1 -α)T -α 0 + T ′ a -T a R m + G s0 T ′ a 1-α -T 1-α a (1 -α)T -α 0 = P d (T j , T ′ a ) + P d (T j = T ′ a , T a ) (a) (b) (c)
Since ∆P d is independent of T j , formulation [START_REF] Esposito | Thermal penetration depth analysis and impact of the BEOL metals on the thermal impedance of SiGe HBTs[END_REF] signifies that the dissipated power dependent junction temperature data for another ambient temperature T ′ a can be obtained just by shifting the original curve along the power axis by ∆P d . From the T j (P d ) at a given T a , one can, therefore, obtain other T j (P d ) curves for multiple different T ′ a values. Using this methodology, we have regenerated all the curves shown in Figs. 2a and 3a by shifting the lowest curve with T a = 300 K and obtained perfect agreements with the corresponding TCAD generated curves. Thus we have verified the methodology which will be useful in the following section while applying our extraction scheme on actual measured data.

IV. EXTRACTION FROM MEASURED DATA

The state-of-the-art SiGe BiCMOS HBT technology from Infineon (B11HFC) having f T = 250 GHz, f max = 370 GHz and six levels of metallization reported in [START_REF] Esposito | Electro-thermal characterization, TCAD simulations and compact modeling of advanced SiGe HBTs at device and circuit level[END_REF] is used to demonstrate our extraction technique. The experimental data presented in this paper corresponds to a maximum of four metal layers (till M4). The transistors are isolated using deep trench isolations (DTI). The test structures had the same FEOL with an emitter dimension of 0.35 × 5 µm 2 but different metallization schemes in the BEOL part. In the first type, we have dummy metal layers of the same width and length directly above the emitter contact but not above the base and collector contacts as shown in Fig. 4a. This test structure will be called E4narr. The number denotes the last layer of the metal stack and E indicates that the metal dummies are on Emitter contact. In the second type, we have two test structures, namely E2wide and E4wide. Here the difference is that the width of the metal dummies on the emitter contact increases progressively as shown in Fig. 4b. Note that the total volume of metal on top of the emitter contact of E4wide structure is more than E4narr. Lastly, we have two other test structures named B2 and B3, where the metal dummies are placed directly on top of the base contacts but not on the emitter and collector contacts as shown in Fig. 4c. All the test structures considered in the study are listed in Table II. In section III, for both synthetic and TCAD data, the BEOL thermal conductivity was kept constant. However, the actual BEOL thermal resistance varies with temperature due to the temperature-dependent thermal conductivity of the BEOL components. Since this dependence is significantly weak compared to that of FEOL, one can assume R m to be independent of temperature especially at temperatures close to T a . Therefore, for the R m extraction, we are going to use only those data points which are close to T a . Since T j (P d ) data for a single T a is available for each kind of devices discussed above, required data for multiple T a values are obtained using the shifting method elaborated in section III-C. Fig. 5a compares the R th (T j ) data corresponding to the structure B3 obtained from measurement (symbols) and that obtained from (5) (dashed line) using the extracted parameters. Note that the extraction is carried out using the data points (T j < 313 K) close to T a . The non-linear regression plot for the same is shown in Fig. 5b. The model (dashed line) and the measurements are in excellent agreement in the lowtemperature range as shown in Fig. 5a. [START_REF] Esposito | Electro-thermal characterization, TCAD simulations and compact modeling of advanced SiGe HBTs at device and circuit level[END_REF]. After extracting α and G s0 and R m = R m0 at low temperature, the junction temperature (T j ) dependent R m can be estimated from the measured R th (T j ) data as

Test α G s0 R m0 m 1 m 2 Structure [mW/K] [K/W] [/K] [/K 2 ]
1 R m (T j ) = 1 R th (T j ) - 1 R s (T j ) (7)
where

R s (T j ) = (1 -α) • T -α 0 G s0 • T j -T a T 1-α j -T 1-α a (8)
is the FEOL component of thermal resistance evaluated using the parameters extracted from the low temperature data. This temperature dependence of R m can be captured with simple polynomial functions. Symbols in Fig. 5c shows the normalized variation of R m (T j ) calculated using [START_REF] Esposito | Innovative SiGe HBT topologies with improved electrothermal behavior[END_REF]. The variation of the extracted R m (T j ) from the measured data is not negligible unlike the ones enforced in the synthetic as well as the TCAD simulation in section III. Possible reasons for such a variation can be as follows. According to our modeling framework, all the regions above the heat source (i.e., the base-collector junction) including the thin SiGe base as well as the mono and poly-emitter portions are considered within the upward heat transfer path contributing to the overall R m along the actual BEOL stacks. The thermal conductivity of these base and emitter regions decrease significantly with temperature and contributes to the rise in R m with temperature. Along with that the thermal conductivity of metals within the BEOL slightly reduces with temperature, causing an increase (although minor) in R m with temperature. Note that in TCAD simulation (presented in section III), we have placed the heat source at the emitter finger and considered constant thermal conductivity for the metal layers; hence obtained a constant R m .

It is clear that the FEOL component can be evaluated from the extracted parameters using [START_REF] Esposito | Influence of the BEOL metallization design on the overall performances of SiGe HBTs[END_REF] and subsequently R m (T j ) can be estimated using [START_REF] Esposito | Innovative SiGe HBT topologies with improved electrothermal behavior[END_REF]. Finally, assuming that the temperature-dependence of R m can be captured with a parabolic function given as

R m (T j ) = R m0 1 + m 1 (T j -T 0 ) + m 2 (T j -T 0 ) 2 , (9)
the parameters m 1 and m 2 are extracted using simple polynomial fit. Fig. 5c presents a simple modeling results for R m (T j ) (normalized w.r.t. R m0 ) obtained by using [START_REF] Gupta | Extraction of true finger temperature from measured data in multifinger bipolar transistors[END_REF]. After including this temperature dependent R m model, the complete model (solid line) and the measured data (symbols) show excellent agreements as demonstrated in Fig. 5a.

Table II for all the devices, α and G s0 values are in good agreement across the devices. Even the parameters of temperature dependent R m are also in the same range. The variation in R m0 is expected because of the deliberate BEOL change made in the test structures. Since stacking more dummy levels reduces the BEOL thermal resistance [START_REF] Esposito | Electro-thermal characterization, TCAD simulations and compact modeling of advanced SiGe HBTs at device and circuit level[END_REF], B3 has smaller R m0 than B2. Use of wider metal layers also reduces R m0 as it increases the volume of metal present in the BEOL structure. Therefore, we expect E4wide to have the least R m0 due to four levels of metal dummies above the emitter contact. However, we see that the extracted R m0 value of E4wide is not consistent. Such an error is mainly due to the sensitivity of the parameters obtained from a nonlinear extraction method. Eventually, we performed a sensitivity study on the extracted R m0 values by fixing the FEOL thermal resistance across the test structures. We calculate different R m0 values for a given test structure with different G s0 values listed in Table II following the relation, 1/R m0 = G th (T 0 ) -G s0 . The resulting R m0 values for all the test structures are represented as a bar chart in Fig. 6. The values of R m0 obtained for the test structures with an average G s0 are also plotted (solid line). Since these structures have the same FEOL, the solid line is a better representation of the actual R m0 variation. Note that these values are in alignment with our expectation, For example, we observe that using identical G s0 for all the test structures yields the least R m0 for the E4wide structure. From the various combinations considered in Fig. 6, the maximum deviation of R m0 from the corresponding average value is found to be less than 5.4%. For such an inverse problem, this range of error is fairly acceptable.

V. CONCLUSIONS

A robust method to extract the parameters related to the FEOL and BEOL components of thermal resistances in SiGe HBTs was presented. The method allows extraction of these components and their temperature dependent parameters. The method, when tested on synthetic or TCAD generated data, gave excellent results. We also applied the method on measured data obtained from devices with the same FEOL but different BEOL stack. The returned values were similar for FEOL and showed an expected trend for BEOL components across the devices confirming the proposed method's accuracy.

Fig. 2 .

 2 Fig. 2. (a) Dissipated power dependent junction temperature at different ambient temperatures: comparison of analytically generated synthetic data and model using the extracted parameters. (b) Ambient temperature dependent G th (T 0 )-G th (Ta): analytical data and fitting by RHS of (4) to extract G s0 and α.

Fig. 3 .

 3 Fig. 3. (a) Dissipated power dependent junction temperature at different ambient temperatures: comparison of 3D TCAD thermal simulation (with α = 1.4) and model using the extracted parameters. (b) Ambient temperature dependent G th (T 0 ) -G th (Ta): TCAD data and fitting by RHS of (4) to extract G s0 and α.

Fig. 4 .

 4 Fig. 4. The test structures considered (a) E4narr -test structure with 4 dummy metal layers of same width and length placed above emitter contact. (b) E4wide -test structure with 4 dummy metal layers placed above emitter contact, whose width increases from first layer to fourth layer. (c) B3-3 level structure where the stacking of dummy metal layers are done directly above the base contact. The images are taken from [23].

6 Fig. 5 .

 65 Fig. 5. (a) Junction temperature dependent thermal resistance: comparison of measured data (symbols) and model (lines) using the extracted parameters. Dashed line represent the constant R m0 model and solid line represent the model with temperature-dependent Rm. (b) Ambient temperature dependent G th (T 0 )-G th (Ta): measured data and fitting by RHS of (4) to extract G s0 and α. (c) Rm estimated using (7) with extracted R m0 .

TABLE II PARAMETERS

 II EXTRACTED FROM MEASURED DATA

  presents the result of extraction for all the devices considered in this study. As expected, since the FEOL is same R m0 for the test structures calculated using different G s0 values listed in TableII. The solid line represents R m0 obtained for each structure with an average G s0 .

		8,000
			using Gs0 of:	E4wide	B2	B3	E2wide	E4narr
		7,500					Using Avg. Gs0
	Rm0 [K/W]	7,000
		6,500
		6,000
		B2	E2wide	E4narr		B3	E4wide
			Test structures
	Fig. 6.

This work was supported in part by ISRO project ELE/17-18/176/ISRO/ANJA, in part by DST, India, under Project EMR/2016/004726 The authors would like to thank Infineon for providing Silicon hardware. The test-structure designs used in this work were a part of the joint European research project DOTSEVEN.


Quentin Vandoolaeghe

Claire Lopez-Piffet

Valérie Bouchart

Raphaël Delépée

The challenge to produce robust data in UHPLC-HRMS/MS based metabolomics for molecular epidemiology studies

Quentin Vandoolaeghe 1,2,3 , Claire Lopez-Piffet 1,3,4 , Valérie Bouchart 1,4,5 , Raphaël Delépée 1,2,3,4 1 Institut National de la Santé et de la Recherche Médicale (INSERM), Unité Mixte de Recherche (UMR) 1086 ANTICIPE, Caen, France; 2 Centre de Lutte Contre le Cancer François Baclesse, Caen, France; 3 Université de Caen Normandie, Caen, France; 4 PRISMM platform, UNICAEN, CLCC François Baclesse, Caen, France; 5 LABÉO, Pôle recherche développement et innovation, Caen, France *quentin.vandoolaeghe@unicaen.fr

The challenge to produce robust data in UHPLC-HRMS/MS based metabolomics for molecular epidemiology studies Utilisation of biobanks samples associated with exposure description and tasks listing for a metabolomic analysis

Context : Agricultural activity and cancer 3 ❑❑

 3 Study cancer biological plausibility, early stages effects, lookout for predictives markers of cancer risks ❑ Holistic approach of cancers occurrence in farmers population Metabolomic Ion annotation tool optimisation ❑ Biobank samples analysis Alseekh S et al. Mass spectrometry-based metabolomics: a guide for annotation, quantification and best reporting practices. Nat Methods. 2021 General Workflow ❖ Experimental design ❖ Sample preparation ❖ Metabolites separation ❖ Ionisation ❖ Detection ❖ Data treatment ❖ Statistical analysis ❖ Metabolites identification ❖ Biological interpretation Many steps to be optimized ❑ Study on blood, hair or saliva ❑ Use of new statistical tools 5


Zawilski B. M., Granouillac F., Brut A., Tallec T., Gibrin H. -Transmission of the data to the CP OBJECTIVES: Build a robust, automatized, "smart" system to optimize data retrieval and processing.

bartosz.zawilski@cnrs.fr

Context & Objectives

loggers' organization

Conclusion Data transmission to CP

➢ Each logger (mainly Campbell-Sci. loggers, UK but also SmartFlux2 from Licor, USA) is connected to a router via internet cables. Other devices such as a camera or local PC are also present on the stations and connected via RJ-45 cables. ➢ Some loggers have its RJ-45 interface others do not have it and the communication with them is insured via their serial port RS-232. To connect it through the internet, a serial server module is used (USR DR-301, China). An internet WiFi relay allows to hard reset each logger via internet commands. ➢ All RJ-45 cables are plugged into the router of the station Lock 500 from the Tosibox with an embedded 4G LTE modem (no subscription, Finland). This router opens a VPN tunnel through the internet.

➢ Some data files, but not all, have to be uploaded every day, around midnight, to the CP.

The uploaded data files have to be checked first to correct missing data lines, if any, and to merge fragmented files if any. Indeed, in some cases, it is possible for loggers to produce several fragmented daily files instead of one complete file and to miss a few data lines. In case of incomplete file transmission to the CP, it will be rejected and flagged as comporting errors or fatal errors. To automatize the process a reliable yet simple-to-use software was developed under Labview from National Instrument, USA. This software is distributed freely on CP for the ICOS community (https://fileshare.icoscp.eu/apps/files/?dir=/ICOS%20RI%20groupdisk/Ecosystem/MSA/Codes%20shared/Files %20Uploader&fileid=1852945) seek for every fragmented file and merge it into one file then scan it determining the right sampling rate and eventually missing data lines. In case of missing lines inserts data lines with the correct timestamp along with agreed virtual data (-9999 for us). Once the file is corrected, it is uploaded on CP. All steps are consigned in a log file to allow accurate tracking.

  Centre d'Etudes Spatiales de la BIOSphère 18 avenue Edouard Belin, 31400 Toulouse, FRANCE ENTIRELY AUTOMATIZED, SECURED DATA RETRIEVAL, PROCESSING AND UPLOADING SYSTEM ➢ In the context of global climate warming, several monitoring scientific stations networks were developed around the world. In Europe, we have an Integrated Carbon Observation System (ICOS) network to which belongs our two cropland stations FR-Lam (class 1) and FR-Aur (class 3). ➢ ICOS stations use an important number of sensors with high sampling rate data recording, especially Eddy Covariance (EC) carbon dioxide measurements setups (10Hz sampling rate). All these data have to be logged and transmitted every day to the ICOS server on Carbon Portal (CP). ➢ Three distinct but interdependent problems arise: -Internal station organization of the loggers. -Collection and transmission of the logged data to the laboratory.

Acknowledgments

Data acquisition at FR-Lam and FR-Aur were mainly funded by the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (CNRS-INSU) through the ICOS and OSR SW observatories. Facilities and staff were also funded and supported by the University Toulouse III -Paul Sabatier, the CNES and IRD (Institut de Recherche pour le Développement). We would also to acknowledge the ETC team and especially Dario Papale and Simone Sabbatini for helping us with data processing and for our Carbon Upload Software sharing on the CP server.

Data transmission from the stations to the laboratory ➢ Tosibox system insure a fast reliable and secured; 256 bits encrypted data transmission via a Virtual Private Network (VPN). In the lab, a PC, a so-called VPN PC, with a matched master key from Tosibox (a USB plug) establishes bidirectional communication between this PC and every device plugged into the router. Official software from Campbell-Sci.; LoggerNet, is retrieving data from Campbell loggers stocking it on the VPN PC. For other data, such as EC data from Smartflux2, an FTP server on VPN PC is embedded and receive files from SmartFlux2 or cameras. All IP addresses can be static with a generic usage SIM card.

➢ Measurements go hand in hand with data collection, verification, and transfer-storage. For this purpose, a robust system should be built on sites and in the lab. Nowadays, solutions derivate from omnipresent IOT technology allow to make an entirely automatized "bulletproof", cost-effective setups.


Abdullah F Alabiad
email: abdullah.alabiad2@femto-st.fr

Blandine E Guichardaz

Fadi I Baida

Control of the light extraction from a photonic crystal nanocavity by coupling with a nanoparticle

Keywords: Photonic Crystal, Light Extraction Enhancement, Directivity, Radiation Angular Spectrum

In this paper, we present a numerical study (using Finite Difference Times Domain (FDTD) method) of the coupling between two resonant structures with different quality factors of their resonances. The first is a Photonic Crystal Cavity (PCC) while the second is a dielectric Nano-Particle (NP). Both structures exhibit optical responses at the same wavelength, which allows a critical coupling between them. The objective of the study is to exploit this coupling to control the light extraction properties of the more resonant structure (the PCC). The challenge is then: first, to model such a near-field coupling through numerical method that can handle the entire structure (PCC+NA) and second, to get physical insight about the interaction allowing a good control of the light extraction. We numerically demonstrate that this coupling strongly depends on the spatial position of the NP relative to the PCC. We also show that the direction of the emission (into the substrate or superstrate) and its intensity can be adjusted by the NP position relative to the PCC. Quantitatively, the critical coupling corresponds to an NP-PCC distance of 120 nm leading to an increase by factor of 3 of the light radiated towards the far field. This is accompanied by a red shift of 0.53 nm of the cavity resonance. This study opens the way to the development of a new kind of sources with both "controllable" wavelength resonance and radiation pattern.

Introduction

The search for techniques to improve the efficiency of light extraction and to manipulate its directivity from sources such that light-emitting diodes (LEDs), Photonic Crystals (PC) or optical nano-antennas (NA), has attracted the attention of researchers in recent years. Some have used a PC as a reflecting layer to improve by a factor of 2.6 [START_REF] Bergenek | Directional light extraction from thin-film resonant cavity light-emitting diodes with a photonic crystal[END_REF] the directivity or the efficiency by 47% [START_REF] Bergenek | Enhanced light extraction efficiency from algainp thin-film light-emitting diodes with photonic crystals[END_REF] and, even more, by 70% [START_REF] Boroditsky | Light extraction from optically pumped light-emitting diode by thin-slab photonic crystals[END_REF] of the light extraction from thin-film LEDs. Others have experimentally demonstrated that they can increase the light extraction efficiency from a PC up to 70% [START_REF] Shanhui Fan | High extraction efficiency of spontaneous emission from slabs of photonic crystals[END_REF][START_REF] Boroditsky | Spontaneous emission extraction and purcell enhancement from thin-film 2-d photonic crystals[END_REF] and up to 8 times higher than non-patterned structure with high directivity of light emitted [START_REF] Ndiaye | Enhanced light-extraction efficiency and emission directivity in compact photonic-crystal based algainp thin-films for color conversion applications[END_REF] by changing its geometrical parameters (lattice constant, pattern shape) . On the other hand, some researchers have used optical nano-antennas to control the radiation pattern of given emitters [START_REF] Svetlana V Boriskina | Q factor and emission pattern control of the wg modes in notched microdisk resonators[END_REF][START_REF] Taminiau | Enhanced directional excitation and emission of single emitters by a nano-optical yagi-uda antenna[END_REF][START_REF] Cp Dettmann | Directional emission from an optical microdisk resonator with a point scatterer[END_REF][START_REF] Devilez | Compact metallo-dielectric optical antenna for ultra directional and enhanced radiative emission[END_REF][START_REF] Rolly | Boosting the directivity of optical antennas with magnetic and electric dipolar resonant particles[END_REF][START_REF] Alexander E Krasnok | Superdirective dielectric nanoantennas[END_REF][START_REF] Rose | Control of radiative processes using tunable plasmonic nanopatch antennas[END_REF][START_REF] Rocco | Controlling the directivity of all-dielectric nanoantennas excited by integrated quantum emitters[END_REF]. One of them [START_REF] Taminiau | Enhanced directional excitation and emission of single emitters by a nano-optical yagi-uda antenna[END_REF] shows a maximum directivity emission up to 3.8 times higher with a rotation change of 90 • compared to the free emitter. In another approach, an enhancement factor up to 260 times for the fluerescence efficiency of a PC was experimentally obtained by coupling it with another resonator such as a nanoparticle [START_REF] Cheng-An Tao | Coupling of nanoparticle plasmons with colloidal photonic crystals as a new strategy to efficiently enhance fluorescence[END_REF][START_REF] Eftekhari | The effect of fluorophore incorporation on fluorescence enhancement in colloidal photonic crystals[END_REF][START_REF] Wu | Manipulating luminescence of light emitters by photonic crystals[END_REF]. In all these studies, the configurations are passive because their geometry is fixed and does not allow to modify the coupling properties. Here we bring a first approach to lift this lock by considering a configuration where one of the elements can move relative to the other. Achieving experimentally this flexibility in the configuration will be discussed in the conclusion. For this purpose, we numerically determine the angular radiation spectrum of light extracted from a Photonic Crystal Cavity (PCC) coupled to an Electric Dipolar (ED) resonant Nano-Particle (NP) as a function of its spatial position and demonstrate the possible control of the light emission properties of this combination. First, we present the properties of the PCC which has already been studied theoretically and experimentally by our team [START_REF] Fadi | Double-way spectral tunability for the control of optical nanocavity resonance[END_REF][START_REF] Grosjean | Extraordinary blueshift of a photonic crystal nanocavity by reducing its mode volume with an opaque microtip[END_REF]. Second, we numerically design a dielectric NP using Mie theory so that its electric dipolar resonance matches the resonance wavelength of the PCC. Third, we calculate, utilizing the FDTD results, the angular spectrum of radiation of the whole system (PCC+NP) both in the substrate and the superstrate with respect to the position of the NP. We demonstrate how to control the light emission direction and energy by moving in the 3 directions the NP with respect to the PCC.

Numerical Studies

In the following, we will first study the optical properties of each resonators (PCC, NP) individually and, second, their coupling. Due to the fact that the cavity mode of the considered PCC occurs for a well given value of the wavelength, we design the NP so it presents an electric dipole (ED) resonance at the same wavelength by modifying its refractive index. Let mention that the calculations are done within the Mie theory in the case of self-suspended spherical NP while an adapted TFSF (Total Field Scattered Field) algorithm [START_REF] Taflove | Computational Electromagnetics: The Finite-Difference Time-Domain Method[END_REF] is used in the case of NP in the vicinity of dielectric substrate. All these codes are home-made codes as the one used to simulate the NP to PCC coupling.

Study of the PCC

The considered PCC was proposed since 2003 by Monat et al. [START_REF] Monat | Two-dimensional hexagonal-shaped microcavities formed in a two-dimensional photonic crystal on an inp membrane[END_REF] to exhibit a cavity mode in the NIR spectral range. It consists on a 2D hexagonal photonic crystal composed of air hole cylinders of radius r = 105 nm with a lattice constant a = 420 nm. The holes are engraved into an InP layer of 300 nm-thickness deposited on a silica layer of 1 µm-thickness on top of silicone substrate. The optical cavity is composed of seven missed and aligned air holes (then called CL7) as shown in figure 1(a). Note that the two air holes at the extremities of the cavity are offset 75 nm towards the outside of the cavity in order to enhance the quality factor of the cavity mode resonance keeping it inside the photonic bandgap of the PC dispersion diagram. The near-field spectrum calculated 15 nm above the PCC center is shown in figure 1(b). The cavity mode occurs at the wavelength λ mode 0 = 1602.92 nm with a quality factor of Q 0 = 4035 as experimentally measured in ref. [START_REF] Grosjean | Extraordinary blueshift of a photonic crystal nanocavity by reducing its mode volume with an opaque microtip[END_REF]. Figures 1(c-e) present the amplitude distributions of the 3 components of the electric field E x , E y and E z respectively while the figure 1(f) presents the total intensity distribution of the electric field. First, we can see that, at the cavity center, the E y -component is maximum while the E x -and E z -components are zero.

Second, it is clear that the electric field components are spatially distributed differently along the major axis of the cavity i.e. a possible coupling to a NP should take into account the position of the latter as we will see later.

The magnetic field components are also given in figure 1(g-i) as well as its intensity in figure 1(j). As expected and because the cavity mode is an in-plane mode, H z is predominant but its value vanishes at the cavity centre similarly to H y while H x is maximum. It may be mentioned here that, in the vertical direction (Oz), the cavity mode is located at the center of the InP membrane.

Study of the NP

As mentioned above, the considered spherical NP is designed so that it exhibits an ED resonance close to the wavelength resonance of the PCC mode. Thus, knowing the operation wavelength to be λ mode 0 = 1602.92 nm, two solutions can be envisaged : (i) Fixing the radius R of the NP and changing its refractive index n p , or (ii) Fixing the index and varying its radius. Nonetheless, if we want the coupling of the NP with the PCC to be local, we should consider the size of the NP smaller than the size of an electric or magnetic field lobe of the PCC mode (see figure 1(f)). This limits the choice to solution (i). From the distributions shown in figures 1(c-f), we estimate the size of this lobe to be about 440 nm (see the scale bar in figure 1(c)). Therefore, we fix the NP radius to be R = 220 nm. Using the Mie theory, we have calculated the spectral position of the ED resonance as a function of the refractive index of the NP (n p ). Figure 2(a) presents the obtained results providing the values of n ED p = 4.85 for the ED to be located at the PCC resonance wavelength. Note that this value of the refractive index does not necessarily correspond to a naturally existing material (even if Re[ε GaAs (λ = 450 nm)] = 4.8427) but it is needed to allow us understanding the physical effects related to the studied phenomenon. Figure 2(b) shows the effective scattering cross section (σ S ) spectrum calculated using the Mie theory [START_REF] Craig | Absorption and scattering of light by small particles[END_REF] (dashed line) in comparison with the one obtained by TFSF-FDTD [START_REF] Taflove | Computational Electromagnetics: The Finite-Difference Time-Domain Method[END_REF] (solid line) algorithm that is described in the Appendix A. Indeed, this comparison helps to validate our FDTD codes which will be used in all the following simulations, especially for the study of the PCC-to-NP coupling. From figure 2(b), we see a fairly good agreement between the two results (see Appendix A for more information). Figures 2(c,d) show the distribution of the normalized electric field amplitude recorded respectively in the xz and yz planes passing through the center of the NP excited at the ED resonance wavelength. (c,d) are the normalized electric field amplitude distributions recorded in the xz and yz planes respectively at the ED excitation. Note that the incident plane wave is propagating along the z > 0-axis and is linearly polarized along the x-direction.

Study of the NP-dielectric substrate coupling

The optical response of a resonator being strongly dependent on its environment [START_REF] Elie M Atie | High optical resonance sensitivity to its environment of a fibered bowtie nano-aperture antenna[END_REF][START_REF] Elie M Atie | Remote optical sensing on the nanometer scale with a bowtie aperture nano-antenna on a fiber tip of scanning near-field optical microscopy[END_REF], it is thus necessary to determine the influence of the presence of the InP substrate (without the cavity) as shown in figure 3(a) placed at a distance D z from the NP on their resonance properties. So, figure 3(b) shows the values of σ S and the spectral position (λ ED ) of the ED excitation as a function of D z . These results were obtained using an adapted TFSF-FDTD code taken into account the presence of the substrate. The illumination consists of a plane wave propagating from the substrate perpendicularly to the air-InP interface (see schematic of the system in figure 3(a)). As expected, a slightly small quasi periodic modification of σ S along the z-direction occurs due to the presence of the interference pattern between the incident field and the small one reflected by the NP itself. This leads to a small modulation of the σ S with a peculiar behavior in the near field where the coupling involves evanescent waves that are diffracted by the NP. We can notice that the oscillations appearing in the two curves are not synchronous even far from the substrate. This can be attributed to a geometrical effect (spherical NP) of the resonance contrary to that of a point dipole. Figures 3(c,d) present the normalized electric field amplitude distributions in the two vertical planes (xz and yz) that cross the center of the NP, when the latter is at D z = 45 nm above the substrate for a ED excitation (λ = 1624 nm). Let notice that, at this distance, the maximum of the electric field amplitude is enhanced compared to the case of the self-suspended NP (see figures 2(c,d)). This increase is accompanied by a large shift in the spectral position of the ED mode (up to 22 nm) when approaching the NP to the substrate, revealing the possibility of a tunable control of the PCC resonance wavelength for small PCC to NP distance as it will be shown in the following.

Study of the PCC-NP coupling and discussion

Let us describe the approach considered in order to achieve the main objective of this study that consists on demonstrating how the presence of NP can modify the pattern of light scattered by the cavity. Obviously, we have almost infinite possibilities for the spatial position of the NP with respect to the PCC. Consequently, we have made the choice to start from the center of the PCC, where the electric field intensity is maximum, and to move the NP along three perpendicular axis parallel to the z, x and y-directions respectively. Therefore, we will divide our study and discussion into two parts: 2.4.1. Part A: Near field, transmission and reflection We begin the study by determine the spectral behavior of the whole system resonance when the NP moves along the z-direction due to the fact that the coupling efficiency is more sensitive to the NP-to-PCC distance. Figures 4(a-c) show the evolution of the normalized spectra for: (a) the near-field "NF " (the normalized electric intensity recorded at 15 nm above the center of the PCC), (b) the transmission "Tr " (the normalized transmitted electric energy of the wave propagating to the far-field region along the z > 0), and (c) the reflection "Refl " (the normalized electric intensity of the wave propagating along the z < 0) as a function of the vertical position (D z ) of the NP. These quantities are calculated according to the equations below:

N F = || - → E N F d (0, 0, z d )|| 2 (1)
T r = +∞ -∞ || - → E t (x, y, z t )|| 2 dx.dy +∞ -∞ || - → E inc (x, y, z t )|| 2 dx.dy (2)
Ref l = +∞ -∞ || - → E r (x, y, z r )|| 2 dx.dy +∞ -∞ || - → E inc (x, y, z r )|| 2 dx.dy (3)
Where :

|| - → V || is the modulus (norm) of the vector - → V , - → E N F d (x d , y d , z d
) is the electric field vector recorded by a detector "d" placed at a distance "z d = 15 nm" above the center of the PCC (x d = 0, y d = 0), -→ E t (x, y, z t ) and -→ E r (x, y, z r ) are the diffracted zero-order of the electric field in the superstrate (z t > 0) and the substrate (z r < 0) respectively.

-→ E inc is the same electric field vector calculated without the presence of the PC nor the NP.

As expected, a red shift of the resonance wavelength occurs for small distances accompanied by an enhancement of the light extraction in both substrate (reflection) and superstrate (transmission) regions. This result clearly demonstrates the possibility to control the light extraction from the PCC by changing the NP position. Nonetheless, contrary to what was obtained for the NP alone, the maximum spectral shift of the PCC-NP resonance does not exceed 2.5 nm instead of 22 nm for the NP-substrate one. This is the signature of the great dominance of the PCC resonance even if the coupling is critical; predictable through the great value of its quality factor compared to that of the ED resonance of the NP alone. To go further, we have also studied the influence of the lateral position of the NP on the PCC resonance behavior. Figures 4(d-f) present the corresponding results when the NP moves along the x-direction parallel to the PCC major axis starting from its center while figures 4(g-i) correspond to a y-displacement of the NP along the perpendicular direction (PCC minor axis). In both cases, the NP vertical position was fixed to the D z = 45 nm corresponding to a detector placed in air at 15 nm between the PCC and the NP. As can be seen, a strong modulation of the transmitted and the reflected intensity is obtained as a function of the position of the NP depending on whether it is above a maximum or minimum of the electric intensity of the cavity mode (see figure 1(f)). Moreover, we see on figures 4(e,f) that the minima of the transmitted and reflected intensities occur simultaneously and correspond to a total extinction for D x = 0.5 a = 210 nm and D x = 1.5 a = 630 nm. The same phenomenon also occurs when the NP moves along the y direction for D y = 0.5 a = 210 nm (see figures 4(h,i)). This implies that the energy is no longer scattered in the direction perpendicular to the structure otherwise the interpretation will be very difficult (along the z-axis) but in other directions which remain to be determined. Moreover, the same figures show that the maxima of transmission and reflection also occur simultaneously meaning that the energy flow is probably mainly directed along the vertical direction.

Part B: Radiation pattern

In order to provide a clear and complete explanation on the modification of the emission, we have calculated the whole angular spectrum (see Appendix B) of the light emitted by the PCC alone as well as in the presence of the NP at different positions along both x and y directions in a plane located at D z = 45 nm. In this part, we will impose that we are in the continuous wave (CW) regime which corresponds to the cavity mode frequency (i.e. without coupling) otherwise the interpretation will be very difficult. In order to fully understand the light emission, we present the 4π sr radiation diagram of the PCC in figure 5 for four cases: in figure 5(a) we see that the angular spectrum of the energy scattered from the isolated PCC presents two lobes located by azimuth angle of ϕ = ±30 • and directed outward from the vertical axis (z-axis) along almost θ = 31 • , while three lobes appear in the substrate side where two lobes oriented outward from the vertical axis along θ = 40.7 • . A third lobe exists and is directed vertically (θ = 0 • ). When the NP is placed at position "A(0,0,45)" (see figure 5(b)), the scattered energy by the whole system (PCC+NP) is mainly directed in the vertical axis in transmission (air) and reflection (substrate) and its efficiency is about 1.79 times higher than that of the isolated PCC. Moving the NP to "B(630,0,45)" and "C(0,210,45)" positions leads to direct the scattered energy in special directions in the transmission and reflection regions, as shown in figure 5(c,d) with a slight enhancement of its efficiency to 1.1 times that of the isolated PCC. To have an overview, we made three videos that show the evolution of the angular spectra of the PCC-NP system according to the position of the NP. The movie "Movie1" shows this evolution when the NP moves vertically away from the center of PCC. The movies "Movie2" and "Movie3" show the evolution of the angular radiation spectra of the system when the NP is moved along the major axis of the PCC (x-axis) and along the minor axis of PCC (y-axis), respectively.

Through figures 4(a-c) and Movie 1, we can see that the PCC-NP coupling, where the NP moves vertically away from the cavity along the z axis, establishes two distinct cases of light extraction : for D z ≤ 210 nm, the majority of the light is extracted vertically, whereas for larger distance, the coupling becomes increasingly weaker and thus the light scatters further away from the z-axis. On the other hand, the lateral displacement of the NP leads us to the different interpretations. In the case of a displacement along the x-axis (movie 2), the light extraction behaves in 3 different aspects: When the NP is positioned above an intensity for the isolated PCC, in (c) for the PCC-NP system when the NP placed 45 nm above the PCC center at position "A", in (d) for the PCC-NP system when the NP at position "B" and in (e) for the PCC-NP system when the NP displaced along the y-axis and placed at position "C". The 2D circular maps (see last two lines) give the angular spectra calculated in air (z>0) and in the substrate (z<0) for the four corresponding cases (isolated PCC, NP at A, at B and at C respectively). maximum (see figure 1(f)), we see that the behavior of the extracted light is almost similar to the case of an NP placed above the center of the cavity (i.e., the majority of the light propagates vertically along the z-axis). When placed above an intensity minimum, the extracted light scatters in a quasi-isotropic shape with an extinction in the vertical direction leading to a "donut" form. In the third case presented by movie 3 (i.e. NP moving along the y-axis), the emitted light behaves differently from what we got in the previous case. The variation of the direction of the extracted light is a bit irregular because the NP will pass over two different media (air and InP) when it moves. This irregularity is illustrated in figures 4(g-i) where there is only one position of the NP (D y = 0.5 a = 210 nm) that causes an extinction of the light emitted in the vertical direction. This position appears clearly in figure 1(f).

Conclusion

In summary, we present the effect of the coupling between a dielectric NP and an optical nanocavity exhibiting a same resonance wavelength. This results obtained demonstrate that this kind of coupling could offer a possible way to shape the angular radiation diagram of the whole structure. The weight and direction of the extracted light from the cavity can be tuned according to the NP position with respect to the cavity allowing a directive and selective emission into the substrate and/or the superstrate. We studied different configurations and showed that the coupling depends greatly on the near electromagnetic field distribution of the isolated nanocavity. When the ED resonant NP is over a maximum of the electric field, the coupling becomes strong and leads to an enhancement of the scattered light. Nevertheless, a particular position of the NP along the z-direction could lead to a critical coupling for which the scattering is maximum. This phenomenon does not exist when the NP approaches a flat substrate. This is illustrated through an additional result that presents the total energy scattered over 4π sr as a function of the NP position. As illustrated in figure 6, this maximum is obtained for D z = 120 nm in the case of our configuration. For other positions of the NP, its possible to direct the energy in a desired direction. For experimental demonstration, the position of the NP can be spatially controlled by sticking it on the apex of a SNOM tip for example, or by using optical tweezers in non contact regime such as in ref. [START_REF] Hameed | Optical manipulation of nanoparticles by simultaneous electric and magnetic field enhancement within diabolo nanoantenna[END_REF]. The numerical tools developed in this work can be easily adapted to the study of any coupling between two optical nano-systems for the determination of their behavior in both near-and far-field regions.

Appendix A. Mie theory and TFSF-FDTD method

In our study, we used two calculation methods in order to determine the effective scattering cross-section (σ S ) of a spherical NP: The first one is the well-known Mie theory that allows to calculate the coefficients related the each electric and magnetic multipoles through a spherical wave expansion of the scattered field. The weight of each component linked to a multipole (dipole, quadrupole, octopole...) is analytically calculated and its value allows to determine the nature of the excited resonance at a given wavelength. The total scattered energy is obtained by summing coherently all the fields and allows determining the scattering and the extinction optical cross sections. The ratio of the latter to the physical cross section of the NP allows determination of the effective extinction, absorption and scattering cross sections [START_REF] Craig | Absorption and scattering of light by small particles[END_REF]. The second method is the so-called "Total Field -Scattered Field" (TFSF) technique associated to the FDTD (Finite Difference Time Domain) algorithm. The latter is based on spatial discretization of the structure under study through parallelepiped cells in the case of Cartesian coordinates. This kind of mesh is not well suitable to model structures with curved surfaces. Thus, in order to minimize an eventual stair-casing artifactual effect, several techniques were developed in the literature [START_REF] Dey | A modified locally-conformal fdtd algorithm for modeling 3-d perfectly conducting objects[END_REF][START_REF] Dey | A conformal finite-difference time-domain technique for modeling cylindrical dielectric resonators[END_REF][START_REF] Yu | A conformal fdtd algorithm for modeling perfectly conducting objects with curve-shaped surfaces and edges[END_REF][START_REF] Nilavalan | Interpolation/averaging in nonorthogonal fdtd algorithm[END_REF]. Nonetheless, in the case of a purely dielectric material without absorption, the averaging technique seems to be quite powerful. In this case, the dielectric constant of each cell, comprising two or more media of different dielectric constants, is calculated as the average value of the dielectric constant over the cell volume. We have tested this technique and demonstrated that it allows a very good agreement with Mie's theory as shown in next paragraph. The TFSF technique is based on a spatial division of the calculation window into two areas: the area of total field (TF) that contains the scatterers and where the calculated electromagnetic field consists on the sum of the incident and the scattered fields. For the second area, which encompasses the first one, only the scattered field (SF) is calculated by subtracting the incident field over all the external interface of the TF area at each temporal iteration. For a simple case where the NP is immersed in a homogeneous medium (see the figure A1(a)), we have compared the TFSF-FDTD results to the Mie theory ones (see figure A1(b)). This figure presents the variations of the relative deviation between the Mie and TFSF results for both the spectral position and the scattering efficiency as a function of the ratio λ/δ; δ being the FDTD cell size. As it can be shown, the results depends greatly on λ/δ and convergence only occurs when the cell size tends to zero. For our simulations, we set the cell size at 15 × 15 × 15 nm 3 corresponding to a λ/δ = 107 for which the deviation from the Mie theory is of only 0.8% (see figure A1(b)). Let notice that in the case of a substrate (figure A1(c)), the TFSF can be improved by subtracting not only the incident field but also the reflected one on the substrate interface [START_REF] Taflove | Computational Electromagnetics: The Finite-Difference Time-Domain Method[END_REF][START_REF] Dmitry | Magnetic dipole radiation tailored by substrates: numerical investigation[END_REF]. We have developed and tested this configuration in view of obtaining the results presented in section 2.3.

Appendix B. Plane Waves Spectrum theory

The angular spectra of the radiation are calculated using the plane wave expansion (PWE) theory. This theory is based on the decomposition into an angular spectrum (i.e. as a function of the Euler angles θ and ϕ) through the known of the electric field -→ E P (x, y, z 0 , λ) on a plane (P) parallel to xy one and located at D z = z 0 (see figure B1). The PWE theory allows determining the electric field associated with each plane wave propagating along the (θ, ϕ) direction through: We note that the integration appearing in equation B.1 runs on the whole plane (P) from -∞ to +∞. This imposes to know the electromagnetic field everywhere in this plane. In practice, care must be taken to consider a sufficiently large calculation window corresponding to a field that cancels at the edges of the window. After the determination of the three components of the angular spectrum, we calculate the far field intensity by considering the asymptotic limit given in [START_REF] George C Sherman | Asymptotic approximations to angularspectrum representations[END_REF][START_REF] Stamnes | Waves in focal regions hilger[END_REF]

Figure 1 :

 1 Figure 1: (a) Schematic of the considered PCC. (b) Near-field intensity spectrum calculated 15 nm in air above the center of the PCC. (c) to (e) and (g) to (i) are the amplitude distributions of the three components of the electric and magnetic field recorded at 15 nm above the PCC in air. (f) and (j) present the intensity distribution of the electric and magnetic field, respectively, recorded on the same plane as previously. Green circles show the position of the photonic crystal air holes.

Figure 2 :

 2 Figure 2: (a) Evolution of the spectral position (λ) of the ED mode as a function of the refractive index of a 220 nm radius NP. (b) Effective scattering cross section (σ S ) spectra corresponding to the NP with a refractive index (n p = 4.85) indicated by the green arrow in (a) for which the ED wavelength is the same as the PCC resonance one. Two different numerical methods are used: the solid line correspond to the TFSF-FDTD method while the dotted line is calculated using Mie theory.(c,d) are the normalized electric field amplitude distributions recorded in the xz and yz planes respectively at the ED excitation. Note that the incident plane wave is propagating along the z > 0-axis and is linearly polarized along the x-direction.

Figure 3 :

 3 Figure 3: (a) Schematic of the considered configuration. (b) Variations of the ED spectral position (λ ED ) and the corresponding scattering efficiency (σ s ) of a NP (refractive index n = 4.85 nm) coupled with an InP substrate (n = 3.17), as a function of the NP-substrate distance D z . (c,d) are the electric field amplitude distributions recorded in the xz and yz planes at the ED excitation for D z = 45 nm.

Figure 4 :

 4 Figure 4: Normalized near-field (first column), transmission (second column) and reflection (third column) spectra of the PCC coupled to the NP (n p = 4.85, R = 220 nm) as a function of its position along the (x, y, z) axes. The dashed green lines indicate the resonance wavelength (λ mode 0

Figure 5 :

 5 Figure 5: (a) Top-view of the PCC schema showing the three specific positions of the NP marked by 3 colored dots. 3D view of the total (4π sr) radiation angular spectra are shown in (b)for the isolated PCC, in (c) for the PCC-NP system when the NP placed 45 nm above the PCC center at position "A", in (d) for the PCC-NP system when the NP at position "B" and in (e) for the PCC-NP system when the NP displaced along the y-axis and placed at position "C". The 2D circular maps (see last two lines) give the angular spectra calculated in air (z>0) and in the substrate (z<0) for the four corresponding cases (isolated PCC, NP at A, at B and at C respectively).

Figure 6 :

 6 Figure 6: Normalized total radiated energy spectrum of the whole system (PCC+NP) as a function of the NP positions along the z-direction (D z ).

Figure

  Figure A1: (a) Schema of a self-suspended NP in air and surrounded by the TF box (red box) and the SF box (blue box). (b) Relative deviation between the wavelength (blue curve) and the normalized scattering cross section (red curve) calculated by the FDTD code and Mie's theory as a function of λ/δ. (c) Schema of a NP suspended in air and placed at a distance D z above an InP substrate. We can see the TF and SF boxes to penetrate into the substrate.

2 y -k 2 xFigure

 22 Figure B1: (a) Schematic of a wave-vector -→ k located by its Euler angles (θ, ϕ) allowing the determination of its coordinates by simple projection. (b) 3D view of the electric field intensity distribution in the (P) plane. Green circles show the position of the photonic crystal air holes.

  by:I(θ, ϕ, λ) = lim

	z→∞	|| -→ E (θ, ϕ, z, λ)|| 2 ∝	-i λ	cos θ|| -→ E (θ, ϕ, z 0 , λ)|| 2	(B.2)
	where || -→ E || is the modulus of the	-→ E field.

Acknowledgments

Computations have been performed on the supercomputer facilities of the "Mésocentre de calcul de Franche-Comté". This work has been partially supported by the EIPHI Graduate School (contract ANR-17-EURE-0002).

Data availability statement

All data that support the numerical results of this study are included within the article and appendixes.

ORCID iDs

F. I. Baida https://orcid.org/0000-0002-7475-525X