G Lefeuvre

P Gorodetzky
email: philippe.gorodetzky@cern.ch

J Dolbeau

T Patzak

P Salin

Absolute measurement of the nitrogen fluorescence yield in air between 300 and 430 nm

Keywords: nitrogen fluorescence, air fluorescence, extensive air showers, ultra-high energy cosmic rays PACS: 29.30.Dn ; 98.70.Sa

.

Introduction

One of the current challenges in high energy particle physics is to find the origin and nature of Ultra-High Energy Cosmic Rays (UHECR, E > 10 18 eV). The measurement of their energy spectrum could confirm their interaction with the cosmological background (GZK theory). But this is a difficult task: the UHECR are not detectable themselves, but only the air shower they induce while going through the atmosphere.

Experiments which have tried to solve this problem still have non consistent results for E > 5•10 19 eV. They use different detecting techniques. On one hand, AGASA [1,2] has an array of ground detectors separated by about 0.7 km. Muons and electrons of the showers arriving at the ground are sampled. The energy reconstruction is difficult because of the lack of knowledge on hadronic cross sections. This method relays on Monte-Carlo particle distributions in the shower approximation to find the energy of the incoming particle. On the other hand, HiRes [3] uses fluorescence telescopes to detect the continuous development of the shower. The accuracy of this method can still be improved if the uncertainties on the fluorescence itself are lowered. Today, the Auger experiment uses both methods but still calibrates energy extracted from the ground detectors by the fluorescence (also on the ground, [4]). In the future, space-based experiments like JEM-EUSO will look for the fluorescence signal from above [5].

Since 1964, several authors [6,7,8] have measured the fluorescence yield of each emission band of nitrogen (see the spectrum in fig. 5). The fluorescence yield is defined as the number of photons produced when an electron goes through one meter of air. In 2002, the cosmic ray community [9,10,11] decided to start an active campaign of fluorescence measurements for cosmic ray physics. Some were made at high energy in electron accelerators, the energy ranging from 80 MeV to 28 GeV [12,13]. Others were made at low energy, around 1 MeV with a 90 Sr radioactive source, or much lower with electron guns ([14, 15, 18, 21]). These experiments confirmed the hypothesis that the fluorescence yield, up to some tens of MeV, is proportional to the energy deposit of the particle, dE/dX. The difference between energy loss and energy deposit is of no importance from 0.1 to 10 MeV [16,17], which includes our range. Therefore, the absolute scale of the yield can be set with a measurement performed with electrons from a 90 Sr source. Now, the relative variation of the fluorescence yield with altitude is generally considered as quite well known [18]. The parameterization of the yield as a function of altitude is deduced from the relationship between pressure and altitude in an atmospheric model [19,20]. These parameterizations are fairly consistent with one another, but can only give the variation of the yield with altitude and not the absolute value of the yield. Moreover, measurements of the absolute fluorescence yield have never reached a precision better than 13 to 20 % [13,21].

The experiment presented in this paper has been realized with a double purpose: • to improve significantly the precision on the fluorescence yield: it is now 5 % ;

• to achieve a first continuous measurement of the fluorescence spectrum in the range from 300 to 430 nm, which had never been done before with a radioactive source.

Experiment

Principle

Two measurements have been carried out, named "integral" and "spectral" measurements. Electrons from a strontium source are used to excite nitrogen. Three detectors are needed, all of them based on Photonis photomultiplier tubes (PMT). One detects the electrons (electron-PMT), the other two, the photons (photon-PMTs).

On one hand, the nitrogen emits fluorescence light through bands emitting lines from 300 to 430 nm. Thus cosmic ray experiments like HiRes and Auger use band pass filters of that range in front of their detectors. This is what is done here in the integral measurement.

On the other hand, fluorescence yield experiments generally use narrow filters to separate the lines. But the overlap between some spectral filter bandwidth and the close proximity of some nitrogen bands make the separation difficult in reality. For this reason, the spectral measurement has been performed with an optical grating spectrometer, allowing a continuous check and a good separation of the lines.

Description of the setup

This setup has been conceived with the aim of improving the measurement accuracy. The fluorescence chamber is a cross-shaped stainless steel chamber (a 6-ports accelerator pipe cross). It is schematically shown in fig. 1. Electrons are emitted by a strong 90 Sr source (370 MBq) which is placed on the top part of the cross. On the bottom part, a plastic scintillator (truncated cone, 20 and 28 mm diameters, 20 mm thick and its angle is the same that the one made by the electron trajectories at the cone periphery). The electrons path is hence defined as the inside of a cone (100 mm height), with the source at its top. To homogenize the light between the scintillator and the PMT, a kaleidoscopic hexagonal light guide is used to carry the light. The phototube is a Photonis XP2262 [22] equipped with an active voltage divider (the rate is high, of the order of 2.5 MHz, and very stable). The consequence of the high counting rate is a high X-ray noise level, due to the interacting electrons with the surrounding matter, including the PMT. To shield to a maximum these X-rays, a 10 cm diameter lead cylinder has been placed inside the chamber with the 90 Sr source at its centre. A vertical cone is dug in the lead so the electrons can reach the scintillator. This cone is 30 mm diameter at its base, slightly wider than the electrons cone defined by the scintillator. Such a cone minimises the number of interactions between electrons and the lead: either they are totally stopped (around the source), or they touch in a grazing way the cone in their way to the scintillator. A GEANT Monte-Carlo simulation shows that half of the useful electrons do reach the scintillator without an interaction with the lead, the other half touching only once. An horizontal cylinder (40 mm diameter) is also dug in the lead, centred on the optical axis, to let the photons go to the photon-PMTs. This structure is represented in fig. 2.

The effective fluorescence volume is thus, on a first approximation, the intersection of the vertical cone with the horizontal cylinder. In a more precise way, all geometrical efficiencies have been determined with Monte-Carlo simulations. They take into account all possible interactions between electrons and their surroundings and give the exact photon solid angles as well as the exact useful electron path length (the length used in the yield expressed in photons per meter).

Two kinds of interactions can occur: with the gas and with the lead. First, part of the energy of the source electrons is lost by ionization, i.e. creation of secondary electrons. If the secondaries go too far away from the useful fluorescence volume, they can either reach the lead or enter the horizontal cylinder. Both modify the detection efficiencies (geometry wise) of the photon-PMTs calculated. But 99 % of these secondaries have an energy smaller than 5 keV. Their range in air at atmospheric pressure is less than 2 mm. This has been checked through a geometrical simulation not to be enough to have an influence on the solid angles. The effect of these δ rays on the yield expressed in photons per deposited energy is explained later in paragraph 3.1

The second type of interaction concerns the electrons from the source that do not reach directly the scintillator. Some of them are scattered by the lead and then bounce back to the scintillator. This has three consequences, studied by a GEANT Monte-Carlo simulation. The first is the increase in the electron counting rate, due to the larger solid angle available for the electrons: it is doubled by this effect. The second is a modification in the electron energy spectrum, which contains more low energy electrons. This has been taken into account in the calculation of the electrons average energies. The third is a slight change in the useful electron path length, also taken into account.

Some electrons can also bounce from the scintillator back into the fluorescence volume. This has been simulated to be a 10 -9 effect on the photon detection, hence totally negligible. For the integral measurement, fluorescence photons are detected and counted by a Photonis XP2020Q (2 in. diameter) placed on the left hand-side of the setup [23]. This PMT will be named integral-PMT and it is equipped with a fused silica window. As the fluorescence efficiency is known to be low (around 4 photons per meter), it has to work in single photoelectron mode, meaning at a very high gain. It is polarized positively (photocathode at ground) in order to reduce its own noise level to around 300 Hz (instead of 3 kHz with a negative polarity). This effectively removes the tiny discharges in the silica between the photocathode and the outside world. The optical filter on the entrance window is a Schott-Desag BG3 (2 mm thick, 34 mm diameter) [START_REF]Optical filters/Band Pass filter/BG3[END_REF]. It is the same as JEM-EUSO intends to use, and very similar to those used by Auger. It is glued with Epotec 301-2 [25] which has the same refraction index than both the window and the filter. For the spectral measurement, photons are analyzed by an optical grating spectrometer (Jobin-Yvon H25) and counted by another Photonis XP2020Q, attached to the spectrometer. This PMT will be named spectral-PMT. The light rays incoming the spectrometer have to be inside its numerical aperture (f = 250 mm, NA = f/4), so a silica converging lens of 150 mm focal length and 46 mm diameter [START_REF][END_REF] is inserted at the right place between the fluorescence volume and the entrance slit of the spectrometer. The optical image of this volume is thus on the latter. The optical magnification is 1/7. Entrance and exit slits of the spectrometer have the same dimensions: 2 x 7 mm, the maximum size available, inducing a resolution (measured with lasers) of 6 nm FWHM. This part of the setup is shown in fig. 3. The main uncertainty in this kind of experiments, where very few photons are counted, arises from the absolute efficiency of the photon-PMTs themselves. The manufacturers provides this value with an uncertainty of 15 to 20 % (at 1σ), which is not precise enough. For this reason, we designed a new way to measure the absolute efficiency of the photon-PMTs in the single photoelectron mode very accurately, based on the comparison with a NIST photodiode. The relative efficiency map of the photocathodes of the PMTs working in single photoelectron mode was measured with a 377 nm LED every 3 mm in X and Y with a relative precision better than 0.5%. This precision was needed to be better than the point to point variations in efficiency (about 2%), in order to control these variations. A measurement of the absolute efficiency of one of the map points with an accuracy of 1.7%, transforms that relative efficiency map in an absolute efficiency map. A dedicated paper, following a patent, will present this measurement. This result was used to determine the flattest region of the photocathode (see fig. 4). The effective detecting area of the integral-PMT was limited by a 20 mm diameter diaphragm where the efficiency stays roughly constant when the diaphragm size is increased. That way, a small error on the diaphragm diameter has a negligible consequence on the PMT efficiency (it has however on the solid angle, but this is easy to control). The absolute spectral efficiency of the detector {PMT + filter + diaphragm} is the dashed line superimposed on the fluorescence spectrum in fig. 5. Its value for photons at 377 nm is (18.9 ± 0.3)%, and, when convoluted with the fluorescence bands [14] and the relative variations given by Photonis [22], is (17.8± 0.4)%. The spectral-PMT is illuminated on a 10 x 18 mm area (smaller area than the 20 mm circle for the integral PMT, hence all photons at the exit slit of the spectrometer reaching different parts of this area would here also have roughly the same probability to produce photoelectrons).

