ELECTRICITE DE FRANCE L Viricel

GROUPE Des Laboratoires

France C Andrieu

ELECTRICITE DE FRANCE J C Segura

SEPTEN France A Rocher

ELECTRICITE D E France

J Thomazet

M H Clinard

Framatome Anp

COMMISSARIAT A L'ENERGIE France F Dacquait

Atomique

France

PRIMARY WATER CHEMISTRY OPTIMIZATION FOR EXTENDED FUEL CYCLE OPERATION : RESULTS OF THE « DUO EXPERIMENTATION » AFTER THREE CYCLES

The primary coolant conditioning in French nuclear power plants is essentially based on the boron-lithium coordinated chemistry, with a target pH of 7.2 at 300°C and a maximum lithium concentration of 2.2 mg/kg.

In 1996, EDF 1300 MWe units began operating 18-month fuel cycles, increasing boron concentrations at the beginning of the cycles. Since today the maximum lithium concentration in normal operation is 2.2 mg/kg, extended cycle operation results in a decrease in the pH at the beginning of the cycles, which may possibly lead to deposits in RCS, and particularly on the fuel cladding, and increased dose rates. It has to be noted that today, the fuel assemblies maximum burnup is set at 52 GWd/tU.

One solution is to adjust the pH by increasing the lithium content at the beginning of the cycles, which is easy to implement and does not require any modification on the units. Hence, EDF is testing a "modified" chemistry regime in the « DUO experimentation » during 4 fuel cycles, with a maximum authorized lithium content of 3.5 mg/kg at the beginning of the cycles in the Cattenom 2 pilot unit. The Golfech 1 reference unit implements a standard boron-lithium coordination pH 300 7.2. The major goal of the experimentation is to assess the impact of elevated lithium concentrations at the beginning of the cycles on fuel cladding oxide behavior, mass transport and dose rates. This paper presents the results of the first three cycles of the DUO experimentation.

Oxide thickness measurements results show that operating with EDF modified chemistry, hence lithium contents up to 3.5 mg/kg at the beginning of the cycles, is encouraging. No influence of elevated lithium concentrations was observed on oxide thickness up to a burnup of about 50 GWd/tU, which is slightly lower than the maximum burnup a fuel assembly is allowed to achieve ; moreover, advanced fuel cladding alloys show a good behavior.

Concerning contamination and dose rates, the 58 Co surface activities measured in Golfech 1 primary system decreased from the last 12-month cycle to the first 18-month cycle, whereas the 60 Co surface activity increased after the first two conversion cycles, probably due to an abnormal cobalt input during a conversion cycle. In Cattenom 2, the 58 Co and 60 Co surface activities in the steam generators increased from the last 12-month cycle to the third conversion cycle, whereas they remain overall stable or even decrease in the crossover legs and in the hot legs. Overall, the DUO experimentation results are comparable to French plants usual behavior.

While any beneficial impact of a constant pH of 7.2 at 300°C throughout the cycle on dose rates is not proved given the French plants results, the easy implementation of increased lithium concentrations at the beginning of the cycles may help face issues due to new French fuel management strategies, leading to fuel assemblies maximum burnups of 52 GWd/tU.

As the DUO fuel assemblies were re-loaded for a fourth cycle in Cattenom 2, burnups of about 55 to 60 GWd/tU should be achieved. The same surveillance program as during the three previous cycles is applied ; hence, it should provide extensive data, especially on the fuel cladding behavior up to 60 GWd/tU and on contamination when a modified chemistry regime is implemented.

Today, given the pH improvement obtained by operating EDF modified chemistry with increased lithium concentrations up to 3.5 mg/kg at the beginning of the fuel cycles, generalizing this chemistry regime to the French plants is under consideration.

INTRODUCTION

The primary coolant conditioning in French nuclear power plants is essentially based on the boron-lithium coordinated chemistry, with a pH target value of 7.2 at 300°C and a maximum lithium concentration of 2.2 mg/kg [START_REF] Nordmann | Primary and secondary water chemistry of French PWR units, Fourth International Seminar on Primary and Secondary Side Chemistry of NPP[END_REF]. Some years ago, the target pH value was changed from 7.0 to 7.2 to decrease dose rates [START_REF] Anthoni | Effects of pH of primary coolant on PWR contamination, Water Chemistry of nuclear reactor system[END_REF].

In 1996, EDF 1300 MWe units began operating 18-month extended fuel cycles, increasing fuel uranium oxide enrichment from 3.2 to 4 %. This change induced :

a significant increase in the boron concentration at the beginning of the fuel cycles (progressively for the first two conversion cycles) ; a decrease of the pH at the beginning of the cycles, since today the maximum lithium concentration in normal operation is 2.2 mg/kg. However, the twenty 1300 MWe units operating extended cycles are not the only EDF units concerned by higher boron concentrations at the beginning of the cycles. Other French nuclear reactors can be impacted if new core designs are implemented (Mixed OXide, ...).

It is well known that operating at pH 300 below 6.9 for a significant period of time makes it difficult to keep corrosion products as low as possible, possibly leading to :

deposits in RCS, and particularly on the fuel cladding, which may participate in axial offset anomalies occurrence [START_REF] Robin | Modified lithium chemistry for extended fuel cycles[END_REF]- [START_REF]PWR primary water chemistry[END_REF] ; increase in dose rates. One solution is to adjust the pH by increasing the lithium concentration at the beginning of the cycles. This solution is easy to implement and does not require any modification on the units. A new maximum lithium concentration, consistent with materials behavior and especially with fuel cladding, has to be determined.

EDF has decided to test a "modified" chemistry regime in the « DUO experimentation » during 4 fuel cycles, one major goal being assessing of the corrosion behavior of several fuel cladding materials. This paper deals with the results of the controls and expertises performed during the first three cycles of the DUO experimentation.

DUO EXPERIMENTATION FEATURES

This experimentation deals with the specific surveillance, from the first conversion cycle to the second 18-month cycle, of two 1300 MWe units operating extended fuel cycles. The reference unit and the pilot unit for this experimentation were chosen to be Golfech 1 and Cattenom 2 respectively. They had to have similar design (steam generator tubes bundle in 600 TT alloy, fuel grids assemblies, no hot spots, no flux anomalies, ...).

Golfech 1 implements a classic boron-lithium coordination pH 300 7.2 (Figure 1). Cattenom 2 operates EDF modified chemistry at the beginning of the cycles, with a maximum authorized lithium concentration of 3.5 mg/kg (Figure 2).

The major goal of the experimentation is to assess the impact of elevated lithium concentrations at the beginning of the cycles on fuel cladding oxide behavior, mass transport and dose rates.

SURVEILLANCE PROGRAM OF THE DUO EXPERIMENTATION

A specific surveillance program has been implemented during the four cycles of the DUO experimentation (Table 1). During each cycle, the primary coolant surveillance deals with its chemical and radiochemical features. During shutdown transients and outages, it concerns primary coolant activity, surface activity and fuel oxide thickness measurements.

These controls aim at assessing the performance of an extended cycles modified chemistry when elevated initial boron concentrations are used compared to a classic chemistry, which may lead to operate, at the beginning of the cycles, at a less favorable pH for fuel crud deposits and dose rates.

Field

DUO EXPERIMENTATION : RESULTS AFTER THREE CYCLES

Chemistry

During the first three cycles of the DUO experimentation in Golfech 1 and Cattenom 2, despite load following, the boron-lithium coordinated chemistry regime matched the specifications (Figures 3 and4). Data out of the normal operating area are essentially due to shutdowns, physical testing, rapid load variations and makeups for tritium dilution in RCS. No specific pollution was observed during these cycles in Golfech 1 nor in Cattenom 2.

(2)

(2) (3) (3)
(2)

(3) (3) (4) (3) (4) (3) (4)
(3) ( 2)

24
(4) (5) (2) (3) (3)
(5)

The initial boron concentrations at hot zero power and full power are reported in table 2. The Cattenom 2 pilot unit maintained a lithium concentration between 3.5 and 2.2 mg/kg for about three months during its third conversion cycle ; operation with a lithium content between 3.5 and 2.5 mg/kg lasted about one month (Table 3).

GOLFECH

Cycle #

Operating time with lithium between 3.5 and 2. In the Golfech 1 reference unit, pH was around 6.8 at 300°C at the beginning of the first and second conversion cycles, and 2 to 3 days were needed to reach 6.9. At the beginning of the first 18-month cycle (third conversion cycle), pH was around 6.7 and about 10 days were needed to reach 6.9 ; overall, it took about 42 weeks to reach a pH of 7.2 at 300°C. In Cattenom 2, pH at the beginning of the cycles was always above 6.9 at 300°C, and 32 weeks were needed during the third conversion cycle to reach a pH of 7.2 at 300°C. The modified chemistry in Cattenom 2 allows a better control of the pH just before criticality than does the standard chemistry regime in Golfech 1. Feed-and-bleeds due to load following may result in primary water chemistry transient changes because of system equilibrium disrupting, dilution, oxygen ingress, etc … Particularly, water and boron makeups may generate oxygen ingress in RCS which may impact the coolant-materials interactions. during power increase at the beginning of the DUO experimentation cycles or during power decrease before refueling shutdown. These results agree with the values typically observed in French plants. In Golfech 1, the concentration in soluble nickel is higher at the beginning of the second conversion cycle and of the first real 18-month cycle than it was for the first conversion cycle. For both units, no significant differences can be observed at the end of the cycles, from the second to the third conversion cycle, except for a slight increase in the chromium content in Golfech 1 and a significant decrease at 0 % power in Cattenom 2. Cr in Golfech and 25 MBq/t in Cattenom. Some higher activities could sometimes be observed, most of them related to significant makeups or load reductions to 0 % power.

Unit

The activities measured during shutdowns are presented in table 6. In Golfech 1, the cobalt 58 peak activity was similar for the second conversion cycle and for the first real 18-month cycle ; it was somewhat smaller than previously observed. In Cattenom 2, an increase was noted at the end of the first conversion cycle compared to the last 12-month cycle ; then, it decreased continuously until the third conversion cycle. In Cattenom 2, the total released activity increased during the first two conversion cycles, but at the end of the third one, it decreased and reached the initial level before extended cycle operation. So today, no significant effect of the cycles length is observed on the released activity. Moreover, the two units behave similarly, despite the differences in the chemistry regime implemented for the last three cycles : they have shown an overall decreasing trend in the 58 Co peak activity for several years.

A B C D E F G Q H I J K L M N O P A B C D E F G Q H I J K L M

