Stergios G Yiantsios
email: yiantsio@auth.gr

Stefania K Serpetsi

Frédéric Doumenc

Béatrice Guerrier

Surface deformation and film corrugation during drying of polymer solutions induced by Marangoni phenomena

Keywords: Marangoni effect, Free convection, Drying, Instability

Solutal Marangoni phenomena occurring during drying of polymer solutions are investigated through numerical simulations. Bénard-Marangoni convection and free surface deformation are considered, taking into account the strong dependence of viscosity with solute concentration. Simulations show that free surface deformation which is resisted by capillary pressure reaches an overall magnitude that is proportional to the capillary number. Decreasing Ca also inhibits the interface deformation by favoring the formation of plumes that enhance mixing of the solute. For the range of capillary numbers corresponding to experiments, free surface deformation is not significant. However, Bénard-Marangoni convection induces non-uniformity of solute concentration and the varying viscosity leads to accumulation of solute at the periphery of convective cells that is observed up to the end of drying. Though simulations have been performed with dimensionless numbers significantly milder than the experimental ones, they are in agreement with experimental observations presented here and in the literature, suggesting that wrinkles may be viewed as the signature of Bénard-Marangoni convection.

Introduction

Marangoni phenomena have a wide range of manifestations, i.e. in crystal growth [START_REF] Schatz | Experiments on thermocapillary instabilities[END_REF], heat or mass transfer enhancement [START_REF] Scriven | The Marangoni effects[END_REF], biological systems [START_REF] Jensen | The spreading of heat or soluble surfactant along a thin liquid film[END_REF], thin film coatings and paints, etc. [START_REF] Abbasian | Study of different planforms of paint's solvents and the effect of surfactants (on them)[END_REF][START_REF] Bassou | Role of Bénard-Marangoni Instabilities during Solvent Evaporation in Polymer Surface Corrugations[END_REF][START_REF] Yiantsios | Marangoni flows during drying of colloidal films[END_REF][START_REF] Yiantsios | A mechanism of Marangoni instability in evaporating thin films due to soluble surfactant[END_REF][START_REF] Serpetsi | Stability characteristics of solutocapillary Marangoni motion in evaporating thin films[END_REF]. They have also been proposed as a possible cause for the development of organized structures in thin films with potentially promising properties for organic photovoltaics and photodiodes [START_REF] Heriot | An interfacial instability in a transient wetting layer leads to lateral phase separation in thin spin-cast polymer-blend films[END_REF]. The most common factors inducing changes in surface tension, and thus Marangoni phenomena, are gradients of temperature or concentration of dissolved species and surfactants [START_REF] Colinet | Nonlinear Dynamics of Surface-TensionDriven Instabilities[END_REF]. Theoretically, the mechanisms leading to hydrodynamic instabilities have been elucidated first by Pearson [11] for the surface-tension-driven flow (Bénard-Marangoni convection). Another instability mechanism due to surface deformation can also be present: it is known as the long-wave instability [START_REF] Scriven | On cellular convection driven by surface-tension gradients: effect of mean surface tension and surface viscosity[END_REF][START_REF] Reichenbach | Linear perturbation analysis of surface tension driven convection at a plane interface[END_REF][START_REF] Goussis | On the thermocapillary instabilities in a liquid layer heated from below[END_REF]. Earlier theoretical and experimental studies of thermocapillary phenomena are reviewed in [START_REF] Schatz | Experiments on thermocapillary instabilities[END_REF][START_REF] Davis | Thermocapillary instabilities[END_REF]. Aspects particular to thin-film flows can be found in the reviews of Oron et al. [START_REF] Oron | Long-scale evolution of thin films[END_REF] and Craster and Matar [START_REF] Craster | Dynamics and stability of thin liquid films[END_REF].

In the specific case of evaporating thin films, thermocapillary phenomena arise spontaneously because temperature gradients develop due to the energy transported through the film that is required to provide the latent heat for evaporation [START_REF] Mancini | Pattern formation without heating in an evaporative convection experiment[END_REF][START_REF] Moussy | Influence of evaporation on Bénard-Marangoni instability in a liquid-gas bilayer with a deformable interface[END_REF]. In multicomponent evaporating films surface tension gradients develop not only due to the effects of temperature but of species concentration as well. Several experimental studies confirm the influence of Marangoni phenomena, with stagnant or spin-coated films of various solutions where the volatile solvent invariably had lower surface tension than the solute [START_REF] Bassou | Role of Bénard-Marangoni Instabilities during Solvent Evaporation in Polymer Surface Corrugations[END_REF][START_REF] Birnie | Surface tension evolution during early stages of drying of sol-gel coatings[END_REF][START_REF] Bormashenko | On the mechanism of patterning in rapidly evaporated polymer solutions: Is temperature-gradient-driven Marangoni instability responsible for the large-scale patterning?[END_REF][START_REF] Strawhecker | The critical role of solvent evaporation on the roughness of spin-cast polymer films[END_REF]. Theoretical and numerical studies have also been devoted to the study of such configurations [START_REF] Yiantsios | Marangoni flows during drying of colloidal films[END_REF][START_REF] Yiantsios | A mechanism of Marangoni instability in evaporating thin films due to soluble surfactant[END_REF][START_REF] Serpetsi | Stability characteristics of solutocapillary Marangoni motion in evaporating thin films[END_REF][START_REF] De Gennes | Instabilities during the evaporation of a film: Non-glassy polymer + volatile solvent[END_REF][START_REF] Ha | Onset of Marangoni instability of a two-component evaporating droplet[END_REF][START_REF] Machrafi | Time-dependent Marangoni-Bénard instability of an evaporating binary-liquid layer including gas transients[END_REF][START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF][START_REF] Doumenc | Free convection in drying binary mixtures: solutal versus thermal instabilities[END_REF].

The present study is based on two previous and complementary approaches concerning solutal isothermal problems. On the one hand, theoretical analyses including linear stability and nonlinear simulations are presented in Serpetsi and Yiantsios' work [START_REF] Serpetsi | Stability characteristics of solutocapillary Marangoni motion in evaporating thin films[END_REF]. They account for a deformable free surface but assume constant viscosity. Distinct modes of instability similar to those of Pearson [START_REF] Pearson | On convection cells induced by surface tension[END_REF] and Scriven and Sternling [START_REF] Scriven | On cellular convection driven by surface-tension gradients: effect of mean surface tension and surface viscosity[END_REF] were found. The former mode leads to cell development with minimal interfacial deformation, whereas the latter is oscillatory and induces significant interfacial deformation. Moreover, it was found that systems where the solute decreases the surface tension are also susceptible to instability in the form of a long-wavelength, non-oscillatory deformational mode. A related phenomenon of solutocapillary motion arises in the leveling of evaporating paints where a characteristic oscillatory behavior is observed. This has been explained by Overdiep [START_REF] Overdiep | The leveling of paints[END_REF] as a result of the synergistic action of capillary pressure and the concentration induced surface tension gradients.

On the other hand, linear stability analysis and direct numerical simulations have been performed by Doumenc et al. [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF][START_REF] Doumenc | Free convection in drying binary mixtures: solutal versus thermal instabilities[END_REF], taking into account the variation of viscosity with solute concentration, but disregarding the deformation of the interface. They focus on the onset of Bénard-Marangoni convection and point out the effect of varying viscosity. The aim of the present study is to bring together the two approaches and to get a view of the evolution of the interface deformation and concentration field during the drying, considering variable viscosity.

