Wuyang Zhong

Patrick Perre

Fanny Duval

Julien Lemaire
email: julien.lemaire@centralesupelec.fr

A predictive model of organic acids separation by chromatography with strong anionic resins in sulfate form

Keywords: Chromatography, Modeling, Strong anionic resin, CE/SE method, Langmuir adsorption, ion-exchange

Organic acids commonly have quite symmetrical chromatography profiles at low pH (< 1.5) with strong anionic resins, but a significant tailing can be observed with succinic and citric acids. Classical adsorption models, like the Langmuir model, fail to predict this behavior, which can have a major influence on mean retention times and profile shapes, therefore on chromatography performances.

A new retention model was developed to better predict organic acid separation with strong anionic resin. This model combines a refined Langmuir adsorption model and an ionexchange model. Organic acid adsorption is assumed to be due to hydrogen bonding with sulfate and hydrogen sulfate counter-anions on the resin. The adsorption capacity depends mostly on molecular size: up to sixteen formic acid molecules could be adsorbed per counter-anions, meanwhile only two succinic acid or one citric acid molecules could be adsorbed.

This adsorption model was then embedded in a generic and accurate modeling approach (continuous column with mass balance equations solved by the conservation element/solution element (CE/SE) method). All parameters of this column model were identified by fitting the simulation to experimental results (equilibrium curves and pulse tests). Then, the column model was validated with original experimental results from a binary mixture pulse test (formic and succinic acids). Results show that the simulations are much more predictive for multi-component pulse tests, both in terms of profile shape and retention time, which cannot be captured without considering ion-exchange.

Introduction

Chromatographic separation is a well-established technology in many domains, such as the food and feed industry [START_REF] Jun | Method for detecting paraffin in food by gas chromatography and gas chromatography-mass spectrometry and application thereof, Worldwide applications[END_REF], pharmaceutical production [START_REF] Shane | Methods for monitoring integrated continuous pharmaceutical manufacturing processes[END_REF], protein and enzyme purification [START_REF] Qiao-Le | Model-based process design of a ternary protein separation using multi-step gradient ion-exchange SMB chromatography[END_REF], and heavy metal separation [START_REF] Ammann | Speciation of heavy metals in environmental water by ion chromatography coupled to ICP-MS[END_REF]. Chromatography can separate or purify target components, even with very close properties like stereoisomers, without using a large amount of energy, chemicals, or solvent. For example, chromatography is often used to separate sugars or organic acids from hemicellulose hydrolysates or fermentation broths [START_REF] Kaifei | Ion exchange separation for recovery of monosaccharides, organic acids and phenolic compounds from hydrolysates of lignocellulosic biomass[END_REF]. These applications often require high purity, high recovery rate, and low operational costs for economic reasons. Thus, chromatography is the ultimate solution when conventional separation techniques such as distillation, liquid-liquid extraction, membrane filtration, or crystallization are not sufficiently selective or consume too much energy, chemicals, or solvent.

Preparative chromatography consists of pushing a mixture (generally liquid) through a column filled with a porous solid phase using a mobile phase called the eluent. Its principle is to separate different compounds according to their affinities towards both phases that affect their travel speed and retention time inside the column. The solid phase can be composed of dense or microporous resin beads. Hundreds of resin types exist on the market to separate ions or organic molecules, depending on their structures and retention mechanisms. Retention of a molecule on a resin surface often implies adsorption by Van Der Waals forces owing to its hydrophobic matrix.

Ion-exchange resins are covered with counter-anions or counter-cations, making them hydrophilic [START_REF] Bolin | Preparation of hydrophobic interaction chromatographic packings based on monodisperse poly (glycidyl-methacrylate-coethylenedimethacrylate) beads and their application[END_REF]. For instance, molecules such as sugars or organic acids can be separated by different adsorption mechanisms (H bond, van der Waals force, chelation) with the counterion fixed on the resin surface, and ionic compounds can be separated by ion-exchange or Donnan exclusion effect. For example, strong anionic resins (coated with quaternary ammonium groups) or weak anionic resins (coated with amine groups) are often used to separate organic acids by interacting with their acid group.

However, the application of chromatography at an industrial scale is severely limited by its low productivity, which impacts capital expenditures, and by how it dilutes separated compounds, which impacts both capital and operational costs because further concentration steps are required. To tackle these limiting factors, several dozen multi-column chromatography processes have been developed over the last three decades, but their development and tuning are difficult. Modeling the elution profile inside columns is a valuable tool to make predictive simulations and optimize the performance of multi-column chromatography systems. Further investigations and model development are still needed to improve prediction accuracy, evidenced by continuing research activity in this field [START_REF] Shanshan | Preparative separation of five polyphenols from the fruits of Sorbus pohuashanensis Hedl. by high-speed counter-current chromatography[END_REF][START_REF] Ton | Charge-based separation of synthetic macromolecules by non-aqueous ion exchange chromatography[END_REF][START_REF] Ali | Mathematical modeling and optimization of industrial scale ELUXYL simulated moving bed (SMB)[END_REF][START_REF] Cheol-Yeon | Development of a simulated moving bed process for ultra-high-purity separation of ribose from a low-selectivity sugar mixture in microalgal hydrolyzate[END_REF].

For instance, our previous work found that outlet profiles of some organic acids have uncommon tailing that cannot be explained by usual Langmuir adsorption models [START_REF] Julien | Purification of organic acids by chromatography with strong anionic resins: Investigation of uptake mechanisms[END_REF]. This phenomenon was explained theoretically by considering the small fraction of organic acid dissociated into anions (< 1% at pH below 1.5). Historically, this fraction was neglected, but it can be strongly retained on resin ion-exchange sites [START_REF] Julien | Purification of organic acids by chromatography with strong anionic resins: Investigation of uptake mechanisms[END_REF]. A new model was then proposed to better predict the behavior of organic acids with strong anionic resins. In addition, Langmuir adsorption mechanisms were described in more detail, based on Honeywell Universal Oil Products (Honeywell UOP) patent assumptions [START_REF] Kulprathipanja | Separation of lactic acid from fermentation broth with an anionic polymeric absorbent[END_REF]. Mono-component simulation of an organic acid profile was performed using a discontinuous column model (Martin and Synge plate model [START_REF] Guiochon | Fundamentals of preparative and non-linear chromatography[END_REF]) and a finite difference method to solve differential equations. Uncommon tailing of succinic and citric acid profiles was better predicted. This model still needs to be adjusted and then validated with more experimental data and improved for multi-component simulation.

Chromatography simulation is based on solving mass balances in the mobile phase (convection-diffusion equations) and local exchanges with the stationary phase. It generally results in a non-linear and coupled system of partial differential equations. The classical numerical solving methods are finite element [START_REF] Juke | Optimal design of batch and simulated moving bed chromatographic separation processes[END_REF][START_REF] Kaczmarski | Use of orthogonal collocation on finite elements with moving boundaries in the simulation of non-linear multi-component chromatography. Influence of fluid velocity variation on retention time in LC and HPLC[END_REF][START_REF] Hørsholt | A Discontinuous-Galerkin Finite-element method for simulation of packed bed chromatographic processes[END_REF], finite difference [START_REF] Dünnebier | Computationally efficient dynamic modelling and simulation of simulated moving bed chromatographic processes with linear isotherms[END_REF][START_REF] Kaczmarski | Fast finite difference method for solving multi-component adsorption-chromatography models[END_REF], or finite volume [START_REF] Lim | Accuracy, temporal performance and stability comparisons of discretization methods for the solution of partial differential equations (PDEs) in the presence of steep moving fronts[END_REF] for the spatial discretization, with an explicit or implicit time scheme. The CE/SE (conservation element/solution element) method, proposed by NASA in 1995 [START_REF] Chang | The method of space-time conservation element and solution element-A new approach for solving the navier-stokes and euler equations[END_REF], has been applied for simulating explosions [START_REF] Wei | Study on the mechanism of the deflagration to detonation transition process of explosive[END_REF][START_REF] Shen | Numerical simulation of liquid-fueled detonations by an eulerian-lagrangian model[END_REF], crystallization [START_REF] Orkomi | Simulation and control of multidimensional crystallization processes[END_REF] and chemical reactions [START_REF] Qamar | The space-time CESE method for solving special relativistic hydrodynamic equations[END_REF][START_REF] Sadighi | Dynamic simulation of a pilot scale vacuum gas oil hydrocracking unit by the space-time CE/SE method[END_REF]. It was first used to simulate chromatography with the Langmuir adsorption model in 2004 [START_REF] Young-Il | A fast and accurate numerical method for solving simulated moving bed (SMB) chromatographic separation problems[END_REF]. The CE/SE method allows simultaneous treatment of space and time integration. In the case of packed bed chromatography [START_REF] Lan | Performance evaluation of the conservation element and solution element method in SMB process simulation[END_REF], global trials of the CE/SE method showed a better trade-off between accuracy, stability, and calculation speed than traditional methods using an upwind scheme. Compared to conventional approaches, CE/SE method was proved, at equivalent accuracy, to save the CPU time by 77% for the pseudo-linear system of glucosefructose separation, and by 96% for non-linear enantiomer separation [START_REF] Chuanyi | Combination of spacetime conservation element/solution element method and continuous prediction technique for accelerated simulation of simulated moving bed chromatography[END_REF].