On their way from the fluorescence chamber to any of the photon-PMTs, integral or spectral-PMT, the 20 cm long tubes are baffled to prevent reflected photons to reach the photocathode.

Gas

The gas used is either nitrogen for the adjustments or Messer dry air for the measurements. Fig. 6 shows the gas setup. It is possible to mix gases and to introduce controlled quantities of different impurities and water vapour. Those measurements including pressure variations are planned for the future.

The gas circulates at a rate of 1 L.h -1 through the chamber to avoid impurities build-up and degassing of the walls. This circulation is controlled by a precise flow meter [27], followed by a pressure controller [28] and an ultra-clean primary pump with an ultimate vacuum of 50 mbar [29]. Internal and external pressures and temperatures are regularly monitored during data taking.

Data acquisition

Hardware data acquisition is basically the same that for other authors, but had to be adapted to high counting rates. Fig. 7 shows its principle.

The lifetimes of excited levels of nitrogen molecules at atmospheric pressure are of the order of a few nanoseconds [8,14,15,18,21]. A measurement based on the time coincidence between an electron and a photon is made to discriminate fluorescence from background photons. Basically, an electron in the scintillator is the trigger for the photon measurements. An electron produces about 0.16 photon in the 4 cm long fluorescence volume. Geometrical efficiencies toward each photon-PMT are very small: 3.69•10 -4 for the integral-PMT and 7.48•10 -6 for the spectral-PMT. As a consequence, only 5.9•10 -5 and 1.2•10 -6 photon respectively reach the PMTs. So, both photon-PMTs work in single photoelectron mode. Their spectra are very "pure" in "1 photo-electron", with a negligible amount of "2". So a discriminator is enough to select the "1 photoelectron" peak. This can be seen in fig. 8.

ADCs (LRS 2249A, 0.25 pC per channel) record the electron and photon spectra. They are used only to check the gain stability of all PMTs during data taking. The "detection inefficiency", D, defined as the proportion of lost fluorescence counts due to the mandatory discriminator threshold on the single photoelectron spectrum (see fig. 8), has also been measured and been found to be 3.76 % of the measured single photoelectron peak. Both ADCs and TDCs could be triggered by the electron signal, but this would produce a very high random rate. We cleaned up the trigger by replacing it by an electron-photon coincidence.

This coincidence is performed by a fast NIM coincidence unit. The photon pulse is 10 ns wide. The electron pulse is 100 ns wide to take into account the delay in photon emission (lifetimes of the levels). This was set to prepare acquisitions at low pressure, where the lifetimes are longer than at atmospheric pressure. Even at very low pressure, 100 ns ensures the loss due to that effect will be negligible. In other terms, a coincidence occurs when the photon arrives with less than 100 ns delay after the electron.

The TDC is stopped by the same electron (see fig. 7), this time as a 10 ns wide pulse and delayed by 150 ns. Therefore, the start has a rate much lower than the stop and this inverts the time axis. Results can be seen in fig. 9. The useful signal is the peak in the spectrum. It can be easily separated from the flat background. The effect on the peak left slope due to the lifetimes of the nitrogen levels cannot be seen on this figure because it has been measured in air at atmospheric pressure. The narrow peak on the far right-hand side is due to random photons arriving some nanoseconds before electrons. These spectra suffer from an important dead time. But the ratio of the peak to total spectrum is not affected. The real number of counts is extracted from scalers measuring the total spectrum with a small and well corrected dead time. These fast scalers (CAMAC CERNSPEC NE003, 25 MHz, and VME V560E, 100 MHz) count all electron, photon and coincidence pulses.

Two thresholds, creating two triggers, are set on the 90 Sr spectrum to study the possible dependence of the fluorescence yield with the electron energy even if it is difficult to imagine an influence other than energy deposit. One is set at ~ 600 keV and another at ~ 1.2 MeV. The corresponding average energies are 1.1 and 1.5 MeV. There are around 1 000 photon pulses per second counted by the integral-PMT (with a background of 300 Hz) and 40 by the spectral-PMT on at the top of the 337 nm line (with a background of 30 Hz). Moreover, the dead times, T M , of both electron scalers have been specifically measured and found to be 1.2 % (scaler rate 2•10 6 Hz) and 0.7 % (scaler rate 6•10 5 Hz). The photon yield per unit length, Y, is derived from the signal portion of the TDC spectrum through the following formula:

Results

Integral measurement

where:

• Π is the number of signal counts in the TDC spectrum (see fig. 9) ; • H is the integral of this spectrum ;

• C is the number of coincidences sent to the TDC measured by a scaler. Only H coincidences have effectively started the TDC due to its dead time ; • 1 + D is the correction to the detection inefficiency ;

• N e is the number of electrons according to the scaler ; • 1 + T M is the correction to the dead time of the scaler itself ; • L moy is the mean length of the electrons path in the fluorescence volume ;

• ε PMT is the efficiency of the integral-PMT ; • T w is the transmission of the fused silica window closing the fluorescence chamber ;

• ε geo is the geometrical efficiency.

Systematic errors are presented in table 1. The main uncertainty of this experiment is due to the high counting rate of the electrons, leading to a non-linear dead time dependence in the TDC module. This effect, which varies from channel to channel, depends on the internal TDC time constants and is not fully understood. It has been evaluated by using all the module's channels, and another TDC module (CAEN V1290N, multihit) to compare their results. All values are found equal within 4% (at 1σ) and this uncertainty has been chosen in a conservative way. In the future, to further reduce this uncertainty, fast flash ADCs will be used, allowing to discriminate pulses with and without pile-up on an event per event basis.

On the contrary, the uncertainty on the efficiency of the integral-PMT is very low. This is due to the photomultiplier efficiency measurement made especially for the fluorescence measurement by the patented method. An example of typical data sample is presented in tab. The pressure and temperature normalizations are then applied. Each excited level has its own lifetime, and the yield can be written with respect to pressure and temperature using the kinetic theory [6]. Thus this normalization has been made for each band using previous parameterizations (which give the same variations for the first kilometers above the ground [20]). Here, Nagano's model [21], who uses the different bands yields is used. The normalized pressure and temperature are those of the US-Standard 1976 model [19] at sea level: 760 mmHg and 15°C. In this model, for an electron energy of 0.85 MeV, one finds that the ratio of the yield at 753.8 mmHg and 295.95° K (the 1.1 MeV conditions) to the yield at 760.0 mmHg and 288.15° K (the US Standard conditions), is 0.9863, and the similar ratio for the 1.5 MeV conditions which are 751.8 mmHg and 296.05° K is 0.9860. These ratios are then used to normalize the measured values.

Finally, the two normalized values are 4.05 and 4.42 ph/m and have, due to this normalization, an added uncertainty of 0.6 %, setting the total relative uncertainty to 5.0 %. They are separated by 8.5%, inside an error bar (± 1σ). Hence, these are averaged to get the fluorescence yield, at 760 mmHg, 15°C and for 0.85 MeV electrons: 4.23 ± 0.21 photons / m. This yield per meter is that of the primary particle, therefore it takes into account its production of δ rays and their fluorescence. This number can also be written in units of photons per deposited energy. For a given energy, these two quantities are strictly proportional, since:

