I Rodriguez-Donis

J Acosta-Esquijarosa

V Gerbaud

E Pardillo-Fondevila

X Joulia

Separation of n-hexane -ethyl acetate mixtures by azeotropic batch distillation with heterogeneous entrainers

Keywords: batch distillation, heterogeneous entrainer, azeotropic distillation

In this article, a systematic study of the separation of the n-hexane -ethyl acetate mixture with an entrainer by heterogeneous azeotropic batch distillation is performed. Based upon the thermodynamic behaviour of the ternary mixtures, potential entrainers partially miscible with one or two original azeotropic components are chosen. In all cases, the entrainer adds a heterogeneous binary or ternary azeotrope that is the lowest boiling point in the ternary diagram. Therefore, it leaves the column by the overhead stream which is subcooled to get two liquid phases in the decanter. The phase with the highest amount of the original component is removed as distillate product whereas the entrainer -rich phase is continuously refluxed to the column. Considering methanol, acetonitrile, water and nitromethane as heterogeneous entrainers, screening was performed based on the 2 composition of the unstable heteroazeotropic mixture, the ratio of both liquid phases in the condensed top vapour and the purity of the distillate product determined by the liquidliquid envelope at the decanter temperature. The process feasibility analysis is validated by using rigorous simulation with the batch process simulator ProSimBatch. Simulation results are then corroborated in a bench experimental column for the selected entrainer, showing several advantages of heterogeneous batch distillation compared to homogeneous systems.

Introduction

Batch processes are again becoming important because of recent expansion of the pharmaceutical and speciality chemical industries in many developing countries. Besides, recycling of liquid waste streams is becoming a key issue in all processing plants to meet stricter environmental regulations. Batch distillation is a widely used separation method in batch processes because of its inherent operational flexibility due to some extent to its time dependent operation. The separation of azeotropic mixtures or close boiling components is a challenging task in many chemical processes as it is impossible using a single conventional column or as a pressure swing distillation is an uneconomical process.

Alternative techniques based on the addition of an auxiliary substance in the original mixture have been developed and are usually known as azeotropic and extractive distillation.

Process synthesis and design of these non-conventional distillation processes proceed in two steps. The first step -process synthesis -is the selection of one or more candidate entrainers along with the computation of thermodynamic properties like residue curve maps that help assess many column features such as the adequate column configuration and the corresponding product cuts sequence. The second step -process design -involves the search for optimal values of batch distillation parameters such as the entrainer amount, reflux ratio, boiler duty and number of stages. The complexity of the second step depends on the solutions obtained at the previous level, because efficiency in azeotropic and extractive distillation is largely determined by the mixture thermodynamic properties that are closely linked to the nature of the entrainer. Hence, we have established a complete set of rules for the selection of feasible entrainers for the separation of non ideal mixtures by homogeneous azeotropic batch distillation [START_REF] Rodriguez-Donis | Entrainer Selection Rules for the Separation of Azeotropic and Close Boiling Temperature Mixtures by Homogeneous Batch Distillation[END_REF]. In contrast to batch homogeneous distillation processes, studies on heterogeneous batch distillation were scarce [START_REF] Düssel | Separation of Azeotropic Mixtures by Batch Distillation Using an Entrainer[END_REF][START_REF] Köhler | Zur diskontinuierlichen Rektifikation azeotroper Gemische mit Hilfsstoffeinsatz[END_REF][START_REF] Skouras | Separation of ternary heteroazeotropic mixtures in a closed multivessel batch distillation-decanter hybrid[END_REF] until a new method for the synthesis of heterogeneous distillation in a batch rectifier was defined and a complete set of rules for the selection of heterogeneous entrainers for the separation of non-ideal mixtures was established [START_REF] Rodriguez-Donis | Heterogeneous Entrainer Selection Rules for the Separation of Azeotropic and Close Boiling Temperature Mixtures by Heterogeneous Batch Distillation[END_REF][START_REF] Rodríguez-Donis | Feasibility of Heterogeneous Batch Distillation[END_REF].

Simulation and experimental results have shown key differences between homogeneous and heterogeneous batch distillation [START_REF] Rodríguez-Donis | Feasibility of Heterogeneous Batch Distillation[END_REF]. Indeed, some advantages of heterogeneous batch processes compared to homogeneous systems are : i) more suitable candidate entrainers and hence, more design alternatives for the separation of non-ideal mixtures, ii) simplified distillation sequences thanks to the liquid -liquid phase split occurring in some parts of the column and in the decanter, iii) the addition of a smaller amount of entrainer to the original mixture, iv) the use of more flexible reflux policies through any combination of the entrainer -rich phase and the distillate -rich phase, v) the still path can cross unstable separatrices allowing the separation of components located in different basic distillation regions. The overall cost reduction generated by these advantages makes heterogeneous batch distillation an attractive option for the separation of non-ideal mixture generated by specialty industries.

The aim of this paper is to study the performance of heterogeneous batch distillation (HBD) for the separation of the mixture n-hexane -ethyl acetate. This azeotropic mixture is often used in pharmaceutical industries for separating bioactive substances by liquid chromatographic processes. The separation of this mixture by batch distillation was first studied by Rodriguez-Donis et al. [START_REF] Rodriguez-Donis | Separation of nhexane-ethyl acetate mixtures by azeotropic distillation[END_REF], by using acetone as a homogeneous entrainer. As the separation by using a single batch distillation process with acetone was not possible, a hybrid process was proposed involving a homogeneous azeotropic distillation column with acetone and a liquid -liquid extraction column with water. To circumvent this complex process, we propose in this paper a simpler process by using heterogeneous azeotropic batch distillation. Based on entrainer selection rules established by Rodriguez-Donis et al. [START_REF] Rodriguez-Donis | Heterogeneous Entrainer Selection Rules for the Separation of Azeotropic and Close Boiling Temperature Mixtures by Heterogeneous Batch Distillation[END_REF], which consider the thermodynamic behaviour of multicomponent mixtures, an initial list of potential entrainers partially miscible with one or both original azeotropic components are selected. From the analysis of residue curve map features, in particular the decanter phases composition and the distillation region separatrices curvature, the sequence of products to be drawn from the rectifier column can be predicted. The number of distillation tasks in the whole process is set by the products targeted purity.

Selection of a suitable entrainer for the separation of n-hexane -ethyl acetate mixtures by heterogeneous azeotropic batch distillation.

Liquid -vapour equilibrium data for the binary mixture n-hexane -ethyl acetate have been scarcely reported in the scientific literature [START_REF] Gmehling | Azeotropic Data[END_REF]. Those authors reported a minimum boiling temperature azeotrope with an average mole fraction of n-hexane equal to 0.6565 and a temperature of 65.15°C at 101.3 kPa. A more detailed experimental study, at atmospheric pressure, of liquid -vapour equilibrium was performed by Acosta et al. [START_REF] Acosta | Vapour-Liquid Equilibrium of the Ternary System Ethyl Acetate + Hexane + Acetona at 101.32 kPa[END_REF].

Their estimation of the azeotropic mixture corresponds to a temperature of 64.85°C with a molar fraction of n-hexane of 0.657 with an average experimental error for the temperature and composition measurements of 0.02 °C and 0.003 mole fraction, respectively.

The candidate homogeneous and heterogeneous entrainers are investigated applying the rules reported by Rodriguez-Donis et al. [START_REF] Rodriguez-Donis | Entrainer Selection Rules for the Separation of Azeotropic and Close Boiling Temperature Mixtures by Homogeneous Batch Distillation[END_REF][START_REF] Rodriguez-Donis | Heterogeneous Entrainer Selection Rules for the Separation of Azeotropic and Close Boiling Temperature Mixtures by Heterogeneous Batch Distillation[END_REF]. The existence of binary and ternary azeotropes is sought by using literature data [START_REF] Gmehling | Azeotropic Data[END_REF] and by simulating liquid -vapour and liquid -liquid -vapour equilibrium by using the BibPhyAddIn Macro in an Excel spreadsheet [10]. Confirming the need for a complex hybrid process in case a homogeneous entrainer like acetone is used [START_REF] Rodriguez-Donis | Separation of nhexane-ethyl acetate mixtures by azeotropic distillation[END_REF], no homogeneous entrainer matching the reported rules has been found whereas some suitable heterogeneous entrainers have been identified; the composition and temperature of their azeotropic mixtures with the original components are displayed in Table 1 [START_REF] Gmehling | Azeotropic Data[END_REF].

As seen in Table 1, the boiling temperature of the possible entrainers is either lower or higher than that of the original components. All candidates form a binary heterogeneous azeotrope with n-hexane. Methanol and acetonitrile add a homogeneous minimum boiling azeotrope with ethyl acetate. Nitromethane forms a zeotropic mixture with ethyl acetate, whereas water is the only component that forms an additional heterogeneous azeotrope with ethyl acetate and also adds a heterogeneous ternary azeotrope to the resulting mixture.

Experimental liquid -vapour equilibrium data are available for all binary entrainer-ethyl acetate mixtures [START_REF] Gmehling | Vapour-Liquid Equilibrium Data Collection[END_REF]. However, for the binary system entrainer -n-hexane that exhibits a significant miscibility gap, the only experimental thermodynamic information is given in Table 1, except for methanol for which liquid -vapour equilibrium data were reported by Gmehling et al. [START_REF] Gmehling | Vapour-Liquid Equilibrium Data Collection[END_REF]. All previous experimental data are used to regress binary interaction parameters for the UNIQUAC and NRTL models with ChemCAD IV minimizing the mean square deviation between the experimental and calculated mole fraction in the liquid -phase (Table 2). Then, ternary liquid -liquid -vapour equilibrium is determined by using these binary coefficients with BibPhyAddIn [10]. As a result of the use of these binary interaction parameters regressed on binary systems to predict the behavior of a ternary system, we have detected the existence of the experimental heterogeneous ternary azeotrope in the mixture n-hexane -ethyl acetate -water, which has the lower boiling point in the system (unstable node).

For the synthesis of heterogeneous batch distillation the liquid -liquid envelope at the decanter temperature is considered in addition to the residue curve map. Therefore, the binary interaction parameters used in predicting liquid -liquid equilibrium are estimated from binary heterogeneous azeotrope or liquid -liquid equilibrium data [START_REF] Gmehling | Azeotropic Data[END_REF]10]. Table 3 shows the calculated purity of original components in each phase split at 25°C for all heterogeneous azeotropes reported in Table 1. The thermodynamic models and binary coefficients used in the calculation of the liquid -liquid -vapour equilibrium, liquidliquid equilibrium at 25°C and the separatrices are reported in Table 2.

