Lucas Joseph

Vincent Padois

Guillaume Morel

Guillaume Morel Online

Online minimization of the projected mass of a robot

Online minimization of the projected mass of a robot for safe workspace sharing with a human Lucas Joseph 1,3 Abstract-In contexts where a human user and a robot share their workspace, the two major sources of danger for the human are related to impacts where kinetic energy is transferred from the robot to the human body and to quasi-static contacts where contact forces are applied to the human body. In order to reduce the robot kinetic energy while maximizing the robot velocity, a solution consists in minimizing the apparent mass of the robot in direction of the human. Using an online measurement of the relative position of a human and a redundant collaborative robot, the work presented in this paper features a control solution which minimizes the robot apparent mass in direction of the human being while ensuring the task assigned to the robot.

Index Terms-kinetic energy, apparent mass, redundant robot, pHRC

I. INTRODUCTION

Human robot collaboration addresses many challenges in terms of safety. When robots perform autonomous motions in an environment shared with humans, one must consider the possible occurrence of undesired contacts. To improve safety, several contributions explore the possibility to avoid contact [START_REF] Fabrizio Flacco | A Depth Space Approach for Evaluating Distance to Objects[END_REF] [START_REF] Khatib | Real time obstacle avoidance for manipulators and mobile robots[END_REF]. However, in constantly evolving environments, contact cannot always be foreseen and avoided. Different control strategies have been proposed to reduce the dangerousness of the contact [3] [4]. If contact has been established, other contributions propose to detect and react to this collision [START_REF] Haddadin | Robot collisions: A survey on detection, isolation, and identification[END_REF].

ISO TS 15066 [START_REF]Robots and robotic devices -Collaborative robots. International Organization for Standardization[END_REF] recommends that during a brief contact with a human, one should prevent the robot from delivering too much kinetic energy. Indeed, this is the energy that will be transferred to the human being. Hence, one must limit the robot kinetic energy at all time to prevent any dangerous contact. From a control point of view this is often dealt with more or less explicitly by monitoring the robot velocity. However, in many applications, the velocity is fixed by the task to realise. Therefore, these solutions reduce the tasks performances for the sake of safety.

In this work, similarly to [START_REF] Mansfeld | Improving the Performance of Biomechanically Safe Velocity Control for Redundant Robots through Reflected Mass Minimization[END_REF], it is proposed instead to minimize the robot projected mass in the direction of an obstacle. Indeed, this mass plays an important part in the transferred energy in case of contact. This mass is a function of the robot configuration and can be minimized, from a control point of view, if the robot is redundant.

II. ROBOT PROJECTED MASS

At each instant, the robot kinetic energy projected in a direction u is

e c,u = 1 2 m u v 2 u ( 1
)
with v u and m u respectively the robot velocity and mass projected in the direction u. This projected mass is expressed as

m u (q) = u T J(q)M -1 (q)J T (q) -1 u. ( 2
)
where q ∈ R n is the joint space configuration, M (q) ∈ R n×n is the joint-space inertia matrix and J(q) ∈ R m×n is the robot Jacobian matrix for a given task requiring m degrees of freedom. The robot projected mass depends on the robot configuration. For a redundant robot (m < n), it is thus possible to find a configuration minimizing this projected mass and thus minimizing the robot kinetic energy in a given direction without influencing the robot main task.

In [START_REF] Padois | Controlling dynamic contact transition for nonholonomic mobile manipulators[END_REF], it is proposed to perform a gradient descent on the projected mass, to determine a new robot configuration which does not influence the main task but minimizes the robot projected mass. This gradient descent can be written q(k + 1) = q(k) -α I n -J(q) + J(q)

T ∇m u (q)∆t. [START_REF] Joseph | Towards X-ray medical imaging with robots in the open: safety without compromising performances[END_REF] with q(k), the new robot configuration at the k th iteration of the gradient descent algorithm, J + the Moore-Penrose pseudo-inverse of J(q), α < 1 a step size and ∇m u the gradient of m u . The gradient descent algorithm returns a new desired configuration, q des , minimizing the robot projected mass without influencing the main task. Given this desired position, it is possible to compute a desired torque to send to reach this configuration τ mu = k p (q desq) -k d q.

(4)

This desired torque is inserted as a regularisation task inside a Quadratic Programming problem as described in [START_REF] Joseph | Towards X-ray medical imaging with robots in the open: safety without compromising performances[END_REF].

In [START_REF] Mansfeld | Improving the Performance of Biomechanically Safe Velocity Control for Redundant Robots through Reflected Mass Minimization[END_REF], this gradient descent is performed to minimize the robot projected mass along the direction of motion of the robot. In this work, it is rather proposed to perform this minimization in the direction of the closest human.

III. SETUP

The algorithm presented in section II is tested on a 7 degrees of freedom LWR 4+ KUKA robot. RGBD cameras are used to capture the position of a human. The RGBD cameras are calibrated and their position towards the robot are known. Three cameras are used to avoid any occlusion during the interaction with the robot. The obtained point clouds are merged to obtain a single point cloud representing the human position relatively to the robot. This position is used to compute the direction of the human towards the robot u. The whole setup is the one described in [START_REF] Meguenani | Energybased control for safe human-robot physical interactions[END_REF] and depicted in Figure 1.

IV. RESULTS

In the following experiment, the robot is required to keep a Cartesian pose while a human moves around it. The experiment shows the robot projected mass with and without the projected mass minimization algorithm. For both experiments, the human-robot distance, d, is recorded. The angle, θ, between the direction u and the unit vector [0, 1, 0] T is computed. During the two experiments, the human performs similar motions around the robot. ), it can be seen that for the same motions of the human around the robot, its projected mass in the direction of the human ranges between 2.5 kg and 2.8 kg which represents an improvement of up to 38%.

The regularisation task can thus be used to minimize the robot perceived mass. Using external sensors it is possible to minimize the perceived mass in the direction of the closest human.

V. CONCLUSION

Working in collaboration with robots requires considering the possible occurrence of a contact. To prevent any dangerous contact, traditional control solutions propose to reduce the robot velocity when a human is too close. In this work it is shown that it is also possible to adapt the robot configuration to minimize its projected mass and thus its kinetic energy in the direction of contact. If a contact occurs, it is the direction where the most kinetic energy will be transferred between the robot and the human. It allows to use the robot at its maximal velocity before reaching its maximum allowed kinetic energy as recommended by ISO TS15066.

1

 1 ISIR Sorbonne Université, CNRS UMR 7222 F-75005 Paris, France lucas.joseph@sorbonne-universite.fr Vincent Padois 1,2 2 Auctus Inria / Bordeaux INP, IMS F-33405 Talence, France vincent.padois@inria.fr Guillaume Morel 1,3 3 INSERM U1150 Agathe-ISIR F-75005 Paris, France guillaume.morel@sorbonne-universite.fr

Fig. 1 .

 1 Fig. 1. Scene observed by the RGBD cameras (extracted from [8])

