Cristian Camilo

Ruiz Vasquez

Noureddine Lebaz

Isabelle Ramière

Sophie Lalleman

Denis Mangin

Murielle Bertrand

Cristian Camilo Ruiz Vasquez

Fixed point convergence and

Keywords: population balance, steady state, oxalic precipitation, fixed point solution, convergence acceleration

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

Introduction

* Corresponding author: camilo.ruiz@etu.univ-lyon1.fr In precipitation processes, the crystal size distribution (CSD) is a decisive product property.

Continuously operated well mixed precipitators present in principle poor performances in this concern [START_REF] Yazdanpanah | The Handbook of Continuous Crystallization[END_REF]. For this reason, the simulation of continuous crystallizers has been the target of many modelling efforts during the last decades and remains a scientific challenge [START_REF] Myerson | Handbook of Industrial Crystallization[END_REF]. Crystallization phenomena are characterized by a strong dependence to the mixing issues: local concentrations involve local crystallization kinetics and macromixing conditions define the transport of the continuous and dispersed phases. This research work appears as a necessary step to further understanding the CSD evolution in crystallizers exhibiting complex hydrodynamics.

The evolution of the CSD in a control volume is tracked by the Population Balance Equation (PBE).

The general statement at the steady state as first proposed in 1964 by Hulburt and Katz is given as [START_REF] Hulburt | Some problems in particle technology[END_REF] : where is the crystal size (m), is the number-based density function (m -4

), is the flowrate (m

3 .s -1
), is the total volume (disperse and continuous phases) (m 3 ), is the total streams, is the crystals growth rate (m.s -1

) and and are the birth and death terms respectively, related to the appearance and disappearance of crystals in the system (through nucleation, aggregation and breakage mechanisms). Hereafter, the CSD dependence over the size will be obviated in order to simplify equations formulation ( ).

The PBE is a well-known mathematical framework used to describe dispersed systems evolution. In precipitation, it is employed to determine CSD from crystallization kinetics (nucleation, growth, breakage and agglomeration) and transport phenomena. It needs to be coupled with the mass, momentum and energy balances to estimate the CSD and the supersaturation. Since several strongly coupled phenomena influence the CSD (hydrodynamics, solid dispersed particles, reaction kinetics, crystallization mechanisms, etc.), a highly nonlinear algebraic partial differential equations system needs to be solved in order to model the precipitation process. Because of the complexity of the mathematical problem, analytical solutions exist only in few cases under highly restrictive hypothesis [START_REF] Mersmann | Crystallization Technology Handbook[END_REF]. Indeed, the numerical solution of the PBE is still an active research field [START_REF] Ramkrishna | Population Balance Modeling: Current Status and Future Prospects[END_REF]- [START_REF] Wang | Recent progress of continuous crystallization[END_REF]. Numerical solution methods can be divided into two types: those predicting some properties of the CSD (moments-based methods) and those predicting the entire CSD (discretization methods mainly) [START_REF] Yazdanpanah | The Handbook of Continuous Crystallization[END_REF]. The first implies converting the original problem into ordinary differential equations via the moments transformation and the most probable crystal size distribution has to be estimated from the computed moments (inverse problem). On the other hand, discretization methods are generally known to have poor computing time performances and being highly dependent on the quantity and distribution of the nodes (discretization grid).

Most of the methods developed for the solution of the PBE aim at predicting the dispersed phase properties in a batch unit or approach the continuous process modelling as a particular case of the transient state analysis [START_REF] Ramkrishna | Population Balance Modeling: Current Status and Future Prospects[END_REF], [START_REF] Omar | Crystal Population Balance Formulation and Solution Methods: A Review[END_REF], [START_REF] Kumar | Comparison of numerical methods for solving population balance equations incorporating aggregation and breakage[END_REF]. The continuous steady state operation is rarely directly addressed in literature. For instance, [START_REF] Blandin | Modelling of agglomeration in suspension: Application to salicylic acid microparticles[END_REF] studied salicylic acid precipitation in a T-mixer modelled by a continuous well-mixed reactor followed by a plug-flow reactor operating at steady state [START_REF] Blandin | Modelling of agglomeration in suspension: Application to salicylic acid microparticles[END_REF]. The PBE was solved with the method of classes and the steady state in the mixing zone was given by the asymptote of the transient regime. [START_REF] Buffo | Simplified volume-averaged models for liquid-liquid dispersions_ Correct derivation and comparison with other approaches[END_REF], [START_REF] Shiea | Numerical Methods for the Solution of Population Balance Equations Coupled with Computational Fluid Dynamics[END_REF].

When considered, steady state PBE direct calculations are often approached under highly simplifying assumptions or via underperforming algorithms. M.J. [START_REF] Hounslow | A discretized population balance for continuous systems at steady state[END_REF] developed a discretization method for the PBE solution in the case of steady state continuous precipitation of Nickel Ammonium Sulphate [START_REF] Hounslow | A discretized population balance for continuous systems at steady state[END_REF]. The considered PBE includes nucleation, size-dependent growth and aggregation but computing speed drastically decreases with the number of nodes used in the discretization scheme. [START_REF] Nicmanis | Finite-element methods for steady-state population balance equations[END_REF]. In this case, the system of equations was solved by an iterative process and validated via numerical case studies and analytical solutions. The numerical performance of the traditional discretization method is improved but the selection of interpolation functions and collocation points remains a mathematical constraint. In order to test this numerical strategy, the oxalic precipitation (equation ( 2)) is considered as application case. In the nuclear industry, this reaction is widely used as the main conversion and purification stage leading to the production of an actinide riche powder from a multicomponent solution. Due to economic and environmental considerations, the study of actinides precipitation is traditionally approached by the investigation of harmless elements such as lanthanides [START_REF] Lalleman | Physical simulation of precipitation of radioactive element oxalates by using the harmless neodymium oxalate for studying the agglomeration phenomena[END_REF], [START_REF] Bertrand-Andrieu | Determination of Nucleation and Crystal Growth Kinetics in Hostile Environment &&num[END_REF].

Hence, the neodymium oxalate precipitation in a MSMPR precipitator is investigated in this study. The reaction kinetics is considered fast enough to neglect its influence on the entire process:

The structure of this paper is as follows: In Section 2, the PBEs are formulated and the corresponding numerical solution methods are described. The fixed point problem solving the agglomeration PBE at steady state is established. A suitable fixed point acceleration method is then introduced in order to improve the performance of the solution algorithm but also to ensure its robustness, especially when standard fixed point iterations would diverge. In Section 3, simulation results are first validated by comparison to analytical solutions and then tested against a range of experimental data. Finally, some concluding remarks and perspectives are pointed out.

Modelling strategy

During continuous precipitation, the CSD and continuous phase composition are the target variables to be predicted since they determine the quality of the produced powder and the recovery of the desired compound. CSD is detailed by the PBE as described in the last section. Furthermore, liquid phase concentrations need mass balances to be evaluated. Thus, the system of equations describing the precipitator includes both mass balances on the continuous and dispersed phases and the PBE.

One may note that crystallization kinetics vary with species concentrations in the continuous phase via the supersaturation. In parallel, the concentrations of the reagents remaining in the liquid phase are determined by the extent of the chemical reaction and further by the amount of solid accounted in the mass balances.

PBE solution

Equation ( 1) assembles all the physicochemical processes modifying the CSD, namely the transport phenomena and the crystallization mechanisms. The precipitator considered in this study has several inlet flows that may contain particles in suspension and one single exit flow ( ). The mean residence time is thus defined as and the outgoing flow rate as . Assuming an MSMPR precipitator, operating at steady state, with nucleation, size independent growth and agglomeration, the PBE can be written as:

with where and represent the particle size (m), is the Dirac delta function, the nuclei size (

) the nucleation rate, represents the agglomeration kernel (m

), and the density associated to the output stream. The term accounts for the particles appearing and disappearing by agglomeration, and respectively. Note that no assumption is made regarding the feed streams, therefore the PBE model solved in this work deals with single or multiple feedings whether they include crystals in suspension or not.

Equation ( 3) is a highly nonlinear integro-differential equation. In this work, the steady state PBE accounting for nucleation, growth and agglomeration is solved by a discretization method. Neodymium oxalate crystals have been identified as loose (open) agglomerates according to microscopic observations: the surface of crystals in contact with the surrounding solution is the same at any stage of the agglomeration process [START_REF] Lalleman | Physical simulation of precipitation of radioactive element oxalates by using the harmless neodymium oxalate for studying the agglomeration phenomena[END_REF]. Thus, the growth mechanism takes place at the surface of every elementary particle forming an agglomerate and is considered as non-affected by the internal diffusion. In other words, growth and agglomeration can be modelled as independent phenomena since the first affects the crystallites size whereas the second modifies the agglomerates size.

Considering this, the steady state PBE solution is performed by describing the whole precipitation process through two distinct PBEs solved successively. The first one accounts for nucleation and growth mechanisms and the second one considers the crystallites agglomeration only, using the CSD obtained as a result of the first PBE. Previously, a similar methodology was implemented in order to determine the neodymium oxalate agglomeration kernel [START_REF] Lalleman | Physical simulation of precipitation of radioactive element oxalates by using the harmless neodymium oxalate for studying the agglomeration phenomena[END_REF].

PBE for nucleation and growth mechanisms

The steady state PBE for a MSMPR undergoing only nucleation and size independent growth phenomena, leading to the crystallites size distribution function , is written as:

Assuming that crystals born at zero size, and that , equation ( 5) can be rewritten as:

with the boundary condition:

Equation ( 6) is an ordinary differential equation which can be easily solved when the inlet population ( ) and the nucleation and growth rates are known. The population fed to the reactor is set as an input variable, the nucleation and growth rates depend on the species concentration in the solution via the supersaturation. The concentration of the precipitating compound and the CSD flowing out the precipitator have also to satisfy the mass balance, both are initially unknown. The nucleation and growth PBE is computed over the discretized size vector by numerical integration ( is defined as the total number of nodes in the discretization grid). It is solved simultaneously with the mass balance by a minimisation iterative gradient-based method (Sequential Quadratic Programming method). The objective function aims to converge to a concentration in solution satisfying both the PBE and the mass balance. The stopping criterion is determined by the relative difference between the concentration calculated at two successive iterations.

PBE for agglomeration

The second steady state PBE concerns only the agglomeration phenomenon. It is expressed as:

From a theoretical point of view, the agglomeration process can be divided into two steps: particle collision and bonds consolidation. The loose agglomerates assumption implies that the amount of matter needed to bond crystallites does not modify the free surface of the initial particles. Thus, the amount of precipitating compound used to build these bonds can be neglected from the mass balance.