From the thermodynamic information given in Table 1, Table 2 andTable 3, the residue curve maps are drawn in figure 1 for each entrainer and the batch distillation task sequence required to perform the separation of the original components is deduced using published rules [START_REF] Rodriguez-Donis | Heterogeneous Entrainer Selection Rules for the Separation of Azeotropic and Close Boiling Temperature Mixtures by Heterogeneous Batch Distillation[END_REF]. The residue curve maps of the ternary mixture n-hexane -ethyl acetateentrainer displayed in Figure 1 indicate the temperature and the topological stability of all singular points. As the heteroazeotropes are unstable nodes, they are recovered in the decanter where the phase rich in one of the original components L W is removed as distillate, whereas the other original component is retained into the still at the process end.

In some cases, the composition of the phase rich in original component defined by the liquid -liquid equilibrium at the decanter temperature may not meet purity requirements; then, further distillation tasks are needed.

One of the key thermodynamic parameter in heterogeneous batch distillation is the decanter split ratioω. It is defined by the position of the heteroazeotrope composition on its liquid -liquid tie line at the decanter temperature or alternatively by the mole ratio of the entrainer -rich phase L R to the overall liquid phase L° into the decanter as follows:

o R L L ⋅ = ω (1)
Then, the molar quantity of entrainer-weak / original component -rich phase W L is

given by equation ( 2):