Indeed, one aspect particularly relevant to evaporating polymer solutions is the patterning of the final dried films, in the form of corrugations, wrinkles, pinholes, etc. [START_REF] Abbasian | Study of different planforms of paint's solvents and the effect of surfactants (on them)[END_REF]. Theoretical studies are virtually non-existent, except from a few that are limited to the long-wavelength approximation and the lubrication regime for the fluid flow and scalar transport [START_REF] Howison | A mathematical model for drying paint layers[END_REF][START_REF] Eres | Three-dimensional direct numerical simulation of surface-tension-gradient effects on the leveling of an evaporating multicomponent fluid[END_REF]. The development of patterns in thin films encompasses a complex and multifaceted array of phenomena attributed to various factors. For some authors, a notable factor is the formation of an elastic skin and the incompatibility of its mechanical properties with those of the substrate or the lateral boundaries. Considerations and analyses along this line have been presented in the literature to explain various manifestations of wrinkling pattern development [START_REF] Cerda | Geometry and physics of wrinkling[END_REF][START_REF] Soumendra | Mechanism of wrinkle formation in curing coatings[END_REF][START_REF] Chen | Wrinkling instabilities in polymer films and their applications[END_REF]. On a different length scale well known are phenomena of phase separation in polymer blends and in evaporating solutions [START_REF] Xue | Phase separation induced ordered patterns in thin polymer blend films[END_REF][START_REF] Xu | Surface textures induced by convection in thin films of polymeric and polymerizable fluids[END_REF][START_REF] Jaczewska | Polymer vs solvent diagram of film structures formed in spin-cast Poly(3-alkylthiophene) blends[END_REF]. However, Marangoni instabilities are also considered a potent factor that may lead to patterned dry films [START_REF] Bassou | Role of Bénard-Marangoni Instabilities during Solvent Evaporation in Polymer Surface Corrugations[END_REF][START_REF] Bormashenko | On the mechanism of patterning in rapidly evaporated polymer solutions: Is temperature-gradient-driven Marangoni instability responsible for the large-scale patterning?[END_REF][START_REF] Xu | Surface textures induced by convection in thin films of polymeric and polymerizable fluids[END_REF][START_REF] Sakurai | Highly organized macroscopic pattern observed on the free surface of a polymer film cast from a solution[END_REF]. Notably, experiments performed by Bassou and Rharbi [START_REF] Bassou | Role of Bénard-Marangoni Instabilities during Solvent Evaporation in Polymer Surface Corrugations[END_REF] with polystyrene/toluene solutions show a strong correlation between the convection cell morphology and the final corrugation of the dried film. From these observations, these authors infer that the final state of the film is the signature of solutal convection which redistributes the solute along the thin film, inducing lateral concentration non-uniformities, which might be finally frozen when most of the volatile solvent is removed by evaporation. Following this track, focus is placed in this paper on solutal Marangoni phenomena that could potentially be responsible for the wrinkling behavior observed in the experiments. Through numerical simulations, we explore the onset and development of solutal Bénard-Marangoni convection and consequently the development of concentration non-uniformities. We also investigate the deformation of the interface that could also results in corrugations if it is significant and continues until the end of the drying. The question addressed in this paper is therefore the relative importance of these two phenomena, according to values of capillary and Marangoni numbers, and the role played by viscosity variation with concentration.

The paper is organized as follows. First an experimental section presents some drying experiments performed with polyisobutylene/toluene solution, showing the drying and corrugation kinetics. Model assumptions, governing equations and numerical approach are developed in Section 3. The last section is devoted to the analysis of simulation results, focusing first on the free surface deformation and then on the solute distribution along the film.

Nomenclature

Experimental

Drying experiments were performed as follows: polymer solution was poured in a glass dish (diameter about 0.11 m) and put on a balance located in an extracting hood at room temperature. The polymer/solvent solution is polyisobutylene/toluene. Polyisobutylene (PIB) was supplied by Sigma-Aldrich with Mw = 5 × 10 2 kg/mol and polydispersity = 2.5. As is well known in polymer solutions, the viscosity strongly depends on solvent content. Fig. 1 displays the viscosity measurements for the PIB/toluene solution used in the experiments. An empirical law interpolated from measurements is used in the simulations of the present study, to describe the variation of viscosity with polymer concentration [START_REF] Doumenc | Free convection in drying binary mixtures: solutal versus thermal instabilities[END_REF]: 𝑙𝑜𝑔 10 (𝜇) = ∑ 𝛼 𝑖 4 𝑖=0 𝑌 𝑖 , with 𝑌 = 𝑙𝑜𝑔 10 (𝜑 𝑃𝐼𝐵 ), 𝜑 𝑃𝐼𝐵 being the polymer volume fraction, and 𝛼 = (8.235, 14.02, 6.575, 1.392, 0.1114). Mutual diffusion coefficient and solvent activity are also concentration dependent, but mainly for the concentrated regime and they can be assumed constant for polymer mass fraction less than about 0.6. There mutual diffusion coefficient is about 10 -10 m 2 /s and solvent activity is close to one [START_REF] Doumenc | Mutual diffusion coefficient and vapor liquid equilibrium data for the system PIB/Toluene[END_REF]. Accordingly, after a short thermal transient regime, the evaporation flux is constant as long as the polymer mass fraction stays less than 0.6 [START_REF] Doumenc | Free convection in drying binary mixtures: solutal versus thermal instabilities[END_REF]. Indeed, due to the thermodynamic behavior of polymer solutions, the evaporation rate of the solution is similar to that of pure solvent even if the concentration of the solute increases at the free surface, as long as the polymer mass fraction does not reach 60% [START_REF] Doumenc | Mutual diffusion coefficient and vapor liquid equilibrium data for the system PIB/Toluene[END_REF]. While both thermal and solutal instabilities may occur during the thermal transient regime [START_REF] Doumenc | Free convection in drying binary mixtures: solutal versus thermal instabilities[END_REF], the temperature is uniform during the constant flux regime and thus only solutal convection is expected. At the end of this constant evaporation flux regime, the flux decreases due to solvent activity reduction. It is during this transition that, for some configurations (initial thickness and initial polymer concentration), corrugations of the surface appear and stay frozen in time. Typical time evolution of evaporation flux and of wrinkle invasion during the transition is shown in Fig. 2. Evaporation flux is derived from mass measurements performed during the drying experiment. Thanks to a camera located above the solution, the percentage of wavy surface is obtained through image processing of surface top views.

A quantitative estimation of the amplitude of the wrinkles was made using a Bruker Nano surface profilometer Dektak XT for the configuration corresponding to an initial thickness of 1.06 mm and an initial polymer volume fraction of 5.07%. Surface visualization is shown in Fig. 3, while Fig. 4 displays typical profiles. Standard deviation estimated on the whole surface is about 15 µm, the dry mean film thickness being about 54 µm. Another set of experiments have been performed by Bassou and Rharbi with PS/toluene solutions. Experimental conditions are close to the ones of the experiments on PIB/toluene. Viscosity also strongly increases with polymer concentration, but unlike PIB which is rubbery at room temperature, PS solution undergoes a glass transition when drying. However, the amplitude of wrinkles is quite similar for those two systems (cf. [START_REF] Bassou | Role of Bénard-Marangoni Instabilities during Solvent Evaporation in Polymer Surface Corrugations[END_REF] for more details).

Mathematical and numerical formulation

Model assumptions and governing equations

The model adopted to describe the phenomena observed experimentally considers a Newtonian and incompressible liquid film containing dissolved polymer. The liquid viscosity is a strong function of polymer concentration. The film is resting on a flat horizontal substrate and the volatile solvent is evaporating in the overlying passive gaseous phase at a constant rate, E', corresponding to the constant evaporation flux regime described in the experimental section. Here E' is defined as a volumetric evaporation rate (mass rate divided by the liquid solvent density) and has dimensions of velocity. It is assumed that surface tension is a linear function of solute concentration and that the solute has a higher surface tension than the solvent. The initial transient thermal regime is disregarded and the problem is assumed isothermal.

The initial mean film thickness, H0, is taken as the length scale and the initial solute volume fraction, ω0, as a measure of concentration variations. Specific volumes of solute and solvent are assumed to be the same, so that no distinction is made between mass fraction and volume fraction and the solution density ρ is constant.

Concentration is scaled by ρω0. If the difference between the surface tensions of solute and solvent is Δσ, then ω0Δσ may be taken as characteristic measure of surface tension variations. Then, the balance of surface tension gradient forces with viscous stresses at the interface may be used to define a velocity scale, U0 = ω0Δσ/μ0, where μ0 is the initial viscosity of the solution. A viscous scaling for the pressure is adopted (P0 = μ0 U0/H0). Then, the equations governing the fluid motion may be written in dimensionless form as follows [START_REF] Serpetsi | Stability characteristics of solutocapillary Marangoni motion in evaporating thin films[END_REF]:

(1) which together with the continuity equation, ∇. 𝒖 = 0, determine the velocity field. Here, Sc is the Schmidt number, defined as Sc = υ0 /D, where the diffusivity, D, of the polymer in the solvent is assumed constant. M is the Marangoni number defined as

𝑀 = 𝜔 0 𝛥𝜎 𝐻 𝑂 𝜇 0 𝐷 . ( 2
)
The ratio M/Sc is essentially a Reynolds number for the flow and μ in Eq. ( 1) is the relative viscosity of the liquid, scaled by the initial one, μ0. Since the Schmidt number is very large the inertia terms may be ignored and the flow can be described by the Stokes equations, as demonstrated in Appendix A with scaling arguments. These equations are supplemented by the no-slip condition at the substrate and the stress condition at the interface