The present work aims to improve the prediction of the outlet profile of different organic acids and to confirm the role of ion-exchange with strong anionic resin. A more detailed Langmuir adsorption model is proposed based on our previous theoretical work [START_REF] Julien | Purification of organic acids by chromatography with strong anionic resins: Investigation of uptake mechanisms[END_REF], assuming that adsorption is due to a hydrogen bond between organic acids and sulfate anions, as mentioned first in a UOP patent [START_REF] Kulprathipanja | Separation of lactic acid from fermentation broth with an anionic polymeric absorbent[END_REF]. Experiments were performed first to identify all model parameters (equilibrium constant, mass transfer coefficients, axial dispersion coefficient), then to validate assumptions, and finally to check model prediction. In addition, a continuous column model, solved by the CE/SE method, was developed to more accurately simulate multi-component separation by chromatography. This new model tackles discrepancies observed in the case of succinic or citric acid separation with strong anionic resin.

Materials and methods

Experimental setup

The separation of organic acids by chromatography with strong anionic resin was studied using a classical chromatography setup. Experimental devices and chemicals used in this study are detailed in Tables 1 and2.

The product or the eluent is fed at the top of column by the volumetric pump. The liquid phase flows downwards in the column, while soluble compounds interact with the resin. The output liquid goes through a pH and conductivity sensor, then a UV absorbance analyzer, and finally is sampled by a fraction collector for offline analysis. All sensors and the fraction collector are connected to a computer to monitor pH, conductivity, UV absorbance, and sample collection.

The eluent was a H 2 SO 4 solution at pH 1.5 (≈ 0.02 mol/L) for all experiments. For information, its ionic strength is 0.0314 mol/L. Product pH was also adjusted at 1.5 to limit pH and conductivity variation, and to keep a similar ionic strength during coumpound elution. The on-line measurement of UV absorbance was used to select samples to analyze by HPLC. The HPLC system (Ultimate 3000, Dionex) is equipped with a refractometer RI-101 (Shodex) at 35°C. The analytical column is an Aminex HPX-87H (Biorad) connected to a guard column (Micro-Guard cartridge Biorad) kept at 45°C within an oven. The mobile phase is a 2 mmol/L H 2 SO 4 solution, prepared with Milli-Q water (Direct 8, Millipore).

The flowrate was set at 0.5 BV/h (BV represents the total bed volume, equal to 0.35 L), corresponding to 3.10 -4 m/s interstitial flow velocity (Re p = 6.6). The internal column diameter was 2.5 cm, and the bed length was 70 cm. A jacketed column linked to a thermostatic bath was used to maintain the temperature at 25°C. The column was filled with the strong anion resin DIAION UMA 150, functionalized with trimethylamine groups (type 1), manufactured by Mitsubishi Chemical, whose characteristics are summarized in Table 3. The resin is supplied in Cl -form, but it was put in sulfate form before running experiments by using the following procedure: 2 BV of 1 mol/L Na 2 SO 4 solution was fed into the column at 1 BV/h, then 8 BV of eluent was fed to wash the bed resin and reach pH equilibrium.

The bed porosity was measured by a pulse test with 10 mL 1 mol/L Na 2 SO 4 solution eluted at 1 BV/h. The porosity was estimated between 0.39 and 0.41, which corresponds to the range given by the resin manufacturer (Table 3). Finally, the extra-column volume (in tubes, valves, sensors, etc.) was measured by performing a similar pulse test without the chromatography column. It was estimated around 0.05 BV (15 mL). All chromatographic profiles were corrected by removing this extra-column volume from the eluted volume.

Determination of transfer and axial dispersion coefficient

A pulse test consists of injecting a small volume of product (0.05 BV) followed by the eluent at a constant flow rate until all compounds leave the column. It is one of the most common methods used in chromatography to study compound retention and dispersion along the column. A pulse test was chosen in our study to determine physical properties such as adsorption isotherm, mass transfer coefficients, bed porosity and extra-column volume.

During pulse tests, component concentrations at column output are plotted as a function of the eluted volume. The curve obtained is usually called a chromatographic profile, which looks like a peak more or less wide and symmetric. This profile is mainly characterized by its mean retention volume and its standard deviation :

where is the eluted volume (in BV), and is the component concentration in the mobile phase at column output (concentration unit does not matter).

The peak width depends almost on axial dispersion and mass transfer resistance while the component passes through the column. The volume of product fed into the column (< 0.05 BV) is small enough to have a negligible influence on the shape of the pulse test profile. Dead volumes were minimized to limit their influence on dispersion. The global dispersion is usually characterized by the height equivalent to a theoretical plate (HETP) as a function of , and the bed length :

The van Deemter curve [START_REF] Van Deemter | Longitudinal diffusion and resistance to mass transfer as causes of nonideality in chromatography[END_REF] describes the relationship between HETP (inversely proportional to the global dispersion) and the interstitial velocity (in m.s -1 ). This curve is estimated experimentally by plotting the HETP of component profiles obtained at different flowrates during the pulse test (Fig. 1).

One of the most known versions of the van Deemter equation reads as:

This equation considers that the global dispersion has three main causes: the mechanical dispersion ( term), the molecular diffusion ( term), and the mass transfer resistance ( term), as depicted in Figure 1. The mathematical expressions to obtain the values of and are important for prediction purposes. The original article by van Deemter [START_REF] Van Deemter | Longitudinal diffusion and resistance to mass transfer as causes of nonideality in chromatography[END_REF] gives the following expressions:

(

) 5
where is a packing factor that equals 1 for a homogeneous packed bed, is the mean particle diameter (in m), is the molecular diffusion coefficient (in m 2 .s -1 ), is the bed tortuosity, is the bed porosity, is the volumetric ratio of mobile phase to stationary phase, is the linear equilibrium constant (concentration ratio between stationary and mobile phase), and is the liquid-side local mass transfer coefficient (s -1 ).

Within narrow operating conditions of preparative chromatography (flowrate between 0.5 and 2 BV/h), van Deemter curves are almost linear over the whole intertitial velocity range. Thus, the most simple version of van Deemter equation was sufficient to describe experimental results and estimate its coefficients.

From experimental points, the range of can be estimated quite accurately by linear regression, then range can be deduced from previous equation. and are difficult to estimate separately because of the lack of experimental points in the non-linear part of van Deemter curve. However, in the present modelling work, only axial disperion coefficient is required, and it can be expressed as the sum of molecular diffusion and mechanical dispersion, according to the following equation: Thus, the range of can be estimated quite accurately from the experimental points by fitting the van Deemter equation with the least squares method, while considering the range of previously estimated by linear regression.

Physical formulation at the bed level 2.3.1 Mass balance