Consequently, the solution of the PBE for agglomerates requires the supersaturation, the composition of the output stream and the crystallites CSD as input information. Equation ( 8) is solved by the fixed pivot technique proposed by [START_REF] Kumar | On the solution of population balance equations by discretization-I. A fixed pivot technique[END_REF] [START_REF] Kumar | On the solution of population balance equations by discretization-I. A fixed pivot technique[END_REF]. The number-based density function ( ) is substituted by the total number of particles ( ), defined by the integration of ( ) over a discrete size interval :

with , the representative size for the k th size range obtained by discretization of the size domain in elements and and the lower and upper boundaries for the k th size range respectively. The right-hand side of equation ( 8) defines the agglomeration rate (m

-4 s -1 ):
The discretized agglomeration rate, still called , is given as proposed by [START_REF] Kumar | On the solution of population balance equations by discretization-I. A fixed pivot technique[END_REF] [START_REF] Kumar | On the solution of population balance equations by discretization-I. A fixed pivot technique[END_REF]:

under the following conditions:

where is defined as the Kronecker delta function, represents the total number of particles in the k th size range (class). The weight coefficient ensures the conservation of two properties (number and mass). In the agglomeration problem, the unknown variable is defined as which refers to the vector resulting from the discretization of the number-based CSD over the crystal size . By substituting equation [START_REF] Jaradat | Effect of phase dispersion and mass transfer direction on steady state RDC performance using population balance modelling[END_REF] in equation ( 8) and by defining the inlet population (before agglomeration process) as , equation ( 11) leads to a system of equations:

where the function represents the agglomeration PBE. The simplest way to solve the fixed point problem in equation ( 14) is to use fixed point iterations (or Picard iterations). Hence the j th iteration is defined as:

Fixed point iterates can be seen as pseudo-transient solutions (pseudo-time steps) before reaching the steady state. However, the convergence of the fixed point sequence strongly depends on the properties of the fixed point function (namely its spectral radius in the neighbourhood of the fixed point solution), which are generally a priori unknown. Hence, the convergence of the fixed point iterations (equation ( 15)) is not guaranteed.

In most precipitation processes, the objective variable varies over several orders of magnitude (between 0 and 10 25 in the case of the neodymium oxalate). For this reason, the following mixed convergence criterion is used as in [START_REF] Michel | A new numerical methodology for simulation of unstable crack growth in time independent brittle materials[END_REF]:

where the index k refers to the k th component of the vector. Thus the fixed point convergence is reached when the difference between two successive iterates (left-hand side term) is less than the total of the tolerances (right-hand side terms). The right-hand side terms include two tolerance constants: the relative tolerance and the absolute tolerance .This convergence criterion turns automatically from a relative error to an absolute error when becomes small. Hence it enables to deal in a unique criterion with high values of N components (for which the relative precision is recommended) and small values of N components (where absolute tolerance is required).

In practice, it is common to represent Equation ( 16) as a boolean:

It is worth noticing that when the peak region (high values of N components) is of interest, the absolute tolerance could remain relatively high, in order to focus on the convergence of the solution in the peak region (with a relative error).

Fixed point convergence acceleration methods

Fixed point iterations are a powerful tool to solve highly non-linear problems, since the algorithm can be applied to any objective function and no information about its derivative is required to implement the method. Nevertheless, in case of convergence, they are often characterized by a slow first order convergence which can become time consuming [START_REF] Brezinski | Numerical Analysis: Historical Developments in the 20th Century[END_REF], [START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF]. As the PBE concerned in this paper is a strongly non-linear and high dimensionality problem, acceleration methods are considered to improve both algorithm performances and robustness. Indeed, as shown in section 3, the fixed point algorithm applied to the agglomeration problem may result in a divergent sequence. Acceleration algorithms are then introduced to ensure fixed point convergence in such cases and to reduce the calculation time (number of iterations) when standard fixed point iterations converge.

In equation [START_REF] Nicmanis | Finite-element methods for steady-state population balance equations[END_REF], the standard iterates are defined by a recurrence relation and can be treated as a series. Then, classical sequence acceleration principles can be applied. In numerical analysis, two types of convergence acceleration methods are commonly studied: the static and the dynamic acceleration procedures. The second is commonly used in practice as it consists in substituting the current fixed point value ( ), by an accelerated iterate which improves the convergence rate [START_REF] Anderson | Iterative Procedures for Nonlinear Integral Equations[END_REF]. Many expressions for the accelerated sequence have been proposed in the literature depending on the characteristics of the initial mathematical problem. In general, it is defined as a function of the previous accelerated and standard iterates [START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF]. Ramiere et al. ( 2015) proposed a general framework to obtain accelerated sequence expressions from a nonlinear hybrid acceleration procedure applied to vector-based fixed point problems [START_REF] Ramière | Iterative residual-based vector methods to accelerate fixed point iterations[END_REF]. In this research work, three acceleration methods were tested to solve the PBE.

The first one is the classical scalar secant method applied independently on each component of the unknown vector . For a scalar fixed point problem, this method can be geometrically interpreted as finding the intersection of the line passing through the two points ( and ( , and the first bisector of the axes ( ). The abscissa of this intersection point defines the new iterate . It is hence mathematically defined as (see for example [START_REF] Bresinski | Accélération de la convergence en analyse numérique[END_REF]):

In this work, Equation ( 18) is applied independently on each scalar :

The second and third acceleration methods tested in this work can be obtained by the application of the secant principle when considering a vector fixed point problem. As each component does not define a scalar fixed point problem, vector sequence acceleration should be applied on the whole vector . The vector secant methods can be directly obtained by considering the (non unique) definition of the inverse of a vector:

where ||.|| represents a given norm, generally the Euclidean norm. Considering equations ( 18) and [START_REF] Kumar | On the solution of population balance equations by discretization-I. A fixed pivot technique[END_REF], two vector secant methods can be derived based on the two possible dot products arising in the numerator. In [START_REF] Ramière | Iterative residual-based vector methods to accelerate fixed point iterations[END_REF], it has been proven that these vector secants can be seen as a minimization process on fixed point residuals [START_REF] Ramière | Iterative residual-based vector methods to accelerate fixed point iterations[END_REF]. For the sake of clarity, the following residual notation is introduced: . The two vector-based secant acceleration methods obtained write:

 The alternate secant method:

 The crossed secant method:

The K-dimensional scalar secant approach (Equation ( 19)) prevented the divergence of the fixed point iterations but did not reduce significantly the number of iterations. The alternate secant method did not improve the fixed point iterations convergence at all, Finally, the crossed secant method was selected to solve agglomerates PBE as it presented really advantageous properties especially in case of divergence of the standard fixed point iterations (see Appendix 1).

As the fixed point sequence and the accelerated one are not mathematically constrained to fit the same limit [START_REF] Ramière | Iterative residual-based vector methods to accelerate fixed point iterations[END_REF], the convergence criterion defined in equation ( 17) must still focus on the standard fixed point residual to guarantee the solution of equation ( 14).

Results and discussion

Validation through analytical solutions

As described previously, the methodology developed in this work implies the formulation of two different population balances to describe precipitation process. As a first step, analytical solutions were employed to validate the numerical approach developed to solve each PBE separately. All the simulations were performed using the same standard personal computer.

PBE for nucleation and growth mechanisms in a cascade of two precipitators

Equation ( 6) may be solved analytically when the reactor is fed with a single flow of clear solution (

). In such case, the CSD flowing out the reactor is determined by [START_REF] Randolph | Transient and steady state size distributions in continuous mixed suspension crystallizers[END_REF]:

In order to test the methodology when the precipitator is fed with a suspension, the CSD flowing through a cascade of two precipitators is studied. Thus the CSD at the output of the second precipitator ( ) is determined analytically as [START_REF] Randolph | Theory of particulate processes: analysis and techniques of continuous crystallization[END_REF]:

Both analytical solutions (equations ( 23) and ( 24)) are compared to the CSD obtained by the numerical method detailed in section 2.1.1 (Figure 1). Simulation points are obtained with a discretization grid containing 1500 logarithmically distributed values between 10 -11 and 10 -4 m for the integration rule. The relative tolerance in the mass balances is set to 10 -6 and the calculation time is less than 10 seconds for both precipitators on an Intel Core i7 machine (1.90 GHz/2.11 GHz) with 32 Go of RAM. Figure 1 shows that the numerical and analytical results are identical with a relative difference less than 10 -4

.

Figure 1.

Comparison between the numerical solution methodology and the analytical solutions: output CSD for two reactors in series considering nucleation and size-independent growth.

PBE for agglomeration with size independent kernel

Steady state agglomeration problems can be solved analytically only when the agglomeration kernel is independent of the crystal size and the inlet population follows a normal distribution. In such case, [START_REF] Hounslow | A discretized population balance for continuous systems at steady state[END_REF] developed an analytical solution based on the crystal volume as internal coordinate and dimensionless variables [START_REF] Hounslow | A discretized population balance for continuous systems at steady state[END_REF]:

where represents the crystal volume, and are the modified Bessel functions of the first kind, of order zero and one respectively, the characteristic time of the agglomeration process, and are characteristics of the inlet population, namely the mean crystal volume and the total amount of particles.

Figure 2 a) presents a comparison between the analytical solution given in equation ( 25) and the numerical solution obtained through the accelerated method described in section 2.2 in terms of the total number of particles. Three values of the parameter are examined in order to verify results coherence and algorithm robustness. In every case, 200 logarithmically distributed points are used to define the discretization grid. The numerical and the analytical solutions match with a relative tolerance lower than 10 -4

.

Figure 2. Comparison between the CSD obtained by the developed numerical methodology and the analytical solution for the agglomeration problem.

In the context of a multicompartmental model, is expected to take values over 100. Additionally, CSD estimation has to be performed between the nuclei size (about 10 -9 m) and the agglomerates maximum size (about 10 -4 m), conditioning the number of intervals in the discretization to be large enough to track the smaller particles and the mean size region simultaneously. Figure 2 b) presents the comparison between the numerical and the analytical solutions for varying between 10 and 1000. Again, 200 logarithmically distributed points are used to define the discretization grid. With such values of t', the accelerated algorithm enables the agglomerates PBE to converge when standard fixed point iterations diverge.

Appendix 1 summarizes the main results obtained with the fixed point algorithm and the tested acceleration methods for two different values of t' (1 and 10 3 respectively), both and were chosen to match the order of magnitude observed in the case of the precipitation of neodymium oxalate. The crossed secant method clearly appears to be the most performant method to accelerate fixed point iteration convergence but also to guarantee sequence convergence.