(3)
Here τ is the viscous stress tensor of the liquid, n the normal vector to the interface pointing towards the gas, ω the dimensionless polymer mass fraction, and h the local dimensionless film thickness. Ca is the capillary number defined as

𝐶𝑎 = 𝜇 0 𝑈 𝑂 𝜎 0 = 𝜔 0 𝛥𝜎 𝜎 0 , ( 4
)
and Bo is a modified Bond number

𝐵𝑜 = 𝜌𝑔 ℎ 0 2 𝛥𝜎 ( 5
)
quantifying the effects of gravity. Such effects are ignored in the present simulations assuming relatively thin films and weak gravity forces compared to the Marangoni forces. The first term on the right-hand-side of Eq. ( 3) corresponds to the capillary pressure, expressed through the curvature, ∇. 𝒏, and the last one corresponds to the Marangoni stresses, expressed through the surface gradient operator, ∇ 𝑠 .

The liquid domain is shrinking due to evaporation and may also be deforming due to the flow. Its evolution is described by the kinematic condition, which written in terms of the instantaneous local thickness, h, has the form [START_REF] Yiantsios | Marangoni flows during drying of colloidal films[END_REF] u and v are the horizontal and vertical component of the velocity. Here, the dimensionless parameter E, defined as E = E'/U0, is the ratio of the time scales for motion due to Marangoni effects and for evaporation. If such effects are to be manifested, E has to be a small parameter.

Obviously, there is coupling between the flow and the distribution of the solute, the transport of which is described by the species conservation equation:

(7)
Hence, the Marangoni number takes the role of a Peclet number for species transport. The above equation is supplemented by the no-flux condition of the solute at the substrate and the condition that the polymer does not evaporate at through the free interface, which is written as (8)

Numerical approach

Galerkin finite-element formulation

The direct numerical solution of the coupled set of the full equations was based on the Galerkin finite element method. Since the flow domain is changing in time and needs to be found as part of the solution, the ALE (Arbitrary Lagrangian-Eulerian) [START_REF] Ho | A Legendre spectral element method for simulation of unsteady incompressible viscous free-surface flows[END_REF] formulation was adopted. Only a few basic features of the formulation and implementation are presented here for reasons of completeness. Periodic conditions in the direction parallel to the substrate were employed. The laterally periodic domains bounded by the substrate and the moving interface were tessellated into elements of quadrilateral shape. Within each element, the velocities and solute concentration were represented by biquadratic shape functions, whereas the pressure was represented by bilinear functions (Taylor-Hood elements, [START_REF] Elman | Finite Elements and Fast Iterative Solvers[END_REF]). The interface was discretized by a set of points equi-spaced in the direction parallel to the substrate, thus representing the local film thickness. In the ALE formulation a grid velocity conforming to the boundary velocity needs to be found. This velocity was taken to have only a vertical component (in the direction towards the substrate) and a magnitude proportional to the vertical velocity of the interface and the vertical position of the grid point. Symbolically, if y is the coordinate of a grid point, h the local thickness and ht its rate of change, then the grid velocity is

𝒗 𝐺 = ℎ 𝑡 𝑦 ℎ 𝒆 𝑦 (9)
If φ is a test function from the solute concentration function space, the weak form of the solute transport equation in the ALE formulation may be written: [START_REF] Colinet | Nonlinear Dynamics of Surface-TensionDriven Instabilities[END_REF] To obtain the above form the diffusive terms have been integrated according to Green's theorem and the interfacial flux condition is employed to introduce the last term in the right-hand-side in place of the interfacial diffusive flux.

Similarly, if v is a test function from the velocity space, the weak form of the momentum equations reads [START_REF] Pearson | On convection cells induced by surface tension[END_REF] where, again, the interfacial stress condition has been employed to introduce the last term in place of the components of the stress tensor. Here, 𝑫𝒖 = 1/2(𝛁𝒖 + 𝛁𝒖 𝑻 ) is the rate of deformation tensor.

A second-order-accurate time integration scheme

A procedure to integrate the above equations together with the kinematic condition with 2nd-order accuracy in time is now briefly outlined. This is largely based on the discussion presented by Lucca and Formaggia [START_REF] Formaggia | Stability analysis of second-order time accurate schemes for ALE-FEM[END_REF], where more information may be found. After discretizing Eq. ( 10) in space according to the approach mentioned above, namely with the Taylor-Hood elements, but retaining the differential form in time one obtains the abstract form: [START_REF] Scriven | On cellular convection driven by surface-tension gradients: effect of mean surface tension and surface viscosity[END_REF] Here ω is the vector of discrete unknown concentrations, M, C, L are the so-called mass, convection, Laplace matrices, and S1 is a matrix containing the interfacial contributions. Since the domain is deforming, these matrices are functions of time. Similarly, the momentum equation may be written [START_REF] Reichenbach | Linear perturbation analysis of surface tension driven convection at a plane interface[END_REF] where u is the vector of unknown velocities, L2 is a matrix analogous to the Laplacian accounting for the viscous terms with variable viscosity, and G, S2 are the pressure gradient and the interfacial stress matrices, respectively.

Suppose now that we have calculated the film thickness, the concentration and the velocity field at an instance in time, say tn. The time advancement scheme starts with finding the new interface position at tn+1 = tn + dt. This is done by integrating the kinematic conditions with a 2 nd -order in time explicit scheme. For this purpose, the spatial derivatives of the film thickness are expressed by finite-differences and 3 rd -order upwinding is applied.

The next step is to advance the concentration field. Integrating equation ( 12) over the time interval (tn, tn+1) and evaluating the mass and diffusive terms implicitly and the rest of the terms explicitly, one may write with 2nd-order accuracy in time: [START_REF] Goussis | On the thermocapillary instabilities in a liquid layer heated from below[END_REF] Thus, the mass and Laplace matrices need to be evaluated at an intermediate level within the time step by interpolating between the domains at tn and tn+1.

The time advancement scheme is completed by solving Eq. ( 13) together with the continuity equation for the velocity un+1 and the pressure. Since no time derivatives are involved here and the concentration and the domain have already been advanced, all the terms refer to the time level tn+1 and there is no truncation error in time. The solution is obtained by a coupled Uzawa algorithm where the conjugate gradient method is employed to find the velocity corrections during the iterations. The curvature and normal vector to the interface that are needed in this substep are obtained by employing cubic splines.

Results and discussion

In the numerical simulations the functional relationship of viscosity vs. concentration was taken to be the same as obtained from the experimental measurements (cf. Fig. 1) and the initial polymer volume fraction was 0.05. The parameter range covered was M = 3000-40000, Ca = 0.05-0.5, while E was kept fixed at a value of 5 × 10 -4 . The duration of most simulations was up to 1000 dimensionless time units, which corresponds to the reduction of the film thickness to half of its original height and doubling of the average polymer concentration in the film. It was found that by that time, the velocities developed had slowed down significantly and the patterning of the solute had reached a quasi-steady distribution.

The above conditions were manageable in terms of computational resources and time, but significantly milder than the experimental ones, which range as M ~ 10 5 -2 × 10 6 , E ~ 10 -5 -10 -4 , and Ca ~ 0.01-0.02. However, under these conditions the high Marangoni numbers demand exceedingly fine spatial resolution, the small capillary numbers demand very small time steps because the interfacial deformation problem is numerically stiff, and, finally, the small evaporation numbers require simulations over extended time spans to cover the evolution of the evaporation process. It is therefore suggested that the results may serve to obtain a qualitative understanding of the main processes and mechanisms rather that to reproduce quantitatively the actual experiments. In this respect, it should be kept in mind that the simulations are 2-dimensional in space. For instance, in the experimental configuration corresponding to Fig. 2, wrinkling can be observed at about 25,000 dimensionless time units.

Free surface deformation 4.1.1. Effect of capillary and Marangoni numbers