o W L L ⋅ - = ) 1 ( ω (2)
As stated by Rodriguez-Donis et al. [START_REF] Rodríguez-Donis | Feasibility of Heterogeneous Batch Distillation[END_REF], the reflux policy to be used is strongly influenced by the split ratio ω at the decanter. If the liquid reflux needed at the top of the column is lower than L R , then the distillation can be performed by using only the reflux of entrainer -rich phase. Otherwise, the separation of original components requires the reflux of a combination of both decanted phases.

Nitromethane shows the simplest residue curve map with one unstable curved separatrix dividing the triangle in two basic distillation regions. Methanol and acetonitrile give rise two binary azeotropic mixtures and three distillation regions that are bounded by two unstable curved separatrices. Water shows the most complicated residue curve maps, due to the presence of a ternary azeotrope and a miscibility gap with both the n-hexane and the ethyl acetate component. In all four cases, the heteroazeotrope (binary or ternary) has the lowest boiling temperature of the system. As it can be seen in Table 3, all entrainers except water provide the n-hexane -rich phase L W as distillate product with a purity better than 0.91. Water is not a desirable entrainer because of the existence of ternary azeotrope whose n-hexane -rich phase has a water purity much lower (0.70). Considering in Table 3 the split ratio ω at the decanter (inside the range 0.3-0.6) and the purity of both liquid phases at 25°C, methanol and acetonitrile are better entrainers than nitromethane. Methanol is a light entrainer that would result in a ternary system with a low boiling temperature that in turn will likely increase the cost of the condenser cooling system and the make-up of fresh entrainer to be added in the next batch. Besides, the mole ratio of the product-to-be nhexane in the methanol rich phase (x II n-hexane in Table 3) is higher than the one calculated for the acetonitrile. Therefore, acetonitrile seems the best entrainer for the separation of nhexane -ethyl acetate by heterogeneous azeotropic batch distillation.

Simulation of the separation of the mixture n-hexane -ethyl acetate by heteroazeotropic batch distillation using acetonitrile.

We study the separation of n-hexane -ethyl acetate mixture by using acetonitrile as a heavy heterogeneous entrainer. The simulation of the process is performed with the batch process simulator ProSimBatch [10]. It enables to evaluate operational parameters like the entrainer amount that are not provided by the feasibility and synthesis analysis The column model consists of usual plate by plate MESH (Material balance, Equilibrium, Summation of fractions and Heat balance) equations which are solved for the whole column, decanter included and taking into account the liquid-liquid phase split. Numerical treatment of the Differential Algebraic Equation (DAE) system and discrete events handling is performed with DISCo, a numerical package for hybrid systems with a DAE solver based on Gear's method. The column technical features and operating conditions are shown in Table 4. A sequence of two operational batch steps, namely an infinite reflux step and a distillate removal step described in detail below are simulated.

The amount of heterogeneous entrainer to be added to the original binary mixture is established from material balances considering that the entrainer composition in the still should be less than 0.01 at the end of infinite reflux batch step 1. Simulation shows that this composition is reached after two hours, which is the end of batch step 1. It should be noted that the initial ternary composition F is placed in distillation region I that is not limited by the vertex of ethyl acetate (see Figure 2). Therefore, in order to reach the vertex of ethyl acetate located in region II at the end of the distillation process, the still path must cross the unstable separatrix connecting the heteroazeotrope with the homoazeotrope nhexane -ethyl acetate. As established by Rodriguez -Donis et al. [START_REF] Rodríguez-Donis | Feasibility of Heterogeneous Batch Distillation[END_REF], such a crossing is possible with different reflux policies, either the reflux of entrainer -rich phase or a combination of both decanted liquid phases. In this example, the distillate removal step (batch step 2) is performed refluxing a combination of entrainer -rich L R and entrainerweak L W phases coming from the decanter, the reflux policy 2.3 of Rodriguez-Donis et al. [START_REF] Rodríguez-Donis | Feasibility of Heterogeneous Batch Distillation[END_REF]. As in Rodriguez-Donis et al. [START_REF] Rodríguez-Donis | Feasibility of Heterogeneous Batch Distillation[END_REF] the acetonitrile -rich phase L' R and n-hexane -rich phase L' W reflux flow rates are a fraction θ and ϕ of the amount of each phase (L R and L W ) contained in L°, the condensed vapour flow rate entering the decanter. Table 4 displays the θ and ϕ values that enable to withdraw the binary heterogeneous azeotrope during the whole distillation process. The process requires that the reflux of the acetonitrile -rich phase is larger (θ >1) than the acetonitrile quantity in the vapour flow rate leaving the column top L°, while the amount of n-hexane to be recycled at the top of the column is determined by the calculated ϕ value. θ > 1 means that drying of the acetonitrile -rich phase in the decanter during operation can occur and care should be taken to avoid it in the definition of an optimal reflux policy.

Figure 2 shows the simulation results, the still path and the column liquid profile, at the end of batch step 1 (infinite reflux) and at the beginning and end of batch step 2 (distillate removal). Table 5 displays significant process values. As can be seen in Figure 2 and Table 5, only 12.2% of acetonitrile is added to the initial azeotropic binary mixture and the resulting ternary mixture lies in the distillation region I. it has been observed [START_REF] Rodríguez-Donis | Feasibility of Heterogeneous Batch Distillation[END_REF] that heterogeneous batch distillation requires less entrainer than homogeneous batch distillation. The liquid profile in the column at the end of the infinite reflux operation links the still composition F R∞ to the binary heteroazeotrope n-hexane -acetonitrile in the first top tray of the column through the vapour line. Then, applying the reflux policy detailed in Table 4, distillate removal starts. Five minutes after the beginning of the distillation step, the column liquid profile changes suddenly and crosses the unstable separatrix due to the recycle of both decanted liquid phases. The reason is that the liquid phase leaving the first tray of the column top has a composition (x T ) inbetween the heteroazeotrope and the acetonitrile -rich phase composition and is located in a distillation region (III) different from region (I) which contains the still composition at infinite reflux F R∞ . Consequently, the simulated still path crosses the unstable separatrix separating distillation regions I and II (full black line in Figure 2). At the end of batch step 2, the still path reaches the ethyl acetate vertex when all n-hexane is drawn as distillate product. The operating parameters θ and ϕ are set to obtain the binary heteroazeotrope at the column top during the whole step 2. Therefore, the split ratio ω is kept constant. The n-hexane -rich phase is drawn as distillate product whereas the acetonitrile -rich phase can be recovered at the decanter at the end of the step 2. Table 5 reports the simulated purity and recovery yield of all components. Simulation computation time is typically less than 4 minutes on a Pentium III -based computer.

Experimental verification of the separation of the mixture n-hexaneethyl acetate by heteroazeotropic batch distillation in a bench scale rectification column.

The technical features of the column and the operating conditions are similar for the experiments and the simulations. The glass SCHOTT distillation column has a total height of 1.7 m and an internal diameter of 0.026 m. The total height is packed with stainless steel wire mesh rings of 3x3 mm. The total packed height is considered to be equivalent to 45 theoretical plates and the total liquid hold goes up of 52.5 ml. The condenser is subcooled at 25°C and total liquid reflux is provided at the top of the column by a solenoid valve. The open/close time of the solenoid valve is defined considering that 70% of the condensed vapour is returned to the column. The 30% remaining condensed vapour is collected into the decanter (70 ml.) where the light phase is drawn as distillate (n-hexane -rich phase) whereas a portion of the heavy entrainer -rich phase is recycled to the top of the column with a peristaltic pump.

Experiments taking into account these operating conditions (Table 4) are carried out in the bench distillation column and the results are presented in Figure 2 and Table 5.

Analysis of the still composition is performed every 0.5 hours. The analysis of purity of the light and heavy phases in the decanter and of the distillate product is made at the end of the process (4.8 h). The acetonitrile, n-hexane and ethyl acetate compositions are determined using a FISONS HRGC chromatograph with a FID detector and with a column HP INNOWax (30 m x 0.53 mm, 0.25µm film). The injector and detector temperature are held at 180°C while the column temperature is set equal to 65°C. Prior calibration is done using internal standards; the samples were prepared by weighing on a SARTORIUS BP 211D balance with a precision of 10 µg. The standard accuracy and reproducibility in the composition for all components are found to be ± 0.1% and 0.3% respectively.

As shown in Figure 2, the still path obtained experimentally (circles) is in excellent agreement with the still trajectory calculated by simulation. The selected reflux policy permits the still path to cross the separatrix into another distillation region than the initial feed region. Hence, the still path is able to reach the ethyl acetate vertex and this component remains pure into the still at the end of the distillate removal step. Such a behavior is not possible with a homogeneous entrainer that gives rise to a similar residue curve map because the distillation process is restricted to the distillation region where the initial composition of the mixture is located. In this case, ethyl acetate could not be obtained as an isolated product.

Table 5 displays the quantity and composition of the main streams of the process. Nhexane is drawn of as distillate product with a molar purity of 96.8%. Another binary batch distillation step is required to obtain n-hexane with a highest purity. Ethyl acetate with 99.6% of purity is recovered in the still at the end of the process, while the acetonitrilerich phase obtained from decanter can be reused in the next batch distillation process. The make-up of fresh entrainer for the next cycle of n-hexane -ethyl acetate separation is approximately 52 %.

Conclusions

Acetonitrile was chosen as a suitable heterogeneous entrainer for the separation of the azeotropic mixture n-hexane -ethyl acetate by heteroazeotropic batch distillation from the comparison of the residue curve maps of several candidate entrainers. We have confirmed the feasibility of the heteroazeotropic batch distillation via simulation and by experimentation in a bench batch distillation column. Good agreement was obtained between simulation and experimental results. The novel process with a single heterogeneous distillation column shows a great potential over a reported hybrid process with a homogeneous distillation column and a liquid -liquid extraction column. Indeed, the consumption of entrainer is low and the separation of all three components can be performed in a single batch distillation step using the defined reflux policy. The ethyl acetate purity is high enough to enable its recycling in a chromatography process whereas the n-hexane -rich phase must be purified in a simple non azeotropic batch distillation process.

FIGURE CAPTION

1. Residue curve maps of ternary systems n hexane -ethyl acetate - heterogeneous entrainer Figure 2. Simulation and experimental results of the separation of n hexane -ethyl acetate by using batch distillation with acetonitrile as a heavy heterogeneous entrainer Figure 1. Residue curve maps of ternary systems n hexane -ethyl acetate - heterogeneous entrainer

	n Hexane (68.7°C) [sn] Ethyl Acetate (77.1°C) [sn] Liquid -Liquid Equilibrium at 25°C		Ethyl Acetate (77.1°C) [sn] Vapour Line Experimental Still Path
	64.8°C [s] 76.1°C [s]	III	50.6°C [un] Liquid -Liquid -Vapour Equilibrium. Vapour Line Still composition (F R∞ ) at R ∞ Liquid -Liquid Still Path (Simulation) 76.1°C Equilibrium at 25°C Initial Still composition (F) Vapour Line [s] 64.8°C [s] 64.8°C [s] Liquid -Liquid -II Column liquid Profile at t = 5 min Vapour Equilibrium. Column liquid Profile at t = 3.8 h
	Ethyl Acetate (77.1°C) [sn] Acetonitrile	Methanol (64.7°C) [sn] I n Hexane Acetonitrile (81.5°C) [sn] Column liquid Profile at R ∞ n Hexane (68.7°C) [sn] 54.4°C [un] 54.4°C
	(81.5°C) [sn]	[un]	(68.7°C) [sn]
	Ethyl Acetate (77.1°C)			Ethyl Acetate (77.1°C)
	[sn]		[s] x T : Liquid composition in the top first tray
	Figure 2.		Liquid -Liquid Equilibrium at 25°C
	61.8°C [s]		Vapour Equilibrium. Liquid -Liquid -	Liquid -Liquid Equilibrium at 25°C
						Liquid -Liquid -
			Vapour Line		Vapour Equilibrium.
			64.8°C		64.8°C
			[s]			[s]
			61.2°C			Vapour Line
			[un]
	Water (100°C)	70.5°C	n Hexane	Nitromethane	n Hexane
	[sn]	[s]	(68.7°C) [sn]	(101.2°C) [sn]	(68.7°C) [sn]

Simulation and experimental results of the separation of n hexane -ethyl acetate by using batch distillation with acetonitrile as a heavy heterogeneous entrainer.

TABLE CAPTION Table 1 . Potential heterogeneous entrainers for the separation of n-hexane-ethyl acetate mixture with thermodynamics data from Gmehling et al. (1994).Table 2 . Binary interaction parameters for NRTL and UNIQUAC models.Table 3 . Purity of n-hexane -rich phase (x I ) and entrainer -rich phase (x II ) at 25°C.Table 4 . Operating conditions of heteroazeotropic batch distillation with acetonitrileTable 5 . Simulation and experimental results (x1: n-hexane; x2: ethyl acetate; x3: acetonitrile)

 CAPTION12345 The mole fraction indicated is that of the lightest component in the mixture ** Heteroazeotrope

		b.p	Binary azeotrope with	Binary azeotrope with	Ternary azeotrope
	Entrainer		n-hexane (68.7°C)	ethyl acetate (77.1°C)
		( ) C °	x	molar	*	( ) C T °	x	molar	*	( ) C T °	x	molar	*	( ) C T °
	Methanol	64.7	0.5100**	50.6	0.708	62.05		-		-
	Acetonitrile 81.5	0.5755**	54.4	0.6730	76.1		-		-
	Water	100.0 0.7860**	61.8	0.7000**	70.5	0.5764**	60.7
											0.2558
	Nitromethane 101.2 0.7271**	62.0			zeotropic		-		-

*

Table 1 . Potential heterogeneous entrainers for the separation of n-hexane-ethyl acetate mixture with thermodynamics data from Gmehling et al. (1994).

 1

	Binary coefficients for NRTL model (cal/mol)
	Binary parameters	A ij	A ji	α ij
	n-hexane -ethyl acetate	720.731	-33.844	0.1521
	n-hexane -acetonitrile	924.616	1212.64	0.2000
	ethyl acetate -acetonitrile	580.783	-253.247	0.1997
	Binary coefficients for UNIQUAC model (cal/mol)
	Binary parameters	A ij	A ji
	n-hexane -ethyl acetate	371.892	-133.73	-
	n-hexane -methanol	1451.225	7.831	-
	ethyl acetate -methanol	770.688	-147.777	-
	n-hexane -water	2292.190	-324.669	-
	ethyl acetate -water	768.058	214.271	-
	n-hexane -nitromethane	1080.270	7.335	-
	ethyl acetate -nitromethane	181.416	-42.973	-

Table 2 . Binary interaction parameters for NRTL and UNIQUAC models

 2

	20

Table 3 . Purity of n-hexane -rich phase (x I ) and entrainer -rich phase (x II ) at 25°C.Table 5 . Simulation and experimental results (x 1 : n-hexane; x 2 : ethyl acetate; x 3 : acetonitrile)

 35

		Simulation Results	Experimental Results
		quantity	molar	quantity	molar
		(mol)	fraction	(mol)	fraction
			x 1 0.705		x 1 0.705
	Initial mixture	12.2	x 2 0.173	12.2	x 2 0.173
			x 3 0.122		x 3 0.122
	n-hexane -rich phase (Distillate + Decanter)	8.6	x 1 0.951 x 2 0.002 x 3 0.047	8.4	x 1 0.968 x 2 0.001 x 3 0.031
	Acetonitrile -rich phase (Decanter)	1.1	x 1 0.107 x 2 0.157 x 3 0.736	0.99	x 1 0.082 x 2 0.190 x 3 0.728
	Final Still Content (ethyl acetate )	1.7	x 1 0.001 x 2 0.992 x 3 0.007	1.6	x 1 0.001 x 2 0.991 x 3 0.008
	Ethyl acetate recovery yield	81.04 %	75.1 %
	n-hexane recovery yield	95.1 %	94.5 %
	Acetonitrile recovery yield	54.5 %	48.4 %


Sylvie Bureau

C Scofield

Barbara Gouble

A Mcglone

H Boldingh

C Clark

Jean-Marc M Audergon

J Stanley

Evaluation of apricot fruit quality diversity in two countries, France and New Zealand

Keywords: Apricot, fruit, genetic diversity, Infrared Spectroscopy, quality traits, MIRS, NIRS, ethylene

à 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.


Imen Ghorbel
email: ghorbel.imen@gmail.com

Fayrouz Haddad
email: fayrouz.haddad@im2np.fr

Wenceslas Rahajandraibe
email: wence.rahajandraibe@im2np.fr

Mourad Loulou
email: mourad.loulou@ieee.org

Ultra-Low-Power Wideband NMOS LC-VCO Design for Autonomous Connected Objects

Keywords: LC-VCO, RF, design methodology, subthreshold, PVT, adaptive body biasing, low power, IoT. I

In this paper, a fully integrated subthreshold LC voltage controlled oscillator (VCO) is presented. A design methodology is also proposed to find the optimal parameters lowering the power consumption. This methodology has been applied to design oscillators for different frequency bands. Furthermore, an adaptive body biasing technique has been used to improve the startup constrains and allows a high immunity to PVT (process, voltage and temperature) variations. A VCO operating in 5 GHz ISM (Industrial, Scientific and Medical) band has been realized with the proposed methodology in 0.13μm CMOS. It consumes only 468 μW from 0.39V supply voltage. This makes it possible to meet the required specifications of autonomous connected objects and IoT applications. The measured oscillation frequency can be tuned from 5.14 GHz to 5.44 GHz. The obtained phase noise is approximately equal to -112 dBc/Hz in post-layout simulation (PLS) and -104.5 dBc/Hz in measure.

INTRODUCTION

The recent years have seen a drastic rise in Internet of Things (IoT) devices and connected products such as home automation systems, smart appliances, wireless sensors, and wearable etc. Objects are becoming intelligent and wirelessly interoperable to enhance transport, logistics, Health care digital, life and digital industry [START_REF] Lannacci | Internet of things (IoT); internet of everything (IoE); tactile internet; 5G -A (not so evanescent) unifying vision empowered by EHMEMS (energy harvesting MEMS) and RF-MEMS[END_REF][START_REF] Gazis | A survey of standards for machine-to-machine and the internet of things[END_REF]. It is envisioned that billions of objects will be equipped with sensors and will be connected to the Internet via heterogeneous access networks. This is due to the development of communication technologies such as radio frequency identification (RFID), semantic web services and wireless sensor networks etc. [START_REF] Razzaque | Middleware for Internet of Things: A Survey[END_REF]. Studies estimate that 212 billion connected IoT objects will be reach in the end of 2020 [START_REF] Al-Fuqaha | Internet of Things: A Survey on Enabling Technologies, Protocols and Applications[END_REF]. Key requirement for these devices is the implementation of low power, low-cost and reconfigurable systems. In fact, an IoT sensor node should be portable and able to work under batteries or even able to harvest ambient energy for several years. The general architecture of an IoT sensor node is shown in Fig. 1. It is composed of a wireless transceiver, a sensor, a power block and a processor unit, or a microcontroller (MCU) [START_REF] Char | Internet of things system design with integrated wireless MCUs[END_REF]. The wireless transceiver allows transmission and reception of the data using radio frequency (RF) connectivity. It is an important part of IoT systems which mostly affects the performance of the system. So, the transceiver RF must meet drastic power constraints.

The voltage controlled oscillator (VCO) is an essential building block in all radio transceivers. It generates a local signal that is used in both sides (reception and transmission) to transport the RF to the required band. This block consumes the significant portion of the power budget. In this work, the main challenge consists in minimizing its power consumption. Many researches have been reported in literature in order to reduce the power consumption of VCOs which can be applied at different levels of abstraction (system, circuit or technology levels). For the circuit level, two main techniques have been presented: reducing the supply voltage and using the current-reuse techniques. This later suffers from asymmetrical outputs in terms of amplitude and phase [START_REF] Yang | Low phase-noise and low power CMOS VCO constructed in currentreused configuration[END_REF][START_REF] Hsu | Design of low phase noise and low power modified current-reused VCOs for 10 GHz applications[END_REF][START_REF] Lin | A Low-power Low-voltage Current-reused Voltage controlled Oscillator with Amplitude-balanced[END_REF]. The first technique requires the use of subthreshold MOS devices which present some drawbacks such as high sensitivity to PVT variations [START_REF] Banchuin | Analysis and comprehensive analytical modeling of statistical variations in subthreshold MOSFET's high frequency characteristics[END_REF][START_REF] Ghorbel | Ultra Low Power RF Cross-Coupled VCO Design in the Subthreshold Regime with High Immunity to PVT Variations in 130nm CMOS Technology[END_REF]. These variations are caused by the exponential property of MOS transistors in the subthreshold region.

The purpose of this work is to design subthreshold VCOs with high immunity to PVT (process, voltage, and temperature) variations. Furthermore, a design methodology saving the power consumption and aiming optimal parameters has been used to satisfy the requirements of IoT applications: autonomy, stability and frequency agility. This paper is organized as follows. Section 2 investigates the subthreshold LC-VCO design. Section 3 describes the proposed design methodology of ultra-low power LC-VCOs. Section 4 presents the simulation results for different frequency bands. The obtained measurement results of an LC-VCO operating at 5.4 GHz are presented and discussed in section 5.

Finally, section 6 draws conclusions.

II.

SUBTHRESHOLD RF LC-VCO DESIGN

Several researches have been carried out in literature to implement ultra-low power analog circuits.

A reduction of the supply voltage is considered as one of the most important techniques to reduce the power consumption, which pushed designers to use special approaches, such as using subthreshold MOS transistors.

In the subthreshold regime, also known as weak inversion, the drain current depends exponentially on the gate source voltage Vgs. This latter is slightly below the transistor's threshold voltage Vth and the concentration of carriers is very low but not negligible. Subthreshold MOS transistors exhibit high transconductance gm for a given bias current and comparable noise performance with the superthreshold ones. However, the subthreshold operation is less robust to transistor mismatches, power-supply noise and temperature compared to the other regimes. It presents a high sensitivity to PVT variations. In order to guarantee stable operation even while working in worst conditions, designers should consider the variations of the transconductance, so that biasing techniques should be used. The adaptive body-biasing technique is an adequate solution, since it improves the immunity to PVT variations and ensures the stabilization of Vth and gm.

In this work, an NMOS cross-coupled topology using the adaptive body biasing technique has been selected. Figure 2 presents the detailed schematic of the VCO topology. It is composed by an LC resonator, also referred as LC-tank, an adaptive body biasing circuit (M3 and M4) and an NMOS cross-coupled pair (M1 and M2) biased in the subthreshold region. This pair is used to provide the negative conductance and to overcome the tank loss. The negative transconductance can be expressed as:

where gm_NMOS is the transconductance of the NMOS transistor. The adaptive body biasing circuit adjusts automatically the threshold voltage of MOS transistors as well as the transconductance. Moreover, this circuit improves the startup constraint and reduces the MOS drain source current. Thus, when the threshold voltage varies, the power consumption can be decreased.

The adaptive body biasing circuit allows to detect the VCO output amplitude (V+ and V-) and to feed it

Cvar Cvar Vdd Vctrl M 3 M 4 V b C 1 M 2 M 1 NMOS adaptive body biasing circuit 2 _ NMOS m active g g - = (1)
back to the body bias in the negative transconductance cell. If the output amplitude of the VCO varies, the bulk voltages of gm_NMOS are automatically adjusted in order to maintain a constant transconductance over PVT variations and to maintain the oscillation. As shown in Fig. 2, the adaptive body biasing circuit is composed by a symmetric NMOS transistor pair and a grounded capacitor C1. This circuit detects the lowest magnitude and stores the voltage Vb into the capacitor. When the output amplitude decreases, Vb increases which lowers the threshold voltage of the negative transconductance cell [START_REF] Park | Design techniques for a low-voltage VCO with wide tuning range and low sensitivity to environmental variations[END_REF]. This causes the increase of gm_NMOS as well as the output amplitude. Initially, the negative transconductance is not sufficient to satisfy the oscillation condition. Then, the adaptive bloc can be used to speed up the oscillation startup of the VCO. In fact, the gate and drain voltages of M3 and M4 are initially set to Vdd/2.

When these transistors are turned on, Vb increases and Vth of M1 and M2 decreases. This causes a higher value of the negative transconductance. In this case, the startup condition is met and the VCO starts oscillating.

III. LC-VCO DESIGN OPTIMIZATION APPROACH

In order to reach an ultra-low power operation suitable for IoT applications and to respect the various constraints, a design methodology using an optimization approach has been used in this work. It allows the design of high performance LC-VCO. This generic design methodology is based on the study and the optimization of VCO's components (inductor L, varactor C and active circuit). The designers should study the relation components/constraints and the impact of the physical and technological parameters of VCO's elements on its performances (power consumption, tuning range, phase noise and area occupation). For example, in order to enhance the power consumption and phase noise (PN) performances, it is necessary to optimize each element of the VCO and especially the LC resonator. Its quality factor (Q) and resistive losses (gp) must be taken into account, since increasing Q and reducing gp can improve both power consumption and PN. The design flow of high-performance LC-VCOs is illustrated in Fig. 3. It is based on three steps, as detailed below: i) the inductor optimization by maximizing the corresponding quality factor and minimizing the resistive losses, ii) followed by the varactor optimization, and iii) finally the active circuit optimization. In our study, the 130nm CMOS technology has been used to design VCOs.