It is worth underlying that a relevant initialization vector ( ) is also required to obtain the best performances. The determination of an optimal initial vector could be a hard mathematical task and return to the solution of the problem itself. Three numerical choices (null initialization, based only on crystallites CSD and based on both the crystallites CSD and the agglomeration dimensionless time ( )) have been tested (see Appendix 1). A null initialization vector ( ) appears to be a simple but robust choice.

Appendix 2, in which the acceleration algorithm is applied to a theoretical case ( and ),

shows that the efficiency (in terms of number of iterations) of the standard fixed point algorithm and the crossed secant method does not depend on the quantity of points used to discretize the internal coordinate (crystals size or volume), except for large agglomeration characteristic times. However, as expected the calculation time increases with the number of unknowns. One may observe that the number of iterations required to converge increases with the agglomeration dimensionless time ( ).

Actually, divergence of standard fixed point iterations is observed for . Also, the number of iterations to converge depends on the values of , and . For this reason, a different number of iterations is observed when on Appendix 1 and Appendix 2.

Application case: oxalate neodymium precipitation

Case study

Lalleman et al. (2012) presented results for neodymium oxalate precipitation [START_REF] Lalleman | Physical simulation of precipitation of radioactive element oxalates by using the harmless neodymium oxalate for studying the agglomeration phenomena[END_REF]. The experiments were carried out in a stirred tank reactor provided with a heating jacket and four stainless steel baffles.

Mixing was performed by a stainless-steel turbine equipped with four 45° pitched blades.

Oxalic precipitation was achieved by mixing an aqueous solution of neodymium nitrate with an aqueous solution of oxalic acid. For the purpose of this study, the reagents concentrations in the reactor obey to the stoichiometric ratio of equation ( 2) and the neodymium oxalate was the only compound found in the outgoing crystals. Mean residence time is about 1 minute and the steady state is typically reached after 15 minutes. Slurry samples were taken from the output stream and analysed with a laser diffraction granulometer in order to obtain the steady state CSD in the output flow.

The nucleation rate for the neodymium oxalate is expressed as [START_REF] Bertrand-Andrieu | Determination of Nucleation and Crystal Growth Kinetics in Hostile Environment &&num[END_REF]:

with the temperature (K) and the driving force, namely the relative supersaturation defined as:

where the activity coefficient is computed from the Bromley correlation [START_REF] Bromley | Thermodynamic properties of strong electrolytes in aqueous solutions[END_REF] with the individual contribution values proposed by Lalleman et al. (2012) for the neodymium electrolytes [START_REF] Lalleman | Determination of the Bromley contributions to estimate the activity coefficient of neodymium electrolytes[END_REF], is the molar concentration (mol.m -3

) and the solubility product (mol

5

.m -15

).

In the case of the oxalic precipitation of neodymium, the growth rate was found to be independent of the crystal size and integration-controlled [START_REF] Bertrand-Andrieu | Determination of Nucleation and Crystal Growth Kinetics in Hostile Environment &&num[END_REF]:

Similarly, the neodymium oxalate agglomeration process was found to be independent of the crystal size [START_REF] Lalleman | Physical simulation of precipitation of radioactive element oxalates by using the harmless neodymium oxalate for studying the agglomeration phenomena[END_REF], the agglomeration kernel is expressed as:

where is the shear rate (s -1 ) and the ionic force of the solution (mol.m -3

).

PBE modelling

In a first time, the developed methodology is tested against a reference experience (Experiment 1 in Table 1)Table 1. . Following the previous study (see also Appendix 1), the initial vector for the fixed point algorithm solving the agglomeration PBE is set to zero. Crystal mean size is correctly predicted under the experimental incertitude (Figure 4 a). The same observation is made concerning the moments, only the 0 th order moment is slightly over predicted by 5% (Figure 4 c)). In contrast, a broader difference is observed when the entire CSD is examined (Figure 4 b)). The gap between the simulated and the experimental CSD arises from the uncertainty associated to the crystallization kinetics. Firstly, nucleation kinetics typically predicts the quantity of elementary crystals over one order of magnitude. Secondly, the methodology employed to determine the neodymium oxalate agglomeration kernel includes the solution of the PBE by the moments method [START_REF] Lalleman | Physical simulation of precipitation of radioactive element oxalates by using the harmless neodymium oxalate for studying the agglomeration phenomena[END_REF]. Such numerical treatment ensures a precise calculation of the moments but generates errors in CSD prediction due to the lack of information about the entire distribution.

Sensitivity analysis

The developed algorithm is tested under several operating conditions. Experimental details and resulting mean sizes are presented in Table 1. Variation of three operating conditions are tested: the temperature, the shear rate via the stirring rate and the total neodymium concentration via the inlet stream concentration. Mean residence time, reagents proportions, and other experimental parameters remain the same as in the reference case (Experiment 1). 5. Calculation time is lower than 100 seconds for all the examined cases.

Also, the number of iterations is about 300 for the worst case and no direct dependency is observed with any of the studied variables. The proposed algorithm seems to be robust.

From Erreur ! Source du renvoi introuvable. Table 1, the mean size of the examined distributions increases with the concentration of Nd and with the temperature. In both cases, the crystallization kinetics are systematically strengthened: a higher concentration implies a higher supersaturation from equations ( 27), ( 29) and [START_REF] Lalleman | Determination of the Bromley contributions to estimate the activity coefficient of neodymium electrolytes[END_REF]. Similarly, the crystallization phenomena occur faster when the temperature is raised. In both cases, larger particles are expected. In contrast, it is not possible to establish a direct link between the shear rate and the mean crystal size, the latter could increase or decrease due to the compromise between mixing performance and mechanical resistance of the crystals. For this reason, a direct tendency is less observable when experiments 4 and 5 are considered. In all the cases the mean size is well predicted considering the experimental uncertainty and the evolution over the examined variables. Erreur ! Source du renvoi introuvable.Table 1.

Conclusions

In this work, a numerical methodology to solve the Population Balance Equation (PBE) including nucleation, size-independent growth and loose agglomeration in a Mixed Suspension Mixed Product Removal (MSMPR) precipitator was developed for the steady state case. It is based on the solution of two coupled PBEs by a discretization method. The PBE formulation takes into account several input streams and no assumption is made about the solid phase entering to the control volume.

The PBE for the crystallites raises in an ordinary differential equation which can be easily solved numerically. Concerning the agglomerates, the PBE is a strongly non-linear and high dimensionality problem. Hence, a fixed point problem is formulated and acrossed-secant acceleration method is included in order to guarantee both convergence and efficiency. No restriction is placed on the location of grid points (quantity or distribution) and a mixed convergence criterion is formulated in order to favour the prediction of the mean size region rather than the furthest portions of the crystal size distribution (CSD). The complete system is solved through the coupling of both PBEs.

The methodology was first validated against analytical solutions. The entire CSD of the crystallites and agglomerates are accurately predicted. The crossed-secant accelerated fixed point improves standard fixed point robustness to solve the agglomeration PBE, especially for high values of the agglomeration characteristic time.

As a second step, the CSD obtained by the developed methodology was compared to a wide range of experimental data. Only the crossed-secant acceleration method achieves agglomeration PBE solution. The mean crystal size and the first seven moments of the CSD are accurately predicted by the numerical method. A slight difference is observed when the entire distribution is compared, which can be explained by the error propagation induced by the nucleation, growth and agglomeration kinetics models and the inherent experimental uncertainty. The convergence performances are independent of the operating conditions when the acceleration method is included on the PBE solution algorithm for agglomeration. The proposed numerical methodology is then has proven to be efficient and robust. Moreover, the computing performances are really satisfactory and the assumptions made are suitable for future coupling with a multicompartmental model or CFD calculations [START_REF] Saikali | Validation of the hydrodynamics in a turbulent un-baffled stirred tank: A necessity for vortex-reactor precipitation studies[END_REF].

Figure 3

 3 presents the behaviour of the agglomeration PBE convergence criterion (equation (17)) over the number of iterations when the developed methodology is applied to a MSMPR precipitator including nucleation, growth and agglomeration mechanisms. The discretization grid contains 1500 logarithmically distributed points between 10 -11 m and 10 -3 m. Standard fixed point iterations (Figure 3 (a)) do not succeed to converge due to the value of the dimensionless agglomeration time. In contrast, the crossed-secant accelerated fixed point achieves the agglomerates CSD prediction after 237 iterations (Figure 3 (b)). Simulation time for the case study is about 80 seconds on an Intel Core i7 machine (1.90 GHz/ 2.11 GHz) with 32 Go of RAM, the relative tolerance ( ) is fixed to 10 -2 while the absolute tolerance is defined as a function of the crystallites CSD:

Figure 3 .

 3 Figure 3. Convergence criterion evolution in the case of a MSMPR modelling including nucleation, sizeindependent growth and loose agglomeration with (a) standard fixed point algorithm b) accelerated fixed point algorithm.

Figure 4 .

 4 Figure 4. Experimental and simulation results in the case of a MSMPR including nucleation, size-independent growth and loose agglomeration: a) Mean crystal size prediction over the number of iterations b) Mass-based CSD and c) CSD moments

Figure 5 .

 5 Figure 5. Comparison between experimental and simulated mean crystal sizes for the conditions reported in

  Mother particles in agglomeration

  Primary and secondary nucleations and size dependent growth were considered in the model. The steady state behaviour of polypropylene particle size was investigated by[START_REF] Luo | Steady-state particle size distribution modeling of polypropylene produced in tubular loop reactors[END_REF] [START_REF] Luo | Steady-state particle size distribution modeling of polypropylene produced in tubular loop reactors[END_REF]. In this case, growth and attrition were both considered in a tubular loop reactor and the PBE was solved by a discretization method. The polymer molecular weight distribution and the particle size distribution were approached by solving the transient state population balances over an extended time scale. A Rotator Disc Contactor column operating at steady state was reported by[START_REF] Jaradat | Effect of phase dispersion and mass transfer direction on steady state RDC performance using population balance modelling[END_REF] where a multivariate PBE was coupled with hydrodynamics and applied to three chemical systems[START_REF] Jaradat | Effect of phase dispersion and mass transfer direction on steady state RDC performance using population balance modelling[END_REF]. Again, the steady state PBE solution is based on asymptotic transient calculations. It is worth noticing that the dispersed phase properties

obtained via transient state estimation are generally employed as initial conditions when coupled CFD-PBE steady state simulations are performed