In Fig. 5 a measure of surface deformation is presented as a function of time for three capillary numbers and various Marangoni numbers. The first and most important observation to be made is that deformation reaches an overall magnitude that is proportional to the capillary number. This parameter has therefore a dominant effect. At the lowest capillary number, the deformation does not exceed a few percent through the course of the simulation. The numerical results therefore are in accord with the intuitive expectation that strong capillary forces prevent interfacial deformation. Hence, the results suggest that film deformation cannot be the dominant reason of patterning and wrinkling. To illustrate the effects of capillary and Marangoni numbers, Fig. 6 shows the evolution of concentration contours through the course of drying for M = 10,000, Ca = 0.5 and 0.1 on the one hand, Ca = 0.1, M = 10,000 and M = 20,000 on the other hand. The immediate effect of evaporation is that vertical concentration gradients develop. This is simply the result of the selective evaporation that is promoting solvent depletion in the liquid near the interface and is initially counteracted only by bulk diffusion. At high Ca number, the small resistance of the interface to deformation promotes a rapid collection of the concentrated top layer and a significant deformation of the interface is observed with the formation of two lenses at t = 200, which are still visible at t = 1000. We also observe the appearance of secondary structures that are consumed continually by the larger and more permanent ones. This is a feature of Marangoni instability observed in several thermocapillary or solutocapillary systems [START_REF] Schatz | Experiments on thermocapillary instabilities[END_REF][START_REF] Schwarzenberger | Pattern formation and mass transfer under stationary solutal Marangoni instability[END_REF].

Decreasing Ca inhibits the interface deformation. The formation of plumes, as for example in the contours for M = 10,000 and t = 300, 400, entrain the concentrated solute from the interface back to the bulk. In this way mixing of the solute takes place, which, in a sense, reduces the driving force and counteracts the deformation process. A more detailed sequence of concentration contours for M = 10,000 in Fig. 7 shows this behavior more clearly. In the interval t = 180-210 the sudden growth in deformation observed in Fig. 5b is related to the nonlinear growth of the instability. Right after t = 220 the sudden drop clearly coincides with the initiation of plume formation. A similar effect is observed for M = 20,000 at an even earlier time (cf. Fig. 5b). A second effect is also evident by examining the sequence of concentration contours for two different Marangoni numbers, M = 10,000 and M = 20,000 in Fig. 6. At M = 20,000 finer cells appear at an earlier time, which again tend to cluster and then form plumes plunging back into the bulk. Now the motion is strong so that more plumes per unit length are able to survive the coarsening process. Therefore, the distance between the final cells is smaller and capillary forces may be more effective in resisting the deformation of the film. Compatible with the above picture is the result for the deformation in a simulation with constant viscosity, as shown in Fig. 5a. The fluid motion is not inhibited since the viscosity does not increase with solute concentration and hence the tendency of plume formation and the resistance to clustering are increased. Both these effects counteract film deformation.

Comparison with order of magnitude analysis

The order of magnitude analysis is based on previous results detailed in [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF][START_REF] Doumenc | Free convection in drying binary mixtures: solutal versus thermal instabilities[END_REF]. It corresponds to steady regime and constant viscosity assumptions. The additional assumptions and scaling derivation are detailed in Appendix A. The amplitude δh of free surface deformation reads

𝛿ℎ ~ 𝐶𝑎 𝑀 1/3 𝐸 2/3 (15)
The results of the simulations clearly support the conclusion that deformation decreases as the capillary number gets smaller, in agreement with the scaling law (Eq. ( 15)). Regarding the overall effect of the Marangoni number, a clear conclusion cannot be drawn because this would require more simulations as the exponent in the power law is only 1/3. However, it may be suggested that film deformation due to its overall small magnitude cannot be the dominant reason of patterning and wrinkling, especially as the capillary number in the experiments is less than 0.02.

Solute Distribution along the film 4.2.1. Effect of Marangoni number

We now turn to the distribution of the solute along the film and examine whether non-uniformities could account for film patterning. For this purpose, the simulations at the lower capillary number of 0.05 where film deformation is small are examined in more detail. The main features obtained from these simulations with increasing Marangoni number are highlighted in Fig. 8, where solute concentration contours are shown again at various instances during film drying. Note that due to increased numerical resolution demands and computational cost limitations, in these simulations the lateral periodic domain size was selected to be 2π. To facilitate comparison with similar plots (Figs. 6,7 and 9) two adjacent periodic cells are drawn for each case.

As may be seen, at the relatively low Marangoni number of 3000 the instability is barely able to develop and induce shearing and deformation of the thickening layer of concentrated solute near the interface. In this way this layer collects into lenses, thus making possible the non-uniformity in the solute distribution. As the Marangoni number is increased, finer cellular structures appear at an ever earlier stage and coalesce into big structures. Beyond that phase, within the large established structures smaller secondary ones continue to appear at large M almost throughout the entire simulation. These secondary structures become constantly consumed by the already established large permanent formations. The process is similar to what is seen in 3D simulations [START_REF] Kollner | Multiscale structures in solutal Marangoni convection: Three-dimensional simulations and supporting experiments[END_REF][START_REF] Schwarzenberger | Pattern formation and mass transfer under stationary solutal Marangoni instability[END_REF]. Up to M = 20,000 the number of the final cells remains the same, whereas at M = 40,000 their number increases. From all the above results it may be suggested that the increase in the strength of the Marangoni forces promotes the instability, the clustering of cells and the non-uniformity of solute distribution. However, again and in direct analogy with the discussion on film deformation there appear to be subtle effects for the same reasons. Stronger plumes result in enhanced mixing of the solute and finer structures are able to survive the process of coalescence. It may be argued here that since the structures finally remaining after cessation of coalescence are more than one in the periodic domains there is no artifact of the size of the lateral domain where periodic conditions are imposed. However, to further test such effects one simulation was repeated with the lateral domain size doubled and it was shown that the dynamics and number of final structures per unit lateral length are not affected.

Effects of the varying viscosity

To examine the effect of the varying viscosity on the distribution of solute along the film, similar simulations were performed where the viscosity was kept constant, independent of solute concentration. The capillary number and the evaporation number in these simulations were the same and the Marangoni number ranged between 5000 and 20,000. The development of the instability may be observed through the evolution of concentration contours during drying, as is shown in Fig. 9.

The main events are similar. A notable difference by comparing Figs. 8 and 9 is that the instability appears at an even earlier stage and the number of final structures is higher after the process of cell coalescence and coarsening. The differences with the variable viscosity simulations are more pronounced at relatively lower Marangoni numbers. Not unexpectedly, when the viscosity is constant the overall magnitude of the velocities is larger and it is sustained for a longer period during drying, as may be seen by comparing the data in Fig. 10. When the viscosity is varying the velocities are smaller, here by a factor in the range 2-5, which depends on the Marangoni number.

A related important difference is in the characteristics of the velocity field, a sample of which is shown in Fig. 11. The large arrows in this figure indicate interfacial points of maximum solute concentration. When the viscosity is constant, the downward motion is strong at the points corresponding to high concentrated regions, thus promoting mixing. The opposite happens when the viscosity increases with concentration. The downward motion almost vanishes at the points of maximum concentration.

The same kind of behavior was observed experimentally by Bassou and Rharbi during drying of PS/toluene solutions. Indeed, using trajectory of probe particles to visualize the flow, they obtained quasi immobile zones between the convective cells.

All the above features suggest that the effect of viscosity is significant in promoting the non-uniformity of the solute distribution along the film. In contrast to the constant viscosity case, where the higher velocities at the points corresponding to maximum concentration promote sustained mixing, the increasing viscosity inhibits convection and mixing of the solute back to the bulk. Indeed, this may also be observed from the solute concentration probability density distributions, shown in Fig. 12, all corresponding to dimensionless time, t = 1000, when roughly half of the solvent has been removed. For constant viscosity these distributions are symmetric around the mean value of 2 (relative to the initial concentration), they are relatively narrow and not significantly affected by the Marangoni number. In contrast, the distributions for the variable viscosity case are skewed and have extended tails towards values larger than the average. This is the result of asymmetry between upward and downward motion promoted by the increasing viscosity.

Patterning of the dried film due to solute lateral non-uniformity

Given all the previous observations, it appears plausible to suggest that the non-uniformity in the solute distribution is a significant effect and could be responsible for the final patterning of the dried films. During the late stages of drying evaporation would certainly become affected by concentration and mainly come from the more dilute parts of the interface. Motion would also be significantly slowed down and the diffusivity of the solvent would decrease dramatically. Whereas these effects may not be easily simulated, it may be argued that diffusion would not be able to restore the concentration uniformity in the remaining film because it would require a relative time scale much larger than O(M).