Figure 5 : Characterized assemblies location map

The zirconia layer thickness measurements were performed using the blade device developed by FRAMATOME ANP, based on Eddy Current measurements [START_REF] Thomazet | In-reactor fuel cladding external corrosion measurement process and results[END_REF]. During each refueling outage of the DUO experimentation, the device is installed on a criss-cross measurement table, allowing blade insertion into each row of rods at a chosen measurement altitude. All analyzed rods, particularly those with claddings made of advanced alloys, are therefore measured during blade insertion and extraction from all four assembly faces leading to more than 2000 measurements collected per campaign. Statistical process using the « Lagos » and « Zircone » software programs is used for this analysis. The measurement accuracy is estimated to be ± 5 µm.

During the first three cycles of the DUO experimentation, neither in core axial offset anomalies nor abnormal cruds on assemblies had to be reported. Hence, only oxide thickness measurements were carried out on fuel rods.

Concerning Zircaloy-4, the maximum oxide thicknesses measured at Golfech 1 and Cattenom 2 at the end of each cycle are plotted on figure 6. Up to a burnup of about 45 GWd/tU, no difference was observed between the maximum oxide thickness measured in Golfech 1 and in Cattenom 2. Since higher burnups were achieved in Cattenom 2 than in Golfech 1 on the characterized assemblies, it was not possible to compared directly the behavior of the Zircaloy-4 fuel cladding in these two units. However, burnups similar to Cattenom 2 were achieved in the Paluel 4 reactor, were the same chemistry regime as the one in Golfech was implemented.

Comparing the maximum fuel oxide thickness measured in Cattenom 2 and Paluel 4, no significant differences could be observed for burnups up to about 50 GWd/tU.

Moreover, the results obtained on the advanced alloys M4 and M5 in Cattenom 2 and Golfech 1 are within the scattering observed on these alloys in French reactors (900 and 1300 MWe) implementing a 2.2 mg/kg maximum lithium chemistry regime (figures 7 and 8 for M5 and M4 alloys respectively).

The modified lithium chemistry does not seem to have an obvious impact on the corrosion behavior of fuel cladding after three fuel cycles, and for average fuel rod burnups up to about 50 GWd/tU.

Legend

DISCUSSION

The results of the DUO experimentation surveillance program after three fuel cycles show a good behavior of the fuel cladding in Cattenom 2, the pilot unit implementing EDF modified chemistry with a maximum authorized lithium concentration of 3.5 mg/kg at the beginning of the extended cycles.

Concerning contamination and dose rates, the 58 Co surface activities measured in Golfech 1 primary system decreased from the last 12-month cycle to the first 18-month cycle, whereas the 60 Co surface activity increased after the first two conversion cycles, probably due to an abnormal cobalt input during a conversion cycle. In Cattenom 2, the 58 Co and 60 Co surface activities in the steam generators increased from the last 12-month cycle to the third conversion cycle, whereas they remain overall stable or even decrease in the crossover legs and in the hot legs. Overall, the DUO experimentation results are comparable to French plants usual behavior.

Moreover, non destructive steam generator tubing examinations proved that no abnormal corrosion due to elevated lithium concentrations occurred on 600 TT alloy, since the beginning of extended cycle operation in Cattenom 2.

From these results, it appears possible to improve the pH at the beginning of extended fuel cycles, beneficial to a better behavior of the RCS. It can be noticed that boron concentrations up to 1800 mg/kg at the beginning of fuel cycles can be managed with EDF modified chemistry, ensuring pH 300 ≥ 6.9.

Studying operating with a constant lithium concentration of 3.5 mg/kg until pH 300 7.2 is reached could also be of interest in order to reach the target pH 300 7.2 as soon as possible, but has not been tested in this DUO Experimentation.

Figure 1 :Figure 2 :

 12 Figure 1 : Classic B-Li coordinated chemistry pH 300 7.2 applied in Golfech 1

Figure 3 :Figure 4 :

 34 Figure 3 : Golfech 1, cycle 8, 3 rd cycle of the DUO Experimentation : B-Li diagram

  : Hot shutdown (2) : Hot standby (4) : Tritium dilution (3) : Physical testing (5) : Load regulation

Figure 6 Figure 8 :

 68 Figure 6 : Zircaloy-4 external corrosion experience

Table 1 : Surveyed parameters

 1

	Surveillance

Table 2 : Boron at hot zero power and full power

 2

			1			CATTENOM 2
		Cycle	B at hot zero	B at full power	Cycle	B at hot zero	B at full power
		#	power (mg/kg)	(mg/kg)	#	power (mg/kg)	(mg/kg)
	Last 12-month cycle	5		970	7		1020
	1 st conversion cycle	6	~ 1450	1020	8	~ 1480	1020
	2 nd conversion cycle	7	~ 1450	1010	9	~ 1580	1100
	3 rd conversion cycle	8	~ 1630	1220	10	~ 1650	1163

Table 3 : Operating time with lithium concentrations greater than 2.2 mg/kg in Cattenom 2

 3

		Operating time with lithium
	5 mg/kg	between 2.5 and 2.2 mg/kg

Table 4 : Golfech 1 and Cattenom 2 makeups to RCS In

 4 Table 4 summarizes the makeup volumes in RCS in Golfech 1 and Cattenom 2. The greater makeup volumes in Golfech 1 compared to Cattenom 2 can be explained by more load following in Golfech. Similarly, the variations in the volumes used during the different cycles on both units is explained by different load variations. table 5 are reported the results of corrosion products measurements performed in Golfech 1 and Cattenom 2

			Makeup water	Boron makeup water
		Cycle #	Volume (m 3 )	Mean dissolved O 2 (µg/kg)	Volume (m 3 )	Mean dissolved O 2 (µg/kg)
				Golfech 1
	1 st conversion cycle	6	12304	37	827	1700
	2 nd conversion cycle	7	9242	49	714	1777
	3 rd conversion cycle	8	14398	57	1381	1923
			Cattenom 2
	1 st conversion cycle	8	9600	45	467	2200
	2 nd conversion cycle	9	8528	51	444	1523
	3 rd conversion cycle	10	9963	54	694	2502

Table 5 : Soluble corrosion products in Golfech 1 and Cattenom 2

 5

		Cycle	Cycle # % Power Ni (ng/kg) Co (ng/kg) Cr (ng/kg)
		Power increase at beginning of cycles
		1 st conversion cycle	6	100	350	6	1 600
	GOL 1	2 nd conversion cycle	7	100	3 170	10	1 190
		3 rd conversion cycle	8	100	2 300	< 20	500
		1 st conversion cycle	8	8 100	11 000 260	< 4 < 4	60 230
	CAT 2	2 nd conversion cycle	9	8 100	1 355 1 012	22 < 4	117 NA
		3 rd conversion cycle	10	8 100	4 780 660	20 6	180 200
		Power decrease at end of cycles
				90	392	< 10	166
		2 nd conversion cycle	7	70 0 (290°C)	215 20 000	< 4 NA	190 800
	GOL 1			0 (175°C) 80	862 000 2 630	NA < 100	900 1 290
		3 rd conversion cycle	8	0 (296°C) 0 (175°C)	16 400 730 000	< 100 4 000	1 320 1 100
				0 (120°C)	960 000	6 500	1 400
				75	890	< 4	690
		2 nd conversion cycle	9	0 (286°C) 0 (175°C)	610 > 200 000	NA NA	1 150 3 210
	CAT 2			0 (120°C) 80	1 584 000 440	NA < 10	6 880 1 070
		3 rd conversion cycle	10	0 (290°C) 0 (175°C)	370 1 060 000	< 10 720	320 970
				0 (120°C)	1 080 000	870	460
	NA = Not Available

Table 6 : Shutdown activities

 6

			Golfech 1		Cattenom 2
		Cycle	58 Co peak	Total released	Cycle	58 Co peak	Total released
		#	activity (GBq/t)	activity (GBq)	#	activity (GBq/t)	activity (GBq)
	Last 12-month cycle	5	160		100 000	7	165	100 000
	1 st conversion cycle	6	170			60 000	8	207	134 000
	2 nd conversion cycle	7	105			70 000	9	197	135 000
	3 rd conversion cycle	8	100 + 33	(1)			10	180	100 000
	(1) : A cold shutdown with oxygenation occurred three months before refuelling shutdown ; the	58 Co peak
	activity was 100 GBq/t for this shutdown and 33 GBq/t for the refuelling shutdown
	Deposits
	Deposited activities were measured on steam generators (SG) tube bundle, hot legs (HL) and crossover
	legs (CL). Cobalt 58 and cobalt 60 representing more than 80 % of the dose rates around the primary system,
	only their activities are reported in tables 7 and 8 for Golfech 1 and Cattenom 2 respectively. The accuracy is about 25 % for 58 Co measurements and 15 % for 60 Co measurements. Moreover, dose rates were measured
	around steam generators and legs (table 9).
						SG (GBq/m²) HL (GBq/m²) CL (GBq/m²)
			Cycle #	58 Co	60 Co	58 Co	60 Co	58 Co	60 Co
	Last 12-month cycle	5		2.41	0.30 1.81 0.71 4.31	0.64
	1 st conversion cycle	6		2.36	0.52 1.49 0.89 3.69	0.89
	2 nd conversion cycle	7		2.42	0.61 1.21 2.03 3.53	1.02
	3 rd conversion cycle	8		1.72	0.60 1.07 1.46 3.12	0.99

Table 7 : Surface activity in Golfech 1

 7

			SG (GBq/m²) HL (GBq/m²) CL (GBq/m²)
		Cycle #	58 Co	60 Co	58 Co	60 Co	58 Co	60 Co
	Last 12-month cycle	7	2.06 0.53	2.35	1.49 1.72 0.73
	1 st conversion cycle	8	2.23 0.57	1.50	1.62 1.52 0.68
	2 nd conversion cycle	9	3.20 0.59	1.41	1.04 1.77 0.63
	3 rd conversion cycle	10	2.64 0.61	1.06	1.50 1.46 0.57

Table 8 : Surface activity in Cattenom 2

 8

Co surface activity in the steam generators was constant from the last 12-month cycles to the second conversion cycle, and then decreased at the end of the first 18-month cycle ; the 58 Co surface activity in the hot legs and crossover legs decreased constantly from the last 12-month cycle to the first 18-month cycle. The 60 Co surface activity results increased from the last 12-month cycle to the second conversion cycle, probably because of an abnormal input of cobalt (i.e. stellites) [START_REF] Bretelle | Extended fuel cycle : results of the first two cycles of the « DUO Experimentation[END_REF] ; then, it stabilized at the end of the first 18-month cycle.

In Cattenom 2, the 58 Co surface activity in the steam generators increased significantly from the first to the second conversion cycle, then decreased at the end of the third conversion cycle. From the last 12-month cycle to the third conversion cycle, the 58 Co surface activity decreased constantly in the hot legs and remained stable in the crossover legs. The 60 Co surface activity increased slightly from the last 12-month cycle to the third conversion cycle in the steam generators and decreased slightly in the crossover legs ; it remained overall stable in the hot legs when comparing the last 12-month cycle to the third conversion cycle. Concerning dose rates, the values of legs are impacted by high surrounding dose rates : 50 % to 70 % of the dose rates in the primary coolant pipes come from the auxiliary system.

Hence, compared to the previous cycles, the results of the third conversion cycle show that :

in Golfech 1, the 58 Co surface activity in the primary system decreased from the last 12-month cycle to the first 18-month cycle, whereas the 60 Co surface activity increased after the first two conversion cycles, probably due to an abnormal cobalt input ; in Cattenom 2, the 58 Co and 60 Co surface activities in the steam generators increased from the last 12-month cycle to the third conversion cycle, whereas they remain overall stable or even decrease in the crossover legs and in the hot legs. Overall, the cobalt 58 and cobalt 60 surface activities measured in the Golfech 1 and Cattenom 2 primary systems are comparable to French plants usual behavior.

Fuel assemblies

In order to assess the impact of the reactor coolant chemistry on the fuel cladding corrosion, oxide thickness measurements were planned on characterized fuel assemblies loaded in the Golfech 1 and Cattenom 2 reactors.

Four characterized fuel assemblies were loaded in Golfech 1 and Cattenom 2. The fuel rods claddings are made of Zircaloy-4, alloy M4 or alloy M5. Each alloy comes from the same tube batch. The chemical composition of the alloys is reported in table 10.

Alloy

Sn (%) Nb (%) Fe (%) Cr (%) V (%) O (%)

Zircaloy-4 1.3 0.2 0.1 0.12 M4 0.5 0.6 0.4 0.12 M5 1 0.12

Table 10 : Chemical composition of fuel rods cladding alloys

The fuel rods were loaded into eight assemblies so as to distribute uniformly all the rods cladding batches into each assembly. Their locations in the assemblies are shown in figure 5. Concerning the advanced alloys, the location of the corresponding rods were chosen in order to be able to measure zirconia layer thickness along four axes 90° apart. However, some of them were also loaded in peripheral location in order to run in-depth visual examinations.

CONCLUSION

Oxide thickness measurements results show that operating with EDF modified chemistry, hence lithium concentrations up to 3.5 mg/kg at the beginning of the cycles, is encouraging. No influence of elevated lithium concentrations was observed on oxide thickness up to a burnup of about 50 GWd/tU, which is slightly lower than the maximum burnup a fuel assembly is allowed to achieve ; moreover, advanced fuel cladding alloys show a good behavior.

Concerning contamination and dose rates, the 58 Co surface activities measured in Golfech 1 primary system decreased from the last 12-month cycle to the first 18-month cycle, whereas the 60 Co surface activity increased after the first two conversion cycles, probably due to an abnormal cobalt input during a conversion cycle. In Cattenom 2, the 58 Co and 60 Co surface activities in the steam generators increased from the last 12-month cycle to the third conversion cycle, whereas they remain overall stable or even decrease in the crossover legs and in the hot legs. Overall, the DUO experimentation results are comparable to French plants usual behavior.

While any beneficial impact of a constant pH of 7.2 at 300°C throughout the cycle on dose rates is not proved given the French plants results, the easy implementation of increased lithium concentrations at the beginning of the cycles may help face issues due to French 18-month extended cycles management, leading to burnups of about 50 GWd/tU.

As going to very high burnups is under investigation, the DUO experimentation has been extended to a fourth cycle. The characterized fuel assemblies were re-loaded for this fourth cycle in Cattenom 2, and burnups of about 55 to 60 GWd/tU should be achieved at the end of this cycle. The same surveillance program as during the three previous cycles is applied ; hence, it should provide extensive data, especially on the fuel cladding behavior up to 60 GWd/tU and on contamination when a modified chemistry regime is implemented.

However, new French fuel management strategies could lead EDF to examine again the possible use of enriched boron 10. According to studies carried out in the Robust Fuel Program, the impact of boron on flux anomalies will also have to be taken into account when defining EDF strategy on using enriched boron 10. Today, given the pH improvement obtained by operating EDF modified chemistry with increased lithium concentrations up to 3.5 mg/kg at the beginning of the fuel cycles, generalizing this chemistry regime to the French plants is under consideration.


G Gaeta
email: giuseppe.gaeta@unimi.it

On some dynamical features of the complete Moran model for neutral evolution in the presence of mutations

We present a version of the classical Moran model, in which mutations are taken into account; the possibility of mutations was introduced by Moran in his seminal paper, but it is more often overlooked in discussing the Moran model. For this model, fixation is prevented by mutation, and we have an ergodic Markov process; the equilibrium distribution for such a process was determined by Moran. The problems we consider in this paper are those of first hitting either one of the "pure" (uniform population) states, depending on the initial state; and that of first hitting times. The presence of mutations leads to a nonlinear dependence of the hitting probabilities on the initial state, and to a larger mean hitting time compared to the mutation-free process (in which case hitting corresponds to fixation of one of the alleles).

Introduction

The Moran model [START_REF] Moran | Random processes in genetics[END_REF][START_REF] Karlin | On a genetics model of Moran[END_REF][START_REF] Donnelly | The transient behaviour of the Moran model in population genetics[END_REF] is a classical tool in studying neutral evolution [START_REF] Kimura | Diffusion Models in Population Genetics[END_REF][START_REF] Kimura | Evolutionary rate at the molecular level[END_REF][START_REF] Kimura | DNA and the neutral theory[END_REF][START_REF] Kimura | Molecular evolutionary clock and the neutral theory[END_REF][START_REF] Kimura | The neutral theory of molecular evolution: A review of recent evidence[END_REF][START_REF] Kimura | Population genetics, molecular evolution, and the neutral theory: selected papers[END_REF][START_REF] Ohta | The nearly neutral theory of molecular evolution[END_REF][START_REF] Ohta | Development of neutral and nearly neutral theories[END_REF][START_REF] Leigh | Neutral theory: a historical perspective[END_REF], i.e. the evolution of genome when several alleles are present but none of these presents any advantage over the others; thus in this framework evolution -and in particular fixation of a single allele, or at intermediate times elimination of some of the different alleles -is due uniquely to random fluctuations in the reproduction rates of individuals carrying the different alleles. Needless to say, the relevance of these fluctuations is intimately related to the fact any population is of finite size.

The Moran model deals with the simple case in which only two alleles are present at a locus, and the population size is fixed; this captures the essen-1 tial features of neutral evolution yet is simple enough to be amenable to exact solution.

On the other hand, the standard Moran model only considers reproduction fluctuations, and does not take into account the presence of mutations.

