A Lemaire

A Perona

M Caussanel

A Dollet

Open-Circuit Voltage Decay Simulations on Silicon and Gallium Arsenide p-n Homojunctions: Design Influences on Bulk Lifetime Extraction

Keywords:

Open-Circuit Voltage Decay (OCVD) method was investigated for few decades as a simple and convenient method to characterize effective lifetime into Ge and Si p-n homojunctions. Minority carrier lifetime (MCL) is an important parameter to optimize device design where being able to investigate one type of carriers is an important goal. The p-n homojunction design is of major importance to reach that purpose by OCVD. We carried out Technology Computer Aided-Design simulations of the OCVD signal. The study focused on Si and GaAs p-n homojunctions. We looked into bulk thickness and doping level influences on bulk lifetime extraction. Previously, those influences have not been quantified for Si and never investigated for GaAs. MCL accurate extraction from bulk required a bulk thickness at least 4 times higher than diffusion length and emitter doping levels at least 2 orders of magnitude higher than bulk. Likewise, the paper shows different accuracy of extraction between GaAs/Si and p-type/n-type bulk in p-n homojunction.

Introduction

Minority carrier lifetime (MCL) is the most relevant parameter to characterize and optimize electronic and optoelectronic devices such as PV cells [START_REF] Schroder | Semiconductor material and device characterization[END_REF]. Lifetime depends on material nature and on various treatments the material or device has been subjected to. Downstream, lifetime allows to know material purity and gives a follow-up of process steps quality (etching, passivation, metallization and annealing). Upstream, it is possible to investigate new materials bulk properties. So it enables to optimize design such as bulk thickness and doping levels of a solar cell. Measured lifetime is always the contribution of various and independent electrons-holes recombination mechanisms which are SRH (Shockley-Read-Hall), Auger and radiative type [START_REF] Schroder | Carrier lifetimes in silicon[END_REF]. They occur on either sides of the p-n homojunction and at contact interfaces. Therefore measured lifetime is always an effective lifetime and thus usually represented by the simple Eq. 1 and 2 where eff is the effective lifetime, SRH the SRH lifetime, Aug the Auger lifetime, rad the radiative lifetime, b the bulk lifetime and s the surface lifetime.

1 eff = 1 b + 1 s ( 1
)
1 b = 1 SRH + 1 Aug + 1 rad (2)
Since extracted lifetime is always effective, bulk MCL investigation is tough. We decided to investigate the Open-Circuit Voltage Decay (OCVD) method [START_REF] Schroder | Semiconductor material and device characterization[END_REF] because of its conveniences. First it is a fast and non-destructive method. Then it requires a simple circuit and enables characterization of full devices. Finally, it allows extraction of multiple parameters ORCID(s):

such as back surface recombination velocity [START_REF] Rose | Minority-carrier lifetime measurements on silicon solar cells using I sc and V oc transient decay[END_REF] (combination with short-circuit current decay), level of injection [START_REF] Green | Solar cell minority carrier lifetime using open-circuit voltage decay[END_REF][START_REF] Jain | Open circuit voltage decay in p-n junction diodes at high levels of injection[END_REF], series resistance [START_REF] Green | Minority carrier lifetimes using compensated differental open circuit voltage decay[END_REF], diffusion length [START_REF] Stutenbaeumer | Comparison of minority carrier diffusion length measurements in silicon solar cells by the photoinduced open-circuit voltage decay (OCVD) with different excitation sources[END_REF] and capacitance [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF] (corresponds to the voltage when space charge region (SCR) capacitance equals diffusion capacitance). On the other hand, it requires a junction (p-n or Schottky). Therefore an appropriate design (thicknesses and doping levels) has to be defined, in order to highlight the behavior of one type of charge carrier (electron for p-type bulk and hole for n-type bulk) and so to measure only bulk lifetime in our case. We have studied Si because it has been widely investigated by OCVD. However design effects have never been quantified for this material. On the other hand, OCVD studies on GaAs are absent. Only one experimental study has been published [START_REF] Bhimnathwala | Lifetime measurements by open circuit voltage decay in GaAs and InP diodes[END_REF] while it is a rather well-known semiconductor. Actually, III/V materials characterized by OCVD are rather scarce [START_REF] Zemel | Carrier lifetime in InP/InGaAs/InP by opencircuit voltage and photoluminescence decay[END_REF][START_REF] Dalapati | Effect of temperature on the intensity and carrier lifetime of an AlGaAs based red light emitting diode[END_REF]. The aim of this work is to define a p-n homojunction design in order to obtain a better usable OCVD signal. We have worked with TCAD Sentaurus which is a powerful tool to simulate design of devices and circuits. In this paper, influences of emitter doping concentration and bulk thickness of Si and GaAs p-n homojunctions on bulk MCL extraction by OCVD are evaluated.

OCVD Theory

OCVD theory was first introduced in 1955, by Gossick [START_REF] Gossick | On the Transient Behavior of Semiconductor Rectifiers[END_REF] and by Lederhandler & Giacoletto [13] to characterize germanium (Ge). The method is based on the interpretation of Open-Circuit Voltage Decay (OCVD) across a p-n junction after an excitation pulse. The excitation pulse is usually obtained by applying a direct constant bias before opening abruptly the circuit with a switch. Initially, a forward current source was employed as excitation and hereafter optical excitation has been also investigated [START_REF] Mahan | Measurement of minority carrier lifetime in solar cells from photo-induced opencircuit voltage decay[END_REF]. Two main circuits can be used: one where switch and characterized diodes 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] Marcelot | Comparison between TCAD simulated and measured carrier lifetimes in CMOS photodiodes using the Open Circuit Voltage Decay method[END_REF] (voltage bias source) and the other where they are in parallel [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). Usually, OCVD signal presents 3 regions as shown in Fig. 1. The region (i) is due to ohmic resistance [START_REF] Green | Minority carrier lifetimes using compensated differental open circuit voltage decay[END_REF]. The region (ii) is the one of interest for MCL extraction and is due to minority charge carriers diffusion/recombination. In the ideal case, region (ii) presents a linear decay at low injection level during periods a few times higher than the lifetime. The inverse slope is used into Eq. 3 to calculate effective MCL, where is the p-n junction ideality factor, B the Boltzmann constant, the temperature, the elementary electric charge and d /d the linear decay slope. Ideality factor depends on the injection level and on the dominant recombination mechanism [START_REF] Sze | Physics of semiconductor devices[END_REF]. In the region (iii), the voltage variation slows down due to the SCR capacitance [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF].

eff = -B 1 d ∕d (3)
( ) = 0 e -∕ eff (4)
( ) = B ln 1 + ( ) 0 (5)
Eq. 3 is obtained after few assumptions [START_REF] Gossick | On the Transient Behavior of Semiconductor Rectifiers[END_REF][START_REF] Lederhandler | Measurement of minority carrier lifetime and surface effects in junction devices[END_REF]. First it was assumed lifetime is constant through the bulk. It means p-n junction must be close to an abrupt junction and that recombination centers are homogeneously distributed. Thus density of excess minority carriers decays exponentially with a time constant eff [START_REF] Nosov | Switching in Semiconductor Diodes[END_REF] represented by Eq. 4 where n 0 represents electrons density at the SCR boundary just before the opencircuit condition. Then, the potential across a p-n junction depends on the logarithm of minority carriers density [START_REF] Sze | Physics of semiconductor devices[END_REF] represented by Eq. 5 where n 0 is electron density at thermal equilibrium. Hence it is clear why the theory pretends to obtain a linear region during the decay. Next, the applied voltage has to be much higher than the thermal voltage ( a ≫ V ℎ ) to obtain a straight region. Finally, low injection condition must be fulfilled otherwise a second linear region appears and the ideality factor shall be adjusted in consequence (assumed to be around 1 in low injection when SRH recombination is dominant). Experimentally, global OCVD decay can either slow down or speed up according to recombination into the emitter [START_REF] Green | Solar cell minority carrier lifetime using open-circuit voltage decay[END_REF] and the SCR [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF], to SCR capacitance (also called depletion capacitance) (cf. Fig. 1 region (iii)) [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF], high injection level [START_REF] Jain | Open circuit voltage decay in p-n junction diodes at high levels of injection[END_REF], low shunt resistance [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF] and opencircuiting element [START_REF] Bellone | An Analog Circuit for Accurate OCVD Measurements[END_REF]. In transient regime, the region (ii) is due to minority charge carriers diffusion/recombination and is represented by the so-called "diffusion capacitance" [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF] while region (iii) is due to the "SCR capacitance" driven by majority charge carriers. According to AC equivalent circuit [START_REF] Salach-Bielecki | Influence of junction parameters on the open circuit voltage decay in solar cells[END_REF] they are in parallel. In consequence linear part is often blurred and so makes MCL extraction tough and unreliable.

Simulation details

Simulations have been carried out to optimize p-n junction design for bulk MCL extraction by OCVD. The simulated circuit is illustrated by Fig. 2 and represents "series" circuit. This circuit has been selected to perform simulations because of its more convenient configuration and its abundance into literature [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] Marcelot | Comparison between TCAD simulated and measured carrier lifetimes in CMOS photodiodes using the Open Circuit Voltage Decay method[END_REF]. 2D simulations have been performed with TCAD Sentaurus software (Synopsys) using Sentaurus Device "mixed-mode" [START_REF]Device Simulation[END_REF]. It allows to simulate an entire circuit while modeling one or several components with the finite element method. The characterized p-n junction, so-called D1 (cf. Fig. 2), has been 2D designed (cf. Fig. 3) using Sentaurus Structure Editor (sde) while other components such as MOSFET and voltage sources (cf. Fig. 2) have been selected from well-known SPICE compact models provided by TCAD Sentaurus.