As a measure of solute lateral non-uniformity, the vertically averaged solute content was calculated as a function of position along the substrate and of time. This vertically averaged polymer content, which is scaled with the initial mean concentration, also represents a measure of the remaining dry film thickness if the solvent were instantaneously removed. The evolution of the standard deviation of this measure of dry thickness with time is shown in Fig. 13. It is clear from this figure that the distribution has reached a steady regime for t > 200 for the cases of constant viscosity. A very different behavior is observed for the cases of variable viscosity, where some regions become more concentrated while others stay diluted (the evaporation rate being constant in the simulations, the mean solute concentration is the same for variable or constant viscosity). No steady regime is observed for variable viscosity all along the drying, which results in an accumulation phenomenon inducing spatial variation of polymer concentration much larger than in the constant viscosity case.

If it is assumed, therefore, that the simulations capture, at least qualitatively, the final solute distribution, then an estimate of the final film thickness non-uniformity may be obtained by the distribution of the vertically averaged solute content at the end of the simulations. The relevant probability density distributions at t = 1000 along the film are shown in Fig. 14. Compared to the constant viscosity case these distributions are wider and skewed towards values smaller than the average.

The trends are also evident by examination of some statistical parameters for various Marangoni numbers, as shown in Table 1. Notably, excursions above the mean are more pronounced and the effect of the varying viscosity is clear.

If we therefore plot the projected final film thickness along the interface we obtain the results shown in Fig. 15. Significant variations appear in the final thickness, of the order of the mean value, and this occurs only when the viscosity is taken to vary with the solute concentration. This is in agreement with orders of magnitude obtained in experiments. Indeed, for PS/toluene solutions, Bassou and Rharbi obtain about 0.25 for the ratio of the average corrugation amplitude to the average dry film thickness. For the PIB/toluene experiments described in Section 2, standard deviation is about 0.28.

Concluding remarks

The purpose of the present paper is to shed light through numerical simulations on the potential mechanisms that cause permanent wrinkling of dried polymer films. Solutal Marangoni phenomena are considered. First, simulations show that free surface deformation induced by capillary pressure reaches an overall magnitude that is proportional to the capillary number, in agreement with scaling analysis. Decreasing Ca inhibits the interface deformation through the formation of plumes that helps mixing of the solute. For the value of capillary number corresponding to experiments, free surface deformation is not significant. The other mechanism investigated is the non-uniformity of solute concentration induced by Bénard-Marangoni convection. Comparisons of simulations with constant and varying viscosity suggest that the effect of viscosity is significant in promoting the non-uniformity of the solute distribution along the film and leads to accumulation of solute that is observed up to the end of the simulation. When the solute concentration is high enough, velocity decreases due to high viscosity and the system is frozen. Resulting corrugation amplitudes are in agreement with experiments. Therefore, though simulations have been performed with dimensionless numbers significantly milder than the experimental ones, they confirm the mechanism proposed by Bassou and Rharbi [START_REF] Bassou | Role of Bénard-Marangoni Instabilities during Solvent Evaporation in Polymer Surface Corrugations[END_REF], i.e. wrinkles may be viewed as the signature of Bénard-Marangoni convection.

It can be shown that the configurations used in the numerical simulations performed in this paper belong to the domain called ''B'' in the scaling analysis performed in [START_REF] Doumenc | Free convection in drying binary mixtures: solutal versus thermal instabilities[END_REF] (cf. Table 3 and Fig. 11 in this reference). Indeed, with the definition of Marangoni and evaporation numbers used in this paper, this region is limited by the following boundaries:

1 << M 2 E << Sc 3/2 and M E 2 << 1. It is easy to check that the values used in the simulations (M = 3000-40,000, E = 5×10 -4 , Sc = 2×10 5 ) obey to these three criteria. In this domain, the hydrodynamic layer is of the order of one and the velocity scales as M 1/3 E 2/3 [START_REF] Doumenc | Free convection in drying binary mixtures: solutal versus thermal instabilities[END_REF].

A.1. Inertia versus viscous terms

Using the above results, it is possible to compare the inertia and viscous terms in the Navier-Stokes equation (1) governing the fluid motion.

u and v being the velocity components in the horizontal and vertical directions, mass conservation with the above velocity scaling reads 𝑢 ~ 𝑣 ~ 𝑀 1/3 𝐸 2/3 (A-1)

In the steady regime and assuming constant viscosity, the ratio between inertia terms (lhs of Eq. ( 1)) and viscous terms (rhs of Eq.( 1)) is therefore of the order of M 4/3 E 2/3 Sc -1 .

According to Fig. 11 and Table 4 of [START_REF] Doumenc | Free convection in drying binary mixtures: solutal versus thermal instabilities[END_REF], with the definition of Marangoni number used here, ''B'' region verifies M 2 E << Sc 3/2 . Thus, the ratio between inertia terms and viscous terms is much less than one, which means that viscous terms are dominant for the configurations studied.

A.2. Scaling analysis for free surface deformation

Free surface deformation is assumed small and small slope approximation is used in the following. No distinction is made between the vertical direction and the normal to the free surface. The curvature is taken as the second derivative of the thickness and is of the order of δh / λ 2 , where δh and λ are the amplitude and wavelength of deformation, respectively. Thus, normal stress balance at the interface reads in dimensionless form ).

The hydrodynamic layer being of the order of one in the ''B'' domain, and using Eq. (A.1) and assumption H.4, the order of magnitude for the normal stress reduces to 𝛤 𝑦𝑦 ~ 𝑢 ~ 𝑀 1/3 𝐸 2/3 . (A-5)

From (A-4) and (A-5) it follows:

𝛿ℎ ~ 𝐶𝑎 𝑀 1/3 𝐸 2/3 .

Fig. 1 .

 1 Fig. 1. Viscosity of PIB/toluene solution as a function of PIB volume fraction. Measurements (triangles, performed with a Low Shear 30 rheometer except for the viscosity of pure PIB) and empirical law used in the simulations (continuous line).

Fig. 2 .

 2 Fig. 2. Wrinkling of the film during the transition from the constant flux regime (t < 7500s) to the last regime (t > 11,000 s). Time evolution of the evaporative flux (continuous line, left axis) and time evolution of the percentage of wavy surface (diamond markers, right axis). After t = 11,000 s, wrinkling is frozen. The remaining flux, controlled by solvent diffusion in the concentrated solution, is very weak. (Configuration: initial thickness = 4 mm, initial polymer mass fraction = 5%, evaporation rate = 0.4 µm/s).

Fig. 3 .

 3 Fig. 3. Wrinkled dry polymer film. Image size is 40 mm × 40 mm. Color scale corresponds to film height, from blue color (holes) to red color (bumps). Zero level corresponds to the mean plane. (Configuration: initial thickness = 1.06 mm, initial polymer mass fraction = 5.07%, evaporation rate = 0.36 µm/s).

Fig. 4 .

 4 Fig. 4. Dry film profiles for 0 < y < 40 mm and, from top to bottom, x = 10 mm, 20 mm, 30 mm (positions corresponding to dashed vertical lines in Fig.3.).

Fig. 5 .

 5 Fig. 5. Effect of the capillary number on interfacial deformation as of function of time for Ca = 0.5 (a), Ca = 0.1 (b) and Ca = 0.05 (c) and various Marangoni numbers.

Fig. 6 .

 6 Fig. 6. Successive snapshots of solute concentration contours. The lateral periodic domain size is 4π. In the arbitrary chromatic scale shown on top, 1 corresponds to the instantaneous maximum concentration and 0 to the minimum.

Fig. 7 .

 7 Fig. 7. Solute concentration contours for M = 10,000 and Ca = 0.1, showing initiation of plumes. Successive snapshots at time instants 180, 190, ... ,270. The lateral periodic domain size is 4π.

Fig. 8 .

 8 Fig. 8. Successive snapshots of solute concentration contours at time instants 100, 200, ... ,1000. The lateral domain size shown is 4π.

Fig. 9 .

 9 Fig. 9. Solute concentration contours for constant viscosity and M = 3000 (left), M = 10,000 (right). Successive snapshots at time instants 100, 200, ... ,1000. The lateral periodic domain size is 4π.

Fig. 10 .

 10 Fig. 10. Maximum velocity for constant (left) and varying viscosity (right) as a function of time for various Marangoni numbers.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Map of instantaneous velocity vectors for Ca = 0.05, M = 5000 at t = 1000, for constant (top) and varying viscosity (bottom). The lateral periodic domain size is 4π.

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Standard deviation of the vertically averaged solute content as a function of time for constant (left) and variable viscosity (right), for various Marangoni numbers, Ca = 0.05.

Fig. 15 .

 15 Fig. 15. Thickness of the dried film on the basis of the solute lateral distribution for varying (top) and constant viscosity (bottom) and various Marangoni numbers, at time t = 1000 and Ca = 0.05.