Two types of unidimensional column model are often used to describe the mass transfer inside a chromatography column [START_REF] Guiochon | Fundamentals of preparative and non-linear chromatography[END_REF]:

 Continuous models, solved by usual finite methods, where the column is discretized into a large number of very small elements. In this case, the axial dispersion depends on the axial dispersion coefficient and the global mass transfer coefficient .  Discontinuous models where the chromatography column is assimilated to a plate column. Each theoretical plate acts as a perfectly stirred reactor. In this case, the axial dispersion becomes an emerging property of the theoretical plate number N.

Discontinuous models, easy to code and to solve, require low CPU time. However, because the axial dispersion is an indirect effect of the theoretical plate number N, this model imposes the same axial dispersion for all compounds. Thus, this model is not appropriate for multi-component chromatography simulation if their axial dispersion coefficients are not similar. That is the reason why a continuous model was chosen in the present work.

At the level of the chromatography column (Fig. 2), the local 1-D mass balance of one component j in the mobile phase, for a section of extension , can be expressed as two mass balances for liquid side and solid side:

where is the concentration of component j in the mobile phase (in mol.m -3 ), is the axial dispersion coefficient (in m.s -2 ), is the interstitial velocity (in m.s -1 ), is the mass variation of component j due to acid-base reactions (in mol.m -3 .s -1 ), is the mass exchange of component j between the mobile phase and the stationary phase (in mol.m -3 .s -1 ), and is the concentration of component j in the stationary phase (in mol.m -3 ).

The

term is due to the acid-base reactions in the mobile phase between SO 4 2-and HSO [START_REF] Ammann | Speciation of heavy metals in environmental water by ion chromatography coupled to ICP-MS[END_REF] -anions, organic acid HA and its conjugate base A -, and H + cations. When the mobile phase goes through the column, organic acids HA can be retained by Langmuir adsorption; meanwhile, A -can be exchanged by SO 4 2-or HSO 4 -anions on resin sites. As the concentration ratio between acid and basic compounds changes in the mobile phase, the acid-base reactions tend to restore equilibrium by transforming the acid into a base or the inverse.

Finally, to avoid any term, mass balances in the mobile phase were written for the total content of each organic acid A ( ), the total sulfate content SO 4 (

) and the total acid content H ( ). These balances are not affected by acid-base reactions. Then, the concentration of each compound can be estimated using the equilibrium constant of each acid-base reaction. Indeed, these reactions are assumed to always be at an equilibrium state because they are considered to be much faster than all other phenomena. The mass exchange between the mobile phase and stationary phase depends on the adopted mass transfer model. With instantaneous mass transfer model, concentrations in the stationary phase are assumed to be in equilibrium with concentrations in the mobile phase. In this case, it is not necessary to estimate because concentrations can be deduced directly from concentrations by solving the equilibrium equations system. With a non-instantaneous mass transfer model, a new equation is needed to estimate , which also depends on equilibrium equations between the mobile and stationary phase.

In a previous work [START_REF] Julien | Purification of organic acids by chromatography with strong anionic resins: Investigation of uptake mechanisms[END_REF], the instantaneous mass transfer model was adopted with a discontinuous model, as the influence of mass transfer resistance can be considered directly in the theoretical plate number that summarizes the global effect of all dispersion mechanisms.

In the present work, a non-instantaneous mass transfer model was chosen to consider compounds with different mass transfer resistances. Mass transfer resistance between the mobile and stationary phase is driven by molecular diffusion and hydrodynamic properties at the interface (summarized in the Sherwood number, depending on Reynolds and Schmidt numbers). The lumped model was chosen to estimate the mass exchange from the global mass transfer coefficient of component j (in s -1 ) between the mobile phase and the stationary phase, due to diffusion inside resin and at its interface [START_REF] Felinger | Comparison of the kinetic models of linear chromatography[END_REF]:

where represents the theoretical concentration of component j in the stationary phase (in mol.m -3 ), which would be in equilibrium with the mobile phase. This equilibrium is calculated from the concentration of all components in the mobile phase by solving all equilibrium equations (adsorption, ion-exchange, and acid-base reactions).

The mass exchange can also be expressed at the liquid side, as:

where is the concentration of the component j at the interface (in mol.m -3 ).

Thus, can be estimated from the value by assuming a linear equilibrium, whose coefficient is estimated from the adsorption isotherm of component j (cf. section 2.3.2).

(11)

Adsorption mechanisms

Counter-anions SO 4 2-or HSO 4 -cover the surface of strong anionic resin in sulfate form (Fig. 3). Their proportion depends on ion-exchange equilibrium and their concentration in the mobile phase. The adsorption mechanism of organic acid on this type of resin is not well documented in the literature. UOP patents [START_REF] Kulprathipanja | Separation of lactic acid from fermentation broth with an anionic polymeric absorbent[END_REF][START_REF] Kulprathipanja | Separation of an organic acid from fermentation broth with an anionic polymeric adsorbent[END_REF] suggested that the molecular form of organic acid was retained only on counter-anions HSO 4

-and SO 4 2-. Lemaire et al. [START_REF] Julien | Purification of organic acids by chromatography with strong anionic resins: Investigation of uptake mechanisms[END_REF] developed this hypothesis in their model, assuming that the interaction was probably hydrogen bonding between the acid hydrogen atom of organic acids and the lone electron pairs of sulfate and hydrogen sulfate oxygen atoms.

A Langmuir model is generally used to describe adsorption on chromatography resins, which involves several assumptions [START_REF] Sircar | Comments on practical use of Langmuir gas adsorption isotherm model[END_REF]:

1. The surface containing adsorption sites is perfectly flat with no corrugation, 2. All adsorption sites are energetically equivalent, 3. Only mono-layer coverage occurs, 4. No interaction occurs between molecules adsorbed on adjacent sites.

As shown in Figure 3, each SO 4 2-or HSO 4 -contains eight lone electron pairs, so up to eight hydrogen bonds can be theoretically observed on each anion. As an improvement compared to previous works [START_REF] Julien | Purification of organic acids by chromatography with strong anionic resins: Investigation of uptake mechanisms[END_REF], this feature was directly integrated into the adsorption equation in the present work. As a first approximation, the adsorption sites on SO 4 2-and HSO 4 -were assumed to be similar. Thus, instead of writing two adsorption equations for SO [START_REF] Ammann | Speciation of heavy metals in environmental water by ion chromatography coupled to ICP-MS[END_REF] 2-and HSO 4 separately, the total concentration of SO 4 2-and HSO 4 -in the stationary phase, called (in mol/L), was used. In the case of a single component model, the Langmuir equilibrium model for component j was described as:

where is the number of lone electron pairs on each counter-anion , is the adsorption constant of component j, and is the number of lone electron pairs occupied by each adsorbed molecule.

In the case of a multi-component model, the competitive Langmuir equilibrium model for m components can be described as:

The value of can be estimated from the value of the ion-exchange resin capacity , which is the concentration of quaternary ammonium groups RMe 3 N + . The resin supplier provides an estimation (ca. 2.3 eq/L where eq corresponds to the equivalent molar charge amount). This value was checked by mixing a given volume of resin in OH -form with a given volume of H 2 SO 4 solution whose concentration is known. The titration of residual acid concentration in the solution enabled estimating the initial OH -concentration in resin (2.1 ± 0.2 eq/L).

Equilibrium curves of additional organic acids, determined by frontal analysis, are currently in progress in our laboratory. This dataset is intended to refine the competitive Langmuir equilibrium model. So far, these additional results are consistent and will be the object of another paper. Equilibrium curves enable us to estimate Langmuir model coefficients ( and ) by fitting the model to experimental points. As the and values are correlated in the linear zone of equilibrium curves, their estimation can be quite inaccurate when it is impossible to obtain sufficient experimental points in the non-linear zone. Consequently, and are given in Table 4 with their uncertainty range. Because represents the number of lone electron pairs occupied by each adsorbed molecule, was assumed afterward to be equal to only rounded values that could be explained physically. When the experimental range of can include two possibilities (for formic and succinic acids), two corresponding ranges of the k s value were identified.