▪ Inductance optimization strategy

The optimization process of the inductor requires a great attention when considering the design of highperformance VCO. The first step is to choose the inductor's type according to the targeted constraints.

In order to satisfy the need of designers, different types of inductors are proposed in the technology design-kit ('Ind_Sym_la', 'Ind_Sym_mf', 'Ind_dif_la', etc.). Furthermore, the used technology provides differential inductance that is also more appropriate to design differential structure [START_REF] Aspemyr | An ultra-Low Voltage, Low Power, Fully Integrated VCO for GPS in 90 nm RF-CMOS[END_REF]. Indeed, the differential inductor results in chip area saving compared to two single-ended inductors. In this work, several S-parameter simulations have been realized using Spectre-RF simulator of Cadence.

The obtained results have demonstrated that the 'Ind_dif_la' inductor type presents a good compromise between quality factor (QL), resistive loss (gL), occupation area and self-resonant frequency (SRF) compared to the other types. The second step is to size the inductance value. A high inductor value with minimal series resistance and maximal quality factor should be chosen since increasing L value reduces both power consumption and PN. Note that increasing L value is limited by its SRF, which should be higher than the oscillation frequency. Finally, designers should adjust both geometrical parameters: the inductor track width (W) and the number of turns (N). In this case, the influence of these two parameters has been studied. In order to decrease the resistive loss (gL), the minimum number of turns (Nmin) should be used. However, it is important to increase N in order to increase the inductance value that is also related to the technological constraints. Besides, the inductor's track width has a significant impact on SRF, QL and gL. Generally, increasing track width allows to reduce the gL and then to increase QL.

However, it presents some drawbacks such as increasing the substrate coupling. This can cause a decrease of self-resonance frequency. Note that this SRF should be higher than the operating frequency to maintain the inductive behavior. So, in order to increase the SRF, designers should use occasionally the minimum track width (Wmin).

▪ Varactor optimization strategy

The varactor is an essential component in LC resonator, since its quality factor affects the power dissipation and the phase noise performances. So, it is necessary to optimize the varactor to assess the performance of VCO. The optimization strategy begins with choosing the varactor type according to constraints of the targeted application. For example, for a high constraint in terms of tuning range, the MOS varactors are recommended. However, for high constraints in terms of consumption and PN, it is preferable to use the diode varactor, since it presents a good compromise between linearity, quality factor (Qv) and effective parallel equivalent conductance. For this design, diode varactor has been used.

Then, the value of C should be minimized in order to reduce PN and power consumption. However, this will contribute to lower the frequency tuning range. Finally, the physical parameters of the varactor have to be adjusted. For example, for a diode varactor, both number and width of P+ bars can be adjusted in order to obtain a good compromise between capacity tuning range and quality factor. These characteristics affect the performance of the oscillator in terms of tuning range and PN.

▪ Active circuit optimization strategy

The active circuit is an essential block of VCO that is used to overcome the tank loss. However, it also brings back parasitic elements which can degrade the phase noise of the VCO and modify the oscillation frequency. In this circuit, using low-leakage transistors with minimum channel length can improve the switching speed and decrease the power consumption. Then, transistor width is determined in order to provide the minimum transconductance. Using the maximum number of fingers (NF) decreases the gate resistance and consequently, improve the performances of VCO in terms of noise and power consumption. These charts (cf. Fig. 4) can be used to facilitate the optimal oscillator choice according to the constraints and requirements of the targeted application. Furthermore, it can be used to estimate the performances of any LC-VCO using the supply voltage reduction technique, in terms of power consumption and phase noise, for different frequency bands. second state combines worst PVT corners (high temperature, low Vdd and slow-slow process). The third state takes the best PVT corners (low temperature, high Vdd and fast-fast process). Table 1 presents the post-layout simulation results and the obtained percentages of variation in the different cases compared to the typical one. It can be noted that the oscillation is preserved even in the worst case. Also, an acceptable immunity to PVT variations is achieved confirming the efficiency of the adaptive body biasing technique. For example, in the best case, the observed variation in frequency, phase noise and power consumption is equal to only 1.2%, 1.7% and 29.3% respectively compared to the typical condition.

IV. SUB-mW RF LC-VCO IMPLEMENTATION AND RESULTS

1) Validation of the design methodology for different frequency bands

3) Measurement results

In order to assess the performances of the RF LC-VCO design, several measurements have been performed, under the following conditions: supply voltage of 0.39V, room temperature, noise and interferences all around the measuring station. This makes the measurement setup of our circuit very complex. The used on-chip measurement station is presented in Fig. 7. The DC and RF probing are positioned into the PADs on the wafer. The measured power consumption is very low and remains constant while changing the control voltage. It is only 468 µW from 0.39 V supply voltage. That confirms the effectiveness of the subthreshold operation and of the design methodology. It can be seen that the oscillation frequency can be tuned from 5.144 GHz to 5.458 GHz with a tuning voltage varying from 0.001 to 0.39V. then observed due to the experimental conditions in a very noisy environment around the test station. In fact, the noise requirements are stringent when measuring RF circuits. So, it is necessary to use shielded chamber (which is not our case). Additional VCO noise measurement will be processed in low-noise and calibrated environment, as a future work, to address the correct phase noise.