The MOSFET is used to switch on and off the circuit (to reach open-circuit condition across the p-n junction). It is controlled by a pulse generator whereas the p-n junction is forward biased by a dc source. MOSFET commutation time has been set up according to our experimental bench (24 ns). Characterized p-n junctions have a constant and abrupt doping profile. Boron (p-type) / phosphorous (n-type) and beryllium (p-type) / tellurium (n-type) are selected doping species for Si and GaAs respectively. Bulk doping concentration has been kept equal to 10 15 cm -3 . Naming emitter the emitter thickness and emitter the diffusion length into the emitter, emitter / emitter ratio has been selected arbitrarily and kept equal to 0.16. Samples width is 500 µm. Voltage source has been set up to ensure injection condition into the ohmic regime: it corresponds to 10 -5 A for Si, 10 -4 A for GaAs p-type bulk and 10 -2 A for GaAs n-type bulk.

Poisson and continuity (electrons and holes) equations have been solved using Boltzmann statistics while the Drift-Diffusion model has been chosen as the charge carrier transport model. No generation has been considered. SRH, radiative and Auger recombination models have been selected for all simulations. SRH recombination rate is represented by Eq. 6 where i is intrinsic charge carrier density, 1 and 1 statistical factors representing the electron and hole densities for the Fermi level equal to the trap energy level. Indexes "n" and "p" always stand for electrons and holes respectively. We have considered the worst case which is when trap energy level is right in the middle of the band gap. SRH recombination doping dependence follows the Fossum law [START_REF] Fossum | Computer-aided numerical analysis of silicon solar cells[END_REF] according to Eq. 7 and 8. Radiative recombination rate is represented by Eq. 9 where rad , the radiative coefficient, is a constant. The TCAD radiative coefficient for Si is by default 0. The Auger recombination rate is represented by Eq. 10 where n,Aug and p,Aug are electron and hole Auger coefficients respectively. The TCAD Auger coefficients for Si is by default doping dependent while it is not for GaAs. All recombination parameters used for simulations are listed in the Tab. 1.

SRH = -i 2
p (TCAD) ( + 1 ) + n (TCAD) ( + 1 ) (6)
n (TCAD) = n,max 1 + ( A ∕ ref ) (7)
p (TCAD) = p,max 1 + ( D ∕ ref ) (8) rad = rad ( -i 2 ) (9) Aug = ( n,Aug × + p,Aug × )( -i 2 ) (10)
Electrical contacts are ohmic and insure an infinite surface recombination velocity (SRV) [START_REF] Aberle | Surface passivation of crystalline silicon solar cells: a review[END_REF]. Actually, TCAD considers electrons and holes with finite maximum velocities ( n = 2.573 × 10 6 cm s -1 and p = 1.93 × 10 6 cm s -1 ) due to thermodynamic considerations. Simulation aim is to extract lifetime through the OCVD signal. Here OCVD response obtained by TCAD simulations has been analyzed according to the following procedure. First step consists in identifying the linear region. Therefore derivative is plotted and so linear region is highlighted through extremum. Then corresponding value is injected into Eq. 3 to get the so-called extracted effective lifetime eff (extracted) . The local ideality factor ( ) [START_REF] Mcintosh | Lumps, Humps and Bumps: Three Detrimental E ects in the Current{Voltage Curve of Silicon Solar Cells[END_REF] is calculated thanks to the dark I-V simulated curve through the following equation:

( ) = B d d (ln ) (11)
where d /d(ln ) is the slope of the ln(I)-V curve at the voltage . In our case, the voltage is where index "ext" stands for extremum. It actually represents the voltage of the extremum identified on the derivative OCVD signal. Finally, extracted value is compared to bulk lifetime value used by TCAD Sentaurus to carry out the simulations and called here TCAD bulk lifetime bulk(TCAD) . Extraction accuracy is directly assessed from this comparison.

Simulation results

Assumptions

OCVD technique provides effective lifetime. Therefore assumptions must be expressed to investigate the behavior of only one type of charge carrier. Following assumptions have been made to extract lifetime from the bulk. First, bulk thickness must be much longer than diffusion length ( bulk ) in order to ensure recombination within and to avoid surface recombination at bulk electrical contact. Secondly total current flowing through the diode must be mainly due to the type of carrier under investigation (electrons or holes). Therefore an asymmetric p-n junction is needed ( emitter ≫ bulk ). Finally, emitter thickness must be much smaller than the diffusion length ( emitter ) to prevent recombination within.

Bulk thickness influence

This section studies bulk thickness influence. Bulk is the layer where minority charge carriers under investigation recombine. As previously introduced, its thickness must play an important role to maximize recombination while avoiding surface recombination at bulk electrical contact. Because contacts are considered ohmic, it is of major importance to prevent charge carriers from reaching them. In the reverse case, surface recombination would be disruptive. Indeed, minority carrier density decays exponentially from SCR region through bulk thickness [START_REF] Mathieu | Physique des semiconducteurs et des composants électroniques-6ème édition: Cours et exercices corrigés[END_REF]. According to exponential decay law, if the bulk thickness is only once the diffusion length 63 % of minority charge carriers should recombine, 99 % if it is 5 times. Thus OCVD signal has been simulated on abrupt silicon (Si) and gallium arsenide (GaAs) p-n homojunctions with various bulk thicknesses. Aug [cm 6 s -1 ] 2.9 × 10 -31 -1.3 × 10 -30 1.03 -9.52 × 10 -31 10 -30 10 -30 ). The presence of 3 distinguishable regions is highlighted in Fig. 4[a] and 5[a]. The region (i) shows a drastic voltage drop until 520 mV for Si and 940 mV for GaAs due to series resistance [START_REF] Green | Minority carrier lifetimes using compensated differental open circuit voltage decay[END_REF]. This voltage is called i-ii because it represents the voltage transition between the regions (i) and (ii). Due to the difference of their energy band gap, i-ii differs for Si and GaAs. The region (ii), which is the region of interest, presents a rather linear decay until about 300 mV for Si and 800 mV for GaAs due to minority carriers recombination. This voltage is called ii-iii and has been clearly explained by Mahan and Barnes [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF]. Once again, this difference between Si and GaAs is expected because they have such a different lifetime. It means diffusion capacitance of GaAs is much lower than Si one, while they have rather the same SCR capacitance. The region (iii) proves a drastic slow down due to this SCR capacitance [START_REF] Mahan | Depletion layer effects in the open-circuitvoltage-decay lifetime measurement[END_REF]. Since total recombination rate is roughly SRH recombination rate for HI and LI, SRH recombination can be considered as the dominant mechanism influencing minority charge carrier behaviour into the bulk. SRH recombination theory exposed in 1952 by Shockley Read [START_REF] Shockley | Statistics of the recombinations of holes and electrons[END_REF] and Hall [START_REF] Hall | Electron-hole recombination in germanium[END_REF] demonstrates that in low injection level bulk lifetime depends only on one charge carrier type (cf. Eq. 12) while in high injection level it depends on both charge carrier type (cf. Eq. 13).

rad [cm 3 s -1 ] 0 2 × 10 -10
bulk SRH (LI) = n [START_REF] Gossick | On the Transient Behavior of Semiconductor Rectifiers[END_REF] bulk SRH (HI) = n + p [START_REF] Lederhandler | Measurement of minority carrier lifetime and surface effects in junction devices[END_REF] These two distinguishable regions are more visible for GaAs and tend to fuse when bulk thickness increases. The higher the thickness the longer the linear decay or plateau duration.

It signifies the extraction becomes more accurate. We have defined plateau duration Δ plateau as time when the d /d extremum varies about ±5 %. Indeed, for Si it is 19.2 µs (≈ 2.1 bulk(TCAD) ) while it is about 71.9 ns (≈ 1.6 bulk(TCAD) ) for GaAs. We observe that extrema position ( extremum ) shifts toward longer times when thickness increases. This delayed phenomenon can be related to the increase of circuit RC constant because of bulk resistance increasing. The main differences between Si Fig. 4 and GaAs Fig. 5 are voltages of transition regime ( i-ii and ii-iii ) and scale of the total decay time. It is shorter for GaAs than for Si. Indeed, bulk lifetime is much lower for GaAs. Fig. 6 shows extracted effective lifetime to TCAD bulk lifetime ratio plotted over bulk thickness to bulk diffusion length ratio for Si and GaAs p-n homojunctions where bulk diffusion length is constant for a given doping level. These two normalized axis allow to compare any materials. It shows results for Si (circle) and GaAs (triangle) p-n junctions in p-type (full symbols) and n-type bulk (empty symbols). Extracted effective lifetimes have been calculated for LI so we have considered the second extremum from derivative d /d Fig. 4[b],[c] and5[b],[c]. In both Si and GaAs case, extracted effective lifetimes increase with the bulk thickness and reach a maximum value for bulk > 4.1 bulk . When thickness is smaller than diffusion length it clearly shows surface recombination influence at bulk electrical contact. The higher the minority charge carrier density at the contact the lower the surface lifetime. Moreover this maximum value gets closer to the TCAD bulk lifetime with the exception of electrons lifetime in GaAs. Indeed, MCL is overestimated around 22.5 % while it is underestimated around 30 % for holes lifetime in both materials and 24 % for electrons lifetime in Si. This trend could be attributed to the higher mobility imbalance between holes and electrons in GaAs.