2 (

 2 𝛤 𝑛𝑛 ~ (𝐶𝑎 -1 + 𝜔) 𝛿ℎ 𝜆 A-3)with 𝛤 𝑛𝑛 the normal stress at the interface. Using 𝛤 𝑛𝑛 ~ 𝛤 𝑦𝑦 , with assumption H4 (λ ~ 1) and in the range of capillary numbers used in the simulations, Eq. (A-3) reads𝛤 𝑦𝑦 ~ 𝐶𝑎 -1 𝛿ℎ (A-4)Disregarding inertia terms in Navier-Stokes equation, Eq. (1) in the vertical direction reads

Table 1 .

 1 Statistical parameters of the dried film on the basis of the solute lateral distribution. The data are normalized to have a mean value of unity for the dry film thickness hdry.

	M		Dry film thickness for variable viscosity		Dry film thickness for constant viscosity
		hdry-min	hdry-max	Standard deviation	hdry-min	hdry-max	Standard deviation
	3000	0.79	1.19	0.14	0.90	1.11	0.050
	5000	0.75	1.33	0.17	0.89	1.14	0.050
	10,000	0.73	1.53	0.22	0.86	1.09	0.049
	20,000	0.72	1.54	0.23	0.85	1.06	0.045
	40,000	0.73	1.40	0.16

Acknowledgements

F.D. and B.G. thank ANR EVAPEC 13-BS09-0010-01 for funding. S.G.Y. thanks the University Paris-Sud for funding as a visiting professor.

Conflict of interest

None declared.

Appendix A.

Based on simplifying assumptions, scaling laws have been obtained by ''solving'' equations of solutal model in terms of order of magnitude. All the equations below must be hence understood as orders of magnitude only since these scalings do not give prefactors. We use results previously obtained in [START_REF] Trouette | Transient Rayleigh-Bénard-Marangoni solutal convection[END_REF][START_REF] Doumenc | Free convection in drying binary mixtures: solutal versus thermal instabilities[END_REF], where scaling analysis on the solutal problem led to a partition in five domains of the different configurations, characterized by the values of the non-dimensional parameters. For each domain, the thicknesses of hydrodynamic boundary layer and solutal boundary layer, the velocity and concentration variation have been derived from scaling analysis and expressed as a function of Marangoni and Evaporation numbers. Based on these results, we analyze here the relative importance of inertia and viscous term, and derive the scaling law for free surface deformation. The following assumptions are used: H1. Below the free surface, we assume the existence of a hydrodynamic and a solutal boundary layer. H2. Time derivative terms are neglected (quasi-steady regime). H3. Solvent concentration variations across a convective cell in the vertical and horizontal directions are of the same order of magnitude. H4. The wavelength of convective structures scales with the layer thickness. Hence the order of magnitude of the dimensionless characteristic length in the horizontal direction is one. H5. The analysis is restricted to fluids characterized by a Schmidt number greater than one. H6. Viscosity and layer thickness are assumed constant and equal to their initial values.


Laure Tosatto
email: laure.tosatto@univ-amu.fr

Joël Fagot
email: joel.fagot@univ-amu.fr

Arnaud Rey
email: arnaud.rey@cnrs.fr

Sequence Learning and Chunk Stability in Guinea Baboons (Papio papio) Apprentissage de séquence et stabilité des chunks chez les babouins de Guinée (Papio papio)

Keywords: sequence learning, associative learning, chunking

Chunking mechanisms, the processes of grouping several items together into a single processing unit, are central to several cognitive processes in human and non-human primates and notably to the acquisition of visuomotor sequences. Individuals segment sequences into chunks to perform visuomotor tasks more fluidly, rapidly, and accurately. Using an operant conditioning device, we previously studied the precise mechanisms by which chunks are formed and reorganized during sequence learning. Eighteen Guinea baboons (Papio papio) repeatedly produced the same fixed sequence of nine movements during 1,000 trials by pointing to a moving target on a touch screen. We found that chunking patterns are reorganized during the course of learning, with chunks becoming progressively fewer and longer. We also identified two forms of reorganization of the chunking pattern: the recombination of preexisting chunks and the concatenation of two distinct chunks into a single one. To understand the conditions under which these reorganizations occur, we study here how the stability of a chunk and the stability of chunk boundaries are related to these reorganizations. Our analyses showed that less stable chunks and less stable boundaries are more likely to produce reorganizations. These results provide new evidence about the finegrained dynamics of chunking mechanisms during sequence learning.

Résumé

Les mécanismes de chunking, processus par lesquels plusieurs items sont regroupés pour former une unité fonctionnelle, sont centraux dans de nombreux processus cognitifs chez les primates humains et non-humains, et particulièrement lors de l'acquisition de séquences visuomotrices. Les individus segmentent les séquences en chunks pour réaliser les tâches visuomotrices plus fluidement, rapidement et précisément. L'utilisation d'un système de conditionnement automatique nous a précédemment permis d'étudier les mécanismes précis par lesquels ces chunks sont formés et réorganisés durant l'apprentissage de séquences. Dix-huit babouins de Guinée (Papio papio) ont répété la même séquence fixe de neuf mouvements durant 1000 essais lors d'une tâche de pointage d'une cible mouvante sur écran tactile. Nous avons constaté que les patterns de chunking sont réorganisés au cours de l'apprentissage, les chunks devenant progressivement moins nombreux et plus longs. Nous avons également identifié deux formes de réorganisation du pattern de chunking : la recombinaison de chunks préexistants et la concaténation de deux chunks distincts en un seul.

Pour comprendre les conditions dans lesquelles ces réorganisations se produisent, nous étudions ici comment la stabilité d'un chunk et la stabilité des frontières entre chunks sont liées à ces réorganisations. Nos analyses montrent que les chunks et les frontières moins stables sont les plus susceptibles de produire des réorganisations. Ces résultats apportent de nouvelles informations sur la dynamique fine des mécanismes de chunking au cours de l'apprentissage de séquences.

Mots-clés : apprentissage de séquences ; apprentissage associatif ; chunking

Introduction

As individuals in constant interaction with the environment, primates evolve in a continuous stream of stimulations that trigger behavioral responses. Most behaviors are produced sequentially and learning rapidly and efficiently long sequences of actions represents an evolutionary challenge. This is the case for humans, especially for language processing and production [START_REF] Conway | Sequential learning in non-human primates[END_REF][START_REF] Lashley | The Problem of Serial Order in Behavior[END_REF][START_REF] Petkov | Structured Sequence Learning: Animal Abilities, Cognitive Operations, and Language Evolution[END_REF], but more broadly in all primate communication systems, where sequences can be long and complex [START_REF] Aychet | Sequential and network analyses to describe multiple signal use in captive mangabeys[END_REF]Habib-Dassetto, Montant & Reboul, this issue;Ouattara, Lemasson & Zuberbülher, 2009).

An important constraint rendering the learning of long sequences difficult is the natural limits of short-term memory [START_REF] Cowan | Evolving conceptions of memory storage, selective attention, and their mutual constraints within the human information-processing system[END_REF][START_REF] Cowan | The many faces of working memory and short-term storage[END_REF]. To bypass these limits, a core cognitive mechanism allowing the compression of information, and increasing short term memory capacity, is the formation of chunks [START_REF] Mathy | What's magic about magic numbers? Chunking and data compression in short-term memory[END_REF][START_REF] Miller | The Magical Number Seven, Plus or Minus Two: Some Limits on Our Capacity for Processing Information[END_REF]. Indeed, chunking is usually defined as the process of associating and grouping several items together into a single processing unit [START_REF] Gobet | Chunking mechanisms in human learning[END_REF][START_REF] Gobet | What's in a Name? The Multiple Meanings of "Chunk" and "Chunking[END_REF].