The measured performances of our VCO are summarized in Table 2 and compared to other similar lowpower VCOs operating around 5 GHz frequency band. It can be noted from Table 2 that our oscillator consumes the lowest power. Also, it achieves a comparable phase noise with the other works in post-layout simulations. A slight degradation of the phase noise in measure is observed due to the very noisy environment in which the measurements were carried out. Note that the phase noise needs to be measured in a low-noise environment soon. The implemented RF LC-VCO has good performances validating the design optimization methodology. applications. Oscillation frequencies vary between 5.14 and 5.44 GHz with compact core area. The post-layout simulation results show a phase noise of -112dBc/Hz @1MHz offset frequency, while the measured phase noise is -104.41 dBc/Hz @1MHz. Note that the phase noise should be done in lownoise environment to address the correct phase noise in the future work.

Fig 1 .

 1 Fig 1. IoT sensor node topology

Fig. 2 .

 2 Fig. 2. Schematic of the implemented LC-VCO

Fig. 3 .

 3 Fig. 3. Design methodology of sub-mW LC-VCO

  Rising the operating frequency presents important challenges in monolithic VCO design. Moreover, the oscillator must ensure high performances in terms of power consumption and PN. The contribution of our work is to design VCOs for different frequency bands and to verify the reliability of the proposed design methodology in high frequencies. Hence, several VCOs working at different frequency bands (from 1.5 to 9 GHz) have been implemented in 0.13µm CMOS technology. A low varactor value has been used (Cv=0.4 pF, Wp+ =1 µm, NP+=4), while, for each frequency band, we have varied both the transistors (W, NF) and the inductor (L, N, W).

Figure 4

 4 Figure4illustrates the obtained results for the different oscillators. Note that each point represents a fixed designed circuit. Once comparing the obtained performances versus frequency, it can be seen that when the frequency increases, the phase noise and power consumption are kept low. That confirms the effectiveness of both subthreshold operation and design methodology.

Fig. 4 .

 4 Fig.4. (a) Power consumption and (b) phase noise of LC-VCOs operating at different frequencies

Fig. 5 .

 5 Fig. 5. (a) Layout and (b) microphotograph of the LC-VCO

Fig. 6 .

 6 Fig. 6. The obtained phase noise histogram of 100 iterations

Fig. 7 .

 7 Fig. 7. The used on-chip measurement station

Fig. 8 .Figure 9

 89 Fig. 8. The output frequency spectrum measured for (a) Vctrl = 0.001V and (b) Vctrl= 0.39V

Fig. 9 .

 9 Fig. 9. The measured phase noise variation for Vctrl = 0.39 V

TABLE I .

 I LC-VCO PERFORMANCES UNDER PVT VARIATIONS Additional analysis of the process and mismatch effect, through Monte Carlo simulations, has been performed with Spectre-RF of Cadence tool taking into account the layout extracted parasitic. The obtained phase noise histogram of 100 iterations together with the fitted normal distribution curve are shown in Fig. 6. It depicts a mean phase noise value of -111.8 dBc/Hz and a standard deviation of only 0.19. These results confirm the robustness of the proposed sub-mW LC-VCO regarding the noise analysis.

TABLE II .

 II Performance comparison of the implemented LC-VCO with similar oscillators operating around 5 GHz This work describes the design and the implementation, in 0.13μm CMOS technology, of subthreshold sub-mW LC-VCOs operating at different frequency bands. These VCOs use the adaptive body biasing technique that allows a high immunity to PVT variations and can improve the startup constraint. of subthreshold VCOs operating in different frequency bands. A 5.4GHz LC-VCO oscillator has been fabricated and measured. This VCO achieves in measure ultra-low power consumption of only 468 μW drawn from 0.39V supply voltage. That makes the presented design and methodology suitable for IoT

		[15]	[16]	[17]	[18]	[19]	This work
	Process (nm)	180	180	180	180	180	130
	Supply (V)	1.1	0.65	1.8	1.2	1.5	0.39
	Frequency (GHz)	5.15	5.23	5.2	5.4	5.2	5.44
	Power (mW)	0.8	2.37	2.59	0.9	3.71	0.468
	Phase noise @1MHz (dBc/Hz)	-114.6	-115.1	-110	-112.7	-116	-104.4 / -112 *
	FOM (dBc/Hz)	-189.9	-185.7	-180	-188	-184	-182.4 / -191

* PLS

V. CONCLUSION

Furthermore, a design methodology, aiming the power consumption optimization in LC VCOs while satisfying design constraints, has been presented in this paper. It has been applied to size the components (2) *


Stéphane Balac

Simulation of rare-earth-doped high-power fiber-lasers using Matlab BVP solver

Keywords: optics, optical fiber-lasers, rare-earth-doped fiber, boundary value problem, collocation method 2000 MSC: 78A60, 78M99, 65L10, 65L60

Numerical simulation of light-wave propagation in double-clad rare-earth-doped fiber-lasers implies dealing with a two-point boundary value problem (BVP) with non separated boundary conditions. We show that this BVP can be solved in a simple way using Matlab BVP solver. However, this requires being able to provide to the BVP solver a relevant initial guess for the solution and the Jacobian matrix of the mapping defining the BVP. We show that when propagation losses and contribution to spontaneous emission are neglected, the BVP is equivalent to an initial value problem (IVP) whose solution, computed by Matlab IVP solver, provides a suitable guess for Matlab BVP solver. We also provide the expression of the Jacobian. This results in a very simple and cost-less Matlab program for the simulation of light-wave propagation in high-power fiber-lasers.

Introduction

In the theoretical study of high-power fiber-lasers, the most commonly used mathematical model is based on the rateequation theory [START_REF] Becker | Erbium-Doped Fiber Amplifiers: Fundamentals and Technology[END_REF][START_REF] Desurvire | Erbium-Doped Fiber Amplifiers: Principles and Applications[END_REF] that describes the distributions of inversion populations and optical powers along the fiber when a stationary regime has been attained. In this model, the power propagation equations consist of several coupled non-linear first order ordinary differential equations (ODE) depending on the position z along the fiber [START_REF] Kelson | Strongly pumped fiber lasers[END_REF]. The main difficulty for solving this system of non-linear ODE is related to the boundary conditions. We have a two-point boundary value problem (BVP), not an initial value problem (IVP), and the boundary conditions (BC) for this BVP are not separated, which means that we do not have independent relations at each endpoint for the unknowns but relations that correlate the values of the unknowns at the two endpoints, see relations (9c) and (9d) in the text below. This mathematical specificity is related to the reflection of the laser signal by Bragg mirrors at the two fiber extremities. Usually, existence and uniqueness for BVPs are much more difficult to investigate than for IVPs since there is no general theory [START_REF] Ascher | Numerical Solution of Boundary Value Problems for Ordinary Differential Equations[END_REF]. As well, it is usually more difficult to solve numerically a BVP than an IVP where the solution can be propagated from one endpoint of the computational domain to the other. For a BVP, the solution is not fully known at the domain endpoints so that the numerical approach used for IVPs does not apply.

They are mainly two approaches for solving a BVP. The first one is to use a shooting method [START_REF] Shampine | Solving ODEs with MAT-LAB[END_REF][START_REF] Ascher | Numerical Solution of Boundary Value Problems for Ordinary Differential Equations[END_REF]. Basically, shooting methods consist in reducing the solving of the BVP to the solving of a sequence of IVPs where the values for the missing boundary conditions at the fiber front end (z = 0) are enforced. After the missing boundary values were chosen arbitrarily at the Email address: stephane.balac@univ-rennes1.fr (Stéphane Balac) first attempt (the initial guess), they are adjusted for the next attempts by comparing the computed back end (z = L) boundary values to the BVP boundary conditions and using a root finding solver to compute the new boundary values at z = 0 for the next iteration. The second approach is to use collocation methods [START_REF] Ascher | Numerical Solution of Boundary Value Problems for Ordinary Differential Equations[END_REF][START_REF] Shampine | Solving ODEs with MAT-LAB[END_REF]. The collocation technique introduces a mesh of points to divide the fiber length into sub-intervals. Over each sub-interval, an approximation formula for the system of ODE is considered. It can be a Runge-Kutta formula or a Lobatto formula as used in Matlab BVP solvers [START_REF] Kierzenka | A BVP solver based on residual control and the Maltab PSE[END_REF]. This results in a global system of non-linear algebraic equations to which the relations given by the boundary conditions are added. Then, this system of non-linear equations is solved to obtain the approximate solution to the BVP.

In the last years, several authors have proposed numerical strategies to solve the BVP stemming from the modeling of light-wave propagation in high-power fiber-lasers. Among the more recent papers, one can quote reference [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF] where the standard shooting approach described above is improved by using some relations (close to relations (A.2) given in Appendix A) between the boundary values of the pump and laser signal powers at the both ends of the fiber. The main drawback of the method identified by the authors is that it is necessary for the initial guessed power laser to be lower than the true (unknown) value for the method to converge fast and provide accurate results. In [START_REF] Shang | A combined algorithm for simulating fiber lasers based on the shooting and relaxation methods[END_REF], the authors also use a shooting method but they combined it with a relaxation method to improve the convergence speed of the shooting method. We can also quote references [START_REF] Hu | Number sequence transition method based on Matlab BVP solvers for high power Yb3+ doped fiber lasers[END_REF][START_REF] Hu | A simple error control strategy using Matlab BVP solvers for Yb3+ doped fiber lasers[END_REF] where a complicated strategy is deployed to solve the BVP using Matlab BVP solver. For an incomprehensible reason, the authors use Matlab BVP solver in conjunction with an iterative method where, starting from guessed initial values and updating these values at each step, convergence to the so-lution to the BVP is expected. Although, according to their authors, these different numerical approaches provide the expected results in terms of simulation of light-wave propagation in high-power fiber-lasers, they appear to be excessively and unnecessarily complicated from a numerical point of view, especially the numerical approaches suggested in references [START_REF] Shang | A combined algorithm for simulating fiber lasers based on the shooting and relaxation methods[END_REF][START_REF] Hu | Number sequence transition method based on Matlab BVP solvers for high power Yb3+ doped fiber lasers[END_REF][START_REF] Hu | A simple error control strategy using Matlab BVP solvers for Yb3+ doped fiber lasers[END_REF]. Moreover, all these numerical approaches are very dependent on the initial guess for the iterative strategies they use and it is well known that fixed-point like methods efficiency is very dependent on the quality of the initial guess, see e.g. [13, chap.14] or [1, chap. 8].

In this paper, we show that it is possible to use Matlab BVP solver, in a very simple way, to solve the BVP stemming from the modeling of light-wave propagation in double-clad rareearth-doped fiber-lasers. Moreover, we show that an initial guess for Matlab BVP solver, close to the solution to the BVP, can be obtained at a very low cost. Indeed, when propagation losses and spontaneous emission are neglected, we show that an equivalent IVP to the BVP at hand can be obtained, for which the solution can be computed simply by using Matlab ode45 solver. This leads to a very simple and efficient Matlab program to simulate light-wave propagation in rare-earth-doped high-power fiber-lasers.

