Arshia Cont
email: cont@ircam.fr

REALTIME AUDIO TO SCORE ALIGNMENT FOR POLYPHONIC MUSIC INSTRUMENTS USING SPARSE NON-NEGATIVE CONSTRAINTS AND HIERARCHICAL HMMS

We present a new method for realtime alignment of audio to score for polyphonic music signals. In this paper, we will be focusing mostly on the multiple-pitch observation algorithm proposed based on realtime Non-negative Matrix Factorization with sparseness constraints and hierarchical hidden Markov models for sequential modeling using particle filtering for decoding. The proposed algorithm has the advantage of having an explicit instrument model for pitch obtained through unsupervised learning as well as access to single note contribution probabilities which construct a complex chord instead of modeling the chord as one event.

INTRODUCTION

Traditionally, score following (or realtime audio to score alignment) was introduced for two main reasons : First for musical accompaniment of human performers with a computer and second, for synchronization of live sound processing algorithms for instrumental electroacoustic composition. Today, audio to score alignment has found applications in wider domains such as audio indexing and pedagogical systems.

Among different systems available the only polyphonic score follower systems that we know of are reported in [START_REF] Raphael | A Hybrid Graphical Model for Aligning Polyphonic Audio with Musical Scores[END_REF] and [START_REF] Dannenberg | Polyphonic Audio Matching for Score Following and Intelligent Audio Editors[END_REF] by Raphael and Dannenberg respectively. Raphael's system in special, has been evaluated and demonstrated for realtime polyphonic performances on the Piano. In both works, the spectral model for a chord is assumed to be fixed and there is no explicit notion of instrument model. In the proposed system, the instrument model is learned explicitly as described in Sections 3.2. Moreover, we have access to individual pitch contribution to a complex chord during the observation process whereas in both systems the model considers a chord as one fixed object when calculating probabilities. For a detailed review of other systems, we refer the curious reader to [START_REF] Cont | Training ircam's score follower[END_REF].

In Section 2 we describe the general architecture of the proposed polyphonic score follower. Due to limited space, we focus mainly on the multiple pitch observation system based on non-negative matrix decomposition (NMD) with sparseness constraint in Section 3. The sequential modeling will be briefly described in Section 4. In the end, we provide results and further discussions on the subject.

GENERAL ARCHITECTURE