dX dE h Y m ! " # " " $ % & ' ( ) * =
where Φ ν is the fluorescence efficiency at the wavelength corresponding to the frequency ν, with the yield per deposited energy being [31]:

! " # $ % & ' = ( ( ( ) h Y E
The energy of the δ rays is naturally included in the dE/dX of the initial electron. Nevertheless, in the context of this experiment, a special care must be taken for the δ rays. A Monte-Carlo simulation has been performed and shows that here, with the geometry described earlier and at atmospheric pressure, only 1 % of the δ rays have more than 5 keV and some of these could be lost in the lead before having produced any fluorescence. The effective dE/dX should thus be reduced. This simulation shows then that almost 67 % of these δ rays with 5keV or more will produce detectable fluorescence. Therefore only 33 % of them are effectively lost and contribute to the reduction of the dE/dX. Finally, this reduction, hence the reduction in yield, is 0.4 %. The central value is shifted from 0.4 % (we measured 4.21 photons/m before applying this correction). The uncertainty stays the same at 5 %. This would be very different at lower pressures, where a more extensive Monte-Carlo simulation would be required and the correction would be increased.

The energy deposited by a 0.85 MeV electron in a meter of the US Standard air is 0.2059 MeV [30]. The fluorescence yield per deposited energy is thus 20.38 ± 0.98 photons per MeV.

The nitrogen yield to air yield ratio has also been measured, (at one electron energy only: 1.1 MeV), and found to be 4.90 ± 0.01, where the low uncertainty is only statistical: all systematic errors compensate, the gas being the only element changed from one measurement to the other. This ratio is compatible with what has been measured at different energies : 6 ± 2 at 28.5 GeV [13], 5.5 ± 0.3 at 1 MeV [21]

Spectral measurement

Spectral measurements are interesting for many purposes. The gas kinetic theory is valid for individual bands. Is the sum of the different bands yield equal to the integrated yield used by the majority of cosmic ray experiments using fluorescence? How do the different bands change yield when temperature and pressure are modified? Up to now, the different bands have been observed either with a grating spectrometer [6,7,14] in a first method, the electrons being produced by an electron gun. Their energy (around 10 keV) is so small that the electron scattering in the gas prevents any measurement of path length, hence cannot give a yield in photons per meter. The second method [8,15,18,21] uses narrow interference filters to analyze the bands. Here, the electrons, as in this work, come from a 90 Sr source. Fig. 1 of [21] illustrates the method complexity. It is difficult to separate overlapping bands in some filters. The only solution to have a good resolution on the bands is to use a grating spectrometer, and the easiest way to have electrons with an energy large enough to know their path is to use a 90 Sr source. The acceptance of a spectrometer is roughly an order of magnitude smaller than that of a filter. So the very first step is to prove the feasibility of such a method. This is what is done here.

Hence, in this experiment, a 90 Sr source 10 times stronger than in [8,15,18,21] was used, The counterpart of the high number of incident electrons is that a very special attention has to be given to pile-up effects and to the lead shield geometry to avoid the introduction of a large X-ray background. Results of the spectral measurements are represented in fig. 10. The whole fluorescence spectrum has been measured, from 300 to 436 nm, in dry air at room temperature and normal pressure, with electrons of average energy of 1.1 MeV. The time allowed to measure this spectrum was short and the spectrometer equipped with an output slit can measure only one wavelength at a time. So it was decided to open the slits at their maximum, setting the 1σ resolution to 3 nm. Thus measurements were made with steps of 3 nm. But even with 2•10 6 electrons detected every second, only 0.16 fluorescence photons per second are recorded in the signal part of the TDC spectrum for the most intense line (337 nm). This is the reason why this spectrum has only been made for the "low" energy threshold, as defined above.

Fluorescence lines are much narrower than the resolution set for the spectrometer as is seen in [14]. The absolute yield is thus given by the height of the curve at a given wavelength, and not by the area integrated over a bandwidth dλ. The sum of the yields of all lines indicated in fig. 10 is 3.9 ± 0.8 photons per meter. This value is in close agreement with the result of the integral measurement. The 20% uncertainty is due most entirely to the bad knowledge of the spectrometer absolute efficiency, especially in the UV.

Superimposed on this spectrum is Ulrich's spectrum [14], convoluted with a 3 nm resolution. As it is a relative yield spectrum, it is shown here normalized to the 337 nm line. Our spectrum is well compatible with what would measure [14] with such a resolution. The first observation is the confirmation of the feasibility of this measurement: all the main lines are indeed observed. Evidence is given that an absolute spectrum can be measured with a basic apparatus, and taking care of reducing the PMT background as explained earlier. The sum of the lines yield is consistent with the previous integral measurement. Moreover, its uncertainty is 20 %, when previous experiments do not give better than 15 %.