To understand the general learning mechanisms shared by primates and involved in sequence learning and chunking, motor sequence learning tasks have been widely used in comparative psychology, where motor sequence learning is defined as the process by which a specific sequence of movements is executed with increased speed and accuracy [START_REF] Willingham | A neuropsychological theory of motor skill learning[END_REF]. In this field, chunking is widely considered as the main motor sequence integration mechanism [START_REF] Diedrichsen | Motor skill learning between selection and execution[END_REF]Wymbs et al., 2012), and is commonly believed to be rooted in elementary associative processes (i.e., [START_REF] Perruchet | The self-organizing consciousness as an alternative model of the mind[END_REF][START_REF] Rey | Regularity Extraction Across Species: Associative Learning Mechanisms Shared by Human and Non-Human Primates[END_REF]. Indeed, several primate species have been shown to spontaneously segment sequences in chunks of 3 or 4 items as indicated by long temporal gaps emerging between successive responses and marking chunk boundaries (e.g., in humans, Homo sapiens, [START_REF] Abrahamse | Control of automated behavior: Insights from the discrete sequence production task[END_REF][START_REF] Bottary | Insufficient chunk concatenation may underlie changes in sleep-dependent consolidation of motor sequence learning in older adults[END_REF]in rhesus monkeys, Macaca mulatta, Scarf et al., 2018[START_REF] Terrace | The Comparative Psychology of Chunking[END_REF]in baboons, Papio papio, Tosatto et al., 2022). Furthermore, as the input's representation might evolve during learning, some studies were also interested in the temporal aspect of chunking and how chunking patterns change with practice. Throughout extended practice, chunks were found, in human and non-human primates, to evolve and grow larger as if more compression of information was possible with increasing familiarity with the sequence (e.g., [START_REF] Acuna | Multifaceted aspects of chunking enable robust algorithms[END_REF][START_REF] Bera | Motor Chunking in Internally Guided Sequencing[END_REF][START_REF] Ramkumar | Chunking as the result of an efficiency computation trade-off[END_REF][START_REF] Wright | Offline Improvement during Motor Sequence Learning Is Not Restricted to Developing Motor Chunks[END_REF]. One mechanism proposed to account for this growth is the concatenation of chunks [START_REF] Verwey | Concatenating familiar movement sequences: The versatile cognitive processor[END_REF][START_REF] Wright | Offline Improvement during Motor Sequence Learning Is Not Restricted to Developing Motor Chunks[END_REF], i.e., independent chunks being executed more fluidly with practice and with a decrease of the temporal gap between them leading to a single and longer chunk [START_REF] Abrahamse | Control of automated behavior: Insights from the discrete sequence production task[END_REF].

We have recently conducted a study on Guinea baboons (Papio papio) on the role of extended practice in the formation and the evolution of chunks [START_REF] Tosatto | The Evolution of Chunks in Sequence Learning[END_REF]. The task was a serial response time (SRT) task where baboons had to point to a moving target on a touchscreen and were repeatedly exposed to the same sequence composed of 9 different locations during 1,000 trials. In accordance with the literature, this study found an increase in chunks size throughout the repeated production of the sequence. However, we also observed two (instead of one) forms of reorganizations of the chunking pattern. Indeed, as in previous studies, we observed concatenations, but we also observed another mechanism, recombinations, i.e., the emergence of a new segmentation pattern -such as two chunks of 3 items becoming a chunk of 4 items followed by a chunk of 2 items. Thus, recombinations seem to also lead to an increase in chunk size but necessitate a previous chunk to be resegmented. These findings show that linear associations between adjacent elements or a hierarchical organization of the sequence do not seem to be the only mechanisms at work during sequence learning, and that chunks may be more flexible than previously thought [START_REF] Diedrichsen | Motor skill learning between selection and execution[END_REF][START_REF] Rosenbaum | Hierarchical control of rapid movement sequences[END_REF][START_REF] Sakai | Chunking during human visuomotor sequence learning[END_REF][START_REF] Fonollosa | Learning of Chunking Sequences in Cognition and Behavior[END_REF].

To explain this variability in the types of reorganizations observed, we hypothesized the role of two main parameters: (1) the internal stability of chunks; (2) the stability of chunk boundaries. Indeed, stable chunks might be less likely to be resegmented from one block to another and survive throughout the task, whereas a less stable chunk might be more likely to disappear (see the Results Section for an operational definition of chunk stability and chunk boundary stability). Similarly, a stable chunk boundary might be less likely to disappear throughout the task than a less stable one.

In this study, we propose to enrich our understanding of these reorganization phenomena by suggesting a method to assess the internal stability of a chunk and the stability of a boundary between two chunks. Considering the large dataset provided by [START_REF] Tosatto | The Evolution of Chunks in Sequence Learning[END_REF] and the many occurrences of chunk reorganizations it contains, using dedicated posthoc analyses, we attempt at shedding light on the link between these parameters and the evolutions of the chunking pattern throughout this experiment.

Method

Fully detailed method and materials can be found in [START_REF] Tosatto | The Evolution of Chunks in Sequence Learning[END_REF].

Participants

Thirteen female and five male Guinea baboons (Papio papio, age range 2.8-23.7 years) from the CNRS primate facility in Rousset (France), living in a social group of 25 individuals, were tested in this study. For practical reasons, we stopped the experiment after

Materials

This experiment was conducted with a computer-learning device based on the voluntary participation of baboons (for details, see [START_REF] Fagot | Automated testing of cognitive performance in monkeys: Use of a battery of computerized test systems by a troop of semi-free-ranging baboons (Papio papio)[END_REF]. Baboons had free access to operant conditioning learning devices equipped with touch screens. The screen was divided into nine equidistant predetermined locations represented by white crosses on a black background, virtually labeled as position 1 to 9 (see Figure 1).

<insert Figure 1 here> A trial began with the presentation of a yellow fixation cross at the bottom of the screen. Once pressed, the fixation cross disappeared and the nine white crosses were displayed, one of them being replaced by the target, a red circle. When the target was touched, it was immediately replaced by the cross. The red circle then replaced the next position in the sequence until it was touched, and a new position was displayed. Reward was provided at the end of a sequence of nine touches. The time elapsed between the appearance of the red circle and the baboon's touch on this circle was recorded as the RT.

Procedure

Baboons were either presented with a single sequence (each baboon was either trained on Sequence 1, N=8, or Sequence 2, N=10, and never received the other sequence) and had to produce it for 1000 successive trials. RTs for each position of the sequence were recorded for all the trials.

Results

Previous results summary

To observe the chunking pattern at different stages of learning, RTs for each of the nine positions and for the 1,000 trials were divided into 10 Blocks of 100 trials, and we identified the chunking pattern for each baboon at each block by computing mean RTs per position over each block of 100 trials. The identification of chunks relied on identifying chunk boundaries by isolating significant increases in the successive responses to each position in the sequence. Practically, for each block of 100 trials, the mean RT for Position n and n+1 were computed and if the mean RT for Position n+1 was significantly higher than the mean RT for Position n, then this difference was taken as marking a boundary between two chunks. More details about this procedure can be found in [START_REF] Tosatto | The Evolution of Chunks in Sequence Learning[END_REF]. Figure 2 illustrates the chunking pattern of a single baboon (Ewine) and its evolution throughout the task.

<insert Figure 2 here> Three main findings were obtained in this previous study. First, in line with the literature, baboons spontaneously segmented the long sequences into short chunks at the beginning of the task. Second, with extended practice, chunks became longer and fewer.

Third, the number of chunks and the increase in chunks' size was due to two types of reorganizations: the recombination of several preexisting chunks and the concatenation of two distinct chunks into a single one. For instance, in Figure 2, at Block 2, Positions 6 to 9 are grouped into two chunks of two items each (6-7/8-9). However, at Block 3, a new segmentation pattern has emerged: Position 6 is no longer grouped with Position 7, as the mean RT for the latter has significantly increased. Position 7 and 8 are now grouped into a chunk, while the mean RT for Position 9 has increased and it is no longer grouped with Position 8. The change in chunking pattern for these Positions between Blocks 2 and 3 is an example of recombination. An example of concatenation can be found on the same figure between Blocks 9 and 10 for Positions 1 to 5. During Block 9, these positions are grouped as two distinct chunks of respectively three and two items. At Block 10, the mean RTs for Positions 4 and 5 decrease, making the boundary between Positions 3 and 4 disappear and leading to a single chunk of five items. Finally, these evolutions of the chunking pattern were accompanied by a significant decrease of the mean response time in all baboons. This result suggests that the formation of chunks and their modification accelerated sequence production, confirming that chunking fosters learning and increases performance fluency.

Post-hoc analyses

Proxies for inferring chunk and boundary stability

Our main goal here is to acquire a better understanding of chunks reorganizations and why chunks either survive from block to block or are reorganized. We hypothesize that the stability of a chunk plays a key role in its survival. One way of looking at the stability of a chunk is to consider the decrease in RT between the successive positions within this chunk. Indeed, if we observe a strong decrease in the successive response times within a chunk, we can qualify this chunk as very stable. For example, if we look at the performance of the baboon Ewine in Block 4 (Figure 2), it appears that in the first chunk (grouping Positions 1-3), the decrease in RTs from one position to the next is relatively strong while the decrease in RTs for the second chunk (grouping Positions 4 and 5) is lower. Following our reasoning, the probability of the first chunk being reorganized would therefore be lower than for the second chunk.

Practically, the large difference in milliseconds between each response time within a chunk will result in a large standard deviation among all RTs. We thus propose to take the within standard deviation (WSD) of responses belonging to the same chunk as a proxy of its stability and assess the link between this value and the survival of the corresponding chunk throughout the blocks of trials. Similarly, we hypothesize that the stability of a boundary between two chunks strongly depends on the difference between the mean RT on the last position of the first chunk and the mean RT on the first position of the next chunk. The larger CHUNK STABILITY IN BABOONS 10 the difference, the stronger the stability of the boundary and the greater the probability the boundary will survive from Block to Block.

Figure 3 provides a schematic illustration of these propositions for a sequence of nine positions. Here the sequence is artificially segmented into three chunks of three positions. The first chunk includes Positions 1-3 and due to the large differences in RTs between each of these positions, the resulting WSD is also quite large, suggesting that the chunk is rather stable. Similarly, the difference between the mean RT in Position 3 (last position of the first chunk) and the mean RT in Position 4 (first position of the second chunk) is large indicating that the chunk boundary is quite stable. Due to smaller WSDs, the second and third chunks are certainly less stable and the same goes for their boundary with a smaller difference between Positions 6 and 7. We therefore predict that the last two chunks should give rise to reorganizations, while the first chunk should not change over the course of repetitions.

<insert Figure 3 here>

Chunk stability

To quantify the survival of a chunk from one block to the next, we used the following method: if a chunk is present in Block n, it is considered as surviving in Block n+1 if either the chunk is exactly the same in Block n+1 or it belongs to a larger chunk. We assessed the effect of WSD on chunk survival using a logistic binomial regression model in which chunk survival is coded 1 and 0 otherwise. Chunks composed of only one position were excluded from this analysis considering they necessarily have a WSD of 0 and a 100% chance of surviving from one block to the next.

The logistic regression was highly significant X 2 (1, N= 421) = 20.8, p< .001 and indicated that the likelihood of chunk survival increased with stability, and that each millisecond added to a chunk's WSD increased its chances of survival in the next block by 3% (OR=1.030; 95% CI [1.015, 1.045], see Figure 4A).

Boundary stability

Here again, we propose that the stability of a chunk boundary is simply reflected by the difference between the mean RT for the last position of the first chunk and the mean RT for the first position of the second chunk. The greater the difference, the greater the probability of survival of the boundary in the next Block should be. To quantify the boundary survivals from one block to the next, we simply assessed for each boundary present at Block n whether it was still marking a chunk boundary at Block n+1.

We used a second logistic regression to assess the relationship between the stability of the chunk boundary and its survival (survival=1; 0 otherwise) in the next block of trial. The logistic regression was highly significant X 2 (1, N= 407) = 70.4, p< .001 and indicated that likelihood of boundary survival increased with its stability, and that each millisecond added to the size of the boundary increased its chances of survival in the next block by 5.7% (OR=1.057; 95% CI [1.038, 1.076], see Figure 4B).

<insert Figure 4 here>

Discussion

There were two main findings in the present study. First, we found that there is a strong link between the internal stability of a chunk and its survival over different stages of learning. Indeed, very stable chunks, as evidenced by their large within standard deviation, are more likely to survive from block to block during a sequence learning task. Conversely, less stable chunks are more likely to be reorganized differently from block to block. Second, we found that the stability of a boundary between two chunks is also strongly linked with its survival. Similarly to chunk stability, very stable boundaries, as evidenced by their larger size, are more likely to survive from block to block. Therefore, both the internal stability of chunks and the stability of chunk boundaries likely determine the survival of a chunking pattern or its reorganization.

These results, although arising from post-hoc analyses, suggest that the role of these parameters -chunk stability and boundary stability -need to be studied further to understand the complex dynamics underlying the evolution of chunking patterns during sequence learning. However, we are not yet able to precisely link these reorganizations -recombination or concatenation -to the variations of these parameters. Indeed, for instance, chunk stability can predict chunk survival, but a chunk can survive (e.g., not being segmented) while still being concatenated with another chunk. Therefore, chunk survival does not indicate the absence of reorganization. Similarly, a boundary can disappear between two chunks and lead to a concatenation, but a new boundary can appear simultaneously on the same chunk and a recombination will be observed. This is the case in Figure 2, between Blocks 2 and 3, where the disappearance of the boundary between Positions 7 and 8 would have led to a concatenation, but two new boundaries have appeared between Positions 6 and 7 and Positions 8 and 9. Therefore, our predicted outcomes (chunk survival and boundary survival) are not themselves predictors of a particular type of reorganization. Understanding what determines either type of reorganizations most likely requires studying the interaction between these parameters (chunk and boundary stability), which we were not able to do on this limited dataset.

Although these parameters seem to clearly play a role in the evolution of the chunking pattern, modeling the evolution of this dynamical system still remains to be done. One way to further explore the role of these parameters and their interaction would be to design an experiment where chunk stability and boundary stability would be manipulated. We could CHUNK STABILITY IN BABOONS 13 test, for example, a situation where two successive chunks have a strong internal stability and are separated by a stable or less stable boundary. As the realization of the second chunk would become more efficient, leading to a lower difference between the last RT of the first chunk and the first RT of the second chunk, we would then predict a higher probability for these chunks to concatenate into a single one, depending on the boundary stability.

Additionally, it is important to note that, despite a variability of the chunking patterns across subjects (see [START_REF] Tosatto | The Evolution of Chunks in Sequence Learning[END_REF], we have identified a central common feature driving the reorganizations of chunks during sequence learning. Furthermore, we argue that these results, obtained with Guinea baboons, can confidently be generalized to many primate species. Indeed, as stated before, chunking mechanisms appear to be shared across primate species and rely on common associative mechanisms, and the precise dynamics of chunk evolution, including concatenations and recombinations, appear to be very similar even between humans and baboons [START_REF] Tosatto | The dynamics of chunking in humans (Homo sapiens) and Guinea baboons (Papio papio)[END_REF].

Note that, as the sequence is mastered, RTs may decrease to a floor performance which should lead to less possibilities for reorganizations. Floor performance for RTs was not reached in the present study but reaching floor performances indeed affect the chunking pattern and its reorganizations. In another experiment, we trained baboons over a longer period of 4,000 trials on a single repeated sequence composed of 5 positions (Tosatto et al., under review). We observed that baboons reached a plateau in terms of mean sequence production speed and, concurrently, a plateau was reached in the growth of chunk size, necessarily decreasing the number of recombinations and concatenations. The same effect was observed in humans, although they needed a shorter training period to reach this plateau [START_REF] Tosatto | The dynamics of chunking in humans (Homo sapiens) and Guinea baboons (Papio papio)[END_REF]. These results suggest that as practice increases, the chunking pattern seems to stabilize and the number of reorganizations decreases.

Finally, the fact that chunking patterns appear flexible and dependent on the stability of associations makes sense from an evolutionary point of view. Indeed, it is important that, during sequence learning and especially motor sequence learning, the cognitive system can maintain a balance between a stable representation of the sequence, to identify or produce it efficiently, and possibilities to reorganize said sequence in memory, to optimize performance and accommodate the behavior to the potential presence of interference (i.e., the presence of an event interrupting the sequence). Likewise, the analyses carried out in this study provide new constraints for modeling the nature and dynamics of fundamental chunking mechanisms with ever greater precision.

  Wymbs NF, Bassett DS, Mucha PJ, Porter MA, Grafton ST. 2012. Differential Recruitment of the Sensorimotor Putamen and Frontoparietal Cortex during Motor Chunking in Humans. Neuron 74: 936-946.

Figure 1 Figure 2

 12 Figure 1

Figure 3

 3 Figure 3

Figure 4

 4 Figure 4

Acknowledgements

This work was supported by the Institut Convergence ILCB (ANR-16-CONV-0002) and the HEBBIAN ANR project (#ANR-23-CE28-0008). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Animal rights

This research adhered to the applicable French rules for ethical treatment of research animals and received ethical approval from the French Ministry of Education (approval APAFIS#2717-2015111708173794 10 v3).

Data from the experiment are available on Open Science Framework at https://osf.io/xcw95/.

Open Practices Statements

Conflict of Interest

The authors declare that they have no conflict of interest.