The document is organized as follows. In Section 2, we introduce the mathematical model for light-wave propagation in rare-earth-doped high-power fiber-lasers giving rise to the BVP with non-separated BC. Then, in Section 3, we explain how to use Matlab BVP solver in an efficient way to solve this BVP. Namely, we discuss the way an initial guess required by Matlab BVP solver can be computed so as to be close to the solution for a faster convergence of Matlab BVP solver and we give the expression of the Jacobian also required by Matlab BVP solver for a better efficiency. Section 4 is dedicated to an overview of the Shipol program written for the simulation of light-wave propagation in rare-earth-doped high-power fiberlasers. Finally, in Section 5, we compare the results provided by Shipol program to the existing literature results for validation purposes.

The mathematical model

We denote by L the length of the fiber and we assume that the axis of the fiber coincides with the z-axis of the reference frame. Using the rate-equation theory under the steady state condition [START_REF] Becker | Erbium-Doped Fiber Amplifiers: Fundamentals and Technology[END_REF][START_REF] Desurvire | Erbium-Doped Fiber Amplifiers: Principles and Applications[END_REF], we obtain a set of coupled equations for the pump and signal laser co-propagative and contra-propagative powers P + p , P - p , P + s , P - s describing light-wave propagation in a linear cavity as follows [START_REF] Kelson | Strongly pumped fiber lasers[END_REF]:

                               ± dP ± p dz (z) = (σ (e) p + σ (a) p ) N 2 (z) -σ (a) p N Γ p P ± p (z) -α p P ± p (z) ± dP ± s dz (z) = (σ (e) s + σ (a) s ) N 2 (z) -σ (a) s N Γ s P ± s (z) -α s P ± s (z) + Γ s σ (s) e P 0 N 2 (z) (1a) (1b)
where -P ± p and P ± s represent respectively the pump power (index p) at pump wavelength λ p and the signal laser power (index s) at laser wavelength λ s in the forward propagating direction (superscript +) and backward propagating direction (superscript -); -Γ p and Γ s are respectively the overlap integral factors between the ion-doping distribution and mode fields of pump and signal laser lights; -σ (a) p and σ (e) p denote respectively the absorption and emission cross-sections at the pump wavelength λ p ; -σ (a) s and σ (e) s denote respectively the absorption and emission cross-sections at the signal laser wavelength λ s ; -α p and α s represent the propagation losses coefficients including background loss and scattering loss at the pump and signal laser wavelengths respectively; -P 0 represents the contribution to spontaneous emission into the propagation laser mode; -N is the doping substance concentration density assumed to be constant along the fiber and N 2 is the upper-level population density.

For continuous-wave laser, the upper-level population density N 2 at position z (in unit m -3 ) is given by [START_REF] Becker | Erbium-Doped Fiber Amplifiers: Fundamentals and Technology[END_REF][START_REF] Kelson | Strongly pumped fiber lasers[END_REF]:

N 2 (z) = σ (a) s Γ s hν s P s (z) + σ (a) p Γ p hν p P p (z) A eff τ + (σ (a) s + σ (e) s ) Γ s hν s P s (z) + (σ (a) p + σ (e) p )
Γ p hν p P p (z)

N

(2) where P p (z) = P + p (z) + P - p (z) is the total pump power and P s (z) = P + s (z) + P - s (z) is the total signal power, τ is the spontaneous emission lifetime, A eff is the effective doping area, h = 6.62607015 10 -34 J • s is the Planck s constant and ν p = c/λ p , ν s = c/λ s , where c = 299 792 458 m • s -1 denotes light velocity in vacuum, are respectively the pump and signal frequencies. The contribution to spontaneous emission is given by P 0 = 2hc 2 δ λ /λ 3 s where δ λ is the spontaneous emission bandwidth.

We define the pump attenuation constant α (a) p and the signal attenuation constant α (a) s (in unit m -1 ) as

α (a) p = N Γ p σ (a) p (3a) α (a) s = N Γ s σ (a) s (3b)
and the saturation powers for the pump and signal (in unit W) as

P sat p = A eff hν p τ (σ (a) p + σ (e) p ) Γ p (4a) P sat s = A eff hν s τ (σ (a) s + σ (e) s ) Γ s (4b)
Using the above defined quantities, the ion population density N 2 given by ( 2) can be expressed as

N 2 (z) = τ A eff α (a) p hν p P p (z) + α (a)
s hν s P s (z) 1 + P s (z)

P sat s + P p (z) P sat p (5)
Substituting N 2 as given by ( 5) into the propagation equations (1), we obtain the following non-linear system of ODE

                           ± dP ± p dz (z) = -(α (a) p + α p ) P ± p (z) + G p (P ± p (z), P ± s (z)) P ± p (z) P sat p ± dP ± s dz (z) = -(α (a) s + α s ) P ± s (z)
+ G s (P ± p (z), P ± s (z))

P ± s (z) P sat s + P 0 P sat,e s (6a) (6b)
where (we recall that P p = P + p + P - p and P s = P + s + P - s )

G p (P ± p (z), P ± s (z)) = α (a) p P p (z) + ν p ν s α (a) s P s (z) 1 + P s (z) P sat s + P p (z) P sat p (7a) G s (P ± p (z), P ± s (z)) = α (a) s P s (z) + ν s ν p α (a) p P p (z) 1 + P s (z) P sat s + P p (z) P sat p (7b)
and where, with reference to (4b), we have set

P sat,e s = A eff hν s τ σ (e) s Γ s (8)
The boundary conditions (BC) are as follows [START_REF] Kelson | Strongly pumped fiber lasers[END_REF]:

                     P + p (0) = R (1) p P - p (0) + P + pump P - p (L) = R (2) p P + p (L) + P - pump P + s (0) = R (1) s P - s (0) P - s (L) = R (2) s P + s (L) (9a) (9b) (9c) (9d)
where P + pump and P - pump are the pump powers available at the front (z = 0) and back (z = L) ends respectively, R (1) s and R (2) s are respectively the front and back ends mirror reflectivities at the signal wavelength λ s and R (1) p and R (2) p are the reflectivities at the pump wavelength λ p . Thus, we have to solve a BVP that consists in the non-linear first order coupled ODE [START_REF] Kelson | Strongly pumped fiber lasers[END_REF] and linear two-point BC [START_REF] Peysokhan | Analytical formulation of a highpower Yb-doped double-cladding fiber laser[END_REF]. The specificity of these BC is that they are not separated due to the reflection at the fiber ends that mixes forward and backward propagating powers. However, this BVP can be solved in a very simple way under Matlab using its BVP solvers contrary to what has been stated in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF][START_REF] Shang | A combined algorithm for simulating fiber lasers based on the shooting and relaxation methods[END_REF][START_REF] Hu | Number sequence transition method based on Matlab BVP solvers for high power Yb3+ doped fiber lasers[END_REF][START_REF] Hu | A simple error control strategy using Matlab BVP solvers for Yb3+ doped fiber lasers[END_REF] and without need of complicated numerical approaches. To this end, we express BVP ( 6)-( 9) in the form

Y ′ (z) = F(Y(z)) ∀z ∈]0, L[ ( 10
)
where the unknown vector Y(z) ∈ R 4 is defined as

Y(z) = P + p (z) P - p (z) P + s (z) P - s (z) ⊤ ∈ R 4
where ⊤ indicates matrix transposition and the mapping F is defined as

F : Y ∈ R 4 -→                      -(α (a) p + α p ) Y 1 + G p (Y) Y 1 P sat p (α (a) p + α p ) Y 2 -G p (Y) Y 2 P sat p -(α (a) s + α s ) Y 3 + G s (Y) Y 3 P sat s + P 0 P sat,e s (α (a) s + α s ) Y 4 -G s (Y) Y 4 P sat s + P 0 P sat,e s                      (11)
where here and throughout the paper Y 1 , . . . , Y 4 denotes the components of Y. The BC (9) can be expressed in matrix form as

M 0 Y(0) + M L Y(L) =                P + pump P - pump 0 0                (12)
where the two matrices M 0 and M L in M 4 (R) are given by

M 0 =                1 -R (1) p 0 0 0 0 0 0 0 0 1 -R (1) s 0 0 0 0                M L =                0 0 0 0 -R (2) p 1 0 0 0 0 0 0 0 0 -R (2) s 1               

MATLAB BVP solvers

Matlab BVP solver bvp4c implements a collocation method based on the three-stage Lobatto IIIa formula, see [START_REF] Kierzenka | A BVP solver based on residual control and the Maltab PSE[END_REF]. Matlab BVP solver bvp5c implements the four-stage Lobatto IIIa formula and slightly differs from bvp4c on some internal implementation choices, see Matlab documentation for details. Matlab BVP solvers efficiency is improved when an initial guess close enough to the solution is provided and when the Jacobian matrix of the function F defining the ODE given by ( 11) is available. We examine now these two topics for BVP ( 6)- [START_REF] Peysokhan | Analytical formulation of a highpower Yb-doped double-cladding fiber laser[END_REF].

Propagation losses in the fiber, although not negligible, are generally small and they contribute lowly to the BVP solution that is mainly driven by other optical phenomena. The same can be observed concerning the contribution to spontaneous emission into the propagation laser mode. This observation leads us to chose as initial guess for Matlab BVP solver, the solution to BVP (6)-( 9) in the case when α p = α s = 0 and P 0 = 0. This idea would be idiotic if solving the BVP under these assumptions was as difficult as solving the complete BVP, but such is not the case. Indeed, we show in Appendix A that when losses and spontaneous emission are neglected, BVP ( 6)-( 9) can be reformulated in an equivalent Initial Value Problem (IVP) that can be solved easily using standard numerical method such as the Runge-Kutta RK4 method as e.g. implemented in Matlab ode45 solver [START_REF] Shampine | The Matlab ODE suite[END_REF]. This solution, though different from the real solution to the BVP is likely to be close to it and it will thus provide a very good initial guess for Matlab BVP solvers. Therefore, using Matlab ode45 solver, we solve the IVP composed of ODE ( 6) and initial conditions (A.9) for P - p (0), (A.10) for P - s (0), (9a) for P + p (0) and (9c) for P - p (0). Then, we provide to Matlab BVP solver bvp4c, as initial guess, the solution Y to this IVP computed at the grid nodes of a subdivision (z j ) j=0,...,J of the interval [0, L].

Let us now consider the calculation of the Jacobian matrix J F of the mapping F introduced in [START_REF] Shampine | The Matlab ODE suite[END_REF]. For all Y ∈ R 4 , we have

J F (Y) =               ∂ 1 F 1 (Y) ∂ 2 F 1 (Y) ∂ 3 F 1 (Y) ∂ 4 F 1 (Y) ∂ 1 F 2 (Y) ∂ 2 F 2 (Y) ∂ 3 F 2 (Y) ∂ 4 F 2 (Y) ∂ 1 F 3 (Y) ∂ 2 F 3 (Y) ∂ 3 F 3 (Y) ∂ 4 F 3 (Y) ∂ 1 F 4 (Y) ∂ 2 F 4 (Y) ∂ 3 F 4 (Y) ∂ 4 F 4 (Y)               (13)
where ∂ j refers to the partial derivation with respect to the j-th variable and F i refers to the i-th component of the mapping F. A straightforward calculation shows that the partial derivatives of F are given for j ∈ {1, . . . , 4} and for Y ∈ R 4 by

∂ j F 1 (Y) = -(α (a) p + α p ) δ 1, j + G p (Y) P sat p δ 1, j + Y 1 P sat p ∂ j G p (Y) ∂ j F 2 (Y) = (α (a) p + α p ) δ 2, j - G p (Y) P sat p δ 2, j - Y 2 P sat p ∂ j G p (Y) ∂ j F 3 (Y) = -(α (a) s + α s ) δ 3, j + G s (Y) P sat s δ 3, j + Y 3 P sat s + P 0 P sat,e s ∂ j G s (Y) ∂ j F 4 (Y) = (α (a) s + α s ) δ 4, j - G s (Y) P sat s δ 4, j - Y 4 P sat s + P 0 P sat,e s ∂ j G s (Y)
where δ i, j denotes the Kronecker symbol (equal to 1 when i = j and 0 otherwise). The partial derivatives of the mappings G p and G s defined in [START_REF] Kierzenka | A BVP solver based on residual control and the Maltab PSE[END_REF] are given by

∂ 1 G p (Y) = ∂ 2 G p (Y) = α (a) p - G p (Y) P sat p 1 ψ(Y) ∂ 3 G p (Y) = ∂ 4 G p (Y) = ν p ν s α (a) s - G p (Y) P sat s 1 ψ(Y)
and

∂ 1 G s (Y) = ∂ 2 G s (Y) = ν s ν p α (a) p - G s (Y) P sat p 1 ψ(Y) ∂ 3 G s (Y) = ∂ 4 G s (Y) = α (a) s - G s (Y) P sat s 1 ψ(Y) where ψ(Y) = 1 + Y 1 + Y 2 P sat p + Y 3 + Y 4 P sat s .

The Shipol program

A Matlab program based on the ideas presented above has been implemented.

It can be downloaded with it documentation from HAL library at the following URL : https://hal.archives-ouvertes.fr/hal-03212156

The program basically consists in the following stages :

1. Computation of a guess (rough approximate solution) by solving ODE (6) with α s = α p = P 0 = 0 under the initial conditions given by (A.9) for P - p (0), (A.10) for P - s (0), (9a) for P + p (0) and (9c) for P - p (0) using Matlab ode45 solver.

2. Solving of the BVP [START_REF] Shampine | Solving ODEs with MAT-LAB[END_REF] with BC ( 12) by Matlab bvp4c solver with the guess computed at stage 1. Efficiency of bvp4c is improved by supplying the Jacobian matrix of the mapping F to bvp4c.

3. Results display. Values of the pump and signal powers in the forward and backward directions (P + p , P - p , P + s , P - s ) at the two fiber ends are shown as well as the value of the total laser signal output power. Graphical display of the pump and signal powers in the forward and backward directions as a function of the position z along the fiber together with the ratio N 2 /N are also provided.

The accuracy of the computed solution can be evaluated using an error estimator based on relations (A.1). We must have P + p (0) P - p (0) = P + p (L) P - p (L) and, when P 0 = 0, we must also have P + s (0) P - s (0) = P + s (L) P - s (L). Therefore, we use as error estimator

E P = P + p (0) P - p (0) -P + p (L) P - p (L) P + p (0) P - p (0) + P + p (L) P - p (L) (14a)
E s = P + s (0) P - s (0) -P + s (L) P - s (L) P + s (0) P - s (0) + P + s (L) P - s (L) (14b)
where E s is relevant only when P 0 ≡ 0. The closer to zero these estimators are, the more accurate the computation can be expected to be.

Numerical validation

We have compared the simulation results provided by the Shipol program to several results available in the literature. (The numerical methods used in these papers are not always detailed and the codes are not available for further comparisons.) The numerical results presented below were obtained using the software Matlab (R2018b) under Linux-Ubuntu OS on an Intel Core i5 computer with 8 Go RAM. We have considered the Yb 3+ -doped fiber-laser investigated in [START_REF] Kelson | Strongly pumped fiber lasers[END_REF], the parameters of which are given in Table 1. The ASE is neglected in [START_REF] Kelson | Strongly pumped fiber lasers[END_REF] so that we set P 0 ≡ 0 in (1b). We have depicted in Fig. 1 the pump and signal powers in the forward and backward directions, i.e. P + p , P - p , P + s , P - s , as a function of the position along the fiber together with the ratio N 2 (z)/N. This figure is identical to [START_REF] Kelson | Strongly pumped fiber lasers[END_REF]Fig. 8]. We have obtained a forward signal power at z = 0 with value P + s (0) = 2.02871 W and a backward signal power at z = 0 with value P - s (0) = 2.07011 W. At the fiber back end z = L, the values are P + s (L) = 10.2465 W and P - s (L) = 0.409861 W. The forward pump power at z = L is 4.15728 W whereas the backward pump power remains zero along the fiber. The CPU time for the simulation was 1.8 s and the error estimator defined in (14b) was E s = 2.7418 10 -12 . (Note that the estimator E p not relevant here since P - p ≡ 0.) Taking into account the ASE in the simulation, we have obtained P + s (0) = 2.02874 W, P - s (0) = 2.07014 W, P + s (L) = 10.2466 W and P - s (L) = 0.409862 W. These values are not significantly different from the one obtained by neglecting the ASE which justify the assumption made in [START_REF] Kelson | Strongly pumped fiber lasers[END_REF]. We want to point out that the numerical method used in [START_REF] Kelson | Strongly pumped fiber lasers[END_REF] to solve the BVP is not specified as it is mentioned only for comparison purposes with a semi-analytical approach valid under strong pumping conditions, such that the signal power is high enough to saturate the gain medium and suppress the spontaneous emission. Pump and signal powers in the forward and backward directions P + p , P - p , P + s , P - s as a function of the position along the fiber for the Yb 3+ fiberlaser the parameters of which are given in Table 1 (left axis). Ratio N 2 (z)/N as a function of the position along the fiber in dotted line and right axis.

Simulation of light-wave propagation in a

As a complement, we have made a simulation for the Yb 3+ doped fiber-laser investigated in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF]. The parameters of the fiber are similar to the one given in Table 1 excepted that the fiber length is L = 150 m, the pump power launched into the front end P + pump is 1000 W and the signal background losses are α s = 3 10 -3 m -1 . We have depicted in Fig. 2 the pump and signal powers in the forward and backward directions P + p , P - p , P + s , P - s as a function of the position along the fiber together with the ratio N 2 (z)/N. We have obtained a forward signal power at z = 0 with value P + s (0) = 105.493 W and a backward signal power at z = 0 with value P - s (0) = 107.645 W. At the fiber back end z = L, the values are P + s (L) = 532.818 W and P - s (L) = 21.3127 W. The forward pump power at z = L is 8.84451 W whereas the backward pump power remains zero along the fiber. The CPU time for the simulation was 0.7 s and the error estimator defined in (14b) was E s = 1.6709 10 -11 . Fig. 2 is to be compared with [8, Fig. 2]. The two figures are significantly different. Since the only difference with the previous simulation that matches the results published in [START_REF] Kelson | Strongly pumped fiber lasers[END_REF] is the fiber length and pump power launched at the fiber front end, we infer that Fig. 2 in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF] do not correspond to the fiber parameter values provided in the text in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF] or that the numerical method presented in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF] is implemented in an erroneous way or that the method itself is erroneous. It should be noted that we have detected an other error in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF] in the results reported for the thulium-doped fiber-laser, see Section 5.3.

Simulation of light-wave propagation in a Nd 3+ doped fiber-laser

We then have compared the results obtained by our numerical method to the one given in [START_REF] Kelson | Strongly pumped fiber lasers[END_REF] for a Nd 3+ -doped high-power double-clad fiber laser characterized by the values given in Table 2. We have depicted in Fig. 3 the pump and signal powers in the forward and backward directions, i.e. P + p , P - p , P + s , P - s , as a function of the position along the fiber together with the ratio N 2 (z)/N. This figure is identical to the one depicted in [6, Fig. 5]. We have obtained a forward signal power at z = 0 with value P + s (0) = 2.50126 W and a backward signal power at z = 0 with value P - s (0) = 2.5523 W. At the fiber back end z = L, the values are P + s (L) = 12.6333 W and P - s (L) = 0.505331 W. The forward pump power at z = L is 2.19664 W whereas the backward pump power remains zero along the fiber. The CPU time for the simulation was 0.6 s and the error estimator defined in (14b) was E s = 2.6426 10 -12 . (Note that the estimator E p is not relevant here since P - p ≡ 0.) For comprehensiveness, we have also made a simulation when the pump is injected at the fiber back end z = L, i.e. P + p (0) = 0 and P - p (L) = 20W, the other parameters are as given in Table 2. The result is depicted in Fig. 4 and it is identical to the one depicted in [START_REF] Kelson | Strongly pumped fiber lasers[END_REF]Fig. 6 2 except that the pump is injected at z = L (left axis). Ratio N 2 (z)/N as a function of the position along the fiber in dotted line and right axis.

Simulation of light-wave propagation in a thulium-doped

fiber-laser Finally, we have considered the Tm 3+ -doped high-power double-clad fiber-laser proposed in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF] characterized by the values given in Table 3. The ASE is neglected in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF] so that we set P 0 ≡ 0 in the simulation. Note that the doping substance concentration given in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF] was N = 4.68 10 26 ions/m 3 and it is changed for N = 8.6 10 25 ions/m 3 to comply with the value of the literature quoted in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF]. We have depicted in Fig. 5 the pump and signal powers in the forward and backward directions (P + p , P - p , P + s , P - s ) as a function of the position along the fiber together with the ratio N 2 (z)/N. This figure is quite similar to the one depicted in [8, Fig. 5]. Note however that when using the doping substance concentration given in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF] (N = 4.68 10 26 ions/m 3 ) rather than the one given in Table 3, the figure we obtain is significantly different from [8, Fig. 5]. We have obtained a forward signal power at z = 0 with value P + s (0) = 74.8308 W and a backward signal power at z = 0 with value P - s (0) = 76.358 W. At the fiber back end z = L, the values are P + s (L) = 377.953 W and P - s (L) = 15.1181 W. The forward pump power at z = L is 2.02216 W whereas the backward pump power remains zero along the fiber. The CPU time for the simulation was 1.3 s and the error estimator defined in (14b) was E s = 7.5682 10 -12 . Note that when the ASE was taken into account in the simulation, we obtained similar results to the one provided above which justify the assumption made in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF] to neglect the ASE.