As a basic description, a score following system takes audio chunks in realtime as input and provides an alignment to the music score in the output. In our system, audio input is represented by This research was supported by the EU project I-MAESTRO (see http://www.i-maestro.org/). feature vectors y k for each realtime frame k. This vector is used to compute the observation likelihood P(y k |x k ) of being on the event x k in the score. The score is defined as a hidden Markov chain of states representing sequential events as described in Section 4 and providing transition priors P(x k-1 |x k ). Current belief of the system is calculated using the observation likelihood, transition priors and previous belief of the system P(x k-1 |y 1:k-1 ).

The system proposed in this paper is recursive in nature. This means that the observation (or acoustic) model and the sequential model are tied together and the way they are applied depends on where we are currently in the score through two different decoding schemes. This characteristic is the main feature of our sequential modeling which is based on Hierarchical Hidden Markov Models [START_REF] Fine | The hierarchical hidden markov model : Analysis and applications[END_REF] described in details in Section 4. Generally speaking, it consists of two main hierarchies, the higher modeling the atoms in the score (notes, chords and silences) and the lower modeling time. A vertical transition switches from the first hierarchy to the second by the introduction of a new event by the performer and a strike or end transition sends the system back to the first hierarchy as shown in Figure 1.

First Hierarchy

NMF Pitch, Attack, Rest and Vertical Transition.

Second Hierarchy

Sustain/Rest and Strike Observations. With respect to this structure there are two separate decoding and belief update processes for each hierarchy level. If the current position is on the first hierarchy, the decoding is done using Sequential Importance Sampling with Resampling (SISR) with a fixed number of particles and using residual resampling [START_REF] Liu | Sequential Monte Carlo methods for dynamic systems[END_REF]. If we are in the second hierarchy, the decoding is done using the regular forward (α) variable that contains the transition probabilities and observations of the second level.

Strike or End Vertical Transition

The observation process can be regarded as two independent systems : one providing attack, sustain and rest probabilities (showing the attribute of the ongoing music signal according to the three categories) and another observation process for polyphonic pitch. The attack, sustain and rest model has been previously discussed along with its training algorithm in [START_REF] Cont | Training ircam's score follower[END_REF] and will not be the subject of this paper.

MULTIPLE PITCH OBSERVATION MODEL

Multiple pitch analysis is currently an important topic in music and speech signal processing among others. It is mostly applied for multiplef 0 recognition and music transcription. However, most of the multiplef 0 analysis systems available are either computationally very heavy or unadaptable for realtime applications. For a detailed overview of different systems and approaches we refer the reader to [START_REF] De Cheveigné | Multiple f0 estimation[END_REF]. In this section we introduce a novel algorithm for multiple pitch analysis of music signals. The proposed algorithm is based on non-negative matrix factorization (NMF) [START_REF] Lee | Algorithms for nonnegative matrix factorization[END_REF]. NMF was first used to find parts representation of image data. In general terms, given a matrix V , the original NMF algorithm is an unsupervised algorithm that learns desired number of templates stored as columns in W and a decomposition matrix H through multiplicative updates which results in V ≈ W H. Non-negativity constraint in the algorithm implies additivity of parts in the representation in order to construct V .

In our application, W is learned offline and contains spectral templates for all pitches of a given instrument. During live performance W is kept fixed and the algorithm tries to learn a decomposition of the ongoing audio spectrum using the fixed templates (i.e. it only learns H). A similar approach has been recently reported by Sha and Saul [START_REF] Sha | Real-time pitch determination of one or more voices by nonnegative matrix factorization[END_REF] where they use regular NMF for both learning and realtime decomposition. Our approach is strictly different in the performance stage with use of Non-negative Matrix Decomposition with sparseness constraints. Later in this section we compare the output of the two and show that the proposed algorithm undergoes a better representation for polyphonic pitch analysis. Sparse nonnegative constraint has been previously introduced by Hoyer [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF]. Our approach is different from Hoyer's with use of different sparse measures adapted for multiple pitch observation and further refinements for realtime decomposition which is not the case in [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF].

We start this section by introducing the signal processing frontend used for the representation of the signal (V ). We continue by introducing the training phase of the pitch observation (learning W ), demonstrating the sparse non-negative matrix decompotision algorithm for realtime performance and how the likelihood probability is calculated from the observation.

Signal processing front-end

The signal processing front end used for this observation is the result of a fixed point analysis of frequency to instantaneous frequency mapping of the ongoing audio spectrum [START_REF] Hideki Kawahara | Fixed point analysis of frequency to instantaneous frequency mapping for accurate estimation of f0 and periodicity[END_REF]. The short-time Fourier transform (STFT) already in use in the observation model is an efficient tool for instantaneous frequency (IF) estimation [START_REF] Abe | Harmonic tracking and pitch extraction based on instantaneous frequency[END_REF]. Given z(ω,t) as the analytical form of the STFT, the mapping

λ (ω,t) = ∂ ∂t arg[z(ω,t)]
can be computed efficiently and in real-time using STFTs [START_REF] Abe | Harmonic tracking and pitch extraction based on instantaneous frequency[END_REF] and the fixed points of this mapping can be extracted using the following criteria [START_REF] Hideki Kawahara | Fixed point analysis of frequency to instantaneous frequency mapping for accurate estimation of f0 and periodicity[END_REF] :

λ (ω * ,t) = ω * and ∂ λ ∂ ω | ω=ω * < 1.
As a result, vector V would be non-negative amplitudes of the fixed-point instantaneous frequency. The same representation has been reported in [START_REF] Sha | Real-time pitch determination of one or more voices by nonnegative matrix factorization[END_REF].

Learning pitch templates

Once V is prepared in realtime, it will be decomposed using pre-learned pitch templates stored in W . The idea here is that the pitch observation is familiar with all pitches that can be produced by an instrument. Moreover, we are interested in instrument specific pitch models to be learned explicitly by the system. For example, for following a piano score, matrix W would contain all 88 pitches as 88 different columns. Learning pitch templates is done off-line and only once for each instrument using databases of instrumental sounds [START_REF] Ballet | Studio online 3.0 : An internet "killer application" for remote access to ircam sounds and processing tools[END_REF]. For each audio file in the database (corresponding to a known pitch), training is an iterative NMF with symmetric kullbackleibler divergence as error measure. In addition, since each audio spectrum contains additional information such as noise besides a harmonic spectrum, we put two constraints on the learned W . First, we decompose V into two vectors (W has two columns) where we only learn one vector and have the other fixed as white noise. This criteria helps the algorithm focus more on the harmonic structure of V . Furthermore, we constrain this harmonic structure in each learning iteration by an envelope. This envelope is constructed from the pitch information of the audio file (taken from the name of the file in the database) and emphasizes frequencies around the fundamental with a decreasing envelope towards the end and close to zero for frequencies less than the fundamental. This constraint improves common octave and harmonic errors.

Equation 1 shows the NMF updates with kullback-leibler divergence for each iteration as described above and adapted from [START_REF] Lee | Algorithms for nonnegative matrix factorization[END_REF].

Here Env is the envelope constraint for each pitch, ⊗ is an element by element multiplication and V is the short-time fixed-point instantaneous frequency representation of the audio file described before.

H aµ ←-Env ⊗ H aµ ∑ i W ia V iµ /(W H) iµ ∑ k W ka W ia ←-Env ⊗W ia ∑ i H aµ V iµ /(W H) iµ ∑ ν H aν (1)
When the training reaches an acceptable stopping criteria, the harmonic spectra in the local W will be saved in the global W and the algorithm continues to the next audio file in the database until it constructs W for all pitches of the instrument.

Sparse non-negative Decomposition for realtime application

The concept of sparse coding refers to a representational scheme where only a few units out of a large population are effectively used to represent typical data vectors. One of the useful properties of NMF is that it usually produces a sparse representation of the data. However this sparseness is more of a side-effect than a goal and one can not control the degree to which the representation is sparse.

Numerous sparseness measures have been proposed and used in the literature. In general, these measures are mappings from R n to R which quantify how much energy of a vector is packed into a few components. For our application there are two sparseness measures working in parallel : an approximation of ℓ ε norm by the tanh function as g(x) = tanh(|ax| b ), where a and b are positive constants with b greater than 1, and an ℓ 2 norm which is crucial for the normalization of the results. This measure is demonstrated in Equation 2 with N being the dimension of the space.

sparseness(x) = √ N -∑tanh(|xi| 2 )/ ∑ x 2 i √ N -1 (2)
In order to obtain H while V (signal processing representation) and W (pitch templates) are known, we use a gradient descent update, H = H -µ H W T (W H -V ), instead of the original NMF multiplicative update and project each vector in realtime to be non-negative and have desired ℓ 2 and ℓ ε norms. This is achieved by first projecting the new H to the ℓ ε space and solving a quadratic equation where the projection has the desired ℓ 2 norm for each iteration. Note that this algorithm can be adapted for other applications.

For the multiple pitch observation model, ℓ 2 measure will be provided in realtime to the algorithm from the log-energy of the spectrum. The ℓ ε constraint represents the sparsity of the output, the bigger the more sparse and vice versa. For this reason, depending on the local activity of future events in the score ℓ ε can be adjusted in realtime.

In order to understand the importance of sparseness constraints, Figure 2 demonstrates three piano chords (shown as the first three chords in Figure 5(a)) using NMD without sparseness constraints and the main portion of Figure 3 shows the same analysis but with the proposed algorithm. Each value on the y-axis corresponds to a specific pitch in the W matrix (in this case Piano) sorted according to their frequency. From the figure, it can be noted that the NMD with sparseness constraint has achieved a more compact representation than the one without. While the activity in the regular NMD spreads over all frequency regions, it can be seen visually that for the sparse representation it is more active in the frequency region that the pitches occur. Regular NMD has a natural tendency to reconstruct a pitch with its octave, dominant and subharmonic whereas NMD with sparseness constraint, due to the imposed constraints tends to enforce main pitches for each chord. Moreover, because of the use of energy as ℓ 2 constraint, we can visually distinguish between three chords and at the onset we can see high activity in the noise template which corresponds to transition states.

Observation probability

The proposed system above can be used independent of the score follower for multiple pitch analysis of polyphonic instruments with further consideration. For a score following application, we are interested in the probability of each chord in the score played by the musician which will serve as observation likelihoods P(y k |x k ) during decoding. Moreover, this probability should consider human errors playing the score (for example playing one or more notes off the written pitch in the score). For this purpose, using the given score, we construct Gaussian mixtures over the pitch templates for each chord. For each note in a chord, a Gaussian will be centered on the index of that note in the overall pitch template. Each Gaussian corresponds to a note and moreover, scores possible human errors if adjacent pitches are played. The NMD decomposition at each frame will be filtered by these gaussian mixtures and therefore the normalized energy in each filtered vector correctly represents the probability of each chord at each frame.

Figure 3 represents the steps mentioned above for the NMD results from previous section and for the first three chords in Figure 5(a). The right portion of Figure 3 demonstrates the Gaussian mixtures for the first three chord and the figure on the bottom shows the chord probabilities extracted from filtering the NMD for each frame. Note that the three chords in analysis in this example have at least two common pitches among each other and hence the shown example is not an ideal analysis.

B0 Db1 Eb1 F1 G1 A1 B1 Db2 Eb2 F2 G2 A2 B2 Db3 Eb3 F3 G3 Bb3 C4 D4 E4 Gb4 Ab4 Bb4 C5 D5 E5 Gb5 Ab5 Bb5 C6 D6 E6 Gb6 Ab6 Bb6 C7 D7 E7 G7 A7 B7 noise 0 0.1 0.2 noise B7 A7 G7 E7 D7 C7 Bb6 Ab6 Gb6 E6 D6 C6 Bb5 Ab5 Gb5 E5 D5 C5 Bb4 Ab4 Gb4 E4 D4 C4 Bb3 G3 F3 Eb3 Db3 B2 A2 G2 F2 Eb2 Db2 B1 A1 G1 F1 Eb1

SEQUENTIAL MODELING

Once the observations are gathered for an ongoing audio frame, they should be used for alignment to a sequential model of the music score. In our model, observations are calculated for pertinent states and serve as likelihood probabilities in a Bayesian framework. The prior information are transition probabilities of the sequential model of the given score. The sequential model used in our system is inspired by Hierarchical Hidden Markov Models [START_REF] Fine | The hierarchical hidden markov model : Analysis and applications[END_REF] which we will call Score Hidden Markov Model (SHMM) throughout this paper.

The SHMM consists of two different hierarchies. The first level consists of high-level musical states that correspond to atoms in the score, that is notes, chords and rests and the second level contains temporal modeling according to the time-value of each event in the score. Since we are not dealing with music transcription we model each chord as one state. The observations pertinent to the first level are attacks, multiple pitch and rest observations. Moreover, the SHMM model described here allows us to have different activations and entries for each high-level state according to where we are in the sequence. Figure 4 shows the note model used in SHMM. The firstlevel in the figure corresponds to a note or chord event in the score and its corresponding observations are attack and multiple pitch. Once this event is detected, SHMM enters automatically into the second hierarchy, the temporal states. For a note/chord state, these states correspond to sustain observation. The transition probabilities and the number of states are set using the time-value of the event according to the score and a temporal modeling scheme. The time modeling used here is the same as [START_REF] Cont | Training ircam's score follower[END_REF]. There is an additional end state in the second hierarchy which automatically sets index to the last first level state observed. Other than passing through the time model, there are strike observations whose function is to halt the temporal model if there is a sudden rest or attack. The rest model is essentially the same but with different observations and only one strike state for attacks.

... The final SHMM of a given score, would be a chain of note/chord and rest models as described above according to the music score.

RESULTS AND DISCUSSION

Figure 5 shows results for a simulation of the realtime audio to score alignment with the score shown in Figure 5(a) as an excerpt of a contemporary music piece ("Pluton" by Philippe Manoury) with the alignment results in Figure 5 In this paper, we concentrated mostly on the multiple pitch feature of the proposed system. Note that this system is independent of the proposed application and can be used freely in other contexts.

Despite learning the templates from a database, the pitch observation has proven to be robust for various performance situations and the observation probabilities maintain the same relative information. Other aspects of this system including particle filtering and results for other instruments were not discussed and will be reported in a separate publication. The realtime version of this system is developed for the Max/MSP realtime graphical programming environment and is currently planned for production of various contemporary music pieces.

Fig. 1 .

 1 Fig. 1. General diagram of score follower's realtime message passing with appropriate observations for each hierarchy.

Fig. 2 .

 2 Fig. 2. Chord analysis for the first three chords of Figure 5(a) using regular NMF without sparseness constraints.

Fig. 3 .

 3 Fig. 3. Chord analysis for the first three chords of Figure 5(a). Figure on the right shows the Gaussian mixtures constructed from the score on the pitch templates. The figure on the bottom shows the probability of chords from the normalized energy of each Gaussian mixture structure multiplied with the raw pitch observation for each frame.

Fig. 4 .

 4 Fig. 4. Score Hidden Markov Model -one note/chord model

  Figure 5 shows results for a simulation of the realtime audio to score alignment with the score shown in Figure 5(a) as an excerpt of a contemporary music piece ("Pluton" by Philippe Manoury) with the alignment results in Figure 5(b) and (c) where Figure 5(c) shows a zoomed view for the last music bar of the score in Figure 5(a). Each index in Figure 5(b) and (c) shows the event (note/chord) number in the order they appear in the score.

  Zoom view of last bar.

Fig. 5 .

 5 Fig. 5. Example of placing a figure with experimental results.

ACKNOWLEDGMENT

The research was supported in part by the European Commission under Contract IST-026883 I-MAESTRO. The author would like to acknowledge the EC IST FP6 for the partial funding of the I-MAESTRO project (http://www.i-maestro.net/).


Jan Sieber
email: j.sieber@exeter.ac.uk

Alain Rapaport
email: rapaport@supagro.inra.fr

Serafim Rodrigues
email: serafim.rodrigues@plymouth.ac.uk

Mathieu Desroches
email: mathieu.desroches@inria.fr

† A Rapaport

§ M Desroches

A new method for the reconstruction of unknown non-monotonic growth functions in the chemostat

We propose an adaptive control law that allows one to identify unstable steady states of the open-loop system in the single-species chemostat model without the knowledge of the growth function. We then show how to use a continuation method to reconstruct the whole graph of the growth function. Two variants, in continuous and discrete time, are presented. The case of two species in competition is also examined.

Introduction

We recall the classical chemostat model [START_REF] Smith | The Theory of the Chemostat[END_REF] for a single species (biomass b) consuming a substrate (mass s):

ṡ = -µ(s)b + D(s in -s) ḃ = µ(s)b -Db (1)
where the dilution rate D (the input ) is the manipulated variable, which takes values in a bounded positive interval [D min , D max ], and µ(•) is a non-negative Lipschitz continuous function with µ(0) = 0 .

We consider here the following scenario: the function µ(•) is unknown and possibly non-monotonic. Our objective is to reconstruct the graph of the function µ(•) on the domain (0, s in ) by varying the input D in time. On-line measurements are only available for the variable s (that is, s is the output ). This setup is realistic for experimental investigations such as in [START_REF] Antonelli | Set point regulation of an anaerobic digestion process with bounded output feedback[END_REF], however, demonstrations in this paper are based entirely on simulations of model system [START_REF] Abed | Stabilization of period doubling bifurcations and implicatons for control of chaos[END_REF].

Remark: We assume that the yield coefficient of the bio-conversion is known. Then we can choose it equal to unity without loss of generality such that µ(s)b appears with the same pre-factor 1 (once positive, and once negative) in both equations of (1)).

The problem of kinetics estimation in biological and biochemical models has been widely addressed in the literature ( [START_REF] Dochain | Adaptive identification and control algorithms for non linear bacterial growth systems[END_REF][START_REF] Dochain | On-line estimation of microbial growth rates[END_REF][START_REF] Dochain | On-line estimation of specific growth rates : an illustrative case study[END_REF][START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF][START_REF] Dochain | Dynamical Modelling, Analysis, Monitoring and Control Design for Nonlinear Bioprocesses[END_REF][START_REF] Perrier | Tuning of observer-based estimators: theory and application to the on-line estimation of kinetic parameters Control Engineering Practice[END_REF][START_REF]Dochain State and parameter estimation in chemical and biochemical processes : a survey[END_REF]), either as a parameter estimation problem (one chooses a priori analytical expression of the function µ(•)), or as an on-line estimation of the kinetics (one aims at determining µ(s(t)) at the current time t). Here, our objective is different: we aim at reconstructing the whole graph of the function µ(•). To our knowledge, this precise problem has not been addressed in the literature.

When the growth function is monotonic, a common way to reconstruct points on the graph of the growth function µ(•) is to design a series of experiments fixing the dilution rate D with different values and wait until the system settles to a steady state (s ⋆ , b ⋆ ) [START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF]. As long as D is less than µ(s in ), it is well known that the dynamics converges to a unique positive equilibrium that satisfies µ(s ⋆ ) = D (see for instance [START_REF] Smith | The Theory of the Chemostat[END_REF]). An alternative approach is to fix a value of s, say s and design an adaptive control law D(•) that stabilizes the system about the steady state (s, s in -s), with the value of D converging to µ(s). This technique has been proposed in [START_REF] Antonelli | Set point regulation of an anaerobic digestion process with bounded output feedback[END_REF] to stabilize such dynamics without the knowledge of µ(•) and under the constraint

D ∈ [D min , D max ].
Unfortunately, these two reconstruction strategies cannot be used to identify non-monotonic growth functions, such as the Haldane law. More precisely, these two techniques require the steady state to be stable in open loop, and consequently cannot reconstruct any part of the graph of a function µ(•) where µ is non-increasing (such as the example, shown schematically in Figure 1. Furthermore, the global convergence of these methods is not satisfied in case of bi-stability, which is present in [START_REF] Abed | Stabilization of period doubling bifurcations and implicatons for control of chaos[END_REF], with non-monotonic growth functions µ (see again [START_REF] Smith | The Theory of the Chemostat[END_REF]). s µ stable unstable In this work, we first propose that it can be useful to introduce a feedback control loop into (1) to identify the growth function µ of the open-loop system (1) (that is, (1) with constant input D). The feedback control law is initially a simple constrained piecewise linear relation between system input D and system output s:

D(s, D, s) = sat [Dmin,Dmax] D -G 1 (s -s) , (2)
where D and s are reference values, and G 1 > 0 is the linear control gain. To ensure realistic values for the input D, we enclose the feedback rule into the saturation function

sat [Dmin,Dmax] (x) =      D max if x > D max , x if x ∈ [D min , D max ], D min if x < D min ,
where the limits D min and D max are the extreme dilution rates that can be achieved experimentally. In sections 3 and 5 we will then explore adaptation rules for the reference values (s, D) which ensure that asymptotically for t → ∞ the input satisfies

D(s, D, s) = D, (3)
or, equivalently, the output satisfies

s = s, (4)
such that the control is asymptotically non-invasive. The result is a new adaptive control law that stabilizes the dynamics about any desired equilibrium point without requiring a-priori knowledge of its location, and whatever is the monotonicity of the growth function. One requirement on the adaptive law is that it should work uniformly well around a local maximum of µ (non-invasive feedback laws such as filtered feedback [START_REF] Abed | Stabilization of period doubling bifurcations and implicatons for control of chaos[END_REF] or time-delayed feedback [START_REF] Pyragas | Continuous control of chaos by self-controlling feedback[END_REF] do not achieve this).

We propose two variants of our method, making adaptation in continuous time (Section 3) or in discrete time (Section 5).

Then we show that a simple continuation method allows us to the reconstruct continuously the graph of the growth function, even in the case of non-monotonic growth functions. Finally we investigate the case of two species that compete for the same common substrate.

Global stability of the simple feedback law

Our starting point is that the feedback law (2) is, within reasonable limits, globally stabilizing. Suppose that we choose the reference value s from an interval [s min , s in ) ⊂ (0, s in ), and that the limits on the input cover the growth function µ on this interval:

D min < µ(s) for all s ∈ [s min , s in ],
(5)

D max > µ(s) for all s ∈ [0, s in ]. (6)
These conditions mean that the graph of µ does not cross the thick parts of the horizontal lines D min and D max bounding the grey area in Figure 2 from below and above.

Proposition 1. Suppose that the reference values (s, D) are chosen from [s min , s in ) × [D min , D max ], that the growth function µ satisfies ( 5)- [START_REF] Boudjellaba | Stability Loss Delay in Harvesting Competing Populations[END_REF], and that the gain G 1 is chosen sufficiently large, that is,

G 1 > -min s∈[0,sin] µ ′ (s), and (7)
G 1 > µ(s in ) - D s in - s . (8)
Then the controlled system (1) with D = D(s, D, s), given in (2), has a stable equilibrium (s eq , b eq ) ∈ [0, s in ) × (0, ∞), which attracts all initial conditions (s(0), b(0)) ∈ [0, s in ) × (0, ∞).

Proof of Proposition 1. If D > 0, and the growth function µ satisfies µ(0) = 0 and, for s > 0, µ(s) > 0 then the set

R = {(s, b) : s ∈ [0, s in ), b > 0}
is positively invariant (that is, trajectories starting in R will stay in R for all positive times). Furthermore, all trajectories starting in R approach the subspace

T = {(s, b) ∈ R : s + b = s in }
with rate at least D min forward in time. This implies that it is sufficient to check if all trajectories in T converge to a unique equilibrium. On T the equation of motion can be expressed as a differential equation for s only:

ṡ = [D(s, D, s) -µ(s)][s in -s]. (9)
First, let us check that the equilibrium at s = s in is unstable. The term

-G 1 (s in -s) is negative such that D -G 1 (s in -s) < D max for all admissible D. Assumption (8) guarantees that D -G 1 (s in -s) < µ(s in ).
Assumption [START_REF] Benoît | Dynamic bifurcations: proceedings of a conference held in Luminy[END_REF] guarantees that also D min < µ(s in ). Hence,

D(s in , D, s) -µ(s in ) < 0
for all admissible (s, D). Thus, the prefactor of s in -s in ( 9) is negative such that the equilibrium at s in is unstable for all admissible (s, D). Since ṡ > 0 at s = 0, there must be other equilibria of ( 9) in (0, s in ), which are given as solutions s eq of D(s eq , D, s) = µ(s eq ). Now let us check indirectly that none of the equilibria can satisfy D min = µ(s eq ).

Assume that ( 9) had an equilibrium s eq with D min = µ(s eq ). Then s eq has to be less than s min due to assumption [START_REF] Benoît | Dynamic bifurcations: proceedings of a conference held in Luminy[END_REF]. However, if s eq < s min , then D -G 1 (s eq -s) > D ≥ D min for all admissible (s, D). Hence D(s eq , D, s) > D min (recall that D min = µ(s eq ) by assumption of the indirect proof) such that D(s eq , D, s) -µ(s eq ) > 0, which means that s eq cannot be equilibrium.

Assumption [START_REF] Boudjellaba | Stability Loss Delay in Harvesting Competing Populations[END_REF] excludes that equilibria with µ(s eq ) = D max exist, hence all remaining equilibria s eq ∈ (0, s in ) must satisfy D -G 1 (s eq -s) = µ(s eq ).

Condition [START_REF] Campillo | Effect of population size in a Prey-Predator model[END_REF] ensures that this equation has a unique solution and that this solution corresponds to a stable equilibrium (which must be in (0, s in ) because the boundaries of (0, s in ) are inflowing for ( 9)).

Proposition 1 ensures that the output s eq of the controlled system (1) with ( 2), after transients have decayed, is a well-defined smooth function of the parameters (s, D) as long as (s, D) are chosen from (s min , s max ) × (D min , D max ). We express this fact by using the bracket notation:

s eq (s, D) = lim t→∞ s(t)
where s is output of ( 1), ( 2).

The function s eq can be evaluated at any admissible point by setting the parameters (s, D) in the definition (2) of the feedback rule, waiting until the transients of (1) have settled, and then reading off the output s. Equilibria of the uncontrolled system can then, according to (4), be found as roots of s eq (s, D) -s. More specifically, we know that, for any admissible s, D = µ(s) if and only if s eq (s, D) = s. [START_REF] Dochain | On-line estimation of microbial growth rates[END_REF] Relation [START_REF] Dochain | On-line estimation of microbial growth rates[END_REF] permits us to identify µ(s) as the unique root of s eq (s, •) -s. Sections 3-5 will explore two strategies to find this root for a range of admissible s efficiently.

An adaptive control scheme

The first strategy is a dynamic feedback that comes on top of the feedback law (2) for D. We treat D not as a parameter but introduce an additional dynamical equation for D, achieving local convergence of the output s to any reference value s ∈ (0, s in ) without the knowledge of the growth function µ. The asymptotic value of D allows then one to reconstruct the value µ(s).

Proposition 2. Fix a number s ∈ (0, s in ) and take numbers D min , D max that fulfill 0 < D min < µ(s) < D max . Then the dynamical feedback law Locally about s = s, the closed loop system is equivalent to the three-dimensional dynamical system

D(s, D) = sat [Dmin,Dmax] D -G 1 (s -s) Ḋ = -G 2 (s -s)( D -D min )(D max -D) (11
   ṡ = -µ(s)b + ( D -G 1 (s -s))(s in -s) ḃ = µ(s)b -( D -G 1 (s -s))b Ḋ = -G 2 (s -s)( D -D min )(D max -D)
This system admits the unique positive equilibrium E ⋆ = (s, s in -s, µ(s)).

For simplicity, we write the dynamics in the variables (z, s, D) coordinates, where z is defined as z = s + b:

   ż = ( D -G 1 (s -s))(s in -z) ṡ = -µ(s)(z -s) + ( D -G 1 (s -s))(s in -s) Ḋ = -G 2 (s -s)( D -D min )(D max -D)
The Jacobian matrix at E ⋆ in these coordinates is

  -µ(s) 0 0 -µ(s) -(µ ′ (s) + G 1 )(s in -s) s in - s 0 -G 2 (µ(s) -D min )(D max -µ(s)) 0  
Its eigenvalues are λ 1 = -µ(s) < 0 and λ 2 , λ 3 as eigenvalues of the sub-matrix

M = -(µ ′ (s) + G 1 )(s in -s) s in - s -G 2 (µ(s) -D min )(D max -µ(s)) 0 Then, one has det(M ) = G 2 (µ(s) -D min )(D max -µ(s))(s in -s) tr(M ) = -(µ ′ (s) + G 1 )(s in -s)
and concludes about the exponential stability of E ⋆ when G 2 > 0 and G 1 > -µ ′ (s). Finally, one obtains from (11) that D or D converges toward the unknown value µ(s).

Reconstruction of the entire growth function

Now, we use a dynamic continuation method, letting s change slowly with time as solution of the simple dynamics ṡ = εs(s in -s) [START_REF] Dochain | Dynamical Modelling, Analysis, Monitoring and Control Design for Nonlinear Bioprocesses[END_REF] to explore the right part of the graph of µ(•) when ε is a small non-negative number, and to explore the left one when ε is a small non-positive number. Recall that continuation denotes a numerical method that belongs to the family of homotopy methods and is used to compute solution branches of parametrized nonlinear equations via a predictor-corrector strategy [START_REF] Allgower | Introduction to Numerical Continuation Methods[END_REF]. For the continuation, the gain G 1 has to be been chosen uniformly large according to (7).

A discrete-time adaptation of the method

In this section, we propose an alternative to the continuous adaptation of D and s: we treat the root problem 0 = s eq (s, D) -s with ordinary numerical root-finders such as the Newton iteration. We present here an approach that combines the two steps of the method (the adaptive control and the continuation) in a discrete time framework.

In an experimental setting one will have to adapt the numerical methods to the lower accuracy of experimental outputs (see [START_REF] Sieber | Experimental continuation of periodic orbits through a fold[END_REF] for a demonstration in a mechanical experiment) but for this paper we restrict ourselves to a numerical demonstration. In the single-species case, one profits from the knowledge of an approximate derivative of s eq with respect to D to make the Newton iteration more efficient. Suppose, we plan to identify the growth function µ in a sequence of points sk = s0 + kδ (where δ > 0 is small). The function values µ(s k ) are the roots Dk of s eq (s k , •) -sk . Then we can approximate

∂s eq ∂ D ( D, sk ) = 1 G 1 + µ ′ (s eq ( D, sk )) ≈ 1 G 1 + D-Dk-1 sk -s k-1
to obtain the iteration Dnew = Ds eq ( D, sk )

G 1 + D-Dk-1 sk -s k-1 , (13)
starting from D = Dk-1 , or (for k > 2)

D = Dk-1 + Dk-1 -Dk-2 sk-1 -sk-2 .
For the initial step (k = 1) the derivative of s eq has to be either guessed or approximated with a finite difference (we used the latter in our numerical simulations). Note that at no point it is necessary to set the internal states s or b of system (1). Only the reference values (s, D) have to be set.

A simplified discrete scheme

The scheme (13) permits one to find µ(s k ) for an a-priori prescribed set of admissible abscissae sk . If one wants to recover only the graph of µ one does not need to prescribe the sequence sk a-priori, thus, avoiding a Newton iteration. Suppose that we know already two points p k-1 = (s k-1 , D k-1 ) and p k = (s k , D k ) on the curve (s, µ(s)). Then we set

(s new,k+1 , Dnew,k+1 ) = p k + δ p k -p k-1 p k -p k-1 , (14)
where δ > 0 is the approximate desired distance between points along the curve (s, µ(s)), and run the controlled experiment with the reference values (s, D) = (s new,k+1 , Dnew,k+1 ) in (2) until the transients have settled to obtain the next point on the curve s k+1 = s eq (s new,k+1 , Dnew,k+1 )

D k+1 = D(s k+1 , Dnew,k+1 , snew,k+1 ) = Dnew,k+1 -G 1 (s k+1 -snew,k+1 ) (15)
This simplified procedure cannot guarantee the identification of µ at prescribed equidistantly spaced values of s but finds µ(s k ) for a (nearly evenly spaced) sequence s k given by the intersections of the lines D = D new,k -G 1 (s -snew,k ) with the graph D = µ(s).

The two species case

The extension of the chemostat model ( 1) considers two species that compete for the same substrate. The two-species model can be written as follows

     ṡ = - 2 i=1 µ i (s)b i + D(s in -s) ḃi = µ i (s)b i -Db i (i = 1, 2) (16)
The two-species model has co-existing equilibria E * i , which correspond to the state where species i is present but the other species is not. The following result shows that feedback stabilization based on input D and output s breaks down for the equilibrium corresponding to the species with the smaller growth rate. Proposition 3. Fix s ∈ (0, s in ) and consider the equilibrium E ⋆ 2 = (s, 0, s in -s).

1. If µ 1 (s) > µ 2 (s), there does not exist feedback D(•) of the form

D = f (s, ξ), ξ = g(s, ξ), (ξ ∈ R k )
with f (s, 0) = µ 2 (s) and g(s, 0) = 0, that stabilizes asymptotically the system ( 16) about E ⋆ 2 .

2. If µ 1 (s) < µ 2 (s), then the feedback (11) exponentially stabilizes the system [START_REF] Perrier | Tuning of observer-based estimators: theory and application to the on-line estimation of kinetic parameters Control Engineering Practice[END_REF] locally about E ⋆ 2 , for any positive constants (G 1 , G 2 ) such that G 1 > -µ ′ 2 (s). Furthermore one has lim t→+∞

D(t) = µ 2 (s)
Proof of Proposition 3. Consider the extended dynamics

         ṡ = - 2 i=1 µ i (s)b i + f (s, ξ)(s in -s) ḃi = µ i (s)b i -f (s, ξ)b i (i = 1, 2) ξ = g(s, ξ) that we write in (z, b 1 , b 2 , ξ) coordinates with z = s + b 1 + b 2 :    ż = f (z -b 1 -b 2 , ξ)(s in -z) ḃi = (µ i (z -b 1 -b 2 ) -f (z -b 1 -b 2 , ξ))b i ξ = g(z -b 1 -b 2 , ξ) 1. At equilibrium E ⋆ 2 , the Jacobian matrix possesses the following form in (z, b 1 , b 2 , ξ) coordinates     -µ 2 (s) 0 0 0 0 µ 1 (s) -µ 2 (s) 0 0 ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆     ,
which has the positive eigenvalue µ 1 (s) -µ 2 (s). This proves that E * 2 is unstable whatever the choice of the feedback D(•).

At equilibrium E ⋆

2 , the Jacobian matrix can be written as follows, in

(z, b 1 , b 2 , D) coordinates    -µ 2 (s) 0 0 µ 1 (s) -µ 2 (s) 0 ⋆ M    with M = -(µ ′ 2 (s) + G 1 )(s in -s) s in - s -G 2 (µ 2 (s) -D min )(D max -µ 2 (s)) 0 Its eigenvalues are λ 1 = -µ 2 (s) < 0, λ 2 = µ 1 (s) -µ 2 (s) < 0, λ 3 and λ 4 with λ 3 .λ 4 = G 2 (µ 2 (s) -D min )(D max -µ 2 (s))(s in -s) λ 3 + λ 4 = -(µ ′ 2 (s) + G 1 )(s in -s)
As in the proof of Proposition 2, one concludes the exponential stability of E ⋆ 2 when G 2 > 0 and G 1 > -µ ′ 2 (s), and the convergence of D(•) toward µ 2 (s).

Consequently, the adaptive control scheme proposed in Section 3 only allows one to reconstruct the larger of the two growth rates at any given s. Nevertheless, the dynamics of [START_REF] Dochain | Dynamical Modelling, Analysis, Monitoring and Control Design for Nonlinear Bioprocesses[END_REF], which defines the dynamic continuation method presented in Section 4, may still be of help in reconstructing the smaller growth rate at least in some cases. Suppose that the graphs µ 1 and µ 2 cross each other at some value s c , say, µ 1 (s) < µ 2 (s) for s < s c but µ 1 (s) > µ 2 (s) for s > s c . If one treats s as a parameter then the equilibria E *

1 and E * 2 undergo an exchange of stability (a degenerate transcritical bifurcation) at s = s c . As the continuation rule (12) lets s drift slowly (with speed ε) the full system exhibits a phenomenon known as delayed loss of stability [START_REF] Boudjellaba | Stability Loss Delay in Harvesting Competing Populations[END_REF] in the context of dynamic bifurcations [START_REF] Benoît | Dynamic bifurcations: proceedings of a conference held in Luminy[END_REF], widely studied in slow-fast systems. Physically this means (say, we are increasing s slowly from below s c ) that concentration b 2 comes close to 0 while s < s c . Then, when s has crossed s c , the concentration b 2 grows exponentially, but still takes some time until it reaches values noticeably different from 0. The value sloss of s at which b 2 becomes noticeably non-zero is in the ideal ODE model independent of the drift speed ε of s. Figure 7 in Section 7 shows this effect: since b 2 is nearly zero the variable D continues to follow the, by now unstable, drifting equilibrium E * 1 . This delay mechanism allows one in principle a reconstruction of a part of the slower growth rate close to the about the bifurcation value of s. This is in line with the philosophy behind continuation methods, which allow one to track both stable and unstable solution branches, hence, to go past bifurcation points.

The numerical simulations show only the case of continuously drifting s. Setting the value of s at discrete times (as one would do in the discrete-time framework of Section 5), achieves the same delay effect.

Figures are given at the end of the paper.

The one species case

For the unknown function, we have chosen the non-monotonic Haldane function µ(s) = s 1 + s + 10s 2 and s in = 1. The grey background curve in Figure 3 shows µ(•), which is clearly non-monotonic on the domain (0, s in ) (s in = 1). Figure 3 shows the stabilization of the control law [START_REF] Dochain | On-line estimation of specific growth rates : an illustrative case study[END_REF] for s in the increasing (left panel) and decreasing (right panel) part of the growth law µ. Figure 4 shows how the continuation [START_REF] Dochain | Dynamical Modelling, Analysis, Monitoring and Control Design for Nonlinear Bioprocesses[END_REF] reconstructs the entire graph of the growth function. For the left panel, ε = +0.01, and for the right panel, ε = -0.01.

Figure 5 shows the output of a simulation with the discrete-time adaptation method proposed in Section 5. The black dots correspond to values at which the control was accepted as non-invasive. Using continuation one achieves small and rapidly decaying transients in every evaluation of s eq (which involves running system (1) with control until transients have settled). This is so because the transients all lie inside the subspace {(s, b) : s + b = s in } after system (1) has run at least once. Second, the initial offset from the equilibrium is always small, because the adjustments of s and D are small.

Figure 6 demonstrates the speed-up using the simplified scheme ( 14)-( 15) (note the times at the abscissae). The difference to Figure 5 is that the values at which the growth function is evaluated are not exactly equidistantly spaced.

The two species case

We have chosen two monotonic Monod functions

µ 1 (s) = s 0.1 + s , µ 2 (s
) = 1.5 s 0.4 + s and s in = 1. Their graphs are shown as faint curves in the background of Figure 7. The simulations depicted in Figure 7 have been made with the drift speed parameter ε equal to ±0.01. One can see that the delayed loss of stability permits us to reconstruct parts of the graphs of the growth curves that do not belong to their supremum envelope.

From a practical view point, the sudden "jump" between the two graphs could inform of the presence of another species, if one belongs that the culture in chemostate was initially pure.

Remark. If the drift speed ε of the continuation is too small the concentration of the initially suppressed species (for example, b 2 in the left panel of Figure 7) will become unrealistically small in the numerical simulations. In a real experiment this would correspond to the species becoming almost extinct. This makes the validity of the ODE model ( 16) with two species questionable as a good representation of reality [START_REF] Campillo | Effect of population size in a Prey-Predator model[END_REF] such that it is difficult to conclude from numerical simulations of the ODE model ( 16) how large the delay effect is in real experiments. This question will the matter of future investigations.

Conclusion

In this work, we have presented a new framework for the functional identification of non-monotonic growth functions in the chemostat. Continuous and discrete time variants have been proposed. Numerical simulations illustrate the potential of the method. Further investigations are required to find out which approach (discrete or continous) is better suited for real experiments.

The approach is more general than the case we have presented here for the chemostat model. We introduce feedback control to reduce a dynamical system to an algebraic equation. Then we apply either numerical methods (in discrete time) or dynamical equations (in continuous time) to identify unstable equilibria of the original uncontrolled dynamical system. We use the chemostat as a conceptually simple example that is still of practical interest.

Another application we plan to explore in the future are regulation problems. For example, one can regulate the single-species chemostat to operate at the substrate concentration s at which the growth rate µ is maximal by following the same recipe. This approach to regulation, which is similar in spirit to the act-and-wait technique for delay compensation [START_REF] Insperger | Act-and-wait concept for continuous-time control systems with feedback delay[END_REF], does not require an a-priori identification of the growth rate µ, and leads to a different algorithm than the methods discussed in the literature [START_REF] Dochain | Extremum Seeking Control and its Application to Process and Reaction Systems: A Survey[END_REF][START_REF] Guay | Adaptive extremum seeking control of stirred tank bioreactors[END_REF]).

Figure 1 :

 1 Figure 1: Domains of stability and instability in open-loop

µ

  ) exponentially stabilizes the system[START_REF] Abed | Stabilization of period doubling bifurcations and implicatons for control of chaos[END_REF] locally about (s, b) = (s, s in -s), for any positive constants (G 1 , G 2 ) such that G 1 > -µ ′ (s).Furthermore one has lim t→+∞ D(t) = µ(s) Note that the assumptions in Proposition 2 (for example, on the gain G 1 ) are weaker than those of Proposition 1 as Proposition 2 is only concerned with local stability and a single reference value s. ( s ) , ( s , b ) u n s t a b l e slopes determining lower bound for gain G 1 grey area = set of admissible ( s, D), s min s (example) d (example) µ ( s ) , ( s , b ) s t a b l e

Figure 2 :

 2 Figure 2: Sketch illustrating shape of growth function µ and set of admissible reference values

Figure 3 :Figure 4 :

 34 Figure 3: Simulations for two different values of s in the (s, D)-projection.

Figure 5 :D

 5 Figure 5: Simulation of discrete-time continuation using (13) where sk = 0.05k. Panel (a) shows the time profile of output s and input D throughout the run. Panel (b) shows the same in the (s, D)-plane. Black dots indicate when convergence (|s eq -sk | < tol) was reached, and the iteration moved on to the next sk . Parameters: tol = 10 -4 , G 1 = 1/(s in -sk ), D max = 0.2, D min = 0.02. Transients were accepted as settled if max s -min s < 0.1 on the interval [t -10, t].

Figure 6 :Figure 7 :

 67 Figure 6: Simulation of simplified discrete-time continuation using (14)-(15) where δ = 0.05. Note the shorter time scales on the x-axis of panel (a) compared to Figure 5. The convention in the panels and the parameters are identical to Figure 5.


Francesco Ancilotto

Manuel ; Barranco
email: manuel@ecm.ub.es

Jesús Navarro

Martí Pi

Jesús M Navarro

Cavitation of electron bubbles in liquid parahydrogen

Keywords: liquid parahydrogen, electron bubbles, density functional theory, capillary approximation liquid parahydrogen, electron bubbles, density functional theory, capillary approximation

Introduction

Excess electrons in cryogenic liquids like helium and hydrogen have been the subject of many theoretical and experimental studies. It is well established that when an electron is injected into such liquids, it looses kinetic energy by ionization and excitation of the liquid atoms, and by the production of elementary excitations in the liquid. When the electron has lost most of its kinetic energy and moves subsonically, it eventually produces a cavity in the fluid and gets trapped inside it. We will refer to these structures as electron bubbles (e-bubbles). Many works have been devoted to measure the mobility of electron bubbles in liquid helium [1][2][START_REF] Eden | [END_REF][4] and hydrogen [5][6][7][8][9].

In the case of helium, there is a wealth of experimental and theoretical work on heterogeneous cavitation produced by excess electrons when the fluid is driven below the saturated vapor pressure by the application of ultrasonic pulses caused by a transducer [10]. At given temperature (T ), there is a limiting (critical) pressure below which electron bubbles explode triggering liquid-gas phase separation. In the case of helium, this process is well known and the agreement between theory and experiment is very satisfactory [11,12].

It was theoretically predicted that electron bubbles should develop in liquid hydrogen [13], and they were experimentally detected shortly after [5]. Mobility experiments [14] disclosed an anomalously high mobility of negative charges in the T range from 17 to 22 K that was attributed to the existence of H -anions in the liquid, as discussed in detail in [7]. It thus appears that in liquid hydrogen two types of negatively charged bubbles may develop when electrons are injected inside it, i.e. ordinary e-bubbles and 'H --bubbles'.

As compared with liquid helium, much less is known on cavitation caused by e-bubbles in liquid hydrogen. We have recently addressed homogeneous nucleation and cavitation in parahydrogen [15] (pH 2 ) using a zero-range finite-temperature density functional [16] (DF) that reproduces the properties of liquid pH 2 in the pressure (P ) and T ranges corresponding to the liquid-vapor phase equilibrium region, i.e., from the triple to the critical point. It was found that the spinodal pressure P sp near the triple point was about -80 bar, whereas the homogeneous cavitation pressure (that at which gas bubbles are created in the liquid at a rate large enough so as to trigger the liquid-vapor phase transition) in the same region of the phase diagram was about -70 bar. This value is very much affected by the presence of electrons (heterogeneous cavitation). Indeed, we have found [15] that in the triple point region the pressure for e-bubble explosion is just P expl ∼ -30 bar. That estimate was obtained within the capillary model that describes the electron bubble as a sharp spherical cavity and neglects the thickness of the cavity inner surface, whose width drastically increases as one approaches the critical point. As a result, in the capillary model approximation P expl eventually becomes more negative than P sp , which is unphysical.

In this work we want to determine P expl more accurately by going beyond the crude capillary model. To this end, we resort to the finite-range density functional (FRDF) of [16], which has proven to accurately describe the temperature dependent, finite width of the liquid-vapor interface, supplementing it with an electron-pH 2 interaction obtained within a Hartree-type model, which has been successfuly used in the case of liquid helium [23]. This work is organized as follows. In Sec. 2 we present the conventional capillary approximation to electron bubbles in cryogenic fluids. The density functional approach for pH 2 is briefly reviewed in Sec. 3, together with the determination of the e-pH 2 Hartree interaction. The results are presented and discussed in Sec. 4. Finally, a summary is presented in Sec. 5.

Capillary approximation

The capillary approximation is found to work fairly well for electron bubbles in liquid helium, especially at low temperatures [11] where the interface separating the electron bubble cavity from the liquid is rather sharp due to the pressure exerted by the electron against the surrounding fluid.

Electron bubbles in hydrogen have been detected long time ago [5]. Their radius have been calculated in [13] using a model characterized by an abrupt interface for the spherical cavity where the electron is localized and taking into account in a simplified way the polarization induced by the electron in the liquid. The model in [13] was supplemented by experimental informations for the electron-H 2 scattering length. The e-bubble radius was found to vary between 10 and 12.5 Å (in the temperature range between 20 K and the critical temperature), that is roughly a factor of two smaller than in liquid 4 He.

In the capillary approximation, the energy of an electron bubble is written as the following sum [START_REF] Konstantinov | [END_REF] whose four terms correspond, respectively, to the electron energy in an impenetrable spherical well potential of radius R, the surface and volume energies and a polarization contribution. In that equation, m e is the electron mass, γ is the surface tension of the liquid-vapor interface, P is the external pressure, and ǫ is the dielectric constant of pH 2 . The minimization of the energy with respect to the bubble radius allows one to obtain the equilibrium radius from the equation

E = π 2 2 2m e R 2 + 4πR 2 γ + 4π 3 R 3 P - e 2 2R ǫ -1 ǫ , (1)
- π 2 4m e R 5 + 2γ R + P + e 2 8πR 4 ǫ -1 ǫ = 0 , (2)
which has a solution at any non-negative pressure. Considering for instance T = 14 K, were the experimental values of the required inputs are P = 0.079 bar (saturated vapor pressure), γ = 2.958 dyne/cm and ǫ = 1.251 [19] one finds R eq = 10.7 Å. It is worth noticing that the polarization energy term has some effect. Had we neglected it, the calculated equilibrium radius would have been ∼ 0.6 Å larger.

If the liquid is submitted to a tensile strength that brings it below saturated vapor pressure, the electron bubble becomes metastable. A potential energy barrier prevents the bubble explosion. This barrier decreases as the pressure does and eventually vanishes at P expl . This is the lowest pressure that can be attained before the electron bubble explodes, that is, expands without limits triggering the phase transition. P expl can be obtained from the conditions dE/dR = 0, and d 2 E/dR 2 = 0, which lead to the following system of coupled equations

- π 2 2 m e R 3 + 8πγR + 4πR 2 P + e 2 2R 2 ǫ -1 ǫ = 0 (3) 3π 2 2 m e R 4 + 8πγ + 8πRP - e 2 R 3 ǫ -1 ǫ = 0 .
An analytical solution is obtained if the polarization term is neglected, and we shall use it to quickly illustrate how the cavitation process is affected by the presence of impurities. One gets R = [5π 2 /(8m e γ)] 1/4 and P = -(8/5)×[8m e /(5π 2 )] 1/4 γ 5/4 . According to this, at T = 14 K liquid pH 2 becomes macroscopically unstable at P expl = -28.1 bar, much higher than the homogeneous cavitation pressure, which is about -70 bar. The P expl curve obtained from Eqs. (3) is displayed in Fig. 5 as a dot-dashed line. It can be seen that it leads to the unphysical result that P expl < P sp for T ≥ 28 K. Analogously to the helium case, this failure is related to the presence of a sizeable gas-like region within the electron bubble [11], which cannot be described within the simple capillary approximation. Any improvement in the method used to address cavitation induced by electrons calls for a more realistic description of the bubble-liquid interface, together with using an electron-pH 2 interaction that reproduces the energy barrier encountered by an electron entering liquid pH 2 , much along what has been done for liquid helium [10,11]. This is the aim of the next Section.

Density functional for parahydrogen

The basic theoretical ingredient in our approach is the total free energy F of the system written in terms of a temperature-dependent free energy density functional

F = drf (ρ, T ) ,
where ρ is the pH 2 number density. All the required thermodynamical quantities can be obtained from f (ρ, T ).

We have recently proposed a free energy density functional for pH 2 in the form [16]

f (ρ, T ) = f ni (ρ, T ) + f c (ρ, T ) . (4)
In this expression, f ni (ρ, T ) is the free energy density for an ideal Bose gas (see e.g. [START_REF] Huang | Statistical Mechanics[END_REF]) and f c (ρ, T ) is the correlation energy density, which incorporates the correlations induced by the interaction. The correlation part f c (ρ, T ) is written as a polynomial in powers of the density

f c (ρ, T ) = a 1 (T )ρ 2 + a 2 (T )ρ 3 + a 3 (T )ρ 4 + a 4 (T )ρ 5 . (5)
The a i (T ) parameters have been fixed so as to reproduce a series of thermodynamics properties along the liquid-gas equilibrium line, namely the equality of pressure and chemical potential of both phases, and the pressure and speed of sound of bulk liquid pH 2 at the saturated vapor pressure [20]. In Fig. 1 are plotted the resulting parameters, transformed into dimensionless quantities by means of the standard Lennard-Jones parameters for pH 2 , namely ǫ = 34.2 K and σ = 2.96 Å. As such, this functional is a zero-range one, and finite range effects have been considered in [16] to account for highly inhomogenous configurations, as those appearing when either strongly repulsive electrons or very attractive impurities are present in the liquid. These effects have been included by introducing two modifications into the zero-range DF. Firstly, the term in ρ 2 has been replaced with

a 1 (T )ρ 2 → 1 2 dr ′ ρ(r)V (|r -r ′ |)ρ(r ′ ) , (6)
where V represents the pH 2 -pH 2 interaction screened at distances shorther than a screening parameter h(T ). It has a Lennard-Jones form if r ≥ h(T ) and V (r) = 0 otherwise. Secondly, a coarse-grained density is introduced in the remaining terms of the functional in the following way

a i (T )ρ 1+i → a i (T )ρ ρi , ( 7
)
where ρ is an average of the density in a small sphere whose radius is of the order of the screening distance. This is a well-known prescription for classical liquids [22], which has been also applied to liquid 4 He. We refer the reader to [16] for the details. It turns out that the finite-range DF (FRDF) so devised accurately reproduces the experimental surface tension of the liquid-vapor interface all over the phase-coexistence region [21].

Once the density functional f (ρ, T ) has been determined, the thermodynamic relationship P = -f (ρ, T ) + µρ (µ being the pH 2 chemical potential) allows one to calculate the equation of state in the negative pressure regime, inaccessible to the experimental determination, and to obtain the spinodal line in bulk liquid pH 2 . In the presence of electrons, the free energy of the system is a functional of the pH 2 density ρ, the excess electron wave function Ψ, and T , which we write as

F [ρ, Ψ, T ] = dr f (ρ, T )+ 2 2 m e dr |∇Ψ(r )| 2 + dr |Ψ(r )| 2 dr ′ ρ(r ′ )V e-H2 (r-r ′ ) , (8)
supplementing the FRDF with two more terms that accounts for the electron kinetic energy and the e-pH 2 interaction. The latter has been written as a function of the local hydrogen density [23] dr |Ψ(r )| 2 V [ρ(r)], where

V (ρ) = 2 k 2 0 2 m e + 2π 2 m e ρ a α -2παe 2 4π 3 1/3 ρ 4/3 . (9)
α = 0.801 Å3 is the static isotropic polarizability of a H 2 molecule [24], and k 0 is determined from the pH 2 local Wigner-Seitz radius r s = (3/4πρ) 1/3 by solving the equation

tan[k 0 (r s -a c )] = k 0 r s . (10)
In the above equations, a c and a α are the scattering lenghts arising from the hardcore and from the polarization potential. Their sum a = a c +a α is equal to the e-H 2 total scattering length, which we take equal to the experimental value a = 0.672 Å [24]. We fix a c in such a way to reproduce the experimental value for the total energy V 0 of a delocalized electron in bulk liquid pH 2 , V 0 = 1.9 eV at 20 K [7]. We find a c = 0.979 Å, and consequently a α ≡ a -a c = -0.307 Å. The above procedure works well for slowly varying density profiles. It lacks however an additional non-local contribution which comes from the long-range polarization energy, namely

U nl (r) = - αe 2 2 dr ′ ρ(r) -ρ(r ′ ) |r -r ′ | 4 . ( 11
)
Such a term is not usually included in DF calculations due to its computational complexity, and we will neglect it in the present work as well.

To obtain P expl we have proceeded as follows. For given P and T values we have first determined the bulk liquid density ρ b from the theoretical equation of state. Next, we have solved the Euler-Lagrange equations that result from the variation of the constrained grand potential density ω(ρ, Ψ, T ) ≡ ω(ρ, Ψ, T ) -ε|Ψ| 2 , where the grand potential density ω(ρ, Ψ, T ) is derived from Eq. ( 8)

ω(ρ, Ψ, T ) = f (ρ, T ) + 2 2 m e |∇Ψ| 2 + |Ψ| 2 V (ρ) -µρ . ( 12
)
It yields where ε is the lowest eigenvalue of the Schrödinger equation obeyed by the electron. These equations are solved in a three-dimensional box using Fast-Fourier techniques [25] imposing that ρ(r) = ρ b at the cell boundaries, and that the electron is in the 1s state. Fixing T , the above equations are solved for decreasing values of P below the saturated vapor pressure P svp (T ). The lower the pressure, the larger the metastable electron cavity. Eventually, a value of the pressure is found for which the procedure is unable to find a metastable e-pH 2 configuration, and the electron cavity expands until it reaches the boundaries of the simulation cell. This pressure value numerically defines P expl .

δf δρ + |Ψ| 2 ∂V ∂ρ = µ (13) - 2 2 m e ∆Ψ + V [ρ]Ψ = εΨ , (14)
As an alternative to the use of the Hartree-type interaction described above, one may take advantage of the large short-range exchange repulsion between the electron and the pH 2 molecule, that dominates the long-range polarization attraction, and consider in Eq. ( 8) a much simpler "contact" interaction

V e-H2 (r -r ′ ) = V 0 ρ b δ(r -r ′ ) (15)
whose intensity scales with the liquid density so that it yields the experimental energy barrier V 0 encountered by an electron entering liquid pH 2 , V 0 = 1.9 eV at T = 20 K [7]. We recall that ρ b is the bulk pH 2 density at such temperature. This type of interaction was already used for liquid 4 He [10], and the results obtained were found to be very similar to those obtained using the Hartree-type interaction term [11]. We have checked that the results obtained using both type of interactions are indeed very similar also in the case of pH 2 .

Results and discussion

Figure 2 shows the calculated e-bubble radii at different temperatures (squares), using the approach described in Sec. 3. The radius is defined as the radial position where the density attains the value ρ b /2, ρ b being again the bulk liquid density at the investigated temperature. Our results are compared in this figure with the estimates from [13]. At low T , where our results are expected to be very similar to the ones obtained within the capillary approximation used in [13], the bubble radius is slightly overestimated, and rather coincides with the value obtained in [13] when the long-range polarization interaction term is omitted. This is in part a consequence of our own neglect of the long-range polarization term Eq. ( 11). At higher temperatures the capillary model of [13] breaks down due to the increasing thickness of the liquid-vapor interface in the bubble cavity, which is instead correctly accounted for in our calculations. The density profiles of the pH 2 at different temperatures are shown in Fig. 3. It can be seen that the surface thickness of the e-bubble increases as T does.

Figure 4 shows some density profiles illustrating the cavitation process. For each T , we show the electron and pH 2 profiles corresponding to the stable e-bubble at saturated vapor pressure P svp , and those corresponding to the nearly critical bubble at P expl .

Finally, we display in Fig. 5 the calculated P expl within the FRDF plus e-pH 2 Hartree potential approach and the capillary model as well, together with the saturated vapor pressure, spinodal and homogeneous cavitation curves. It can be seen how the selfconsistent solution of the liquid and excess electron equations circumvents the limitations of the capillary model at high temperatures and yields sensible results for the e-bubble explosion all over the temperature range where the calculations have been carried out. At lower temperatures, the slight deviations between the FRDF results and the capillary model predictions are due in part to our neglect of the long-range polarization energy term, as discussed previously.

Summary

We have investivated heterogeneous cavitation in liquid pH 2 triggered by the presence of excess electrons. To this end, we have built a Hartree-like electronparahydrogen interaction, and have used a previously devised finite-range density functional for pH 2 that works well in the liquid-vapor equilibrium region.

Our approach allows for a realistic and flexible description of the electron cavity and of the interaction between the electron and the surrounding liquid. It does not impose a priori the density profile of the critical bubble, allowing for a proper description of the process from the saturated vapor pressure down to the spinodal line.

We have found that the critical pressure at which electron bubbles explode is in absolute value about a factor of two smaller than the homogeneous cavitation pressure of liquid pH 2 . We have also used a capillary model to discuss the onset of electron bubble explosion in liquid pH 2 , and have discussed the limitations of this model at high temperatures. Electron bubble radius at different temperatures. The solid lines show the results of [13] corresponding to two experimental estimates for the electron-H 2 scattering length. The dashed lines are the results of [13], where the long-range polarization term has been omitted. The squares show the values calculated using the FRDF with the effective Hartree-type electron-pH 2 interaction described in the text.

  Barranco, Jesús Navarro and Martí Pi

4

 4 Francesco Ancilotto, Manuel Barranco, Jesús Navarro and Martí Pi f (ρ, T )

Figure 1 .

 1 Figure 1. Dimensionless parameters a i (T )/(ǫσ 3i ) of the pH 2 density functional of Eq. (5) as a function of temperature. The values ǫ = 34.2 K and σ = 2.96 Å have been used.

Figure 2 .

 2 Figure2. Electron bubble radius at different temperatures. The solid lines show the results of[13] corresponding to two experimental estimates for the electron-H 2 scattering length. The dashed lines are the results of[13], where the long-range polarization term has been omitted. The squares show the values calculated using the FRDF with the effective Hartree-type electron-pH 2 interaction described in the text.

Figure 3 . 11 Figure 4 .Figure 5 .

 31145 Figure 3. Density profiles of pH 2 bubble cavities at saturated vapor pressure conditions and different temperatures. From top to bottom, T = 20, 22, 24, 26, 28, 30, and 32 K.

Acknowledgments

This work has been performed under Grants No. FIS2008-00421/FIS and FIS2007-60133 from DGI, Spain (FEDER), and Grant 2009SGR01289 from Generalitat de Catalunya.


Patrik Lambert
email: patrik.lambert@lium.univ-lemans.fr

Jean Senellart
email: senellart@systran.fr

Laurent Romary
email: laurent.romary@inria.fr

Holger Schwenk
email: holger.schwenk@lium.univ-lemans.fr

Florian Zipser
email: f.zipser@gmx.de

Patrice Lopez
email: patrice.lopez@inria.fr

Frédéric Blain
email: frederic.blain@lium.univ-lemans.fr

Collaborative Machine Translation Service for Scientific texts

French researchers are required to frequently translate into French the description of their work published in English. At the same time, the need for French people to access articles in English, or to international researchers to access theses or papers in French, is incorrectly resolved via the use of generic translation tools. We propose the demonstration of an end-to-end tool integrated in the HAL open archive for enabling efficient translation for scientific texts. This tool can give translation suggestions adapted to the scientific domain, improving by more than 10 points the BLEU score of a generic system. It also provides a post-edition service which captures user post-editing data that can be used to incrementally improve the translations engines. Thus it is helpful for users which need to translate or to access scientific texts.

Introduction

Due to the globalisation of research, the English language is today the universal language of scientific communication. In France, regulations require the use of the French language in progress reports, academic dissertations, manuscripts, and French is the official educational language of the country. This situation forces researchers to frequently translate their own articles, lectures, presentations, reports, and abstracts between English and French. In addition, students and the general public are also challenged by language, when it comes to find published articles in English or to understand these articles. Finally, international scientists not even consider to look for French publications (for instance PhD theses) because they are not available in their native languages. This problem, incorrectly resolved through the use of generic translation tools, actually reveals an interesting generic problem where a community of specialists are regularly performing translations tasks on a very limited domain. At the same time, other communities of users seek translations for the same type of documents. Without appropriate tools, the expertise and time spent for translation activity by the first community is lost and do not benefit to translation requests of the other communities.

We propose the demonstration of an end-to-end tool for enabling efficient translation for scientific texts. This system, developed for the COSMAT ANR project, 1 is closely integrated into the HAL open archive, 2 a multidisciplinary open-access archive which was created in 2006 to archive publications from all the French scientific community. The tool deals with handling of source document format, generally a pdf file, specialised translation of the content, and user-friendly userinterface allowing to post-edit the output. Behind the scene, the post-editing tool captures user postediting data which are used to incrementally improve the translations engines. The only equipment required by this demonstration is a computer with an Internet browser installed and an Internet connection.

In this paper, we first describe the complete work-flow from data acquisition to final postediting. Then we focus on the text extraction procedure. In Section 4, we give details about the translation system. Then in section 5, we present the translation and post-editing interface. We finally give some concluding remarks.

The system will be demonstrated at EACL in his tight integration with the HAL paper deposit system. If the organizers agree, we would like to offer the use of our system during the EACL conference. It would automatically translate all the abstracts of the accepted papers and also offers the possibility to correct the outputs. This resulting data would be made freely available.

Complete Processing Work-flow

The entry point for the system are "ready to publish" scientific papers. The goal of our system was to extract content keeping as many metainformation as possible from the document, to translate the content, to allow the user to perform post-editing, and to render the result in a format as close as possible to the source format. To train our system, we collected from the HAL archive more than 40 000 documents in physics and computer science, including articles, PhD theses or research reports (see Section 4). This material was used to train the translation engines and to extract domain bilingual terminology.

The user scenario is the following:

• A user uploads an article in PDF format 3 on the system.

• The document is processed by the opensource Grobid tool (see section 3) to extract 3 The commonly used publishing format is PDF files while authoring format is principally a mix of Microsoft Word file and LaTeX documents using a variety of styles. The originality of our approach is to work on the PDF file and not on these source formats. The rationale being that 1/ the source format is almost never available, 2/ even if we had access to the source format, we would need to implement a filter specific to each individual template required by such or such conference for a good quality content extraction the content. The extracted paper is structured in the TEI format where title, authors, references, footnotes, figure captions are identified with a very high accuracy.

• An entity recognition process is performed for markup of domain entities such as: chemical compounds for chemical papers, mathematical formulas, pseudo-code and object references in computer science papers, but also miscellaneous acronyms commonly used in scientific communication.

• Specialised terminology is then recognised using the Termsciences4 reference terminology database, completed with terminology automatically extracted from the training corpus. The actual translation of the paper is performed using adapted translation as described in Section 4.

• The translation process generates a bilingual TEI format preserving the source structure and integrating the entity annotation, multiple terminology choices when available, and the token alignment between source and target sentences.

• The translation is proposed to the user for post-editing through a rich interactive interface described in Section 5.

• The final version of the document is then archived in TEI format and available for display in HTML using dedicated XSLT style sheets.

The Grobid System

Based on state-of-the-art machine learning techniques, Grobid [START_REF] Lopez | GROBID: Combining automatic bibliographic data recognition and term extraction for scholarship publications[END_REF] performs reliable bibliographic data extraction from scholar articles combined with multi-level term extraction. These two types of extraction present synergies and correspond to complementary descriptions of an article. This tool parses and converts scientific articles in PDF format into a structured TEI document5 compliant with the good practices developed within the European PEER project [START_REF] Bretel | Back to meaning -information structuring in the PEER project[END_REF]. Grobid is trained on a set of annotated scientific article and can be re-trained to fit templates used for a specific conference or to extract additional fields.

Translation of Scientific Texts

The translation system used is a Hybrid Machine Translation (HMT) system from French to English and from English to French, adapted to translate scientific texts in several domains (so far physics and computer science). This system is composed of a statistical engine, coupled with rule-based modules to translate special parts of the text such as mathematical formulas, chemical compounds, pseudo-code, and enriched with domain bilingual terminology (see Section 2). Large amounts of monolingual and parallel data are available to train a SMT system between French and English, but not in the scientific domain. In order to improve the performance of our translation system in this task, we extracted in-domain monolingual and parallel data from the HAL archive. All the PDF files deposited in in computer science and physics were made available to us. These files were then converted to plain text using the Grobid tool, as described in the previous section. We extracted text from all the documents from HAL that were made available to us to train our language model. We built a small parallel corpus from the abstracts of the PhD theses from French universities, which must include both an abstract in French and in English. Table 1 presents statistics of these in-domain data.

The data extracted from HAL were used to adapt a generic system to the scientific literature domain. The generic system was mostly trained on data provided for the shared task of Sixth Workshop on Statistical Machine Translation6 (WMT 2011), described in Table 2.

Table 3 presents results showing, in the English-French direction, the impact on the statistical engine of introducing the resources extracted from HAL, as well as the impact of domain adaptation techniques. The baseline statistical engine is a standard PBSMT system based on Moses [START_REF] Koehn | Moses: Open source toolkit for statistical machine translation[END_REF] and the SRILM tookit [START_REF] Stolcke | SRILM: an extensible language modeling toolkit[END_REF]. Is was trained and tuned only on WMT11 data (out-of-domain). Incorporating the HAL data into the language model and tuning the system on the HAL development set, tracted from all documents in HAL, in computer science (cs) and physics (phys). The following statistics are given for the English (En) and French (Fr) sides (Lg) of the corpus: the number of sentences, the number of running words (after tokenisation) and the number of words in the vocabulary (M and k stand for millions and thousands, respectively).

yielded a gain of more than 7 BLEU points, in both domains (computer science and physics). Including the theses abstracts in the parallel training corpus, a further gain of 2.3 BLEU points is observed for computer science, and 3.1 points for physics. The last experiment performed aims at increasing the amount of in-domain parallel texts by translating automatically in-domain monolingual data, as suggested by [START_REF] Schwenk | Investigations on large-scale lightly-supervised training for statistical machine translation[END_REF]. The synthesised bitext does not bring new words into the system, but increases the probability of indomain bilingual phrases. By adding a synthetic bitext of 12 million words to the parallel training data, we observed a gain of 0.5 BLEU point for computer science, and 0.7 points for physics.

Although not shown here, similar results were obtained in the French-English direction. The French-English system is actually slightly better than the English-French one as it is an easier translation direction.

Post-editing Interface

The collaborative aspect of the demonstrated machine translation service is based on a post-editing tool, whose interface is shown in Figure 1. This tool provides the following features:

• WYSIWYG display of the source and target texts (Zones 1+2)

• Alignment at the sentence level (Zone 3)

• Zone to review the translation with alignment of source and target terms (Zone 4) and terminology reference (Zone 5)

• Alternative translations (Zone 6)

The tool allows the user to perform sentence level post-editing and records details of postediting activity, such as keystrokes, terminology selection, actual edits and time log for the complete action.

Conclusions and Perspectives

We proposed the demonstration of an end-to-end tool integrated into the HAL archive and enabling efficient translation for scientific texts. This tool consists of a high-accuracy PDF extractor, a hybrid machine translation engine adapted to the scientific domain and a post-edition tool. Thanks to in-domain data collected from HAL, the statistical engine was improved by more than 10 BLEU points with respect to a generic system trained on WMT11 data.

Our system was deployed for a physic conference organised in Paris in Sept 2011. All accepted abstracts were translated into author's native languages (around 70% of them) and proposed for post-editing. The experience was promoted by the organisation committee and 50 scientists volunteered (34 finally performed their post-editing). The same experience will be proposed for authors of the LREC conference. We would like to offer a complete demonstration of the system at EACL. The goal of these experiences is to collect and distribute detailed "post-editing" data for enabling research on this activity.

Table 3 :

 3 Results (BLEU score) for the English-French systems. The type of parallel data used to train the translation model or language model are indicated, as well as the set (in-domain or out-of-domain) used to tune the models. Finally, the number of words in the parallel corpus and the BLEU score on the in-domain test set are indicated for each domain: computer science and physics.

Figure 1 :

 1 Figure 1: Translation and post-editing interface.

Table 1 :

 1 Statistics for the parallel training, development, and test data sets extracted from thesis abstracts contained in HAL, as well as monolingual data ex-

	Parallel data
	Train cs+phys En 55.9 k 1.41 M 43.3 k
		Fr 55.9 k 1.63 M 47.9 k
	Dev cs	En 1100 25.8 k 4.6 k
		Fr	1100 28.7 k 5.1 k
	phys	En 1000 26.1 k 5.1 k
		Fr	1000 29.1 k 5.6 k
	Test cs	En 1100 26.1 k 4.6 k
		Fr	1100 29.2 k 5.2 k
	phys	En 1000 25.9 k 5.1 k
		Fr	1000 28.8 k 5.5 k
	Monolingual data
	Train cs	En 2.5 M	54 M 457 k
		Fr 761 k	19 M 274 k
	phys	En 2.1 M	50 M 646 k
		Fr 662 k	17 M 292 k

Table 2 :

 2 Out-of-domain development and training data used (number of words after tokenisation).

http://www.cosmat.fr/

http://hal.archives-ouvertes.fr/?langue=en

http://www.termsciences.fr

http://www.tei-c.org

http://www.statmt.org/wmt11/translation-task.html Set Domain Lg Sent. Words Vocab.

Acknowledgements

This work has been partially funded by the French Government under the project COSMAT (ANR ANR-09-CORD-004).