According to the Langmuir adsorption model and estimation from Table 4, formic, succinic, and citric acids seem to occupy a different number of lone electron pairs on oxygen atoms of sulfate or hydrogen sulfate counter-anions. The values are rather consistent with their different molecular size and the number of carboxylic acid groups (Fig. 4).

Citric acid occupies more adsorption sites, most likely because of its large molecular size, which hinders access to free lone pairs, and its three carboxylic acid groups that can make H bonds with up to three lone pairs. = 8 corresponds to only one citric acid molecule per sulfate or hydrogen sulfate counter-anion (all eight lone pairs occupied: H-bond or cluttered). Inversely, = 0.5 or 1 for formic acid, which is the smallest acid and has only one carboxylic acid group. It corresponds to one or two formic acid molecules per lone pairs (none cluttered and only single H-bonds). The formation of formic acid dimers could explain that up to two molecules can occupy the same adsorption site [START_REF] Benni | Catalytic effect of water, water dimer, or formic acid on the tautomerization of nitroguanidine[END_REF]. Succinic acid is an intermediate case because of its medium size and two carboxylic acid groups. = 4 or 8 corresponds to one or two succinic acids per sulfate or hydrogen sulfate counter-anions (four or eight lone pairs occupied). values of formic, succinic, and citric acid are quite similar and correspond to the strength of the H-bond between their carboxylic acid group and lone electron pairs of sulfate or hydrogen sulfate counter-anions. Their value has been refined further in the present work to remove the possible influence of ion-exchange that may have caused overestimation.

Ion-exchange mechanisms

In the present model, ion-exchange was also considered to explain the retention of some organic acids with a particular chromatographic profile (large tailing) with strong anionic resin in sulfate form. Indeed, despite the very low dissociation of organic acids at pH below 2, a small concentration of corresponding organic anions could greatly influence their profile in the case of a very strong affinity for resin ion-exchange sites.

In the mobile phase, organic acid HA and hydrogen sulfate anion HSO [START_REF] Ammann | Speciation of heavy metals in environmental water by ion chromatography coupled to ICP-MS[END_REF] -can dissociate into its corresponding conjugate base A -or SO [START_REF] Ammann | Speciation of heavy metals in environmental water by ion chromatography coupled to ICP-MS[END_REF] 2-and one proton H + . The corresponding acid-base equilibria can be written as follow:

(14) (15)
In the present model, usual pKa values of formic, succinic, and citric acids and hydrogen sulfate at 20 °C (operating temperature) were considered: 3.76, 4.16, 3.13, and 1.99, respectively.

Organic anions

can replace or counter-anions on quaternary ammonium ion-exchange sites. Ion-exchange reactions correspond to the following equilibria:

(16) (17) (18)
where is the ion-exchange coefficient of anions X -by anions Y -on resin sites.

In the stationary phase, the concentration of ion-exchange sites is always conserved:

All corresponding equilibrium constants can be deduced from resin manufacturer data listed in the appendices [START_REF]Ion Exchange Resins Selectivity, DuPont[END_REF][START_REF]DOWEX™ Ion Exchange Resins Technical Information[END_REF], which provides the value of between many anions X - and OH -.

Computational strategy

As detailed in the introduction, several computational methods could be used to solve the system of equations in space and in time. The CE/SE method was chosen in the present work for its ability to deal with space and time integrals simultaneously. The integral is calculated between normal mesh points and intermediate mesh points, and it takes two integrals to move one step forward in space or time.

This method is based on Gauss's divergence theorem [START_REF] Sin-Chung | Courant Number and Mach Number Insensitive CE/SE Euler Solvers[END_REF]. In the case of a 1D space model, it is also called Green's theorem. For a partial differential equation (PDE) like:

Green's theorem tells us that the line integral along a closed curve of a function equals the integral of the derivative over the corresponding plane region :

Figure 5 shows the space and time discretization used and the mesh points where all variables are estimated. By definition, the solution element SE (j,n) is the gray space-time interior region centered at point A (j,n). It includes a horizontal line segment, a vertical line segment, and their immediate neighborhood points (B, D, F). In this solution element, we chose to approximate and values by a first-order Taylor formula from its center point A (j,n):

where and correspond to the space derivative of and at the mesh point (j,n), and and correspond to the time derivative of and at the mesh point (j,n).

A zero-order Taylor formula was chosen to approximate the source term , which means that was treated as a constant in each solution element:

The conservation element CE (j,n) corresponds to the rectangle FADE (blue) and ABCD (green), where the source term was consequently considered equal at mesh points A, B, D, and F, but different at mesh points C ( and E ( .

Finally, Green's theorem applied within the rectangle FBCE (blue + green) by approximating , and leads to the following explicit equation that enabled to be estimated from values of , and at the previous semi-time step: where In the case of chromatography convection-diffusion equation, the CE/SE method can be applied to solve the mass balance in the mobile phase by defining:

(27)
As a consequence, , , , then can be deduced from at the mesh point (j,n):

(Eq 20) (28)
This set of equations prove that only and are needed to solve Eq 25. There are many possibilities to calculate the derivative from . In the present model, it was calculated by using a modified harmonic mean [START_REF] Chuanyi | Combination of spacetime conservation element/solution element method and continuous prediction technique for accelerated simulation of simulated moving bed chromatography[END_REF]: [START_REF] Kulprathipanja | Separation of lactic acid from fermentation broth with an anionic polymeric absorbent[END_REF] where are the partial derivatives on the left and right sides, respectively. The modified harmonic mean equals zero when have opposite signs and equals the traditional harmonic mean when have the same sign.

Thus, only one variable is needed to solve the system, which is calculated via Eq 25 with an explicit scheme. Starting from points at time step n, Eqs 22-30 allow us to move to time step n+1/2. Points at time n+1 are then calculated in the same way as points at time n+1/2.

Meanwhile, the mass balance in the stationary phase can be solved separately by a 1 st order finite-difference method, at each mesh point (j,n):

Boundary conditions are used to obtain values of first and last mesh points, as proposed by Chuanyi et al. [START_REF] Chuanyi | Combination of spacetime conservation element/solution element method and continuous prediction technique for accelerated simulation of simulated moving bed chromatography[END_REF]:

Results and discussion

Influence of Langmuir and ion-exchange coefficients

A tailing can be observed in the chromatography outlet profile of some organic acids with a strong anionic exchange resin. Theoretically, this could be due to a high value of their adsorption affinity or ion-exchange coefficient . Figures 6-A and6-B show the impact of values and on the outlet profile during the single component pulse test simulation of citric acid with parameters of Tables 4, 5 and 7. In Figure 6-A, is fixed at 1, a negligible value, to study only the influence of on the profile shape and position. Then, the influence of is depicted in Figure 6-B. Pulse tests with different values of were simulated. As the objective was to highlight the effect of only on the profile shape, the value of was adjusted for each value of to obtain the same retention volume.

As expected, Figure 6-A shows that the higher is the adsorption coefficient , the higher is the retention volume. However, the value of the coefficient barely changes the profile shape and its relative symmetry. In comparison, Figure 6-B shows that significantly impacts the profile shape by increasing its tailing.

To summarize, the retention volume depends on both and while the profile tail is influenced almost solely by . Indeed, very high values could also lead theoretically to significant profile tailing, but these values would not be consistent with experimental retention volume, which corresponds to a ranging between 0 and 0.3 L/mol. The impact of ion-exchange constant depends on the pH of the liquid phase and the pKa of the organic acid. When the pH is much lower than pKa (pH < pKa -2), organic acids are mainly in molecular form. In this case, a very small proportion is dissociated into the anionic form, and the peak is less sensitive to the value of .

Identification of Langmuir and ion-exchange coefficients

In the present study, was estimated to reproduce the experimental profile tailing obtained for formic, succinic, and citric acids, starting from the resin supplier's estimation. Meanwhile, was adjusted to fit the experimental retention volume. Initially, was set to the experimental value estimated previously from equilibrium curves, considering different possible values (Table 4). Figure 7 shows experimental mono-component pulse tests that were used to refine estimations of and given in Table 6. The operating conditions of these pulse tests are detailed in Table 5.

The mean value of obtained from the van Deemter curve was used for simulation. Only the value of was slightly adjusted to fit better the peak left front (Table 7).

In Table 6, the identified ranges of correspond quite well to the supplier's data for formic and citric acids. ranges identified from pulse test simulations were not influenced by different possible values considered for each organic acid. Indeed, the denominator of the Langmuir equation remained close to 1 because of weak component concentration during pulse test simulations. For formic and succinic acids, the range estimated from pulse test simulation is consistent with its range identified from equilibrium curves, for = 0.5 and 4, respectively. Meanwhile, the range of citric acid identified from pulse tests lies at the bottom of its range estimated from equilibrium curves. Thus, 0.5, 4 and 8 were considered afterward as the unique, consistent values for formic, succinic and citric acids, respectively. For citric acid, the slight overestimation from equilibrium curves can be explained by the strong influence of ion-exchange (illustrated by the significant tailing of the citric acid profile, cf. Fig. 7). Both models can predict the experimental results for formic acid. This is not surprising as its tailing is negligible, and its profile is almost symmetrical. The formic acid profile is marginally influenced by ion-exchange due to its very low (< 30). Moreover, it is much less dissociated into anions than citric acid because of its higher pKa (3.76 vs. 3.13). Thus, at pH 1.5, only 0.5% formic acid is dissociated, making the ion-exchange effect even more negligible.

Analysis of ion-exchange effect on the chromatography profile

In comparison, succinic acid has the highest pKa (4.16) but also the highest according to Table 6 (900 ± 100). Therefore, despite a low proportion dissociated into anions (0.2% at pH 1.5), the strong retention of succinate anion on resin ion-exchange sites can explain its slight profile tailing.

Finally, the citric acid profile has the largest tailing because of its lower pKa (3.13) and quite high . Indeed, citric acid is the most dissociated into anions (2.3% at pH 1.5) and citrate anions are strongly retained by ion-exchange.

Figure 7 confirms that the tailing of succinic and citric acid profiles cannot be reproduced by the Langmuir adsorption model. In contrast, the hybrid model is in good agreement with experimental results. It confirms that organic anions retention by ion-exchange must be considered to make predictive simulation, even if they are poorly dissociated at pH = 1.5.

Multi-component pulse test prediction

To confirm that the hybrid model is more predictive than the model with only Langmuir adsorption, a pulse test with a mixture of formic acid and succinic acid was both simulated and carried out (Fig. 8). Table 8 gives parameter values used in both models, whose ranges were estimated previously from experimental data analysis (Tables 6 and7).

Experimental profiles depicted in Figure 10 correspond to a pulse test performed with 10 mL pH 1.5 solution, composed of 0.65 mol/L formic acid and 0.267 mol/L succinic acid, eluted with pH 1.5 H 2 SO 4 solution at 0.5 BV/h, using the same setup described above.

Figure 8 confirms clearly that the hybrid model predicts the chromatographic profile of both organic acids much better. More specifically, the hybrid model is able to predict the significant tailing of the succinic acid profile, unlike the Langmuir model. Moreover, the mean retention volume of formic acid is better predicted by the hybrid model, even though there is still a slight difference (0.1 BV) compared to experimental points. The slight discrepancies observed can be due to assumptions and simplifications in the present hybrid model. For example, all adsorption sites were considered equivalent. To be more accurate, adsorption sites on the upper side of and counter-anions should be considered more accessible than sites closer to the resin matrix (Fig. 4). In addition, some studies [START_REF] Neil | Adsorption of formic acid on clean and oxygen covered Pt(111)[END_REF][START_REF] Nejad | Concerted Pair Motion Due to Double Hydrogen Bonding: The Formic Acid Dimer Case[END_REF] showed that the original Langmuir competitive equation is not suitable for multi-layer adsorption. Thus, a more accurate model considering different types of adsorption sites might be an interesting approach to better predictions of multicomponent mixture separation by chromatography.

Conclusion

Strong anionic resins, often used in sulfate form in chemical or food and feed industries to separate organic acids by chromatography, require better retention mechanisms to improve the model's prediction. In this work, a hybrid model was developed and used to explain the elution profiles of succinic and citric acids, which display unusual tailing. Compared to classic Langmuir adsorption models, it also considers the effect of organic acid dissociation into anions that can be retained on resin ion-exchange sites. The model was solved with the CE/SE method.

Estimation of adsorption equilibrium coefficient ( and ) were consistent with the adsorption of organic acids by means of hydrogen bonding between the hydrogen atom of the carboxylic acid group and one lone electron pair of oxygen atoms of sulfate or hydrogen sulfate counter-anions. Present results confirm that the tailing of succinic and citric acid profiles is mainly caused by ion-exchange of their very low dissociated fraction, even at low pH (< 1.5). Organic acid dissociation and ion-exchange is, therefore, must be considered to develop a predictive model for succinic and citric acids. The prediction of the hybrid model was successfully validated in the case of a multi-component pulse test.

This hybrid model is currently implemented in a multi-column chromatography simulation program to predict and optimize the performance of simulated moving bed technologies. It is also useful to design new multi-fraction processes.

Appendices

SBA (strong base anion) resins selectivity coefficients [34, 35]

Anion

Figure 7

 7 Figure 7 presents simulation curves of pulse tests for different organic acids, to compare the Langmuir adsorption model with the hybrid model, which considers that organic acid can be retained by Langmuir adsorption but also by ion-exchange.

Figures and TablesFig. 1 .

 1 Figures and Tables

Fig. 2 .

 2 Fig. 2. Mass balance (convection-diffusion equation) and mass transfer model

Fig. 3 .

 3 Fig. 3. Schematical representation of the strong anion resin (type 1) in sulfate form.

Fig. 4 .

 4 Fig. 4. Schematic illustration of organic acid adsorption by hydrogen bonding on SO 4 2-or HSO 4 -counter-anions, which have eight lone electron pairs each. N = 0.5 for formic acid could be explained by dimer conformation, and means that up to 16 molecules could be adsorbed (only 4 of them are represented). N = 8 for citric acid corresponds to one molecule maximum adsorbed on each counter-anion. Indeed, citric acid can form up to 3 hydrogen bonds and the other five free sites could be cluttered. N = 4 for succinic acid corresponds to 2 molecules maximum adsorbed on each counter-anion, because it can form up to 2 hydrogen bonds and could clutter less free sites.

Fig. 5 .

 5 Fig. 5. The SE (solution element, gray part) and CE (conservation element, green and blue part) at position j and time n. Values at point A are calculated from points C and E at n-1/2; values at C or E are obtained from points at time n-1 via the CE and SE at the point at time n-1/2.

Fig. 6 .

 6 Fig. 6. Influence of (Fig.6-A on the top) and (Fig.6-B at the bottom) on organic acid outlet profile.

Fig. 7 .Fig. 8 .

 78 Fig.7.Influence of ion-exchange on outlet profile of formic (Fig.7-A), succinic (Fig.7-B) and citric acids (Fig.7-C)

Acknowledgments

This study was carried out in the Centre Européen de Biotechnologie et de Bioéconomie (CEBB), supported by Région Grand Est, Département de la Marne, Grand Reims and the European Union. In particular, the authors would like to thank Département de la Marne, Grand Reims, Région Grand Est and European Union with European Regional Development Fund (ERDF Champagne-Ardenne 2014-2020) for their financial support to the Chair of Biotechnology of CentraleSupélec.

Citric acid 300 ± 50 250 0.14-0.28 ( = 8) 0.12-0.14 ( = 8)


Omar Darwich

Ahmet-Sefa Ulucan

Quentin Bramas

Anissa Lamani

Anaïs Durand

Pascal Lafourcade

Perpetual Torus Exploration by Myopic Luminous Robots

We study perpetual torus exploration for swarm of autonomous, anonymous, uniform, luminous robots with a common chirality. We consider robots with only few capabilities. They have a finite limited vision (myopic), they can only see robots at distance one or two. We show that the problem is impossible with only two luminous robots and also with three oblivious robots (without light). We design an optimal algorithm for three luminous robots using two colors and with visibility one. We also propose an optimal algorithm with visibility two with four oblivious robots.

Introduction

Swarm robotics has drawn a lot of attention the past decade. Inspired by natural systems, a lot of investigations focused on how to reproduce autonomous behaviors observed in nature within artificial systems. Given a collection of autonomous mobile entities called robots, the main focus is to determine the minimum hypothesis in order for the robots to solve a given task. Robots can evolve either on a continuous 2D plane on which they can freely move or on a discrete universe, generally represented by a graph, where nodes indicate possible locations of the robots and the edges the possibility for the robots to move from one node to another.

In this paper, we assume that the mobile robots are autonomous (i.e. there is no central authority to coordinate their move), anonymous (i.e. they have no identity), uniform (i.e. they all execute the same algorithm) and luminous (i.e. they are endowed with lights of different colors). Moreover, they cannot communicate directly but are endowed by visibility sensors allowing them to sense their environment within a certain distance called visibility range. We assume myopic robots that can only sense at small distances. Robots operate in the well-known LCM model. That is, they operate in cycles which comprise three phases: Look, Compute, and Move. During the first phase (Look), robots take a snapshot of their environment using their visibility sensors. In the second phase (Compute), based on the taken snapshot, they first decide whether to move or remain idle and then whether they change their color. If they decide to move, they compute a neighboring destination. Similarly, they compute a new color if they decide to change it. Finally, in the last phase (Move), they move to the computed destination (if any) and they change their color (if they decided to). We consider the fully synchronous model (FSYNC) in which all robots execute the LCM cycle synchronously and atomically.

In the following, we investigate the case in which the robots have to solve the perpetual exploration problem. In this problem, robots evolve in a discrete universe and have to ensure that each location (node) is visited by at least one robot infinitely often. We are interested in torus shaped networks and focus on optimal exclusive solutions with respect to both the visibility range and the number of robots. Exclusiveness add an additional constraint on robots behavior as they can neither occupy the same node simultaneously or traverse the same edge at the same time.

Related work

The exploration problem is considered as one of the benchmarking tasks when it comes to robots evolving on graphs. Various topologies have been considered: lines [?], rings [?,?,?,?,?], tori [?], grids [?,?,?,?], cuboids [?], and trees [?]. Two variants of the problem has been investigated: (i) the perpetual exploration problem [?,?,?,?], considered in this paper, which requires the robots to visit each node of the graph infinitely often and (ii) the terminating exploration problem [?,?,?,?,?,?,?] which requires the robots to visit each node of the graph at least once and then stop moving.

Most of the investigations consider robots with unlimited visibility range allowing them to observe every node of the system [?,?,?,?,?,?,?,?]. Robots are in this case oblivious (i.e. they cannot remember past actions) and have to solve the terminating exploration problem. Myopic robots have also been considered in both variants of the problem [?,?,?,?,?]. When it comes to the perpetual exploration problem, an additional assumption has an impact on the feasibility of the task and the optimality of the proposed solutions. This assumption endow the robots with a common chirality. In fact, chirality is usually assumed when robots evolve in the continuous 2D Euclidean plan but some investigations have also considered it recently in the discrete universe. On finite grids, it has been shown that two (resp. three) synchronous robots with three colors (resp. one color) are sufficient to solve the problem when robots have visibility one and share a common chirality [?]. The case in which robots have no common chirality was investigated in [?]. It was proven that the problem is not solvable with only two robots having any finite number of colors and a finite visibility range. An optimal solution is also presented using only three robots having visibility range one, using only three colors. The case in which robots are oblivious and visibility range 2 was solved using five robots. In the case of infinite grids, assuming robots with visibility range one and few colors (O(1)), five (resp. six) synchronous robots are necessary and sufficient to solve the problem with (resp. without) the common chirality assumption [?,?]. Finally, in the case of cuboids, it has been shown in [?] that three synchronous robots with a common chirality endowed with five colors are necessary and sufficient to solve the perpetual exploration problem.

Contribution:

We first present two impossibility results: we start by showing that the perpetual torus exploration problem is not solvable with only two robots if the number of colors is finite and their visibility range is limited. We then show that three oblivious robots are not sufficient to solve the PTE problem. Next, we propose two optimal solutions A 2 3 and A 1 4 with respect to both the number of robots and the number of colors for the case of visibility one and two respectively.

Model

We consider a set R of n > 0 robots located on a torus.

A graph G = (V, E) is a (l, L)-torus (or torus for short) if |V | = l × L and for any v (i,j) ∈ V ; i ∈ [0, l -1], j ∈ [0, L -1]: -{v (i,j) , v ((i+1) mod l,j) } ∈ E, and -{v (i,j) , v (i,(j+1) mod L) } ∈ E.
The order on the nodes of G forms a coordinate system. For example node v (i,j) is at coordinate (i, j), or, the node is at column i and row j. For simplicity we note node (i, j) instead of v (i,j) . This order/coordinate is used for the analysis only, i.e., robots cannot access it.

We assume a discrete time where at each round, the robots synchronously perform a Look-Compute-Move cycle. In the Look phase, a robot gets a snapshot of the subgraph induced by the nodes within distance Φ ∈ N * from its position. Φ is called the visibility range of the robots. The snapshot is not oriented in any way as the robots do not agree on a common North. However, it is implicitly ego-centered since the robot that performs a Look phase is located at the center of the subgraph in the obtained snapshot. Robots agree on a common chirality. Then, each robot computes a destination (either Up, Left, Down, Right or Idle) based only on the snapshot it received. Finally, it moves towards its computed destination. We also assume that robots are opaque, i.e., they obstruct the visibility in such way that if three robots are aligned, the two extremities cannot see each other. We forbid any two robots to occupy the same node simultaneously. A node is occupied when a robot is located at this node, otherwise it is empty.

Robots may have lights with different colors that can be seen by robots within distance Φ from them. We denote by Cl the set of all possible colors. For simplicity, we assume that all tore has dimensions l × L where l, L ≥ nΦ + 1.

The state of a node is either the color of the light of the robot located at this node, if it is occupied, or ⊥ otherwise. In the Look phase, the snapshot includes the state of the nodes (within distance Φ, including its current node). During the compute phase, a robot may decide to change the color of its light.

In all our algorithms, we also prevent any two robots from traversing the same edge simultaneously. Since we already forbid them to occupy the same position simultaneously, this means that we additionally prevent robots from swapping their position. Algorithms verifying this property are said to be exclusive. However, to be as general as possible, we do not make this additional assumption in our impossibility results.

In the following, we borrow some of the definitions already presented in [?].

Configurations A configuration C in a torus G(V, E) is a set of pairs (p, c), where p ∈ V is an occupied node and c ∈ Cl is the color of the robot located at p. A node p is empty if and only if ∀c, (p, c) / ∈ C. We sometimes just write the set of occupied nodes when the colors are clear from the context.

Views We denote by G r the globally oriented view centered at the robot r, i.e., the subset of the configuration containing the states of the nodes at distance at most Φ from r, translated so that the coordinates of r is (0, 0). We use this globally oriented view in our analysis to describe the movements of the robots: when we say "the robot moves Up", it is according to the globally oriented view. However, since robots do not agree on a common North, they have no access to the globally oriented view. When a robot looks at its surroundings, it instead obtains a snapshot. To model this, we assume that the local view acquired by a robot r in the Look phase is the result of an arbitrary indistinguishable transformation on G r . The set IT of indistinguishable transformations contains the rotations of angle 0 (to have the identity), π/2, π and 3π/2, centered at r. Moreover, since robots may obstruct visibility, the function that removes the state of a node u if there is another robot between u and r is systematically applied to obtain the local view. Finally, we assume that robots are self-inconsistent, meaning that different transformations may be applied at different rounds.

It is important to note that when a robot r computes a destination d, it is relative to its local view f (G r ), which is the globally oriented view transformed by some f ∈ IT . So, the actual movement of the robot in the globally oriented view is f -1 (d). For example, if d = Up but the robot sees the torus upside-down (f is the π-rotation), then the robot moves Down = f -1 (Up). In a configuration C, V C (i, j) denotes the globally oriented view of a robot located at (i, j).

Algorithm An algorithm A is a tuple (Cl , Init, T ) where Cl is the set of possible colors, Init is a mapping from any considered torus to a non-empty set of initial configurations in that torus, and T is the transition function V iews → {Idle, Up, Left, Down, Right} × Cl , where V iews is the set of local views. When the robots are in Configuration C, a configuration C ′ obtained after one round satisfies: for all ((i, j), c) ∈ C ′ , there exists a robot in C with color c ′ ∈ Cl and a transformation f ∈ IT such that one of the following conditions holds:

- ,c). We denote by C → C ′ the fact that C ′ can be reached in one round from C (n.b., → is then a binary relation over configurations). An execution of Algorithm A in a torus G is then a sequence (C i ) i∈N of configurations such that C 0 ∈ Init(G) and ∀i ≥ 0,

((i, j), c ′ ) ∈ C and f -1 (T (f (V C (i, j)))) = (Idle, c), -(((i -1) mod l, j), c ′ ) ∈ C and f -1 (T (f (V C ((i -1) mod l, j)))) = (Right, c), -(((i + 1) mod l, j), c ′ ) ∈ C and f -1 (T (f (V C ((i + 1) mod l, j)))) = (Lef t, c), -((i, (j -1) mod L), c ′ ) ∈ C and f -1 (T (f (V C (i, (j -1) mod L)))) = (U p, c), or -((i, (j +1) mod L), c ′ ) ∈ C and f -1 (T (f (V C (i, (j +1) mod L)))) = (Down
C i → C i+1 .
Definition 1 (Perpetual Torus Exploration). An algorithm A solves the Perpetual Torus Exploration (PTE) problem if in any execution (C i ) i∈N of A and for any node (i, j) ∈ V of the torus and any time t, there exists t ′ > t such that (i, j) is occupied in C t ′ . Notations. ⃗ t (i,j) (C) denotes the translation of the configuration C of vector (i, j).

Impossibility Results

Lemma 1. Let A be an algorithm using a set R of n > 0 robots. If A solves the exploration problem (perpetual or with termination) for any torus then, there exists a tori such that for any execution (C i ) i∈N of A on this torus, there is a configuration C i such that the distance between the two farthest robots is at least 2Φ + 3.

Proof. We proceed by contradiction. Assume, there is an algorithm A that solves the PTE problem and let 0 < B be the farthest any of the robots will be from each other, in any torus. Let (C i ) i∈N be the execution of A on a very large torus l, L ≫ B. When all robots are at distance at most B, then the occupied positions are included in a square sub-grid of size B × B. Since the number of possible configurations included in a subgrid of size B × B is finite, there must be two indices t 1 and t 2 , when the positions and colors of the robots in the corresponding sub-grids are the same, formally, such that C t2 = ⃗ t (i,j) (C t1 ) and t 1 < t 2 for a given translation ⃗ t (i,j) . By making the adversary choose the same rotation, the movements done by the robots in configurations C t1 and C t2 are the same as each robot has the same globally oriented view in both configurations, only their positions on the torus change. Thus C t2+1 = ⃗ t (i,j) (C t1+1 ) and so on so forth, so that ∀x, C t2+x = ⃗ t (i,j) (C t1+x ). We obtain that the configurations are periodic with period p = t 2 -t 1 , up to translation.

Suppose, that the torus being explored is of dimensions l×L with l = 3np 3 max(|i|, 1) and L = 3np 3 max(|j|, 1). The dimensions of the torus are proportional to the non-null scalar components of translation ⃗ t (i,j) i.e., i3np 3 ≡ 0 mod l and j3np 3 ≡ 0 mod L. This means that,

( ⃗ t (i,j) ) 3np 3 (C t1 ) = ⃗ t (i3np 3 ,j3np 3 ) (C t1 ) = ⃗ t (0,0) (C t1 ) = C t1 .
Since translation ⃗ t (i,j) is performed in p rounds, after p × 3np 3 = 3np 4 rounds, all robots will retake their initial positions, so the whole configuration is periodic with period 3np 4 . In this setting, a node is visited infinitely often if and only if it is visited between round t 1 and t 1 + 3np 4 . Now we have to prove that some nodes are left unvisited between round t 1 and t 1 + 3np 4 . Between time t 1 and t 1 + 3np 4 , each robot visits at most 3np 4 nodes, hence all the robots visit at most n × 3np 4 nodes after t 1 . However, there are at least 9n 2 p 6 ≤ l × L nodes in the torus. Hence, there exist some nodes which are not visited infinitely often, which is a contradiction.

Note that we only proved there are some nodes that are not perpetually visited. Nevertheless, observe that at most nt 1 nodes are visited before t 1 and we can increase arbitrarily the chosen period p by a factor f ∈ N * without changing the result (in particular t 1 does not depend on f ). By taking f ≥ 1 such that 9n 2 (f p) 6 -3n 2 (f p) 4 > nt 1 , we have that the number of visited nodes (before or after t 1 ) is nt 1 + 3n 2 (f p) 4 and is smaller than the number of nodes in the torus (9n 2 (f p) 6 ), hence there is at least one node that is never visited. This implies that the impossibility also holds for nonperpetual algorithms as well (where each node must be visited at most once).

We restate the following lemma proven in [?].

Lemma 2. A robot with self-inconsistent compass and that sees no other robot, either stays idle or the adversary can make it alternatively move between two chosen adjacent nodes.

Theorem 1. It is impossible to solve the exploration problem (perpetual or with termination) with two myopic robots equipped with self-inconsistent compasses that agree on a common chirality.

Proof. By Lemma ??, there is a torus and a configuration where the two robots are at distance 2Φ + 3 from each other. In this case the two robots are isolated. By Lemma ??, the two robots will remain idle or the adversary can make them alternatively move between two nodes, never being in vision from each other and never visiting another node.

Theorem 2. It is impossible to solve the exploration problem (perpetual or with termination) with three anonymous, oblivious and myopic robots equipped with self-inconsistent compasses that agree on a common chirality.

Proof. By Lemma ??, there is a torus and a configuration where the distance between the two farthest robots is 2Φ + 3 from each other. We have one of the two following possibilities, (i) there are three isolated robots, or (ii) there is an isolated robot and two robots in vision from each other.

In the first case, it is easy to see that the three isolated robots cannot explore the torus because, by Lemma ??, they have to stay idle or the adversary can make them alternatively move between two nodes, never being in vision from each other and never visiting another node.

In the second case, the two robots that see each other cannot travel together in a direction (because they have the same view). All they can do is get either closer to each other or further from each other. Formally, there is a point P at the middle of the two robots and, if they stay in vision, they will always be at the same distance from that point. The two robots can explore a subgrid Φ × Φ centered at a given middle point. This point is at distance at least 3Φ 2 + 2 from the isolated robots. If the two robots in vision gets isolated from one another, they will be at distance Φ 2 + 1 from the middle point. In this case, the closest robot to the originally isolated robot will be at distance Φ + 1. Now the three robots are isolated, and, as in the first case, they cannot explore the torus.

5 Visibility range one: A 2

3

We present an algorithm, denoted by A 2 3 , which assumes a visibility range one and uses three robots and two colors. By Theorem ??, A 2 3 is optimal w.r.t. the number of robots, and by Theorem ??, A 2 3 is also optimal w.r.t. the number of colors. Animations are made available online [?] to help the reader visualize the algorithm.

The idea of the algorithm is to make the robots alternate between exploring a row and exploring a column. To explore the whole torus, robots move so that all the nodes of the torus are explored. More precisely, after exploring row r i and column c j , the robots will proceed at exploring row r i-1 mod L and then column c j-1 mod l and so on.

Initially the robots are co-linear with respectively color L, F , F3 . The line of the torus on which they are located is considered as a row. The robot with color F which does not sense the leader moves up changing its color to L while the two other robots move along their current row in the following manner: the robot initially with color L moves away from the one with color F and the remaining robot just follows it.

To explore a row (resp. column), one robot stays idle while the two others travel in a straight line along the nodes of the row (resp. column) being explored until they reach the idle robot. The idle robot is located on an neighboring row (resp. column). The idle robot has color L and is called the landmark. The two robots traveling together on a straight line have different colors. One robot, called the follower, has color F and the other robot, called the leader, has color L. To explore a row (resp. column), the two robots have to be next to each other on that row (resp. column). The follower always follows the leader and leader always moves away from the follower. This is done by executing the rules presented in Fig ? ?. 4 The tricky part of this algorithm is how robots switch from exploring a column to exploring a row and vice versa. When exploring a column, the robot left behind (aka the landmark) is on the right side of the traveling group. When the leader of the traveling group reaches the landmark, it moves away from the landmark on its current row and updates its color to F . Meanwhile, the follower continues to follow the leader. In the next round, the three robots are aligned on the same row. The landmark then moves away from the follower and remains on its row followed by the follower. These two robots become the new traveling group. Whereas the leader, moves to the next row so that it becomes on the left side of the traveling group. That is, the landmark and the leader switch their roles and the new traveling group proceed at the exploration of the row on which there are located The traveling group are now exploring a row, when they reach the landmark again, the landmark is this time, on the right side. The robots proceed to move to the next column to be explored. More precisely, when the leader reaches the landmark, it continues forward on its current row and changes its color to F . The follower also continues to move towards the leader. After one round, the robots will be in a L-shaped form with the two robots colored F and the one robot colored L. In the next round, the two robots on the left form the new traveling group and they both move to explore the new column. The robot on the right, moves down and changes its color to L, it becomes the new landmark. The set of rules relevant to this sequence is in It is important to note that every node on a column/row is visited during the exploration of that column/row. Also, the landmark moves two nodes to the left and one node up when going from exploring a column to exploring a row. And, it moves one node to the right and two nodes downward when going from exploring a row to exploring a column. This means that between two consecutive columns (rows) exploration, the landmark moves one node to the left and one node downward. Theorem 3. A 3 2 solves the PTE problem with three robots and two colors.

Proof. By induction on l × L, where l is the number of columns and L is the number of rows of the torus.

We have validated the base case, for torus of size 4 × 4, using our simulation tool. Such a checking is easy since, from a given initial configuration, there is only one possible execution (the algorithm is well-defined and the execution is synchronous). So, we just have to execute the algorithm until reaching an already encountered configuration from which all the nodes have been visited.

We assume now that A 3 2 solves the PTE problem in all tori x × y with 4 ≤ x ≤ l and 4 ≤ y ≤ L for some values l, L ≥ 4 and show that A 3 2 solves also the PTE problem in a torus of size l × (L + 1) and (l + 1) × L.

Consider first the torus of size l × (L + 1). Then, it is easy to see that after adding one row, our algorithm still solves the PTE problem. Indeed, when robots are traveling upward (i.e. they are exploring a row), they move in a straight line periodically until they reach the landmark, so adding one row just increases by one the number of

Fig. 1 :

 1 Fig. 1: Rules for moving straight.

Fig. 3 :

 3 Fig. 2: Rules for switching from moving upward to sideward.

Fig. 5 :

 5 Fig. 4: Rules for switching from moving upward to sideward.

Table ? ?

 ? summarizes our contribution:

	Visibility # Robots # Colors	Algorithm
	finite	2	finite Impossible (Thm. ??)
	finite	3	1	Impossible (Thm. ??)
	1	3	2	A 2 3
	2	4	1	A 1 4

Table 1 :

 1 Summary of our results.

Note that any reachable configuration can be an initial configuration

In all figures, colored letters inside nodes indicate the color of the robots occupying the nodes. Moreover when a colored letter is given next to a node, it indicates which color the robot will take in the next round.

times they perform their periodic movement. And, when robots are traveling sideward (i.e. they are exploring a column), they visit all the nodes of the corresponding column. Now, for the torus of size (l + 1) × L. The same argument from the torus of size l × (L + 1) could be used. When robots are traveling sideward, they will perform an extra step for the added column. And, when they travel upward, they will revisit the same nodes visited during the exploration of rows.

6 Visibility range two: A 1 4

We present an algorithm, denoted by A 1 4 , which assumes a visibility range two and uses four oblivious robots. A 1 4 is optimal w.r.t. the number of colors. By Theorem ??, A 1 4 is optimal w.r.t. the number of robots, for oblivious robots. Animations are made available online [?] to help the reader visualize the algorithm.

The idea of the algorithm is to make the robots explore the torus rows by rows in a given direction. This is achieved as follows: Three robots, referred to as the traveling group, move to explore three adjacent rows at the same time, and one robot is left behind to be used as their landmark. When the traveling group reaches the landmark, all four robots perform a three rounds sequence to move to the next rows to be explored.

When exploring the rows, the traveling group form a > shape. That is, two robots are located on the same column separated by one empty node, denoted u. And, on the right of u, the third robot is placed. The three robots move to the right until they sense the landmark. Note that the direction is pointed by the third mentioned robot. The landmark is left behind so that the traveling group knows when they are done exploring the current rows and have to move to the next ones. Note that the landmark is on the same row as the top most robot. When that robot is one node away from the landmark it goes down, same for the landmark since they have the same view. The bottom robot keeps going right because it does not see the landmark. And, the center robot stays idle. After one round, the robots form a T-shape. The rules executed by the robots are presented in Fig ?

Fig. 7: Rules executed when robots initiate rows change.

From the T-shape, the robot move to create a reverse L shape i.e. the two robots in the center of the T-shape move down while the robot on the right goes left. Fig ? ?  presents the rules executed during this process. 8: Rules for the creation of the reverse L shape.

Within the reverse L shape, three robots are co-linear (the ones located on the long side). Among these robots, the one in the middle moves to the right to recreate the > shape while all the other robots remain idle. Refer to the rule presented Fig ??. That is, after three rounds the robots changes rows and the > shape is built again. Robots form initially the reverse L shape.

Theorem 4. A 4 1 solves the PTE problem with four oblivious robots.

Proof. By induction on l × L, where l is the number of columns and L is the number of rows of the torus. Similar to the proof of Theorem ??, we have validated one base case for l, L = 9, using our simulation tool.

We assume now that A 4 1 solves the PTE problem in all tori x × y with 9 ≤ x ≤ l and 9 ≤ y ≤ L for some values l, L ≥ 9. We should show that A 4 1 solves the PTE problem in the tore of size l × (L + 1) and (l + 1) × L.

Consider first the torus of size (l + 1) × L. When we add one column the traveling group will have to perform an extra round to reach the landmark again as the robots perform a periodic movement when traveling until they observe the landmark. Now, consider the torus of size l × (L + 1). When we add a row. The robots will have to perform an extra row exploration: an additional three round sequence to change rows followed by the row exploration.

Conclusion

We presented two optimal solutions for the PTE problem with respect to both the number of robots and the number of colors when robots share a common chirality and have visibility one and two respectively. Indeed, we have shown that three robots endowed with two colors are necessary and sufficient to solve the problem when robots have visibility one and four oblivious robots are necessary and sufficient to solve the problem when robots have visibility two. One direct open question is to extend the study to consider (L, l)-tori such that l, L < nΦ + 1. Ad-hoc solutions might be needed in this case as robots observe the same robots on different sides. Another interesting extension would be to investigate the case in which robots are completely disoriented, i.e., they do not have a common chirality. We conjuncture that three robots remain sufficient to solve the problem with an additional color in the case where robots have visibility one and an additional robot might be needed in the case of oblivious robots with visibility two.