M i n o r i t y c a r r i e r l i f e t i m e o f : e l e c t r o n s i n p -S i h o l e s i n n -S i e l e c t r o n s i n p -G a A s h o l e s i n n -G a

Doping asymmetry influence

In this section is investigated doping asymmetry influence in p-n homojunction. As previously introduced, a heavily doped emitter is required to extract bulk lifetime accurately. A p-n + or n-p + homojunction configuration is essential to extract accurately by OCVD electrons or holes lifetime respectively. This configuration should lower the weight of minority charge carriers which are not under investigation.

Figure 8:

[a] Simulated Open-Circuit Voltage Decay signals, [b] OCVD signal derivatives plotted over time in GaAs p-n homojunctions for emitter (n-type) doping level equal to 10 15 cm -3 up to 10 18 cm -3 and bulk (p-type) doping level equal to 10 15 cm -3 . bulk ≈ 4.1 bulk .

p-n homojunctions respectively. Bulk doping level is kept constant and equal to 10 15 cm -3 while emitter doping level varies from 10 15 cm -3 up to 10 18 cm -3 . Bulk thickness over bulk diffusion length is about 4.1. We still observe 3 different regions in Fig. 7[a] and 8[a]: regions (i) still show a drastic voltage drop where i-ii is roughly 560 mV and 960 mV for Si and GaAs respectively. It increases slightly with emitter doping level because of built-in voltage rising; regions (ii) also present a rather linear decay due to minority charge carriers diffusion/recombination where ii-iii goes up from 300 mV to 350 mV for Si and from 800 mV to 850 mV for GaAs. This is due to the asymmetric junction which reduces the absolute SCR thickness and thus faintly increases the SCR capacitance. Therefore it reduces the region (ii) duration and hence the extraction accuracy but it is actually compensated by the rise of i-ii . Two extrema are appearing on Fig. 7[b] and 8[b] by increasing emitter doping level. The first is due to HI and the second to LI. This effect is more evident for GaAs. Therefore it reduces LI plateau length in the case of GaAs. Indeed, Δ plateau varies with emitter doping from 71.9 ns to 63.3 ns for GaAs and from 19.2 µs to 24.1 µs for Si. This difference could be attributed to higher injection level in the case of GaAs homojunction. Thus injection level is an important parameter because the longer the plateau duration the better the lifetime extraction. Derivative plateau is shifted to a longer time for both Si and GaAs. However there is no significant emitter doping level influence for emitter ∕ bulk >10 2 in Fig. 7 and8. Fig. 9 shows the extracted effective lifetime to TCAD bulk lifetime ratio plotted over bulk thickness to bulk diffusion length ratio where main results are summarized in Tab.2. Local ideality factors extracted to calculate eff(extracted) through Eq. 3 are given in Tab.2. They are rather close to unity for Si while they are higher for GaAs. It expresses the level injection moving which is more prominent in the case of GaAs. Fig. 9[a],[c] shows results for Si and [b], [d] for GaAs p-n homojunctions. Moreover full symbols represent p-type bulk (Fig. 9[a] and [b]) while empty symbols n-type bulk (Fig. 9[c] and[d]). Emitter doping level has been increased by several orders of magnitude from 10 15 cm -3 to 10 18 cm -3 while bulk doping level has been kept constant and equal to 10 15 cm -3 . Fig. 9 analysis points out that all extracted effective lifetimes for both materials increase and tend to reach a constant value with increasing bulk thickness. Furthermore we can clearly see the influence of doping asymmetrization. and so more accurate when doping asymmetry increases. Actually, there is no significant variation for more than 2 orders of magnitude between emitter and bulk doping levels in Si. This put forward the minimum doping ratio required to get wished asymmetry effect. However, there are tiny nuances for lifetime extraction in GaAs. On the one hand, there is still a significant variation between emitter ∕ bulk =10 2 and 10 3 for holes lifetime extraction Fig. 9[d]. Moreover, extractions for bulk < 1.3 bulk have not been carried out because linear regions were not identifiable. On the other hand, there is no massive effect of the asymmetry when emitter ∕ bulk >10 1 for electrons lifetime extraction Fig. 9[b]. This discrepancy is due to the higher imbalance mobility between holes and electrons in GaAs ( n ∕ p =18) than in Si ( n ∕ p =3). Indeed, this induces that GaAs p-n junction is intrinsically asymmetric when emitter = bulk . Therefore, the higher the imbalance mobility, the higher the required asymmetric ratio.

Lifetime extraction becomes closer to TCAD bulk lifetime

Conclusion

OCVD simulations have been performed with TCAD Sentaurus to investigate design influence on bulk lifetime extraction. Extraction has been carried out by identifying linear voltage decay through d /d evolution. Bulk thickness and doping levels asymmetry between the bulk and the emitter have been considered as the major influence design parameters through study of Si and GaAs p-n homojunctions. Electrons (in p-type bulk) and holes (in n-type bulk) lifetime extractions have been carried out successfully in Si and GaAs. OCVD simulations results enable to extract the bulk lifetime with the accuracy of -6.5 % and -5.3 % in Si homojunction while +47.3 % and -2.3 % in GaAs homojunction for electrons and holes respectively. This results have been obtained with the help of simulated dark I-V curves in order to get the local ideality factor. These results have been reached with a ratio bulk ∕ bulk ≈ 4.1, emitter ∕ bulk =10 3 and emitter ∕ emitter ≈ 0.16. Initially it has been illustrated by simulation that the thicker the bulk the better the extraction. First because Δ plateau increases and secondly because extracted lifetime gets closer to TCAD bulk lifetime and so becomes more accurate. Indeed, this configuration avoids surface recombination at bulk contact. Since electrical contacts have been considered ohmic, their OCVD signal contribution is not negligible. Indeed, ohmic contact represents infinite surface recombination velocity. Moreover it has been illustrated that plateau is shifted to a longer time, surely due to circuit RC constant rise. Finally, it has been demonstrated by simulation that emitter doping level should be at least 2 orders of magnitude higher than bulk doping level to get the best extraction accuracy. However, this ratio of 2 may be lower/higher depending on the mobility imbalance ( n ∕ p ) of the semiconductor studied. Difference of MCL extraction accuracy between electrons and holes is lower in Si than in GaAs. Electrons lifetime extraction is overestimated while holes one is underestimated in GaAs. This could be attributed to the higher mobility and SRV discordance between holes and electrons in GaAs. Our study has well confirmed and quantified our assumptions. It has also initiated research on GaAs OCVD behaviour.

Figure 1 :

 1 Figure 1: Typical behaviour of OCVD signal in silicon showing usual regions.

1 Figure 2 :

 12 Figure 2: OCVD "series" circuit simulated by TCAD Sentaurus mixed-mode.

Figure 3 :

 3 Figure 3: Typical 2D structure of the characterized p-n junction designed in TCAD and used for simulations.

Fig. 4

 4 shows [a] OCVD signals, [b] OCVD signals derivatives and [c] normalized OCVD signal derivatives obtained

Figure 4 :

 4 Figure 4: [a] Simulated Open-Circuit Voltage Decay signals, [b] OCVD signal derivatives and [c] Normalized OCVD signal derivatives plotted over time in Si p-n junctions ( emitter = bulk =10 15 cm -3 ) for various p-type bulk thicknesses normalized by electron diffusion length. OCVD signal derivatives have been normalized by their respective extremum value.

  Fig. 4[b] and 5[b] show 2 extrema in region (ii) that appear at 2 different injection levels: high injection (HI) for a short time and low injection (LI) for a longer time. Fig. 4[c] and 5[c] are plotted to see further extrema. OCVD signal derivatives have been normalized by the second extremum.

Figure 5 :

 5 Figure 5: [a] Simulated Open-Circuit Voltage Decay signals, [b] OCVD signal derivatives and [c] Normalized OCVD signal derivatives plotted over time in GaAs p-n junctions ( emitter = bulk =10 15 cm -3 ) for various p-type bulk thicknesses normalized by electron diffusion length. OCVD signal derivatives have been normalized by their respective extremum value.

1 Figure 6 :

 16 Figure 6: Normalized extracted lifetime plotted over normalized bulk thickness in Si (circle) and GaAs (triangle) p-n homojunctions from OCVD simulations. Full symbols represent electrons in p-type bulk and empty symbols holes in n-type bulk. bulk(TCAD) is the bulk lifetime used by TCAD and bulk is the bulk diffusion length.

Fig. 7

 7 and 8 show [a] OCVD signals obtained by simulation and [b] OCVD signal derivatives for Si and GaAs

Figure 7 :

 7 Figure 7: [a] Simulated Open-Circuit Voltage Decay signals, [b] OCVD signal derivatives plotted over time in Si p-n homojunctions for emitter (n-type) doping level equal to 10 15 cm -3 up to 10 18 cm -3 and bulk (p-type) doping level equal to 10 15 cm -3 . bulk ≈ 4.1 bulk .

  e s l i f e t i m e ( n -b u l k ) N e m i t t t e r / N b u l k c t r o n s l i f e t i m e ( p -b u l k )

Figure 9 :

 9 Figure 9: [a] [b] normalized extracted lifetime plotted over normalized bulk thickness in Si p-n homojunction and [c] [d] in GaAs p-n homojunction from OCVD simulations. [a] [c] represent electrons (full symbols) in p-type bulk and [b] [d] holes (empty symbols) in n-type bulk. bulk(TCAD) is the bulk lifetime used by TCAD and bulk is the bulk diffusion length. Each symbol type stand for specific emitter doping level.

Table 1

 1 Summary of recombination model parameters implemented in TCAD for Si and GaAs.

		Si		GaAs
	charge carriers	electrons	holes	electrons	holes
	max [s]	10 -5	3 × 10 -6	5 × 10 -8	5 × 10 -8
	ref [cm -3 ]	10 16	10 16	10 16	10 16

Table 2

 2 Summary of extracted results from TCAD simulations with bulk ∕ bulk ≈ 4.1 and emitter ∕ bulk =103

			bulk type	Si	GaAs
	bulk [µm]	n-doped p-doped	730 230	126 29
	bulk [µm]	// //	179 57	31 7
	bulk(TCAD) [s]	// //	9.09 × 10 -6 2.73 × 10 -6	4.55 × 10 -8 4.55 × 10 -8
	eff(extracted) [s]	// //	8.49 × 10 -6 2.58 × 10 -6	6.70 × 10 -8 4.44 × 10 -8
	error [%]	// //	-6.5 -5.3	+47.3 -2.3
	extremum [ bulk(TCAD) ]	// //	3.89 5.78	3.40 3.98
	Δ plateau [ bulk(TCAD) ]	// //	2.65 2.38	1.39 1.68
	(	) [-]	// //	1.10 1.05	1.32 1.21

Acknowledgment

This work was supported financially by PPCM and COS-MIC teams from PROMES-CNRS laboratory.


A Legarra
email: andres.legarra@toulouse.inra.fr

V Ducrocq

Keywords: genetic evaluation, genomic selection, mixed model Computational strategies for national integration of phenotypic, genomic, and pedigree data in

published or not. The documents may come

Computational strategies for national integration of phenotypic, genomic, and pedigree data in a single-step best linear unbiased prediction

Andres Legarra, Vincent Ducrocq

INTRODUCTION

Genetic evaluations in dairy cattle nowadays include massive SNP marker data. This will be referred to as genomic evaluation hereinafter. Methods include adaptations of BLUP evaluations, sophisticated regressions based mainly on Bayesian models, and nonparametric techniques. Results have generally shown a superiority of genomic evaluations over pedigree-based ones.

In pedigree-based models using BLUP, large data sets are used. Solution of mixed-model equations (MME) usually takes O(n) computing time; that is, central processing unit (CPU) time roughly proportional to n, the number individuals in the pedigree, using iteration on data techniques that reduce storage needs to O(n) elements [START_REF] Schaeffer | Computing strategies for solving mixed model equations[END_REF][START_REF] Misztal | Indirect solution of mixed model equations[END_REF][START_REF] Strandén | Solving large mixed linear models using preconditioned conjugate gradient iteration[END_REF]. Genomic evaluations that use mixed models, such as Bayesian regressions, Lasso, or Elastic Net, can be efficiently computed by several strategies, with an O(pm) cost, where p is the number of genotyped individuals and m is the number of SNP effects [START_REF] Legarra | Technical note: Computing strategies in genome-wide selection[END_REF][START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF][START_REF] Shepherd | Genomic selection and complex trait prediction using a fast EM algorithm applied to genome-wide markers[END_REF]. Iteration on data techniques can reduce storage to O(m) variables.

However, all of these strategies work for models where the own phenotype of the genotyped individual is being analyzed. In dairy cattle, daughter yield deviations (DYD;[START_REF] Vanraden | Derivation, calculation, and use of national animal model information[END_REF] are used.

Unfortunately, bias in computation of DYD resulting from genomic preselection of the best animals may become a problem [START_REF] Mäntysaari | Use of bivariate EBV-DGV model to combine genomic and conventional breeding value evaluations[END_REF]Patry and Ducrocq, 2011a,b). Other problems of DYD are their variable accuracy, difficult definition, or computation for some nonlinear traits (e.g., longevity) or species (e.g., sheep or swine).

A solution to this problem is the one-or single-step method, suggested simultaneously by [START_REF] Misztal | Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information[END_REF], [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF], [START_REF] Aguilar | Hot topic: A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF], and [START_REF] Christensen | Genomic prediction when some animals are not genotyped[END_REF]. We propose here to call this method single-step genomic BLUP (SSGBLUP). The method uses the property that, at a given locus, the covariance of the gene frequency (g = {0, 1/2, 1}) for the ith allele of 2 individuals j and k is described by the expression Cov(g ij , g ik ) = θ jk p i (1 -p i ) [START_REF] Cockerham | Variance of gene frequencies[END_REF], where θ jk is the coancestry of the 2 individuals. Thus, the (unobserved) genotype of an individual can be modeled using the numerator relationship matrix [START_REF] Gengler | A simple method to approximate gene content in large pedigree populations: Application to the myostatin gene in dual-purpose Belgian Blue cattle[END_REF][START_REF] Christensen | Genomic prediction when some animals are not genotyped[END_REF]. This will equally be true for any linear function of genotypes (such as additive breeding values, which are a linear function of genotypes at causal loci). Therefore, the expression can be applied to breeding values, markers, or haplotypes, and warrants generality of the derivations of SSGBLUP to models based on markers, haplotypes, or even nonparametric genomic covariances.

Consider the model y = Xb + Wu + e, [1] where y is a vector of phenotypes at a trait; b, u, and e are fixed effects, random effects, and residuals, respectively; and X and W are incidence matrices.. For clarity and without loss of generality, we will consider here that R I = σ e 2 , where R is the covariance matrix of e, I is an identity matrix, and σ e 2 is the residual variance.

The MME for single-trait SSGBLUP are

′ ′ ′ ′ +                     = ′ ′          - X X X W W X W W H b u X y W y α u 1  , [2]
with terms as previously defined, and

H A A A A G A - - - = + -           1 11 12 21 22 1 22 1 ,
where A and G are, respectively, pedigree-based (numerator) and genomic relationship matrices, partitioned in genotyped and ungenotyped individuals,

σ u Var u 2 =
( ) is the polygenic genetic variance and α σ σ

u e u
= 2 2 . The model is easily extended to a multiple trait or a model including maternal or random regression genetic effects by including a more general structure in R, substituting σ u 2 by a covariance matrix, and using Kronecker products. The G matrix is usually formed as in [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF]; that is, G = ZDZ′, where Z is a centered incidence matrix of SNP, D i,i = m/2p i q i or D i,i = 1/2∑p i q i , and p i is the allelic frequency of the ith (i = 1 . . . m) locus. An improved G may require an appropriate scaling of G so that G and A 22 are proportional to the same variance component σ u 2 [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] Vitezica | Bias in genomic predictions for populations under selection[END_REF]. Computation strategies for SSGBLUP are described in [START_REF] Aguilar | Hot topic: A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF][START_REF] Aguilar | Efficient computations of genomic relationship matrix and other matrices used in the single-step evaluation[END_REF]. The cost of inversions of G and A 22 is a cubic function O(p 3 ) of the number of genotyped animals (p).

In our opinion, 3 major barriers for widespread use of SSGBLUP are (1) this (cubic) computing cost; (2) the need to modify the existing software, which covers the wide variety of evaluation models currently used (including heterogeneity of variances, multiple-trait models, threshold models, test-day models, and survival models, among others); and (3) its linear form, which precludes the use of nonlinear estimators including genomic information. We address here 2 of these 3 problems. We show how to write 2 sets of equations equivalent to Equation 2, computable in O(n) + O(pm), without explicit inversion and we present an alternative derivation of Equation 2, leading to an iterative solution procedure that allows the use of already written BLUP software, which cannot otherwise incorporate genomic information. Both allow iteration on data techniques. As for the third problem, a partial solution is to precompute an improved G matrix [START_REF] Zhang | Best linear unbiased prediction of genomic breeding values using a trait-specific marker-derived relationship matrix[END_REF][START_REF] Legarra | Improved Lasso for genomic selection[END_REF]. We outline a more general procedure to include Bayesian estimators, based on the iterative approach.

MATERIALS AND METHODS

We propose to solve Equation 2 by iteratively considering 2 systems: the first one as close as possible to the regular MME, apart from a correction of the right-hand side, and the second one isolating the block corresponding to genotyped animals. In other words, an iterative solution of Equation 2 is obtained, alternating between a regular genetic evaluation and a regular genomic evaluation, whatever the model and the methodology chosen for either one. Such strategy should lead to only minimal modifications of the existing BLUP software used for national evaluations and should be compatible with most genomic evaluation methods.

Two different derivations of these 2 systems will be presented. Although their starting point is quite different, they lead to equivalent results.

Extended SSGBLUP MME Introducing Additional Unknowns

The first technique is the derivation of an equivalent model by including additional unknowns. This form of data augmentation can be seen as the opposite of the well-known absorption technique, and aims at simplifying computations. For instance, [START_REF] Henderson | Best linear unbiased prediction of breeding values not in the model for records[END_REF] showed that including all animals in a pedigree in the set of unknowns allowed the use of his simple rules for inversion of A [START_REF] Henderson | A simple method for computing the inverse of a numerator relationship matrix used in prediction of breeding values[END_REF]. Similarly, including a permanent environment effect in a simple repeatability model leads to an easier solution than explicitly using the residual covariance structure between repeated records of a same animal.