Parameter

Conclusion

We have shown that the boundary value problem (BVP) stemming from the modeling of light-wave propagation in highpower fiber-lasers can be solved in a simple way by Matlab BVP solver. This approach renders meaningless the numerical methods proposed in [START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF][START_REF] Shang | A combined algorithm for simulating fiber lasers based on the shooting and relaxation methods[END_REF][START_REF] Hu | Number sequence transition method based on Matlab BVP solvers for high power Yb3+ doped fiber lasers[END_REF][START_REF] Hu | A simple error control strategy using Matlab BVP solvers for Yb3+ doped fiber lasers[END_REF] to solve this BVP. The effectiveness of the method relies on the efficiency of Matlab BVP solver program and on the practicability to provide an initial guess close to the solution to the BVP at a very low cost as shown in the paper. Moreover we have provided the expression of the Jacobian of the mapping defining the BVP that further improves the overall effectiveness of the simulation. This results in a very simple and efficient Matlab program for the simulation of light-wave propagation in high-power fiberlasers. We hope that the author Matlab program distributed under CeCILL-C free software license, will benefit to the greatest number of researchers in the optics community. The same calculation can be done with the three other equations in (A.5a)-(A.5b). Taking into account relations (A.2), it remains the following two equations for the unknowns P - p (0) and P - s (0): where we have set η = R (1) s R (2) s -R (1) s -R (1) s R (2) s + 1 log(R (1) s R (2) s )α a s L + α a p L (A.8)

±

                                 log         P + pump + R (
Φ(P - p (0)) = (1 -R (
From (A.7), we deduce that P - p (0) is solution to the following algebraic equation with unknown X 1 4 P - pump 2 + R (1) p R (2) p X 2 + P - pump R (2) p Xe -A R (1) p R (2) p = e -A R (2) p P + pump + 1 2 P - pump .

Figure 1 :

 1 Figure1: Pump and signal powers in the forward and backward directions P + p , P - p , P + s , P - s as a function of the position along the fiber for the Yb 3+ fiberlaser the parameters of which are given in Table1(left axis). Ratio N 2 (z)/N as a function of the position along the fiber in dotted line and right axis.

Figure 2 :

 2 Figure 2: Pump and signal powers in the forward and backward directions P +p , P - p , P + s , P - s as a function of the position along the fiber for the Yb 3+ doped fiber-laser considered in[START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF].

Figure 3 :

 3 Figure 3: Pump and signal powers in the forward and backward directions P +p , P - p , P + s , P - s as a function of the position along the fiber for the Nd 3+ fiberlaser the parameters of which are given in Table2(left axis). Ratio N 2 (z)/N as a function of the position along the fiber in dotted line and right axis.

Figure 4 :

 4 Figure 4: Pump and signal powers in the forward and backward directions P +p , P - p , P + s , P - s as a function of the position along the fiber for the Nd 3+ fiberlaser the parameters of which are given in Table2except that the pump is injected at z = L (left axis). Ratio N 2 (z)/N as a function of the position along the fiber in dotted line and right axis.

Figure 5 :

 5 Figure 5: Pump and signal powers in the forward and backward directions P +p , P - p , P + s , P - s as a function of the position along the fiber for the Tm 3+ doped fiber-laser.

  2) p )φ(P - p (0)) + (1 -R (1) p )P - p (0) -(P + pump + P - pump )

Table 1 :

 1 Parameter values for the Yb 3+ doped fiber-laser investigated in[START_REF] Kelson | Strongly pumped fiber lasers[END_REF].

	Yb 3+ doped

Table 2 :

 2 ]. Parameter values for the Nd 3+ doped fiber-laser investigated in[START_REF] Kelson | Strongly pumped fiber lasers[END_REF].

	Parameter	Notation	Value
	Fiber length	L	50 m
	Pump wavelength	λ p	0.808 10 -6 m
	Signal wavelength	λ s	1.060 10 -6 m
	Pump power at front end Pump power at back end	P + pump P -pump	20 W 0 W
	Pump overlap factor	Γ p	0.01
	Signal overlap factor	Γ s	0.8
	Doping substance concentration	N t	2 10 24 ions/m 3
	Pump absorption cross section	σ (a) p	2 10 -24 m 2
	Pump emission cross section	σ (e) p	0 m 2
	Signal absorption cross section	σ (a) s	0 m 2
	Signal emission cross section	σ (e) s	2.5 10 -24 m 2
	Pump background losses	α p	4.5 10 -3 m -1
	Signal background losses	α s	0 m -1
	Front mirror reflectivity at λ p	R (1) p	0
	Output mirror reflectivity at λ p	R (2) p	0
	Front mirror reflectivity at λ s	R (1) s	0.98
	Output mirror reflectivity at λ s	R (2) s	0.04
	Lifetime of the upper-level atoms	τ	4 10 -4 s
	Effective mode area	A eff	1 10 -11 m

Table 3 :

 3 Parameter values for the Tm 3+ doped fiber-laser investigated in[START_REF] Liu | An improved shooting algorithm and its application to high-power fiber lasers[END_REF].

		Notation	Value
	Fiber length	L	25 m
	Pump wavelength	λ p	0.790 10 -6 m
	Signal wavelength	λ s	1.973 10 -6 m
	Pump power at front end Pump power at back end	P + pump P -pump	1000 W 0 W
	Pump overlap factor	Γ p	0.01
	Signal overlap factor	Γ s	0.752
	Doping substance concentration	N t	8.6 10 25 ions/m 3
	Pump absorption cross section	σ (a) p	5 10 -25 m 2
	Pump emission cross section	σ (e) p	0 m 2
	Signal absorption cross section	σ (a) s	1 10 -26 m 2
	Signal emission cross section	σ (e) s	2.5 10 -25 m 2
	Pump background losses	α p	3 10 -3 m -1
	Signal background losses	α p	5 10 -3 m -1
	Front mirror reflectivity at λ p	R (1) p	0
	Output mirror reflectivity at λ p	R (2) p	0
	Front mirror reflectivity at λ s	R (1) s	0.98
	Output mirror reflectivity at λ s	R (2) s	0.04
	Lifetime of the upper-level atoms	τ	1.45 10 -5 s
	Effective mode area	A eff	1.39 10 -11 m

  Let us consider equation (A.5a) for P + p . We divide both sides by P + p (z) under the assumption that the forward pump power P + p never cancels and we integrate over the interval [0, L]:

													0	L	1 P + p (z)	dP + p (z) dz	dz = -	1 P sat p	0	L	dP + p dz	(z) -	dP -p dz	(z) dz
													-	ν p ν s P sat p	0	L	dP + s dz	(z) -	dP -s dz	(z) dz -α a p L.
													It follows that
													log	       	P + p (L) pump + R (1) P + p P -p (0)	        = -	1 p P sat	P + p (L) -P + pump -R (1) p P -p (0)
													-P -pump -R (2) p P + p (L) + P -p (0) -	ν p p ν s P sat	P + s (L) -R (1) s P -s (0)
													-R (2) s P + s (L) + P -s (0) -α a p L.
		dP ± p dz	(z) = -	dP ± p dz	(z) -	dP -p dz	(z)	P ± p (z) P sat p
			-	ν p ν s	dP ± s dz	(z) -	dP -s dz	(z)	P ± p (z) p P sat	-α a p P ± p (z)	(A.5a)
	±	dP ± s dz	(z) = -	ν s ν p	dP ± p dz	(z) -	dP -p dz	(z)	P ± s (z) P sat s
			-	dP ± s dz	(z) -	dP -s dz	(z)	P ± s (z) s P sat	-α a s P ± s (z)	(A.5b)

Acknowledgments

The author would like to thanks Prof. Thierry Chartier, from Foton Laboratory, ENSSAT-University of Rennes 1 (France), for having drawn attention to this numerical problem and for fruitful discussions on fiber-lasers.

Appendix A. An equivalent initial value problem in the absence of losses

We consider in this appendix the special case when the losses and the contribution to spontaneous emission into the propagation laser mode can be neglected. That is to say, we consider BVP ( 6)- [START_REF] Peysokhan | Analytical formulation of a highpower Yb-doped double-cladding fiber laser[END_REF] under the assumptions that α s = α p = 0 and P 0 = 0. In this special case, the solutions P ± p , P ± s to the set of ODE (6) exhibit the following features: There exist two real numbers C p and C s such that for all z ∈ [0, L]

This property can be easily proved by considering the mapping z ∈ [0, L] → P + p (z) P - p (z) and showing that its derivative is zero. We have

We then sum equation (6a) for P + p preparatorily multiplied by P - p and equation (6a) for P - p preparatorily multiplied by P + p to show that the derivative is zero. Note that the properties (A.1) are true even when the losses are not zero. However, when P 0 0, only the relation on P ± p remains valid. We deduce from (A.1) and from BC (9a) and (9b), the following relationships between the boundary values at the fiber ends:

p (A.2a)

Note that when R (2) p = 0, (A.2a) has to be changed for

Let us now show how we can express the two quantities P - p (0) and P - s (0) in terms of the fiber-laser characteristic parameters. By difference between (1a) considered for P + p and P - p on the one hand and by difference between (1b) considered for P + s and P - s the other hand, one can easily show that N 2 given by ( 2) can be expressed as

Substituting this expression of N 2 into equations (1a)-(1b), we obtain

Solving this algebraic equation, we obtain that its non-negative solution reads

p + P - pump e A e 2A -R (1) p R (2) p (A.9)

Then, from (A.6b) , we have (1) s R (2) s ) -

Finally, under the assumption that α s = α p = 0 and P 0 = 0, the BVP ( 6)-( 9) is equivalent to the IVP composed of the ODE (6) and initial conditions (A.9) for P - p (0), (A.10) for P - s (0), (9a) for P + p (0) and (9c) for P - p (0). Note that in [START_REF] Kelson | Strongly pumped fiber lasers[END_REF] the authors also propose an IVP, different of the one proposed here, equivalent to the BVP ( 6)-( 9) under the following additional assumptions: σ (a) s ≪ σ (a) p , σ (e) p = 0 and R (1) p = R (2) p = 0. We did not need these additional assumptions here. As well, in [START_REF] Peysokhan | Analytical formulation of a highpower Yb-doped double-cladding fiber laser[END_REF] the authors proposed an analytical solution to the BVP in the absence of losses but also require additional assumptions.