It should be mentioned in this regard that the original paper by Moran [START_REF] Moran | Random processes in genetics[END_REF] also took mutations into account and provided results for this case, see below for detail. Most of the literature dealing with the Moran mode (including that considering generalizations of the Moran model, see e.g. [START_REF] Baake | Ancestral processes with selection: Branching and Moran models[END_REF][START_REF] Chalub | On the stochastic evolution of finite populations[END_REF][START_REF] Chalub | Fitness potentials and qualitative properties of the Wright-Fisher dynamics[END_REF][START_REF] Muirhead | Modeling Multiallelic Selection Using a Moran Model[END_REF][START_REF] Nowak | Evolutionary dynamics in structured populations[END_REF][START_REF] Traulsen | Coevolutionary Dynamics: From Finite to Infinite Populations[END_REF][START_REF] Traulsen | Stochastic evolutionary dynamics on two levels[END_REF][START_REF] Traulsen | Stochastic dynamics of invasion and fixation[END_REF][START_REF] Vogl | A nearly-neutral biallelic Moran model with biased mutation and linear and quadratic selection[END_REF][START_REF] Dalmau | Convergence of a Moran model to Eigen's quasispecies model[END_REF][START_REF] Cerf | The distribution of the quasispecies for a Moran model on the sharp peak landscape[END_REF][START_REF] Cerf | Critical population and error threshold on the sharp peak landscape for a Moran model[END_REF][START_REF] Cerf | Critical population and error threshold on the sharp peak landscape for the Wright-Fisher model[END_REF][START_REF] Dalmau | The Wright-Fisher model for class-dependent fitness landscapes[END_REF][START_REF] Cerf | Quasispecies on class-dependent fitness landscapes[END_REF][START_REF] Cerf | The quasispecies equation and classical population models[END_REF][START_REF] Saakian | Finite population size effects in quasispecies models with single-peak fitness landscape[END_REF][START_REF] Saakian | Allele fixation probability in a Moran model with fluctuating fitness landscapes[END_REF]), however, disregards the possibility of mutations; we will hence refer to the Moran model without mutations as the standard Moran model.

Taking into account mutations does of course raise immediately a substantial problem. That is, if mutations are possible, fixation becomes impossible. The question to be asked concerns then on the one hand the equilibrium distribution of the process (this was determined by Moran [START_REF] Moran | Random processes in genetics[END_REF]) and on the other the dynamics of the model, in particular, given an initial state, the probability first hitting either one of the extreme states, i.e. those corresponding to a uniform population and which will therefore also called pure states (all other states will be called mixed states); and the hitting times for these. That is, we are led to study a problem of stopping times, or of exit from a domain; in this case, the domain is that of "mixed states", and exit from this domain corresponds to hitting of the "pure states", representing the border of the set of allowed states.

In the present paper, we consider a version of the Moran model, in which mutations are taken into account, and evolution -in the same sense as in the standard Moran model -is the result of fluctuations in the mutation rates as well as in reproductive rates for the two alleles; this will hence be dubbed as the "complete" Moran model (as opposed to the "standard" Moran model, see above). As already mentioned, this is not a new model, but corresponds to the full model considered by Moran in 1957. We will then consider the two problems mentioned above (first hitting either one of the pure states and hitting times), see below for more detail. 1Remark 1. When mutation is present, the Moran process is an ergodic Markov chain, hence it has a unique equilibrium distribution; this is attributed by Moran to Wright (see eq.( 9) in Moran's paper [START_REF] Moran | Random processes in genetics[END_REF]), and is given by

P (x) = Γ(β 1 + β 2 ) Γ(β 1 ) Γ(β 2 ) (1 -x) β1-1 x β2-1 , (1)
where

x = k/n , β i = n α i ,
with n the population size and α i the mutation rates for the allele i = 1, 2, i.e. the probability that a descendant of an individual with allele i has the different allele. ⊙

When discussing neutral evolution one is of course interested in the case where mutations do not favor any of the alleles, hence -in the framework where only two alleles are present, as in the Moran model -they are symmetric in the two alleles. We will thus consider, in the notation of Remark 1, the case α 1 = α 2 = α, which also entails

β 1 = β 2 = β.
The complete Moran model we consider is amenable to a limited random walk, like many biological models [START_REF] Berg | Random Walks in Biology[END_REF][START_REF] Codling | Random walk models in Biology[END_REF]; at difference with the standard Moran model, however, it does not correspond to a symmetric random walk (even for symmetric mutation rates).

We will state more precisely the questions we will investigate:

1. How the introduction of mutations affects the first hitting probabilities (that is, the probability that one "pure" state, in which all the population shares the same allele, is reached before the other one) for the two alleles, and more precisely their dependence on the initial state k, 0 < k < n;

2. How the introduction of mutations affects the mean hitting time (this is the time at which either one of the "pure" states is first reached) for a process starting from the initial state k, 0 < k < n.

It turns out that for both questions one could obtain analytical formulas in terms of hypergeometric functions; these explicit formulas are however rather involved and of little practical use (they are given in the Appendix). On the other hand, the two questions mentioned above may also be discussed in terms of exact (analytical, not numerical) computations for arbitrary but given population number n; these will allow to get to the point with minimal mathematical complications.

Remark 2. Needless to say, the "complete" model considered in this note (like the standard Moran model usually considered in the literature) is far too simple to be directly applicable to concrete current biological problem; it should instead be seen as a way to get a qualitative -and partially quantitative -understanding of the balance among reproductive versus mutational fluctuations, and their role in neutral evolution. ⊙ Remark 3. Finally, we note that the fixation problem can of course also be considered when one of the two alleles is advantageous, in particular in the case of beneficial mutations; this framework is outside of the scope of the present work, and is discussed e.g. in the review paper by Patwa and Wahl [START_REF] Patwa | The fixation probability of beneficial mutations[END_REF]. ⊙

The standard Moran model

In the standard Moran model [START_REF] Moran | Random processes in genetics[END_REF][START_REF] Karlin | On a genetics model of Moran[END_REF][START_REF] Donnelly | The transient behaviour of the Moran model in population genetics[END_REF], we have a haploid population of n individuals, each of them carrying an allele A or B of the only gene under study. At each time step, two individuals from the population are chosen (they may coincide). One is allowed to reproduce, and the other is eliminated; the one which is reproducing transmits the allele A or B to its heir. Denoting by n A = k, n B = n -k the number of individuals carrying the alleles A and B respectively this dynamics means that at each step, depending on which alleles are carried by the individuals chosen for reproduction and for elimination, the number k can either remain unchanged, k = k, or change to k = k ± 1. Denoting by p k the probability for k → k = k + 1 and by q k the probability for k → k = k -1, it is easily seen that

p k = q k = k n 1 - k n . (2)
Note that for k = 0 or k = n, there will be no further change in the composition of the population: all its members carry the same allele (B for k = 0, A for k = n), and this is definitely fixed.

The fixation time T (k, n), i.e. the mean time to fixation for a process starting from state k in a population of n individuals, is computed to be

T (k, n) = n k m=1 n -k n -m + n-1 m=k+1 k m . (3)
For large n we have

T (k, n) ≈ T (k, n),
where

T (k, n) := -n 2 1 - k n log 1 - k n + k n log k n . (4)
3 The Moran model in the presence of mutations

The standard Moran model aims at illustrating how fixation emerges as a result of fluctuations in the reproductive rate of the two alleles; on the other hand, in general there will also be mutations making that descendants of individuals with one allele will actually have the other allele. Even in the case where the mutation rates are the same for the A → B and the B → A mutations, there will be fluctuations in mutations rates -exactly as in reproduction rates -and these can also lead to reaching transiently a "pure" state (fixation is not possible, since mutations can halter a population sharing the same allele). More precisely, in the neutral case of interest here and in the presence of mutations, reaching a pure state will result from a combination of reproduction and mutation fluctuations. Thus we would like to understand how the combination of the reproductive and mutation fluctuations affects the dynamics, and in particular the hitting of the extreme, pure states, in neutral evolution; in particular, how the hitting probabilities for the two pure states and the mean hitting time respond to the presence of mutations.

We will thus consider a simple Moran-like model in which mutations are also present; this will be dubbed the "complete" Moran model (as mentioned in the Introduction, this was actually considered by Moran in his seminal 1957 paper).

The process undergone by this model will be described as follows. We have a homozygous population of n individuals, and the gene looked upon can exist in two alleles, A and B; the state of the system is labeled by the number k of individuals with the A allele (correspondingly, there will be n -k with the B allele).

At each step, several operations are performed, i.e.:

(i) A single individual is chosen, with equal probability for each of the n individuals;

(ii) The chosen individuals reproduces, giving birth to a single descendant; this has the same allele of the parent with probability (1 -µ), and the other allele with probability µ ;

(iii) Then again a single individual is chosen from the pre-existing population, again with uniform probability, and this is eliminated (the newborn is excluded from this choice, but the chosen parent can coincide with the eliminated individual).

We stress that the process is quite obviously a Markov chain, and actually an ergodic one. Note that for µ = 0 we recover exactly the standard Moran model (which is not ergodic).

It is clear that the population size remain constant, by construction; moreover at each step the state k can only be changed to k = k ± 1, or remain equal (unless k = 0, in which case it cannot pass to k -1; or k = n, in which case it cannot pass to k + 1). We thus have to evaluate the probabilities for these three events. Actually we are only interested in evaluating such probabilities for k ̸ = 0 and k ̸ = n: in fact once either one of the pure states is reached we stop following the process.

We note that (for k ̸ = 0, n) the process can jump to k = k + 1 in two ways: either the reproducing individual is of type A, mutation does not take place in the reproduction, and the eliminated individual is of type B; or the reproducing individual is of type B, mutation does occur in reproduction, and the eliminated individual is again of type B. The probability for these to occur are easily evaluated, and the transition probability

P (k → k + 1) := p k is just their sum.
Similarly, the process can jump (for k ̸ = 0, n) from k to k = k -1 in two ways: either the reproducing individual is of type B, mutation does not take place in the reproduction, and the eliminated individual is of type A; or the reproducing individual is of type A, mutation does occur in reproduction, and the eliminated individual is again of type A. The probability for these to occur are easily evaluated, and the transition probability P (k → k -1) := q k is just their sum.

In explicit terms,

p k = k n (1 -µ) n -k n + n -k n µ n -k n = 1 - k n µ + k n -2 µ k n ; (5)
q k = k n µ k n + n -k n (1 -µ) k n = k n 1 -µ + k n -2 µ k n . (6)
Note again that for µ → 0 we recover the probability transitions of the standard Moran model. At difference with the standard Moran model, we have in general (the exception being for k = n/2, if n is an integer) p k ̸ = q k , i.e. a non-symmetric random walk.

In this case we have for the drift

δk = k -k ⟨δk⟩ = p k -q k = µ 1 -2 k n ; (7)
thus for µ ̸ = 0 this dynamics drives the process towards equilibrium, and we expect that hitting the borders requires larger times than in the case of the standard Moran model. In the following, the ratio ξ k = q k /p k will be of interest; it follows from ( 5) and ( 6) and simple algebra that this is given by

ξ k = q k p k = - k(k(2µ -1) -µn + n) (k -n)(-2µk + k + µn) . (8)
Writing k/n = x ∈ [0, 1], this also reads

F (x) = x x -1 1 + 1 2µx -x -µ ; (9)
this is plotted as a continuous function in Figure 1 (this kind of representation is justified for large n, and will be used again in the following).

It should be noted that we always have F (1/2) = 1 and moreover, unless µ is large, F (x) is substantially different from the unit function only for x ≃ 0 or x ≃ 1; in fact, expanding F in a µ series, we get

F (x) = 1 - (1 -2x) x (1 -x) µ + O µ 2 . ( 10
)
Thus for a small mutation rate µ ≪ 1 the properties of the model -or at least those controlled by the parameter ξ k , which we will see is a central one -are expected to be quite similar to those of the standard Moran model, except near the pure states k = 0 and k = n (corresponding to x = 0 and x = 1).

0.2 0.4 0.6 0.8 1 0.8 0.9

1.1

1.2

Figure 1: The function F (x) given by ( 9), plotted here for µ = 0 (black, constant and equal to unity), µ = 10 -3 (blue), µ = 10 -2 (green), and µ = 10 -1 (red). For small but nonzero µ, the curve flattens on F = 1 (the behavior of the standard Moran model) except near the borders α = 0 and α = 1; that is -going back to k and n -except for k ≃ 0 and k ≃ n.

The behavior in this region is of course essential for actually reaching a pure state, so this deviation in a limited region will suffice to provide rather different hitting properties, in particular concerning the mean hitting time.

We also note that the F (x) with smaller µ will be more flat in the central region, and steeper at the borders, see Figure 1. They become of course singular in the limit µ → 0. Remark 4. Note that one could in principles also consider a model in which only mutational fluctuations -and not reproductive ones -are at play; this would be a "purely mutational" Moran model. This has some mathematical interest, but is less interesting biologically, in that in the limit of vanishing mutation rate it yields a trivial model in which whatever initial composition of the population remains fixed. ⊙

Hitting probability: a general discussion

We want to evaluate the hitting probability π(k 0 ) and χ(k 0 ), depending on the initial state k 0 , for the two alleles; these are respectively the probability π(m) that if the initial state is k 0 = m, the system hits the state k = n before passing through the state k = 0; and the probability χ(m) that it hits the state k = 0 before passing through the state k = n. We note preliminarily that by their definition, these quantities satisfy

π(0) = 0 , π(n) = 1 ; χ(0) = 1 , χ(n) = 0 . (11)
We should also note that, in principles, the system could also oscillate indefinitely along the (finite) lattice; thus we are not guaranteed apriori that π(m) + χ(m) = 1.

It is well known that actually random walks on a finite one-dimensional lattice always hit the border (in the case of absorbing boundaries, as for the standard Moran model, this means reaching fixation), as the measure of forever oscillating paths is null. It will be useful to give a simple scheme to prove this, and to actually compute π(k) and χ(k).

We consider a general random walk on the lattice Z n+1 = {0, 1, ..., n}, with transition probabilities p k for k → k + 1, q k for k → k -1, and r k for k → k, with p k + q k + r k = 1, with p k ̸ = 0 ̸ = q k for all k but k = 0 and k = n; π(k) and χ(k) are as above. The situation is trivial for k 0 = 0 (by definition, π(0) = 0, χ(0) = 1) and for k 0 = n (by definition, π(n) = 1, χ(1) = 0) , so we will from now on assume k 0 ̸ = 0, n.

In order to compute π(k) and χ(k), 1 ≤ k ≤ n -1, one filters the process according to its first step. We will first focus on π. Recalling that p k +q k +r k = 1 we write

(p k + r k + q k ) π(k) = p k π(k + 1) + r k π(k) + q k π(k -1) . (12)
Introducing the function (defined for 0 ≤ k ≤ n -1)

η(k) := π(k + 1) -π(k) , (13)
and writing for ease of notation

ξ k := q k / p k , (14)
eq. ( 12) is promptly rearranged as the recursion relation

η(k) = ξ k η(k -1) . (15)
Thus, starting from η 0 = η(0) = π(1) -π(0), we have

η(k) = k h=1 ξ h η 0 . (16)
Now we note that, due to [START_REF] Baake | Ancestral processes with selection: Branching and Moran models[END_REF],

π(k) = π(0) + k-1 h=0 η(h) ; (17)
recalling also that π(0) = 0, specializing to k = n, and using ( 16), we get

π(n) = n-1 k=0 η(k) = n-1 k=0 k h=1 ξ h η 0 . (18)
Using now π(n) = 1, we get in the end

η 0 = n-1 k=0 k h=1 ξ h -1 . ( 19
)
This yields immediately

η(k) =   k h=1 ξ h n-1 k=0 k h=1 ξ h   , (20)
and hence

π(k) =   k-1 h=0 h m=1 ξ m n-1 h=0 h m=1 ξ m   . (21)
In order to evaluate χ(k), we could proceed along the same lines; but we can as well -and more economically -use a symmetry argument. Mapping (k; p k , q k ) → (n -k, q k , p k ) we have just reversed the role of the A and B alleles. This shows at once that we will get

χ(k) =   n-1 h=k h m=1 ξ m n-1 h=0 h m=1 ξ m   . (22)
Note this implies immediately that π(k) + χ(k) = 1 for all k.

Remark 5. Albeit the probabilities π(k) and χ(k) can sometimes -and in particular for the model discussed here -be expressed in terms of hypergeometric functions, in concrete computations it may be convenient to fix n and µ, and use this scheme to provide exact values of the quantities of interest through a symbolic manipulation program rather than going through complex and rather obscure analytic expressions. ⊙ Example 1. For the standard Moran model, see [START_REF] Karlin | On a genetics model of Moran[END_REF], and more generally for any model in which p k = q k , we have ξ k = 1 for all k ̸ = 0, n; then [START_REF] Traulsen | Stochastic dynamics of invasion and fixation[END_REF] and [START_REF] Vogl | A nearly-neutral biallelic Moran model with biased mutation and linear and quadratic selection[END_REF] grant that π(k) grows linearly from π(0) = 0 to π(n) = 1, each increment η(k) being of course given by η(k) = η 0 = 1/n.

Hitting time: a general discussion

We now come to consider the mean hitting time T (k), i.e. the time (measured in time unit such that there is one step per unit) needed to reach either one of the end states k = 0 or k = n. Note that, obviously, T (0) = 0 = T (n). In our discussion we consider n as fixed. We denote the average time spent in state j before hitting the border by a random walk starting in site k as τ (j|k) (this is defined for j ̸ = 0, n); by definition,

τ (j|0) = 0 = τ (j, n) . (23)
Obviously we have

T (k) = τ (1|k) + ... + τ (n -1|k) = n-1 j=1 τ (j|k) . (24)
Now we have to evaluate τ (j|k); we will again (as in the discussion of hitting probability) filter the random walk w.r.t. its first step, which takes the system from state k to state k + 1 with probability p k , to state k -1 with probability q k , and leaves it in state k with probability r k = 1 -p k -q k . Thus -noticing that if the walk starts from k = j by the time we make a first step it will have spent a unit time in site j -we have

τ (j|k) (p k + r k + q k ) = δ jk + p k τ (j|k + 1) + r k τ (j|k) + q k τ (j|k -1) . ( 25
)
We write as above ξ j = q j /p j and also, for ease of notation,

ρ j := 1/p j . (26)
Moreover we introduce

θ(j, k) := τ (j|k + 1) -τ (j|k) . (27)
Then the relation [START_REF] Cerf | Critical population and error threshold on the sharp peak landscape for the Wright-Fisher model[END_REF] reads

θ(j, k) = ξ k θ(j, k -1) -ρ j δ jk . (28)
Thus the θ(j, k) can be determined recursively. We will write

β j := τ (j|1) -τ (j|0) = θ(j, 0) . ( 29
)
With this, we easily get

θ(j, k) =      β j k i=1 ξ i -ρ j k i=j+1 ξ i for j < k , β j k i=1 ξ i -ρ j for j = k , β j k i=1 ξ i for j > k . (30)
Note that here the ρ j = 1/p j are given, while β j are constants yet to be determined.

We next observe that, by [START_REF] Cerf | Quasispecies on class-dependent fitness landscapes[END_REF] and recalling [START_REF] Cerf | The distribution of the quasispecies for a Moran model on the sharp peak landscape[END_REF],

τ (j|k) = θ(j, 0) + ... + θ(j, k -1) . (31)
Recall that, again by [START_REF] Cerf | The distribution of the quasispecies for a Moran model on the sharp peak landscape[END_REF], we have τ (j|n) = 0 as well; thus specializing [START_REF] Berg | Random Walks in Biology[END_REF] to

k = n we get n-1 ℓ=0 θ(j, ℓ) = 0 . (32)
This allows to determine β j , as we now discuss. Using [START_REF] Berg | Random Walks in Biology[END_REF], the relation ( 28) is promptly converted into an equivalent one for the τ . Using τ (j|0) = 0 and (31), we get

τ (j|k) =                        1 + k-1 m=1 ( m i=1 ξ i ) β j -1 + k-1 m=j+1 m i=j+1 ξ i ρ j (j < k), 1 + k-1 m=1 ( m i=1 ξ i ) β j -ρ j (j = k), 1 + k-1 m=1 ( m i=1 ξ i ) β j (j > k). (33)
In order to determine β j , we should once again recall [START_REF] Cerf | The distribution of the quasispecies for a Moran model on the sharp peak landscape[END_REF]. Specializing (33) to the case k = n (which is always in the class j < k, as j = 1, ..., n -1) and requiring this to vanish finally determines β j as

β j =   1 + n-1 m=j+1 m i=j+1 ξ i 1 + n-1 m=1 ( m i=1 ξ i )   ρ j . (34)
Now, inserting this in [START_REF] Patwa | The fixation probability of beneficial mutations[END_REF], one is able to determine explicitly -modulo the problem of evaluating the sums and products appearing in our formulas -the τ (j|k); and hence, using [START_REF] Cerf | Critical population and error threshold on the sharp peak landscape for a Moran model[END_REF], also the mean hitting time T (k).

Needless to say, the evaluation of sums and products can be trivial or extremely difficult (or impossible), depending on the expression of ξ k for the model at hand. Remark 6. These formulas are substantially simplified for symmetric random walks, i.e. for ξ k = 1. In this case, in fact, we get

θ(j, k) = β j -ρ j for j ≤ k , β j for j > k . (35)
This in turn yields

τ (j|n) = n β j -(n -j) ρ j ; (36)
hence requiring this to vanish we just get

β j = n -j n ρ j . (37)
In the end, it results

T (k) = k n-1 j=1 β j - k-1 j=1 (k -j) ρ j . (38)
Combining this with (37), we get

T (k) = k-1 j=1 k n -j j -(k -j) ρ j + n-1 j=k k n -j j ρ j . (39)

⊙

Example 2. For the simple symmetric random walk in [0, n], p = p k = q k = q = p (0 < p ≤ 1/2), and hence all the ξ i appearing in the formula reduce to unity, while ρ j = 1/p j = 1/p. Then

β j = n -j n 1 p , τ (j|k) = j (n -k) n p ,
and in conclusion

T (k) = (n -k) (n -1) 2 p .
Example 3. For the simple non-symmetric random walk in [0, n], p = p k , q k = q (0 < p + q ≤ 1), all the ξ i are equal to ξ = q/p. In this case we get

β j = 1 p 1 -ξ n-j 1 -ξ n , τ (j, k) = 1 p 1 -ξ j 1 -ξ n 1 -ξ n-k 1 -ξ ξ k-j .
In conclusion

T (k) = 1 p [nξ n (1 -ξ) -ξ(1 -ξ n )] 1 -ξ k-n (1 -ξ) 2 (1 -ξ n ) .
6 Hitting probability: the complete Moran model

We will now apply the general discussion of Section 4 to the complete Moran model discussed in Section 3, i.e. to the case ξ k is given by ( 8). We will proceed as mentioned in Remark 5, i.e. analytically with given values for n and µ, and consider the dependence of the emerging results both on the total population number n and on the mutation rate µ.

Dependence on the total population number

The dependence of the form of π(k) on the total population number is quite substantial. This is shown in Figure 2 and Figure 3. In particular, it appears from these that for small n the hitting probability grows basically as in the standard Moran model (i.e. linearly) from π(0) = 0 to π(n) = 1, while for larger n the behavior becomes more and more different; in particular, in the limit of large n the probability of first hitting each of the two pure states is substantially 1/2 unless the initial state is very near to k = 0 or to k = n. For small µ the behavior is very similar to the linear one holding for the standard (µ = 0) Moran model, while as µ is increased the behavior is more and more flattened towards a constant 1/2 probability unless the process starts very near to the borders. [Note that the red line is perfectly horizontal in its core, but may appear tilted due to an optical effect related to the nearby diagonal lines.]

Dependence on the mutation rate

From the biological point of view, the most interesting question is the dependence on the mutation rate µ. This is shown in Figure 4. This shows a transition between the behavior of the standard Moran model, with the first hitting probability of an allele growing linearly with its initial weight in the population, to a situation in which the probability is essentially the same for the two alleles unless the initial state is extremely near to the borders. In other words, nowunless the initial state is extremely strongly skewed towards one of the alleles -first hitting one or the other of the alleles is a random process with equal probabilities. This is easily understood in terms of our previous discussion: in this case the drift will lead the process towards equilibrium (k = n/2), so hitting the pure states will occur as the result of a large fluctuations from equilibrium, which can of course take place in both directions.

This observation also suggest that the mean hitting time will substantially increase in the presence of mutations (unless the mutation rate is negligible); this point will be discussed in the next Section.

The leading control parameter

We have considered the dependence of the hitting probability π(k) (and intrinsically also of χ(k) = 1-π(k)) on the two parameters appearing in the complete Moran model; that is, the population number n and the mutation rate µ. Actually, one can consider a combination of these, i.e.

σ := n µ , (40)
which plays the role of leading control parameter. By this we mean that the behavior of the model depends to a large extent on σ, independently on how this is obtained in terms of n and µ. For small σ the behavior of the model is similar to that of the standard Moran model, while for large σ the behavior is different, showing the flattening of the π(k) function in the core and its steepening at the borders. Equivalently, we can say that the mutation rate is small or large by comparing it with µ 0 = 1/n. We stress that a mutation rate µ 0 means that there is one mutation once each individual has reproduced once (on the average); thus it is a very natural measure unit for the mutation rate.

Note that there is no simple scaling relation in the complete Moran model, as appears from the functional form of F (x) defined in [START_REF] Kimura | Population genetics, molecular evolution, and the neutral theory: selected papers[END_REF]. This is also shown in Figure 5, where we plot f (x) := π(xn) in order to compare the situation for different values of n and µ yielding the same value for σ. Together with an increase of the hitting time, there is a variation in the shape of T (k), see also Figure 7.

Hitting time: the complete Moran model

We can now apply our general formalism developed in Section 5 to the complete Moran model, i.e. to the case where ( 5) and ( 6) hold. We proceed in the same way as in our general discussion above, and will work with given population number n and mutation rate µ, which allow to determine exactly all the quantities of interest.

Dependence on the total population

We will first consider the situation where the mutation rate is fixed, and look at how the mean hitting time varies with the population number. The result of computations for µ = 0.01 and different population numbers are shown in Figure 6. Needless to say, an increase in n leads to larger hitting times. Together with this, we observe that the shape of the function T (k) also changes. In order to compare the results for different n, we consider interpolating functions Θ n (k) for T (k) with total population n, and normalize these so that their maximum is one; in other words, we consider increased, the function T (k) gets a flatter form in the core, and a steeper one near the borders.

ϑ n (x) = Θ n (x * n)/Θ n (n/2) , (41)

Dependence on the mutation rate

The main motivation for this work was to consider the effect of introducing mutations in the standard Moran model; this is considered as a paradigmatic model for neutral evolution, and also used to get a rough estimate of the timescales for allele fixation. Albeit in the complete Moran model fixation is not possible, and is substituted by reaching (temporarily) a "pure" state, the dependence of the mean hitting time on the mutation rate is specially interesting from the point of view of Evolutionary Dynamics (see the discussion in Sect. ?? below).

Our discussion of the hitting probabilities for the two alleles suggest that on the one hand the presence of mutations will lead to an increase of the mean hitting time; on the other hand we have seen that when µ increases the hitting probabilities are essentially independent of the initial state (and equal for the two alleles), unless the initial state is very near to the borders. This suggests that the process spends a long time wandering back and forth in the core of the state space, until a large fluctuation leads to fixation. This picture suggests that the shape of the mean hitting time function T (k) will also be changed as µ is increased, and in particular that its dependence on the initial state k will be not so relevant except when k is near to the borders. The increase in µ leads to a change is the shape of the curve, which becomes wider, and above all to a raising in the curve, i.e. to an increase in the hitting time.

The exact evaluation of T (k) for a population of fixed size (n = 100) and different values of µ, whose results are reported in Figure 8, confirms these expectations.

Reversal time in the presence of mutations

In the absence of mutations, once the process reaches either one of the pure states the situation is fixed forever. If mutations are possible (but staying within the simple frame in which only two alleles can exist) the allele which at some point is absent from the population will appear again due to a mutation; and the process will at some point, after a time T from hitting the pure state, reach the opposite pure state.

The typical time scale for this can be estimated, at first order, as

T = [µ π(1)] -1 T (1) , (42)
where, as before, π(k) is the probability to fix A starting from the state k 0 = k and T (k) is the mean time needed to fixation starting from state k 0 = k.

A comparison between T and the maximal mean hitting time T M (which n π(1)

T

Discussion and conclusions

We have considered a modification of the standard Moran model for neutral evolution, actually already considered in the original paper by Moran [START_REF] Moran | Random processes in genetics[END_REF]. This model takes into account mutational as well as reproductive fluctuations, and is hence called the "complete" Moran model. The presence of mutations forbids that an allele is fixed forever, so the questions to be investigated concern (a) the probability π(k) to first hit one of the pure states or the other if initially there are k individuals with the first allele in a population of n individuals, and (b) the time T (k), under these initial conditions, for a pure state to be first hit. Note that the hit condition corresponds, in the standard (mutation-less) Moran model, to fixation of an allele.

As for point (a), we have observed that the main difference is that while in the standard Moran model π(k) grows linearly from π(0) = 0 to π(1) = 1, in the complete model the behavior is nonlinear. More precisely, there is a sharp increase for k ≃ 0 followed by a substantially flat behavior and by a new sharp increase for k ≃ n. In particular this means that unless the initial conditions are strongly skewed towards one of the pure states, each of them will have a probability 1/2 to be selected for first hit by fluctuations. This effect gets quickly more important with the growth of the mutation rate µ for constant n, and also increases with the growth of the population size n for constant µ.

As for point (b), we have observed two effects. On the one hand, the presence of mutations slows down the hitting of pure states, i.e. for a given population size n and initial condition k, the corresponding mean hitting time T (k) increases (as obvious) with increasing µ. On the other hand, the shape of the function T (k) is also modified as µ increases; here too, for larger µ the function gets steeper near the borders k = 0 and k = n and flatter in the core.

As for the methods employed, we have used exact recursive relations to evaluate π(k) and T (k) at given n. General formulas providing these quantities in terms of special functions (in particular, hypergeometric ones) can be obtained, and are given in the Appendix, but the recursion formulas appear to be more manageable for our purposes.

A Explicit expressions for the complete Moran model

As mentioned above, the hitting probability and mean hitting time for the complete Moran model can be explicitly computed; the resulting formulas are a bit complex and involve special functions, in particular hypergeometric functions and the digamma function [START_REF]Handbook of Mathematical Functions[END_REF].

A.1 Hitting probability

Recalling that the general expression for π(k) is given by eq.( 21), where ξ m := q m /p m , see eq. ( 14), and that for the complete Moran model p m and q m are given by eqs. ( 5) and ( 6) respectively, we obtain -by a direct computation with Mathematica -that in this case the function π(k) can be expressed in terms of the Hypergeometric function the series expansion for small µ yields

T (k) ≈ π 2 6 kn 2 -γ e kn + n 2 ψ (0) (1 -n) -ψ (0) (k -n) -kψ (1) (n) -kn ψ (0) (k) + ψ (0) (1 -n) -ψ (0) (k -n) + π 2 2 kn 2 -γ e n(n + k) + 1 2 kn 3 ψ (2)
(1) -ψ (2) (n) -kn 2 ψ (1) (k) -2ψ (1) (n)

+ n 2 ψ (0) (1 -n) -ψ (0) (k) -ψ (0) (k -n) + kn ψ (0) (k -n) -ψ (0) (k) -ψ (0) (1 -n) µ + O µ 2 . (45)

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: Plots of π(k) for µ = 0.01 and different values of n. For small n the values of π(k) grow nearly linearly (as for µ = 0), while for higher values of n the function π(k) gets flatter in the central region and steeper near the borders.

Figure 5 :

 5 Figure 5: The value of f (x) = π(xn) for α ∈ [0, 1] for σ = nµ = 1 obtained in two different ways: (a) n = 1000 and µ = 0.001 (blue points); (b) n = 100 and µ = 0.01 (red points). The functions obtained are similar but not equal. See the discussion in the text.

Figure 6 :

 6 Figure 6: The mean hitting time T (k) for the complete Moran model with µ = 0.01 and different values of n.Together with an increase of the hitting time, there is a variation in the shape of T (k), see also Figure7.

Figure 7 :

 7 Figure 7: The functions ϑ n (x) defined in eq.(41) for the different values of n considered in Figure6: n = 10 (black curve), n = 25 (blue), n = 50 (orange), n = 100 (red). For higher n, the function shape becomes flatter in the core and steeper near the borders.

Figure 8 :

 8 Figure 8: Dependence of mean hitting time T (k) on the mutation rate for constant population size. Here we consider a population of n = 100 individuals, and plot the function T (k) for different values of the mutation rate µ; in detail:for µ = 0 (black points), corresponding to the standard Moran model; for µ = 0.001 (red points); for µ = 0.005 (green points); and for µ = 0.01 (blue points). The increase in µ leads to a change is the shape of the curve, which becomes wider, and above all to a raising in the curve, i.e. to an increase in the hitting time.

2 F 1

 1 (a, b; c; z) := ∞ k=0 P(a, k) P(b, k) P(c, k)z k k! , where P(a, k) is the function 2 P(a, k) := Γ(a + k) Γ(a)with Γ the Euler Gamma function.In order to have a lighter notation, we will setF(b; c; z) := 2 F 1 (1, b; c; z)A.3 Mean hitting time for small mutation rateFinally, let us consider the series expansion for small mutation rate µ. Reverting to the usual polygamma notation[START_REF]Handbook of Mathematical Functions[END_REF], i.e. withψ (k) (z) := d k Ψ(z) dz k = d k+1 dz k+1 log[Γ(z)] ,

Table 1 :

 1 Comparison of reversal time T and of the maximal mean hitting time T M = T (n/2) in the complete Moran model for µ = 0.01 and different values of the total population n. In all cases, T is two orders of magnitude larger than T M . corresponds to T (n/2)) is provided in Table1. It appears from this that T /T M ≈ 10 2 in all considered cases. Remark 7. As stressed in the Introduction, the complete Moran model is an ergodic Markov chain, admitting a single equilibrium state; in particular, no fixation is possible. One can however consider, from the mathematical point of view, a related model with absorbing boundary condition. This corresponds to stopping the process once a pure state is reached. In biological terms, this would correspond to the -extremely unrealistic -hypothesis that mutations are not possible if the system is in a pure state, and in this sense such a model should be seen as a mathematical one, possibly with real-world applications but surely not relevant in the context we have been discussing and which provides motivation for the present work. It is obvious that all the computations presented above go unaltered in such a case, albeit with a different meaning: hitting probabilities become fixation probabilities, and hitting time becomes fixation time.⊙

			1)	T	T M
	10	0.11	30.63	27845	68
	25	0.06	131.85	219750	496
	50	0.06	548.51	914183	2558
	100	0.10 5290.58 5290580 19977

These problems were already considered in the literature, but the available results are mostly of numerical nature (we will not discuss them); we will instead adopt an analytic approach, and provide exact formulas -albeit in recursive form. See however the Appendix for closed-form expressions in terms of special functions.

This is also known in the literature as Pochhammer symbol; a more convenient name for it is the rising factorial, and it is associated to Stirling numbers of the first kind[START_REF] Knuth | Two notes on notation[END_REF].

Acknowledgements

This work was started while I was on sabbatical leave from Università degli Studi di Milano, enjoying the hospitality and the relaxed work atmosphere of SMRI; in particular I thank L. Peliti for sharing with me his notes on Knuth's "Two notes on notation" [35]. This work was supported by "iniziativa specifica MMNLP" of INFN; I am also a member of IAMP and of ISNMP. My work is supported by GNFM-INdAM.

and moreover write ν = 1 -n ; η = µ -1 µ .

We then have

A.2 Mean hitting time

Recalling that the general expression for T (k) is given by eq.( 39), where ρ j is given by eq.( 26), and that for the complete Moran model p j is given by eq.( 5), we obtain -again by a direct computation with Mathematica -an explicit formula. In order to simplify notation, we will denote by Ψ(z) the digamma function -which is often denoted as ψ (0) (z) -that is,

with Γ(z) the usual Euler Gamma function. Moreover, we denote by γ e the Euler-Mascheroni gamma constant, γ e ≈ 0.577216 , and write for short

With these notations, it results

This might be somewhat simplified using the relations

and its asymptotic behavior for large n can be inferred from Ψ(z) ≈ log(z) -1 2z

which holds for |z| → ∞ for |arg(z)| < π -ε.


Angelo Koudou

Yves Ismaël

Ngounou Bakam

Denys Pommeret

Lancaster

Yves Ngounou

Lancaster copulas

Keywords: Copula, Lancaster distributions, Maxiset theory, Minimax theory, Non-parametric estimation, Orthogonal polynomials

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

Introduction

Consider a d-random vector X = (X 1 , • • • , X d ) T with joint distribution function H and marginal distribution functions F 1 , • • • , F d , that we assumed to be continuous. According to Sklar's Theorem [START_REF] Sklar | Fonction de répartition dont les marges sont données[END_REF]), there exists a unique d-variate function C such that

H(x 1 , • • • , x d ) = C(F 1 (x 1 ), • • • , F d (x d )).
(1)

The function C is called the copula (function) associated to X. The copula is a joint distribution function on I d = [0, 1] d with uniform margins and satisfying

C(u 1 , • • • , u d ) = H(F -1 1 (u 1 ), • • • , F -1 d (u d ))
, where, for j = 1, • • • , d, F -1 j (u) = inf{x; F j (x) ≥ u}, is the quantile function of F j . Assuming that for j = 1, • • • , d, F j is differentiable, we can express the joint density h of X as

h(x 1 , • • • , x d ) = c (F 1 (x 1 ), • • • , F d (x d )) d j=1 f j (x j ), (2)
where for j = 1, • • • , d, f j is the marginal density of X j and where

c = ∂ d C ∂x 1 • • • ∂x d , (3)
is called the copula density of X. Copulas and density copulas have a large spectra of applications as described for instance in [START_REF] Joe | Dependence modeling with copulas[END_REF]. For a general account on copulas, see, e.g. [START_REF] Roger | An introduction to copulas[END_REF]. Several copula families or classes are based on the Sklar's Theorem and have been derived from certain families of multivariate distribution. Various constructions have been proposed, such as, for instance, Archimedean copulas (see [START_REF] Hofert | Archimedean copulas in high dimensions: Estimators and numerical challenges motivated by financial applications[END_REF]), Vines (see [START_REF] Aas | Paircopula constructions of multiple dependence[END_REF]), Elliptical copulas (see [START_REF] Frahm | Elliptical copulas: applicability and limitations[END_REF]) which are simply the distribution functions of componentwise transformed elliptically distributed random vectors, or the Extreme-value copulas (see [START_REF] Gudendorf | Extreme-value copulas[END_REF]) which relies to the copulas of random vectors distributed according to multivariate extreme value distributions.

In this paper we propose a novel copula class called Lancaster copulas. This class is based on the class of distributions satisfying an orthogonality condition involving orthogonal polynomials with respect to their marginal laws. That is the copula of random vectors whose joint distributions are in the class of Lancaster distributions. The reader is referred to Lancaster [START_REF] Ho Lancaster | The structure of bivariate distributions[END_REF][START_REF] Ho Lancaster | Correlations and canonical forms of bivariate distributions[END_REF][START_REF] Oliver | Joint probability distributions in the meixner classes[END_REF]) the synthesis by [START_REF] Efoevi | Problèmes de marges et familles exponentielles naturelles[END_REF][START_REF] Efoévi | Probabilités de lancaster[END_REF] and the characterization by [START_REF] Pommeret | A characterization of lancaster probabilities with margins in a multivariate additive class[END_REF]; [START_REF] Koudou | Theory & methods: A construction of lancaster probabilities with margins in the multidimensional meixner class[END_REF] for more details. To recapitulate this class of Lancaster distributions we generally follow the standard notation in that area, as given by these authors.

If X and Y are two random variables with probability measures µ and ν, and if (P n ) n∈N (respectively (Q n ) n∈N ) is a µ (respectively ν) complete orthonormal basis of polynomials, the the joint density of (X, Y ), says Consider bivariate distributions σ(dx, dy) with margins µ and ν defined on spaces X and Y respectively. Also, we denote by (X, Y ) a random vector with distribution σ. If there exists bi-orthonormal polynomials, (P n ) n∈N and (Q n ) n∈N , associated to µ and ν, respectively, such that

E (P n (X)P k (X)) = δ n,k ; E (Q n (Y )Q k (Y )) = δ n,k and E (P n (X)P k (Y )) = ρ n δ n,k ,
for n, k ≥ 0, then the class of bivariate distributions σ are known as Lancaster probabilities and the sequence ρ n = E (P n (X)P k (Y )) is known as the Lancaster sequence of σ.

These types of distributions are widely used in various areas of probability and statistics. For instance we find them in the literature concerning stationary Markov processes ( see [START_REF] Jakob | Crosscorrelation functions of amplitude-distorted gaussian signals[END_REF], [START_REF] Barrett | An expansion for some second-order probability distributions and its application to noise problems[END_REF]Lampard (1955), Brown (1958), [START_REF] Wong | On polynomial expansions of second-order distributions[END_REF], [START_REF] Wong | The construction of a class of stationary markoff processes[END_REF], etc.) and also in canonical analysis (see [START_REF] Dauxois | Une extension de l'analyse canonique. quelques applications[END_REF]Pousse (1975) andBuja (1990)).

They are also use for disjunctive kriging in geostatistics (see [START_REF] Droesbeke | Méthodes bayésiennes en statistique[END_REF] and the reference given there). There are also a few more recent references: [START_REF] Griffiths | Stochastic processes with orthogonal polynomial eigenfunctions[END_REF] on reversible Markov processes whose proper functions are orthogonal polynomials, [START_REF] Diaconis | Gibbs sampling, exponential families and orthogonal polynomials[END_REF] where new techniques are proposed to calculate convergence speeds to stationary measurement in the context of bivariate Gibbs sampling and where it is shown that these techniques give good results if the target measurements are Lancaster distributions. In the area of discrete probabilities, [START_REF] Diaconis | Exchangeable pairs of bernoulli random variables, krawtchouck polynomials, and ehrenfest urns[END_REF] provide an interpretation of Lancaster probabilities of binomial margins using the generalized Ehrenfest ballot box model and [START_REF] Robert | Orthogonal polynomial kernels and canonical correlations for dirichlet measures[END_REF] for Dirichlet measurements.

The plan of the paper is as follows: Section 2 introduces the Lancaster copula and their densities. Some examples of these functions are presented in 3. Non-parametric estimators are proposed in Section 4. provides an extension to the multivariate case. The proofs and some mathematical and statistical properties are given in Appendix.

Lancaster copulas and their densities

Bivariate case

Let µ and ν be two probability measures on X and Y with density f X and f Y with respect to Lebesgue such that there exists an open set Θ satisfying e θx µ(dx) < ∞ and e θy ν(dy) < ∞, for all θ ∈ Θ, respectively.

Assume that (X, Y ) is a random vector with margins µ and ν and with joint distribution σ. We say that σ is a Lancaster probability with margins µ and ν if it satisfies the following decomposition:

σ(dx, dy) = n∈N ρ n P n (x)Q n (y)µ(dx)ν(dy),
where (P n ) n and (Q k ) k are orthonormal polynomials sequences with respect to µ and ν. Here we assume that µ, ν and σ are absolutely continuous and σ has density h. This yields the expansion,

h(x, y) = f X (x)f Y (y) n∈N ρ n P n (x)Q n (y) (4)
Identifying this expression ( 12) with (2) from Sklar's Theorem ??, we obtain the Lancaster copula density given by

c(u, v) = n∈N ρ n P n (F - X (u))Q n (F - Y (v)), (5)
where F -denotes the pseudo inverse of F . By integration, we get the Lancaster copula function as follows

C(u, v) = n∈N ρ n F -1 X (u) 0 P n (x)µ(dx) F -1 Y (v) 0 Q n (y)ν(dy) (6)
Remark 1. The Lancaster copula function ( 6) and the Lancaster copula density function ( 5) are completely specified by their marginal measures and the correlations ρ n between P n (X) and Q n (X). In terms of canonical analysis, the variables P n (X) and Q n (Y ) are the canonical variables. So, is is important to note that this Lancaster sequence (ρ n ) n∈N characterizes the copula and contains the polynomial correlations between all the marginal random variables.

Examples of Lancaster Copulas

Bivariate normal Lancaster copula

Let (X, Y ) ∼ N (0, Σ) follow a bivariate normal distribution with mean vector 0 and covariance Σ = 1 ρ ρ 1 . For n ∈ N, let

H n (x) = (-1) n n! √ 2π exp 1 2 x 2 d dx n exp - 1 2 x 2 , x ∈ R
denote the nth Hermite polynomial. It is well-known that the polynomials {H n : n ∈ N} are orthonormal with respect to the standard normal distribution µ and form a complete orthonormal basis for the Hilbert space L 2 (µ). [START_REF] Ov Sarmanov | Probabilistic properties of bilinear expansions of hermite polynomials[END_REF] consider series expansions of the form

ϕ X,Y (x, y) = ϕ X (x)ϕ Y (y) ∞ n=0 ρ n n! H n (x)H n (y), (x, y) ∈ R 2 (7)
as the Lancaster expansion of (X, Y ) where the sequence of real numbers

(ρ n ) n∈N satisfies ∞ n=0 ρ 2 n < ∞; ϕ X,Y denoted the density function of (X, Y )
and ϕ X , ϕ X their marginal densities. We also denote by Φ the probability distribution function and ϕ -1 the associated quantile function. [START_REF] Ov Sarmanov | Probabilistic properties of bilinear expansions of hermite polynomials[END_REF] proved that for the expansion (7) to be nonnegative, and therefore to be a valid probability density function, it is necessary and sufficient that the sequence (ρ n ) n∈N be the moment sequence of a random variable supported on the interval [-1, 1] and set ρ n = ρ n , if n even 0, otherwise Thereby, the Lancaster copula density associated to (X, Y ) is given by

c(u, v) = ∞ n=0 ρ n n! H n Φ -1 X (u) H n Φ -1 Y (v) , (u, v) ∈ I 2
and its Lancaster copula function is as follows

C(u, v) = n∈N ρ n n! Φ -1 X (u) 0 H n (x)µ(dx) Φ -1 Y (v) 0
H n (y)ν(dy)

Bivariate gamma Lancaster copula

The Lancaster expansion for a bivariate gamma distribution, which was derived by Sarmanov (1970a,b) can be stated as follows (see [START_REF] Kotz | Continuous multivariate distributions[END_REF], page 437 -438). Let α > -1 and for n ∈ N, the classical Laguerre polynomial is defined by

L (α) n (x) = x -α e -x n! d n dx n e -x x n+α = (α + 1) n n! n k=0 (-n) k (α + 1) k x k k! , x > 0 where (α) j = Γ(α+n) Γ(α) = α(α + 1) • • • (α + j -1
) denotes the rising factorial. By standardizing the classical Laguerre polynomial, we obtain the orthonormal version

L (α) n (x) = L (α) n (x)/q k where q 2 k = L (α) n 2 = (α + 1) n n!
Let λ ∈ (0, 1) and let α and β satisfy α ≥ β > 0. Sarmanov (1970a,b) derived for certain bivariate gamma random variables (X, Y ) the joint probability density function,

f XY = f X (x)f Y (y) ∞ n=0 a n L (α-1) n (x)L (β-1) n (y), x, y > 0 where a n = (β) n (α) n 1/2 λ n (8)
The corresponding marginal density functions are

f X (x) = 1 Γ(α) x α-1 exp(-x) and f Y (y) = 1 Γ(β) y β exp(-y),
which we recognize as the density functions of one-dimensional gamma random variables with index parameters α and β respectively. [START_REF] Kotz | Continuous multivariate distributions[END_REF] proved that if α = β then the density function ( 8) reduces to the Kibble-Moran bivariate gamma density function with Corr(X, Y ) = λ and ( 8) represents the Lancaster expansion for (X, Y ). Therefore, the Lancaster copula density associated to the bivariate gamma distribution (X, Y ) as follows

c(u, v) = ∞ n=0 λ n L (α) n F -1 X (u) L (α) n F -1 Y (v) , (u, v) ∈ I 2 , α > 1
and its Lancaster copula function is given by

C(u, v) = n∈N λ n F -1 X (u) 0 L (α) n (x)µ(dx) F -1 Y (v) 0 L (α) n (y)(dy), (u, v) ∈ I 2 , α > 1
where µ is the gamma measure.

Bivariate beta Lancaster copulas

Let considering the triangular bivariate beta distribution (X, Y ) with

σ(dx, dy) = a + b B(a, b) x a-1 y b-1 1 A (x, y) where A = {(x, y) : x, y ≥ 0, x+y ≤ 1, 1 ≤ a, b < ∞}
and B the beta function. This of joint distribution can be found in [START_REF] Buja | Remarks on functional canonical variates, alternating least squares methods and ace[END_REF], page 1049. The marginal distribution of X is β a,b+1 i.e µ(dx) ∝ β a,b+1 (dx) which has the Jacobi polynomials with paramter p = a + b and q = a on the interval (0, 1) as orthogonal polynomials where β is the univariate beta distribution. Similarly, ν(dy) ∝ β a+1,b (dy). Thus, the Lancaster sequence is given by

ρ n = (-1) n √ ab (a + n)(b + n)
and so an explicit expression of ( 5) and ( 6) follows without difficulty. Bivariate quadratic natural exponential families Lancaster proposed a construction of Lancaster distributions associated to natural eponential families.

An estimation procedure

For any positive integer N we can define the following N -th order approximations

c [N ] (u, v) = N n ρ n P n (F -1 X (u))Q n (F -1 Y (v))
and

C [N ] (u, v) = N n ρ n F -1 X (u) 0 P n (x)µ(dx) F -1 Y (v) 0 Q n (y)ν(dy) If we observe an m-sample (X 1 , Y 1 ), • • • , (X m , Y m )
, of iid random data, having joint distribution function F XY and joint density f XY of (X, Y ), we can estimate the ρ n by

ρ n = 1 m m i=1 P n (X i )Q n (Y i ).
Therefore, a N -th order non-parametric estimator of the Lancaster copula density c is given by

c [N ] (u, v) = N n ρ n P n ( F -1 X (u))Q n ( F -1 Y (v)) (9)
where F -1 Z is the empirical quantile estimation of the random variable Z. By integration, we get a non-parametric estimator of the Lancaster copula function as follows

C [N ] (u, v) = N n ρ n F -1 X (u) 0 P n (x)µ(dx) F -1 Y (v) 0 Q n (y)ν(dy) (10)
Remark 1.

• Marginal quantiles can be estimate through kernel method in the following way

F -1 Z (t) = inf{z : F Z (z) ≥ t} with F Z (z) = 1 m m i=1 K( z -Z i h m )
where k is a kernel on R, that is, a symmetric density function, K(z) = z -∞ k(s)ds and h m is a sequence of strictly positive smoothing parameter which controls the smoothness of the estimator, with lim m→∞ h m = 0 and lim m→∞ mh m = ∞.

• The estimator of C and c have the same simple form whether the margins are known or not.

Properties of such estmators are summarized in Appendix. In particular, we have a truncation parameter selection, asymptotic properties and minimax results in the bi-orthogonal case.

Multivariate extensions

To extend the definition of the Lancaster copula in multivariate case, we introduce the following multivariate notations. Given

x = (x 1 , • • • , x d ) ∈ R d and n = (n 1 , • • • , n d ) ∈ N d , we write x n = x n 1 1 × • • • × x n d
d and we call the integer |n| = n 1 + • • • + n d the order of n. A |n|th degree polynomial, say P n , has the following form

P n (x) = k∈N d |k|≤|n| α k x k ,
where at least one real α k is non null for |k| = |n|. Let X and Y be two random variables on R d with probability measures µ and ν respectively. Let (P n ) n∈N d and (Q n ) n∈N d be two basis of orthonormal polynomials with respect to the distributions of X and Y , respectively.

Here and subsequently, P n and Q n denote the polynomials of the |n|th and |k|th degrees.

The joint distribution σ is called a Lancaster probability if

E (P n (X)Q k (Y )) vanishes for n ̸ = k.
that is

P n (x)Q k (Y )σ(dx, dy) = ρ n δ n,k
where δ nk is the Kronecker's delta and σ the joint distribution of (X, Y ). ρ n is called a Lancaster sequence. It is well known that the joint probability measure σ has a polynomial expansion under the condition that a belongs to the Hilbert space L 2 (µ ν). Then we can write

σ(dx, dy) = n∈N d ρ n P n (x)Q n (y)µ(dx)ν(dy) (11)
We assume that µ, ν and σ are absolutely continuous with respect to Lebesgue measure and we denote their corresponding probability density functions by f X ,f Y and f XY respectively. This yields the expansion,

f XY (x, y) = f X (x)f Y (y) n∈N d ρ n P n (x)Q n (y) (12)
Identifying ( 2) and ( 12) we obtain the relation

c(u, v) = n∈N d ρ n P n F -1 X (u) Q n F -1 Y (v) , (13)
where

u = (u 1 , • • • , u d ), v = (v 1 , • • • , v d ) ∈ I d = [0, 1] d and F -1 Z (z) = (F -1 Z 1 (z 1 ), • • • , F -1 Z d (z d ))
We deduce an expression of the copula C, as follows The lancaster copula density can be approximated by the following truncated sum

C(u, v)) = n∈N d ρ n F -1 X 1 (u 1 ) 0 • • • F -1 X d (u d ) 0 P n (x)µ(dx) F -1 Y 1 (v 1 ) 0 • • • F -1 Y d (v d ) 0 Q n (
c [N ] (u, v) = |n|≤N ρ n P n F -1 X (u) Q n F -1 Y (v) ,
where N is an arbitrary integer. The square norm of the error of such approximation can be easily expressing using the orthogonality of the polynomials. We obtain

(c -c [N ] ) 2 (dx, dy) = |n|>N ρ 2 n .
We can also the estimate the Lancaster sequence ρ n by the following empirical correlations

ρ n = 1 m m i=1 P n (X(i)) Q n (Y (i)) .
Hence, we get the following unbiased estimation of c c [N ] (u, v)

= |n|≤N ρ n P n F -1 X (u) Q n F -1 Y (v) ,
The same reasoning makes it possible to obtain an analog result for the Lancaster copula function. All properties in bivariate case can also be extended in the multivariate case.

Discussion

  y)ν(dy)5.1. EstimationLet assume ((X(1),Y (1)) , • • • , ((X(m), Y (m)) are i.i.d with Lancaster distribution such that ((X(k), Y (k)) are distributed as (X, Y ). X(k) denote the vector (X 1 (k), • • • , X d (k)) and Y (k) denote the vector (Y 1 (k), • • • , Y d (k)).