Table 1 .

 1 Experimental conditions, mean crystal size and numerical performances for the simulation of a continuous crystallizer including nucleation, size-independent growth and loose agglomeration.

	Experiment	[Nd] 0 (mol m -3 )	T (K)	ϒ (s -1 )	Experimental d 43 (µm)	Simulated d 43 (µm)	Number of iterations (Agglomeration	Calculation time (s)
							PBE)
	1	62.0	293.15	362	45	43	237	76
	2	201.7	293.15	362	53	53	272	92
	3	284.9	293.15	362	59	56	274	83
	4	142.2	293.15	665	48	47	319	97
	5	142.2	293.15 1024	41	46	250	84
	6	144.7	303.15	362	55	53	303	92
	7	144.7	313.15	362	54	57	317	91
	8	144.7	323.15	362	60	60	306	89
	Experimental and numerical volume-based mean sizes (d 43 ) obtained for the experiments described in
	Table 1 are plotted in Figure

Acknowledgments

The authors thank ORANO for its financial support.

Abbreviations


Arnaud Casteigts
email: casteig@site.uottawa.ca

Serge Chaumette
email: serge.chaumette@labri.fr

Afonso Ferreira
email: afonso.ferreira@sophia.inria.fr

Characterizing Topological Assumptions of Distributed Algorithms in Dynamic Networks ⋆

Keywords: Dynamic networks, distributed algorithms, evolving graphs, local interactions, topological assumptions

Besides the complexity in time or in number of messages, a common approach for analyzing distributed algorithms is to look at their assumptions on the underlying network. This paper focuses on the study of such assumptions in dynamic networks, where the connectivity is expected to change, predictably or not, during the execution. Our main contribution is a theoretical framework dedicated to such analysis. By combining several existing components (local computations, graph relabellings, and evolving graphs), this framework allows to express detailed properties on the network dynamics and to prove that a given property is necessary, or sufficient, for the success of an algorithm. Consequences of this work include (i) the possibility to compare distributed algorithms on the basis of their topological requirements, (ii) the elaboration of a formal classification of dynamic networks with respect to these properties, and (iii) the possibility to check automatically whether a network trace belongs to one of the classes, and consequently to know which algorithm should run on it.

Introduction

The past decade has seen a considerable research effort devoted to the design of distributed algorithms and protocols targeting dynamic network topologies. It appears, however, that most of the assumptions considered when examining algorithm requirements still relate to static properties such as the size, density, or geometry of the target network. Assumptions that really relate to the network dynamics, when any, are generally stated using non-formal expressions, such as "nodes are expected to move slowly enough to..", or "nodes cannot leave the network for a long time", or by assuming a given mobility model (whose concrete topological implications remain unclear). The dynamic properties highlighted in such a way make the fair comparison of algorithm requirements rather complicated and often ambiguous.

Our work aims at providing general formalisms and methods for studying fundamental properties of dynamic networks, and more particularly their impact on distributed systems. As an illustrative example, let us consider the broadcasting of an information within the dynamic network depicted by Figure 1. The possibility to complete the broadcast in this scenario clearly depends on which node is the initial emitter: a and b may succeed, while c cannot. Why? How can we formulate this intuitive property that the topology evolution must have with regard to the emitter and other nodes? How can we formally prove it as a necessary condition to obtain broadcast completion? While rather simple, such a characterization might be difficult to obtain with usual graph formalisms and computation models. This paper introduces a theoretical framework dedicated to such kind of analyses. This framework is intended to serve as a general basis for studying fundamental properties of distributed algorithms in dynamic networks. Contrary to the work in [AAD + 06], where the authors first make a strong topological assumption (all pairs of nodes repeatedly meet during the execution) and then characterize the solvable problems in this context, we consider the exact opposite approach by studying, for a given solution (i.e., algorithm), the necessary and/or sufficient conditions it requires on the topology. To the best of our knowledge, this is the first attempt in such a direction.

The strength of the proposed framework lies in its basic components, which are an appropriate combinatorial model to represent dynamic topologies (evolving graphs [START_REF] Ferreira | Building a reference combinatorial model for MANETs[END_REF]), and a very high-level interaction model to describe distributed operations (local computations, with the associated formalism of graph relabellings [START_REF] Litovsky | Graph relabelling systems and distributed algorithms[END_REF]). The next section is devoted to the presentation of these existing components. In Section 3 we combine them to set up the new analysis framework. This framework is then applied in Section 4 to the analysis of three basic algorithms (one propagation and two enumeration algorithms), whose intuitively apparent properties are here formally characterized. Based on the analysis results, Section 5 shows how algorithms can be compared on the basis of their topological requirements, and reciprocally how dynamic networks can be classified according to the algorithms they support. Finally, we discuss the possibility to check automatically the inclusion of a given network trace to one of the classes. Section 6 concludes with some avenues for further research.

Related work

We describe here the formalisms and theoretical tools that compose the proposed analysis framework: Local Computations to abstract the communication model, Graph Relabelling Systems to describe local computation algorithms, and Evolving Graphs to express properties on dynamic topologies. Their comprehension is required to ensure a clear understanding of the following sections, where they are combined together.

Abstracting communications through local computations and graph relabellings

Distributed algorithms can be expressed using a variety of communication models (e.g. mailbox, shared memory, and message passing). Whereas a vast majority of algorithms is designed in one of these models (predominantly the message passing model), the very fact that one of them is chosen implies that the obtained results (e.g. positive or negative characterizations and associated proofs) are limited to the scope of this model. This problem of diversity among formalisms and results, already pointed out twenty years ago in [START_REF] Lynch | A hundred impossibility proofs for distributed computing[END_REF], led researchers to consider higher abstractions when studying fundamental properties of distributed systems.

Local computations and Graph relabellings were jointly proposed in this perspective in [START_REF] Litovsky | Graph relabelling systems and distributed algorithms[END_REF]. These theoretical tools allow to represent a distributed algorithm as a set of local interaction rules that are independent from the effective communications. Within the formalism of graph relabellings, the network is represented by graph whose vertices and edges are associated with labels that represent the algorithmic state of the corresponding nodes and links. An interaction rule is then defined as a transition pattern (preconditions, actions), where preconditions and actions relate to the label values. Since these interactions are local, each transition pattern must involve a limited and connected subset of vertices and edges. Figure 2 shows different scopes for the transition patterns, which are not necessarily the same for preconditions and actions.

More formally, let the network topology be represented by a finite undirected loopless graph G = (V G , E G ), with V G representing the set of nodes and E G representing the set of communication links between them. Two vertices u and v are said neighbors if and only if they share a common edge {u, v} in E G . Let λ : V G ∪ E G → L * be a mapping that associates every vertex and edge from G with one or several labels from an alphabet L (which denotes all the possible states these elements can take). The state of a given vertex v, resp. edge e, at a given time t is thus denoted by λ t (v), resp. λ t (e). The whole labelled graph is represented by the pair (G, λ), noted G. The dashed elements represent entities (vertices or edges) that are considered by preconditions but remain unaffected by actions. The reader is referred to [START_REF] Chalopin | Local computations in graphs: The case of cellular edge local computations[END_REF] for a comparative study of some of these models.

According to [START_REF] Litovsky | Graph relabelling systems and distributed algorithms[END_REF], a complete algorithm can be given by a triplet {L, I, P }, where I is the set of initial states, and P is a set of relabelling rules (transition patterns) representing the distributed interactions. The Algorithm 1 below (A 1 for short), gives the example of a one-rule algorithm that represents the general broadcasting scheme discussed in the introduction. We assume here that the label I (resp. N ) stands for the state informed (resp. non-informed). Propagating the information thus consists in repeating this single rule, starting from the emitter vertex, until all vertices are labelled I. 4Algorithm 1 A propagation algorithm coded by a single relabelling rule (r 1 ). Regarding the organization of collaborations between nodes, it is important to note that the algorithm specification does not stipulate how the nodes must collaborate, i.e., the way they select each other to perform a common computation step. From the abstraction level of local computations, this underlying synchronization is seen as an implementation choice, which implies that local computation algorithms may not be deterministic at this level. As discussed later on, characterizing sufficient conditions will require additional assumptions on this underlying layer (which is not the case for necessary conditions).

Expressing dynamic network properties using Evolving Graphs

In a different context, evolving graphs [START_REF] Ferreira | Building a reference combinatorial model for MANETs[END_REF] have been proposed as a combinatorial model for dynamic networks. The initial purpose of this model was to provide a suitable representation of fixed schedule dynamic networks (FSDNs), in order to compute optimized communication schemes such as shortest, fastest and foremost paths. In such a context, the evolution of the network was known beforehand. In the present work, we propose to use evolving graphs for a different purpose, namely to express topological properties in dynamic networks. It is important to keep in mind that the analyzed algorithms are never supposed to know the evolution of the network ahead of time.

