T K N Nguyen

B Dierre

F Grasset

A Renaud

S Cordier

P Lemoine

N Ohashi

T Uchikoshi
email: uchikoshi.tetsuo@nims.go.jp

Formation Mechanism of Transparent Mo6 Metal Atom Cluster Film Prepared by Electrophoretic Deposition

Keywords: Metal atom cluster, molybdenum halides, electrophoretic deposition

Among the solution deposition techniques, electrophoretic deposition (EPD) is a fairly rapid and low cost two-step process well-known for ceramic shaping and conductive surface coating.

Recently we have reported the fabrication of transparent inorganic thin film composed of octahedral molybdenum clusters and shown that the stabilization of the cluster solutions during EPD is very important to obtain high quality films. In this study, we selected the Cs2Mo6Br i 8Br a 6 or ((n-C4H9)4N)2Mo6Br i 8Br a 6 clusters compound ( i = inner ligands, a = apical ligands) as precursors in a Methyl Ethyl Ketone (MEK) solution to fabricate the inorganic Mo6 cluster transparent film on a conductive indium tin oxide (ITO) glass substrate. For the first time, the mechanism for the deposition process of metal clusters is proposed based on results obtained by the combination of complementary techniques such as infrared spectroscopy, X-ray analysis and electron microscopy. The Mo6 film on the ITO glass substrate exhibiting a Br-rich layer is quickly deposited during the first stage of the EPD followed by a multilayer structure consisting of two types of cluster compounds, [Mo6Br i 8Br a 4(H2O) a 2] and (H3O)2[Mo6Br i 8Br a 4(OH) a 2].

Introduction

The octahedral molybdenum halide cluster has been extensively studied since the 1970s after the discovery in 1971 of the AyMo6Q8 (Q = chalcogen, A = cation) compounds, the so-called "Chevrel phases" 1,2 . More recently, their high potential use in nanoarchitectonic devices has been reported [3][4][5] . The [Mo6Br i 8Br a 6-xL a x]

x-2 cluster unit, based on the [Mo6X i 8]

x+ core family (X i =Br,

Cl, I at face-capping inner position) capped by 6 halogens or functional groups at the terminal apical position (L a = donor ligand; Br, Cl, I, OH, H2O, CF3COOH,…) 2,5 , has recently been studied in applications of pH-sensor 6 , singlet oxygen sensitization 7 , photocatalytic reactions 4,8 , solar cells [9][10][11] or biotechnology [12][13][14][15][16][17][18] . The synthesis of hybrid nanocomposites based on metal atom clusters has also attracted the interest in optoelectronic applications 4 , but only a few studies have reported the preparation of transparent thin films [19][20][21] .

The surface of glass is often coated for intrinsic reasons (safety, convenience or better stability)

or specific functionalization (ultraviolet or heat protection), water-repellent coating, shatter protection, easy-to-clean and self-cleaning coatings [22][23][24][25] . As a method for these applications, electrophoretic deposition (EPD) would be an attractive process for preparing transparent, smooth and thickness-controlled thin coating. In addition, this process can be easily scalable to an industrial level. The EPD process consists of two steps: i) the electrophoresis of charged particles moving toward the electrode under an electric field, and ii) the deposition of the particles on the electrode surface by coagulation [26][27][28] . Our group has recently reported the fabrication of transparent inorganic thin films composed of octahedral molybdenum metal clusters deposited on indium tin oxide (ITO) glass slides by EPD 21 . In the study, Cs2Mo6Br14 cluster compound was dissolved in various solvents to prepare the EPD solutions. The films exhibit a high transparency in the visible light wavelength range and good absorption at wavelengths lower than 400 nm and higher than 1100 nm, which is of significant interest for ultraviolet (UV) and/or near infrared (NIR) cut filter applications, such as windows and solar concentrators [29][30][31][32][33][34][35][36][37][38][39] . Moreover, a luminescent emission in the wavelength range of 650 to 900 nm resulting from the Mo6 cluster units could be observed and depends on the EPD solutions and conditions. The efficiency and lifetime of the luminescent emission of the molybdenum halide cluster depend on the lowest excited state determined by the prominent improvement of the halides 40 . The electronic interaction between the ligands of the

[Mo6X i 8L a 6]
2-clusters in the network and different kinds of cations depends on the fabrication process of the metal cluster for films, crystals, powders or composites 4,9,41 .

From the results of the previous report, Methyl Ethyl Ketone (MEK) was evidenced as the stably dispersing medium in order to remain the original luminescent characterization of the Cs2Mo6Br14 cluster precursors 21 . In this study, we used TBA2Mo6Br14 (TBA = Tetra Butyl Ammonium, ((n-C4H9)4N)) and Cs2Mo6Br14 as the raw materials. The main objectives are to (i) select the best cluster precursor in order to obtain good quality Mo6 cluster films and (ii) evaluate the stability of the Cs2Mo6Br14 and TBA2Mo6Br14 cluster compounds in Methyl Ethyl Ketone (MEK). Indeed, in MEK solution, the Br -ligands may be partially replaced and released out of the octahedral cluster during the dissolving and depositing procedures. The deposition mechanism to fabricate a transparent Mo6 cluster film by the EPD process and the structure of the obtained film were studied for the first time as far as we know by using complementary tools and techniques; i.e., color 3D laser microscopy, field-emission scanning electron microscopy (FE-SEM), X-ray diffraction (XRD), scanning transmission electron spectroscopy (STEM), transmission electron spectroscopy-energy dispersive spectroscopy (TEM-EDS), Fourier transform infrared spectroscopy (FTIR), X-ray fluorescence spectroscopy (XRF) and X-ray photoelectron spectroscopy (XPS).

Experimental Preparation of cluster films by EPD

The Cs2Mo6Br14 powder (CMBP) was synthesized from MoBr2 and CsBr agents by a solid-state method at high temperature, while the (TBA)2Mo6Br14 powder (TMBP) was prepared from CMBP and TBABr agents by solution chemistry 42 . Both powders were separately dissolved in reagent grade methyl ethyl ketone (MEK) (99 %, Nacalai Tesque, Inc., Kyoto, Japan) at the concentration of 5g/L with agitating by a magnetic stirrer for 20 minutes. ITO glass, washed with distilled water and acetone during ultrasonication, was connected to a Keithley Model 2400

Series Source meter as the anodic substrate and a stainless steel sheet as the cathode. Based on our previous study 21 , the voltage from 7.5 to 25 V and deposition time up to 40s were applied.

The best quality Cs2Mo6Br14 films (CMBF) and TBA2Mo6Br14 films (TMBF) were obtained by applying voltages at 13 and 17 V, respectively. The Mo6 cluster-based films deposited on the ITO glass substrates, whose surface area were about 1×1.5 cm 2 , were characterized after drying in air for 24h.

Characterization methods

The zeta potential and electric conductivity of the suspensions were measured by a zeta-potential analyzer (Malvern Instrument, Ltd., Zetasizer Nano Z). The particle size of the cluster in an MEK suspension was measured by a dynamic light scattering (DLS) technique using a Nano Q V2.5.9.0 analyzer (Cordouan Technologies) at the wavelength of 657 nm. A high resolutioncolor 3D laser microscope with a 408 nm laser was used for the measurement of the thickness of the films formed at different deposition times. The surface morphology of the thin films was observed by an FE-SEM (Hitachi S4800). The crystallographic structure of the films was determined by XRD (SmartLab, RIGAKU, 40 kV and 30 mA) in the 2θ angle range from 5° to 55° and at the scan speed of 1°/min with Cu Kα radiation (λ = 1.54 Å). The chemical state and composition were verified by coupling several techniques, such as FTIR (Thermoscientific Nicolet 4700) in the wavenumber range from 4000 to 400 cm -1 and XRF measurement (EZX Primus II, Rigaku). XPS spectra of the powder and film were measured by a PHI Quantera SXM (ULVAC-PHI) using Al Kα radiation at 20 kV and 5 mA, neutralization by Ar + , pass energy of 55 eV and the take-off angle of 45°. High resolution observation of the cluster-deposited films was performed by a TEM (JEOL JEM 2100F) equipped with an EDS analysis device.

Results and discussion

Characterization of the solutions:

The CMBP and TMBP cluster-containing suspensions in the MEK solvent presented an electric conductivity of approximately 0.3 mS/cm and a negative zeta potential in the range of 20-25 mV as listed in Table . 1. These values are slightly different compared to those in our previous study.

It is related to the stirring time, which was optimized in order to achieve a good colloidal liquid and a thick film by the EPD. The dissolution of TMBP in MEK is quicker than for CMBP due to the better compatibility of the TBA + cation in an organic solvent. Nevertheless, the hydrodynamic diameters of the CMBP and TMBP clusters in the MEK solution (1g/L) obtained by DLS were in the same range of 11±1 and 13±1 nm, respectively (Table 1). It could correspond to an aggregate of less than 1000 [Mo6Br i 8Br a 6]

2-cluster units 42 .

Structural and chemical characterization of the film:

In order to discuss the deposition process, the chemical composition and structure were analyzed for the CMBF and TMBF films prepared at different voltages and deposition times. The Br/Mo atomic ratios evaluated by XRF of the Mo6 cluster-based films obtained by the EPD as a function of the deposition time and applied voltage are shown in Figure 1. At 13V, which was used for the CMBF, was applied to TMBF, a good homogeneity was not obtained because the film was easily broken by the re-dissolution in the MEK solution after stopping the applied voltage. The Br/Mo atomic ratio as a function of the deposition time were determined at 13 V for CMBF and 17 V

for TMBF, which are the optimal voltages to attain the best uniformity in the films, respectively, for those as a function of the applied voltage at 10s for both CMBF and TMBF.

The atomic ratios of the Br and Mo atoms were estimated to be 2.20 (±0.11) and 2.25 (±0.13) for the CMBP and TMBP, respectively, which were in fairly good agreement with the theoretical value of 2.33 (14 Br atoms/6 Mo atoms). In addition, the Cs + atoms were detected in CMBP solution, but it was not possible to measure any signal showing the existence of Cs atoms in all of the CMBF films, suggesting the absence of the Cs + cations in the electrodeposited films.

The curves in Fig. 1a show a decrease in the Br/Mo atomic ratio of the films from 7.2 (±0.3) (TMBF) or 4.9 (±0.2) (CMBF) to 2.8 (±0.2) in the first 10s when the applied voltage increased from 7.5 to 25V. These values are much higher than the theoretical value (2.33) and the experimental values determined for CMBP and TMBP, suggesting that the Br-atom rich Mo6

cluster layers are preferentially deposited on the ITO glass for the first 10s, especially at the low applied voltages.

In a similar way, the Br/Mo atomic ratio decreases to the theoretical index (2.33) for the higher deposition time (at a fixed applied voltage) and reaches the theoretical value after 30s (TMBF) or 40s (CMBF) of the deposition. These data also indicate that a Br-rich layer containing Mo6clusters is deposited during the first seconds of the EPD process, followed by the deposition of the Mo6-cluster units. One explanation of these two steps during the deposition process can be expressed as follows: During the dissolution process, the free Br -anions are generated by the substitution of the solvent molecules to form [Mo6Br i 8Br a 6-x(solvent) a x]

x-2 cluster units, then move toward the electrode faster than the bigger [Mo6Br i 8Br a 6-x(solvent) a x]

x-2 cluster units under the impact of an electric field.

Considering the Br/Mo atomic ratio of CMBF and TMBF versus the applied voltage and deposition time, the EPD parameters were selected to obtain the Mo6-cluster films with a Br/Mo atomic ratio close to the theoretical one (2.33); i.e., the deposition time was longer than 20s, and the voltage fixed at 13V for CMBF and 17V for TMBF.

The thickness of the Mo6 cluster films prepared from the CMBP and TBMP suspensions in the MEK solution at the different deposition times of 10, 20, 30 and 40 seconds was measured using a color 3D microscope as presented in Table .2. The average values were recorded at least on three samples and three positions on each films. It can be seen that CMBF reached a maximum thickness of 1.89 ± 0.08 μm for the deposition time of 40 s (applied electric field of 13 V), while TMBF presented a maximum thickness of 1.68 ± 0.07 after 30 s of the deposition. For higher deposition at 40 s, the thickness of the TMBF sample decreased. Such phenomena was previously reported for CMBF 21 and was explained by the fact that the deposited Mo6 film on the ITO glass with a high thickness plays the role as a resistance layer and then reduces the influence of the electric field on the movement of further cluster particles. Moreover, the diffusion of cations to the outside cluster layer will break up the balance of anion cluster particles with the charge on the electrode. Consequently, the Mo6 cluster units will be released out of the layer by the re-dissolution in the MEK solution after 40 s for the TMBF film. Considering these results, the use of CMBP for fabricating the Mo6 cluster thin film seems to be better for controlling the thickness of the Mo6 cluster-deposited films.

The surface morphology of the Mo6 cluster films was characterized by SEM on the 20s deposition time samples. The SEM images of the CMBF20s and TMBF20s films show the evidence of homogeneous and flat morphology area of the film even at low magnifications of 1000 and 50000 times (Fig. 2a). However, the surface of the film seems not to be stable under high vacuum because many hairline cracks with 500 nm widths were observed. Such cracks are probably generated on the smooth surface of the films by the quick evaporation of the absorbed water molecules or residual solvent, especially during the evacuation in the SEM chamber. The cross-sectional SEM images of the high resolution images of the CMBF20s and TMBF20s films shown in Fig. 2b exhibit a layered structure formed on the glass matrix, the ITO layer (0.2 µm) and the Mo6-cluster layer (~1.2 µm). These values are in good agreement with the color 3D microscope results. A larger view of the cross-section of the CMBF40s film confirmed the homogeneity of the deposition process and also the thickness could reach 2 µm as shown in x-2 cluster units during the EPD process. The intensity of O-H vibrational band in the wavenumber range between 3600 cm - 1 (ν3) and 3200 cm -1 sharply increases, continuously including a stretching vibration (ν3) at 3600 cm -1 , a stretching vibration (ν1) at 3450 cm -1 and an overtone of the bending vibration (2ν2) at 3250 cm -1 . The TMBF30s and CMBF30s films have similar FTIR spectra consisting of the asymmetric O-H stretch vibrational band (2ν2) at 3250 cm -1 , combined bending modes (ν2) and vibrations (νL) of the H-O-H mode at 2304 cm -1 , and the H-O-H bending mode at 1589 cm -1 , 1400 cm -1 (ν2) and 792 cm -1 (νL) 43 . The band assigned to O-H of the H-O-H stretch vibration (ν3) at 3527 cm -1 in the CMBP sample, representing the adsorbed free H2O molecules, shifted to the O-H of the H-O-H stretch vibration (2ν2) at 3250 cm -1 in the CMBF30s film, the signal of the hydrogen bond [44][45][46] . For this reason, we assumed that the Br apical ligands were exchanged by the H2O molecules or OH -groups, consequently, new hydrogen bonds were formed in the Mo6 cluster film prepared by the EPD process. In addition, the obvious increase in the intensity of the O-H vibrational band in the FTIR spectrum originated from the increase of the H3O + cation absorption in the film, which would be generated on the surface of the ITO glass anode by the oxidization reaction of the H2O molecules contaminated in the suspension in order to neutralize the negative charges of the [Mo6Br i 8Br a 6-xL a x]

x-2 cluster molecules 21 .

The Br/Mo atomic ratio of the CMBF (4.9±0.2) after the first 10 seconds is lower than that of the TMBF (7.2±0.3) at 7.5 V (Fig. 1a), which suggests that the Br apical ligands of CMBP are less labile than TMBP probably due to a limited exchange by the H2O molecules or OH -groups in the MEK solution. Thus, the CMBP suspension was selected as the precursor solution to prepare the Mo6 cluster film (CMBF30s) that was investigated for the chemical elements through the XPS and EDX measurements. In order to elucidate the difference in the structure and chemical components of the [Mo6Br i 8Br a 6-xL a x]

x-2 cluster unit, the XPS spectra of the O 1s (Fig. 5a), Cs 3d and 4d (Fig. 5b), Mo 3d (Fig. 5c) and Br 3d region (Fig. 5d, 5e and 5f) of the Cs2Mo6Br14 powder (CMBP) and the Mo6 film (CMBF30s) were recorded and analyzed. The results are presented in Figure 5 and Table 3.

The XPS spectra of O 1s evidences a peak at 533.2 eV for both samples (Fig. 5a), which is assigned to either H2O molecule or OH -group. However, the intensity of the signal of O 1s atoms in the CMBF30s film is higher than that in CMBP. This result supports the substitution of apical Br atoms by H2O molecules or OH -groups during the dissolution process and the absorption of the H3O + cations during the deposition on an ITO glass. The XPS spectra of Cs 4d evidences the peaks at 738 eV (3d3/2), 724 eV (3d5/2), 76.1 eV (4d3/2) and 78.4 eV (4d5/2) only for the CMBP sample (Fig. 5b), confirming (i) the presence of Cs atoms as counter cation in the precursor and

(ii) the absence of Cs atoms in the film after the EPD process. The Mo 3d XPS spectra show equivalent peak intensity at 229.5 (3d5/2) and 233 eV (3d3/2) for both the CMBP and CMBF30s samples (Fig. 5c). Consistently with the XRF and EDS results, it confirms the presence of octahedral Mo6 clusters in the film after the deposition.

The peaks in the Br 3d region of the CMBP and CMBF30s samples have relatively similar values (Fig. 5d). The peaks at 70.4 eV and 68.5 eV are assigned to the inner Br atoms and terminal apical Br atoms, respectively. While the peak intensities corresponding to the inner Br atoms of the CMBP and CMBF30s samples are equivalents, those of the Br atoms at the face-capping apical position are clearly different, as confirmed by their deconvolutions (Fig. 5e and5f). Indeed, a reduction of the Br 3d5/2 peak intensity at 68.5 eV of apical Br atoms is clearly observed for the CMBF30s sample. Similar feature was also observed for the Mo6 clusters immobilized on graphene oxide 8 . It can be interpreted as a new valence linking between the Mo6 octahedral clusters and OH -functional groups by the removal of apical Br atoms from the cluster unit.

The atomic ratio of CMBP sample estimated from XPS spectra of Br 3d and Mo 3d is 14.1/6 (≈ 2.35). This value is very close to the theoretical index (2.33), while that in the CMBF30s sample is lower: 12.1/6 (≈ 2). It indicates that in average each octahedral [Mo6Br i 8Br a 6]

2-cluster unit deposited in the film lost 2 Br atoms (Table .3). This ratio is inferior to that obtained from the XRF measurements (Fig. 1), probably due to the free Br -anions lost during the preparation of the XPS measurement. Consequently, only the Br atoms directly linking with the Mo6 clusters are then detected by XPS.

From XPS and FTIR analyses, it could be conclude that the dissolution in the MEK solvent and the EPD process induces the substitution of two Br apical ligands from the [Mo6Br 2] cluster units. Both of them are known to exhibit a stable octahedral structure [49][50] .

The chemical analysis, shape and size diameter of the first deposited cluster units inside the EPD film (CMBF30s) were clearly characterized by the TEM-EDX and STEM images (Fig. 6). The EDX-TEM mapping (Fig. 6a) shows homogenous distribution of Br and Mo atoms, while the presence of Cs atoms is not significantly recognized. The later point supports the absence of Cs atoms previously suggested by XRF and XPS analyses. The TEM images indicate that the CMBF30s sample consists of many crystallized nanoparticles of average diameter about 6 nm, homogeneously distributed near the substrate (Fig. 6b and6c). Considering that a single cluster unit presents an average volume of 1.2 nm 3 42 , it can be deduced that a nanoparticle of 6 nm diameter is formed by less than 90 [Mo6Br i 8Br a 6-x(OH) a x] 2-or [Mo6Br i 8Br a 4(H2O) a 2] cluster units, the number of cluster units per nanoparticle being dependent on the cluster unit density. Considering that in MEK solution, the cluster units are linked together by hydrogen bonds, it leads to the formation of bigger negative charged spherical particles, which are accumulated on the electrode surface through the effective impact of the electric field to form the multilayers.

Therefore, the structure in the first 100 nm layer is proposed to be an arrangement of 2] cluster units. A schematic representation is shown in Figure 7. It is important to note that in the proposed model, the nature of the cations, namely Cs + or TBA + , of the initial powder precursors does not seem to play a key role at the surface of the electrode. These cations seem to play only an important role in the dissolution of metal cluster in the solvent, and consequently, on the rate of the deposition and thickness.

Conclusions

This study describes on the preparation, and chemical, structural and optical analyses of smooth 2-paired with two (H3O) + , and (iii) the exchange of Cs + and TBA + cations by H3O + in order to keep a neutral charge. Finally, the films present a high transmittance in the visible range and strong absorption in UV and NIR ranges.

Figure 2c at

  Figure 2c at 40s.

  nanoparticles of 6 nm diameter. Consequently, the micron-sized layer is formed by the closed-pack stacking of different types of the species: (H3O) + cations, free Br -anions, spherical nanoparticles containing a large number of units.Based on the variation in the Br/Mo atomic ratios as a function of the applied voltage and deposition time and the modifications of the chemical linking observed by FTIR, XPS and TEM, we propose the following mechanism for the fabrication of the Mo6 cluster film by the EPD process. In the first seconds of the EPD process, high mobile Br -anions (originated from OH -to Br -substitution in MEK solution) move toward the ITO glass anode in the electric field to form a Br -rich sublayer. The [. Due to lower pH close to the electrode, these later are quickly protonated to form [units interact together by hydrogen bonds to form approximately 6-nm size nanoparticles. The outer layers would mostly contain the (H3O)2units which are simply packed by the EPD process (high velocity and collision of clusters). In summary, the Mo6 film is a nanocomposite material with a multilayer structure made by the stacking of Br -anions,

  Mo6 cluster-based films deposited by the EPD process from the Cs2Mo6Br14 and (TBA)2Mo6Br14 cluster compounds in the MEK solution. The combination of X-ray analyses, electron microscopy and optical characterizations allow us to propose for the first time a mechanism for the multilayer structure deposition process of the metal atom clusters. The deposition of a Brrich layer during the first seconds of the EPD process is proved through the determination of the Br/Mo atomic ratio in relation to the deposition time and applied voltages. Moreover, three important points should be pointed out: (i) the exchange of 2 Br apical ligands by OH -or H2O, (ii) the existence of two kinds of clusters in the Mo6 film, i.e., [Mo6Br i

Figure 1 .

 1 Figure 1. The Br/Mo atomic ratio in the CMBF and TBMF films as a function of (a) the applied voltage for a deposition time of 10s, and (b) the deposition time for an applied voltage of 13V and 17V, respectively.

Figure 2 .

 2 Figure 2. The FE-SEM images of a) the surface and b) the cross-section of the Mo6 cluster films fabricated by EPD: CMBF20s at 13 V (top) and TMBF20s at 17 V (bottom). c) FE-SEM image of the cross-section of the CMBF40s film.

Figure 3 .

 3 Figure 3. The XRD patterns of ITO glass, CMBP, TMBP, CMBF30s and TMBF30s samples.

Figure 4 .

 4 Figure 4. The FT-IR spectra of CMBP, TMBP, CMBF and TMBF samples.

Figure 5 .

 5 Figure 5. The XPS spectra of O 1s (a), Cs 3d and 4d (b), Mo 3d (c), and Br 3d region (d), deconvolution spectra of Br 3d region of CMBP (e) and deconvolution spectra of Br 3d of CMBF (f).

Figure 6 .

 6 Figure 6. a) The EDX-TEM mapping of Cs, Mo and Br elements and b) TEM image (left) and STEM mode image (right) of the Mo6 cluster nanoparticles from film by EPD. c) Spherical Mo6 cluster nanoparticles included in the CMBF30s.

Figure 7 .

 7 Figure 7. Schematic representation of the multilayered structure of Mo6 cluster thin film

Figure 1 .

 1 Figure 1. The Br/Mo atomic ratio in the CMBF and TBMF films as a function of (a) the applied voltage for a deposition time of 10s, and (b) the deposition time for an applied voltage of 13V and 17V, respectively.

Figure 2 .Figure 3 .

 23 Figure 2. The FE-SEM images of a) the surface and b) the cross-section of the Mo6 cluster films fabricated by EPD: CMBF20s at 13 V (top) and TMBF20s at 17 V (bottom). c) FE-SEM image of the cross-section of the CMBF40s film.

Figure 4 .

 4 Figure 4. The FT-IR spectra of CMBP, TMBP, CMBF and TMBF samples.

Figure 5 .

 5 Figure 5. The XPS spectra of O 1s (a), Cs 3d and 4d (b), Mo 3d (c), and Br 3d region (d), deconvolution spectra of Br 3d region of CMBP (e) and deconvolution spectra of Br 3d of CMBF (f).

Figure 6 .

 6 Figure 6. a) The EDX-TEM mapping of Cs, Mo and Br elements and b) TEM image (left) and STEM mode image (right) of the Mo6 cluster nanoparticles from film by EPD. c) Spherical Mo6 cluster nanoparticles included in the CMBF30s.

Figure 7 .

 7 Figure 7. Schematic representation of the multilayered structure of Mo6 cluster thin film

Table 1 .

 1 Physical properties of the CMBP and TMBP cluster suspensions in MEK solutions.

	Suspension	Conductivity	Zeta potential	DLS Particle size (nm)
	(5g/L in MEK)	(mS/cm)	(mV)	(1g/L in MEK)
	CMBP	0.33±0.01	-21.0±1.2	11±1
	TMBP	0.30±0.03	-26.4±1.7	13±1

Table 2 .

 2 Thickness (average) of the films deposited from the CMBP and TMBP cluster suspensions in MEK solutions for deposition times.

	Deposition time	Thickness of films fabricated from MEK solution (μm)
	(Second)	Cs2Mo6Br14 solution (13V)	(TBA)2Mo6Br14 solution (17V)
	10	CMBF10s	0.81 ± 0.03	TMBF10s	0.49 ± 0.05
	20	CMBF20s	1.05 ± 0.09	TMBF20s	1.23 ± 0.09
	30	CMBF30s	1.54 ± 0.09	TMBF30s	1.68 ± 0.07
	40	CMBF40s	1.89 ± 0.08	TMBF40s	1.29 ± 0.12

Table 3 .

 3 The percentage of the component elements in the CMBP and the CMBF30s films prepared at 13 V estimated by the peak analysis of the XPS measurements.

	Element	O 1s	C 1s	Mo 3d Cs 4d	Br 3d	Cs: Mo: Br (Wt.% at.)
	CMBP	4.9%	22.3% 20.2% 5%	47.6%	1.5: 6.0: 14.1
	CMBF30s 13.7 % 37.2% 16.1% -	32.5%	0: 6.0: 12.1

Acknowledgement

This study was carried out as part of the France-Japan international collaboration framework through the Laboratory for Innovative Key Materials and Structures (UMI 3629-LINK Center).

The study was financially supported by the National Institute for Materials Science (NIMS), Saint-Gobain (France), CNRS and Université de Rennes 1 (UMR 6226). The authors wish to thank David Lechevalier and Dr. Mari Kono of Saint-Gobain KK (Tokyo, Japan) involved in LINK and related activities. We would like to thank Dr. Hirohito Ohata at NIMS for his help with the XPS measurements, Dr. Yoshio Matsui for his help with the STEM measurements and Dr. Chika Matsunaga at AIST for her valuable comments about the EPD experiments. The authors thank the ScanMAT platform at the University of Rennes 1 for the FE-SEM image (Figure 2c).


Antoine Lemaire

Arnaud Perona
email: arnaud.perona@univ-perp.fr

Matthieu Caussanel

Hervé Duval

Alain Dollet

Open-circuit voltage decay: moving to a flexible method of characterisation

published or not. The documents may come

Introduction

Minority carrier lifetime measurement is of major importance to optimize electronic and optoelectronic devices such as solar cells [START_REF]Handbook of photovoltaic science and engineering[END_REF]. It enables to tailor solar cells characteristics such as bulk thickness and doping levels. Moreover, when lifetime can be measured between process steps, it also gives a follow-up of these various steps quality. Indeed, lifetime depends on material nature and on the various treatments the material or device have been subjected to [START_REF] Schroder | Semiconductor material and device characterization[END_REF]. That's why it informs on material purity combined with interface quality. Measured lifetime always results from several electrons-holes recombination mechanisms such as SRH (Shockley-Read-Hall), Auger and radiative [3] and thus is called effective lifetime τ ef f . Nowadays, photoconductance decay (PCD) [START_REF] Yablonovitch | A contactless minority lifetime probe of heterostructures, surfaces, interfaces and bulk wafers[END_REF][START_REF] Sinton | Quasi-steady-state photoconductance, a new method for solar cell material and device characterization[END_REF] is commonly used for silicon solar cells. It is a contactless method which allows really fast characterization. However, the method expects equipments (light flash and filter) specific to the characterized semiconductor and does not allow to measure short lifetimes (< 1 µs). Time-resolved photoluminescence (TRPL) [START_REF] Ahrenkiel | Measurement of minority-carrier lifetime by time-resolved photoluminescence[END_REF] allows lifetime measurement down to the ps scale. It is usually applied to III-V and II-IV semiconductors for that reason. This method is usually slow when applied to materials having poor quality and/or bad radiative properties. Furthermore, it is an expensive method compared to PCD.

By comparison, Open-Circuit Voltage Decay (OCVD) is fast, lowcost (as the PCD) and enables characterization of COTS devices. Moreover, it allows characterization of any semiconductors, no matter its band gap: wide/narrow and direct/indirect. However, due to several drawbacks it has been scarcely employed to measure short lifetime. It requires electrical contacts plus a junction (p-n or Schottky) and its signal is often blurred by the space charge region (SCR) effects (capacitive and generations/recombinations (G/R) current) [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF]. These effects make the method rather subjective for measuring τ ef f . A well-known method, the "Suns-Voc" [START_REF] Sinton | A QSS open circuit voltage method for solar cell characterization[END_REF], has been developed on the basis of the QSS-PCD [START_REF] Sinton | Quasi-steady-state photoconductance, a new method for solar cell material and device characterization[END_REF] to avoid SCR capacitive effect. However, the bulk doping level of the junction must be known.

Our approach is to take advantage of capacitive effect and to the best of our knowledge, this has never been investigated. The SCR capacitance of an asymmetric junction (n-p + or p-n + ) depends on the doping level of the lowest-doped region and built-in voltage, whereas SCR G/R current can be first approximated to a simple shunt resistance. The flexibility of the OCVD makes it extremely attractive since none of the methods faced here can compete. Thus, according to theory these quantities can be extracted by assuming a proper fitting model. Therefore, this is the purpose of this paper.

OCVD Theory

OCVD theory was initially exposed in 1955 by Gossick [START_REF] Gossick | On the Transient Behavior of Semiconductor Rectifiers[END_REF] and Lederhandler & Giacoletto [START_REF] Lederhandler | Measurement of minority carrier lifetime and surface effects in junction devices[END_REF]. It has been first developed to characterize germanium (Ge) but it was quickly applied to silicon (Si). In theory, it is suitable to any semiconductors. Its working principle is based on the interpretation of the open-circuit voltage decay across a junction after an excitation pulse. First, the junction is forward biased either electrically or optically [START_REF] Mahan | Measurement of minority carrier lifetime in solar cells from photo-induced open-circuit voltage decay[END_REF]. Then, open-circuit condition is provided by a pulse generator or by a switch like a MOSFET. Finally, the voltage decay across the open-circuited junction is recorded over the time. There are two types of circuits using a switch: the so-called series circuit where switch and characterized junction are in series [START_REF] Green | Minority carrier lifetimes using compensated differental open circuit voltage decay[END_REF][START_REF] Caussanel | Doping-Type Dependence of Damage in Silicon Diodes Exposed to X-Ray, Proton, and He + Irradiations[END_REF][START_REF] Bellone | An Analog Circuit for Accurate OCVD Measurements[END_REF][START_REF] Lemaire | Mesure par OCVD de la durée de vie des porteurs minoritaires dans des jonctions en GaSb, en GaAs et en Si : simulations et expérimentations[END_REF] (voltage bias source) and the so-called parallel where they are in parallel [START_REF] Lemaire | Mesure par OCVD de la durée de vie des porteurs minoritaires dans des jonctions en GaSb, en GaAs et en Si : simulations et expérimentations[END_REF][START_REF] Deshmukh | Measurement of solar cell ac parameters using the time domain technique[END_REF][START_REF] Ariel | A computerized method for carrier lifetime measurement in PN junctions at high and low-level injection[END_REF] (current bias source). The simplicity of the circuit makes the method low-cost and easy to operate. Usually, OCVD signal exhibits 3 regions (see Fig. 1). Region (i) is an abrupt voltage drop mainly due to series resistance [START_REF] Green | Minority carrier lifetimes using compensated differental open circuit voltage decay[END_REF].

Region (ii) is due to diffusion/recombination of minority charge carriers into quasi neutral regions (QNR). This is the region of interest for lifetime extraction. Ideally, region (ii) presents a linear decay whose slope is inversely proportional to effective lifetime τ ef f :

τ ef f = - ηk B T q 1 dV /dt (1)
where η is the p-n junction ideality factor, k B the Boltzmann constant in J K -1 , T the temperature in K, q the elementary electric charge in C and dV /dt the voltage decay rate in V s -1 . η depends on injection level and dominant recombination mechanism [START_REF] Sze | [END_REF]. This linear behaviour eases greatly τ ef f extraction and makes the method really simple in terms of analysis. However, the linear OCVD decay has an intrinsic limit, depending on materials [START_REF] Lemaire | Mesure par OCVD de la durée de vie des porteurs minoritaires dans des jonctions en GaSb, en GaAs et en Si : simulations et expérimentations[END_REF]. Indeed, for short lifetimes or highly-doped junctions, there isn't region (ii) any more [START_REF] Lemaire | Mesure par OCVD de la durée de vie des porteurs minoritaires dans des jonctions en GaSb, en GaAs et en Si : simulations et expérimentations[END_REF]. Furthermore, for intermediate lifetimes or doping levels, the method becomes rather subjective because it depends on where you clearly identify the linear region. In transient regime, minority charge carriers stored into QNRs behave like a capacitance. This diffusion capacitance [START_REF] Castañer | Investigations of the OCVD transients in solar cells[END_REF] (in F cm -2 ) is defined by:

(i) (ii) (iii) t (µs) V 1 (V)
C D = q ηk B T τ ef f J(V ) (2)
where J(V ) is the current density of the p-n junction in A cm -2 derived from the Shockley equation. Region (iii) is due to SCR capacitance which typically decreases decay rate [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF]. Experimentally, it makes region (ii) usually tough to identify and thus unreliable. Moreover, the lower the lifetime, the more difficult the region (ii) identification. SCR capacitance depends on material and junction design (doping levels) and is defined [START_REF] Sze | [END_REF] according to the following formula:

C SCR = W SCR = q N A N D 2(N A + N D )(V bi -V ) (3)
where C SCR is the SCR capacitance in F cm -2 , the semiconductor permittivity in F cm -1 , W SCR the thickness of the SCR in cm, N A (N D ) the acceptors (donors) doping level of the p-doped (n-doped) QNR in cm -3 , V bi the built-in potential in V and V the voltage across the SCR in V. V bi is calculated [START_REF] Sze | [END_REF] through:

V bi = k B T q ln N A N D n 2 i ( 4
)
where n i is the intrinsic carriers density in cm -3 . When a p-n junction is asymmetric (n-p + or p-n + ), (3) may be written in the following form:

C SCR = q N l 2(V bi -V ) (5)
where N l is the doping level of the lowest-doped QNR in cm -3 . SCR capacitance actually acts in parallel with diffusion capacitance, thus global capacitance is always the sum of them. Nevertheless, for a voltage higher than V bi , SCR region and its capacitance disappear. In this case, OCVD signal is fully due to diffusion capacitance. As the voltage across the diode decreases over the time, the inner electric field E i and SCR capacitance reappear. Then, SCR capacitance decreases OCVD decay rate when its value becomes higher than diffusion capacitance. SCR capacitance makes decay rate lower because E i induces a strong increase of instantaneous resistance [START_REF] Mahan | Measurement of minority carrier lifetime in solar cells from photo-induced open-circuit voltage decay[END_REF].

However, decay rate after region (ii) can also be increased with low shunt resistance [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF]. Shunt resistance effect is often attributed to G/R into the SCR [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF] and perimeter recombinations in case of small area samples like solar cells dedicated to CPV [START_REF] Ochoa | 3-D modeling of perimeter recombination in GaAs diodes and its influence on concentrator solar cells[END_REF][START_REF] Espinet-Gonzalez | Analysis of perimeter recombination in the subcells of GaInP/GaAs/Ge triple-junction solar cells[END_REF]. Therefore region (iii) is filled with information. For a given material, (5) clearly shows that SCR capacitance depends on N l , V bi and V . However, to the best of our knowledge, region (iii) has never been exploited as it could be. Thus, applying a fitting model to OCVD measurements theoretically allows to extract N l and V bi . Shunt resistance R sh is also extractable when it has a strong influence on OCVD signal. In addition, being able to fit the entire signal allows more objective and accurate lifetime determination when region (ii) identification is difficult.

Modelling

OCVD signal modelling was carried out with Python programming language [START_REF] Python | Welcome to Python[END_REF] which provides many powerful libraries that allow data treatment, numerical computation, optimization and so on. Moreover it works on almost all operating systems (Windows, Linux/UNIX, Mac OS X and so on). We used Python 3.7.1 with Anaconda 4.2.0 distribution on Spyder environment. We developed a model and first compared it to Mahan and Barnes (M&B) previous work [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF]. Then we studied the influence of several parameters such as τ ef f , N l , V bi and R sh on OCVD signal.

One-diode model

The model we developed is based on the lumped element model of one asymmetric diode in transient regime [START_REF] Colinge | Physique des dispositifs semi-conducteurs[END_REF] which is represented Fig. 2.

Our work has been inspired by M&B one's [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF] who simulated OCVD

• D 1 j D 1 • R sh j R sh C SCR j C SCR C D j C D j tot V 1
Fig. 2: Lumped element model of one asymmetric diode in transient regime.

signal in 1981. They used two different models: a one-diode model to simulate SCR capacitance effect and a two-diode model to simulate effect of the G/R current into the SCR. The second diode from the model #2 is similar to a shunt resistance having a dynamic behaviour.

In our model, the G/R current into the SCR was assigned to a constant resistance so-called shunt resistance. In first approximation, this has no major influence on the global behaviour because the G/R current into the SCR is known to occur at low voltage over a short range. However we applied a dynamic behaviour (voltage dependence) to the SCR capacitance while M&B considered a constant one. Moreover, we simulated influences of SCR capacitance and leakage current into the SCR together whereas they did it separately: at least they didn't show influence when both have significant effects. Therefore, they didn't show their mutual influence. Unfortunately, we didn't take into account the series resistance because it requires a current flowing at t = 0, so it does not fulfil the open-circuit condition. Details of the models and their differences are summarized in the Table 1.

The lumped element model of OCVD represented in Fig. 2

leads to the following equation expressing the total current density in

J tot = J D1 + J Rsh + J C SCR + J C D (6)
where J D1 , J Rsh , J C SCR et J C D are the current densities through diode D1, shunt resistance R sh , SCR capacitance C SCR and diffusion capacitance C D respectively. J D1 corresponds to the Shockley equation mentionned in [START_REF] Schroder | Semiconductor material and device characterization[END_REF]. Diffusion capacitance in an n-p + (pn + ) junction deals with hole diffusion (electron diffusion) into the n-doped (p-doped) QNR. By expressing all these current densities in a p-n + junction, we can rewrite (6) as:

J tot = Qn(t) τn + V 1 (t) R sh + C SCR (t) dV 1 (t) dt + dQn(t) dt (7)
with:

Qn(t) = Q n0 exp V 1 (t) ηV th -1 (8)
Q n0 = q n 2 i N A √ Dnτn ( 9
)
where Qn is the electron density into p-doped QNR in C cm -2 , Q n0 the electron density into p-doped QNR at thermal equilibrium in C cm -2 , τn the electron lifetime in s, V 1 the OCVD voltage in V, R sh the shunt resistance in Ω cm 2 and C SCR the SCR capacitance in F cm -2 . N A is the acceptor doping level. It corresponds also to the doping level of the lowest-doped region, so N A = N l . As the p-n junction is asymmetric we assumed τn is equal to τ ef f [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF]. During OCVD measurement, it is assumed that J tot = 0 because of the opencircuit condition. Then, simulation of the V 1 decay is performed by solving the following 1 st order ordinary differential equation (ODE):

0 = Q n0 τn exp V 1 (t) ηV th -1 + V 1 (t) R sh + dV 1 (t) dt Q n0 ηV th exp V 1 (t) ηV th + C SCR (t) (10)
We used the function "odeint" [24] provided by Python-SciPy library called "scipy.integrate" to solve the ODE. The resolution method is Adams/BDF [25].

Comparison to Mahan and Barnes (M&B) model

First modelling results have been obtained using M&B parameters [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF]. They are presented in the following Table :   Table 2 Parameters of M&B models [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF] at 300 K. signal without major leakage current into the SCR and so the signal is solely governed by diffusion and SCR capacitances. C SCR decreases decay rate into region (iii). On the other hand, the signal increases it abruptly without major SCR capacitance effects but with major leakage current into the SCR. When SCR capacitance and shunt resistance are omitted, the signal is only due to diffusion capacitance and it is linear over the entire decay. Therefore, region (iii) is logically vanished. OCVD decay without diffusion capacitance represents the case when lifetime is so small that region (ii) is completely fused with region (iii) [START_REF] Lemaire | Mesure par OCVD de la durée de vie des porteurs minoritaires dans des jonctions en GaSb, en GaAs et en Si : simulations et expérimentations[END_REF]. These results are consistent regarding the theory and similar to those obtained by [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF].

D N l n i r V bi τ ef f [cm 2 V -1 s -1 ] [cm -3 ] [cm -3 ] [-] [V] [µs]
t/τ ef f (-) V 1 (V) w/o R sh w/o C SCR w/o R sh and C SCR w/o C D
As we previously explained, the work of M&B didn't consider the dynamics of SCR capacitance. Then, we compared our model to their one-diode model with two different applied voltages Va (0.9 V and 0.6 V) below the built-in voltage V bi (0.93 V) in Fig. 4