Figure 2

 2 depicts the results of this experiment. During the first part of the experiment (Figure a.), the robot projected mass ranges between 2.6 kg and 4.3 kg. In the second part of the experiment (Figure b.

Fig. 2 .

 2 Fig.2. Reconfiguration of the robot to minimize its perceived mass in the direction of a human. a. shows the projected mass for a given configuration of the robot when the human moves around it. b. shows the projected mass for the same human motion but with the active reconfiguration of the robot to minimize the projected mass.


domain: spi phys

text:
Gregory R Brady
email: gbrady@stsci.edu

Christopher Moriarty

Peter Petrone

Iva Laginja

Keira Brooks

Tom Comeau

Lucie Leboulleux

Rémi Soummer

Phase-Retrieval-Based Wavefront Metrology for High Contrast Coronagraphy

Keywords: Phase retrieval, coronagraph, wavefront-sensing, deformable mirror, high-contrast imaging

We discuss the use of parametric phase-diverse phase retrieval as an in-situ high-fidelity wavefront measurement method to characterize and optimize the transmitted wavefront of a high-contrast coronagraphic instrument. We apply our method to correct the transmitted wavefront of the HiCAT (High contrast imager for Complex Aperture Telescopes) coronagraphic testbed. This correction requires a series of calibration steps, which we describe. The correction improves the system wavefront from 16 nm RMS to 3.0 nm RMS.

INTRODUCTION

The next generation of space telescopes for direct imaging and spectroscopy of exoplanets will require coronagraphs with high-contrasts, on the order of 10 -10 . To achieve this, iterative speckle nulling procedures will be used with science data to eliminate residual light passing through the coronagraph by perturbing a deformable mirror (DM) and correcting the end-to-end system wavefront on the scale of a few picometers. However, most of these procedures assume a linear approximation when relating the measured intensity to the wavefront phase, which becomes invalid for phase distributions with an amplitude larger than a small fraction of the wavelength. [START_REF] Give'on | Broadband wavefront correction algorithm for high-contrast imaging systems[END_REF][START_REF] Bottom | Speckle nulling wavefront control for Palomar and Keck[END_REF][START_REF] Miller | Spatial linear dark field control: stabilizing deep contrast for exoplanet imaging using bright speckles[END_REF] Thus we desire a method for decreasing the coronagraph system's wavefront error so that it is within the valid range of the linear approximation.

In this work we describe a procedure for measuring the wavefront of a well-aligned system using a phase retrieval algorithm based on an explicit optimization of parameters describing the phase and other pertinent system characteristics. Phase retrieval algorithms accurately model the propagation from the pupil plane to image plane and do not make use of a linear approximation. The wavefront measured using phase retrieval is used to derive a correction that is applied to the DM. Phase retrieval is not sufficiently sensitive to perform measurements on the picometer scale and is used to pre-correct the system wavefront error to the order of a few nanometers prior to speckle nulling.

In addition, we will describe a series of phase-retrieval-based calibration steps that are necessary to ensure that the wavefront correction is properly registered to the DM. Finally, we describe our final measurement of the system wavefront demonstrating dramatic improvement, bringing the wavefront error well within the valid range of the linear phase approximation.

PARAMETRIC PHASE RETRIEVAL ALGORITHM

Phase retrieval algorithms typically infer the phase in the exit pupil of an optical system using one or more measurements of the intensity patterns that the system produces in a plane near its nominal focus. [START_REF] Fienup | Reconstruction of an object from the modulus of its Fourier transform[END_REF] A parametric phase retrieval algorithm specifies a system to be considered using a set of parameters and uses an optimization algorithm to find values of these parameters that gives a minimum value for an objective function (also termed an error metric). The parameters used typically specify the phase in the exit pupil using a polynomial expansion, such as Zernike polynomials, or another basis set, such as family of delta-functions to specify phase at individual pixel locations. [START_REF] Fienup | Phase retrieval algorithms: a comparison[END_REF][START_REF] Fienup | Phase retrieval for undersampled broadband images[END_REF] Other parameters can make the algorithm more robust by allowing it to correct for inexactly known system parameters, such as the distance from the exit pupil to the measurement plane or planes, [START_REF] Brady | Application of Phase Retrieval to the Measurement of Optical Surfaces and Wavefronts[END_REF][START_REF] Thurman | Amplitude metrics for field retrieval with hard-edged and uniformly illuminated apertures[END_REF] the centration of the measured spots on the detector array, [START_REF] Thurman | Amplitude metrics for field retrieval with hard-edged and uniformly illuminated apertures[END_REF][START_REF] Brady | Optical wavefront measurement using phase retrieval with transverse translation diversity[END_REF] the amplitude distribution in the pupil plane, [START_REF] Brady | Nonlinear optimization algorithm for retrieving the full complex pupil function[END_REF][START_REF] Thurman | Complex pupil retrieval with undersampled data[END_REF] or the sampling rate of the image data. [START_REF] Jurling | Phase retrieval with unknown sampling factors via the two-dimensional chirp z-transform[END_REF] Further, these algorithms have the ability to jointly optimize over multiple images, allowing for a more robust estimate of the wavefront. This ability to compensate for systematic uncertainties make a parametric, optimization-based algorithm an attractive choice when compared to iterativetransform-type algorithms such as the Gerchberg-Saxton-Misell algorithm. [START_REF] Fienup | Phase retrieval algorithms: a comparison[END_REF][START_REF] Gerchberg | Practical agorithm for determination of phase from image and diffraction plane pictures[END_REF][START_REF] Misell | A method for the solution of the phase problem in electron microscopy[END_REF][START_REF] Dean | Phase retrieval algorithm for JWST flight and testbed telescope[END_REF] Iterative transform algorithms are less flexible and less robust (i.e. more prone to stagnating in local minima) and rely on properties of the algorithm alone for their convergence characteristics.

Objective Function

Here we employ a bias and gain-insensitive objective function to characterize the agreement between our modeled and measured intensity patterns, 8, 16

Φ = 1 - 1 K K k p,q W k (p, q)G k (p, q) G k (p, q) 2 p,q W k (p, q)G 2 k (p, q) p,q W k (p, q) G 2 k (p, q) , (1)
where p and q are sample indices, G k (p, q) is one of K measured intensity patterns, G k (p, q) is the corresponding calculated intensity pattern and W k (p, q) is a weighting function used to ignore the contribution of bad pixels or to favorably weight certain portions of the intensity pattern with, for example, better Signal-to-Noise ratio (SNR). Throughout this paper we use the convention that variables topped with a circumflex, such as G k (p, q), are being estimated while other values are known through measurement or convention. If G k (p, q) and G k (p, q) are identical, the value of the objective function is zero.

Phase Model

We typically model the field, E p (m, n), in the exit pupil using one of two methods. The first method is by simply using sampled representations of the amplitude A (m, n) and phase φ (m, n) so that we have

E p (m, n) = A (m, n) exp i φ (m, n) (2)
at sample indices (m, n). We term this representation as point-by-point.

The second method is by representing the phase by the coefficients α j of a basis set such as the Zernike polynomials Z j (m, n) so that we have

E p (m, n) = A (m, n) exp   i j α j Z j (m, n)   . (3)
It is common practise to initially optimize over a few low order Zernike terms before continuing the optimization with higher order Zernikes or point-by-point phase values.

Propagation Model

We employ a two step propagation to calculate the estimated field in each of the measurement planes from the field in the pupil plane. First, Fresnel diffraction is used to calculate the field E f (p, q) in the paraxial focal plane,

E f (p, q) = 1 N exp iπ λf p 2 ∆ 2 p + q 2 ∆ 2 q N m,n E p (m, n) exp -i2π N (mp + nq) , (4)
where λ is the wavelength, f is the distance between the exit pupil and paraxial focal plane, N is the size of the square arrays representing the pupil and focal plane fields, and (∆ p , ∆ q ) are the spacing between samples (in units of length) in the paraxial focal plane. This is typically equivalent to the pixel spacing of the detector used to record the image data. The second step uses angular spectrum propagation to propagate the small distance between the paraxial focal plane and each of the slightly defocused measurement planes. We calculate the angular spectrum of plane waves U f (m, n) in the paraxial focal plane using a Fourier transform

U f (m, n) = 1 N N p,q E f (p, q) exp -i 2π N (pm + nq) . (5)
The angular spectrum is then propagated to each of the K measurement planes,

U k (m, n) = U f (m, n) exp i2π z k 1 λ 2 -m 2 ∆ 2 m -n 2 ∆ 2 n , (6)
where (∆ m , ∆ n ) are frequency domain sample spacings (in units of inverse length) given by

∆ m,n = 1 N ∆ p,q (7)
and z k is the distance to the k th measurement plane. Note that we have included a circumflex in z k as we typically include this as a free parameter in the optimization to compensate for possible positioning errors. The propagated field in the spatial domain is computed from the angular spectrum using an inverse Fourier transform,

E k (p, q) = 1 N N m,n U k (m, n) exp i 2π N (mp + nq) . (8)
The intensity in the measurement plane is then determined from the field,

I k (p, q) = E k (p, q) 2 . (9)
A final step is often performed on this intensity image to shift the intensity spot laterally to model the effect of a detector array not centered at the origin. This is done on a sub-pixel scale using the Fourier shift theorem. First, we perform the forward transform

f k (m, n) = 1 N N p,q I k (p, q) exp -i 2π N (mp + nq) , (10)
then we multiply by the linear phase Fourier filter

g k (m, n) = exp -i2π N (m x k + n y k ) f k (m, n) , (11)
where ( x k , y k ) is the amount by which the image is shifted in pixel units. Note that these are also typically optimized for. Inverse transforming this gives us our final modeled intensity distribution,

G k (p, q) = 1 N N m,n g k (m, n) exp i 2π N (mp + nq) , (12)
which, along with the measured intensity data G k (p, q), is used to calculate the error metric given in Equation 1.

Optimization

We typically employ gradient-based algorithms to minimize the objective function, such as the conjugate gradient [START_REF] Press | Numerical Recipes 3rd Edition: The Art of Scientific Computing[END_REF] or limited-memory BFGS [START_REF] Broyden | Quasi-newton methods and their application to function minimisation[END_REF][START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF][START_REF] Zhu | Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization[END_REF] algorithms. These algorithms are most useful when efficient analytic expressions for the derivatives of the objective function with respect to the system parameters are available, rather than computing the derivatives using finite differences. This is particularly true in the case of point-bypoint estimates of the phase, where N × N derivatives must be calculated and N is typically a few hundred to a few thousand. Efficient expressions for all of the parameters discussed have been derived and are compiled in the appendix of Reference 8. To reiterate, our phase retrieval code has the capability to optimize over

• Lateral spot position in each phase-diverse image, x k , y k • Axial location of each phase-diverse image, z k • Pupil amplitude distribution, A (m, n)
• Pupil phase expressed as point-by-point phase values, φ (m, n)

• Pupil phase expressed as coefficients of Zernike polynomials (or other basis), α j .

These can be used jointly or individually, with the output of a preceding optimization serving as the input for subsequent optimizations, as is necessary for the problem at hand. In the code used in this work, the following sequence is typically used:

1. Preprocess image data to calculate data weight masks, W k (p, q), and initial guesses for lateral positions from image centroids.

2. Optimize over lateral positions alone, using guess from centroids.

3. Optimize jointly over lateral position (using result from step 2), axial position and Zernike coefficients.

4. Fix values for lateral and axial positions at values determined in step 3, optimize over point-by-point phase values using Zernike-derived phase as initial guess.

5. Optimize over lateral and axial position and point-by-point phase values using results from step 4.

6. Optimize over lateral and axial position, point-by-point phase using results from step 5 as initial values, and pupil amplitude using assumed pupil as inital guess.

This multiple-step process allows us to determine the unknown parameters of the system robustly. If we were to proceed directly to step 6 it is likely that the optimization would stagnate in a local minima. Steps 4 to 6 are necessary only when point-by-point phase values are desired and may be omitted when Zernikes are all that are needed.

Data Weights

The data weights W k (p, q) are typically derived from the measured data in an attempt to weight the data suitably for the optimization to be done. Typically they are formed by thresholding the input images so that only the central core of the image is retained. However, there are regions near the image center where the intensity is below the threshold value. Furthermore this data is important. To retain these regions a dilation operation is performed to fill in the mask. The resulting weight mask is highly suitable for retrieving the phase over low spatial frequencies, particularly for determining Zernike coefficients. This type of mask is often useful to reject noise from low intensity regions of the images, which is often where the SNR is unacceptably low. This type of mask is used in steps 2 and 3 above.

In our application we need to work with data that has significant spatial frequency content, where a data mask retaining only the central core would not recover the phase information of interest. Further, this high spatial frequency information typically has much less energy than the central core. If a uniform weighting is used across the frame the central core will still dominate and the high-frequency information will effectively be ignored. To recover this information we have found that it is necessary to weight the region away from the central core much more heavily than the central core. For example, we have achieved good results with a weight of 1 in the central core and a weight of 500 outside it. We employed this type of mask in steps 4 to 6 above. This weighting requires data that has good SNR in the high-frequency region, which we typically obtain by registering and summing hundreds of individual frames.

EXPERIMENT

Our goal in this work is to measure, and subsequently correct, the phase of the HiCAT optical system up to the occulting aperture (or Focal Plane Mask, FPM) of the coronagraph. A sketch of the system is shown in Figure 1.

The system uses a single-mode optical fiber collimated by an Off-Axis Parabolic (OAP) mirror, labeled as OAP1 in Figure 1, to illuminate the pupil of a telescope simulator. The telescope simulator consists of a pupil amplitude aperture that is imaged onto an Iris AO hexagonal-segment deformable mirror (IRIS-AO DM) by OAP2 and a parabolic mirror. This DM allows us to simulate the segments of a segmented aperture telescope. The DM pupil plane is again reimaged by the parabola and OAP3 onto a reflective apodizer needed for an Apodized Pupil Lyot Coronagraph (APLC). The apodizer is the first step in controlling the light in the dark zone of the coronagraph. The apodizer is also used as a tip/tilt mirror to stabilize the beam under servo control. The apodizer pupil plane is imaged by OAP4 and a toroidal mirror, TOR1, onto a Boston Micromachines continuous face-sheet deformable mirror, DM1. DM1 is used to correct for residual static aberrations of the system, the aberrations that we are measuring in this work. The beam reflected from DM1 is directed onto a second identical deformable mirror, DM2. DM2 is a short distance away and is not exactly conjugate to the other pupils of the system. Together DM1 and DM2 are used to suppress light in the dark zone of the coronagraph. After DM2, a second toroidal mirror, TOR2, is used to focus the light to a focus at which the FPM is placed when the coronagraph is operating. The FPM is responsible for blocking the majority of the star's light. However, between TOR2 and the FPM, a high-quality mirror can be inserted to intercept the converging beam and direct it to a camera on a motorized stage. This camera is used to collect phase-diverse phase retrieval data sets, the results of which are our topic here. A very high-quality point-source image is formed at the FPM, allowing it to block star light most effectively. Beyond the FPM, a spherical mirror is used to re-image the pupil to the Lyot stop of the coronagraph. The Lyot stop blocks light diffracted by the hard edge of the FPM that falls outside of the system pupil. Finally, fold mirrors and lenses are used to form a coronagraphic image with our desired dark zone or a pupil image on the respective cameras.

We have made a small number of notable changes from the nominal arrangement described above to help ensure good wavefront correction and facilitate phase retrieval measurements using our current software. First, the nominal pupil aperture has been replaced with a 22 mm diameter aperture to cover the entire area of the Boston Micromachine DMs. This allows us to better correct aberrations that are near the edge of the apodizer pupil (circumscribed diameter of 19.725 mm, slightly undersized compared to the DM diameter for alignment tolerancing purposes). Because of the apodized profile, we use an 18mm aperture mask as a full-aperture proxy to evaluate the wavefront over the apodizer. Secondly, the IRIS-AO DM and apodizer have been replaced with surrogate flat mirrors. This allows us to use our current phase retrieval code which does not currently handle the case of a segmented, apodized pupil. Code is in development to handle this situation.

Pupil Calibration

In order to apply the measured phase distributions it is necessary to perform a geometric calibration so that the wavefront can be applied to the DM. This was achieved by applying calibration phase patterns to the DM and recovering them using phase retrieval. The first step of the calibration is to apply a letter "F" pattern to the DM, which allows us to determine how the DM and phase retrieval results are flipped, left-right and up-down, and rotated with respect to each other. In this step we also confirmed the polarity of the retrieved phase as compared to the voltages applied to the DM. The DM command and retrieved phase are shown in Figure 2, with the determined flips and rotations applied.

The second calibration step is to determine the location of the center of the DM in the array of data calculated by the phase retrieval algorithm. To facilitate this, four actuators around the center of the DM are poked, and the associated phase is retrieved. This is illustrated in Figure 3. The intersection of two diagonal lines connecting the actuators unambiguously locates the center of the DM.

Then, the mapping of the recovered wavefront to the individual DM actuators is determined. This is done by putting a "checkerboard" pattern on the DM in which every fourth actuator is poked, illustrated in Figure 4. Sixteen data sets were recorded where each actuator in a 4 × 4 cell is poked. From this, the location of each actuator can be unambiguously determined in the array of wavefront data by comparison with the corresponding phase retrieval result. The phase retrieval results corresponding to Figure 4 are shown in Figure 5. Mapping the pupil in this way calibrates for any pupil imaging distortion between the DM plane and the exit pupil of the

Poked Column Varies

Uncorrected System Wavefront Measurement

In order to measure the system wavefront, we first applied our best-known flat maps to DM1 and DM2, determined before their installation using a Fizeau interferometer. Thus the system was as well-corrected as possible without any additional correction using phase retrieval. The recorded intensity data is shown in the top row of Figure 6. Each of these images is the sum of 200 registered images. This summing process results in an image with high dynamic range and excellent SNR. The phase recovered from this data using the first 45 Zernike polynomials is shown in Figure 7 and the corresponding modeled intensity distributions are shown in the bottom row of Figure 6. The Peak-to-Valley (PV) wavefront error is 120 nm and the Root-Mean-Squared (RMS) wavefront error is 16 nm. retrieved measured Figure 6. Intensity patterns measured in the system before any wavefront correction was performed using phase retrieval. The images in the top row are the minimally processed intensity measurements from the phase retrieval camera. The images in the bottom row are the intensity distributions modeled by our phase retrieval algorithm using the recovered wavefront. The distances labeled across the top of the figure indicate the axial position of measurement plane with respect to the nominal focus plane. Figure 7. Initial wavefront in system exit pupil before applying corrections. The PV wavefront error is 120 nm and the RMS wavefront error is 16 nm. The color scale is in units of nm.

System Wavefront and Point-Spread Function with Corrections Applied

The conjugate of the wavefront in Figure 7 was applied to the DM using the calibration and mapping determined above. The recorded PSF data is shown in the top row of Figure 8. The phase recovered from this data using the first 45 Zernike polynomials is shown in Figure 9. The PV wavefront error is 35.2 nm and the RMS wavefront error is 5.5 nm. This is a significant improvement over the uncorrected case. This correction is a strong indication of the precision and accuracy of our phase retrieval results and demonstrates the validity of our methodology for applying the correction to the DM. Corrected system wavefront, using correction derived from data in Figure 7, in exit pupil of system in the case where an oversized, 22 mm diameter aperture was used. The PV wavefront error is 35.2 nm and the RMS wavefront error is 5.5 nm. The color scale is in units of nm.

Final Wavefront Over 18 mm Aperture

After the first wavefront correction, the oversized 22 mm aperture was removed and the nominal 18 mm pupil aperture was inserted. A subsequent data set was taken. Example intensity images are shown in the top row of Figure 10. The retrieved wavefront is shown in Figure 11. The PV wavefront error is 23.9 nm and the RMS wavefront error is 3.0 nm. An image of the in-focus PSF of the system, resulting from the average of 1000 images is shown in Figure 12. It is difficult to distinguish this from a plot of a perfect theoretical result. 12. In-focus PSF over 18 mm aperture. The data here is the sum of 1000 registered camera frames, resulting in a very low SNR image with high dynamic range. The data is plotted on a logarithmic scale with a compressed dynamic range so that very dim outer rings are visible. The PSF is almost indistinguishable from a perfect theoretical result for a circular pupil with no aberrations and, qualitatively, is consistent with our measured wavefront error of 3 nm RMS.

SUMMARY

We have described the parametric approach to phase retrieval, including details of our forward model and description of the parameters and corresponding derivatives used as inputs to an optimization algorithm. We then described its application to the measurement of the transmitted wavefront error of a coronagraphic testbed system and its subsequent correction using a deformable mirror. We included details of the calibrations necessary to register the measured wavefront maps to the DM actuators, including confirming the orientation and polarity of the data, determining the center point of the DM and calibrating for distortions. The measured wavefront error improved from 16 nm RMS before correction to 3.0 nm RMS after correction, which is summarized in Table 1. The fact that we are able to measure a low amplitude wavefront, derive and apply a correction from this measurement, and realize a significantly improved transmitted wavefront, demonstrates the precision and accuracy of our phase retrieval algorithm and correction methodology.

Figure 1 .

 1 Figure1. HiCAT testbed layout. The testbed includes a telescope simulator (IRIS-AO 37 segment deformable mirror combined with a pupil mask to produce the central obstruction and support structures); an Apodized Pupil Lyot Coronagraph (APLC) combining an apodizer, a hard-edge focal plane mask and a Lyot Stop; Imaging camera including a coronagraphic focal plane, a pupil viewing camera; and a phase retrieval arm to measure the wavefront at the focal plane mask using a high-quality removable fold mirror to minimize non-common path errors.

Figure 2 .

 2 Figure 2. Letter F pattern (a) DM command and (b) retrieved phase. Color scale is in nm.

Figure 3 .

 3 Figure 3. (a) DM command and (b) retrieved phase for the case where four known actuators are poked to determine the location of the center of the DM aperture with respect to the phase retrieval phase map. Color scale is in nm.

  system.

Figure 4 .

 4 Figure 4. Sixteen checkerboard DM commands which drive, in turn, all sixteen actuators in a four by four cell. Note the shift of which actuator row and column is active between each of the 16 commands.

Figure 5 .

 5 Figure 5. Retrieved phase maps for the 16 checkerboard patterns applied to the DM. Color scale is in nm.

Figure 8 .

 8 Figure 8. Intensity patterns measured in the system after a correction based on the wavefront shown in Figure 7 was applied. The images in the top row are the minimally processed intensity measurements from the phase retrieval camera. The images in the bottom row are the intensity distributions modeled by our phase retrieval algorithm using the recovered wavefront. The distances labeled across the top of the figure indicate the axial position of measurement plane with respect to the nominal focus plane.

Figure 9 .

 9 Figure 9. Corrected system wavefront, using correction derived from data in Figure7, in exit pupil of system in the case where an oversized, 22 mm diameter aperture was used. The PV wavefront error is 35.2 nm and the RMS wavefront error is 5.5 nm. The color scale is in units of nm.

Figure 10 . 5 10 15 Figure 11 .

 1051511 Figure 10. Corrected PSFs, measured and retrieved over 18 mm aperture

Figure

  Figure 12. In-focus PSF over 18 mm aperture. The data here is the sum of 1000 registered camera frames, resulting in a very low SNR image with high dynamic range. The data is plotted on a logarithmic scale with a compressed dynamic range so that very dim outer rings are visible. The PSF is almost indistinguishable from a perfect theoretical result for a circular pupil with no aberrations and, qualitatively, is consistent with our measured wavefront error of 3 nm RMS.

Table 1 .

 1 Summary of wavefront errors before and after correction

	Case	RMS Wavefront Error (nm)
	Initial Aligned System (DM1 and DM2 flattened)	16
	Correction Applied over 22 mm Aperture	5.5
	Correction Applied over 18 mm Aperture	3.0

ACKNOWLEDGMENTS

This work is supported by the National Aeronautics and Space Administration under Grants NNX12AG05G, NNX14AD33G issued through the Astrophysics Research and Analysis (APRA) program (PI: R. Soummer) the JWST Telescope Scientist Investigation, NASA Grant NNX07AR82G (PI: C. Matt Mountain) and the STScI Directors Discretionary Research Funds.


Amanda J Ross

Patrick Crozet

Allan G Adam

Dennis W Tokaryk

Observations of the low-lying Ω=1/2 states of nickel deuteride, NiD

Keywords: NiH, NiD, Fourier transform laser-induced fluorescence, transition metal hydride Highlights

Vibrational levels of the low-lying 2 Σ + and 2 Π 1/2 states have been located in 58 NiD • Analysis validate predictions from a multi-isotope supermultiplet fit • A new electronic Ω' = 0.5 state has been identified, with T' v ~ 16724.5 cm -1 above the first rotational level of the ground state.

.

Introduction

This investigation of dispersed fluorescence in NiD was motivated by attempts to use isotopic substitution to improve a 'global' description of the three low-lying electronic states ( 2 Δ, 2 Π, 2 Σ + ) of Ni 1,2 H, either through 'supermultiplet' models [1], or potential energy curve fits [START_REF] Havalyova | 25th International Conference on High Resolution Molecular Spectroscopy[END_REF].

The three lowest-lying electronic states of NiH show strong effects of spin-orbit interactions, resulting in irregular vibrational spacings in all but the 2 Δ 5/2 ground state. The upper states accessible to visible laser excitation are associated with f-d excitation on Ni atoms, the strongest transitions being to upper states with significant 2 Δ 5/2 or 2 Φ 7/2 character.

Resolved laser-induced fluorescence from Ω' = 3/2, 5/2 and 7/2 excited states [START_REF] Abbasi | [END_REF] has recently located vibrational levels of 2 Δ 5/2 , 2 Δ 3/2 and 2 Π 3/2 in NiD, (corresponding quite closely to the three low-lying electronic states of NiCN reported by Kingston, Merer and Varberg [4]) but provided no information about the Ω" = ½ states. The low-lying 2 Π 1/2 and 2 Σ + states of NiH were observed via collisional energy transfer between close-lying excited states, in particular to an Ω'=1/2 upper state [5], but there was no evidence of an equivalent process in the NiD spectra discussed in ref. [START_REF] Abbasi | [END_REF].

Survey spectra.

The analysis of ref. [START_REF] Abbasi | [END_REF] was severely hampered by lack of information on the two Ω=1/2 components of the Ni + (3d 9 )-D -supermultiplet, so the objective here was to determine as many rovibrational term values as possible for the low-lying 2 Σ + and 2 Π 1/2 states. The hybrid Hund's case (a) -Hund's case (b) Hamiltonian proposed by Abbasi in [START_REF] Abbasi | [END_REF] predicted their first rotational levels (v=0) to be respectively ~2073 and 3507 cm -1 above the lowest level of the electronic ground state, but the predictive powers of this model were so far untested. We reasoned that whilst ΔΩ = ΔΛ (and particularly ΔΛ=0) transitions are the strongest in the visible systems of NiH and NiD, some ΔΩ = ± 1 systems are observed. A (weak) 2 Π 3/2 ← 2 Δ 5/2 pump transition could conceivably yield (even weaker) 2 Π 3/2 → 2 Σ + fluorescence. We surveyed the excitation spectrum of NiD in an experiment similar to those performed by Guo et al. on NiH [6]. NiD was produced in the collision-free environment of a molecular-jet laser-ablation source at University of New Brunswick, using CD 4 as a precursor. A pulsed dye laser (YAG-pumped Hyperdye HD-500, repetition rate 10 Hz, resolution approximately 0.1 cm -1 ) was scanned from 16000 to 18580 cm -1 , detecting through a 0.25 m monochromator (zero order) and using λ > 610 or 640 nm filters to reduce laser scatter, since the signals of primary interest -to excited electronic states -were expected to occur more than 1500 cm -1 below the laser line. Dispersed fluorescence from each resonance was then recorded at low resolution with a 0.5 m monochromator and red-sensitive Hamamatsu R-928 photomultiplier tube, to see whether there were any features likely to correspond to emission to Ω" = ½ states. The only excitations to show any promise were those exciting the first rotational levels of E [17.5]1.5 (v=1), using [T v ' in 1000 cm -1 ]Ω notation. This upper state had already been identified in ref. [START_REF] Abbasi | [END_REF] from collisionally-induced emission to 2 Π 3/2 , but no E→ 2 Σ + transitions had been found. Direct excitation of E[17.5]1.5 (and D[17.2]1.5) had also been attempted, but gave emission signals too weak to record in Fourier transform resolved fluorescence. In the collision-free ablation source environment, with all population concentrated in the lowest rotational levels, the strong feature around 15450 cm -1 (see Figure 1) was readily assigned as E → 2 Σ + (v=0) emission, since it falls, as predicted, in the gap between transitions to v=1 of the ground state and to v=0 of the 2 Π 3/2 state. The peaks around 14100 cm -1 were not immediately assignable; on energetic grounds, they could plausibly come from transitions either to v=1 of 2 Σ + , or to v=0 of the 2 Π 1/2 state.

Dispersed fluorescence at higher resolution

The low-resolution dispersed fluorescence spectrum gave convincing evidence of transitions at least to v=0 of the 2 Σ + state, but with insufficient precision to test or refine existing predictions. Given earlier failure to record laser-induced fluorescence following E!X excitation with the Lyon experiment, we attempted a two-laser experiment with the ablation source, counter-propagating output from two pulsed dye lasers (rhodamine 6G to pump E!X, and DCM to dump E" 2 Σ + ), and looking for depletion while monitoring fluorescence with the monochromator set 2800 cm -1 below the laser, i.e. at the strongest band in Figure 1. This was unsuccessful, because issues with electrical noise from the lasers made it impossible to find the signals needed to adjust timing of the laser pulses. Not having cw laser resources to perform a high-resolution depletion experiment with the collision-free source at UNB, we were forced to revert to the Lyon experiment to take this further, using the sputter source described in detail in ref. [7].

We needed an approximately three-fold increase in signal/noise ratio compared with spectra reported in [START_REF] Abbasi | [END_REF], to produce workable interferograms. To achieve this, we worked with a higher enrichment (15%) of D 2 in argon ("Mélange Crystal" from Air Liquide). We also passed the laser beam as close as possible to the annular anode, where the concentration of NiD radicals is highest. There are disadvantages to doing this: background scatter is higher, collisions transfer excited state population to other electronic states very efficiently, and the molecular fluorescence at the longer-wavelength end of the spectrum is dwarfed by argon atomic emission (see Figure 2).

Even with maximum amplification from the Si-avalanche detector and (Bomem DA3) interferometer's amplifier, the maximum signal at zero-path difference rarely exceeded 3% of the analog-digital converter's ± 10 V full-scale input. Nevertheless, by averaging typically 50 scans (recording time roughly 1 hour) at close to Doppler resolution (instrumental resolution 0.04 cm -1 in most cases), we could identify direct fluorescence to v=0 levels of the 2 Σ + quite easily -and made use of the argon lines to verify wavenumber calibration from the Fourier transform spectrometer.

Figure 2 illustrates fluorescence from v'=1 of the E state, with the ring dye laser tuned to excite just one of the parity components (J'= 1.5,f) accessed with the pulsed laser in Figure 1. We recorded 16 spectra following excitation of a selection of rotational levels (1.5 ≤ J' e,f ≤ 11.5 ) of E [17.5]1.5 (v=1). The R and Q branches were relatively strong in transitions to v=0 of 2 Σ + ; the P branches were stronger in transitions to v=0 2 Π 1/2 , with the R branches dying out quickly with increasing J (see Table 1). The transitions to v=0 of 2 Π 1/2 were convincingly identified and assigned by subtracting spectra of different upper-state parity at given J, as illustrated in Figure 3.

This approach removed all the 'common to both' collisionally induced emission, and greatly reduced the atomic argon lines, so that some very weak features 12600-12650 cm -1 were tentatively assigned (and later confirmed) as P f transitions to 2 Π 1/2 , v=1, with no hint of any branch to levels of the other parity. The strongest features in the lower part of Figure 2 are rotationally relaxed E[17.5]1.5 → 2 Π 3/2 (1-0) and (1-1) bands. The (1-2) band of this system, near 12000 cm -1 , is hardly perceptible in Fig 2, but it locates previously unobserved rotational levels of 2 Π 3/2 , some 5500 cm -1 above the electronic ground state. Other bands have been identified as E[17.5]1.5 → 2 Δ 3/2 (v"= 0,1,3) and E[17.5]1.5 → 2 Δ 5/2 (1-0). These bring no new information on the energy level patterns in NiD, but establish upper state energies with respect to the known low-lying 2 Δ Ω states. Table 1 lists assigned transitions from the E[17.5]1.5 state. These transitions were included in the multistate term-energy fit of ref. [START_REF] Abbasi | [END_REF] to generate T v,J parity energies for the low-lying Ω" = ½ levels. The supermultiplet model of Gray et al. established that the low-lying Ω = ½ levels are almost 50:50 mixtures of the 2 Σ + and 2 Π 1/2 states of NiH. This is true for NiD, too, and although the expression for energy levels of a 2 Σ + state (eq. 1) can be used to reproduce term energies of both states, it returns un-physical parameter values. As expected, though, the parity splittings at J=0.5 are opposite in sign for the two states: E(v=0,f) -E(v=0,e) = 19.77 cm -1 in 2 Σ + , and -9.07 cm -1 for

2 Π 1/2 . T v,F 1 N ( ) = T v + B v N N + 1 ( ) -D v N N + 1 ( ) ⎡ ⎣ ⎤ ⎦ 2 + ... + 1 2 γ 1 N + 1 2 γ 2 N 2 + ... T v,F 2 N ( ) = T v + B v N N + 1 ( ) -D v N N + 1 ( ) ⎡ ⎣ ⎤ ⎦ 2 + ... - 1 2 γ 1 N + 1 ( ) + 1 2 γ 2 N + 1 ( ) 2 + ... ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ (1)
The information on the hitherto unobserved levels in NiD provide a stringent test on the predictions made by Abbasi et al. [START_REF] Abbasi | [END_REF] from a multi-isotope fit of term energies to a modified version of the supermultiplet model originally proposed by Gray and Field [1]. Retaining spin-orbit interaction terms derived from a Ni + 3d 9 configuration, ref. [START_REF] Abbasi | [END_REF] imposed simple mass-dependence of rotational and vibrational terms so that a single (though extensive) set of parameters could be used to represent vibronic energies of the spin-orbit coupled 2 Σ + , 2 Π and 2 Δ states of 58,60,62 NiH and 58 NiD up to 6000 cm -1 above the electronic ground state. This non-standard model was not able to reproduce experimental precision for the observed term values, but was intended to guide future experimental work, indicating to within a few cm -1 where transitions should be expected. It served extremely well! The f parity term values determined for the v=0 level of the 2 Σ + state from our spectra are within 1 cm -1 of predictions from the multi-isotope deperturbation fit, likewise for levels of e symmetry in the v=0 level of the 2 Π 1/2 state (f parity overestimated by 6-8 cm -1 ). The predictions placed T J ,e and T J+1 f 7 cm -1 apart for J=6.5, 2 Π 1/2 v=1, initially making it hard to decide whether the weak features now labelled P f (J) in the lower part of Figure 2 were truly P f (J) lines, or Q fe (J+1). The assignments given in Table 1 were confirmed by re-recording fluorescence from the B 2 Δ 5/2 state v=1 with the 15% D 2 in argon mixture, and examining collisionally-induced fluorescence carefully.

Collisional energy transfer from the B 2 Δ 5/2 state v=1 populates another Ω'=1.5 state, D[17.6] 1.5. Some of the weakest features in the spectrum could now be assigned as R f lines in the (1-0) and (1-1) D[17.6] 1.5→ 2 Π 1/2 bands. Better still, the secure lower-state combination differences from the E[17.5]1.5 → 2 Σ + (1-0) spectra provided J and parity assignments for R e , P f , Q ef and Q fe branches from a 'new' excited Ω' = ½ state with T' v ~16700 cm -1 (hereafter labelled [16.7]0.5), terminating on v" = 0 of the low-lying 2 Σ + state. These transitions are illustrated in Figure 4. The 'unassigned' feature mentioned in ref. [START_REF] Abbasi | [END_REF] was the [16.7]0.5 -2 Σ + v"-0 R e branch of this system (filled circles in Figure 4), easily spotted because it happens to fall between emission lines from stronger systems. The widely-spaced (~30 cm -1 ) P f transitions were nearly twice as strong, but less noticeable because they are interspersed with direct B-X 1 (1-2) fluorescence (the strongest lines of the P and Q branches of B-X 1 (1-2) go off the scale of Figure 4 by a factor ~8). and Q fe transitions to 2 Σ + v"=1, but only P f lines were seen to 2 Σ + , v"=2. Conversely, the R f and P e branches were strong in emission to v"=0 of the W 2 2 Π 1/2 state. R f lines are assigned with confidence (via an upper-state perturbation producing line -extra-line pairs) to v"=1, while the P e assignments are based on some features barely emerging from noise that happen to match predictions for R e transitions remarkably well; these are less secure. Given that Q branches appear in the [16.7]0.5 → 2 Σ + bands, and not in [16.7]0.5 → 2 Π 1/2 , it seems reasonable to suggest that this new upper state has some 2 Π 1/2 character. It is too low in energy to be the NiD analogue of the [17.9]0.5 state observed [5] in NiH. Figure 4 indicates an upper state vibrational assignment v'=0; this comes from the isotope shift < 0.1 cm -1 for the first rotational lines of the R and P branches in 58 NiD and 60 NiD. The isotope splitting is expected to be ~ 0.3 cm -1 for Δv = ± 1. The missing, or very weak, branches in the [16.7]0.5 → 2 Σ + and the [16.7]0.5 → 2 Π bands are explained by destructive interference between parallel µ ! and perpendicular µ ⊥ transition moments to the spin-orbit-mixed lower states ( 2 Π 1/2 and 2 Σ + respectively). The strong branches are enhanced by constructive interference. These effects result from ΔΛ = ± 1 mixings produce opposing intensity changes for Δ J = +1 and Δ J = -1 lines, because transition amplitude phases are of opposite sign for the P and R lines of ΔΛ = ± 1 'perpendicular' transitions [START_REF] Lefebvre-Brion | The spectra and dynamics of diatomic molecules[END_REF]. The situation in NiD is illustrated schematically in Figure 5. Unexpectedly strong or weak P and R branches result from combinations of µ ⊥ and µ || transition moments, for example,

for P J ( ) : c 1 µ ! - 1 2 c 2 µ ⊥ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ 2 and for R J ( ) : c 1 µ ! + 1 2 c 2 µ ⊥ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ 2 ,
where signed eigenvector 'parentage' coefficients c 1 and c 2 multiply µ ⊥ and µ || to determine the relative intensities of P(J) and R(J) lines. Similar patterns were seen in NiH from a higher-lying Ω' = 0.5 state, and the phenomenon was discussed in some detail in ref. [5.] Figure 5 also shows that the parity splittings in the 16.7]0.5 state are even larger than in the lower 2 Σ + state.

Figure 5 (in colour online). Heavy lines indicate the stronger branches (arising from constructive interference between transition dipole moments) observed in emission from the [16.7]0.5 excited state to 2 Π 1/2 and 2 Σ + lower states. The Q branches to 2 Σ + are weaker (broken lines). Rotational stacks are drawn to scale within each state, indicating J and parity for e and f symmetry levels.

All assigned transitions from [16.7]0.5 (listed in table 2) were included in the multi-state termvalue fit for 58 NiD, currently comprising ~2250 transitions and returning an unweighted root mean square deviation of 0.0056 cm -1 . Term energies for the 2 Σ + state were reduced to effective band parameters (Table 3) using expression (1); the new data come from optical transitions shown schematically in Figure 6. A more restricted fit to eq (1) was used for [16.7]0.5 terms, eliminating f-parity data with J ≥ 9.5 (clearly affected by a local perturbation); the resulting parameters are shown in the last column of Table 3.

Figure 6 (in colour online). Schematic view of transitions defining previously unobserved levels of the low-lying supermultiplet states of 58 NiD. Upwards arrows indicate laser pump transitions, with the Ω'=0.5 state populated by collisional energy transfer. Black downwards arrows connect to previously observed levels of the supermultiplet states. Thick lines denote strong bands.

Table 4 lists previously unreported term energies for levels of the 'supermultiplet states' of NiD extracted from the term-value fit. This includes not only levels of the 2 Σ + and 2 Π 1/2 states, but also v=2 of 2 Π 3/2 . Comparisons with predictions from Abbasi's work are also given. The agreement for the Ω" = 0.5 v"=0 levels of 2 Σ + and 2 Π 1/2 is remarkable. The predictions furthest from experiment were for 2 Π 3/2 v=2, where the rotationless energy term was overestimated by 14

cm -1 .
The relative positions of the observed ligand-field state termvalues for NiH and NiD are illustrated in Figure 7. The plot emphasizes the modest difference in rotationless energies at v=0 for 2 Π 1/2 and for 2 Σ + , contrasting with a very large shift between NiH and NiD for the 2 Π 3/2 state. The proximity between v ≥ 1 of the 2 Δ 3/2 state and 2 Π 1/2 explains the abnormal parity splittings in the 2 Δ 3/2 state, noticed in [START_REF] Abbasi | [END_REF] to be greater than in the 2 Π 3/2 state.

Conclusions

Finding the first rotational levels of the lowest Ω=1/2 state (v=0 2 Σ + ) of NiD within 1 cm -1 of predictions is a triumph for Abbasi's multi-isotope Dunham parameter-based fit, with its various constraints and assumptions (off-diagonal spin-orbit, L-and S-uncoupling matrix elements derived from Ni + atomic properties, Born-Oppenheimer breakdown effects neglected, save for calculated mass-dependent zero-point energy shifts). Looking at the columns of differences between model predictions and current observations, the model overestimated the rotationless energy of v=2 in the 2 Π 3/2 state, but reproduced the rotational structure quite well. On the contrary, the prediction for first rotational terms of v=0 and v=1 of the 2 Π 1/2 state were remarkably good, within 2 cm -1 , but the parity splittings increase more quickly with J than the model predicted. At J=7.5 the energy term for e symmetry was underestimated by 5 cm -1 , with the f parity term overestimated by around 7 cm -1 . Additional information from spectra recorded after selective excitation of the lower abundance 60 NiD will, we hope, provide sufficient information to make a potential-curve fit option viable. There remain some unassigned features in the spectra not (yet) connected to known states of NiD, or even compatible with predictions, for which we cannot even offer tentative assignments. The remarks of A. Merer and co-workers' on spectra of CaF [START_REF] Kay | [END_REF], " … Rigorous assignment is frustrated only when multiple rotational transitions to a common upper state are unobservable because of systematic intensity effects", are highly pertinent in this case, too, as mixing between states causes 'allowed' branches to vanish in the electronic spectra of NiD.

Table 1 Compilation of assigned fluorescence transitions (cm -1 ) from E[17.5] Ω =1.5. Estimated uncertainties 0.005 cm -1 for unblended lines. Comparison of experimentally determined term energies for levels of 2 Σ + (v=0) , 2 Π 1/2 (v=0,1) and 2 Π 3/2 (v=2) with predictions from the multi-isotope fit of ref. [START_REF] Abbasi | [END_REF].

E[17.5]1.5 → 2 Δ 5/2 v"=0 J" R e R f P e P f Q ef Q fe 2.
J

Fig 1 .

 1 Fig 1. (in colour online) Low-resolution dispersed fluorescence in NiD. The laser (17507.5 cm -1 ) excited P(2.5) E [17.5]1.5 ! X 1 2 Δ 5/2 (1-0). Boxes indicate the transitions of particular interest, to previously unobserved low-lying electronic states.

Figure 2 .

 2 Figure 2. (in colour online) Fourier-transform record (0.04 cm -1 resolution) of fluorescence from E [17.5] Ω =1.5, v=1, J=1.5f in 58 NiD, using a sputter source to produce NiD. The vertical arrows show enlargements of the regions with transitions to v=0 of 2 Π 1/2 and 2 Σ + states, where direct fluorescence from J'=1.5f can be distinguished from a baseline of many collisionally induced, (rotationally unselective) transitions.

Figure 3 (

 3 Figure 3 (in colour online). Taking the difference between spectra differing only in upper state parity for the pump transition (E, v=1, J = 5.5e and 5.5f in this example) makes direct fluorescence easier to distinguish from collisionally-induced relaxation. View (a) shows transitions to 2 Σ + v=0, and view (b) illustrates weaker transitions to v=0 of 2 Π 1/2 , surrounded by P and R transitions to v=2 2 Δ 3/2 (denoted by Δ symbols), with much smaller parity splittings in the lower state.

Figure 4 (

 4 Figure 4 (in colour online). Collisionally-induced fluorescence in 58 NiD, following laser excitation of 58 NiD B 2 Δ 5/2 v=1, J' e,f =2.5, showing transitions from the [16.7]0.5 upper state to the low-lying 2 Σ + state, v"=0. Panel a) illustrates the distinctive, widely-spaced P f (▼) and slightly weaker R e (•) branches). Panel b) highlights the crowded Q-branches between them; Q ef with red line markers, Q fe in blue, with markers extending below the experimental spectrum for clarity.

Figure 7 (

 7 Figure 7 (in colour online) Term value plot for NiH and NiD in the region of the lowest vibrational levels of the 2 Π and 2 Σ + states, large labels denoting new observations for 2 Π 3/2 (squares), 2 Π 1/2 (diamonds) and 2 Σ + (triangles) in NiD. Open and filled symbols distinguish between e and f parity.

Table 2

 2 Compilation of assigned fluorescence transitions (cm -1 ) from [16.7] Ω' =0.5 to lower states 2 Δ 3/2 , 2 Σ + and 2 Π 1/2 . Estimated uncertainties 0.005 cm -1 for unblended lines, 0.015 cm -1 for particularly weak branches(#). Asterisks indicate extra-lines with J'=9.5, f.

	[16.7]0.5 → 2 Δ 3/2 v"=0
	J"	R e	R f	P e	P f	Q ef	Q fe
	1.5
	2.5					overlap /	15713.221
	3.5					HeNe	15698.278
	4.5						15682.353
	5.5						15665.503
	6.5						15647.794
	7.5						15629.288
	8.5						15610.056
	9.5						15590.405
	9.5						*15589.297
	10.5						15569.255
	11.5						15548.050
	[16.7]0.5 → 2 Δ 3/2 v"=1
	J"	R e	R f	P e	P f	Q ef	Q fe
	1.5			14431.445

Table 3 .

 3 Effective parameters for vibrational levels v = 0-2 of the low-lying 2 Σ + state of 58 NiD, fitting term values to eq. 1. Although these fitting parameters (quoted with 1 standard deviation in units of last digit in parenthesis) have little physical meaning, they reproduce the term energies with 0.014 cm - 1 unweighted root mean square deviation. The last column gives parameters fitting term values J ≤ 12.5 e and J ≤ 8.5 f of the hitherto unobserved [16.7] 0.5 upper state to eq. This fit returns an unweighted rms deviation of 0.02 cm -1 .

	Parameter (cm -1 )	2 Σ + , v=0	2 Σ + , v=1	2 Σ + , v=2	[16.7] 0.5
	T v	2073.333 (6)	3353.833 (8)	4606.791(23)	16724.53(1)
	B v	3.75575(36)	3.55749 (66)	3.4067 (10)	3.3800(3)
	D v	2.0486(52) × 10 -4 4.713(131) × 10 -4	1.7288(1280) × 10 -4	1.73(33)× 10 -5
	H v	1.612(195) × 10 -7 1.429 (73) × 10 -6	3.153(440) × 10 -7
	γ v1	-12.289(3)	-9.182 (3)	-10.0	28.828(4)
	γ v2	8.372(56) × 10 -3	2.534 (12) × 10 -2		-0.01468(4)
	γ v3	-2.33(26) × 10 -6	-3.543(87) × 10 -5

Table 4

 4

Acknowledgements. Financial support for this work came from NSERC Canada, and from CNRS France as an equipment grant from CSAA/INSU. AR thanks UNB for arranging a Harrison McCain visiting fellowship in 2017, and DT is grateful to support from the CNRS as chercheur invité, for 2 months spent at ILM in 2018.


Dominique Larchey-Wendling
email: dominique.larchey-wendling@loria.fr

Jean-François Monin
email: jean-francois.monin@univ-grenoble-alpes.fr

Simulating Induction-Recursion for Partial Algorithms

We describe a generic method to implement and extract partial recursive algorithms in Coq in a purely constructive way, using L. Paulson's ifthen-else normalization as a running example.

. We claim that the use of (1) IR, which is still absent from Coq, and (2) an informative predicate (of sort Set or Type) for the domain, preventing its erasing at extraction time, can be circumvented through a suitable bar inductive predicate.

⇒ ω(α, nm y, nm z) | ω(ω(a, b, c), y, z) ⇒ nm(ω(a, nm(ω(b, y, z)), nm(ω(c, y, z)))) We illustrate our technique on L. Paulson's algorithm for if-then-else normalization [START_REF] Giesl | Termination of Nested and Mutually Recursive Algorithms[END_REF][START_REF] Bove | Modelling general recursion in type theory[END_REF] displayed in Fig. 1. For concise statements, we use ω to denote the ternary constructor for if then else expressions, and α as the nullary constructor for atoms. As witnessed in the third match rule ω(ω(a, b, c), y, z), nm contains (two) nested recursive calls, making its termination depend on properties of its semantics. This circularity complicates the approach of well-founded recursion and may even render it unfeasible.

Our method allows to show these properties after the (partial) function nm is defined, as proposed in [START_REF] Krauss | Partial and Nested Recursive Function Definitions in Higher-order Logic[END_REF], but without the use of Hilbert's εoperator. We proceed purely constructively without any extension to the existing Coq system and the recursive definition of Fig. 1 can be extracted as is from the Coq term that implements nm.

We start with the inductive definition of the graph G : Ω → Ω → Prop of nm (Fig. 2) and we show its functionality. 1 Then we define the domain/termination predicate D : Ω → Prop as a bar inductive predicate with the three rules of Fig. 3. We show a dependent induction principle for D (see Fig. 5). The term d_nm_ind states that any dependent property P : ∀e, D e → Prop contains D as soon as it is closed under the simulated constructors d_nm_[012] of D. The assumption ∀ e D 1 D 2 , P e D 1 → P e D 2 restricts the principle to proof-irrelevant properties about the dependent pair (e, D e ). This is exactly what we need to establish properties of nm. Then we can show partial correctness and termination as in [START_REF] Giesl | Termination of Nested and Mutually Recursive Algorithms[END_REF] in this example, nm happens to always terminate on a normal form of its input. In a more relational approach, these properties can alternatively be proved using nm_spec and induction on G x n x . Though our approach is inspired by IR definitions, in contrast with previous work, e.g. [START_REF] Bove | Another Look at Function Domains[END_REF], the corresponding principles are established independently of any consideration on the semantics or termination of the target function (nm), i.e. without proving any properties of D/nm a priori. This postpones the study of termination after both D and nm are defined together with constructors and elimination scheme, fixpoint equations and proof-irrelevance. Moreover, our domain/termination predicate D is non-informative, i.e. it does not carry any computational content. Thus the code obtained by extraction is exactly as intended.

G α α G y n y G z n z G (ω α y z) (ω α n y n z ) G (ω b y z) n b G (ω c y z) n c G (ω a n b n c ) n a G (ω (ω a b c) y z) n a Figure 2: Rules for the graph G : Ω → Ω → Prop of nm. D α D y D z D (ω α y z) D (ω b y z) D (ω c y z) ∀ n b n c , G (ω b y z) n b → G (ω c y z) n c → D (ω a n b n c ) D (ω (ω a b c) y z)
2 i.e. nm_pirr : ∀ e D 1 D 2 , nm e D 1 = nm e D 2 .

Our Coq code is available under a Free Software license [START_REF] Larchey-Wendling | The If-Then-Else normalisation algorithm in Coq[END_REF]. We have successfully implemented other algorithms using the same technique: F91, unification, depth first search as in [START_REF] Krauss | Partial and Nested Recursive Function Definitions in Higher-order Logic[END_REF], quicksort, iterations until 0, partial list map as in [START_REF] Bove | Partiality and recursion in interactive theorem provers -an overview[END_REF], Huet&Hullot's list reversal [START_REF] Giesl | Termination of Nested and Mutually Recursive Algorithms[END_REF], etc. The method is not constrained by nested/mutual induction, partiality or dependent types. On the other hand, spotting recursive sub-calls implies the explicit knowledge of all the algorithms that make such calls, a limitation that typically applies to higher order recursive schemes such as e.g. substitutions under binders. Besides growing our bestiary of examples, we aim at formally defining a class of schemes for which our method is applicable, and more practically propose some automation like what is done in Equations [START_REF] Sozeau | Equations: A Dependent Pattern-Matching Compiler[END_REF].

  type Ω = α | ω of Ω * Ω * Ω let rec nm e = match e with | α ⇒ α | ω(α, y, z)

Figure 1 :

 1 Figure 1: L. Paulson's if-then-else normalisation algorithm.

Figure 3 :

 3 Rules for the bar inductive definition of D : Ω→Prop.There, we single out recursive calls using G but proceed by pattern-matching on e following the recursive scheme of nm of Fig.1. Then we define nm_rec : ∀ e (D e : D e), {n | G e n} as a fixpoint using D e to ensure termination. However, the term nm_rec e D e does not use D e to compute: the value n satisfying G e n is computed by pattern-matching on e and recursion, following the scheme of Fig.1.Finally, we define nm e D e := π 1 (nm_rec e D e ) and get nm_spec e D e : G e (nm e D e ) using the second projection π 2 . Extraction of OCaml code from nm outputs exactly the algorithm of Fig.1, illustrating the purely logical (Prop) nature of D e . In order to reason on D/nm we show that they satisfy the IR specification given in Fig.4: the constructors of D are sufficient to establish the simulated constructors d_nm_[012], while nm_spec allows us to derive the fixpoint equations of nm. Us-ing g_nm_fun, we get proof-irrelevance of nm. 2Inductive Ω :Set := α : Ω | ω : Ω → Ω → Ω → Ω.Inductive D : Ω → Prop := | d_nm_0 : D α | d_nm_1 y z : D y → D z → D(ω α y z) | d_nm_2 a b c y z D b D c : D ω a (nm (ω b y z) D b ) (nm (ω c y z) D c ) → D ω (ω a b c) y z with Fixpoint nm e (D e : D e) : Ω := match D e with | d_nm_0 → α | d_nm_1 y z D y D z → ω α (nm y D y ) (nm z D z ) | d_nm_2 a b c y z D b D c D a → nm ω a (nm (ω b y z) D b ) (nm (ω c y z) D c ) D a end.

Figure 4 :

 4 Figure 4: IR spec. of D : Ω → Prop and nm : ∀e, D e → Ω.

  Theorem d_nm_ind (P : ∀e, D e → Prop) : ∀ e D 1 D 2 , P e D 1 → P e D 2 → P _ d_nm_0 → ∀ y z D y D z , P _ D 1 → P _ D z → P _ (d_nm_1 y z D y D z ) → ∀ a b c y z D b D c D a , P _ D b → P _ D c → P _ D a . . . . . . → P _ (d_nm_2 a b c y z D b D c D a ) → ∀ e D e , P e D e .

Figure 5 :

 5 Figure 5: Dependent induction principle for D : Ω → Prop.

i.e. g_nm_fun : ∀ e n 1 n

, G e n 1 → G e n 2 → n 1 = n 2 .


Maxime Forster
email: maxime.forster@insa-lyon.fr

Erwann Fourmond

Roland Einhaus

Hubert Lauvray

Jed Kraiem

Mustapha Lemiti

Ga co-doping in Cz-grown silicon ingots to overcome limitations of B and P compensated silicon feedstock for PV applications

Keywords:

In this paper, we investigate gallium co-doping during crystallisation of boron and phosphorus compensated Si. It is shown that the addition of gallium allows to obtain a fully p-type ingot with high resistivity despite high B and P contents in the silicon. Segregation of doping impuri-ties is consistent with theory with a slight adjustment to take Ga evaporation into account. We have determined that 10 to 15% of P atoms are electrically inactive. Relatively high minority carrier lifetimes are measured in this material.

1 Introduction On the way to reduce the cost of solar energy, one promising path is the usage of low cost Si feedstock, or UMG-Si, obtained by simplified purification processes [START_REF] Einhaus | PHOTOSIL Simplified production of solar silicon from metallurgical silicon in[END_REF]. One drawback of UMG-Si is that it contains large amounts of both B and P, these impurities being the most difficult to remove through the metallurgical route. In such a compensated Si, electrical properties are differing from solely B-doped Si. Previous results [START_REF] Dubois | [END_REF] have shown the minority carrier lifetime to be higher and thus to give improved solar cell performances when the compensation level (Cl) was high i.e. for low net doping NA-ND (high resistivity). In addition, high Cl seems to have a positive effect on the so-called light-induced degradation observed in B-doped and oxygen contaminated Si [3,4]. However, because of the lower segregation coefficient of phosphorus (0.35) compared to that of boron (0.8), the net doping is not uniform along ingot height and an inversion of the Si polarity, from p-to n-type may occur during the crystallisation, reducing the material yield for industrial solar cell fabrication. Simulations [START_REF] Enebakk | 3rd Intenational Workshop on Crystalline Silicon Solar Cells[END_REF] using the Scheil's law for the dopants distribution and simple models for the carrier mobility calculation show that a minimum resistivity of 0,5Ω.cm required for solar cell fabrication and an ingot yield of at least 90% implies the use of Si feedstock containing less than 0.45ppmw of B and 1.0ppmw of P. Galli-um co-doping has recently been proposed [START_REF] Kraiem | Proceedings 34th IEEE Photovoltaic Specialists Conference[END_REF][START_REF] Macdonald | Cuevas in: 19th Workshop on Crystalline Silicon Solar Cells & Modules: Materials and Processes[END_REF][START_REF] Kraiem | [END_REF] as a potential solution to control the Cl along the ingot height even when using Si containing large quantities of B and P. It relies on the low segregation coefficient of Ga which allows counterbalancing the increase of P concentration over that of B during crystal growth. This thus permits to obtain a low net doping and a p-type resistivity along the whole ingot height. This method has been experimentally validated [START_REF] Kraiem | Proceedings 34th IEEE Photovoltaic Specialists Conference[END_REF] on a mc-Si ingot grown with Ga co-doped UMG-Si. A deeper understanding of the effect of the dopants distribution was however difficult to reach due to the presence of other impurities and crystal defects that are present in mc-Si. In this work, we explore the influence Ga co-doping on a Cz-Si crystal, grown with high purity solar grade Si to which was added known concentrations of B, P and Ga. The study focuses on the segregation effects during crystal growth and on the electrical properties dependence with the doping impurities concentrations, such as dopant activation and minority carrier lifetime.

2 Experimental details For this study, an <100> oriented Si ingot of about 160mm diameter was grown by the Cz method, using 15kg of Hemlock solar grade poly-Si free of doping impurities to which was added known concentrations of B, P and Ga (Table 1). The concentrations of B and P correspond to a UMG-Si from which B has been efficiently removed by purification but which still contains relatively high amounts of P. The added Ga concentration is calculated using the Scheil equation in order to obtain a fully p-type ingot and a relatively high resistivity. B and P were introduced through highly doped Si wafers in which the doping impurity concentration was evaluated by means of 4-points probe resistivity measurements correlated to the Irvin curves. Ga was introduced in the form of high purity Ga particle grains. All the dopant sources were precisely weighed and put into the crucible together with the Si nuggets before fusion. The crystal was grown with an average pulling rate of v = 0.7 mm/min, a crystal rotation rate of wcryst. = 1.26 rad.s -1 and an inverse rotation of the crucible wcruc. = 0.84 rad.s -1 under an Ar pressure of 20 torr. The obtained ingot was then shaped into a 125x125 mm² pseudosquare brick and into 240 μm thick wafers which were KOH etched to remove saw damages. SIMS analyses were carried out on wafers originating from different ingot heights in order to determine the doping impurities vertical distribution. 4-points probe resistivity and Hall effect measurements were performed to study the electrical activation of dopants. Minority carrier lifetimes were measured by QSSPC after chemical polishing and iodine ethanol passivation.

Results and discussion

3.1 Doping impurities distribution Figure 1 shows the evolution of resistivity, measured by 4-points probe, along the ingot height. The resistivity is plotted against the solidified fraction fS until 0.85 since 15% of the initial Si melt mass remained in the crucible after crystallisation or was removed with the ingot tail. The measured resistivity is p-type along the whole ingot height and is comprised between 0,5 and 5,5 Ω.cm, this range being very well suited to the fabrication of good performance solar cells. SIMS analyses (Fig. 2) show that this p-type resistivity is to be attributed to the addition of Ga. Indeed, one can observe that P concentration is at each height over that of B, suggesting that without Ga co-doping, the ingot would be fully n-type. To study the axial distribution of doping impurities, the segregation coefficients keff of each impurity as well as their initial concentrations in the melt C0 were extracted by fitting the Scheil's equation ( 1) to the experimental SIMS data,

1 0 ) 1 (      eff k S eff S f C k C ( 1
)
where CS is the concentration in the crystal at the solidified fraction fS. In order to compare those results with theoretical values, B, P and Ga effective segregation coefficients were calculated for the growth conditions used in this work, with the equation ( 2) derived from the BPS theory [START_REF] Burton | [END_REF],

) / exp( ) 1 ( 0 0 0 D v k k k k eff       ( 2
)
where k0 is the equilibrium segregation coefficient, δ is the thickness of the diffusion layer in the Si melt near the solid to liquid interface and D is the diffusion coefficient of the impurity in the liquid phase. All the calculated keff and C0 are reported in Table 1 together with the k0 taken from [10]. The segregation coefficients of B and P extracted by fitting the Scheil's law are in good agreement with the BPS theory. However, the extracted initial concentrations in the melt are slightly below the values expected from the doping sources that were introduced in the Si melt. The segregation coefficient of Ga derived from the fit to the Scheil's law is one order of magnitude higher than that calculated with the BPS theory. Such an increase of the effective segregation coefficient has already been observed in Ga-doped Si crystal co-doped with B [11]. This increase was explained, deriving the Thurmond and Struthers theory [12], by a potential Ga-B pair formation showing a better accommodation to the Si lattice and thus leading to an enhanced Ga incorporation during crystal growth. However, this interaction is unlikely to explain the increase of the Ga segregation coefficient observed here, since it is expected to occur at B concentrations higher than 1×10 19 cm -3 and to reach a maximum keff of only 0.011 for very high B concentrations. Another plausible explanation for the observed difference is the evaporation of Ga from the melt due to the low pressure conditions. Considering in such case that the quantity of evaporated Ga atoms per unit of time is given by ( 3), the Scheil's law has to be modified into the form (4),

L L Ga C S E dt dn    (3) 1 0 ) 1 (         C L eff S v S E k S eff S f C k C (4)
E being the evaporation rate, CL the Ga concentration in the melt, SL the free surface of the Si melt and SC the surface of the crystal transversal section. The parameters obtained using this modified Scheil law (Fig. 3) are shown in Table 2. Table 2 keff and C0 used for the fit of Ga distribution to the Scheil's law with and without taking evaporation into account.

C0 added (cm -3 )

C0 Scheil (cm -3 ) keff BPS keff Scheil E (cm/s)
Fit without E 5.79×10 18 3.54×10 17 0.0086 0.105 0

Fit with E 5.79×10 18 4.63×10 18 0.0086 0.0086 2.7×10 -5

Doping impurities activation

In order to evaluate the activation of dopants in this material, we executed Hall effect measurements on samples taken out from different ingot heights. The aim was to compare the majority carrier density to the net doping deduced from SIMS analyses. The difficulty to interpret Hall effect measurements arises from the fact that the measured Hall carrier density pH differs from the true carrier density p by a factor close to unity (4), called the Hall factor rH, which depends on the main scattering mechanism.

H H p r p   (5)
rH dependence with doping impurities concentration and temperature have already been studied for uncompensated p-type Si [13,14]. Although no data was found in literature for rH in highly compensated Si, recent results suggest that it is only weakly depending on the compensation level [15]. In this work, rH will be assumed to be equal to 0.7 at ambient temperature for all the samples, as predicted for uncompensated p-type Si of similar doping impurities concentration. The validity of this assumption is discussed below. Fig. 4 displays the evolution of the carrier density p together with the net doping NA-ND deduced from SIMS, both plotted against the solidified fraction. These results show that the measured carrier density p is systematically higher than the net doping NA-ND. One has to remember here that SIMS measurements provide the total dopant concentration, which can differ significantly from the activated dopant concentration located in substitutionnal position. Two hypothesis can explain the larger value of p compared to NA-ND: 1) the true Hall factor is smaller than that used to deduce p from pH; 2) a certain quantity of P atoms are electrically inactive (interstitial position, precipitates or complexes formation).

The resistivity measurements (Fig. 1) allow us to prefer the hypothesis 2). Indeed, assuming all dopants to be activated implies the carrier density to be inferior to NA-ND and thus the majority carrier mobility to be superior to a minimum μmin calculated by [START_REF] Enebakk | 3rd Intenational Workshop on Crystalline Silicon Solar Cells[END_REF]: [START_REF] Kraiem | Proceedings 34th IEEE Photovoltaic Specialists Conference[END_REF] Such a calculation leads to some values for μmin up to 690 cm²/V.s, values that are far over the theoretical limits predicted for uncompensated p-type Si of similar doping levels. It is thus very probable that the actual net doping, taking into account only the electronically activated dopants, is higher than that deduced from SIMS measurements, suggesting that 10 to 15% of the P atoms are inactive. As the P concentration in this material is significantly lower than the solubility limit, it can be inferred that the deactivation of P atoms is to be attributed to complexes or precipitates involving the contribution of other dopants or impurities.

      q N N D A ) ( 1 min

3.3

Lifetime measurements QSSPC lifetime measurements were carried out to evaluate the impact of doping impurities on the electrical quality of Si. Results show that relatively high lifetimes are obtained in the top part of the ingot considering the quality of the Si feedstock. A continuous decrease of the minority carrier lifetime is observed with the increase of the solidified fraction (Fig. 5) which can be attributed, as proposed by Dubois et al. [START_REF] Dubois | [END_REF], to the increase of the majority carrier density and thus of the recombination strength of defects' energy states in the band gap of the Si. The very low carrier lifetimes measured in the bottom part of the ingot may also be attributed to the apparition of high densities of dislocations in this part of the ingot as was observed during crystal growth by the disappearance of the 4 facets on the edge of the ingot, representative of the dislocation-free <100> oriented mono-crystal. 3 Conclusion Gratefully to Ga co-doping we have managed to control the net doping along the ingot height during crystallisation of B and P compensated Si. This method may thus be used to increase the B and P concentration limits allowed in UMG-Si to produce high quality solar cells. The axial distribution of doping species was studied and showed a good agreement with theory. A slight adjustment of the Ga segregation profile can be attributed to an evaporation of Ga during crystal growth with a rate calculated to be equal to E = 2.7×10 -5 cm/s. It was shown that 10 to 15% of the P atoms are electrically inactive. Solar cells are currently under processing in order to evaluate the photovoltaic potential of this material.

Figure 1

 1 Figure 1 Resistivity averaged on 4-points probe measurements done at 5 different points of each wafers taken out from different

Figure 2 B

 2 Figure 2 B, P and Ga distributions along the ingot height measured by SIMS. The continuous lines represent the fit to the Scheil's law.

Figure 3

 3 Figure 3 Axial distribution of Ga fitted with the Scheil's law taking into account evaporation from the melt. Dashed line displays the distribution of Ga expected from the Scheil's law without taking evaporation into account.

Figure 4

 4 Figure 4 Evolution of the carrier density p measured by Hall effect and of the net doping NA-ND measured by SIMS with ingot height.

Figure 5

 5 Figure 5 Evolution of the minority carrier lifetime measured for each sample at the injection level Δn = 0.1×p. The samples surfaces are passivated using an iodine ethanol solution.

Table 1

 1 Concentrations of B, P and Ga added to the Si melt before crystallisation.

	Element	C0 added to the melt (cm -3 )	C0 from Scheil (cm -3 )	k0	keff from Scheil	keff from BPS
	Boron	2.67×10 16	1.99×10 16 0.80 0.771 0.818
	Phosphorus 1.78×10 17	1.58×10 17 0.35 0.316 0.367
	Gallium	5.79×10 18	3.54×10 17 0.008 0.105 0.0086

Acknowledgements

The authors would like to thank SILTRONIX for the crystallisation of the Cz-Si ingot that was used in this work and Maxime Levaslot from IUT Lannion for the Hall effect measurements.