The second observation concerns the discrepancies, around 316 nm and beyond 375 nm. The extraction of the absolute yield involves the efficiency of the spectrometer. Its spectral efficiency curve is provided by the manufacturer with a low accuracy in the UV (which is true also for Ulrich's results). Furthermore two effects induce a very important loss in our spectrometer:

• the ageing of each optical element (mirrors and grating). The spectrometer efficiency can be reduced to a value around 15 % after ten years of use ; • the extensive use with intense UV light before this experiment.

At 400 nm, the absolute efficiency of the spectrometer was measured with an accuracy of about 2% to be only 15 %, instead of the 61 % given by the manufacturer (who confirms in a private communication such a low value compatible with ageing). The method used is very similar to that taken to determine the absolute efficiency of the photon-PMTs (comparison to a NIST photodiode). This value of 15 % has been used to calculate the absolute yields in the entire spectral range. But there is no reason that this loss is constant with respect to wavelength. The 20% uncertainty on the yield arises from this unknown but limited variation. The question of getting the absolute efficiency of spectrometers in the UV is quite challenging and is the object of specific attentions by the community of atomic / spectral physics.

It is unfortunate that this spectrometer method does not yield yet an accuracy better than the "integrated yield" method. It will if the spectrometer resolution can be made high enough to totally separate the bands, which is possible if it is equipped with a CCD readout to minimize the experiment duration. Then, the absolute efficiency of the spectrometer will have to be determined with a high accuracy in the UV, not an easy task according to Ulrich, but possible through our patented method.

In the future, two measurements will be done: • calibrate in an absolute way the old Jobin-Yvon;

• use the new spectrometer able to measure the full spectrum at once with a 0.1 nm resolution and equipped with a light intensifier to measure the fluorescence yield.

Conclusion

The absolute fluorescence yield of nitrogen in dry air at atmospheric pressure has been measured. The precision of the measurement is improved by a factor of three, which has an immediate impact on the cosmic ray energies found by HiRes which uses Bunner's [6] yield. Their energies are increased by 22%, hence their spectrum is much closer to AGASA's (which incidentally have been recently lowered the energy of their points by 10% [START_REF] Teshima | Ultra High Energy Cosmic Rays Observed by AGASA[END_REF]).

The first continuous fluorescence spectrum of nitrogen excited by electrons from a 90 Sr source was also measured.

Next steps are to introduce impurities in the gas, such as argon, water vapour and pollutants. A pressure study of the total yield will be made. On another hand, the spectral measurement will be improved thanks to a new spectrometer. Papers will follow to account for these measurements, which will provide an overall and realistic view of the fluorescence phenomenon.

Figure 1 .

 1 Figure 1. Schematic view of the fluorescence chamber. The strontium source and the scintillator are inside the lead shield. The internal structure of the lead is shown in fig. 2 (side view).

Figure 2 .

 2 Figure 2. Internal structure of the lead shield. Inside the hollow cone is the smaller cone of the path of the useful electrons defined by the source and the scintillator (thick dotted line).

Figure 3 .

 3 Figure 3. Schematic view of the optical components for the spectral measurement.

Figure 4 .

 4 Figure 4. Absolute efficiency of the integral-PMT (dashed line) and spectral-PMT (bold line) as a function of the size of the effective detecting surface size. A 20 mm diameter diaphragm is used to limit the effective detecting area of the integral-PMT to the flattest region of the photocathode. Using a diaphragm of 20 mm diameter (314 mm 2 ), the efficiency is 18.7%. The spectral-PMT is used in its central region (~ 10×18 mm) where the efficiency including the filter is 19.3%. The uncertainty is 1.7% of the efficiency values. Note that the zero suppression on the vertical axis.

Figure 5 .

 5 Figure 5. Fluorescence spectrum (vertical lines) with relative intensities and absolute response of the photon detector (PMT + filter + diaphragm, dashed curve) (source of the fluorescence spectrum : [14]). The integral-PMT efficiency for this spectrum is 17.8%.

Figure 6 .

 6 Figure 6. The gas circulation scheme.

Figure 7 .

 7 Figure 7. The data acquisition principle

Figure 8 .

 8 Figure 8. Single photoelectron spectrum, with and without the pedestal.

Figure 9 .

 9 Figure 9. Example of TDC spectrum for the integral measurement

2 .

 2 The low (1.1 MeV) and high (1.5 MeV) energy measurements give respectively 3.95 and 4.34 photons per meter. Statistical uncertainties are 0.2 % and 0.8 %. The energy normalization is made at 0.85 MeV using the dE/dX ratios with values interpolated from NIST data [30] and yields respectively 4.05 and 4.41 photons per meter.

Figure 10 .

 10 Figure 10. Spectrum of the absolute fluorescence yield of nitrogen in air between 300 and 436 nm. The bold line is the result of the present experiment, and the dashed line is Ulrich's spectrum [14] as if analyzed by our spectrometer. Discrepancies could be explained by the ageing of our spectrometer.

Table 1 .

 1 Systematic uncertainties of the experiment.

Table 2 .

 2 Data sample for the low energy measurement.

Acknowledgements

We would like to thank Bernard LEFIEVRE for his help with GEANT simulation of the setup, and François LELONG and Jean-Paul RENY for their help in building the bench.


domain: spi math

text:
Judit Chamorro-Servent

Laura Bear

Josselin Duchateau

Mark Potse

Rémi Dubois

Yves Coudière

Possible improvements of the method of fundamental solution to solve the ECGI problem

This study received financial support from the French Government as part of the « Investissement d'avenir » program managed by the National Research Agency (ANR), Grant reference ANR-10-IAHU-04 from the Conseil Régional Aquitaine as part of the project « Assimilation de données en cancérologie et cardiologie ».

 Removing the Neumann conditions reduces the forward and inverse computational cost.  MFS with new distribution of sources made the problem less ill-conditioned and less sensitive to the choice of the regularization parameter.  The HNC are related with the small SVs, i.e. with higher frequencies, this can be seen in the respective noisy reconstructed potentials.  MFS with new distribution of virtual sources and without HNC has better CC and less REs. However, the REs are still high. This indicates that a good similarity of the potential pattern is achieved but not a good amplitude. Finally, best ATs maps were achieved by standard MFS without HNC.  Future work: The amplitude can be improved with better choice of the regularization parameter according to the distance of the virtual sources and/or by using spatio-temporal regularization.

(1)

  ,in . Under:

Dirichlet conditions: = on Homogeneous Neumann conditions: 0 on

T T n T div          
Geometrical meshes used: The body is cut at the top and bottom of the torso and the arms.  Despite all the success of ECGI, the understanding and treatment of many cardiac diseases is not feasible yet without an improvement of the inverse problem solution.

 Can we improve the forward problem to make the inverse problem less ill-conditioned?

 In other fields it has been shown that the ill-posedness of MFS problem and the efficacy of its solution depends on the placement of the virtual sources [3].  The center of the heart used in [2] to inflate the torso surface (in order to place its respective virtual sources) can be a misfit reference, for example when structural diseases are present.

 Therefore, we studied the effect of changing the torso virtual placement: i. For each electrode on the torso, we find the location of closest point on the heart surface ii. We construct the location of virtual source related to the torso electrode as iii. We use SVA to find the ratio that makes our problem better Conditioned, i.e. that improves the SV decay.

      2 2 2 2 0 1 2 2 2
1 1 , for 0 (MFS without HNC) and 0.5

(Standard MFS) 2               T J a B a B a a (4)    min min 2 s.t. min , 1, , d i k H h d x h k M       min i i T i d y x R x h    i y i x T R                    
                                                                                                                     ^, 1 4 i T j T E i j i j x y f x y x y     0 1 , ^, M E j j j E j T E
x a a f x y  The MFS matrix can be split in two matrices:

 We can remark that:

-The contribution of the submatrices B 1 and B 0 in the whole solution is related to their respectives norms. -Due to the meshing of the torso, non-physiological conditions are applied to the top, bottom of the torso and the arms. -Accurate computation of the normals at torso surface is required.

Do we need to enforce the homogeneous Neumann conditions in the MFS inverse problem?

•To evaluate the contribution of both boundary conditions:

The solution of (4) was found by using Tikhonov regularization and was chosen by using Creso method as in [2].  Study of the effect of changing the torso virtual sources placement • The new distribution of torso virtual sources reduces the ill-posedness of the problem

Study of placement of virtual sources

Method of evaluation

COMPARISON OF ELECTROGRAMS

Table 3. CCs in the interval (beat) outlined by the orange box and REs for potentials below (figure 5).

  Barr R, et al. Relating epicardial to body surface potential distributions by means of transfer coefficients based on geometry measurements. IEEE Trans Biomed Eng Jan 1977; BME-24 (1): 1-11. [2] Wang Y, Rudy Y. Application of the method of fundamental solutions to potential-based inverse electrocardiography. Annals of Biomedical Engineering Aug 2006; 34(8): 1272-1288. [3] Karageorghis A, et al. A survey of applications of the MFS to inverse problems. Inverse problems in science and Engineering 2011; 19(3): 309-336. [4] Potse M, et al. Cardiac anisotropy in boundary-element models for the electrocardiogram. Medical & Biological Engineering & Computing. 2009, 47(7): 719-729.



  Robust calculations of the inverse ECGI problem may require accurate specification of boundary conditions at the torso and cardiac surfaces [1]. The aim of the inverse problem is to find  The classical formulation of the ECGI inverse problem with the method of fundamental solution (MFS) involves a linear system [2].



  To test the effect of HNC -Five activation patterns (1 single site pacing and 4 single spiral waves) were simulated using [4]: i. The currents from the monodomain problem are used to compute the BSPM. The simulations provide both the theoretical, in-silico and every 1ms. ii. were computed by solving a static bidomain problem in a realistic, heterogeneous torso model.  To test the ill-posedness of the problem with standard and new distribution of sources -Real data from a patient presenting scar was employed.  To compare the four combinations: Standard MFS, Standard MFS without HNC, MFS with new distribution of sources, and MFS with new distribution of sources and without HNC -Experimental data (LBBB) from: for epi /endo surfaces 108 electrodes at sock Heart and sock meshes  Study of Boundary conditions: Respective contribution of B0 and B1 • The depends only on the geometry and location of sources (for which we used a fixed rule defined in [2]). We found , indicating predominance of B 0 compared to B 1 in the standard MFS.

Figure 1 .

 1 Figure 1. Activation maps for electrograms, single site pacing.

Figure 3 .

 3 Figure 3. Reconstructed potentials for different parameters of regularization.(Data from patient with scar)

Figure 4 .

 4 Figure 4. Activation times maps. From left to right: signals recorded at sock, Standard MFS, MFS with new distribution of sources, Standard MFS without HNC, MFS with new distribution of sources and without HNC.

Figure 5 .

 5 Figure 5. Potentials reconstructed by the 4 combination of methods against potentials recorded on the sock.

Figure 2 .

 2 Figure 2. DPC plot of standard MFS against MFS with new distribution of sources for two instances of time different. (Data from patient with scar)

Table 1 .

 1 CC and RE for electrograms, five simulations.

Table 2 .

 2 CCs and REs for Ats maps above (figure4).


domain: info.info-au math spi sdv

text:
Francis Mairet

Olivier Bernard

Robustness of Bioprocess Feedback Control to Biodiversity *

Keywords: Process Systems Engineering Biotechnology, stability analysis, nonlinear systems, bioreactor, multispecies, coexistence

The design of control laws for bioprocesses are generally based on simplified single-species models. Biodiversity (made of a mixture of different species or strains) is nonetheless inherent in any artificial ecosystem in a bioreactor. Here we propose to define and study the robustness to biodiversity of bioprocess control laws: given a control law designed for one species, what happens when additional species are present? We illustrate our approach with a well used control law which regulates substrate concentration using measurement of growth activity. Depending on the properties of the additional species, the control law can lead to the required objective, but also to an undesired monospecies equilibrium point, coexistence, or even a failure point. Finally, we show that, for this case, the robustness can be improved by a saturation of the control. Robustness to biodiversity is a difficult issue which should be better accounted for in the control design.

Introduction

Biodiversity is inherent in any bioprocess, even when the use of a single microorganism is initially targeted. This diversity, desired or endured, can involve a lot of different species (e.g. several hundreds of archaea and bacteria species in anaerobic digestion [START_REF] Rivière | Towards the definition of a core of microorganisms involved in anaerobic digestion of sludge The[END_REF]). For monospecific cultures, different strains can be initially present or resulting from natural mutation within the process (e.g. in pharmaceutical biotechnology). Nevertheless, theoretical developments in automatic control for bioprocesses are generally based on the single-species dogma, i.e. one macroscopic reaction involves one species (see e.g. [START_REF] Julio | Integrated nonlinear optimization of bioprocesses via linear programming[END_REF][START_REF] Dochain | Automatic control of bioprocesses[END_REF][START_REF] Neli | Nonlinear adaptive control of a model of an uncertain fermentation process[END_REF]). A few works have shed light on the importance of biodiversity for the modelling and control of bioprocesses. For example, the effect of multispecies has been shown in simulation for anaerobic digestion [START_REF] Ramirez | Modeling microbial diversity in anaerobic digestion through an extended ADM1 model[END_REF]. It has also been tackled in the development of optimal strategies (via simulation for the start-up of an anaerobic digester [START_REF] Sbarciog | Some Considerations About Control of Multi-species Anaerobic Digestion Systems[END_REF] and analytically for fed-batch operation [START_REF] Gajardo | Minimal time sequential batch reactors with bounded and impulse controls for one or more species[END_REF]).

Classically, robustness of control laws are considered with respect to parameter uncertainties or noise measurements (see e.g. [START_REF] Tao | Adaptive extremum seeking control of continuous stirred-tank bioreactors[END_REF][START_REF] Mailleret | Robust Nonlinear Adaptive Control for Bioreactors with Unknown Kinetics[END_REF][START_REF] Jean-Luc | Robust control for an uncertain chemostat model International Journal of Robust and Nonlinear Control[END_REF][START_REF] Neli | Nonlinear adaptive control of a model of an uncertain fermentation process[END_REF]), but to our knowledge, never to the presence of other species having their own dynamics. Indeed, given a control law designed for one species, what happens when two or more species are actually present? To tackle this problem, we introduce the concept of robustness to biodiversity. A control law is said to be robust to biodiversity for a given set of species if the performance of the control law (in a sense that will be defined later) is not reduced by the presence of any species from this set.

To define the robustness to biodiversity of a control law, one should study the asymptotic behavior of a system including not only a single species, but an assemblage of species competing for a substrate. This has been widely done in mathematical ecology in order to study the competition of species [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF][START_REF] Frédéric | Competition between phytoplankton and bacteria: exclusion and coexistence Journal of mathematical biology[END_REF]. When a constant dilution rate is applied, the principle of competitive exclusion for Monod-type models, states that at most one species will survive. On the other hand, [START_REF] Leenheer | Feedback control for chemostat models[END_REF][START_REF] Gouzé | Feedback control for nonmonotone competition models in the chemostat Nonlinear Analysis[END_REF][START_REF] Mazenc | Further results on stabilization of periodic trajectories for a chemostat with two species[END_REF] have proposed control laws in order to obtain coexistence.

Finally, recent developments have been proposed for the selection of species [START_REF] Mairet | Driving Species Competition in a Light-limited Chemostat in 9th[END_REF][START_REF] Bayen | Optimization of the separation of two species in a chemostat[END_REF][START_REF] Francis | The photoinhibistat: operating microalgae culture under photoinhibition for strain selection[END_REF], with the objective to design a control strategy in order to select a species of interest. These papers propose some useful tools which will be used to study robustness to biodiversity.

After defining the concept of robustness to biodiversity, we illustrate our approach with a control law proposed in [START_REF] Mailleret | Robust Nonlinear Adaptive Control for Bioreactors with Unknown Kinetics[END_REF] in order to regulate substrate concentration. Finally, we propose a modification of the control law in order to increase its robustness.

Framework and definition

We consider a generic model of micro-organism growth limited by one substrate in a chemostat [START_REF] Dochain | Automatic control of bioprocesses[END_REF]. Denote x 1 and s the concentrations of biomass and substrate respectively. The chemostat model is given by:

     ṡ = u(t)(s in -s) -k 1 µ 1 (s)x 1 ẋ1 = (µ 1 (s) -u(t)) x 1 (1)
where u(t) is the dilution rate (which will be our control signal), µ 1 (s) the specific growth rate, s in the input substrate concentration, and k 1 the pseudo yield coefficient.

A species is defined by a parameter vector θ ∈ S ⊂ R n θ + , gathering stoechiometric and kinetic parameters (i.e. parameters of the growth rate µ ). We will refer to S as a set of species.

Denoting ξ the state vector, we consider a control law u(ξ) which globally stabilizes System (1) towards a set-point ξ * = (s * , (s in -s * )/k 1 ) associated to a scalar productivity criterion P(ξ) (often maximising this criterion). For example, P(ξ) can be related to the biomass productivity at the reactor outlet (P(ξ) = ux 1 ) or to the substrate removal rate (P(ξ) = u(s in -s)). Note that the productivity is often zero for u = 0 or for the species washout (x 1 = 0 and

s = s in ).
By definition, the robustness to biodiversity of this control law guarantees that, in presence of multispecies, the system will be stabilized at a point ξ with the same or better productivity P(ξ) than for the set-point in monoculture ξ * .

Definition 1. The control law u(ξ) is said to be (S, n, P) robust iff, for any n species x i , i = 2, ..., n + 1 such that θ i ∈ S, the control law globally stabilizes the system:

     ṡ = u(ξ)(s in -s) - n+1 i=1 k i µ i (s)x i ẋi = (µ i (s) -u(ξ)) x i , i = 1, 2, ..., n + 1, (2)
towards a point ξ such that P( ξ) ≥ P(ξ * ), where ξ * = (s * , (s in -s * )/k 1 , 0, ..., 0).

A didactic example

To highlight the concept of robustness to biodiversity, we consider the control law proposed in [START_REF] Mailleret | Robust Nonlinear Adaptive Control for Bioreactors with Unknown Kinetics[END_REF]. This must be seen as a simple didactic example to introduce the ideas and show their relevance for bioprocess management.

Control law design

In the following, we assume the following hypotheses on the growth rate and the measurement: Hypothesis 1. (i) The specific growth rates µ i are nonnegative functions C 1 of s with µ i (0) = 0.

(i)The measurement of total growth is available:

y(ξ) = n+1 i=1 l i µ i (s)x i
where l i is the associated yield coefficient for species i.

In anaerobic digestion, the methane flow rate measurement is directly related to the growth rate, given the low solubility of methane (see e.g. [START_REF] Mailleret | Robust Nonlinear Adaptive Control for Bioreactors with Unknown Kinetics[END_REF][START_REF] Karafyllis | A vector Lyapunov function characterization of input-to-state stability with application to ro-bust global stabilization of the chemostat[END_REF]).

For other bioprocesses, the total growth can often be estimated using observerbased estimator [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF][START_REF] Perrier | Tuning of observerbased estimators: theory and application to the on-line estimation of kinetic parameters Control Engineering Practice[END_REF] using for example the measurement of gaseous O 2 or CO 2 .

Given a set-point s * ∈ (0, s in ), we consider the feedback law proposed by [START_REF] Mailleret | Robust Nonlinear Adaptive Control for Bioreactors with Unknown Kinetics[END_REF]: Theorem 1. Under Hypotheses 1, the feedback control law

u(ξ) = γy(ξ) (3)
with γ = k1 l1(sin-s * ) globally stabilizes System (1) towards the positive set point (s * , x * 1 ), where

x * 1 = (s in -s * )/k 1 .
Proof. See [START_REF] Mailleret | Robust Nonlinear Adaptive Control for Bioreactors with Unknown Kinetics[END_REF].

To illustrate our presentation, we will consider the productivity criterion defined as the treated organic load (for waste-water treatment):

P(ξ) = u(ξ)(s in -s) (4)

Adding another species

Robustness of Control law (3) with respect to fluctuations in s in has been considered in [START_REF] Karafyllis | A vector Lyapunov function characterization of input-to-state stability with application to ro-bust global stabilization of the chemostat[END_REF]. Here, we will study its robustness to biodiversity. For sake of simplicity, we restrict our analyses to the presence of only one additional species.

Recalling that the dynamics of M := s + k i x i is given by

Ṁ = u(ξ)(s in -M ),
we additionally consider initial conditions (s, x 1 , x 2 ) in the attractive invariant manifold {(s,

x 1 , x 2 ) | s + k 1 x 1 + k 2 x 2 = s in }. Under this hypothesis, System
(2) with Control law (3) becomes:

ẋi =   µ i (s) -γ 2 j=1 l j µ j (s)x j   x i , i = 1, 2. ( 5
)
with s = s in -2 i=1 k i x i . System (5) admits the following equilibria:

• Substrate depletion (s = 0): E 0 = {(x 0 1 , x 0 2 ) ∈ R 2 + | k 1 x 0 1 + k 2 x 0 2 = s in }.
This corresponds to an equilibrium line.

• Washout (s = s in ): E w = (0, 0).

• Monospecies 1(s = s * ): E 1 = (x * 1 , 0), where x * 1 = 1 γl1 • Monospecies 2 : E 2 = (0, x * 2 ), where x * 2 = 1 γl2 iff k 2 l 2 < γs in (C1)
The substrate concentration at equilibrium is s = s * 2 := s in -k2 l2γ . Note that (C1) guarantees that s * 2 > 0.

• Coexistence : this equilibrium point is possible if the two growth functions have an intersection point:

∃s c ∈ (0, s in ) | µ 1 (s c ) = µ 2 (s c ). (C2) E c = (x 1 , x2
) exists iff x1 and x2 are positive, i.e. iff:

x1 = l 2 (s c -s * 2 ) k 2 l 1 -k 1 l 2 > 0 and x2 = l 1 (s c -s * ) k 1 l 2 -k 2 l 1 > 0. ( C3
)
These conditions impose that the substrate concentrations of monospecies equilibria are located on either side of the intersection point of the growth functions:

-if l2 k2 > l1 k1 , then s * < s c < s * 2 . -if l1 k1 > l2 k2 , then s * 2 < s c < s * (recalling that s * 2 ≤ 0 if E 2 does not exist).
For each equilibrium E • , we denote P • the corresponding productivity.

Local stability analysis

The Jacobian matrix J of System (5) in the invariant manifold is:

J =    B 1 -γ(y(ξ) + x 1 A 1 ) -k 2 µ ′ 1 (s)x 1 -γx 1 A 2 -k 1 µ ′ 2 (s)x 2 -γx 2 A 1 B 2 -γ(y(ξ) + x 2 A 2 )   
where

A i = l i µ i (s) -k i (l 1 µ ′ 1 (s)x 1 + l 2 µ ′ 2 (s)x 2 ) and B i = µ i (s) -k i µ ′ i (s)x i . 3.3.1 Local stability of E 0 .
The eigenvalues of the Jacobian matrix J evaluated at E 0 are λ 1 = 0 and

λ 2 = µ ′ 1 (0)x 0 1 (γs in l 1 -k 1 ) + µ ′ 2 (0)x 0 2 (γs in l 2 -k 2 ). Given that γs in l 1 > k 1 , E 0 is unstable if γs in l 2 > k 2 , i.e. if E 2 exists (C1). If the second term is negative, recalling that E 0 = {(x 0 1 , x 0 2 ) ∈ R 2 + | k 1 x 0 1 +k 2 x 0 2 = s in },
there will be a subset of E 0 where λ 2 < 0, which corresponds to locally stable equilibria.

Stability of E w .

At E w , the Jacobian matrix J becomes:

J =    µ 1 (s in ) 0 0 µ 2 (s in )    so E w is unstable. 3.3.3 Local stability of E 1 and E 2 .
The eigenvalues of the Jacobian matrix J evaluated at

E 1 are λ 1 = -µ 1 (s * ) and λ 2 = µ 2 (s * ) -µ 1 (s * ). Thus, E 1 is locally stable iff µ 1 (s * ) > µ 2 (s * ). (C4)
Similarly, E 2 is locally stable iff

µ 2 (s * 2 ) > µ 1 (s * 2 ). (C5) 3.3.4 Local stability of E c .
At the equilibrium point E c , the trace and determinant of J writes:

tr(J) = -µ 1 (s c ) -γ x1 x2 [µ ′ 2 (s c ) -µ ′ 1 (s c )] (k 2 l 1 -k 1 l 2 ) det(J) = µ 1 (s c )γ x1 x2 [µ ′ 2 (s c ) -µ ′ 1 (s c )] (k 2 l 1 -k 1 l 2 )
Thus, E c is locally stable iff

[µ ′ 2 (s c ) -µ ′ 1 (s c )] (k 2 l 1 -k 1 l 2 ) > 0. (C6)

Robustness to biodiversity of Control law (3)

We now consider the following assumption on the growth functions:

Hypothesis 2. The specific growth rates µ i (s) are assumed to be of Monod type

µ i (s) = μi s s+Ki
, where μi and K i are respectively the maximum growth rate and the half-saturation constant.

Remark that only the ratio of the yield coefficients α i := l i /k i matters to compute local stability and productivity at steady state. As a consequence, a species can be defined by the set of three parameters

θ i = (μ i , K i , α i ),
S := [μ -, μ+ ] × [K -, K + ] × [α -, α + ].
Now, we can define the largest set 1 S R ⊂ S such that Control law ( 3) is

(S R , 1, P) robust. Proposition 1. Given initial conditions (x 0 1 , x 0 2 ) ∈ R * + 2 such that k 1 x 0 1 +k 2 x 0 2 <
s in (i.e. s 0 > 0), the asymptotic behavior of System (5) is given in Table 1.

Proof. First, note that System ( 5) is dissipative: all positive trajectories lie in

the bounded set {(x 1 , x 2 ) ∈ R 2 + | k 1 x 1 + k 2 x 2 ≤ s in }.
If (C1) does not hold, a subspace of E 0 is locally stable (so the control law is not robust to any P which cancels when u cancels, and in particular for ( 4)).

Now, we assume that (C1) holds, and we consider the subset {θ ∈ S |

γs in α 2 > 1}. Let z = l 1 x 1 + l 2 x 2 .
In the coordinates (z, x 1 ), System [START_REF] Ramirez | Modeling microbial diversity in anaerobic digestion through an extended ADM1 model[END_REF] becomes:

     ż = y(ξ)(1 -γz) ẋ1 = (µ 1 (s) -γy(ξ)) x 1 . (6)
Given that E 0 and E w are repulsive, y(t) cannot tend to zero, and thus, z tends to γ -1 , which shows the absence of periodic solutions or cycles. Therefore, any trajectory converges to an equilibrium point. The previous studies of equilibria and local stability allow us to conclude for the different cases.

The asymptotic behavior of System ( 5) is illustrated on Fig. 1 for different cases. The control law is robust whenever all the stable equilibria lead to a productivity greater or equal to P 1 (see Table 1 and Appendix).

Results

In order to restrict our analysis to a more realistic case, we consider a subset of S more relevant from an ecological point of view. We consider that each parameter represents a trait that is involved in the fitness of the species. A

1 note that we ignore subsets of null measure corresponding to non-generic cases, since they will never appear in practice (for example the cases s * = sc, or s * 2 = sc). [START_REF] Shoval | Evolutionary trade-offs, Pareto optimality, and the geometry of phenotype space[END_REF]. In our example, we consider three archetypes x μ, x K , x α , with

θ μ = (μ + , K + , α -), θ K = (μ -, K -, α -)
, and θ α = (μ -, K + , α + ). Note that the best value for K is K -. These archetypes define the subset S:

S := {θ = aθ μ + bθ K + cθ α , ∀(a, b, c) ∈ [0, 1] 3 | a + b + c = 1}
Thus, for any triplet (a i , b i , c i ) ∈ [0, 1] 3 | a i + b i + c i = 1, the parameter vector of the species x i (a i , b i , c i ) is given by:

θ i =       a i μ+ + (1 -a i )μ - b i K -+ (1 -b i )K + c i α + + (1 -c i )α -       (7)
The morphospace, i.e. the space of trait values, is represented on a ternary plot in Figure 2. Each vertex of the triangle represents an archetype. We consider Control law (3) designed for a species x 1 (0.3, 0.4, 0.3) located near the center of the morphospace. After discretizing the morphospace, we test for each additional species which conditions (C1, C4, C5) hold in order to determine the asymptotic behavior of the system (see Proposition 1 and Table 1). This allows us to define S R (represented in color on Figure 2), a subset of S, such that Control law ( 3) is (S R , 1, P) robust. First, one can see that the presence of an additional species can increase the productivity. On the other hand, the control law is not robust for a large subsets of S, corresponding mainly to an additional species with a smaller yield coefficient than species x 1 (α 2 < α 1 ).

Two situations may occur: the productivity at a stable equilibrium point (E 2 or E c ) is smaller than the productivity P 1 (the white area), or there is a reactor shutdown (s and y converge towards zero, the black area). A small decrease in productivity can actually be tolerated (although it is not considered in the present definition of robustness to biodiversity). On the other hand, the reactor shutdown is much more problematic and represents a real drawback of Control law (3). This motivates new control law design in order to avoid as far as possible such situation, i.e. to increase the robustness to biodiversity of the control law.

4 Robustification of the control law

Saturated control law

In order to avoid reactor shutdown, we propose to saturate the dilution rate when growth activity is low. The feedback law proposed by [START_REF] Mailleret | Robust Nonlinear Adaptive Control for Bioreactors with Unknown Kinetics[END_REF] becomes:

u(ξ) =      γy(ξ) if y(ξ) > y γy if y(ξ) ≤ y (8)
where the design parameter y must verify:

y < min{µ 1 (s in ), µ 1 (s * )} γ . (9)
This condition guaranties that the washout is repulsive, and that there is no non-trivial equilibrium with y ≤ y. x μ(1, 0, 0), x K (0, 1, 0), x α (0, 0, 1). The diamond represents species x 1 (0.3, 0.4, 0.3). The set of species S R , such that Control law ( 3) is (S R , 1, P) robust, is represented in color. The color map represents the relative increase of productivity (with respect to P 1 ). The black area corresponds to reactor shutdown (Case A).

Theorem 2. Under Hypotheses 1, the feedback control law (8) globally stabilizes System (1) towards the positive set point (s * , x * 1 ).

Proof. Consider M = s + k 1 x 1 whose dynamics is:

Ṁ = u(ξ)(s in -M ). ( 10
) It follows that the set Ω = {(s, x) ∈ R 2 + | M ≤ max{M (0), s in }} is positively invariant.
In open loop with a constant dilution rate D < µ 1 (s in ), it is well known that System (1) has a non-trivial equilibrium point which is globally asymptotically stable (GAS) [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF]. Given Condition [START_REF] Mailleret | Robust Nonlinear Adaptive Control for Bioreactors with Unknown Kinetics[END_REF], for a constant dilution rate γy, the system has a GAS equilibrium (s, (s in -s)/k 1 ), where s is the smallest root of γy = µ 1 (s). At this point, we have:

y(s, (s in -s)/k 1 ) = α 1 γy(s in -s) = y s in -s s in -s * .
Given that s < s * , we get y(s, (s in -s)/k 1 ) > y, so this point is located outside the set Ω = {(s, x) ∈ Ω | y(ξ) ≤ y}. Thus, under Control law ( 8), all trajectories with initial conditions in Ω will leave this set. Now, in Ω := Ω \ Ω, System [START_REF] Rivière | Towards the definition of a core of microorganisms involved in anaerobic digestion of sludge The[END_REF] rewrites:

     ṡ = γy(ξ)(s * -s) ẋ1 = γy(ξ)(x * 1 -x). (11)
Given that y(ξ) is lower bounded and that y(s * , x * 1 ) > y, (s * , x * 1 ) belongs to Ω and it is the unique stable equilibrium point in this set. Thus, System (1) under Control law (8) has a unique stable equilibrium point (s * , x * 1 ) in the positively invariant set Ω. Equation [START_REF] Jean-Luc | Robust control for an uncertain chemostat model International Journal of Robust and Nonlinear Control[END_REF] eliminates the possibility of periodic orbit, so

Poincaré-Bendixson theorem allows to conclude that (s * , x * 1 ) is GAS.

The global stability of the control law relies on the transition between regions: every trajectory starting in Ω will enter in Ω. Such transitions between regions where a constant dilution rate is applied has supported the design of an hybrid control law in [START_REF] Mairet | Hybrid Control of a Bioreactor with Quantized Measurements IEEE Transactions on Automatic Control[END_REF] where quantized measurements y are considered.

Stability analysis with multispecies

We will now characterize how the saturation of the control law increases its robustness to biodiversity. In line with Section 3, we consider the presence of one additional species with initial conditions (s, x 1 , x 2 ) in the attractive invariant manifold {(s,

x 1 , x 2 ) | s + k 1 x 1 + k 2 x 2 = s in }.
Given that the saturation level y is typically small, we assume for sake of simplicity that γy < µ 2 (s in ). We divide

the positively invariant set Λ = {ξ ∈ R 2 + | k 1 x 1 + k 2 x 2 ≤ s in } into two regions: Λ := {ξ ∈ Λ | y(ξ)
≤ y} and Λ := Λ \ Λ. We get the following system:

ẋi = (µ i (s) -u(ξ))x i , i = 1, 2. ( 12
)
with s = s in - 2 i=1 k i x i and u(ξ) =      γy(ξ) if ξ ∈ Λ γy if ξ ∈ Λ
For Λ, the existence and stability of the equilibria defined in Section 3 remain unchanged whenever they belong to this set, i.e. this concerns E 1 , and respectively E 2 and E c iff :

y(0, x * 2 ) > y, (C7) and y
(x 1 , x2 ) > y. (C8)
Note that E 0 does not exist anymore. Now in Λ, recalling the classical results of microbial competition in the chemostat [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF], we can state that the washout E w is unstable given that u = γy < µ 1 (s in ). Additionally, one equilibrium with only species x 2 can exist. Let s 2 the root of γy = µ 2 (s), with 0 < s 2 < s in .

The equilibrium E 2 = (0, x 2 ), where

x 2 = (s in -s 2 )/k 2 , exists if y(0, x 2 ) ≤ y,
or equivalently, if (C7) does not hold. It is locally stable iff:

µ 2 (s 2 ) > µ 1 (s 2 ). ( C9
)
We now want to characterize the global behavior of the system.

Lemma 1. For γy < µ 2 (s * 2 ), all trajectories of System (12) converge towards an equilibrium point.

Proof. Givent that γy < µ 2 (s * 2 ), we get:

y(0, x * 2 ) = l 2 µ 2 (s * 2 )x * 2 = µ 2 (s * 2 ) γ > y,
and , whenever E c exists:

y(x 1 , x2 ) = µ 1 (s c ) γ > min{µ 1 (s * 1 ), µ 2 (s * 2 )} γ > y, so E 2 and E c (if it exists) belong to Λ. Recalling z(ξ) = l 1 x 1 + l 2 x 2 , consider the function V (ξ) = 1 2 (z(ξ) -γ -1 ) 2 .
In Λ, we have:

V (ξ) = -γy(ξ)(z(ξ) -γ -1 ) 2 ≤ 0.
In Λ, if z(ξ) > γ -1 , we obtain:

V (ξ) = (y(ξ) -γyz(ξ))(z(ξ) -γ -1 ) ≤ y(1 -γz(ξ))(z(ξ) -γ -1 ) ≤ -γy(z(ξ) -γ -1 ) 2 ≤ 0.
Let z m = max ξ∈Λ|z(ξ)<γ -1 z(ξ). For any z c ∈ (z m ; γ -1 ), we can deduce that the

set Λ c = {ξ ∈ Λ|z(ξ) ≥ z c } is positively invariant. Let S = {ξ ∈ Λ c | V (ξ) = 0} = {ξ ∈ Λ c |z(ξ) = γ -1 }
No solution can stay identically in S, other than the equilibria E 1 , E 2 , or E c (if it exists). Applying the LaSalle theorem, every solution starting in Λ c approaches one of these equilibrium points as t → +∞. Finally, in Λ \ Λ c , there is only one equilibrium E w , which is a repeller located on the boundary. This eliminates the possibility of periodic or homoclinic orbit, so any solution starting in Λ will converge towards an equilibrium points, either E 1 , E 2 , or E c .

From extensive numerical simulations, the lemma seems to hold also for γy > µ 2 (s * 2 ). In the following, we will assume that it is true, so all trajectories of System ( 12) converge towards an equilibrium point. For Monod growth rates (Hypothesis 2), we can deduce the global behavior of System [START_REF] Frédéric | Competition between phytoplankton and bacteria: exclusion and coexistence Journal of mathematical biology[END_REF], given in Table 2. The saturation of the control has increased the robustness of the control law, see Fig. 3 in comparison with Fig. 2. In particular, reactor shutdown -the most problematic case -does not occur anymore. The analytic study of robustness for n species is obviously more delicate and is currently under investigation.

Simulation study

We illustrate our approach with numerical simulations. System (2) was simulated under Control laws (3) and [START_REF] Tao | Adaptive extremum seeking control of continuous stirred-tank bioreactors[END_REF]. The number of species n, the species characteristics θ i ∈ S and the initial conditions ξ(0) were drawn randomly. In conditions. Finally, Figure 5 illustrates the case where the additional species allows to increase the productivity of the system.

Discussion

The robustness of control strategies with respect to parametric uncertainty has been classically studied for bioprocesses. Here, we have defined the robustness to biodiversity. These two approaches are clearly not equivalent. Actually, the n species can be seen as one species with time-varying kinetic and stoechiometric parameters. Such time-varying aspect is generally not considered in parametric uncertainty. Moreover, the robustness to biodiversity also involves species competition/selection/coexistence. These phenomena, which affect bioprocess productivity, are not considered in the robustness with respect to parametric uncertainty. [START_REF] Tao | Adaptive extremum seeking control of continuous stirred-tank bioreactors[END_REF] for another set of species and initial conditions. Same legend as Fig. 4. One additional species outcompetes species x 1 and allows to improve the productivity.

Link with species competition

The classical theory of species competition [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF] does not apply when the chemostat is operated with a closed loop control, and some surprises can arise. For example, if we consider a slow growing species x 2 (with µ 2 (s) < µ 1 (s), ∀s ∈ (0, s in )), we expect that this species will rapidly be outcompeted (as it happens in open loop), so it should not affect the asymptotic behavior of the system.

Actually, with Control law (3), the substrate depletion equilibrium E 0 can become locally stable in presence of x 2 (if α 2 γs in < 1), so this species can cause reactor shutdown (see Figure 1, Case A).

Coexistence of two species

The stable coexistence of two species limited by one substrate is possible whenever their specific growth rates intersect. In open loop, this is not feasible in practice since the dilution should be chosen exactly equal to the rate at which growth curves intersect. On the other hand, [START_REF] Leenheer | Feedback control for chemostat models[END_REF] and [START_REF] Gouzé | Feedback control for nonmonotone competition models in the chemostat Nonlinear Analysis[END_REF] have proposed feedback controls for the stable coexistence of two species, based on measurement of biomass concentrations. Here, we have shown that the feedback law (3) can also lead to stable coexistence, whenever the set-point s * is chosen accordingly.

Assuming without loss of generality µ ′ 1 (s c ) > µ ′ 2 (s c ), the coexistence point E c is globally asymptotically stable iff α 1 < α 2 and 0 < s * < s c < s * 2 < s in (see Fig. 1, case D). If α 1 > α 2 , Control law (3) does not generate stable coexistence.

Challenges

Maintaining a high bioreactor productivity despite invasive species (or new individuals resulting from natural mutations) turns out to be a challenging problem to ensure process robustness. New questions arise now, such as detecting and preventing the apparition of "bad species". Observers could be set-up to early identify new competing species with reduced productivity capability. Control strategies, such as extremum-seeking [START_REF] Tao | Adaptive extremum seeking control of continuous stirred-tank bioreactors[END_REF], could be adapted in order to consider multispecies and maintain conditions favorable for the settlement of enhanced species with productivity increase. As demonstrated in this paper, such control strategies could use this opportunity to continuously improve the reactor performance.

Conclusion

We have introduced the concept of robustness to biodiversity for bioprocess control laws. We have illustrated our approach with a control law proposed in [START_REF] Mailleret | Robust Nonlinear Adaptive Control for Bioreactors with Unknown Kinetics[END_REF]. Depending on the characteristics of the additional species, some counterintuitive results may appear such as coexistence, or even reactor shutdown when a slow-growing species is introduced. The saturation of the control law allows to avoid this phenomena, and so it increases the robustness. As demonstrated in our example, this framework can be used to design robust control laws and therefore better tame biodiversity within a biotechnological process.

Figure 2 :

 2 Figure 2: Robustness of Control law (3) on a ternary plot of the morphospace S. Each vertex represents an archetype: x μ(1, 0, 0), x K (0, 1, 0), x α (0, 0, 1). The diamond represents species x 1 (0.3, 0.4, 0.3). The set of species S R , such that Control law (3) is (S R , 1, P) robust, is represented in color. The color map represents the relative increase of productivity (with respect to P 1 ). The black area corresponds to reactor shutdown (Case A).

Figure 4 ,Figure 3 :

 43 Figure 4, we show an example where Control law (3) leads to a reactor shutdown, while the trajectory converges towards the working equilibrium point with Control law[START_REF] Tao | Adaptive extremum seeking control of continuous stirred-tank bioreactors[END_REF]. We have shown previously that the saturation of the control law allows to increase the robustness for one additional species with initial conditions in the invariant set. Our simulations tends to show that the new control law is actually robust for n > 2 additional species with any initial

Figure 5 :

 5 Figure 5: Simulation of System (2) with Control law[START_REF] Tao | Adaptive extremum seeking control of continuous stirred-tank bioreactors[END_REF] for another set of species and initial conditions. Same legend as Fig.4. One additional species outcompetes species x 1 and allows to improve the productivity.

  and we

		Case A										Case D										Case E
	2.4										2.4											3
	2.2										2.2											2.8
	2										2											2.6
																						2.4
	1.8										1.8											2.2
	1.6										1.6											2
	1.4										1.4											1.8
	1.2										1.2											1.6
																						1.4
	1										1											1.2
	0.8										0.8											1
	0.6										0.6											0.8
	0.4										0.4											0.6
																						0.4
	0.2										0.2											0.2
	0										0											0
	0	1	2	3	4	5	6	7	8	9	10	0	1	2	3	4	5	6	7	8	9	10	0	1	2	3	4	5	6	7	8	9	10
	10										10											10
	9										9
	8										8											8
	7										7											7
	6										6											6
	5										5											5
	4										4											4
	3										3											3
	2										2											2
	1										1											1
	0										0											0
	0	1	2	3	4	5	6	7	8	9	10	0	1	2	3	4	5	6	7	8	9	10	0	1	2	3	4	5	6	7	8	9	10

9 Figure 1: Asymptotic behavior of System (5), Case A, D, and E (see Proposition 1 and Table

1

). Top: Specific growth rates as a function of substrate concentration. Bottom: Phase portraits with some trajectories (in purple). Filled circle: stable equilibria, open circles: unstable equilibria. Blue: species x 1 , green: species x 2 , red: coexistence. consider

Table 1 :

 1 Equilibria and robustness for Control law (3) -System (5) (see Propospecies with the best values for each trait will outcompete all the other species in any environmental conditions, and such super mutant should not exist. Actually, one may assume trade-offs between the different traits. For m traits, we consider that the m archetypes (or specialist, i.e. a species with the best value for one trait and the worst values for the others) define a Pareto front where all the species lie

	sition 1)
	Case	Conditions	E 0 E w	E 1	E 2	E c	Robustness c
	A	¬C1	LS uns	-	no	-	F
	B	C1.C4.(¬C5)	uns uns GAS uns	no	T
	C	C1.(¬C4).C5	uns uns uns GAS	no	T iff P 2 ≥ P 1
	D	C1.(¬C4).(¬C5) a uns uns uns	uns GAS T iff P c ≥ P 1
	E	C1.C4.C5 b	uns uns	LS	LS	uns T iff P 2 ≥ P 1
	LS: locally stable, GAS:globally stable, uns: unstable, -:whatever, T/F: True/False.
	a: or equivalently C1.C2.C3.C6, b: or equivalently C1.C2.C3.(¬C6), c: see Appendix

Table 2 :

 2 Equilibria and robustness for Control law (8) -System[START_REF] Frédéric | Competition between phytoplankton and bacteria: exclusion and coexistence Journal of mathematical biology[END_REF]

	Case	Conditions	E w	E 1	E 2	E 2	E c	Robustness a
	B'	C4.(¬C5).C7	uns GAS no or uns no or uns	no	T
	C'	(¬C4).C5.C7	uns uns	GAS	no	no	T iff P 2 ≥ P 1
	C"	(¬C4).C5.(¬C7)	uns uns	no	GAS	no	T iff P 2 ≥ P 1
	D'	(¬C4).(¬C5).C8	uns uns no or uns no or uns GAS T iff P c ≥ P 1
	D"	(¬C4).(¬C5).(¬C8) uns uns	no	GAS	no	T iff P 2 ≥ P 1
	E'	C4.C5.(¬C8)	uns GAS	no	uns	no	T
	E"	C4.C5.C7	uns	LS	LS	no	uns	T iff P 2 ≥ P 1
	E"'	C4.C5.(¬C7)C8	uns	LS	no	LS	uns	T iff P 2 ≥ P 1
	LS: locally stable, GAS:globally asymptotically stable, uns: unstable, T/F: True/False, a: see Appendix.

  Figure 4: Simulation of System (2) without (middle) and with (right) saturation of the control law (Eq. (3) and (8) respectively). Top: Blue thick line: species x 1 , thin lines : species x i , i > 1, red dashed line: substrate s. Bottom: productivity (blue line) and set-point (black dashed line). Left: ternary plot of the morphospace S (see legend of Fig.2). The black dots represent the additional species drawn randomly in the robust set of Control law[START_REF] Tao | Adaptive extremum seeking control of continuous stirred-tank bioreactors[END_REF].

												3.5
												3			8
												2.5
															6
												2
												1.5			4
												1
												0.5			2
												0
												0	1 000	2 000	3 000	4 000	5 000	6 000	0
															0	1 000	2 000	3 000	4 000	5 000	6 000
												0.8			0.8
												0.6			0.6
												0.4
															0.4
												0.2
															0.2
												0
												0	1 000	2 000	3 000	4 000	5 000	6 000	0	1 000	2 000	3 000	4 000	5 000	6 000
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
													8
													6
													4
													2
													0
													0	1 000	2 000	3 000	4 000	5 000	6 000	7 000	8 000
													1
													0.5
													0
													0	1 000	2 000	3 000	4 000	5 000	6 000	7 000	8 000
	0	0.1	0.2	0.3		0.4	0.5	0.6	0.7	0.8	0.9	1

Appendix

We give the conditions on parameters θ 2 that should be fulfilled to avoid any productivity loss at alternative steady states (see Tables 1 and2). For E 2 , this yields:

, we have:

For P 2 , using s 2 = K2γy μ2-γy , we finally get: