Renaud Ansart
email: renaud.ansart@ensiacet.fr

Pablo García-Triñanes

Benjamin Boissière

Hadrien Benoit

Jonathan P K Seville

Olivier Simonin

Dense gas-particle suspension upward flow used as heat transfer fluid in solar receiver: PEPT experiments and 3D numerical simulations

Keywords:

A dense particle suspension, also called an upflow bubbling fluidized bed, is an innovative alternative to the heat transfer fluids commonly used in concentrated solar power plants. An additional advantage of this technology is that it allows for direct thermal storage due to the large heat capacity and maximum temperature of the particle suspension. The key to the proposed process is the effective heat transfer from the solar heated surfaces to the heat transfer fluid, i.e. the circulating solid suspension. In order to better understand the process and to optimise the design of the solar receiver, it is of paramount importance to know how particles behave inside the bundle of small tubes. To access to the particle motion in the solar receiver, two different techniques are carried out: experimental using positron emission particle tracking (PEPT) and 3D numerical simulation via an Eulerian n-fluid approach with NEPTUNE_CFD code. Both numerical predictions and PEPT measurements describe an upward flow at the centre of the transport tube with a back-mixing flow near the wall which influences the heat transfer mechanism. Comparisons between experiment and computation were carried out for the radial profiles of the solid volume fraction, and vertical and radial time-averaged and variance velocities of solid, and demonstrating the capability of NEPTUNE_CFD code to simulate this peculiar upflow bubbling fluidized bed.

Introduction

Conventional Heat Transfer Fluids (HTF) used in solar power plants have many drawbacks, in particular a limited working temperature domain. A solution to overcome these drawbacks uses solid particles as the HTF. In a new patented concept [START_REF] Flamant | PCT Extension[END_REF], it is proposed to use a dense particle suspension, also called an upflow bubbling fluidized bed, as the HTF, flowing upwards in an array of tubes at the focus of the solar receiver. The solids can also be used as an energy storage medium. The present study reports work carried out within a FP7 EC project CSP2, http://www.csp2-project.eu, aimed at developing this alternative HTF in order to extend the working temperature. The use of high temperature heat transfer fluids allows operating at high temperature thermodynamic cycles such as super critical and combined cycles in associated steam turbines. Dense particle suspensions will enable operating temperatures over 1000 • C limited only by the construction materials of the receivers, compared with 560 • C for the most efficient molten salts currently used, thus increasing the plant efficiency and decreasing the cost per produced kWh. In addition, this new HTF has no lower temperature limitation. Flamant et al. [START_REF] Flamant | Dense suspension of solid particles as a new heat transfer fluid for concentrated solar thermal plants: On-sun proof of concept[END_REF] have demonstrated the capacity of dense gas-particle suspensions to transfer concentrated solar power from a tubular receiver to an energy conversion process by acting as a HTF. A mean value of heat transfer coefficient corresponding to 400 W/m 2 K has been measured for standard operating conditions at low temperature and much higher values, up to 1100 W/m 2 K, were measured at temperatures up to 750 • C [START_REF] Benoit | On-sun demonstration of a 750 • C heat transfer fluid for concentrating solar systems: dense particle suspension in tube[END_REF]. The solar receiver consists of one or more multi-tube exchangers, which are the absorbing modules for the solar radiation. The number of exchangers depends on the power of the receiver, and also on the heliostat field configuration. The walls of the tubes are heated by the concentrated solar radiation, and solid particles in dense suspension circulate inside them (Fig. A.1). Boissière et al. [START_REF] Boissière | Experimental hydrodynamic study of gas-particle dense suspension upward flow for application as new heat transfer and storage fluid[END_REF] have studied the hydrodynamics of upward flow of a gas-solid dense suspension at ambient temperature, in a vertical two tube bundle of small diameter tubes (diameter 34 mm), which have their bases immersed in a slightly pressurized fluidized bed. The results obtained confirmed that it is possible to ensure stable upward flow of dense gas-solid suspensions in a bundle of tubes in parallel. Operating conditions for stable upward flows of the suspension and an equal repartition of the total solid flow rate between the tubes were determined experimentally. Boissière et al. [START_REF] Boissière | Experimental hydrodynamic study of gas-particle dense suspension upward flow for application as new heat transfer and storage fluid[END_REF] demonstrated that the suspension at the inlet is under minimum fluidizing conditions at any solid flow rate. Hence, the gas velocity in the tube is the sum of the aeration velocity and the minimum fluidization velocity.

The key to the proposed process is effective heat transfer from the solar heated surfaces to the heat transfer fluid, i.e. the circulating solid suspension. In order to understand this fully and so optimise the design of the solar receiver, it is of paramount importance to know how particles behave inside the bundle of small tubes within it. Indeed, the radial movement induces particle renewal at the wall, and determines the amount of heat removed from the wall to the bulk of the bed. To access to the particle motion in the solar receiver, two different techniques (experimental and numerical) are employed in this study. Measurement of particle motion in opaque suspensions such as this is difficult since conventional optical methods cannot be applied. Positron emission particle tracking (PEPT) is used here since it is the only non-intrusive method that is capable of giving detailed information on the particle circulation within the tubes.

The two CFD approaches commonly used for exploring fluidized bed systems are Eulerian-Lagrangian and Eulerian-Eulerian models. In the Eulerian-Lagrangian approach, particles are modelled as discrete elements and the Newtonian equations of motion for each individual particle are solved, with inclusion of the effects of particle collisions and forces acting on the particles due to the gas [START_REF] Xie | Eulerian-Lagrangian method for three-dimensional simulation of fluidized bed coal gasification[END_REF]. The particle-particle and particle-wall interactions may be accounted for in a deterministic manner rather than stochastic. This method is computationally expensive when the number of real particles to be tracked is high (more than 10 10 particles in our geometry). The silicon carbide particles used are less than 100 lm in diameter, so considering the solar receiver dimensions more than several billion of particles have to be considered. As the computation demands are considerable and limiting, the Eulerian-Lagrangian approach is limited to simulation of lab-scale fluidized beds with large particles. In the Eulerian-Eulerian model, the two phases are treated as interpenetrating continua and separated but coupled mass, momentum and energy Eulerian transport equations are written for the fluid and particle phases. The kinetic theory of granular flow and frictional theory are used to describe the rheology of the particle phase. This model is the most frequently used approach for predicting the dynamic behavior of fluid particle systems on a large scale up to and even including commercial-scale reactors [START_REF] Fede | 3D unsteady numerical simulation of the hydrodynamic of gas phase polymerization pilot reactor[END_REF].

The aim of this study is to examine the capability of the mathematical models implemented in NEPTUNE_CFD code to predict the behavior of an upflow bubbling fluidized bed by comparison with very accurate experimental measurements and to enhance the understanding of the physical mechanisms involved. The goal is to achieve the couple the hydrodynamics and the heat transfer in order to optimise the design of the solar receiver.

Experimental description: PEPT technique

Positron emission particle tracking (PEPT) is a technique for tracking a single radioactive tracer particle. It allows non-invasive observation of the motion of a single particle within a dense bed of similar particles. PEPT uses positron emitting radioisotopes which have the unique attribute that their decay leads to simultaneous emission of a pair of back-to-back c-rays. From the detection of a small number of such pairs the tracer position can be determined by triangulation. As a result accurate tracking is possible even in very dense systems involving significant c-ray attenuation and scattering.

Using the positron camera at the Positron Imaging Centre (University of Birmingham) a SiC tracer particle moving in the upflow bubbling fluidized bed can be located by PEPT to within 0.5 mm (in 3 dimensions). The PEPT technique employs an iterative algorithm [START_REF] Parker | Positron emission particle tracking-a technique for studying flow within engineering equipment[END_REF] to discard outlying data arising from scattered radiation. Once the parameters of the algorithm are selected, the tracking of a tracer particle for an extended period of time in a closed, circulating system builds up an integrated picture of particle behavior at each point in space and this allows the visualisation of the particle behavior and reconstruction of maps of axial and radial movement, circulation times and pseudo-density (using the PEPT "occupancy" function). The resulting locations are of the form (t k , x k , y k , z k ).

The tracers can be produced either by direct irradiation of the sample in a cyclotron, converting oxygen in the sample directly to 18 F, or by irradiation of water, which is then exchanged with, or attached to, molecules on the surface of the tracer. SiC hardly adsorbs 18 F; therefore in a separate project a technique was developed for depositing a very thin layer of Al 2 O 3 on SiC tracer particles under ambient conditions using gas-phase deposition in a fluidized bed [START_REF] Valdesueiro | Enhancing the activation of silicon carbide tracer particles for PEPT applications using gas-phase deposition of alumina at room temperature and atmospheric pressure[END_REF]. The core-shell particles produced were identical in properties that affect the hydrodynamic behavior in the bed (density, shape and size). The density comparison between uncoated SiC and the sample coated with 40 cycles is detailed in [START_REF] Valdesueiro | Enhancing the activation of silicon carbide tracer particles for PEPT applications using gas-phase deposition of alumina at room temperature and atmospheric pressure[END_REF]. A density difference of 0.8% is measured, which is a negligible difference. A single tracer particle was then mixed with other SiC particles and added to the Dispenser Fluidized Bed (DiFB). The tracer is supposed to be representative of all the particles of the bed, these are treated as statistically indistinguishable and therefore the trajectories are considered equivalent. In order to obtain sufficient experimental data it was necessary to operate this rig in a continuous circulation mode so that sufficient passes of the radioactive tracer could be observed. A Geiger counter was used to detect the particle as it entered the uplift transport (length 1.1 m) to check that the tracer particle was circulating in the experimental loop. A sketch of the single 30 mm ID tube in the circulating experimental rig is shown in of the set-up during the PEPT experiments, with the two positron "cameras" on each side. The dense upward-flowing suspension rises up the transport tube, which has additional aeration, and terminates in a cyclonelike disengaging zone. From here the disengaged solids fall under gravity into a downcomer. It is necessary to feed particles from the downcomer into the pressurized DiFB in a controlled way using a purpose-designed venturi nozzle [START_REF] García-Triñanes | Hydrodynamics and particle motion in upward flowing dense particle suspensions: application in solar receivers[END_REF]. In the experiments, the PEPT camera was able to locate and track particles over a height of 500 mm, from the top of the DiFB, including the region situated just before the aeration port. In this way it was possible to obtain a distribution of all locations in the transport tube over a typical 2-h experimental run (half life of 18 F = 110min).

For purposes of comparison with the simulation results a region situated 10 cm above the aeration port will be considered here. The operating parameters including aeration flow rate of the uplift transport tube, air flow in the venturi and fluidisation flow rate were selected to ensure steady flow of solid in the tubes, so that the results have high statistical significance, as reported in previous experiments [START_REF] García-Triñanes | Hydrodynamics and particle motion in upward flowing dense particle suspensions: application in solar receivers[END_REF]. The flow in the vertical tubes is greatly influenced by the fluidisation flow rate of the DiFB and the amount of aeration introduced directly to the tube. Table A.1 summarizes the operating conditions for the experiments carried out in this study, which were kept constant.

The position data were generated in a Cartesian coordinate system, which is convenient for viewing projections on the coordinate planes. Fig. A.3 shows an example of a tracer particle moving upwards with time, where y is the vertical coordinate. As in related work on bubbling beds by Stein et al. [START_REF] Stein | Solids motion in bubbling gas fluidised beds[END_REF], vertical motion consists of a series of vertical "jumps" which are associated with bubble motion, interspersed with quiescent periods in which the particle may descend in the downward-flowing dense phase between bubbles. The information obtained from PEPT yields a continuous trajectory of a single labelled particle while present in the field of view.

During this PEPT experiment, the tracking precision of a 68 lm tracer reaches 0.85 mm in 3D with a location frequency of 10 Hz for a tracer moving at 0.2 m/s. Acrylic perspex was chosen specifically for the construction of the apparatus for this experiment because of its low gamma ray attenuation. The length of the tube is 500 mm and the precision was found to be at least 0.6 mm [START_REF] García-Triñanes | Hydrodynamics and particle motion in upward flowing dense particle suspensions: application in solar receivers[END_REF]. (Tracking can still be performed in apparatus of higher attenuation, but precision is then reduced.).

Determination of solid volume fraction

Let us define a cylindrical grid of annular cells, the occupancy of each cell ' Oc ′ being defined as the fraction of the overall experimental duration t T which the tracer spends in each cell and is given by:

Oc = npass n=1 t cell,n t T (1)
where t cell,n is the single-pass residence time of the particle for pass n in a given cell. n pass is the total number of passes through the cell during the experiment. If we assume that all the particles are statistically indistinguishable (ergodicity) then Oc can be written as the ratio of N cell the mean number of particles in the cell to N T the total number of particle in the tube:

Oc = N cell N T (2)
The time outside the field of view is also tracked and it represents on average between 3 and 4.5 times the time spent inside the transport tube depending on the experimental conditions.

The mean solid volume fraction of particle is calculated as follows:

a p = N cell • m p q p • v cell = N T • Oc • m p q p • v cell = Oc • M T q p • v cell (3)
where m p is the mass of one particle, M T is the overall mass of the powder contained in the system and v cell represents the volume of the cell considered.

The calculation of the time-average value for each corresponding element in the circumference gives the plot presented in Fig. A.8.

Determination of time-averaged Eulerian solid velocity

To diminish the effect of the individual location error on the velocity determination, the local instantaneous particle velocities are calculated using an interpolation from six subsequent particle locations. The particle velocity in direction i at the kth particle track location, u p,i (t k ), is calculated from successive values of location, P k,i the ith component of the point in space at time t k :

u p,i (t k ) =0.1 P k+5,i -P k,i t k+5 -t k + 0.15 P k+4,i -P k-1,i t k+4 -t k-1 + 0.25 P k+3,i -P k-2,i t k+3 -t k-2 + 0.25 P k+2,i -P k-3,i t k+2 -t k-3 + 0.15 P k+1,i -P k-4,i t k+1 -t k-4 + 0.1 P k,i -P k-5,i t k -t k-5 (4)
The 6-point method used for PEPT introduces some smoothing (particularly in chaotic/turbulent systems) but reduces the overall error in velocity calculations. The instantaneous tracer particle velocity can be calculated according to Eq. ( 4). For overlapping sets, the tracer particle velocity is obtained by comparing 11 consecutive locations [START_REF] Stewart | Simulated and measured flow of granules in a bladed mixer-a detailed comparison[END_REF] and performing a weighted rolling average by replacing each data point with the average of the neighbouring six estimates of the resulting velocities.

During each loop, the tracer particle left the detectable range, because the PEPT detectors only covered approximately 500 mm of the transport tube. Therefore, the readings immediately before exit and after return of the particle were discarded for the purposes of velocity determination. The velocity is subsequently assigned to the cell considering the interpolation over 11 consecutive locations of the tracer particle which reduces the effect of PEPT measurement error on the results.

Finally, the time-averaged particle velocity is given by:

U p,i = 1 n meas nmeas i=1 u p,i (t k ) (5)
where n meas is the number of measurements in the cell. The time-variance of solid phase velocity is defined as follows:

s 2 p,i = 1 n meas nmeas n=1 (u ′ p,i ) 2 (6)
where u p,i ′ = u p,i (t k ) -U p,i is the solid fluctuating velocity. The locations in PEPT occur at random time intervals due to the random nature of radioactive decay, and due to the inefficiencies in detection. The average value of velocity error for PEPT is reasonable due to the six-point-averaging method, however the standard deviation of the error for this particular set of data is 0.2 m/s. The location precision has been dominated by errors in the z direction (0.6 mm) because of the low number of detection elements along the axis [START_REF] García-Triñanes | Hydrodynamics and particle motion in upward flowing dense particle suspensions: application in solar receivers[END_REF]. This filter will smooth the data by removing high-frequency events however it may be too restrictive in the case of the radial component.

Numerical simulation description: Euler n-fluid approach

Three-dimensional numerical simulations are carried out using the code NEPTUNE_CFD. This Eulerian n-fluid unstructured parallelized multiphase flow software has been developed in the framework of the NEPTUNE project financially supported by CEA (Commissariat à l'Energie Atomique), EDF (Electricite de France), IRSN (Institut de Radioprotection et de Sureté Nucléaire), and AREVA-NP [START_REF] Méchitoua | An unstructured finite volume solver for two-phase water/vapour flows modelling based on elliptic oriented fractional step method[END_REF]. The modelling approach for poly-dispersed fluidparticle flows is implemented by the Institut de Mécanique des Fluides de Toulouse (IMFT) in the NEPTUNE_CFD V1.08 version. The numerical solver has been developed for High Performance Computing [START_REF] Neau | High Performance Computing (HPC) for the Fluidization of Particle-Laden Reactive Flows[END_REF][START_REF] Neau | Neptune_CFD High Parallel Computing Performances for Particle-Laden Reactive Flows[END_REF].

Mathematical models

The Eulerian n-fluid approach used is a hybrid method [START_REF] Morioka | Modeling of gas and solid particles two-phase flow and application to fluidized bed[END_REF] in which the transport equations are derived by phase ensemble averaging for the continuous phase and by use of the kinetic theory of granular flows supplemented by fluid effects for the dispersed phase. The momentum transfer between gas and particle phases is modelled using the drag law of Wen and Yu [START_REF] Wen | Mechanics of fluidization[END_REF], limited by the Ergun [START_REF] Ergun | Fluid flow through packed columns[END_REF] equation for dense flows [START_REF] Fede | 3D numerical simulation of a lab-scale pressurized dense fluidized bed focussing on the effect of the particle-particle restitution coefficient and particle-wall boundary conditions[END_REF][START_REF] Gobin | Fluid dynamic numerical simulation of a gas phase polymerization reactor[END_REF]. The collisional particle stress tensor is derived in the frame of the kinetic theory of granular media [START_REF] Boëlle | Second-order prediction of the prediction of the particle-phase stress tensor of inelastic spheres in simple shear dense suspensions, In Gas-Particle flows[END_REF]. In the present study the gas flow equations are treated as laminar because the gas Reynolds stress tensor in the momentum equation is neglected compared to the drag term. For the solid phase, a transport equation for the particle random kinetic energy, q 2 p , is solved. The gas-particle turbulent correlation is negligible. The effects of the particle-particle contact force in the very dense zone of the flow are taken into account in the particle stress tensor by the additional frictional stress tensor [START_REF] Srivastava | Analysis of a frictional kinetic model for gas/particle flows[END_REF].

Numerical parameters 3.2.1. Geometry

The geometry is that used by Boissière et al. [START_REF] Boissière | Experimental hydrodynamic study of gas-particle dense suspension upward flow for application as new heat transfer and storage fluid[END_REF] simplified to access to the gas pressure drop in a longer tube. Only one tube of the solar receiver was simulated and the volume of the pressurized fluidized bed was divided by two. The tube diameter for the simulation is 34 mm instead of 30 mm for the PEPT experiments. As the two tube diameters are very close, the authors believe that this difference does not have a significant effect on the hydrodynamics.

Mesh

The mesh (Fig (20 cm X 20 cm), a height of 40 cm and is equipped with a lateral solid entrance. A tube of 2.05 m height and 34 mm in internal diameter has its lower end submerged in the DiFB to a distance of 10 cm from the gas distributor. Secondary air injection, hereafter termed "aeration", is set at 57 cm from the bottom of the tube.

Phase properties

The properties of the SiC particles used for the experimental measurements are presented in Table A.2. The experiments were carried out at ambient temperature and the fluidization gas was air. the non-sphericity of particles in the drag term Loth [START_REF] Loth | Drag of non-spherical solid particles of regular and non regular shape[END_REF] proposed a correction to the Wen and Yu drag law using the coefficients (C shape and f shape ):

Re * p = C shape Re p f shape , (7)
C * d = C d C shape . ( 8
)
where Re p = a g qg | vr | p dp lg is the particle Reynolds number. The spheroidal shape factor f shape is inversely proportional to the change in terminal velocity since the drag dependence is linear in creeping flow (Re p ). It ranges from 1 for a sphere to 1.19 for a tetrahedron. C shape can be thought of as a Newton-drag correction. It ranges from 1 for a sphere to 4.5 for a tetrahedron. In the creeping flow regime Re p ≪ 1, as encountered in this study, the parameter to take into account is f shape . According to the Wen and Yu drag law, Fig. A.6 shows that a value of f shape of 3.2 over range of Loth value recommended.

Hence for the numerical simulations, the classical drag laws without shape correction have been retained by assuming a monodispersed mean diameter of 40 lm for the silicon carbide particles. This mean diameter predicts a void fraction for homogeneous bed expansion in good agreement with the experimental measurement.

The gas density is a function of gas pressure and obtained by following the ideal gas law with a constant temperature of 293K.

The phase properties used for the simulation are summarized in Table A.3.

Boundary conditions

The geometry is composed of three inlets:

• The fluidization gas distributor, which is a wall for the particles and where an air flow rate of 2.268 kg/h is imposed. This air flow rate corresponds to a superficial gas velocity of 1.66cm/s ≈ 3 × U mf , where U mf is the minimum fluidization velocity. • The lateral injection of solid in the DiFB with a solid flow rate imposed of 25.5kg/s/m 2 (relative to the tube surface) with a kg/s which corresponds to a gas velocity of 4 cm/s in the tube as for the experimental rig. • The aeration is a gas inlet set at 56.4 cm from the lower end of the tube where a gas mass flow rate of 5.3875 • 10 -5 kg/s is imposed.

The superficial gas velocity in the tube is equal to 4.5 cm/s= 9U mf (sum of both aeration and gas coming from DiFB as discussed above).

The geometry is also composed of two outlets with an imposed pressure equal to atmospheric pressure: the top of the tube is a free outlet for air and particles, and an outlet placed at the top of the DiFB imitates the regulation valve. The pressure regulation valve is modeled as a pressure loss implemented on the cells of the valve conduit introduced so as to keep a constant pressure inside the DiFB equal to the atmospheric pressure plus 240 mbar.

The wall-type boundary condition is no-slip for the two phases.

Simulation progress

The numerical simulations have been performed on parallel computers with 140 cores. A numerical simulation is divided into two steps: a transitory step to reach a predicted constant total mass of particles in the column (solid mass flow rate injected equal to solid mass flow rate leaving the geometry) and an established regime during which the statistics are computed for 400 s. The time-average of the solid flow variable w p is defined by:

w p (X) = n w p (X, t n )a p (X, t n )Dt n n a p (X, t n )Dt n (9)
and the time-variance by:

w ′ p (X) 2 = w p (X, t n ) -w p (X) 2 (10)
where a p and q p are respectively the volume fraction and the density of the particle, and Dt the time step.

Results and discussion

The numerical results for gas pressure drop along the tube are in very good agreement with the experimental measurements obtained by Boissière et al. [START_REF] Boissière | Experimental hydrodynamic study of gas-particle dense suspension upward flow for application as new heat transfer and storage fluid[END_REF] (Table A.4). As can be noticed, the numerical simulations predict perfectly the decrease of the gas pressure due to the aeration. Below the aeration, the mean solids volume fraction is equal to 43% and this decreases to 36% above the aeration. This decrease of the solids volume fraction can be observed in solid is injected into the DiFB from the right side of the geometry. The height of solid in the DiFB is controlled by the pressure imposed at the surface of the bed.

Experimental vs. simulation results

The PEPT measurements and the numerical predictions are compared as radial profiles of flow variables extracted at 10 cm above the aeration level. This position, imposed by the experimental mockup, is close to the aeration and we can presume that the flow is not fully stabilized. As the simulated tube is about two meters in height, the influence of the distance from the aeration is also discussed. The numerical predictions are plotted at 10 cm and 50 cm above the aeration. As the diameters of tubes used in the experiment and in the simulation are not exactly the same (30 mm for the experiment and 34 mm for the numerical simulation) the radial position is normalized by the tube radius. Moreover, owing to the axisymmetry of the tube receiver, the time-averaged flow may be assumed to obey cylindrical symmetry. So, spatial averaging of the time-averaged variables was performed in the azimuthal direction.

Fig. A.8 shows the radial profile of the solids volume fraction. PEPT measurements and numerical predictions both present a minimum at the center of the tube and a maximum near the wall. These results are in relatively good agreement with a slight over-prediction from the model. When the distance from the aeration increases, the solid volume fraction decreases. Indeed, the size and velocity of the bubbles increase as a function of the height in the tube which increases the bed expansion. Fig. A.9 shows the vertical component of the solids phase velocity and the associated solids mass flux defined as the product of this vertical velocity, the solids volume fraction and the solids density. The numerical results and the experimental data for the vertical component of solid flux are similar, with a maximum value at the center and a negative minimum close to the wall; this is due to bubble rise in the center of the tube. A rotating vortex can be observed near the wall. This trend is the same as that normally observed in dense bubbling fluidized beds [START_REF] Fotovat | Sand-assited fluidization of large cylindrical and spherical biomass particles: experiments and simulation[END_REF]. In the present study, the time-averaged solids flow direction reverses at approximatively 78% of the tube radius. This downward flow has a strong influence on the behavior of the heated dense particle suspension inside the solar receiver. Indeed, the particle suspension will be heated prior to reaching the irradiated area by the hot flow coming from the upper zone. This is confirmed by the temperature measurements conducted by the authors' partners in an experimental solar receiver [START_REF] Benoit | On-sun demonstration of a 750 • C heat transfer fluid for concentrating solar systems: dense particle suspension in tube[END_REF]. It can be noticed, that the downward solid mass flux is slightly overestimated by the numerical simulations. The extrema of the solid mass flux, as does the velocity, increase as a function of the height in the tube due to the increase of the bubbles velocity. ponents of solids velocity ( U 2 p,x + U 2 p,z , x and z are the directions perpendicular to the gravity direction), which corresponds to the time-averaged radial velocity by assuming that the azimuthal component of solids velocity is equal to zero by axisymmetry of the geometry, and the associated solids mass flux. The numerical predictions and experimental results are in good agreement in the central part of the tube. However, a discrepancy can be observed near the wall. Indeed, the experimental measurements do not tend to zero at the wall probably due to the error location reported in the horizontal z direction as mentioned before. The simulations seem more realistic as far as the horizontal radial velocity is concerned. It should also be noted that the comparison is made very close to the aeration point, where the flow is not fully established. Thus, the norm of horizontal solids velocity is not equal to zero. It is also noteworthy that the numerical value of this velocity decreases very fast a few tens of centimeters above the aeration point and tends to zero. The horizontal velocity is three times lower than the vertical one in the center of the tube; the flow is not isotropic and predominantly vertical. The horizontal solids mass flux decreases, as does the velocity, from the center to equal zero at the wall. However, the simulation results seem to decrease faster than the experimental data near the wall. This radial movement of the particle is very important for heat exchange between the wall and the particle suspension, in order to avoid "hot spot" developing at the wall.

Fig.

A.13 presents radial profiles of vertical and horizontal components of time variance of solids phase velocity. In the numerical approach the solids velocity time variance is the sum of a predicted

part (< u p,x ′ u p,x ′ >, < u p,y ′ u p,y ′ >, < u p,z ′ u p,z ′ >) which corresponds
to the collective motion of particles and a modelled uncorrelated part (q 2 p ). This last part is less than 10 -5 m 2 /s 2 and could be negligible compared to the fluctuations of the predicted solid velocity. The PEPT data and numerical simulations are in good agreement for the vertical variance. Thus, the numerical simulations are able to predict the vertical mixing of the particles well. At the center of the tube the standard deviation of the vertical component of solids velocity is of the same order (≈ 0.2 m/s) as the time-averaged vertical solids velocity. This means that the instantaneous particle velocity can be negative at the center leading to non-zero downward PEPT measurements show a variance of the horizontal velocity much lower than in the vertical direction. This anisotropy of the particle velocity fluctuations is not represented by the numerical results which overestimate the radial mixing. In contrast, the numerical simulations show that the particle velocity fluctuations are almost isotropic at the center of the tube and tend to be anisotropic near the wall.

Discussing the plausible reasons of discrepancy

As discussed above, the mathematical models that underpin the numerical simulation approach of this study are able to predict the behavior of an upflow bubbling fluidized bed in the solar receiver. Nonetheless, some discrepancies have been observed such as an overestimation of the downward solids mass flux by the numerical simulation.

To improve the simulation approach several tracks can be followed. First of all, it is important to consider the differences between the experimental and numerical geometries. The tube diameter for the numerical simulations was 34 mm in the simulation instead of 30 mm for the experimental solar receiver and its length was 2.06 m instead of 1.1 m. A mesh refinement was tested on a simplified fluidization column with the same diameter, the same particle properties and showed no effect on the bed expansion. However, a refinement of the mesh, especially in the radial direction, could improve the prediction of the downward flow near the wall.

Clearly, the treatment of the particles in the simulation is grossly simplified, making use of a single diameter of an (effectively) spherical particle to describe a size distribution of extremely nonspherical particles. As discussed above, the particle diameter was arbitrarily decreased in order to match the bed expansion measured by experiments but this method can be only a stopgap. One drawback of this diameter reduction is that the particle-particle collisions are not well described. Many studies were recently carried out in the literature [START_REF] Clift | Bubbles, Drops, and Particles[END_REF][START_REF] Hölzer | New simple correlation formula for the drag coefficient of non-spherical particles[END_REF][START_REF] Loth | Drag of non-spherical solid particles of regular and non regular shape[END_REF][START_REF] Tran-Cong | Drag coefficients of irregularly shaped particles[END_REF] to take into account the non-sphericity of particles in the formulation of drag coefficients. A 3D characterization by experimental measurement of irregularly shaped particles used in the study (Fig. A.5) is tricky but the way to model it still remains a challenge.

To reduce the over-prediction of downward solids mass flux near the wall, one possibility would be to increase the friction of the particles at the wall through adjustment of the wall-particle boundary condition, but the wall-particle boundary condition friction is already at the maximum through the no-slip condition used. Another possibility would be to increase the particle-particle friction. The model used in this study [START_REF] Srivastava | Analysis of a frictional kinetic model for gas/particle flows[END_REF] contains several parameters but the authors have retained the use of recommended value. Nonetheless, other models to describe frictional stress could be used. For example, Farzaneh et al. [START_REF] Farzaneh | The crucial role of frictional stress models for simulation of bubbling fluidized beds[END_REF] have recently shown that the results of a visco-plastic model proposed by Jop et al. [START_REF] Jop | A constitutive law for dense granular flows[END_REF] are in better agreement with the experiments for prediction of the bed hydrodynamics and the movement of the fuel particles than the models of Sundaresan [START_REF] Srivastava | Analysis of a frictional kinetic model for gas/particle flows[END_REF] and Schaeffer [START_REF] Schaeffer | Instability in the evolution equations describing incompressible granular flow[END_REF] conventionally used for CFD simulation.

Conclusions

The positron emission particle tracking technique and an Eulerian n-fluid approach were used in this work to shed light on dense particle suspension motion in an upflow bubbling fluidized bed applied as a heat transfer fluid in solar receiver. The objectives were not only to measure the properties of the system but also the test the utility of the computational model, which might then be used in future design work for this technology. Comparisons of the time averaged radial profiles of solid volume fraction, vertical and radial velocity and velocity variance of particles present satisfactory agreement and demonstrate the capability of NEPTUNE_CFD software to predict this peculiar upflow bubbling fluidized bed.

Although a strong assumption was made on drag term, the simulation was successful in predicting the axial gas pressure drop and the effect of a secondary air injection on this profile. Simulation also succeeded in reproducing the radial evolution of the vertical and horizontal solid velocity and the time variance of the solids vertical velocity which controls the axial mixing. Both numerical predictions and PEPT measurements describe an upward flow at the centre of the transport tube due to rising bubbles with a back-mixing flow near the wall which will strongly influence the solar to particles heat transfer mechanism. When compared to the experimental data, simulation predicts a slightly higher solid downward flow near the wall. The observed discrepancies between the experimental and simulation results were mainly attributed to the simplified treatment in the model of the nature of the particles.

Future numerical simulations will investigate the coupling between the hydrodynamics of the dense particle suspension and the heat transfer to reproduce the experiments conducted in the solar furnace. Interphase momentum transfer (kg

Nomenclature

• m -2 • s -2 ) K p granular diffusivity (m 2 • s -1 ) K col p collisional granular diffusivity (m 2 • s -1 ) K kin p kinetic granular diffusivity (m 2 • s -1 ) m p particle mass (kg • m -3 ) n p particle number density(m -3 ) P g mean gas pressure (N • m -2 ) q 2 p mean particle random kinetic energy (m 2 • • • -2 ) Re p particle Reynolds number (-) U k,i mean velocity of phase k (m • s -1 ) U mf minimum fluidization velocity (m • s -1 ) u k,i ′ fluctuating velocity of phase k (m • s -1 ) V d drift velocity (m • s -1 ) V f
superficial gas velocity (m • s -1 ) V r relative gas-particle velocity (m

(kg • m -1 • s -2 ) S col k,ij collisional stress tensor of the phase k (kg • m -1 • s -2 ) S fric k,ij frictional stress tensor of the phase k (kg • m -1 • s -2 ) S kin k,ij kinetic stress tensor of the phase k (kg • m -1 • s -2 ) t c
collisional timescale (s) t F gp mean gas-particle relaxation timescale (s) Subscripts g gas p particle

In the following equations, subscript k = g refers to the gas phase and k = p or q to the particle phase. a p q p in the particle transport equation represent n p m p where n p is the number density of particle centers and m p is the mass of a single particle: a p = npmp qp is an approximation of the local volume fraction of particle p. Hence, gas and particle volume fractions, a g and a p have to satisfy: a g + a p = 1 (A.1)

Mass transport equation:

∂ ∂t

a k q k + ∂ ∂x j a k q k U k,j = 0 (A.2)
where q k is the density of k-phase and U k,i is the i-component of its velocity.

Momentum transport equation

a k q k ∂U k,i ∂t + U k,j ∂U k,i ∂x j = -a k ∂P g ∂x i + a k q k g i + q=g,p I q→k,i - ∂S k,ij ∂x j (A.3)
where P g is the mean gas pressure, g i is the gravity i-component. I g→p are the interphase momentum transfer between particle and gas without the mean gas pressure contribution and S k,ij is the effective stress tensor of phase k.

Interphase transfer modelling:

According to the particle to gas density ratio, the dominant forces between the gas phase and particles are the drag and Archimedes' force, so the mean momentum gas-particle transfer term may be written: V rp,i =< v r,i > p = U p,i -U g,i -V dp,i (A.9)

I
V dp,i is the drift velocity which can appear due to turbulence [START_REF] Simonin | Eulerian prediction of the fluid/particle correlated motion in turbulent two-phase flows[END_REF] or sub grid effect [START_REF] Parmentier | A functional subgrid drift velocity model for filtered drag prediction in dense fluidized bed[END_REF] which is assumed to be negligible (due to the inertia of the particles).

Fig. A. 1 .

 1 Fig. A.1. Schematic view of a module of the solar receiver with dense upward solid flow.

  Fig. A.2 (left). On the right a photo depicts the disposition

Fig. A. 2 .

 2 Fig. A.2. Schematic drawing of the experimental rig (left) and location of the single tube circulating apparatus positioned relative to the detectors (right).

Fig. A. 3 .

 3 Fig. A.3. Example of individual trajectory against time (Fluidization velocity= 0.096m/s, Aeration velocity=0.04m/s).

Fig. A. 4 .

 4 Fig. A.4. 3D mesh for the numerical simulation 1,650,000 cells.

  . A.4) is composed of 1,650,000 hexahedra, based on the O-grid technique with approximately D r = 1.2 mm and D z = 1.5 mm. The DiFB has a horizontal sectional area of 0.04 m 2

Fig. A. 6 .

 6 Fig. A.5. SEM photograph of silicon carbide particles.

Fig. A. 7 .Fig. A. 8 .

 78 Fig. A.7. Instantaneous field of solid volume fraction. On the left: an overall view of the geometry, at the center: a cross section view at the center of the column and on the right: enlargement of the cross sectional view of the bottom part.

  Fig. A.7. TheTable A.4 Time-averaged gas pressure drop at the wall on 50 cm height. Exp. data [3] Numerical results Below the aeration 136 mbar/m 135 mbar/m Above the aeration 115 mbar/m 117 mbar/m Vertical solids mass flux.

Fig. A. 9 .

 9 Fig. A.9. Radial profiles of time-averaged vertical velocity and mass flux of solid phase at 10 cm above the aeration point.

Fig. A. 10

 10 shows the solids velocity vectors with the field of solid volume fraction in background. An upward flow in the center of the tube and a downward more concentrated solid flow near the wall are clearly observed.

Fig. A. 11

 11 presents the upward and absolute values of downward mean vertical solids mass flux. At 10 cm above the aeration, the downward solids mass flux is very low but not zero at the center of tube, at about 25 kg/s/m 2 , and very high at the wall, at about 180 kg/s/m 2 . This trend of downward vertical solids mass flux is the reverse of the upward solids mass flux with a maximum at the center of the tube and minimum near the wall. With the increase of the height in the tube, the downward solid mass flux decreases at the center of the tube while the upward solid mass flow increases. At the wall, the reverse phenomenon is observed. The ratio between the upward and the downward solid mass flux does not change with the increase of the distance from the aeration.

Fig. A. 10 .Fig. A. 11 .

 1011 Fig. A.10. Time-averaged solid vector velocity with selected pathline of particle, and solid volume fraction for a section of the riser tube between 45 cm and 50 cm above the aeration point predicted by numerical simulations.

Fig. A. 12

 12 Fig. A.12 shows the norm of the time-averaged horizontal com-

  above aeration Simulations 50 cm above aeration Experiments 10 cm above aeration (a) Norm of time-averaged horizontal components of solid velocity. above aeration Simulations 50 cm above aeration Experiments 10 cm above aeration (b) Horizontal solids mass flux.

Fig. A. 12 .

 12 Fig. A.12. Radial profiles of time-averaged horizontal solid velocity and mass flux of particle at 10 cm above the aeration point.

Fig. A. 13 .

 13 Fig. A.13. Radial profiles of time variance of solids phase velocity at 10 cm above the aeration point.

Table A .1

 A Experimental conditions. *Between y=150 mm and 250 mm above the aeration, where y is the vertical coordinate with respect to the bottom of the PEPT camera.

	U mf	U mb	Fluidization velocity	Aeration velocity	Fluidization flowrate	Aeration flowrate	Venturi nozzle flowrate	Net solids mass flux*
	(mm/s)	(mm/s)	(m/s)	(m/s)	(Nm 3 /h)	(NL/h)	(NL/h)	(kg/s/m	2 )
	5	8	0.096	0.04	4.0	190	48	27.2

  g→p,i = -I p→g,i = -a p q p

								gp t F	V rp,i	(A.4)
	1 t F gp	=	3 4	q g q p	| v r | p d p	C d (Re p )	(A.5)
	C d (Re p ) =	C d,WY min C d,WY , C d,Erg else a g < 0.7 if a g ≥ 0.7	(A.6)
				
	C d,WY =	 	24 Rep 1 + 0.15Re 0.687 p 0.44a -1.7 g	a -1.7 g	Re p < 1000 Re p ≥ 1000	(A.7)
	C d,Erg = 200	(1 -a g ) Re p	+	7 3	, Re p = a g	q g | v r | p d p l g	(A.8)

  Fr = 0.05, 0 = p/4, n = 2, m = 5, a p,min = 1 -a g,mf = 0.43 (A.23)

	Particle stress modelling: S p,ij = S kin p,ij + S col p,ij + S frict p,ij	(A.10)	K kin p =	2 3	q 2 p	5 9	t F gp (1 + a p g 0 0 c) 1 +	5 9	t F gp	n c p t c	-1	(A.26)
																				0 c =	3 5	(1 + e c)	2 (2e c -1), n c =	(1 + e c )(49 -33e c ) 100	(A.27)
	S kin p,ij + S col p,ij = P p -k p	∂U p,n ∂x n	d ij -l p	∂U p,i ∂x j	+	∂U p,j ∂x i	-	2 3	∂U p,n ∂x n	d ij (A.11)	K col p = a p g 0 (1 + e c )	  6 5	K kin p +	4 3	d p	2 3	q 2 p p	 	(A.28)
	P p = a p q p [1 + 2a p g 0 (1 + e c )]	2 3	q 2 p	(A.12)
	g 0 = 1 -	a p a p,max	-2.5ap,max	, e c = 0.9, a p,max = 0.64	(A.13)
	l p = a p q p m kin p + m col p			(A.14)
	m kin p =	1 2	t F gp	2 3	q 2 p (1 + a p g 0 V c) × 1 +	t F gp 2	s c p t c	-1	(A.15)
	V c =	2 5	(1 + e c )(3e c -1), s c =	1 5	(1 + e c )(3 -e c )	(A.16)
																			
	m col p =	4 5	a p g 0 (1 + e c )	 m kin p + d p	2 3	q 2 p p		(A.17)
	k p = a p q p	4 3	a p g 0 (1 + e c )d p	2 3	q 2 p p	(A.18)
	1 t c p	= 24	a p g 0 d p		2 3	q 2 p p					(A.19)
	S	frict p,ij = P frict p	∂U p,i ∂x j	+	∂U p,j ∂x i	-	2 3	∂U p,n ∂x n	d ij	(A.20)
	m	frict p	=	√ 2P f sin0 a s q p 2I 2D + 8/3	d 2 p p q 2	(A.21)
	P f = Fr	(a p -a p,min [ap,max -a p]	n m			(A.22)
	I 2D =		∂U p,i ∂x j		+	∂U p,j ∂x i		∂U p,i ∂x j	-	2 3	∂U p,i ∂x i	2	(A.24)
				Particle random kinetic energy transport equation:
	a p q p		∂q 2 p ∂t	+ U p,j	∂q 2 p ∂x j		=	∂ ∂x j	a p q p K kin p + K col p	∂q 2 p ∂x j	(A.25)
																			-S kin p,ij + S coll p,ij	∂U p,i ∂x j	-	a p q p gp t F	2q 2 p	(A.25)
																			-	1 2	1 -e 2 c	a p q p t c p	2 3	q 2 p -a p q p m p frict	4 3	q 2 p d 2 p
																				(A.25)

f d ij -a p q p m

Acknowledgements

This work was developed in the frame of the CSP2 European project. The authors acknowledge the European Commission for cofunding the "CSP2" Project -Concentrated Solar Power in Particles -(FP7, project N 282 932). This work was granted access to the HPC resources of CALMIP under the allocation P1132 and CINES under the allocation gct6938 made by GENCI.

Appendix A. NEPTUNE_CFD mathematical equations

The Eulerian n-fluid approach used is a hybrid method, in which the transport equations are derived by phase ensemble averaging for the continuous phase and by using kinetic theory of granular flows supplemented by fluid and turbulent effects for the dispersed phase, thanks to a joint fluid-particle probability density function (PDF) approach.


Jair Montoya-Martínez

Jean-François Cardoso

Alexandre Gramfort

Jean-Franc ¸ois Cardoso

Alexandre Gramfort Ltci

Jean-Franc ¸ois Cardoso

Caveats with stochastic gradient and maximum likelihood based ICA for EEG

Keywords: Independent component analysis (ICA), maximum likelihood, stochastic gradient method, infomax, electroencephalography (EEG), neuroscience

Introduction

Independent Component Analysis (ICA) is a multidimensional statistical method that seeks to uncover hidden latent variables in multivariate and potentially high-dimensional data. In the ICA model we consider here, the observations x satisfy x = As, where s are referred to as the sources or independent components, and A is the mixing matrix considered unknown [START_REF] Jutten | Blind separation of sources, part I: An adaptive algorithm based on neuromimetic architecture[END_REF]. In the following, we assume as many sources as sensors: A is a square matrix. This model is usually described as a latent linear stochastic model, where x and s are random variables (r.v.) in R N , and A ∈ R N ×N is a nonsingular matrix. The goal of ICA is, given a set of observations of the r.v. x, to estimate the hidden sources s and the unknown mixing matrix A. In order to accomplish this task, the key assumption in ICA is that the components s 1 , s 2 , . . . , s N are mutually statistically independent [START_REF] Cardoso | Blind signal separation: statistical principles[END_REF], a plausibe assumption if each individual source signal is thought to be generated by a process unrelated to any other source signal.

In neuroscience, and in particular when working with Electroencephalography (EEG) data, ICA is extremely popular. It is used for artifact removal as well as estimation of brain sources. Linear ICA is justified by the fact that EEG data are linear mixtures of volume-conducted neural activities [START_REF] Makeig | Independent component analysis of electroencephalographic data[END_REF]. Each brain source is thought to represent near-synchronous local field activity across a small cortical patch [START_REF] Nunez | Electric fields of the brain: the neurophysics of EEG[END_REF], which can be modeled as an electrical current dipole (ECD) located within the brain [START_REF] Scherg | Use of prior knowledge in brain electromagnetic source analysis[END_REF].

To solve the ICA estimation problem, we need to estimate a linear operator W ∈ R N ×N , such that ŝ = Wx = WAs ≈ s, where ŝ is an estimation of the sources. In the context of EEG, the estimated mixing matrix A = W -1 gives us information about how the estimated sources are seen at the sensor level. Indeed each column of A can be visualized as a scalp map (topography). This helps the EEG users to identify plausible brain sources which correspond to ECDs. For such sources, topographies are spatially smooth and exhibit a dipolar pattern.

A common approach to tackle the ICA problem is to cast it as a maximum likelihood estimation problem [START_REF] Pham | Blind separation of mixture of independent sources through a quasimaximum likelihood approach[END_REF]: given a probability density function (p.d.f.) p s (s), associated with the sources, and a set X = {x(1), x(2), . . . , x(T )} = {x j } j=T j=1 , containing independent and identically distributed (i.i.d) samples of the r.v. x, one wants to find the unmixing matrix W that maximizes the log-likelihood function (W, X ) (for the sake of simplicity, we make the dependence on X implicit in ):

(W) = T j=1 N i=1 log p si (w i x j ) + log|det W| (1)
where w i denotes the i-th row of W. One of the most popular algorithm in the EEG community is Infomax [START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF] and it can be shown to follow this likelihood approach [START_REF] Cardoso | Infomax and maximum likelihood for blind source separation[END_REF].

In order to maximize the log-likelihood function [START_REF] Amari | Natural gradient works efficiently in learning[END_REF], we have at our disposal two different families of optimization methods: batch methods, such as gradient descent, which use at each iteration the entire set of observations X , and stochastic methods which access at each iteration only one observation x j , or a small group of observations B = {x j } ⊂ X (also known as mini-batch). When using stochastic gradient (SG) methods, the gradients used as update directions, with one sample or a mini-batch, are affected by 'noise' [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF]. The consequence is that unless so-called step size annealing strategies are employed, SG will not reach a minimum of the minimized function [START_REF] Schaul | No more pesky learning rates[END_REF]. On the contrary, gradient descent (GD), which is a non-stochastic batch method, does guarantee a decay of the minimized function at every iteration and does reach points with zero gradients (cf. Prop. 1.2.1 in [START_REF] Bertsekas | Nonlinear Programming[END_REF]). Convergence rates can be up to linear for strongly convex functions with Lipschitz gradients. However, one update of the parameters by GD requires a full pass on the whole dataset while SG does already reduce the cost function after accessing a fraction of it. That is why when working with many samples, which is the case for EEG, SG exhibits a rapid convergence during the early stages of the optimization procedure, yet this convergence then slows down and the cost function reaches a plateau well before a point with zero gradient is reached. In other words, plain SG will stop too early if a high numerical precision solution is needed.

Infomax uses SG to maximize the log-likelihood function [START_REF] Amari | Natural gradient works efficiently in learning[END_REF]. In particular, it uses a mini-batch SG method in combination with a step size annealing policy, which is applied after one pass on the full data (Infomax considers one iteration as one pass on the full data). As in any stochastic method, the Infomax solver needs an initial step size.

In the first part of the paper, we explore the impact of algorithm initialization, the initial value of the step size, jointly with the annealing policy and the stopping criterion used by the standard Infomax implementation. We explain theoretically why the commonly used initialization of Infomax produces highly dipolar sources. We then explain the observation that Infomax can eventually waste a lot of computation time without converging, or worse can report convergence while the norm of the gradient is still high. Finally, using public EEG data and an alternative optimization strategy we investigate the impact of convergence on source dipolarity, highlighting specificities of EEG.

Infomax: description of the optimization algorithm

The maximum likelihood problem tackled by Infomax can be written as the following minimization problem:

W = argmin W L(W) (2)
where L(W) = -(W)/T denotes the normalized negative log-likelihood function.

In order to solve the problem (2), Infomax uses the relative gradient [START_REF] Amari | Natural gradient works efficiently in learning[END_REF][START_REF] Cardoso | Equivariant adaptive source separation[END_REF]]

L (W) = 1 T T j=1 φ φ φ(y j )y j -I N (3)
where y j = Wx j , where φ φ φ(y j ) = [φ 1 (y j1 ), . . . , φ N (y j N )] , and where

φ k (y j k ) = -p s k (y j k )/p s k (y j k ) = tanh (y j k /2)
In order to solve (2), the reference Infomax implementation, included for example in the EEGLAB software [START_REF] Delorme | EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis[END_REF], uses a mini-batch stochastic gradient method, whose iterative expression can be written as follows:

W k+1 = W k -α j∈B k φ φ φ(y j )y j -I N W k (4)
where, before each pass on the full data, the set of samples {x j } j=T j=1 is randomly permuted, and then, during the full pass, each mini-batch B k is created by taking, sequentially, a subset of samples {x j } of size |B k |. Once the pass on the full data is completed, the stopping criterion is checked and the annealing policy is applied to determine whether or not the step size α should be decreased. In this policy, the step size is never increased. Finally, this process is repeated until the stopping criterion is fulfilled or until the maximum number of iterations is reached.

Let us denote ∆ k = W k+1 -W k . The stopping criterion used by the standard Infomax implementation is ∆ k 2 F < tol, where • F is the Frobenius norm, tol is by default 10 -6 if N ≤ 32, and 10 -7 otherwise. This implementation uses the following heuristic for its annealing policy: if the angle between matrices ∆ k and ∆ k-1 is larger than 60 • , i.e., arccos (Trace(∆ k ∆ k-1 )/( ∆ k F ∆ k-1 F )) > π/3, it decreases the step size by 10% (α ← 0.9α), otherwise the step size remains the same.

Regarding the stopping criterion ∆ k 2 F < tol, it is important to notice that by Eq. ( 4), ∆ k is proportional to the step size α so that the algorithm will stop if the gradient or the step size is small. Even if the stoppping criterion is not met, the step size may have become small enough (even using the standard default values) to prevent any significant update. Example of such behaviors on EEG data are given in Section 4.

3 Assessing the performance of ICA using dipolarity Dipolarity metric. ICA can be seen an unsupervised learning method, therefore, it is in general difficult to assess its performance in real life scenarios where the ground truth is unknown. In order to help to mitigate this issue, when using EEG data, Delorme et al. [START_REF] Delorme | Independent EEG sources are dipolar[END_REF] proposed to use the physics underlying the propagation of the electromagnetic field throughout the head. Physics states that the signals measured on the scalp can be modeled as linear mixtures of the electrical activities generated by ECDs located inside the brain. To assess the biological plausibility of ICA sources, Delorme et al. [START_REF] Delorme | Independent EEG sources are dipolar[END_REF] proposed to take each column of the estimated mixing matrix A, which can be represented as a topography, and compute how well it can be modeled by a single ECD. The following metric is defined by:

dipolarity( A j ) = (1 -A j -Āj 2 2 / A j 2 2 ) × 100 (5)
where A j , Āj denote respectively the j-th column of the estimated mixing matrix and the corresponding topography obtained by fitting a single dipole. Taking into account ( 5), the "Near-Dipolar percentage (ND%)" of an ICA decomposition is defined in [START_REF] Delorme | Independent EEG sources are dipolar[END_REF] as: ND%( A) = {#j : dipolarity( A j ) > τ }/N , that is, the percentage of returned components whose topographies can be modeled by a single ECD with more than a specified dipolarity threshold τ (specified as percentage of explained/residual variance). Following [START_REF] Delorme | Independent EEG sources are dipolar[END_REF], we will consider an ICA source to be biologically plausible when its dipolarity is larger than τ = 90.

Relationship between initialization and dipolarity. ICA leads to nonconvex optimization problems: solutions found by algorithms necessarily depend on their initialization. In this section we discuss connections between initialization and dipolarity.

Learning of a separating matrix W starts with some initial value W 0 . While it is generally possible to start with the identity matrix W 0 = I N , it is a sound and common practice to start with some whitening matrix, that is, with a matrix W 0 such that W 0 Σ x W 0 = I N where Σ x = Cov(x) denotes the covariance matrix of x. There are infinitely many such matrices; two popular choices are 'PCA' and 'sphering', which can be defined in terms of the eigen-value decomposition Σ x = UDU :

W pca = D -1/2 U , W sph = UD -1/2 U .
The topographies (the columns of matrix A 0 = W -1 0 ) associated with the three aforementioned initializations W 0 = I N , W sph , W pca , are displayed on Fig. 1. Of course, the topographies associated with W 0 = I N (Fig. 1(a)) are 'quasi-dipolar' in the sense that activating only one channel could be interpreted as the effect of a single source located just beneath the scalp. Much more striking is the fact that the sphering W 0 = W sph produces topographies which all look dipolar. Fig. 1(b)) shows 6 of them, randomly selected. Nothing similar is observed in Fig. 1(c) after PCA W 0 = W pca . See also Fig. 4, which shows the dipolarity index for all components after sphering or PCA, sorted in decreasing order.

If the dipolarity criterion is to be used for assessing the biological plausibility of a source, one has to understand why a simple sphering would produce dipolar topographies. An explanation can be provided by the observation that somehow sphering is the

Relationship between initialization and dipolarity

ICA leads to nonconvex optimization problems, therefore solutions necessarily depend on the initialization used. To illustrate this issue, let us analyze the following scenario. Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal to the identity matrix I N . The topograhies associated to the columns of A 0 are shown in Fig. 1(a). Each column of A 0 is putting a value of one at each sensor, and zero in the other ones. 'smallest' whitening transform. Indeed, if a whitening matrix is close to the identity, then it should not modify much the 'quasi-dipolar' patterns of Fig. 1(a) and this is what seems to happen upon observation of Fig. 1(b). So, in which sense would sphering be the 'smallest whitening transform'? An answer is provided by Theorem 1 of Eldar et al. [START_REF] Eldar | MMSE whitening and subspace whitening. Information Theory[END_REF], which implies that, among all whitening matrices W, the sphering matrix is the one with the minimal mean-squared difference E x -Wx 2 . In other words, among all white random vectors Wx, W sph x is the closest to x with closeness measured in the mean-squared sense. In other words, sphering is the whitening transform which moves the data the least. In terms of matrix norms, we can write the mean-squared difference

E x -Wx 2 as E (I N -W) x 2 = Trace (I N -W) Σ x (I N -W) ,
so that, the sphering matrix is the closest to the identity in the matrix norm M 2 Σ = Trace MΣ x M . For later reference, we note that sphering is the default initialization used by the current Infomax implementation in the EEGLAB package [START_REF] Delorme | EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis[END_REF].

Numerical experiments

Comparison of EEGLAB and MNE implementations. Our numerical experiments were conducted with the Infomax implementation of MNE-Python [START_REF] Gramfort | MNE software for MEG and EEG data[END_REF]. We checked that this implementation matches the reference Infomax implementation in EEGLAB [START_REF] Delorme | EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis[END_REF] by reproducing the Infomax results published in [START_REF] Delorme | Independent EEG sources are dipolar[END_REF] based on 13 anonymized EEG datasets (publicly available at http://sccn.ucsd.edu/wiki/BSSComparison [START_REF] Delorme | Independent EEG sources are dipolar[END_REF][START_REF] Onton | Frontal midline EEG dynamics during working memory[END_REF]). This comparison is presented in Table 1. It shows the average across the 13 EEG datasets used in [START_REF] Delorme | Independent EEG sources are dipolar[END_REF]. In this table, "MIR" stands for Mutual Information Reduction [START_REF] Delorme | Independent EEG sources are dipolar[END_REF], whereas "ND 90%" denotes the percentage of ICA components with dipolarity larger than τ = 90. In order to fit a single ECD to a topography, we used the same four-sphere model used in [START_REF] Delorme | Independent EEG sources are dipolar[END_REF] for forward computation. The radius of each sphere was equal to 71, 72, 79 and 85 mm, and their corresponding conductivities relative to the cerebrospinal fluid were equal to 0.33, 1.0, 0.0042 and 0.33, respectively.

In the course of this comparison, we found that EEGLAB does not constrain the fitted dipoles to be located inside the brain, whereas MNE-Python does. For this reason, EEGLAB tends to report an artificially high number of dipolar components as compared to MNE-Python (second row of Table 1). However, we checked that EEGLAB and MNE-Python agree in the number of high dipolar components for dipoles located inside the brain (third row of Table 1). Even better, we checked that they agree on the locations of those dipoles. Step size Evaluation of the stochastic gradient approach. We proceed to evaluate the performance of the SG method used by Infomax. We use subject kb77 from the same study [START_REF] Delorme | Independent EEG sources are dipolar[END_REF][START_REF] Onton | Frontal midline EEG dynamics during working memory[END_REF]. The EEG dataset is composed of 306600 samples of 71 channels sampled at 250 Hz. The algorithm is evaluated by monitoring the step size, as well as the Frobenius norm of the relative gradient after each pass on the full data. We consider the following scenarios, differing by the initial step size α 0 and tolerance tol:

-Infomax a: EEGLAB defaults: α 0 = 6.5 × 10 -4 /log(N ) ≈ 1.5 × 10 -4 and tol = 10 -7 . -Infomax b: α 0 = 10 -3 as in [START_REF] Delorme | Independent EEG sources are dipolar[END_REF], and EEGLAB default value for tol = 10 -7 .

-Infomax c: same α 0 as in scenario "Infomax b" and tol = 10 -14 .

The left panel of Fig. 2 displays the evolution of the relative gradient norm L (W) F across passes on the full data. As we can see in this figure, in none of the scenarios the relative gradient ever go to zero. Scenarios a and b, which use the default tolerance, either reach a plateau or stop early. Yet, reducing the tolerance in scenario c, reveals that it is not sufficient to push convergence compared to b: the trajectory plateaus just after the stopping point for b. However, the two plateaus of cases a and c are of different natures, as revealed by the right panel of Fig. 2, which displays the step size trajectories. For a, the step size remains constant after pass ∼ 83, hence we observe the known plateau of SG methods with fixed step size, whereas for b and c the annealing policy drives the step size to zero exponentially fast, therefore preventing the algorithm to make any further progress.

Fig. 3 shows the evolution of dipolarity of the components across passes on the full data. One can see in Fig. 3(a) that when starting with sphering most of the lines are red which means that they exceed at some point the value of 90. This is due to initialization as explained previously. When initializing Infomax with PCA, see Fig. much less components reach this high dipolarity threshold. In both cases, we observe that the dipolarity stops evolving after approximately 70 passes on the dataset, which is consistent with the plateaus of Fig. 2. The two plots on the right column show the same dipolarity metrics, but this time using the quasi-Newton method known as LBFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno), taking as initialization the unmixing matrices estimated by Infomax. In Fig. 3(b), we can see that two highly dipolar sources according to Infomax leave the region of high dipolarity. In other words, pushing the convergence with LBFGS reduces here the number of highly dipolar components as quantified in [START_REF] Delorme | Independent EEG sources are dipolar[END_REF]. To evaluate the convergence of LBFGS towards a stationary point, we computed the Frobenius norm of the relative gradient at the end of the iterations. While this norm was about 10 -4 after SG, it is about 10 -7 after LBFGS, which confirms that LBFGS does push significantly the convergence. The dipolarity of the components after sphering or PCA, as well as following Infomax and LBFGS in the same four cases as Fig. 3 is presented in Fig. 4. This plot is an extra evidence that simple sphering already yields almost only highly dipolar sources. PCA, on the contrary, contains far less dipolar sources. This is also in line with Fig. 1. This plot also reveals that Infomax followed by LBFGS reaches almost identical dipolarities. This suggests that LBFGS manages to wash out the effect of initialization by converging to the same local minimum.

To gain further evidence, we ran a number of checks. to permutation). Fig. 5 reports that it is far from it, demonstrating that LBFGS deviates non-trivially from the output of Infomax. One can also see that the change operated by LBFGS is larger in the PCA case. In other words, Infomax following PCA brings the estimate further away from a stationary point than the sphering initialization. This conclusion only holds if the same stationary point is reached in both settings. Evidence for this is presented in Fig. 6, where one can see that for this subject (kb77) the estimated unmixing matrix obtained in the PCA condition is close to the inverse of the mixing matrix obtained following sphering (up to a permutation). To assess if there is convergence to the same local minimum when SG is followed by LBFGS, we ran the same computation on all the subjects. Fig. 6 shows that for 8 out of 13 subjects the result perfectly replicates.

Conclusion

We explored the annealing policy, the initial step size and the stopping criterion used by the SG Infomax. We reported results where this algorithmic choices lead Infomax to stop before reaching an accurate stationary point, despite a high number of iterations and long computation times. We explained theoretically why the sphering initialization used by Infomax produces highly dipolar sources. By further pushing the convergence using a quasi-Newton method, we showed that the initialization influences the output of Infomax, hence overestimating the number of highly dipolar sources. This observation could explain why practitioners tend to avoid dimensionality reduction when Infomax is used on EEG. Indeed sphering cannot be used in this case. This paper should be seen as an instantaneous picture on current usage of ICA for EEG data. Given the massive use of such techniques, we hope that it will motivate the development and dissemination of better optimization schemes in this scientific community.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

3. 2

 2 Relationship between initialization and dipolarityICA leads to nonconvex optimization problems, therefore solutions necessarily depend on the initialization used. To illustrate this issue, let us analyze the following scenario. Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal to the identity matrix I N . The topograhies associated to the columns of A 0 are shown in Fig.1(a). Each column of A 0 is putting a value of one at each sensor, and zero in the other ones.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

3. 2

 2 Relationship between initialization and dipolarityICA leads to nonconvex optimization problems, therefore solutions necessarily depend on the initialization used. To illustrate this issue, let us analyze the following scenario. Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal to the identity matrix I N . The topograhies associated to the columns of A 0 are shown in Fig.1(a). Each column of A 0 is putting a value of one at each sensor, and zero in the other ones.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

3. 2

 2 Relationship between initialization and dipolarityICA leads to nonconvex optimization problems, therefore solutions necessarily depend on the initialization used. To illustrate this issue, let us analyze the following scenario. Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal to the identity matrix I N . The topograhies associated to the columns of A 0 are shown in Fig.1(a). Each column of A 0 is putting a value of one at each sensor, and zero in the other ones.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

3. 2

 2 Relationship between initialization and dipolarityICA leads to nonconvex optimization problems, therefore solutions necessarily depend on the initialization used. To illustrate this issue, let us analyze the following scenario. Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal to the identity matrix I N . The topograhies associated to the columns of A 0 are shown in Fig.1(a). Each column of A 0 is putting a value of one at each sensor, and zero in the other ones.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

Fig. 1 .Fig. 1 .Fig. 1 .Fig. 1 .

 1111 Fig. 1. Topographies associated to different initialization.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

3. 2

 2 Relationship between initialization and dipolarityICA leads to nonconvex optimization problems, therefore solutions necessarily depend on the initialization used. To illustrate this issue, let us analyze the following scenario. Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal to the identity matrix I N . The topograhies associated to the columns of A 0 are shown in Fig.1(a). Each column of A 0 is putting a value of one at each sensor, and zero in the other ones.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

3. 2

 2 Relationship between initialization and dipolarityICA leads to nonconvex optimization problems, therefore solutions necessarily depend on the initialization used. To illustrate this issue, let us analyze the following scenario. Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal to the identity matrix I N . The topograhies associated to the columns of A 0 are shown in Fig.1(a). Each column of A 0 is putting a value of one at each sensor, and zero in the other ones.

Fig. 1 .Fig. 1 .Fig. 1 .

 111 Fig. 1. Topographies associated to different initialization.

Fig. 1 .Fig. 1 .Fig. 1 .

 111 Fig. 1. Topographies associated to different initialization.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

3. 2

 2 Relationship between initialization and dipolarityICA leads to nonconvex optimization problems, therefore solutions necessarily depend on the initialization used. To illustrate this issue, let us analyze the following scenario. Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal to the identity matrix I N . The topograhies associated to the columns of A 0 are shown in Fig.1(a). Each column of A 0 is putting a value of one at each sensor, and zero in the other ones.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

3. 2

 2 Relationship between initialization and dipolarityICA leads to nonconvex optimization problems, therefore solutions necessarily d on the initialization used. To illustrate this issue, let us analyze the following sce Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal identity matrix I N . The topograhies associated to the columns of A 0 are sho Fig.1(a). Each column of A 0 is putting a value of one at each sensor, and zero other ones.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

3. 2

 2 Relationship between initialization and dipolarityICA leads to nonconvex optimization problems, therefore solutions necessarily depend on the initialization used. To illustrate this issue, let us analyze the following scenario. Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal to the identity matrix I N . The topograhies associated to the columns of A 0 are shown in Fig.1(a). Each column of A 0 is putting a value of one at each sensor, and zero in the other ones.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

3. 2

 2 Relationship between initialization and dipolarityICA leads to nonconvex optimization problems, therefore solutions necessarily depend on the initialization used. To illustrate this issue, let us analyze the following scenario. Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal to the identity matrix I N . The topograhies associated to the columns of A 0 are shown in Fig.1(a). Each column of A 0 is putting a value of one at each sensor, and zero in the other ones.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

3. 2

 2 Relationship between initialization and dipolarityICA leads to nonconvex optimization problems, therefore solutions necessarily depend on the initialization used. To illustrate this issue, let us analyze the following scenario. Let A 0 2 R N ⇥N denote the mixing matrix at iteration zero. Let A 0 be equal to the identity matrix I N . The topograhies associated to the columns of A 0 are shown in Fig.1(a). Each column of A 0 is putting a value of one at each sensor, and zero in the other ones.

Fig. 1 .

 1 Fig. 1. Topographies associated to different initialization.

Fig. 1 .Fig. 1 .Fig. 1 .Fig. 1 .

 1111 Fig. 1. Topographies associated to differ

Fig. 2 .

 2 Fig. 2. Left: Frobenius norm of the relative gradient vs iterations. Right: Step size α vs iterations.

Fig. 3 .

 3 Fig.3. Evolution for scenario c of dipolarity during Infomax (left) followed by LBFGS (right). Initialization with Sphering (top row) and PCA (bottow row). A line is colored red if it exceeds the value of 90% during the iterations.

Fig. 5 Fig. 4 .

 54 Fig. 4. Dipolarity of the components sorted in decreasing order. Plain lines correspond to sphering initialization while dashed lines correspond to PCA initialization.

Fig. 5 .

 5 Fig. 5. LBFGS further transforms the Infomax components by a matrix T . Left shows log 10 (|T |) using sphering initialization (for display, the sources are sorted to have the largest Tij in the lower left corner). Middle shows the update following PCA initialization. Right plots the rows of the same matrices after sorting each row (red is sphering and black is PCA).

Fig. 6 .

 6 Fig. 6. Changing the initialization changes the local minimum otherwise the transform T linking the sources obtained with the two initializing whiteners (PCA and sphering) would be the identity. Each plot shows log 10 |T |. Left: Infomax, right: LBFGS. By completing the convergence process, we find that in 8 cases out of 13, the resulting sources do not depend on initialization.

Table 1 .

 1 Comparison of EEGLAB and MNE-Python Infomax implementations.

	Metric	EEGLAB MNE-Python
	MIR	43.092901 43.092938
	ND 90%	43.445287 31.744312
	EEGLAB loc. in ND 90% 22.751896 22.751896

Jair Montoya-Martínez, Jean-Franc ¸ois Cardoso, and Alexandre Gramfort

Acknowledgments. This work was funded by the Paris-Saclay Center for Data Science with partial support from the European Research Council (ERC-YStG-676943).