Distinguishing between phenotypes of nongenotyped and genotyped animals, we can write Equation 1 as

y y X X b W 0 0 W u u 1 2 1 2 1 2 1 2           =           +                      +           e e 1 2
.

Regular SSGBLUP. Consider then an extended form of system [2]:

′ ′ ′ ′ ′ + ′ ′ + + X X X W X W W X W W A A W X A WW A 1 1 2 2 1 1 1 1 11 12 2 2 12 2 2 22 α α α α α u u u u u ( ( ) ˆG A b u u X y W y - - -                               = ′ ′ 1 22 1 1 2 1 1 1 2 2 ′               W y , [3]
where, again, subscripts 1 and 2 refer to nongenotyped and genotyped animals, respectively, and X′ = (X 1 ′X 2 ′). Complex models (e.g., random regression or maternal models) would include terms such as W 1 ′W 2 , but this will not be used here (although our derivations are completely general). In Equation 3, u was split in ungenotyped (u 1 ) and genotyped (u 2 ) animals.

Unsymmetric Extended SSGBLUP. An equivalent system of equations is

′ ′ ′ ′ ′ + ′ ′ + X X X W X W 0 0 W X W W A A 0 0 W X A WW A 1 1 2 2 1 1 1 1 11 12 1 2 12 2 2 2 α α α α u u u u 2 2 22 1 α α α α α α u u u u u u I I 0 0 I A 0 0 0 I 0 G b u -                         ˆû X y W y W y 0 0 2 1 1 2 2 - -                       = ′ ′ ′           ϕ γ              , [4]
where φ and γ are new vectors of unknowns, with the same size as u 2 . Vector γ corresponds to vector s in [START_REF] Strandén | Technical note: Derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit[END_REF]. The lack of symmetry of these equations will be used later to describe the convergence of iterative methods to the solution. Absorption of the last 2 equations into the first 3 equations leads to system [3]. Equation 4 can be solved, for instance, by successive underrelaxation (SUR; [START_REF] Broyden | On convergence criteria for the method of successive over-relaxation[END_REF] which is a special case of successive overrelaxation (SOR) with relaxation factors lower than 1; we will return to this point later) or the biconjugate gradient stabilized method (BiCGstab;[START_REF] Van Der Vorst | Iterative Krylov Methods for Large Linear Systems[END_REF][START_REF] Misztal | Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information[END_REF]. If one changes the sign of the bottom row of Equation 4, this creates a symmetric system of equations. In that case, a negative sign is associated with G, which is unusual; this does not prevent γ from being a variable with real domain. It is, in fact, a linear transformation of u; actually,

γ = G -1 u 2 , ϕ = - A u 22 1 2 , with variances Var Var u u γ ϕ ( ) = ( ) = - - G A 1 2 22 1 2 σ σ , .

Extended SSGBLUP MME Introducing Additional Unknowns and Marker Effects

For models where G = ZDZ′, another set of equations explicit on marker solutions ĝ is as follows:

′ ′ ′ ′ ′ + ′ ′ + X X X W X W 0 0 W X W W A A 0 0 W X A WW A 1 1 2 2 1 1 1 1 11 12 1 2 12 2 2 2 α α α α u u u u 2 2 1 22 1 1 α α α α α α u u u u u u I G Z 0 0 I A 0 0 0 Z G 0 D - - ′ -                   - - -        - -                       = ′ ′ ′ ˆb u u g X y W y W y 0 1 2 1 1 2 2 ϕ 0 0                      
.

[5]

These equations show that ĝ = DZ′G -1 û 2 , which is the best linear predictor of marker effects g given u 2 [START_REF] Henderson | Sire evaluation and genetic trends[END_REF][START_REF] Strandén | Technical note: Derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit[END_REF]. Similar systems of equations are being used by some scientists [P. M. VanRaden, US Department of Agriculture, Bethesda, MD, personal communication; D. L. Johnson, Livestock Improvement Corp. (LIC), Hamilton, New Zealand, personal communication; N. Gengler, Université de Liège (ULg), Gembloux, Belgium, personal communication]. A reduced, but possibly more cumbersome, system of equations including u 1 and g (but not u 2 ) is as follows:

′ ′ ′ ′ ′ + ′ ′ ′ ′ ′ X X X W X W Z 0 W X W W A AZ 0 Z W X ZA ZW W 1 1 2 2 1 1 1 1 11 12 2 2 12 2 α α α u u u 2 2 22 1 2 22 1 Z ZA Z D Z 0 0 Z A b u + ′ + ′                     - α σ α α α u e u u u ˆĝ X y W y Z W y 0 --ϕ ϕ                   = ′ ′ ′ ′                  1 1 2 2   .
In contrast to the previous one, the left-hand side of these equations is positive definite and symmetric. Its derivation is presented in Appendix A.

Extended SSGBLUP MME Derived from an Equivalent Model

In model 1, one can decompose u into a strictly polygenic part u* and a deviation due to the genomic information that we will call d = u -u*:

y = Xb + W (u* + d) + e.
[6]

Let us assume (or impose) that Var

u u d A 0 0 G A * .               = -               σ 2
Although G is guaranteed to be semipositive definite and A is positive definite, G -A (although full rank) has no guarantees whatsoever; in fact, more often than not, the eigenvalues of G -A span negative and positive values. [START_REF] Thompson | Sire evaluation[END_REF] and [START_REF] Thompson | Estimating genetic parameters using an animal model with imaginary effects[END_REF] suggested including imaginary effects to deal with negative variances. These equivalent models were found to considerably reduce computing time for estimation of genetic parameters with a multivariate reduced animal model [START_REF] Besbes | Estimation of genetic parameters of egg production traits of laying hens by restricted maximum likelihood applied to a multiple-trait reduced animal model[END_REF]. Developing this idea, in Appendix C (which we suggest reading after Appendix B) we show that mixed model equations are identical whether G -A is positive definite or not, as far as it is invertible. Distinguishing between phenotypes of nongenotyped and genotyped animals, we can write Equation 1as

y y X X b W 0 0 W u d u d 1 2 1 2 1 2 1 1 2           =           +           + + * * 2 2 1 2           +           e e .
For nongenotyped animals, the deviation d 1 due to the genomic information is modeled as a linear function of d 2 , obtained by regressing on genomic contribution d 2 from genotyped individuals:

d A A d 1 1 2 22 1 2 = - . There- fore, y y X X b W 0 0 W u u 1 2 1 2 1 2 1 2           =           +                * *       +           +           - W A A W d e e 1 1 12 22 2 2 1 2 .
From above, u 2 * and d 2 are assumed to be indepen-

dent cov , , * u d 0 2 2 ( ) =       var , u G 2 2 ( ) = σ u and var . d G A 2 2 2 2 ( ) = - ( ) σ u Also, var . d A A G A A A 1 1 2 22 1 22 22 1 21 2 ( ) = - ( ) - -
σ u From this model, MME can be derived (details are provided in Appendix B), and are

′ + + - X X W W X W 0 W X W W A A 0 W X A WW A 1 1 2 2 1 1 1 1 11 12 2 2 21 2 2 22 ' ' ' ' ' ' α α α α u u u u α α α α u u u A 0 0 A A G A 1 1 1 1 22 22 22 22 - - - - - + - ( ) { }                                          = ′             ˆ' ' b u u d X y W y W y 0 1 2 2 1 1 2 2          .
[7]

Journal of dairy Science Vol. 95 no. 8, 2012

STraTegIeS For SIngLe-STeP genoMIc eVaLuaTIon 4633

In Appendix B, we show that the absorption of the last equation of system [7] into the first ones leads to system [3]: models [1] and [6] are equivalent.

The last equation of [12] is

A G A d A u 1 1 22 22 1 22 2 - - - + - ( )         = 2 ˆˆ or ( ) ˆˆˆˆ, * G A d A u d A u 1 1 1 - - - - 22 22 2 2 2 2 2 2 = ( -)=
which becomes, after some manipulations,

ˆˆˆ. d u A G u 2 = - - 2 2 2 1 2 [8] This result implies that ˆˆ. * u A G u 2 2 2 1 2 = - Note also that ˆd I A G Zg 2 2 2 1 = - ( ) - if u 2 = Zg.

Iterative Solution of SSGBLUP MME: Iterative SSGBLUP

Splitting the MME in [4] into 3 parts, one gets the following block iterative solution algorithm:

(1) Solve

′ + +    X X X W X W W X W W A A W X A WW A 1 1 2 2 1 1 1 1 11 12 2 2 21 2 2 22 ' ' ' ' ' ' α α α α u u u u                             = ′         '' b u u X y W y W y 1 2 1 1 2 2          + -               0 0 α α u u ˆφ γ [9]
for ˆ, b û 1 , and û 2 . Equivalently, splitting the MME in [7] leads to the following:

(2) Solve

(1) Solve

′ + +    X X X W X W W X W W A A W X A WW A 1 1 2 2 1 1 1 1 11 12 2 2 21 2 2 22 ' ' ' ' ' ' α α α α u u u u                             = ′         '' b u u X y W y W y 1 2 1 1 2 2          +                 - 0 0 A d 1 α u 22 2 ˆ; [11] (2) Solve A G A d A u 1 1 1 22 22 2 2 2 2 - - - + -       = ( ) ˆˆ
for ˆ, d 2 which becomes, after some manipulations comparable to those leading to [7] and using [8]:

(2′) Compute

A d A G u 1 1 22 22 1 2 - - - = - ( ) = - ˆˆˆˆ, 2 ϕ γ [12]
which is identical to It is well known [START_REF] Coleman | Large Sparse Numerical Optimization[END_REF] that convergence of (block) Gauss-Seidel or SOR for Bx = t, where B is the matrix of coefficients of the MME and x is the vector of unknowns at each iteration, does exist if and only if B is positive definite when B is symmetric and, therefore, [9-10] will converge if [7] is positive definite. It can be shown that this is the case if G -A 22 is itself positive definite. This is rarely the case.

However, a less-known result is that, for B not symmetric, if B+B′ is positive definite (this is the case of [4]) then, a positive, real number w g exists so that for 0 < w < w g , the iterative solution using SOR/SUR with a relaxation parameter w will converge to the solution [START_REF] Broyden | On convergence criteria for the method of successive over-relaxation[END_REF][START_REF] Varga | A note on a perturbation analysis of iterative methods, with an application to the SSOR and ADI iterative method[END_REF]. Usually w g will be lower than 1. Therefore, because G -A 22 is usually not positive definite, updates of iterative SS-GBLUP using [9-10], need to be done by block SUR [START_REF] Broyden | On convergence criteria for the method of successive over-relaxation[END_REF], where new solutions are a weighted average of the solutions to [9-10] and the former solutions. Therefore, SUR iterations should be as follows. At iteration t,

ˆˆ* b u u b u u b u 1 2 1 2 1               =               + - ( ) t ω ω 1 1 2 1 ˆ, u               - t where ˆ* b u u 1 2               is a solution to [9], ˆˆˆ, * ϕ ϕ ϕ t t = + - ( ) - ω ω 1 1 and ˆˆˆ, * γ γ γ t t = + - ( ) - ω ω 1
1 where ˆ* ϕ and ˆ* γ are solutions to [10]. As for the system of Equation 5, this can be implemented iteratively as follows:

(1) Solve

′ + +    X X X W X W W X W W A A W X A WW A 1 1 2 2 1 1 1 1 11 12 2 2 21 2 2 22 ' ' ' ' ' ' α α α α u u u u                             = ′         '' b u u X y W y W y 1 2 1 1 2 2          +               0 0 - α α u u ˆˆ; ϕ γ (

Simulated Example

We used the QMSim simulator [START_REF] Sargolzaei | QMSim: A large-scale genome simulator for livestock[END_REF] to simulate 28,800 individuals, from which 1,800 were genotyped for 30,000 polymorphic loci. The scenario is roughly as in [START_REF] Vitezica | Bias in genomic predictions for populations under selection[END_REF]. Genotyped individuals were the last 2 generations, half of them with performance records and half without. Values of ω = (0.2, 0.5, 0.6, 0.7, 1.0) were tried. In all cases, equations were stored in memory using hash techniques, but conclusions about convergence are the same as for iteration on data techniques. The number of nonzero elements in traditional MME was 214,911. The number of nonzero elements in G and A 22 was 1,800 2 = 3,240,000. Convergence levels were set to 10 -12 in all STraTegIeS For SIngLe-STeP genoMIc eVaLuaTIon 4635 iterative solvers (PCG, BiCGstab, and SUR), although they are not strictly comparable across methods. Table 1 summarizes the different computing options.

RESULTS AND DISCUSSION

Simulated Example

Table 1 reflects the computing times and number of operations for each method. Figure 3 shows the progression of convergences for each method. All methods reached the same solutions (correlation >0.9999) once converged. The iterative SSGBLUP did not converge for ω = 0.6 or higher, whereas convergence for ω = 0.5 took 83 iterations (213 iterations for ω = 0.2). Overall, and according to results in Table 1, all solvers are competitive. However, this will not be so for other numbers of total and genotyped individuals and number of SNP. We will discuss this point later.

The need for integrating all data in genomic evaluation implementations is well acknowledged [START_REF] Vanraden | Invited review: Reliability of genomic predictions for North American Holstein bulls[END_REF][START_REF] Wolc | Breeding value prediction for production traits in layer chickens using pedigree or genomic relationships in a reduced animal model[END_REF]. In addition to being theoretically appealing, implementation of the single-step methodology provides good results in both simulated [START_REF] Christensen | Genomic prediction when some animals are not genotyped[END_REF][START_REF] Vitezica | Bias in genomic predictions for populations under selection[END_REF] and dairy cattle, poultry, and pig real data [START_REF] Aguilar | Hot topic: A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF][START_REF] Chen | Effect of different genomic relationship matrices on accuracy and scale[END_REF][START_REF] Christensen | Genomic predictions on pigs using the single step method[END_REF]. It is crucial to use single step to get unbiased estimates if early selection based on genomic proofs becomes a reality (Patry and Ducrocq, 2011a). Yet, the need for new software as well as a lack of understanding of features of the method seem to somehow hamper its implemen-tation. We herein provide an alternate derivation and techniques (expansion of MME, iterative solutions, and algorithms that can work for iteration on data) that may help practitioners. The expansion was used by [START_REF] Henderson | Best linear unbiased prediction of breeding values not in the model for records[END_REF] to include the sparse inverse of the numerator relationship matrix in the MME. The alternate derivation follows the logic of unknown parent groups (e.g., [START_REF] Quaas | Modified equations for sire models with groups[END_REF], with the added complexity that, in the latter, group effects were added to the founders, whereas here deviations d due to SNP effects are added "in the middle" of the pedigree. We prove in Appendix C that this derivation is correct regardless of whether G -A 22 is positive definite or not.

Distributional Properties

In all symmetric transformations, the distributional properties of [1] are kept. For instance, reliabilities of û can be computed from the inverse of the left-hand side of [4] (in its symmetric counterpart) or [5] or [7]. Further, a Gibbs sampler can be implemented by using [9-10] or [9, 10, and 13]. This Gibbs sampler can eventually update G without the need to calculate its inverse.

Theoretical Convergence

As presented above, let B be the left-hand side coefficients of [4]. Because B + B′ is positive definite, then for some 0 < w < w g , the iterative solution using SUR will converge to the solution [START_REF] Broyden | On convergence criteria for the method of successive over-relaxation[END_REF][START_REF] Varga | A note on a perturbation analysis of iterative methods, with an application to the SSOR and ADI iterative method[END_REF] Niethammer, 1993). Usually w g is lower than 1. Therefore iterative SSGBLUP [9-10], which is a block variant of SUR [START_REF] Coleman | Large Sparse Numerical Optimization[END_REF], will converge to a solution. As for the BiCGstab used in unsymmetric extended SSGBLUP, this is a robust algorithm (van der [START_REF] Van Der Vorst | Iterative Krylov Methods for Large Linear Systems[END_REF]. Note that only (block) stationary methods were considered here. Nonstationary methods (e.g., variants of PCG) could lead to faster convergence.

Different Definitions of G

An overlooked aspect about Equation 2 is that the derivation is general: other matrices can be used instead of the VanRaden (2008) G, as long as they describe additive covariances (not necessarily relationships) among individuals. For instance, in the French dairy cattle evaluation, gametic effects based on haplotype segments are used [START_REF] Guillaume | Short communication: Correlations of marker-assisted breeding values with progeny-test breeding values for eight hundred ninety-nine French Holstein bulls[END_REF]. If the incidence matrix of gametic effects in locus i is Z i , and the covariance matrix among gametic effects is Λ i , then G = κ 0 A 22 + ∑κ i Z i Λ i Z i ′, where κ i are fractions of the total genetic variance.

Computing

In unsymmetric extended SSGBLUP or iterative SS-GBLUP, no explicit inverses of A 22 and G are needed. These are expensive. Computation of A 22 has a cost of O(pn) or O(p 2 ) if a sub-pedigree is extracted [START_REF] Aguilar | Efficient computations of genomic relationship matrix and other matrices used in the single-step evaluation[END_REF]; computation of G has O(mp 2 ) cost. Inverting either matrix costs O(p 3 ). Direct computation of G -1 costs O(mp 2 ) if the [START_REF] Sherman | Adjustment of an inverse matrix corresponding to changes in the elements of a given column or a given row of the original matrix[END_REF] formula is used. This is not tenable for very large p.

However, in such a case, equations can be solved by iteration on data using Jacobi or PCG methods, without explicit computation or storage of either A 22 or G. In both, the kernel and more expensive part of the solving method consists in repeated multiplications of the form t = Bx (in BiCGstab, this computation is done twice per round). For classical MME (i.e.,[9]), the PCG method using iteration on data has been described by [START_REF] Strandén | Solving large mixed linear models using preconditioned conjugate gradient iteration[END_REF], and consists in reading data and pedigree files and adding contributions. Modifications for unsymmetric extended SSGBLUP [4] would consist of adding the A 22 φ and Gγ contributions. The first one can be done, without explicit computing of A 22 , by reading twice the pedigree file, with cost O(n) [START_REF] Misztal | Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information[END_REF]. The second one can be done, if G has a form similar to G = ZDZ′ (e.g., [START_REF] Guillaume | Short communication: Correlations of marker-assisted breeding values with progeny-test breeding values for eight hundred ninety-nine French Holstein bulls[END_REF][START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] Calus | Effects of the number of markers per haplotype and clustering of haplotypes on the accuracy of QTL mapping and prediction of genomic breeding values[END_REF] by computing Gγ = Z[D(Z′γ)] [START_REF] Strandén | Technical note: Derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit[END_REF], at a cost O(mp) if D is diagonal or very sparse, which is usually the case. If mp is very large, the product can be computed by iteration on data. In iterative SSGBLUP, solving G u ˆγ = 2 and A u 22 2 ˆφ = (i.e., [10]) can be achieved by PCG with the same requirements. Indeed, O(mp) is the lowest bound attainable, because this is the cost of reading the raw genotypic data.

In regular SSGBLUP equations (i.e., [3]), inversions of G and A 22 can be avoided as well. The multiplication t = Bx, needed in PCG:

′ ′ ′ ′ ′ + ′ ′ + + X X X W X W W X W W A A W X A WW A 1 1 2 2 1 1 1 1 11 12 2 2 12 2 2 22 α α α α α u u u u u ( ( ) Ĝ A b u u - - -                               1 22 1 1 2
can be obtained by summing

′ ′ ′ ′ ′ + ′ ′ +    X X X W X W W X W W A A W X A WW A 1 1 2 2 1 1 1 1 11 12 2 2 12 2 2 22 α α α α u u u u                             b u u 1 2 and α u G u A u - - - ( ) 1 2 2 2 1 2 ˆˆ,
where ˆˆ, ϕ = when needed). This is equivalent to the cost of running a 2-step genetic evaluation where pedigree-based and genomic-based genetic evaluations are run in succession. In addition, solving of regular MME (as in [9]) can be done using any existing, already optimized, software. Table 2 summarizes the number of operations and memory needed for each system in the cases of iteration on data. The case of regular SSGBLUP without inversion of G and A 22 is as iterative SSGBLUP [1]. It can be seen that, whereas all methods are in principle feasible, its effectiveness for large p depends critically on the number of iterations needed to achieve convergence of global equations and of blocks in the case of iterative SSGBLUP and of regular SSGBLUP without inversion. This deserves further investigation with real, large data sets.

Estimation of SNP Effects

As for Equation 5, this is explicit in SNP effects. Note that ĝ = DZ′G -1 û 2 ; that is, assuming multivariate normality, ĝ is the best predictor of g, given û 2 , or, equivalently, the solution to the system û 2 = Zg, with D the covariance matrix of g and no residual terms. There are ways of solving for g, the marker effects, which do not involve explicit construction or inversion of G [START_REF] Legarra | Technical note: Computing strategies in genome-wide selection[END_REF][START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] Strandén | Technical note: Derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit[END_REF] and also nonlinear methods (BayesB and Bayesian Lasso, among others; [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF][START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] De Los Campos | Predicting quantitative traits with regression models for dense molecular markers and pedigree[END_REF]. Note that even if there is no residual (û 2 = Zg), g has a posterior distribution because g has a prior distribution; therefore, the different g where û 2 = Zg (and, therefore, equally likely) have different a posteriori probabilities. Thus, the same software used for nonlinear genomic evaluation methods can be used for single step, provided it iteratively uses current values of û 2 as input data, possibly (for the sake of practicality) fitting a very small residual variance. Use of nonlinear models for g within [5] will be an approximation because [5] was derived on the pretence of multivariate normality. Matrix G involved in computation of γ = G -1 Zg should be updated if the variance explained by each SNP changes (as in BayesA, BayesB, or the Bayesian Lasso). Whereas individual variances of SNP effects depend much on the prior assumptions, use of estimated G, which is a function of all of them, seems to be a robust alternative [START_REF] Zhang | Best linear unbiased prediction of genomic breeding values using a trait-specific marker-derived relationship matrix[END_REF][START_REF] Legarra | Improved Lasso for genomic selection[END_REF], possibly because the weight of the prior fades out.

A Special Case

Situations exist where it is desired to base genomic evaluations on a reference population for genomic selection differing from the complete list of genotyped animals in the national evaluation. This is the case when there are some phenotypes from genotyped animals that we would like to exclude from the genomic evaluation part: in most European countries and Canada, own records from bull dams are excluded for the estimation of marker effects, because of fear of preferential treatment. Even if cows are not subject to preferential treat-ment, the bias in genotyping (genotyped animals having higher breeding values than the average) might lead to bias of genomic proofs [START_REF] Vitezica | Bias in genomic predictions for populations under selection[END_REF]. On the other hand, there may also be some extra phenotypes from genotyped animals that need to be included in the genomic evaluation part but are not routinely included in regular evaluations. A typical example is phenotypes (usually deregressed international EBV) from foreign bulls in multinational reference populations. A variant of the single-step approach would consist of including a usual genomic evaluation before the second step (i.e., Equation 10) of the iterative procedure to get a û 2 estimate based on the desired information. This estimate would be plugged into [10]. The convergence in practice to a proper solution needs to be studied.

Current Limits

In practice, to our knowledge, current software (as implemented by [START_REF] Aguilar | Efficient computations of genomic relationship matrix and other matrices used in the single-step evaluation[END_REF] and hardware has allowed the use of single step with computed inverses of G and A 22 for up to 17,000 genotyped individuals and 18 traits [START_REF] Tsuruta | Multiple-trait genomic evaluation of linear type traits using genomic and phenotypic data in US Holsteins[END_REF]. This is possibly enough for many animal populations (e.g., in other dairy cattle or sheep populations). However, in the Holstein population of large countries (United States, Germany, and France, among others), the number of genotyped animals is already or will soon exceed 50,000 or 100,000 animals. Therefore, at some point, explicit inverses will no longer be possible to be computed and stored. Also, with high-density SNP chips, computation of G itself will become a challenge; again, using approaches that do not compute G explicitly, (e.g., solving [9-10] or [4] by BiCGstab) the cost is O(mp), that is, the same order of magnitude as reading the data. Further, iteration on data are available so memory is not a real challenge. Extensive testing in challenging situations (complex 9), G u ˆˆ, γ = 2 and A 22 φ = û 2 (Equation 10), respectively. 2 G = genomic relationship matrix; A 22 = pedigree relationship matrix. models and large number of genotypes) will be needed to find out the feasibility and accuracy of the solutions we propose in this work.

CONCLUSIONS

Single-step evaluations including records, pedigree, and genotypes can, in principle, be carried out with little additional computing costs or changes compared with 2-step procedures. Convergence of the method is guaranteed from existing theory. This will facilitate its introduction in routine evaluation.

( )

          - ′       exp 0.5 .
After some algebra, and using A A A A 0

( ) ′ - ( )           - - - u A A Zg A u A A Zg 1 1 2 22 1 1 1 1 1 2 22 1 σ u 2 -- ′       = - ′ ′ - 0.5 g D g u 1 exp . 0 5 1 g g A AZ Z A Z A A Z D ( ) ′ ′ - ( ) +     - - - - -- - 11 2 1 2 2 21 2 2 1 2 22 1 2 1 σ σ σ σ σ u u u u u                                         u g 1 .
This is the kernel of a normal distribution and, thus, p (u 1 , g) = N (0, T), where

T A AZ Z A Z A A Z D - - - - - -- - ′ ′ - ( ) +     1 11 2 1 2 2 21 2 2 2 2 22 1 2 1 = σ σ σ σ σ u u u u u               .
In fact,

T I 0 0 Z A A A A A A A A - - - - ′                  1 2 11 12 21 22 1 21 11 12 22 1 =σ u                             +                - I 0 0 Z 0 0 0 D 1 ,
Journal of dairy Science Vol. 95 no. 8, 2012

STraTegIeS For SIngLe-STeP genoMIc eVaLuaTIon 4641 which shows that T is always positive definite. Therefore, the MME can be constructed as

′ ′ ′ ′ ′ + ′ ′ ′ ′ ′ - X X X W X W W X W W A AZ Z W X ZA ZW 1 1 2 2 1 1 1 1 11 12 2 2 12 2 2 α α α σ u u u u W W Z Z A A Z D b u g 2 22 22 1 1 2 1 + ′ - ( ) +                        - - α σ u e ˆ        = ′ ′ ′ ′               X y W y Z W y 1 1 2 2 , [A1]
which, by including additional unknowns φ ( ) as in expression [4], can be written as

′ ′ ′ ′ ′ + ′ ′ ′ ′ ′ X X X W X W Z 0 W X W W A AZ 0 Z W X ZA ZW W 1 1 2 2 1 1 1 1 11 12 2 2 12 2 α α α u u u 2 2 22 1 2 22 1 Z ZA Z D Z 0 0 Z A b u + ′ + ′                     - α σ α α α u e u u u ˆĝ X y W y Z W y 0 --ϕ ϕ                   = ′ ′ ′ ′                  1 1 2 2  
.

[A2]

These equations are nevertheless more cumbersome than [4], [5], or [7].

1 1 b u u d } ˆ* *                                          1 2 2       = +               - Xy W y W y A A W y W y 1 1 1 2 2 22 21 1 1 2 2 ' ' ' ' { }            . [B2]
Now, mimicking the development of [START_REF] Quaas | Modified equations for sire models with groups[END_REF] for groups of unknown parents, let

S I 0 0 0 0 I 0 A A 0 0 I I 0 0 0 I 1 =                   - 12
1 - - = - -                   1 12 22 .
Inserting S -1 S between the coefficient matrix and the vector of unknowns and premultiplying both sides of the system by S -T , we get

                 =      -                                   ˆ* * b u u d 1 2 2 and ′ + - + - X X W W X W 0 W X W W A A A A A A W 22 1 1 1 2 2 1 1 1 1 11 12 11 12 12 2 ' ' ' ' ' ( ) α α α u u u X X A W W A AA A A 0 A A A A 22
- + + -                               ) ( ) } ˆ 1 1 2 2 1 1 2 2 ˆ' ' u d X y W y W y 0                   = ′                     
The rules for the inverse of a partition matrix [START_REF] Searle | Matrix Algebra Useful for Statistics[END_REF] -. Therefore, the system above simplifies to

′ + + - X X W W X W 0 W X W W A A 0 W X A WW A 1 1 2 2 1 1 1 1 11 12 2 2 21 2 2 22 ' ' ' ' ' ' α α α α u u u u α α α α u u u A 0 0 A A G A 1 1 1 1 22 22 22 22 - - - - - + - { }                     ( ) ˆ' ' b u u d X y W y W y 0 1 2 2 1 1 2 2                   = ′                      .
Absorbing the last equation of system, we get Applying the rules for the inverse of a partition matrix [START_REF] Searle | Matrix Algebra Useful for Statistics[END_REF] to matrix K leads to [C3]

K G A A A A G A A G

We have Now, this shows that if we have a model [C2] such that var(d) is not positive definite, there is an equivalent model [C3] in which all variance matrices are positive definite, which leads to the same MME [C4]. The only requisite is that the covariance matrix be nonsingular (i.e., invertible). In other words, there is no need to care whether G -A is positive definite or not.

E E E E E y Xb Z u ZV V e Z       = +       +             +       + + - - ' ' ' i i -       +                                    - + - - i u V D 0 0 D 1 α                                     = - + - + ˆ* ˆ' ' ' b u X y Z y V δ δ i iV Z y -                                     ' ' . Using d V V = -       + - i δ, V V V V I + - + -       -           = i i
STraTegIeS For SIngLe-STeP genoMIc eVaLuaTIon

2

  [10]. In other words, the solution of the equivalent systems [2],[4], and [7] can be achieved by iterating between system [9] and[10], where [9] can be viewed as the regular MME used for most classical genetic evaluations, except for a correction term γ of the right-hand side for genotyped animals. This contribution of genomic information is proportional to the difference between 2 similar expressions: ˆφ = - If A 22 and G are of large size, an iterative algorithm to solve A ˆγ = 2 can be implemented.

  Either solve ĝ = DZ′G -1 û 2 and ˆγ = - in Figure1. Only 2 individuals in this pedigree are genotyped: 4 and 6. We simulated by gene dropping 20 biallelic markers. As a result, matrices G and A 22 have values , and the heritability is h 2 = 0.1; the left-hand side of the mixed model Equation [4] is shown in Figure2to show the sparsity pattern. Equation [5] and regular SSGBLUP [3], unsymmetric extended SSGBLUP [4], iterative SSGBLUP [9-10] (with ω = 1) provided exactly the same results, namely û = (-0.12, 0.00, -0.17, 0.03, -0.11, -0.19) with ˆ. , . . 05, 0.05, -0.07, 0.10, -0.03, -0.03) when the genomic information is ignored.

  Matrix G was created and tuned to fit A 22 as in Vitezica et al. (2011) BLUP FST . Matrices G and A 22 were positive definite; G -A 22 was not: about half of its eigenvalues were negative. Three ways of estimating breeding values were used: regular SSGBLUP [3], unsymmetric extended SSGB-LUP [4], and iterative SSGBLUP [9-10]. In regular SSGBLUP, G and A 22 and their inverses were created using preGSf90 software (Aguilar et al., 2011); resolution was by preconditioned conjugate gradient (PCG) as implemented in the Blupf90 package. In unsymmetric extended SSGBLUP, G and A 22 were read from file and stored in memory; resolution was by BiCGstab as implemented by I. Aguilar (INIA Las Brujas, Canelones, Uruguay) in the Blupf90 package. In iterative SS-GBLUP, 2 solution loops were used: the outer loops updated the right-hand side in [9] and [10]; the inner loops solved [9] and [10]. Equation 9 was solved by PCG; as for [10], 2 strategies were used: one was to solve this mimics the iteration on data situation. The other one was to apply (17) using ˆφ = -

Figure 1 .

 1 Figure 1. Example pedigree.

Figure 2 .

 2 Figure 2. Left-hand side of the mixed model equations according to Equation 4 for the example.

Figure 3 .

 3 Figure 3. Convergence of regular single-step genomic BLUP (SSGBLUP), unsymmetric extended SSGBLUP, and iterative SSGBLUP in the simulated data. Color version available in the online PDF.

  the latter can be constructed by solving the systems of equations G u ˆγ = 2 and A the regular SSGBLUP Equation [3], the unsymmetric expanded SSGBLUP using [4] or [5], or the iterative SSGBLUP system [9-10] has a cost O(n) for the least squares, A -1 , and A 22 parts; and O(mp) for the G part (times a certain number of iterations to solve the equations G u ˆγ = 2 and A u 22 2

  u = u* + d = u* + V + δ + -iV -δ -.

Table 1 .

 1 Computing options and performances to solve the single-step model in a simulated example

	Item 1	Regular SSGBLUP 2	Unsymmetric extended SSGBLUP	Iterative SSGBLUP [1]	Iterative SSGBLUP [2]
	Needs inverted G -1 and A 22 1 -	Yes	No	No	Yes
	Solving method	PCG 3	BiCGstab 4	SUR 5 (outer loop)	SUR (outer loop)
	Number of iterations	62	240	PCG (Equations 9 and 10) 83 6	PCG (Equation 9) 83 6
	Computing time, s	2	47	286	19
	1 G = genomic relationship matrix; A 22 = pedigree relationship matrix.
	2 SSGBLUP = single-step genomic BLUP.

3 Preconditioned conjugate gradient. 4 Biconjugate gradient stabilized method. 5 Successive underrelaxation. 6 Outer loop.

Table 2 .

 2 Number of operations and memory needed for iteration on data with different solvers 1

			Unsymmetric
	Item 2	Regular SSGBLUP 3	extended SSGBLUP	Iterative SSGBLUP [1]	Iterative SSGBLUP [2]
			No. of operations
	G -1 and A 22 1 -needed Computing G -1 and A 22 1 -Solving (1 iteration)	Yes O(mp 2 ) + O(pn) + O(p 3 ) + O(p 3 ) O(p 2 + n)	No 0 O(pm + n + n)	No 0 O(t M n + t G pm + t A n)	Yes O(mp 2 ) + O(pn) + O(p 3 ) + O(p 3 ) O(t M n +2p 2 )
		O(p 2 + n)	Memory O(n + 2p)	O(n + 2p)	O(p 2 + n + 2p)

1 m = number of SNP; n = number of individuals in pedigree; p = number of genotyped individuals; t M , t G , and t A : number of iterations needed to achieve convergence of regular genetic evaluation (Equation

Journal of dairy Science Vol. 95 no. 8, 2012 STraTegIeS For SIngLe-STeP genoMIc eVaLuaTIon

Journal of dairy Science Vol. 95 no. 8, 2012

ACKNOWLEDGMENTS

The authors thank the French National Research Agency (ANR, Paris, France) project AMASGEN and APISGENE for financing. A. Legarra thanks the European FEDER funds for project GENOMIA through POCTEFA (http://www.poctefa.eu). Ignacy Misztal (University of Georgia, Athens), Ignacio Aguilar (INIA, Las Brujas, Uruguay), Gregor Gorjanc (University of Ljubljana, Ljubljana, Slovenia), Paul VanRaden (USDA, Beltsville, MD), Dave Johnson (LIC, Hamilton, New Zealand), and Daniel Gianola (University of Madison, WI) are thanked for fruitful comments. Use of code by Ignacio Aguilar for the BiCGstab solver and the unsymmetric storage for the Blupf90 software is particularly acknowledged. The authors acknowledge the editor and reviewers for suggestions and requirements largely improving the final manuscript. This project was partly supported by the Toulouse bioinformatics platform (France; http://bioinfo.genotoul.fr/).

APPENDIX A

A possibility in SSGBLUP MME is to specifically fit marker (or gametic) effects. Assume u 2 = Zg, where g are marker effects. Consider, as in [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF] ( )

The joint distribution of u 1 and g can be written as

APPENDIX B

Here, we provide details about the equivalent derivation of the MME. We start from the model

According to the assumptions in the main text, the corresponding MME are

APPENDIX C

In this appendix, we show how the problem of lack of positive definiteness of G -A can be circumvented. Consider the model

We want to consider an equivalent model where

where A is the numerator relationship matrix.

Matrices D + and D -are positive diagonal matrices of size n+ and n-. Decompose the matrix of eigenvectors accordingly:

.

We have

(note: d and δ + are vectors of real numbers and δ -is a vector of imaginary numbers). Assume δ δ

+ Therefore, the model with u* and δ is equivalent to the model with u* and d.

The model with u* and δ leads to usual MME, where the variance of δ is the positive definite matrix

Regular MME are usually derived for variables with real domain. However, in this case, we deal with δ -being imaginary random effects. Because δ + and δ -are uncorrelated and, further, they are either real or imaginary (but not complex numbers), they vary according to the so-called circular symmetric complex normal distribution, for which, roughly speaking, regular expressions of multivariate normality for real numbers apply [START_REF] Picinbono | Second-order complex random vectors and normal distributions[END_REF]. Indeed, some derivations of BLUP [START_REF] Henderson | Sire evaluation and genetic trends[END_REF] do not involve normality or even probability distribution functions, but only covariance matrices. These MME are