(ii) (iii) V a = 0.6 V V a = 0.9 V t (µs) V 1 (V) C SCR (
10 -5 Vii-iii V bi Voltage (V) Capacitance (F cm -2 ) C D C SCR (V(t)) C SCR (0.6 V)
C SCR (0.9 V) Fig. 5: Diffusion and SCR capacitances over the voltage according to M&B parameters [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF].

varying influences on region (iii) depending on Va. The dynamics of SCR capacitance has a greater influence when Va is higher. Indeed, the closer Va from V bi , the higher the SCR capacitance. This is clearly shown in Fig. 5 where diffusion and SCR capacitance evolutions over the voltage are presented.

V ii-iii represented in Fig. 5 is the voltage for which SCR capacitance is equal to diffusion capacitance. It can be seen as the transition voltage between region (ii) and region (iii). The lower the voltage, the lower the slope of the voltage-dependent SCR capacitance C SCR (V (t)). Its most important variation occurs around the built-in voltage that is, when SCR is really thin.

The difference appearing Fig. 4 in region (iii) is the motivation for taking the dynamics of SCR capacitance into account since we wanted to extract accurately N l , V bi and R sh from this region.

Parametric study

This study allows to better understand OCVD behaviour. We studied the influences of the effective lifetime (τ ef f ), doping level of the lowest-doped region (N l ), shunt resistance (R sh ) and built-in voltage (V bi ) on the OCVD signal of a silicon p-n + junction. The default parameters of the model are reported in Table 3. N h is the doping Table 3 Default parameters of silicon at 300 K and for Va = 0.8 V.

N l N h n i r R sh τ ef f [cm -3 ] [cm -3 ] [cm -3 ] [-] [Ω cm 2 ] [µs] 1 × 10 15 1 × 10 19 1 × 10 10 11.7 1 × 10 19 1
level of the highest-doped region and r the relative electrical permittivity of silicon. All simulations were performed for a p-n + junction with a ratio N h /N l > 100. The built-in voltage V bi was calculated from (4) and shunt resistance was overestimated on purpose to omit its influence. The diffusion coefficient Dn was calculated through electron mobility µn in cm 2 V -1 s -1 using the following equation:

Dn = k b T q µn (11)
Electrons mobility was calculated according to the doping-dependent model of Masetti [START_REF] Masetti | Modeling of carrier mobility against carrier concentration in arsenic-, phosphorus-, and boron-doped silicon[END_REF] expressed here:

µ M = µ 1 exp -Pc N + µc -µ 2 1 + N Cr α - µ 3 1 + Cs N β (12)
The parameters for [START_REF] Green | Minority carrier lifetimes using compensated differental open circuit voltage decay[END_REF] are provided Table 4. electrons holes

µc [cm 2 V -1 s -1 ] 1417 470.5 µ 1 [cm 2 V -1 s -1 ] 52.2 44.9 µ 2 [cm 2 V -1 s -1 ] 52.2 0 µ 3 [cm 2 V -1 s -1 ]
43. Simulated OCVD signal for the various lifetimes τ ef f is presented in Fig. 6. We normalized x axis by τ ef f in order to plot all the results on the same figure. We know from ( 7) and ( 2) that diode current and diffusion capacitance are lifetime dependent. Thus, the slope of region (ii) would be as low as the lifetime is high according to [START_REF]Handbook of photovoltaic science and engineering[END_REF]. Because of the x axis normalization, the three slopes are the same. On the other hand, the most important influence is related to V ii-iii : the lower the lifetime, the higher V ii-iii . Thus for lower lifetime, duration of region (ii) is reduced for a given doping level N l and a given applied voltage Va. This can be explained by the reduction of diffusion capacitance with lower lifetime [START_REF] Schroder | Semiconductor material and device characterization[END_REF]. If the diffusion capacitance is lower, the SCR capacitance influence appears at higher voltage and so earlier, reducing region (ii) duration. It means that the lifetime extraction accuracy is reduced with decreasing lifetime. Fig. 7 represents simulated OCVD signals for several doping levels N l . We observe that V ii-iii increases either with increasing N l or with decreasing τ ef f we saw previously. According to (5) and (2), N l which represents majority carriers density, influences both capacitances (diffusion and SCR). Increasing N l reduces the minority carriers density (i.e. decreases the diffusion capacitance) with a dependence to 1/N l on the one hand, and reduces the SCR thickness (i.e. increases the SCR capacitance) with a dependence to √ N l on the other hand. These double and imbalanced influences makes the effect of N l stronger on the signal, than its of τ ef f . Indeed, a change of τ ef f by two orders of magnitude has less influence than one oder of magnitude of N l . Therefore, the lifetime extraction accuracy is also reduced as well with increasing N l since this reduces the region (ii) duration.

(ii) (iii) t/τ (-) V 1 (V) τ ef f = 10 ns τ ef f = 1 µs τ ef f = 100 µs
(iii)

t (µs) V 1 (V) N l = 1 × 10 15 cm -3 N l = 1 × 10 16 cm -3 N l = 1 × 10 17 cm -3
Fig. 7: Simulated OCVD signals for various doping levels N l

Shunt resistance:

To the best of our knowledge, the effect of the shunt resistance on the OCVD signal has never been investigated jointly with the presence of a SCR capacitance. This investigation has been conducted here thanks to modelling (see Fig. 8). First, as M&B have previously demonstrated [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF], we observe that V 1 decay rate increases with decreasing shunt resistance, i.e. with increasing leakage current. Indeed, shunt resistance is a bypass for charge carriers stored into QNR (diffusion capacitance) and SCR which makes recombinations faster (see Fig. 8). OCVD signal for R sh = 1 × 10 19 Ω cm 2 can be considered free from shunt resistance influence because of its high value. However, region (iii) is influenced for R sh = 1 × 10 4 Ω cm 2 . On the other hand, region (ii) and region (iii) are both affected for R sh = 1 × 10 2 Ω cm 2 and R sh = 1 × 10 3 Ω cm 2 . Indeed, the shunt resistance decrease and the SCR capacitance increase, have opposite influences on the OCVD signal. Thus region (ii) can also be affected as soon as the diode instantaneous resistance is roughly equal to the shunt resistance. Therefore, its duration may be extended when both effects balance each other. This is shown in Fig. 8 for R sh = 1 × 10 3 Ω cm 2 . In this specific case, the lifetime extraction accuracy is enhanced. This is precisely from that phenomenon that Green has created his alternative circuit [START_REF] Green | Minority carrier lifetimes using compensated differental open circuit voltage decay[END_REF]. He wanted to increase the accuracy of region (ii) identification by cleverly playing with these antagonistic effects. But this required to find a suitable configuration for each characterized samples. Finally, when the shunt resistance

(ii) (iii) t (µs) V 1 (V) R sh = 1 × 10 2 Ω cm 2 R sh = 1 × 10 3 Ω cm 2 R sh = 1 × 10 4 Ω cm 2 R sh = 1 × 10 19 Ω cm 2
Fig. 8: Simulated OCVD signals for various shunt resistances R sh has a strong influence (R sh = 1 × 10 2 Ω cm 2 ), V ii-iii is shifted up and the lifetime extraction accuracy is reduced.

Built-in voltage:

The built-in voltage V bi depends on N l and N h for a given material according to [START_REF] Yablonovitch | A contactless minority lifetime probe of heterostructures, surfaces, interfaces and bulk wafers[END_REF]. Therefore, we varied N h from 3.4 × 10 18 cm -3 to 1.1 × 10 20 cm -3 in order to change V bi (from 0.805 V to 0.894 V) but we kept N l = 1 × 10 15 cm -3 in order to work with an asymmetric p-n junction where N h /N l > 1000. The variation of V bi does 9: Simulated OCVD signals for various built-in voltages V bi not significantly affect the OCVD signal (see Fig. 9). V bi defines the value of the voltage below which a SCR capacitance appears. However, since the diffusion capacitance is much higher than the SCR capacitance for this particular value of the lifetime, its influence remains negligible. It means that V bi extraction from region (iii) would be rather complicated if diffusion capacitance is too high.

(ii) (iii) t (µs) V 1 (V) V bi = 0.805 V V bi = 0.845 V V bi = 0.894 V Fig.

Results in silicon diode

This section is dedicated to the fitting procedure of an experimental OCVD signal for a silicon diode (Vishay 1N4007). The aim is to validate the extraction of τ ef f , N l , V bi and R sh and so, to highlight the flexibility of OCVD method. We assumed the diode is an asymmetric junction. We decided first to apply our model to a silicon diode because this semiconductor has well know properties on the one hand, and it has a long lifetime that guarantees clear (ii) and (iii) regions on the other hand. Regarding the silicon diode, we have no information about its design (thickness layers, doping levels and so on). Thus we extracted the experimentally required variables from dark I-V and C-V curves in order to compare them to those obtained by OCVD modelling.

Experimental details

OCVD characterization was carried out using the "series circuit" presented Fig. 10. The OCVD bench is composed of a 1 A current/voltage source (Keithley 2400), a 2 GHz bandwidth oscilloscope (Tektronix TDS784D), an n-type enhancement MOSFET (Semelab D2081UK) and a function generator (MB Electronique Wavetek FG3B) to drive the MOSFET. The MOSFET has been selected for its low commutation time (< 30 ns). As the design of the silicon diode (Vishay 1N4007) is unknown, we recorded a dark I-V characteristic using the Keithley 2400 and a C-V characteristic using a C-V plotter (EG&G Model 410) at 1 MHz. The dark I-V curve allows to determine the ideality factor and the shunt resistance R sh while the C-V curve allows to extract N l and V bi from an asymmetric p-n junction [START_REF] Schroder | Semiconductor material and device characterization[END_REF][START_REF] Van Opdorp | Evaluation of doping profiles from capacitance measurements[END_REF]. Both curves are presented in been calculated through the local ideality factor η(V ) [START_REF] Mcintosh | Lumps, Humps and Bumps: Three Detrimental Effects in the Current-Voltage Curve of Silicon Solar Cells[END_REF] with the following equation: capacitances are clearly distinguishable and similar to those obtained from theory presented in Fig. 5. It is possible to extract N l and V bi by fitting the C-V curve with linear regression of 1/C 2 SCR [START_REF] Schroder | Semiconductor material and device characterization[END_REF][START_REF] Van Opdorp | Evaluation of doping profiles from capacitance measurements[END_REF] given by:

η(V ) = q k B T dV d (ln I) (13
1 C 2 SCR = - 2V q N l A 2 + 2V bi q N l A 2 ( 14
)
where A is the area of the silicon diode in cm 2 . Fig. 12 shows results of 1/C 2 and its linear regression over the range [-0.7: 0 V]. We obtained N l = 2.19 × 10 12 cm -3 and V bi = 0.694 V for a diode area of 3.14 mm 2 . All extracted parameters have been reported in Table 5.

As we didn't have any information about the nature of the junction, we performed modelling on the p-n + junction and the n-p + junction. We determined diffusion coefficient D through [START_REF] Mahan | Measurement of minority carrier lifetime in solar cells from photo-induced open-circuit voltage decay[END_REF] and calculated mobility µ M from [START_REF] Green | Minority carrier lifetimes using compensated differental open circuit voltage decay[END_REF] for each junction assumption. We considered Voltage (V)

1/C 2 (F -2 )
experimental data linear regression Fig. 12: Squared inverse SCR capacitance of the silicon diode (Vishay 1N4007) over the voltage.

Table 5 Extracted parameters of the silicon diode (Vishay 1N4007) from dark I-V (η and R sh ) and C-V (N l and V bi ) characteristics.

η R sh N l V bi [-] [Ω cm 2 ] [cm -3 ] [V]
1.56 1.57 × 10 6 2.19 × 10 12 0.694 electrons (holes) mobility for a p-n + junction (n-p + junction). This feature of charge carriers in semiconductors opens the possibility for our model to identify the junction type (p-n + ou n-p + ). But it also leads to a lack of accuracy if we do not choose the right junction type.

Fitting process details

The main purpose of our model is to extract various parameters such as τ ef f , N l , V bi and R sh from an experimental curve. This would make the OCVD a flexible method of measurement for semiconductors. Therefore, the model requires an optimization algorithm to properly fit the experimental curve where fitting variables would be τ ef f , N l , V bi and R sh . We used an iterative optimization algorithm to minimize the following RMSE equation:

RM SE = 1 n n i=1 V mod (t i ) -Vexp(t i ) 2 (15)
where n is the number of measurements of the simulated signal, Vexp(t i ) the experimental OCVD voltage at time t i and V mod (t i ) the simulated OCVD voltage at time t i . OCVD voltage corresponds to V 1 in Fig. 2. RMSE is the Root Mean Square Error between Vexp and V mod over the entire time range. The optimization algorithm we used is provided by Python-SciPy library called "scipy.optimize" with the function "minimize" [START_REF] Scipy | minimize(method='Nelder-Mead') -SciPy v1.3.0 Reference Guide[END_REF]. We selected "Nelder-Mead" algorithm [30] which is the default one. We didn't work on algorithms during this work. Only 3 fitting variables were given to the algorithm for each simulation.

Extraction results

Presented below is an extraction of variables R sh , τ ef f , N l and V bi thanks to OCVD modelling. We compared these values obtained to those derived from the dark I-V and C-V characteristics. These values were taken as reference for comparison and are given in Tables 5 and6. Lifetime is estimated from OCVD region (ii) which is linear and long enough to consider reliability of region (ii) by using [START_REF]Handbook of photovoltaic science and engineering[END_REF]. We present plots for a p-n + junction only (see Fig. 13 and Fig. 14) First, we show results on Fig. 13 where τ ef f , N l and V bi are the fitting variables while R sh is kept constant. We took a large value of the shunt resistance (R sh = 1 × 10 9 Ω cm 2 ) in order to avoid its influence on OCVD signal. Initial values of variables were selected arbitrary but are consistent for silicon diodes: τ ef f = 7 µs, N l = 1 × 10 15 cm -3 and V bi = 0.75 V. We obtained RM SE = 0.715 % (the same for a n-p + junction). However there are differences concerning the extracted variables depending on the junction type. Indeed, column Si 1 of Table 6 shows that the extracted values of N l are slightly different for the two kinds of junction while τ ef f are the same. It means, τ ef f does not depend on the junction type. N l and τ ef f values are consistent (around 5 × 10 11 cm -3 and 9 µs respectively) while V bi is far too high (1.16 V). Actually, it can not exceed 0.738 V that is obtained from (4) for N l = 5 × 10 11 cm -3 and N h = 5 × 10 20 cm -3 (roughly the intrinsic limit before being an alloy). We knew that the extraction of V bi would be difficult due to its low influence on OCVD signal (see subsection 3.3.4). Finally, it is not possible to conclude on the junction type even if N l in the case of the p-n + junction is slightly closer from N l extracted by C-V.

We performed a new simulation in order to carry out a fit with a consistent value of V bi . Therefore, we fixed V bi = 0.694 V (value extracted by C-V). The three fitting variables were R sh , N l and τ ef f . The initial values were those extracted previously (given in column Si 1 of Table 6). We chose R sh as a new fitting variable because the experimental OCVD signal shows a smooth increase of the slope in region (ii) around 0.35 V (see Fig. 13). As discussed previously, the increase of the decay rate can be explained by a shunt resistance effect (see subsection 3.3.3). The result for a p-n + junction is presented in Fig. 14. The use of R sh as new fitting variable does not enable to reduce drastically the RMSE whatever the junction type. The extracted variables are presented in column Si 2 of Table 6. Extracted lifetime τ ef f does not change notably by varying R sh instead of V bi . Regarding N l , the extraction provides results slightly higher than those obtained previously: 9.11 × 10 11 cm -3 and 6.51 × 10 11 cm -3 for a p-n + junction and an n-p + junction respectively. R sh values are also not so different for both types of junction. The n-p + junction presents the closest shunt resistance value (4.24 × 10 5 Ω cm 2 ) from that obtained from the dark I-V curve (1.57 × 10 6 Ω cm 2 ). Nevertheless, the junction type is also complicated to define. Thus, these results show how important it is to know the nature of the junction, otherwise, the reliability of the method could be affected. Moreover, many semiconductors such as III-V, show more imbalanced mobility (between electrons and holes) as compared to silicon [START_REF] Levinshtein | Handbook series on semiconductor parameters[END_REF]. That would probably lead to increase the gap in N l between both junction types.

(ii) (iii) t (µs) V 1 (V)
The results of column Si 2 in Table 6 demonstrates that our model is able to properly fit the experimental OCVD curve of a silicon diode (RM SE < 0.57 %) and to extract rather suitable value of R sh , N l and τ ef f with the support of a dark I-V curve (i.e. to determine η).

Conclusion and discussion

The theoretical OCVD signal presents a region with a linear decay which enables extraction of minority carrier effective lifetime τ ef f in simple way. This region is often difficult to identify because of SCR capacitance and/or generation/recombination (G/R) current into the SCR. There is sometimes no linear region at all when lifetime is too short for instance. The quasi-steady stade mode [START_REF] Sinton | Quasi-steady-state photoconductance, a new method for solar cell material and device characterization[END_REF] has been applied to the OCVD (know as "Suns-Voc" [START_REF] Sinton | A QSS open circuit voltage method for solar cell characterization[END_REF]) in order to avoid detrimental SCR capacitance effect. Our approach has been to take advantage of this SCR capacitance effect since it also allows, theoretically, extraction of N l , V bi and R sh . Therefore, OCVD signal modelling is mandatory to overtake this difficulty. A first OCVD modelling has been carried out a few decades ago to explain influences of SCR capacitance and SCR G/R [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF]. M&B did not derive τ ef f or other parameters from their modelling, although the OCVD signal is full of information.

We developed a model under Python programming language in order to simulate the OCVD signal. Our simulation is based on a lumped element model of one asymmetric diode in transient regime. We took M&B work as reference to develop our model. We went further by taking into account the dynamic of the SCR capacitance on the one hand and on the other hand, the joint influences of the SCR capacitance and the SCR G/R. We used this model to study the influences of τ ef f , N l , R sh and V bi on the OCVD signal. We first concluded that the lower τ ef f and the higher N l , the lower the extraction accuracy. We also concluded that the low impact of V bi on OCVD signal makes its extraction unreliable.

Then OCVD signal fitting was carried out on an experimental curve of a silicon diode through an optimization algorithm provided by SciPy libraries. We reached RM SE < 0.57 % for both junction types (n-p + and p-n + ). We have been able to extract a reliable value of τ ef f (i.e. no variation whatever is the junction type). We obtained values of N l and R sh rather close from those obtained from C-V and dark I-V characteristics. The extracted N l values are rather similar for both junction types with a constant ratio N l(p-n + ) /N l(n-p + ) ≈ 1.4. This gap is probably due to the small difference between holes and electrons mobilities in silicon. However, this difference is usually higher in semiconductors such as III-V. It means this gap should be higher with III-V semiconductors impacting then the reliability of the extraction for unknown junction types. On the whole, the extracted values of N l and R sh were consistent but didn't allow to rule on the junction type.

We pointed out how important it is to know the right junction type in order to increase the reliability of N l and R sh extraction. By going further, this work could be complemented by assuming a symmetric junction through (3) or a different doping profile. Indeed, our model is based on the assumption of an abrupt junction. A non-abrupt junction would induce a different dynamics of the SCR capacitance over the time [START_REF] Chawla | Transition region capacitance of diffused p-n junctions[END_REF]. Then, the model could be improved by taking into account a second level of injection (two-diode model) in relevant case and different mobility models. Furthermore, it would be interesting to take into account the series resistance effect in order to simulate region (i). Finally, our model is only based on diode effects while it is clear that the MOSFET also affects the signal. Moreover, the circuit itself has different influences whether it is a "series" or "parallel" configuration. Finally, the method has been subjected to such a high lifetime junction which provides clear linear decay but this is not always the case. Therefore, today, we don't know if the flexibility of the method can be applied to OCVD signal having no linear decay but it opens the perspective to. We can clearly consider, the presented method is well suitable for Si solar cells since it has such a high lifetime, well known properties and an asymmetric design. However, its operation still has to be demonstrated for short lifetime junction whose clear linear decay is not guaranteed.

The first results presented here are promising and even open up the perspective for measuring short lifetime by OCVD with the support of a dark I-V curve and assuming a known design. Instead of measuring lifetime through the linear region (subjective method) it allows fast extraction by a fitting model (objective method). Since short lifetime usually does not provide linear region on OCVD signal, this work is of major importance to apply OCVD to semiconductors such as III-V. These materials are of major interest for PV applications such as multijunction solar cells. On the one hand, our approach makes the method less simple in terms of analysis, because semiconductor properties (µ, n i and r ) have to be known but it enhances greatly OCVD capabilities on the other hand. Furthermore, this is the first time, to the best of our knowledge, OCVD method demonstrates its measurement flexibility. It is a low-cost and a low-tech method so it represents an important breakthrough to achieve widespread use.

Fig. 1 :

 1 Fig. 1: Typical OCVD signal of a silicon diode under moderate injection obtained by "series circuit".

20 1 × 1 Fig. 3

 113 Fig.3shows four different OCVD signals. First, we note that there is no region (i) because no series resistance has been taken into account. The curve without (w/o) shunt resistance represents OCVD

Fig. 3 :

 3 Fig. 3: OCVD signals simulated according to M&B parameters [7].

  . For the model of M&B, we calculated the SCR capacitance at Va. Results shows

Fig. 6 : 3 . 3 . 2

 6332 Fig. 6: Simulated OCVD signals for various lifetimes τ ef f

Fig. 11 . 1 VaFig. 10 :

 11110 Fig. 10: OCVD "series" circuit used for silicon diode characterization.

Fig. 11 :

 11 Fig. 11: Dark I-V and C-V characteristics of the silicon diode (Vishay 1N4007).

Fig. 14 #

 14 fitting variables with random and consistent initial values of variables * fitting variables with Si 1 results as initial values of variables because the simulated signals are about the same for both p-n + and n-p + junctions.

dV 1 /Fig. 13 :

 113 Fig.13: OCVD signal: experiment and modelling. Modelling signal has been obtained with τ ef f , N l and V bi as fitting variables.

Fig. 14 :

 14 Fig. 14: OCVD signal: experiment and modelling. Modelling signal has been obtained with τ ef f , N l and R sh as fitting variables.

Table 1

 1 Comparison between M&B models[START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF] and ours.

		M&B #1	M&B #2	Lemaire et al.
	Assignement of C SCR effect	Capacitance (constant)		Capacitance (voltage-dependent)
	Assignement of G/R current effect		Diode (voltage-dependent)	Resistance (constant)
	Model	one-diode	two-diode	one-diode
	A cm -2 :

Table 4

 4 Parameters of Masetti mobility model in silicon at 300 K[START_REF] Masetti | Modeling of carrier mobility against carrier concentration in arsenic-, phosphorus-, and boron-doped silicon[END_REF].

Table 6

 6 Extracted parameters of the silicon diode (Vishay 1N4007) from OCVD fitting and I-V/C-V characteristics. N l [cm -3 ] 2.19 × 10 12 5.76 × 10 11# 4.11 × 10 11# 9.11 × 10 11 * 6.51 × 10 11 *

	dark I-V and C-V	p-n +	Si 1	n-p +	p-n +	Si 2	n-p +
	RM SE [%] τ ef f [µs] 8.78 V bi [V] 0.715 0.569 --8.93 # 0.694 1.159 # 0.694 1.57 × 10 6 1 × 10 9 3.58 × 10 5 * 4.24 × 10 5 * R sh [Ω cm 2 ] η [-] 1.56 1.56 1.56 n [-] -20 000 20 000
	Figure	Fig. 13

*

This paper is a preprint of a paper accepted by Circuits, Devices and Systems journal and is subject to Institution of Engineering and Technology Copyright. When the final version is published, the copy of record will be available at IET Digital Library"

This paper is a preprint of a paper accepted by Circuits, Devices and Systems journal and is subject to Institution of Engineering and

Technology Copyright. When the final version is published, the copy of record will be available at IET Digital Library"

Acknowledgments

Antoine Lemaire would like to thank Clément Lacroix, Hanany Tolba and Julien Nou for their support in Python programming and advices for the optimization algorithm. He is also grateful to Thomas Pettigrew.


domain: spi info math phys

text:
Léon Reboul

Michel Kieffer

Hélène Piet-Lahanier

Sébastien Reynaud

Cooperative guidance of a fleet of UAVs for multi-target discovery and tracking in presence of obstacles using a set membership approach

Keywords: Cooperative guidance, membership-set estimation, target search and tracking

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Cooperative guidance of a fleet of UAVs for multi-target discovery and tracking in presence of obstacles using a set membership approach

Léon Reboul * Michel Kieffer * * Hélène Piet-Lahanier * Sébastien Reynaud * Observation of multiple targets by a fleet of cooperative UAVs constitutes a challenging task involved in a wide range of applications including crowd monitoring, search operations, surveillance, or tactical patrolling. The determination of the UAV trajectories must fulfill the two following requirements. All the targets contained within the search area must be detected and the detected target trajectories must be tracked with sufficient precision. These usually time-critical missions may be more efficiently fulfilled when the agents are cooperating. Fleets of UAVs have already demonstrated their flexibility to mission requirements, robustness to faults, and ability to achieve the global objective in a shorter time.

In this paper, our aim is to control a fleet of UAVs so as to search and track an a priori unknown number of targets moving in an urban-type environment. Most of the methods presented in the surveys Robin and Lacroix (2016); [START_REF] Khan | Cooperative robots to observe moving targets[END_REF], addressed either the target detection or the target tracking problem, but seldom tackled simultaneously the two aspects. Examples of methods addressing both aspects can be found, e.g., in [START_REF] Adamey | Collaborative multi-msa multi-target tracking and surveillance: a divide & conquer method using region allocation trees[END_REF]; [START_REF] Dames | Distributed multi-target search and tracking using the phd filter[END_REF]. In [START_REF] Adamey | Collaborative multi-msa multi-target tracking and surveillance: a divide & conquer method using region allocation trees[END_REF], a grid in the configuration space is considered. A Bayesian filter updates the probability attached to each grid point to be a target location. An allocation of the UAVs to the region of interest is performed using a binary tree search. The approach in [START_REF] Dames | Distributed multi-target search and tracking using the phd filter[END_REF] defines search patterns using Lloyd's algorithm, combined with a probability hypothesis density filter to evaluate the occurrence probability of

This work has been performed with DGA support.

targets, allowed to enter or leave the search area. Target discovery and tracking is performed combining information from sensors mounted on several UAVs. This problem is formulated as a path planning problem to determine the future UAV sensor locations that are reachable and decrease the target state estimation uncertainty. Various methods have been considered in this context, see, e.g., [START_REF] Yu | Cooperative path planning for target tracking in urban environments using unmanned air and ground vehicles[END_REF]; [START_REF] Capitan | Decentralized multi-robot cooperation with auctioned pomdps[END_REF]; [START_REF] Morbidi | Active target tracking and cooperative localization for teams of aerial vehicles[END_REF].

In all of the aforementioned approaches, a stochastic framework is used to describe the measurement errors and model uncertainties. As pointed out in [START_REF] Gu | Active persistent localization of a three-dimensional moving target under set-membership uncertainty description through cooperation of multiple mobile robots[END_REF], the resulting performance may prove sensitive to the a priori assumptions on the probability density functions (pdfs) describing the process and measurement noises.

As the definition of suitable pdfs might prove tedious, a set-membership description of uncertainties is suggested in [START_REF] Gu | Active persistent localization of a three-dimensional moving target under set-membership uncertainty description through cooperation of multiple mobile robots[END_REF]. The only assumption made on the noises and uncertainties is that their realizations remain within known bounds. Using this description, one no longer searches for a single point estimate associated with a posterior density function but for sets guaranteed to contain the target states at each time step.

In this paper, we assume that each UAV is equipped with a sensor able to detect and localize targets in some compact subset of the search area. A distributed set-membership estimator has been presented in [START_REF] Reynaud | A set-membership approach to find and track multiple targets using a fleet of uavs[END_REF] that accounts for information provided by each UAV, but also for the absence of detected target in the sensed subset. Compared to [START_REF] Reynaud | A set-membership approach to find and track multiple targets using a fleet of uavs[END_REF], this paper addresses the target discovery and tracking problem in urban-type environment, where buildings and obstacles are described by convex polytopes. The information provided by UAV sensors has thus to account for visibility issues. A control input for each UAV is derived by predicting the impact of future measurements on the set estimates of target states. The control inputs minimizing the estimation uncertainty are selected.

The paper is organized as follows. Section 1 describes the multi-target multi-UAV detection and tracking problem.

In Section 2, the set membership estimator introduced in Reynaud et al. ( 2018) is recalled. Section 3 describes the way obstacles modify the field of view of UAVs. Section 4 presents a criterion to evaluate the impact of future measurements on target state estimates and the cooperative guidance scheme that defines for each UAV the trajectory to follow so as to minimize an upper bound of the estimation uncertainty. In Section 5, the control input design algorithms are evaluated on simulations. Section 6 concludes this paper.

PROBLEM FORMULATION

We consider N u identical UAVs flying over some area of interest, assumed to be a compact subset of a plane. Buildings and obstacles are present in the environment of the UAVs. They are represented by a set of convex polytopes P , = 1, . . . , N o , each of which is defined by its set of vertices V = {z ,1 , . . . , z ,n } , where z ,m ∈ R 3 . N t targets are potentially moving within the area of interest. They are assumed not being able to enter the obstacles. N t is fixed but not known a priori. Each UAV is equipped with a sensor providing observations of some subset of the area of interest. The objective of the fleet of UAVs is to detect and track the targets using these observations.

The dynamics of targets and UAVs are described by discrete-time dynamical systems of the form

x T j,k+1 = g T k x T j,k , v j,k (1)
x

U i,k+1 = f U k x U i,k , u i,k (2)
where x T j,k ∈ R nT and x U i,k ∈ R nU are the state vectors of the jth target and the ith UAV, respectively. Time is sampled with a constant period T and k represents the discrete time instant (t = kT ). The inputs v j,k are unknown target state perturbations belonging to some known box common to all targets [v k ] and u i,k is the control input for UAV i, whose values belong to the set U of admissible control inputs. At time k = 0, all initial target states x T j,0 are only assumed to belong to some a priori known compact set Z 0 ⊂ R nT .

Representation of information

While evolving along their trajectories, at each time instant k, the UAVs observe subsets

F i (x U i,k ) ⊂ R nT , i = 1, . . . , N u of the area of interest. F i (x U i,k ) is called the field of view (FoV) of UAV i at time k. More details about F i (x U i,k
) are provided in Section 3. Let L i,k be the list of the targets detected by UAV i at time k. Assuming that the probabilities of non-detection and of false alarm are null, one has

x T j,k ∈ F i x U i,k ⇔ j ∈ L i,k . (3)
If the jth target state belongs to F i (x U i,k ), the ith UAV collects a noisy observation of x T j,k described as

y i,j,k = h i x U i,k , x T j,k + w i,j,k (4
) where h i is the observation equation and w i,j,k is the measurement noise, bounded in the known box [w k ].

Each UAV i maintains a set D i,k of indices of targets already detected at time k. The availability of F i (x U i,k ) and of the detection or non-detection information are used to evaluate a set of set estimates Z i,k = {Z i,j,k } j∈D i,k containing the state values of the already detected targets. Each set Z i,j,k has to contain all possible values of x T j,k that are consistent with the information available to UAV i up to time k. Moreover, UAV i also maintains a set Z i,k containing the possible state values of not detected targets.

Target state estimation uncertainty

The target state estimation uncertainty can be evaluated using the previously introduced set estimates. The choice considered in [START_REF] Reynaud | A set-membership approach to find and track multiple targets using a fleet of uavs[END_REF] is

Φ Z i,k , Z i,k = 1 max {1, |D i,k |} j∈D i,k φ (Z i,j,k ) + αφ Z i,k where φ (Z i,j,k ) is the volume of the set Z i,j,k , |D i,k | is the cardinal number of D i,k
, and α some parameter to adjust the relative importance of the state estimation uncertainties of detected and of not yet detected targets. The (target state) estimation uncertainty at time k is defined as the average of Φ

Z i,k , Z i,k over all UAVs Φ k = 1 N u Nu i=1 Φ Z i,k , Z i,k .
(5)

Communications

Assume that the topology of the fleet is described by an

undirected graph G = (N U , E). N U = {1, 2, ..., N U } is the set of nodes and E ⊂ N U × N U the set of edges of the network. The set of neighbors of UAV i is N i = {j ∈ N U | (i, j) ∈ E, i = j}.
UAV i may exchange information with all its neighbors, which allows each of the neighboring UAVs to update their set estimates.

SET ESTIMATION

At time k, the information available at UAV i consists in the sets

D i,k , Z i,k = {Z i,j,k } j∈D i,k and Z i,k .
The way these sets are updated at time k + 1 considering additional information available has been detailed in [START_REF] Reynaud | A set-membership approach to find and track multiple targets using a fleet of uavs[END_REF]. The update rules alternating a prediction step followed by a correction step are briefly recalled here.

Prediction step

The prediction step is only used under the hypothesis of moving targets. The predicted set of detected targets is simply evaluated as D i,k+1|k = D i,k . For each target j ∈ D i,k+1|k , the set of possible target state values at time k + 1 consistent with Z i,j,k , with the dynamics (1), and the bounded state perturbation is

Z i,j,k+1|k = g Z k (Z i,j,k , [v k ]) , for all j ∈ D i,k+1|k (6)
Similarly, the predicted set Z i,k+1|k has to contain all possible state values of targets not already detected and evolving with dynamics (1)

Z i,k+1|k = g Z k Z i,k , [v k ] . (7)

Correction step from measurements

Assume that at time k + 1, UAV i has evaluated L i,k+1 from the observation of F i (x U i,k+1 ) and, for each j ∈ L i,k+1 , has access to y i,j,k+1 obtained using (4). The set of detected targets is then updated as D i,k+1|k+1 = D i,k+1|k ∪ L i,k+1 . Then, three cases have to be considered.

If Target j has been previously detected, i.e., if j ∈ D i,k and if x T j,k+1 ∈ F i (x U i,k+1
), only Z i,j,k+1|k is updated using the new measurement y i,j,k+1 and the measurement noise bound [w k+1 ] by

Z i,j,k+|k+1 = Z i,j,k+1|k ∩ h -1 i,k+1 (y i,j,k+1 -[w k+1 ]) , (8) where h -1 i,k+1 (y i,j,k+1 -[w k ]) is the pre-image of the box y i,j,k+1 -[w k ] by the function h i,k+1 (z) = h k+1 (r i,k+1 , z). If Target j has just been detected at time k + 1, i.e., if j / ∈ D i,k , but x T j,k+1 ∈ F i (x U i,k+1
), one knows that x T j,k+1 ∈ Z i,k+1|k and that x T j,k+1 has to be consistent with Z i,k+1|k , y i,j,k+1 , the measurement equation ( 4), and the measurement noise bound [w k+1 ], thus

Z i,j,k+1|k+1 = Z i,k+1|k ∩ h -1 i,k+1 (y i,j,k+1 -[w k+1 ]) . (9)
If Target j has been already detected, but is not detected at time k + 1, i.e., j

∈ D i,k , but x T j,k+1 / ∈ F i (x U i,k+1 ), Z i,j,k+1|k+1 = Z i,j,k+1|k \ F i (r i,k+1 ) , ( 10
) where B \ A = {x ∈ B | x ∈ A}.
To evaluate the set containing the state of targets still to be detected, one has to account for the fact that all targets in F i (x U i,k+1 ) have been processed. Thus, one has

Z i,k+1|k+1 = Z i,k+1|k \ F i x U i,k+1 . (11)

Correction step after communications

During each time step, UAV i communicates to its neighbors the sets D i,k+1|k+1 , Z i,k+1|k+1 , and Z i,k+1|k+1 and receives the corresponding sets from its neighbors. Using this additional information, UAV i updates its estimates.

The set of detected targets is evaluated taking the unions of all detected targets

D i,k+1 = D i,k+1|k+1 ∈Ni D ,k+1|k+1 . (12)
Consider the subset N j i,k+1 of neighbors of UAV i which have already detected target j at time k +1 and the subset N j i,k+1 of neighbors which have not detected target j. For all targets j which have been already detected by UAV i, the update equation is

Z i,j,k+1 = Z i,j,k+1|k+1 ∈N j i,k+1 Z ,j,k+1|k+1 ∈N j i,k+1 Z ,k+1|k+1 .
For all targets j which have not yet been detected by UAV i, the update equation is

Z i,j,k+1 = Z i,k+1|k+1 ∈N j i,k+1 Z ,j,k+1|k+1 ∈N j i,k+1 Z ,k+1|k+1 .
The subset of the state space still to be explored becomes

Z i,k+1 = Z i,k+1|k+1 ∈Ni Z ,k+1|k+1 . (13)

EVALUATION OF THE FIELD OF VIEW

In the frame F i attached to the ith UAV, one assumes that the FoV is defined as a conical combination of N f vectors x F i,1 , . . . , x F i,N f . When expressed in the state space of the UAV, in absence of obstacles, the FoV is

F i (x U i,k ) = α i,1 x F i,1 x U i,k + • • • + α i,N f x F i,N f x U i,k | α i,1 0, . . . , α i,N f 0} , (14)
where

x F i,n (x U i,k
) are the coordinates of x F i,n expressed in the global frame.

The area of interest contains obstacles. Due to occlusions, these obstacles will limit the FoV F i (x U i,k ) of UAV i. To simplify the analysis, UAVs are assumed to be represented by points in R nU and one only evaluates the impact of obstacles on the intersection of the FoV with the area of interest, which is also a compact subset of a plane. A point z in the area of interest will be visible by UAV i if and only if it belongs to the intersection of F i (x U i,k ) with the plane containing the area of interest, and the segment linking z and x U i,k does not intersect any obstacle polytope. If the segment linking z and x U i,k intersects the th polytope P , it can be seen to belong to the shadow of P , assuming that UAV i behaves as a point light source.

The impact of obstacles on the FoV requires thus to evaluate the shadow on the area of interest created by the obstacles when illuminated by UAV i. For that purpose, one adopts the approach introduced in Blinn (1988). The shadowed areas have then to be removed from the intersection of F i (x U i,k ) with the area of interest. Consider the set of vertices V = {z ,1 , . . . , z ,n } of the th obstacle P . Let z U i,k denote the projection of the state x U i,k of the ith UAV in the space R 3 of obstacles. Since P is convex, the shadow S (x U i,k ) on the area of interest of P for an UAV with state x U i,k is the convex hull of the intersections of the lines (z U i,k , z ,m ), m = 1, . . . , n with the plane containing the area of interest. These lines are defined as

z = z ,m + α z U i,k -z ,m , (15)
with α ∈ R. Let (z 1 , z 2 , z 3 ) T denote the components of the vector z. Assuming that the area of interest belongs to the plane defined by z 3 = 0, the intersection of (z U i,k , z ,m ) with this plane is obtained when α = -z ,m,3 /(z ,m,3 -z U i,k,3 ) in (15). The intersection F z3=0 i (x U i,k ) of the FoV with the area of interest, accounting for all obstacles is then obtained as

F z3=0 i x U i,k = F z3=0 i x U i,k \ No =1 F z3=0 i x U i,k ∩ S x U i,k
.

COOPERATIVE CONTROL DESIGN

Our aim is to determine for each UAV i, the sequence of control inputs that minimizes the predicted estimation uncertainty knowing the sets D i,k , Z i,k , and Z i,k . This determination is performed adopting the distributed Model Predictive Control (MPC) formalism introduced, e.g., in [START_REF] Morari | Model predictive control: past, present and future[END_REF]; [START_REF] Rochefort | Model predictive control of cooperative vehicles using systematic search approach[END_REF].

At time k, once the communications have been performed with its neighbors, each UAV i has access to D i,k , Z i,k , and Z i,k . If T h = hT denotes the time horizon of the MPC scheme, the sequence of control inputs should be designed for each UAV so as to minimize

Φ k+h = 1 N u Nu i=1 Φ Z i,k+h , Z i,k+h . (16)
Nevertheless, this is very difficult to evaluate in a decentralized manner. Several simplifications will thus be considered. First, one assumes that the communication graph does not change over the MPC time horizon. Second, UAV i assumes that the estimates at time k + κ of all its neighbors in N i are equal to D i,k+κ , Z i,k+κ , and Z i,k+κ , κ = 0, . . . , h. Then, instead of minimizing ( 16), UAV i will try to design its control inputs so as to minimize

Φ i,k+h = 1 N i + 1 ∈N i∪{i} Φ Z ,k+h , Z ,k+h = Φ Z i,k+h , Z i,k+h , (17)
where N i = N i ∪ {i}, since the estimates have been assumed to be equal. The design of the control inputs for UAV i has to account for the control inputs that will be evaluated by its neighbors. If the control inputs evaluated by UAV i at time k + κ are denoted u ,k+κ , κ = 1, . . . , h, ∈ N i , then the problem to be solved by UAV i is ( u ,k+1 , . . . , u ,k+h

) ∈N i = arg min Φ Z i,k+h , Z i,k+h , (18)
where the minimization is performed with respect to (u ,k+1 , . . . , u ,k+h ) for all ∈ N i .

Prediction of Φ

The solution of ( 18) requires evaluations of future values Φ i,k+κ of Φ i,k , knowing only D i,k , Z i,k , and Z i,k . As the potential locations of undetected targets are unknown, this is still very difficult. A bound on Φ i,k+1 has been introduced in [START_REF] Reynaud | A set-membership approach to find and track multiple targets using a fleet of uavs[END_REF], which depends only on the prediction of sets Z i,k and Z i,k from 6 and 7, and using the dynamics of the UAVs. This bound is adapted here, considering further approximations, to get approximations of Φ i,k+h in a recursive way.

At time k, determining whether future control inputs will help to identify new targets is very difficult. Thus, one assumes that the predicted values D P i,k+κ of D i,k+κ satisfy

D P i,k+κ = D i,k , κ = 1, . . . , h.
Assume that predicted values Z P i,k+κ-1 and Z P i,k+κ-1 of Z i,k+κ-1 and Z i,k+κ-1 have been evaluated, that the predicted states x U,P ,k+κ-1 of the UAVs in N i are known at time k + κ -1. Then, considering the control inputs u ,k+κ , ∈ N i , one is able to evaluate x U,P ,k+κ using x U,P ,k+κ-1 , u ,k+κ , and (2) for all ∈ N i . The evaluation of Z P i,k+κ|k+κ-1 requires considering a prediction step similar to that in (7) starting from Z P i,k+κ-1 . One assumes to simplify that Z P i,k+κ|k+κ-1 = Z P i,k+κ-1 , neglecting in the MPC approach the potential increase of the uncertainty related to the targets still to be detected.

Then Z i,k+κ can be easily predicted combining (11) and (13) as

Z P i,k+κ = Z P i,k+κ|k+κ-1 \ ∈N i F x U,P ,k+κ
and consequently,

φ Z P i,k+κ = φ   Z P i,k+κ-1 \ ∈N i F x U,P ,k+κ   . (19)
The sets Z P i,j,k+κ are much more complex to evaluate from Z P i,j,k+κ-1 , since they will depend on the state of the UAVs and the actual state of Target j. Nevertheless, Z P i,j,k+κ|k+κ-1 can again be predicted from Z P i,j,k+κ-1 using (7). A coarse approximation then is to consider Z P i,j,k+κ = Z P i,j,k+κ|k+κ-1 , neglecting information which may be provided by UAV measurements. With this approximation,

φ Z P i,j,k+κ = φ g Z k+κ-1 Z P i,j,k+κ-1 , [v k+κ-1 ] , (20
) which may be further approximated as

φ Z P i,j,k+κ = (1 + µ) φ Z P i,j,k+κ-1 , (21)
where µ depends on g Z k and [v k ], but can be taken as a constant. The main advantage of (21) compared to 20 is that it is not necessary to compute g Z k+κ-1 (Z P i,j,k+κ-1 , [v k+κ-1 ]), only the measures φ(Z P i,j,k+κ-1 ) of the sets Z P i,j,k+κ-1 are evaluated.

Using ( 19) and ( 21), one is able to evaluate recursively (17) for a given sequence of control inputs (u ,k+1 , . . . , u ,k+h ) for all ∈ N i , and to solve (18) considering any local minimization technique. As a result, it could be possible to obtain hippodrome type trajectories with UAVs circling above a detected target. To spur further explorations of the zone containing undetected targets, the criterion ( 18) is supplemented with

Analysis of uncertainty criterion

Φ add = ∈Ni d max x U , Z
where d max is the maximum distance between the current positions of the UAVs and points in the set Z . This term allows to pull the trajectories out of the area corresponding to local minima.

SIMULATIONS

Two scenarii with static and moving targets are considered. ImageSp and SIVIA have been implemented using the Intlab library [START_REF] Rump | Intlab-interval laboratory[END_REF].

Scenario with static targets

The targets are located in a plane, z j,k ∈ R 3 , with z j,k,3 = 0, j = 1, . . . , 6. The initial search area is a square of 500 × 500 m 2 . Three UAVs fly at a constant height of 100 m above the terrain at a constant speed of 15ms -1 . Their control inputs consist in the heading angle. The UAVs are assumed to communicate at each time instant. They are equipped with an optical sensor able to detect targets within its FoV. Its opening angles are equal to π 4 in both azimuth and elevation. When a target is detected at time k, one assumes that the measurement equation provides its actual location with an uncertainty bounded in [-10 m, 10 m] for both components.

Figure 3 presents the evolution of Z from k = 1 (top left) to k = 100 (bottom left). Figure 3 also presents the corresponding Z at k = 100 (bottom right), when all targets have been detected. Figure 4 presents the resulting trajectories of the three UAVs.

Scenario with moving targets

The targets are now assumed to move with a maximum speed of V T max = 5m.s -1 . The other parameters are the same as those considered in the previous section. Figure 5 presents the evolution of Z k for k = 10 (top left) and k = 15 (top right) and of Z at k = 15 (bottom right) and k = 40 (bottom left). One observes that due to the evolution of the targets, the size of the components of Z increases with k, especially when the targets are not redetected. Videos presenting the evolution of Z k and Z with time are available on https://drive.google.com/drive/folders/ 10EzWYZVGu-8J2qhcv8XkcL1-P97ZpEgj?usp=sharing. The trajectories of the three UAVs are presented on Figure 6. The resulting performances are evaluated using the decrease of uncertainty function Φ k . Figure 7 presents the evolution of the measures of Φ k when the targets are static (left), and when the targets are moving for respectively 3 and 4 UAVs(right).

In the static case, the global uncertainty monotonically decreases. When the targets are moving, Φ k appears to tend to a limit which is a compromise between the decrease of Z k due to the UAVs search and the growths of Z k and Z k during the prediction step to account for the displacement of the detected and not yet detected targets with the only knowledge of their maximal speed. As illustrated in Figure 7 (right), the localization uncertainty depends on the number of UAVs involved. This evaluation could provide a mean to determine the lower limit of the number of UAVs required to guarantee a certain level of uncertainty, given the size of the search zone, the characteristics of the UAVs in terms of speed, the size of the FoV, and the bound on the target speed.

CONCLUSIONS

A cooperative distributed guidance scheme has been presented for the coordinated control of a fleet of UAVs aiming to search and track an a priori unknown number of targets spread over a delimited geographical area. The trajectories of the UAVs are evaluated so as to minimize a criterion combining a measure of the uncertainty on the locations of detected targets and of the potential targets not yet detected. This criterion relies on the evaluation of a measure of the sets guaranteed to contain the targets, which are recursively updated while the UAVs are scanning the search zone. Evaluation of the variations of these sets accounts for potential obscuring of the UAVs FoV by obstacles located on the zone. Simulations considering several UAVs show that the proposed control input optimization is able to provide good localization and tracking performance for multiple targets. Further extensions of the proposed approach include taking into account false-alarm and non-detection in the estimation scheme and evaluating the search behavior when the connectivity graph evolves when the UAVs move.

Figure

  Figure 1 illustrates the effect of an obstacle on F z3=0 i (x U i,k ).

Figure 2

 2 Figure 2 presents the variations of the criterion on a time horizon T h = 2T for two cooperating UAVs. The control inputs are limited to variation of heading angle. As illustrated, the uncertainty criterion presents several

Fig. 3 .Fig. 4 .

 34 Fig. 3. Initial Z 1 (top left), Z 15 (top right), Z 100 (bottom left), and Z 100 (bottom right)

Fig. 5 .Fig. 7 .

 57 Fig. 5. Z 10 (top left), Z 15 (top right), Z 15 (bottom left), and Z 40 (bottom right)