An evolving graph is a structure in which the changing connectivity of a dynamic network is recorded (see Figure 3). More formally, let S T = t 0 , t 1 , ..., t n be a sequence of dates in T (usually R + ). Except for t 0 and t n , all these dates correspond to a topological event that modifies the network. Let S G = G 0 , G 1 , ..., G n-1 be the corresponding sequence of graphs, with each G i the graph corresponding to the period [t i , t i+1 [. Finally, let us denote by G (alone) the union graph of all G i , called the underlying graph. Then the triplet G = (G, S G , S T ) is the corresponding evolving graph. As shown in Figure 3, this graph can be represented by the underlying graph G whose edges and vertices (only edges here) are associated with their presence interval indices. Henceforth, we will use the notations V G and E G to denote V (G) and E(G), the sets of all vertices and edges that exist at some point of the network life. Note that whereas used as undirected in this paper, evolving graphs were initially introduced as directed, and considered also bandwidth restrictions on edges, which is not used here. Further definitions on evolving graphs (given an evolving graph G = (G, S G , S T )).

Predecessor of a date: for any date d in T [ti,ti+1[ , with t i , t i+1 ∈ S T , we say that t i is the predecessor of d in S T , and we note pred(d) = t i .

Journey: given a sequence of couples J = {(e 1 , σ 1 ), ..., (e i , σ i ), ..., (e k , σ k )} composed of edges from E G and dates from the continuous domain T, J is called a journey if and only if σ 1 , σ 2 , ..., σ k is non-decreasing and for all i in 1..k, e i ∈ E(G pred(σi) ), that is, the edge e i exists at time σ i . Less formally, a journey can be thought of as a path over-time from one vertex to another. A journey from a vertex u to a vertex v is noted J (u,v) .

Discrete Journey: a discrete journey is a journey so that every date of the sequence σ1, ..., σ k is in S T , instead of T. It allows to represent in a single entity all the possible journeys occurring on the same sequence of edges during the same sequence of intervals. This also allows to consider such entities as subgraphs of the evolving graph G, and to note J ⊆ G. The point is that every normal journey {(e1, σ1), ..., (ei, σi), ..., (e k , σ k )} can be associated with a discrete journey {(e1, pred(σ1)), ..., (ei, pred(σi)), ..., (e k , pred(σ k ))} ⊆ G, and every discrete journey implies an infinity of normal journeys for the corresponding edges and intervals.

Strictness of a discrete journey: a discrete journey is said strict, noted J strict , if its sequence of dates σ 1 , σ 2 , ..., σ k is strictly increasing.

To give a few examples on the graph of Figure 3, -J (a,e) ={(ac,

σ 1 ∈ [t 0 , t 1 [), (ce, σ 2 ∈ [t 2 , t 3 [)}
is a normal journey from a to e ; -J strict(a,e) ={(ac, 0), (ce, 2)} is a discrete (and strict) journey from a to e ; -J (a,e) ={(ac, 0), (cd, 0), (de, 3)} is a discrete (non-strict) journey from a to e ; -J strict(a,e) ={(ac, 0), (cd, 1), (de, 3)} is a discrete (and strict) journey from a to e.

Note that journeys are naturally oriented, in the sense that a journey from one vertex to another does not imply the existence of a journey in the reverse direction (e.g. from e to a). From this point on, unless said explicitly, we will only consider discrete journeys, and denote them by the sole term journey.

The proposed analysis framework

As a recall of the previous section, the algorithmic state of the network is given by a labelling on the corresponding graph G, then noted G. As another recall, we denote by G i the graph covering the period [t i , t i+1 [ in the evolving graph G = (G, S G , S T ), with G i ∈ S G and t i , t i+1 ∈ S T . Note that the notation G was used here with two different meanings: the first as the generic letter to represent the network, the second to denote the underlying graph of G. Both notations are kept in the following, while preventing the text from ambiguities.

Putting the pieces together: relabellings over evolving graphs

For an evolving graph G = (G, S G , S T ) and a given date index i | t i ∈ S T , we denote by G i the labelled graph (G i , λ ti+ǫ ) representing the network state just after the topological event of date t i , and by G i[ the labelled graph (G i-1 , λ ti-ǫ ) representing the network state just before it. We note,

Event ti (G i[ ) = G i .
A number of distributed operations can occur between two consecutive events. Hence, for a given algorithm A and two consecutive dates t i , t i+1 ∈ S T , we denote by R A [t i ,t i+1 [ the relabelling sequence induced by A on the graph G i during the period [t i , t i+1 [, and have,

R A [t i ,t i+1 [ (G i ) = G i+1[ .
For the sake of simplicity, we authorize the notation r i (u, v) ∈ R A [t,t ′ [ to denote the fact that a rule r i is applied on the edge (u, v) during [t, t ′ [. A complete execution sequence from t 0 to t last is given by an alternated sequence of relabelling steps and topological events, noted,

X=R A [t last-1 ,t last [ •Event t last-1 •..•Event t i •R A [t i-1 ,t i [ •..•Event t 1 •R A [t 0 ,t 1 [ (G 0 )
The combined formalism is summed up on Figure 4. As mentioned at the end of Section 2.1, the execution of a local computation algorithm is not necessarily deterministic, and may depend on the way nodes select one another at a lower level. Hence, we denote by X A/G the set of all possible execution sequences of an algorithm A over an evolving graph G.

time start t0 G0 G 1[ R [t 0 ,t 1 [ z }| { G0 Evt 1 t1 G1 G 2[ R [t 1 ,t 2 [ z }| { G1 Evt 2 t2 Evt last-1 t last-1 G last-1 G last[ R [t last-1 ,t last [ z }| { Gt last-1 end t last . . . . . .

Characterizing the topological assumptions of an algorithm

Below are some proposed methods and additional concepts to characterize the requirement of an algorithm in terms of topology dynamics. More precisely, we use the new combined formalism to define the notions of topology-related necessary and sufficient conditions, and discuss how they can be proved.

Objectives of an algorithm. Given an algorithm A and a labelled graph G, the state one desires to reach can be given by a logic formula P on the labels of vertices and/or edges. In the case of the propagation scheme (Algorithm 1 Section 2.1), this could be that all nodes are informed,

P 1 (G) = ∀v ∈ V (G), λ(v) = I ,
Now, if the objective (noted O) is to reach such state at some point, then it can be simply expressed as P to be satisfied on the very last labelled graph of G (e.g. O A1 = P 1 (G last ) in the example). Whereas not covered in the examples, one could also consider algorithms whose objectives are to maintain a state (e.g. selfstabilizing algorithms), and express it for example as

O A = ∀G i ∈ S G , P(G i+1[ ).
Necessary conditions. Given an algorithm A, its objective O A , an evolving graph G and an evolving graph property C N . The property C N is a (topology-related) necessary condition for O A if and only if

∀G, ¬C N (G) =⇒ ¬O A Proving this result comes to prove that ∀G, ¬ CN (G) =⇒ ∄X ∈ X A/G | P(G last ).
Sufficient conditions. Symmetrically, an evolving graph property C S is a (topologyrelated) sufficient condition for A if and only if

∀G, C S (G) =⇒ O A Proving this result comes to prove that ∀G, C S (G) =⇒ ∀X ∈ X A/G , P(G last ).
Discussion. No topology of any kind can guarantee, alone, that the nodes will effectively communicate and collaborate with each other. Hence, the characterization of any sufficient condition necessarily requires to make additional assumptions on the collaboration of nodes. We propose below a generic such assumption for the pairwise interaction model (depicted on Figure 2(c)). This assumption may or may not be considered as realistic depending on the expected rate of topological changes.

Progression Hypothesis 1 (P H 1 ). For every given time interval

[t i , t i+1 [, with t i in S \{t last } T
, every vertex will be able to apply at least one relabelling rule with each of its neighbors, provided the rule preconditions are already satisfied at time t i (and still satisfied at the time the rule is applied).

First applications of the proposed framework

This section illustrates the proposed framework by the analysis of three basic algorithms, namely the propagation algorithm previously given, and two enumeration algorithms (one centralized, the other decentralized). The results obtained here are used in the next section to highlight some implications of this work.

Analysis of the propagation algorithm

We want to prove that the existence of a journey (resp. strict journey) between the emitter and every other node is a necessary (resp. sufficient) condition to achieve O A1 (complete the propagation). The point here is to show how these intuitive conditions can be formally established.

Condition 1 ∀v ∈ V \{emitter} G
, ∃J (emitter,v) ⊆ G (It exists a journey between the emitter and every other vertex).

Lemma 1 ∀v ∈ VG | λt 0 (v) = N, ∀σ ∈ T [t 0 ,t last [ , λσ(v) = I =⇒ ∃u ∈ V \{v} G , σ ′ ∈ T [t 0 ,σ[ | λ σ ′ (u) = I, ∃J (u,v) ⊆ G (If a
non-emitter vertex has the information at some point, it implies the existence of an incoming journey from a vertex that had the information before)

Proof. ∀v ∈ VG | λt 0 (v) = N, ∀σ ∈ T [t 0 ,t last [ , (λσ(v) = I =⇒ ∃v ′ ∈ V \v G | r 1(v ′ ,v) ∈ R A 1 [t 0 ,σ[ ) (If a
non-emitter vertex has the information at some point, then it has applied rule r1 with another vertex)

=⇒ ∃v ′ ∈ V \v G , σ ′ ∈ T [t 0 ,σ[ | λ σ ′ (v ′ ) = I, (v ′ , v) ∈ E(G pred(σ ′ ) )
(An edge existed at a previous date between this vertex and a vertex labelled

I) By repetition, =⇒ ∃v ′′ ∈ V \v G , σ ′′ ∈ T [t 0 ,σ[ | λ σ ′′ (v ′′ ) = I, ∃J (v ′′ ,v) ⊆ G (A
journey existed from a node that had the information to the considered node) Proposition 1 Condition 1 (C1) is a necessary condition on G to allow Algorithm 1 (A1) to reach its objective OA 1 .

Proof. (using Lemma 1). Following from Lemma 1 and the initial states (I for the emitter, N for all other vertices), we have OA 1 =⇒ C1, and then

¬C1 =⇒ ¬OA 1 Condition 2 ∀v ∈ V \{emitter} G , ∃J strict(emitter,v) ⊆ G
Proposition 2 Assuming the progression hypothesis (P H1, defined in the previous section), Condition 2 (C2) is sufficient on G to guarantee that A1 will reach OA 1 .

Proof. (1): By

P H1, ∀ti ∈ S \(t last ) T , ∀(u, u ′ ) ∈ E(Gi),(λt i (u) = I =⇒ λt i+1 (u ′ ) = I) By iteration on (1): ∀u, v ∈ VG, (∃J strict(u,v) ⊆ G) =⇒ (λt 0 (u)=I =⇒ λt last (v)=I) Now, because λt 0 (emitter) = I, we have C2(G) =⇒ ∀X ∈ X A/G , P1(G last )

Analysis of a centralized enumeration algorithm

Like the propagation algorithm, the distributed algorithm presented below assumes that one distinguished vertex is given a different initial state. This vertex, called the counter, is in charge of counting all the vertices it meets during the execution (its successive neighbors in the changing topology). Hence, the counter vertex has two labels (C, i), meaning that it is the counter (C), and that it has already counted i vertices (initially 1, i.e., itself). The other vertices are labelled either F or N , depending on whether they have already been counted or not, respectively. The counting rule is given by r 1 in Algorithm 2, below.

Algorithm 2 Enumeration algorithm with a pre-selected counter. N } ((C,1) for the counter, N for all other vertices) alphabet: {C, N, F, N * } rule r1:

C, i N C, i + 1 F
Objective of the algorithm. Under the assumption of a fixed number of vertices, the algorithm reaches the desired state when all vertices are counted, which corresponds to the fact that no more vertices are labelled N :

P 2 = ∀v ∈ V (G), λ(v) = N
The objective of Algorithm 2 is then to satisfy this property at the end of the execution (O A2 = P 2 (G last )). We want to prove here that the existence of an edge at some point of the execution between the counter node and every other node is a necessary and sufficient condition.

Condition 3 ∀v ∈ V \{counter} G
, ∃ti ∈ S T | (counter, v) ∈ E(Gi), or equivalently with the notion of underlying graph, ∀v ∈ V \{counter} G , (counter, v) ∈ EG Proposition 3 For a given evolving graph G representing the topological evolutions that take place during the execution of A2, Condition 3 (C3) is a necessary condition on G to allow A2 to reach its objective OA 2 .

Proof.

¬C3(G) =⇒ ∃v ∈ V \{counter} G | (counter, v) / ∈ E(G) =⇒ ∃v ∈ V \{counter} G | ∀ti ∈ S \{t last } T , r1(counter, v) / ∈ R A 2 [t i ,t i+1 [ =⇒ ∃v ∈ V \{counter} G | ∀X ∈ X A 2 /G , λt last (v) = N =⇒ ∄X ∈ X A 2 /G | P2(G last ) =⇒ ¬OA 2
Proposition 4 Assuming the progression hypothesis (P H1), C3 is also a sufficient condition on G to guarantee that A2 will reach its objective OA 2 .

Proof. C3(G) =⇒ ∀v ∈ V \{counter} G , ∃ti ∈ S T | (counter, v) ∈ E(Gi) by P H1, =⇒ ∀v ∈ V \{counter} G , ∃ti ∈ S T | r1(counter, v) ∈ R A 2 [t i ,t i+1 [ =⇒ ∀v ∈ V \{counter} G , λt last (v) = N =⇒ ∀X ∈ X A 2 /G , P2(G last ) =⇒ OA 2

Analysis of a decentralized enumeration algorithm

Contrary to the previous algorithm, Algorithm 3 below does not require a distinguished initial state for any vertex. Indeed, all vertices are initialized with the same labels (C, 1), meaning that they are all initially counters that have already included themselves into the count. Then, depending on the topological evolutions, the counters opportunistically merge by pairs (rule r 1 ) in Algorithm A 3 . In the optimal case, at the end of the execution, only one node remains labelled C and its second label gives the total number of vertices in the graph. A similar counting principle has been used in [AAD + 06] to monitor a flock of birds for fever, with the role of counters being played by sensors that have measured a high temperature level.

Algorithm 3 Decentralized enumeration algorithm. initial states: {(C, 1)} (for all vertices); alphabet: {C, F, N * } rule r1:

C, i C, j C, i + j F
Objective of the algorithm Under the assumption of a fixed number of vertices, this algorithm reaches the desired state when exactly one vertex remains labelled C:

P 3 = ∃u ∈ V G | ∀v ∈ V \{u} G , λ(u) = C, λ(v) = C, and O A3 = P 3 (G last )
For the sake of simplicity, we introduce one additional definition: the destination set of a vertex v in an evolving graph G is the set of all the vertices that can be reached from v by a journey, noted Dest G (v). Note that v ∈ Dest G (v) through an empty journey. We want to prove here that the existence of a journey from every vertex to at least one common destination vertex is a necessary condition for this algorithm.

Condition 4 ∃v ∈ VG | ∀u ∈ VG, v ∈ DestG(u) Lemma 2 ∀u ∈ VG | λt i (u) = C, ∃u ′ ∈ DestG(u) | λt j≥i (u ′ ) = C
(Whatever the C-labelled vertex considered at some point, there will be at a later point of the execution at least one vertex labelled C among its destination vertices)

Proof. (by contradiction). The application of r1 is the only operation that can suppress a counter, while preserving the other counter in the pair. If Lemma 2 was false, then it would imply either that both counters have been discarded by r1 at some point, or that the relabelling sequence has occurred from a C-labelled vertex towards a vertex that is outside of its destination set. Both are impossible.

Proposition 5 Condition 4 (C4) is necessary for A3 to reach its objective OA 3 .

Proof. (using Lemma 2). ¬C4(G) =⇒ ∄v ∈ VG | ∀u ∈ VG, v ∈ DestG(u) (no vertices are destination for all the others).

=⇒ ∀v ∈ VG | λt last (v) = C, ∃u ∈ VG | v /
∈ DestG(u) (Whatever the final counter, there is a vertex that could not reach it by a journey).

Now, thanks to Lemma 2, =⇒ ∀v ∈ VG | λt last (v) = C, ∃v ′ ∈ V \{v} G | λt last (v ′ ) = C (There are at least two final counters). =⇒ ¬P3(G last ) =⇒ ¬OA 3
The characterization of a sufficient condition for A 3 is left open. We believe such a condition exists, but would be satisfied on a very few specific graphs.

Applications of the analysis results

This section presents some applications of the new framework. In particular, we show how the previously characterized conditions can be used to define evolving graph classes, some of which are included in others. This leads to a de facto classification of dynamic networks according to the algorithms they support. The relations between classes can be used in turn to compare algorithms on the basis of their topological requirements. Finally, we propose a method to check a given network trace for inclusion in each introduced class.

From conditions to graph classes

From C1 = ∀v ∈ V \{emitter} G
, ∃J (emitter,v) ⊆ G, we derive two classes of evolving graphs. F 1 is the class in which at least one vertex can reach all the others by a journey. If an evolving graph does not belong to this class, then there is no chance for A 1 to succeed whatever the initial emitter. F 2 is the class where every vertex can reach all the others by a journey. If an evolving graph does not belong to this class, then at least one vertex, if chosen as an initial emitter, is guaranteed to fail to inform all the others using A 1 .

From C2 = ∀v ∈ V \{emitter} G
, ∃J strict(emitter,v) ⊆ G, we derive two classes of evolving graphs. F 3 is the class in which at least one vertex can reach all the others by a strict journey. If an evolving graph belongs to this class, then there is at least one vertex that could, for sure, inform all the others using A 1 (assuming the progression hypothesis). F 4 is the class of evolving graphs in which every vertex can reach all the others by a strict journey. If an evolving graph belongs to this class, then the success of A 1 is guaranteed for any vertex as initial emitter (again, if the progression hypothesis is assumed).

From C3 = ∀v ∈ V \{counter} G , (counter, v) ∈ EG, we derive two classes of graphs. F 5 is the class of evolving graphs in which at least one vertex shares, at some point of the execution, an edge with every other vertex. If an evolving graph does not belong to this class, then there is no chance of success for A 2 , whatever the vertex chosen for counter. Here, if we assume the progression hypothesis, then F 5 is also a class in which the success of the algorithm can be guaranteed for one specific vertex as counter. F 6 is the class of evolving graphs in which every vertex shares an edge with every other vertex at some point of the execution. If an evolving graph does not belong to this class, then there exists at least one vertex for which, if it is chosen as the counter, the failure of A 2 is guaranteed. Again, if we consider the progression hypothesis, then F 6 becomes a class in which the success is guaranteed whatever the counter.

Finally, from C4 = ∃v ∈ VG | ∀u ∈ VG, v ∈ DestG(u), we derive the class F 7 , which is the class of graphs such that at least one vertex can be reached from all the others by a journey. If a graph does not belong to this class, then there is absolutely no chance of success for A 3 .

Relations between classes

Since all implies at least one, we have: F2 ⊆ F1, F4 ⊆ F3, and F6 ⊆ F5. Since a strict journey is a journey, we have: F3 ⊆ F1, and F4 ⊆ F2. Since an edge is a (strict) journey, we have: F5 ⊆ F3, F6 ⊆ F4, and F5 ⊆ F7. Finally, the existence of a journey between all pairs of vertices (F2) implies that each vertex can be reached by all the others, which implies in turn that at least one vertex can be reach by all the others (F7). We then have: F2 ⊆ F7. Although we have used here a non-strict inclusion (⊆), the inclusions described above are strict (⊂). This can be easily proved by finding for each inclusion a graph that belongs to the parent class but is outside the child class. Figure 5 summarizes all these relations.

F1 : ∃u ∈ VG | ∀v ∈ V \{u} G , ∃J (u,v) ⊆ G F2 : ∀u, v ∈ VG, ∃J (u,v) ⊆ G F3 : ∃u ∈ VG | ∀v ∈ V \{u} G , ∃J strict(u,v) ⊆ G F4 : ∀u, v ∈ VG, ∃J strict(u,v) ⊆ G F5 : ∃u ∈ VG | ∀v ∈ V \{u} G , (u, v) ∈ EG F6 : ∀u, v ∈ VG, (u, v) ∈ EG F7 : ∃u ∈ VG | ∀v ∈ V \{u} G , u ∈ DestG(v) F6 F4 F5 F2 F3 F7 F1
Fig. 5. A first classification of dynamic networks, based on evolving graph properties that result from the analysis of three distributed algorithms (arrows denote inclusion).

Comparison of algorithms according to topological assumptions

Let us consider the two enumeration algorithms given in Section 4. To have any chance of success, A 2 requires the evolving graph to be in F 5 (and a fortunate choice of counter) or in F 6 (for possibly any vertex as counter). On the other hand, A 3 requires the evolving graph to be in F 7 . Now, both classes F 5 (directly) and F 6 (transitively) are included in F 7 . As a consequence, there are some topological scenarios (i.e., G ∈ F7 \F 5 ) for which A 2 has no chance of success, while A 3 has some. Such observation allows to claim that A 3 is more general than A 2 with respect to its topological requirements. Hence, two algorithms can be fairly (and formally) compared on the basis of their topological requirements. In the particular case of these two enumeration algorithms, however, the claim could be balanced by the fact that a sufficient condition is known for A 2 , while no one is known for A 3 . The choice for the right algorithm may thus depend on the target mobility context: if this context is expected to induce topological scenarios in F 5 or F 6 , then A 2 could be preferred, otherwise A 3 should be considered. More generally, it is however important to realize that a large gap may exist between necessary and sufficient topology-related conditions, and other topological properties (resp. evolving graph classes) could offer intermediate probabilities of success, which was not investigated for the given algorithms in this initial work.

Checking network traces for inclusion in the classes

We consider here the problem of checking automatically whether a given evolving graph belongs to one of the classes listed before. While having potentially a large scope of applications, this could allow in particular to help decide which algorithm is relevant to a given mobility context, by checking how the corresponding topological traces distribute over the classes. Below is a sketch of solution for each class met so far. The point is that all solutions can rely on common static graph properties, provided a few transformations. The transitive closure of an evolving graph G is the graph H = (V, AH ), where AH = {(vi, vj) :

∃J (v i ,v j ) ⊆ G)}.
A transitive closure is by nature a directed graph, as illustrated in Figure 6, since journeys are oriented entities. As explained in [START_REF] Bhadra | Complexity of connected components in evolving graphs and the computation of multicast trees in dynamic networks[END_REF], the computation of transitive closures can be done efficiently (in O(|VG|.|EG|.(log|S T |.log|VG|)), by building the tree of shortest journeys for each node in the network. We extend this notion to the case of strict journeys, with Hstrict = (V, AH strict ), where Given an evolving graph G, its underlying graph G, its transitive closure H, and the transitive closure of its strict journeys H strict , the inclusion of G in each class can be checked as follows:

AH strict = {(vi, vj) : ∃J strict(v i ,v j ) ⊆ G)}.
-G ∈ F1 ⇐⇒ H contains an out-dominating set of size 1.

-G ∈ F2 ⇐⇒ H is a complete graph.

-G ∈ F3 ⇐⇒ Hstrict contains an out-dominating set of size 1.

-G ∈ F4 ⇐⇒ Hstrict is a complete graph.

-G ∈ F5 ⇐⇒ G contains a dominating set of size 1.

-G ∈ F6 ⇐⇒ G is a complete graph.

-G ∈ F7 ⇐⇒ H contains an in-dominating set of size 1.

We expect most of the future classes to be possibly checked with similar approaches. This is however not a certainty.

Conclusion

This paper introduced a set of tools and methods dedicated to the analysis of distributed algorithms in dynamic networks. This new framework allows to characterize assumptions that a given algorithm requires in terms of topological evolution during its execution. It was illustrated by the analysis of three basic algorithms, and the analysis results were used to highlight potential implications of this work, including the possibility to compare algorithms on the basis of their topological requirements, and a sketch of classification of dynamic networks according to the corresponding properties. The problem of checking whether a given evolving graph belongs to the introduced classes was finally discussed.

Analyzing the requirement of algorithms is not a novel approach. It appears however that no proper transposition was previously done in the context of dynamic networks, where the usual practice is to liken dynamic topologies to static graphs. This is particularly striking in the recent field of population protocols [START_REF] Angluin | The computational power of population protocols[END_REF], where a common assumption is that all pairs of nodes interact repeatedly. In the light of the classification shown is this paper, such scenarios actually represent a subset of the most specific class among those discussed (namely, F 6 ). We think the framework proposed here could help characterize weaker assumptions for most population protocols.

The algorithms studied in this paper are simple. An interesting question for further research is whether the framework will scale to more complex algorithms, which remains unclear at this stage. We hope it could suit the study of common problems such as electing, naming, or building spanning structures (note that electing and naming may not have identical assumptions in a dynamic context). Another prospect is to investigate how intermediate properties could be explored between necessary and sufficient conditions, for example to guarantee a desired probability of success. Finally, as more properties are characterized and the classification grows, new insights may follow in the study of mobility models, based on checking generated traces for inclusion in the classes. Ultimately, this could answer questions like what kind of problems can be solved within a given mobility model, such as the well-known random way point model [START_REF] Bettstetter | The node distribution of the random waypoint mobility model for wireless ad hoc networks[END_REF], or in more realistic pedestrian and vehicular contexts.

Fig. 1 .

 1 Fig. 1. A basic dynamic scenario, where a node (b) moves during the execution.

Fig. 2 .

 2 Fig.2. Different powers of local computations; the scope of preconditions is depicted in white (on left sides), while the scope of actions is depicted in black (on right sides). The dashed elements represent entities (vertices or edges) that are considered by preconditions but remain unaffected by actions. The reader is referred to[START_REF] Chalopin | Local computations in graphs: The case of cellular edge local computations[END_REF] for a comparative study of some of these models.

Remark 1 .

 1 initial states: {I, N } (I for the initial emitter, N for all other vertices) alphabet: {I, N } preconditions(r1): λ(v0) = I ∧ λ(v1) = N actions(r1): λ(v1) := I graphical notation : Although the three algorithm examples provided in this paper consider pairwise interactions (and more specifically the model depicted on Figure 2(c)), the concepts and discussions developed in Sections 3 and 5 are not dedicated to it. Note that models such as those of Fig. 2(b) and 2(b) reflect well a wireless computing environment where nodes update their states according to those of their neighbors.

period t0 → t1 period t1 → t2 period t2 → t3 period t3 → t4Fig. 3 .

 3 Fig. 3. Example of an evolving graph covering a period of time from t0 to t4

Fig. 4 .

 4 Fig. 4. Graph Relabellings and Evolving Graphs -Combined formalism.

  initial states: {(C, 1),

Fig. 6 .

 6 Fig. 6. Example of transitive closure of an evolving graph

Detecting such a final state is not part of the given algorithm. The reader interested in termination detection as a distributed problem is referred to[START_REF] Godard | Termination detection of distributed algorithms by graph relabelling systems[END_REF].

⋆ Partially supported by A.N.R. grant No ANR-05-SSIA-0002-01. ⋆⋆ A.Ferreira is currently on leave as Head of Science Operations at the COST Office, Brussels, BE. http://www.cost.esf.org.


Thanh-Tung Dang

Michael Spence

Suzanne K Thomas

Matthew Carnie

Thanh-Tuân Bui
email: tbui@cyu.fr

T.-T Bui

Tetraphenylpyrrolo[3,2-b]pyrroles: Synthesis, Thermal, Optical, Electrochemical Properties and Photovoltaic Applications

Keywords: tetraarylpyrrolo[3,2-b]pyrrole, tetraphenylpyrrolo[3,2-b]pyrrole, molecular glass, hole transporting material, perovskite solar cell

We reported two tetraphenylpyrrolo [3,2-b]pyrroles in which the central core is functionalized by electron-rich functions: methoxy and di(4-methoxyphenyl)aminyl. The targeted compounds have good solubility in common organic solvents, high thermal stability and glass forming capacity. The optical and electrochemical properties of these compounds suggest use in a diverse range of applications. When utilized as the hole transporting layer of the planar perovskite solar cells, power conversion efficiencies ranging from 12.8 to 14.5% were obtained.

Introduction

In the recent years there has been a huge interest in designing small molecules for organic optoelectronic applications such as organic photovoltaics, 1 photodetectors, 2 organic lightemitting diodes, 3 etc. Each application requires molecules possessing specifically suitable thermal, morphological, optical, electrochemical, and electrical properties and processabilities. In the last decade, photovoltaic solar cells based on halide perovskite materials has been emerged as one of the most promising renewable energy conversion technology. Structurally, the perovskite solar cells (PSC) are constructed around a halide perovskite thin layer as light harvester. This active layer is sandwiched between two charge transporting layers of different nature: hole and electron transportation. Finally two current collectors complete the device structure. Hole transporting materials (HTM) extract and transport hole while on the opposite site electron is extracted and transported through the electron transporting materials (ETM). In order to ensure the good operation of the device, both ETM and HTM must have frontier orbital energy levels compatible with the valance band and conduction band energy of the perovskite materials. The ability to form a perfect charge transporting material/perovskite interface is important to have efficient charge extraction. In addition, high charge transporting capacities are of highly desired. The charge transporting materials can also protect perovskite materials from moisture and slow down its degradation by water vapor presented in the surrounding environment, thus effectively increase the device operation stability. 4,5 The HTM can be classified into two types which are inorganic and organic materials. Comparing to satisfied inorganic HTM in PSCs, 6 organic HTM show some advances in large-scale commercial production thanks to the unlimited molecular design possibility to tune the characteristics of targeted materials (HOMO LUMO energy levels, hole mobility, solubility, thermal stability, etc.). Notably, the availability of highly scalable synthetic procedures for these targeted molecular scaffolds are primordial in order to envisage a large-scale application. In this context, 1,2,4,5-tetraarylpyrrolo[3,2b]pyrroles (TAPP) are very promising candidates. 7 They can be easily prepared in one step by a straightforward multicomponent reaction of widely commercially available starting materials: an aromatic aldehyde, an aromatic amine, and butane-2,3-dione. 8 This synthetic procedure can be scaled for preparing the TAPP in multigram quantities and isolating them by simple filtration. 9,10 The electrochemical, electrical, and optical properties of these derivatives can be easily tuned by molecular structural design. 7,11 They can be found in diverse applications such as nonlinear optics, 12 fluorescent linkers in metal-organic frameworks, 13 light absorbing molecules in dye-sensitized solar cells, 14 electron donating molecules for solution-processed organic solar cells, 15 hole transporting molecules for perovskite photovoltaics. 16 In this studies, we designed, synthetized and characterized TPPP. The chemical structure of targeted molecules, named PP-2TPA and PP-4TPA, are shown in Figure 1. In these molecules, the TPPP central core is substituted by electron-rich functionalized groups: methoxy and di(4-methoxyphenyl)aminyl. We have commonly used these groups in designing novel molecules for diverse emerging technological applications such as hole transporting materials for hybrid photovoltaics, solution-processed electroluminescent devices, or as photoinitiators for photopolymerization and 3D printing. In addition, methoxy chains favor intermolecular interaction in solid state and increase charge transporting properties. 17

Results and Discussion

Materials Synthesis

Thermal properties

In comparison to polymeric materials, small molecules have features such as facile purification, well-defined structure and better solubility in common organic solvents, making them suitable for both solution-processing and vacuum processing techniques. However, small molecules tend to recrystallize and aggregate in thin film, thus a stable amorphous state is highly desired. Both compounds have satisfying thermal stability as revealed from ATG analysis (Figure S1 in ESI). Their thermal decomposition temperatures, determined at offset weight loss, are higher than 400°C, which is high enough for optoelectronic applications (T d (PP-2TPA) = 420°C and T d (PP-4TPA) = 418°C). Their glass forming properties are revealed by DSC analysis (Figure S2 in ESI). During the first heating cycle of the asprepared molecules, both compounds show melting transition. After cooling down from the melted states, both compounds formed a glassy state as evidenced from the subsequent heating cycles. As anticipated, thanks to the tridimensional structure of the central core, both PP-2TPA and PP-4TPA possess relatively high glass transition temperatures (T g = 134°C and 122°C, for PP-2TPA and PP-4TPA respectively).

Optical properties

The optical properties of these compounds have been investigated by absorption and emission spectroscopy in dilute Toluene solution (Figure 3) and their optical properties are gathered in Table 1. As can be observed, both compounds absorb and emit visible light in the same spectral range. In comparison with the tetraphenyl-1,4-dihydropyrrolo[3,2-b]pyrrole backbone, both absorption and emission of 2TPA and PP-4TPA have similar behaviors with important redshift. This is ascribed to the electron donating effect of periphery di(4methoxyphenyl)aminyl groups. 18 The nature of the substituent at the N-atom of the pyrrolopyrrole core has no influence to their corresponding optical characteristics. The The electrochemical properties of these compounds have been investigated in 0.1M tetrabutylammonium hexylfluorophosphate (n-Bu 4 NPF 6 ) in dichloromethane solution by using cyclic (CV) and differential pulse (DPV) voltammetries. Their CVs are given in Figure 4 and redox properties are listed in Table 1. The CVs of these molecules show multiple reversible oxidation processes. The first two oxidation processes are ascribed to the triphenylamine functions while the third at higher potential is issued from the pyrrolo[3,2b]pyrrole core. However, no reduction was observed. From the optical bandgap and the HOMO energy levels, the E LUMO can be calculated and reported in Table 1.

Potential application in photovoltaic devices

In a recent work we have reported a triphenylamine-functionalized thieno[3,2-b]thiophene derivatives as effective hole transporting materials for perovskite based solar cells (PSC). 20,21 In comparison with the published molecule, 2,5-bis(para-di(4methoxyphenyl)aminylphenyl)thieno[3,2-b]thiophene (TT-2,5-TPA), the thermal, optical and electrochemical of compounds A-D are similar to that of TT-2,5-TPA (E HOMO = -5. 13 eV, E LUMO = -2.48 eV). Encouraged by the success of TT-2,5-TPA in PSC, we have performed preliminary investigation of PP-2TPA and PP-4TPA as charge transporting molecules in hybrid photovoltaic cells. In this study, these two compounds were tested in planar PSC with device structure of Glass-ITO / ETM (SnO 2 ) / Perovskite (CH 3 NH 3 PbI 3 ) / HTM (molecule PP-2TPA or PP-4TPA) / Au (Figure 5). The PSC device made from Spiro-OMeTAD HTM were used as control devices. The photocurrent density-voltage curves of the best PSC devices are gathered in Figure 5. The device made from Spiro-OMeTAD delivered a photovoltaic performance of 18.98 % with open circuit voltage (V oc ) = 1.07 V, short-circuit current density (J sc ) = 22.10 mA.cm -2 and fill factor (FF) = 0.80. The concentration of the HTM solutions were optimized. Under these optimized conditions, the PSCs utilizing TT-2TPA exhibited the photovoltaic performance with the efficiency of 12.81%. This low efficiency was resulted from low J sc (20.72 mA.cm -2 ), V oc (0.95 V) and FF (0.65). However, replacing TT-2TPA in the hole transporting layer by TT-4TPA provided devices within increased power conversion efficiency due to improved V oc and FF. The average value of device performances are also given in Figure 5 with PCE Spiro- OMeTAD > PCE PP-4TPA > PCE PP-2TPA . Despite having similar HOMO and LUMO energy levels, TT-4TPA has superior photovoltaic property compared to TT-2TPA. Even further investigations have to be realized, we ascribe this difference to a better solubility of TT-4TPA leading to a more homogenous HTM layer and a better HTM/perovskite interface (Figure S3 in ESI). The presence of four triphenylamine moieties in TT-4TPA creates a bulky tridimensional structure and enhances its solubility and thin film morphology. 22

Conclusion

To summarize, two novel tetraphenylpyrrolo[3,2-b]pyrroles have been designed and synthesized. These molecules have been characterized in order to investigate their thermal, optical and electrochemical properties. They have high thermal stability with glass forming properties. They absorb in the visible/near UV spectral domain. Their frontier molecular orbital energies are calculated by a combination of experimental optical and electrochemical characterizations. These energy levels are found to be suitable to be employed as hole transporters in hybrid perovskite solar cells. Additional characterization must be performed to explore the full potential of these molecules, however the preliminary photovoltaic data are very promising.

Experimental

The detailed descriptions of general conditions for materials synthesis, optical, electrochemical and thermal characterization, perovskite-based photovoltaic device fabrication and measurements are given in the Electronic Supporting Information. The synthesis of 1,2,4,5-Tetrakis(4-bromophenyl)-1,4-dihydropyrrolo[3,2-b]pyrrole (1b) has been reported elsewhere. 23 Synthesis of 2,5-bis(4-bromophenyl)-1,4-bis[4-[(methoxy)oxy]phenylpyrrolo[3,2-b]pyrrole (1a): In a 50 mL Schlenk flask equipped with a reflux condenser and magnetic stir bar, glacial acetic acid (30 mL) was placed followed by addition of p-methoxyaniline (2 eq., 2.463 g), 4bromobenzaldehyde (2 eq., 3.7 g) and TsOH (0.1 eq, 172 mg). The mixture was stirred at 90°C for 30-45 min. After that time, butane-2,3-dione (1 eq., 0.87 mL) was slowly added via syringe and the resulting mixture was stirred at 90 °C for 4 h. The mixture was then cooled to room temperature. Synthesis of PP-2TPA: To a Schlenk flask was charged 1a (1eq., 751.2 mg), di(pmethoxyphenyl)amine (2.5 eq., 687 mg), Pd(OAc) 2 (0.2 eq., 53.76 mg), NaOtBu (10 eq., 1.15 g) and put in vacuum line for 5 mins then fill argon back (repeat 3 times). The mixture of starting materials was added toluene (40 mL) and PtBu 3 (0.4 eq., 0.48 mL) then heat 95 °C for overnight. After cooling down to room temperature, the solution reaction was filter by funnel with celite and washed the funnel by CHCl 3 . The organic solvent was removed under high pressure rotavapor. The crude product was dissolve in minimum of hot toluene with high rotation of stir bar, then large amount of MeOH was added. The precipitate of pyrrolopyrrole product was formed and isolated by funnel. The yield is 0.75 g (67.6 %). 1 H NMR (DMSO-d6, 250Hz): 7.2 (d, J = 7.7 Hz, 4H), 6.9-6.8 (broad, 24H), 6.6 (d, J = 7.9 Hz, 4H), 6.1 (s, 2H pyr ), 3.75 (s, 6H), 3.7 (s, 12H). HRMS (ESI+): Calculated for [M] •+ : 924.3881, found: 924.3886

Synthesis of PP-4TPA: To a Schlenk flask was charged 1b (1eq., 0.722 g), di (pmethoxyphenyl)amine (2.5 eq., 0.572 g), Pd(OAc) 2 (0.2 eq., 44.8 mg), NaOtBu (10 eq., 0.96 g) and put in vacuum line for 5 mins then fill argon back (repeat 3 times). The mixture of starting materials was added toluene (40 mL) and PtBu 3 (0.4 eq., 0.4 mL) then heat 95 °C for overnight. After cooling down to room temperature, the solution reaction was filter by funnel with celite and washed the funnel by CHCl 3 . The organic solvent was removed under high pressure rotavapor. The crude product was dissolve in minimum of hot toluene with high rotation of stir bar, then large amount of CH 3 OH was added. The precipitate of pyrrolopyrrole product was formed and isolated by funnel. The yield is 0.85 g (65 %). 1 H NMR (DMSO-d6, 250 Hz): 7.0-8.6 (multi, 48H), 6.1 (broad, 2H pyr ), 3.7 (s, 24H). HRMS (ESI+): Calculated for [M+H] + : 1319.5641, found: 1319.5643

Acknowledgement

The works at CY Cergy Paris Université is funded by French National Research Agency (Agence Nationale de la Recherche, ANR) through two projects: ESPOIR2 ANR JCJC (ANR-17-CE05-0006) and FUTUR -CY Initiative of Excellence (ANR-16-IDEX-0008).

Figure 1 .

 1 Figure 1. Molecular structures of tetraphenylpyrrolo[3,2-b]pyrroles (PP-2TPA and PP-4TPA)

Figure 2

 2 Figure2outlines the synthesis of PP-2TPA and PP-4TPA in a two-step synthetic procedure. The formation of the central TPPP core has been realized by following well described classical methods, as mentioned above, leading to intermediate 1a and 1b by one step condensation reaction of three components under acetic acid as catalyst. Mechanistically, para-bromobenzaldehyde firstly reacts with amine to form Schiff base. In the second step, cyclocondensation of Schiff base with 2,3-butanedione followed by an oxidation providing the intermediate pyrrolopyrrole (1a and 1b). The advances of this reaction are possible for low-cost starting materials, gram-scale and reproducibility. Notably, the work-up is easy by filtration, washing, and then recrystallization. Purification by column chromatography is not needed in this step. The presence of bromine atoms in their structures offer the possibility of further functionalization by Pd-catalyzed C-N coupling reactions. In the next reaction, twofold and fourfold Buchwald-Hartwig coupling reactions providing desired compounds in satisfying yields (> 65%). The isolation and purification of these products (PP-2TPA and PP-4TPA) based on the different solubility of product and starting materials in MeOH and Toluene. This synthetic strategy is highly efficient due to the gram-scale reaction and easy work-up.

Figure 2 .

 2 Figure 2. Synthesis of PP-2TPA and PP-4TPA

  optical bandgap ( ) of PP-2TPA and PP-4TPA (documented in Table 1) were calculated from their corresponding absorption onset in solution according the equation: = hc/ in which h and c are the Planck constant and the speed of light in vacuum, respectively. Both compounds have close optical bandgaps.

Figure 3 .

 3 Figure 3. Absorption and Fluorescent Spectra of PP-2TPA and PP-4TPA

Figure 4 .

 4 Figure 4. Cyclic and differential pulse voltammetry curves of PP-2TPA and PP-4TPA

Figure 5 .

 5 Figure 5. Photovoltaic device structure, J-V output and photovoltaic characteristics of the champion devices of Spiro-OMeTAD, PP-2TPA and PP-4TPA (The average values and standard deviations (calculated from 8 cells for Spiro-OMeTAD and PP-4TPA and 6 cells for PP-2TPA) are given in bracket)

Table 1 .

 1 Thermal and optoelectrochemical properties of PP-2TPA and PP-4TPA

	Compound	T d /T g	/			/		E HOMO	E LUMO
		°C	nm	nm	eV	V/Fc/Fc	+	eV	eV
	PP-2TPA	213/134	386/449	462	2.76	-0.22/0.05/0.66	-4.88	-2.12
	PP-4TPA	416/125	388/442	462, 482	2.81	-0.25/0.02/0.65	-4.85	-2.04

Data availability statement:

The data of the current study are available from the corresponding author on reasonable request.

Author Contribution: Materials syntheses were formed by Thanh-Tung Dang. Thermal, Optical, and Electrochemical data were collected by Thanh-Tuân Bui. Photovoltaic data collection and analysis were performed by Michael Spence and Suzanne K. Thomas. Thanh-Tuân Bui and Matthew Carnie designed and supervised the study. The first draft of the manuscript was written by Thanh-Tuan Bui and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Conflicts of interest