Alain Rapaport
email: rapaport@supagro.inra.fr

Ihab Haidar
email: ihabhaidar@yahoo.com

Frédéric Gérard

Effects of spatial structure and diffusion on the performances of the chemostat

Keywords: Chemostat model, interconnection, diffusion, global stability

Given hydric capacity and nutrient flow of a chemostat-like system, we analyse the influence of a spatial structure on the output concentrations at steady-state. Three configurations are compared: perfectly-mixed, serial and parallel with diffusion rate. We show the existence of a threshold on the input concentration of nutrient for which the benefits of the serial and parallel configurations over the perfectlymixed one are reversed. In addition, we show that the dependency of the output concentrations on the diffusion rate can be non-monotonic, and give precise conditions for the diffusion effect to be advantageous. The study encompasses dead-zone models.

Introduction

The chemostat is a popular apparatus, invented simultaneously by Monod [START_REF] Monod | La technique de la culture continue: théorie et applications[END_REF] and Novick & Szilard [START_REF] Novick | Szilard Description of the chemostat[END_REF], for the so-called continuous culture of micro-organisms. It has the advantage to study bacteria growth at steady state, in contrast to batch cultivation. In the classical experiments, the medium is assumed to be perfectly mixed, that justifies mathematical models described by systems of ordinary differential equations [START_REF] Smith | The theory of chemostat, dynamics of microbial competition[END_REF]. The chemostat model is also used in ecology for studying populations of micro-organisms, such as lake plankton or wetlands ecosystems. In natural ecosytems, or in industrial applications that use large bioreactors, the assumption of perfectly mixed medium is questionable. This is why spatial considerations have been introduced in the classical model of the chemostat, such as the gradostat model [START_REF] Lovitt | The gradostat: a bidirectional compound chemostat and its applications in microbial research[END_REF] that is a series of interconnected chemostats (of identicalvolumes). Segregated habitats are also considered in lakes, where the bottom can be modeled as a dead zone and nutrient mixing between the two zones is achieved by diffusion rate [START_REF] Nakaoka | Competition in chemostat-type equations with two habitats[END_REF]. The consideration of dead zones is also often used in bioprocesses modelling [START_REF] Levenspiel | Chemical reaction engineering[END_REF][START_REF] Hu | Cellular Bioprocess Technology, Fundamentals and Frontier[END_REF][START_REF] Grobicki | Hydrodynamic characteristics of the anaerobic baffled reactor[END_REF][START_REF] Roux | Adaptive non-linear control of a continuous stirred tank bioreactor[END_REF][START_REF] Roca | Hydraulic model of a gas-lift bioreactor with flocculating yeast[END_REF][START_REF] Valdes-Parada | An approximate solution for a transient two-phase stirred tank bioreactor with nonlinear kinetics[END_REF][START_REF] Saddoud | A mathematical study of an enzymatic hydrolysis of a cellulosic substrate in non homogeneous reactors[END_REF].

Series of chemostats, instead of single chemostat, have shown to potentially improve the performances of bioprocesses, reducing the total residence time [START_REF] Hill | Minimum tank volumes for CFST bioreactors in series[END_REF][START_REF] Luyben | Optimal design for continuously stirred tank reactors in series using Michaelis-Menten kinetics[END_REF][START_REF] Harmand | Optimal design of two interconnected bioreactors-some new results[END_REF][START_REF] Harmand | Optimal design of two interconnected enzymatic reactors[END_REF][START_REF] Harmand | Towards a unified approach for the design of interconnected catalytic and auto-catalytic reactors[END_REF] or allowing species persistence [START_REF] Stephanopoulos | Effect of inhomogeneities on the coexistence of competing microbial populations[END_REF][START_REF] Rapaport | Coexistence in the design of a series of two chemostats[END_REF]. These properties have of course economical impacts for the biotechnological industry, and there is a significant literature on the design of series of reactors and comparison with plug-flow reactors (that can be seen as the limiting case of an arbitrary large number of tanks of arbitrary small volumes) [START_REF] Van't Riet | Basic bioreactor design[END_REF][START_REF] De Gooijer | Bioreactors in series: an overview of design procedures and practical applications[END_REF][START_REF] Nelson | Evaluating the performance of a cascade of two bioreactors[END_REF][START_REF] De Gooijer | Optimal design of a series of continuous stirred tank reactors containing immobilised growing cells[END_REF][START_REF] Doran | Design of mixing systems for plant cell suspensions in stirred reactors[END_REF][START_REF] Dramé | Multiple steady state profiles in interconnected biological systems[END_REF][START_REF] Dramé | Multiple stable equilibrium profiles in tubular bioreactors[END_REF]. Sometimes a radial diffusion is also considered in plug-flow reactors [START_REF] Grady | Biological Wastewater Treatment, 3nd Edition[END_REF], but surprisingly, configurations of tanks in parallel have been much less investigated [START_REF] Levenspiel | Chemical reaction engineering[END_REF]. One can argue that knowing input rates and volumes of tanks in parallel, there dynamical characteristics can be studied separately, and there is no need of devoting a specific study for these configurations. This is no longer the case if one considers a passive communication between the tanks, through a membrane for instance. In saturated soils or wetlands, a spatial structure could be simply represented by separated domains with diffusive communication. This consideration is similar to patches models or islands models, commonly used in ecology [START_REF] Macarthur | The Theory of Island Biogeography[END_REF][START_REF] Hanski | Metapopulation ecology[END_REF], or lattice differential equations [START_REF] Scheel | Lattice differential equations embedded into reaction-diffusion systems[END_REF]. For instance, a recent investigation studies the influence of such structures on a consumer/resource model [START_REF] Gravel | Source and sink dynamics in metaecosystems[END_REF]. Consumer/resource models in ecology are similar to chemostat models, apart the source terms that are modeled as constant intakes of nutrient, instead of dilution rates that one rather met in liquid media.

In this paper, we propose to bring new insight on parallel configurations of chemostats with communication, in a spirit different than the one usually taken in bioprocesses design. One usually chooses a target for the output concentration of substrate, and looks for minimizing the total volume, or equivalently the residence time, among all the configurations that provide the same desired output at steady state. Here, we fix both the total hydric volume and the input flow and study the input-output map at steady-state, investigating the role of the spatial structure on the performances of the system. The performance is here measured by the level of substrate that is degraded by the system, and collected at the output. We draw precise comparisons between the three configurations: perfectly mixed, serial and parallel (with diffusion rate) with the same total hydric volume and flow rate. This set of configurations is far to be exhaustive, being limited to two compartments only, but it is a first attempt to grasp this input-output map of a structured chemostat, and study how a spatial structure can modify this map, and what are the key parameters. We believe that this study is of interest for the modelling of ecosystems such as saturated soils for which it is not easy to know the spatial structure, and where one has only access to input-output observations of the substrate degradation.

The paper is organized as follows. In Section 2, we present the three configurations under investigation and give the equations of the models. The main part of the paper is devoted to the analysis of the steady states, given in Section 3. The proofs of the global stability of the equilibriums are postponed to the Appendix, for lightening the presentation. Finally, discussion and numerical simulations are given in Section 4.

The models

The flow rate is labeled Q and V is the total capacity of the system. The three simple patterns we analyze are depicted on Figure 1: 2

-one single compartment of volume V -two compartments of volume V 1 , V 2 such that V = V 1 + V 2 connected in cascade -two compartments of volume V 1 , V 2 such that V = V 1 + V 2 in parallel with a diffusive connection. S = S + (1-) S α Q Q (1-α) S in α out α 2 1 out S = S 2 out S = S Q Q S in d + S in rV (1-r)V rV V Q Q Q Q (1-r)V
We recall the dynamical equations of resource (nutrient) and biomass concentrations, respectively denoted by S i and X i in a compartment i of volume V i fed from a compartment i -with a flow rate Q i and connected by diffusion rate d to a compartment i d (see Figure 2).

Ṡi = - µ(s i ) y X i + Q i V i (S i--S i ) + d V i (S i d -S i ) Ẋi = µ(s i )X i + Q i V i (X i--X i ) + d V i (X i d -X i ) d i Q Q i i i d i - Figure 2: Possible interconnections of a compartment.
For sake of simplicity of the analytical analysis, we assume that the growth function µ(•) is a linear function of the resource concentration:

µ(S) = mS
In Section 4, we shall consider Monod growth function and show that the qualitative results of our study are not changed. The yield coefficient y of the bio-conversion is kept equal to one (this is always possible by choosing the unit measuring the biomass). It is convenient to write dimensionless concentrations: for each concentration C i in the compartment i (C i can denote S i or X i ), we define

c i = m V Q C i , and
r i = V i V .
We shall also consider that the time t is measured in units such that Q = V . Finally, we assume that the input concentration S in is large enough to avoid the (trivial) wash-out equilibrium to be the only steady-state in each compartment.

3 Steady-state analysis of the three configurations

Configuration with one compartment

The dynamical equations of the configuration with a single compartment are

ṡ = -sx + s in -s ẋ = sx -x
The non-trivial equilibrium is (1, s in -1) under the condition s in > 1. Then, one has

s ⋆ out = 1 .
Remark. This is a well known property from the theory of the chemostat that the output concentration at steady state is independent of the input concentration, provided this latter to be large enough (i.e. s in ≥ 1).

Serial connection of two compartments

The dynamical equations of the model with two compartments in series (see Figure 1), assuming r to be different to 0 and 1, are

                   ṡ1 = -s 1 x 1 + 1 r (s in -s 1 ) ẋ1 = s 1 x 1 - 1 r x 1 ṡ2 = -s 2 x 2 + 1 1 -r (s 1 -s 2 ) ẋ2 = s 2 x 2 + 1 1 -r (x 1 -x 2 ) (1) with r = V 1 /V .
Proposition 1. When s in > 1/r, there exists an unique equilibrium (s

⋆ 1 , x ⋆ 1 , s ⋆ 2 , x ⋆ 2 )
of (1) on the positive orthant. One has necessarily s ⋆ 1 = 1/r and s ⋆ 2 < min(1/r, 1/(1 -r)). Furthermore, one has

s ⋆ out < 1 ⇐⇒ s in > 1 + 1/r .
Proof. One can readily check that there exists a non-trivial equilibrium (1/r, s in -1/r) for the first compartment exactly when s in > 1/r. Furthermore, this equilibrium is unique. Then, any equilibrium for the overall system (1) has to be (s ⋆ 2 , s in -s ⋆ 2 ) for the second compartment, with s ⋆ 2 solution of the equation

s 2 (s in -s 2 ) = 1 1 -r (1/r -s 2 ) ( 2
)
with s ⋆ 2 < 1/r. One can easily verify that there exists a unique s ⋆ 2 solution of ( 2) on (0, 1/r). Graphically, s ⋆ 2 is the abscissa of the intersection of the graphs (see Figure 3) of the polynomial function

φ(s 2 ) = s 2 (s in -s 2 )
and the affine function

l(s 2 ) = 1 1 -r (1/r -s 2 ) . 1/r s in 1/r-s 1-r in s 2 2 2 s* 2 s (s -s ) 2 Figure 3: Graphical determination of s ⋆ 2 .
Remark that

s in > 1/r implies the inequality φ(1/(1 -r)) > l(1/(1 -r)), from which one deduces s ⋆ 2 < 1/(1-r). Finally one can compare s out = s ⋆
2 with the value obtained in the configuration of one compartment:

s ⋆ out < 1 ⇐⇒ φ(1) > l(1) ⇐⇒ s in > 1 + 1/r .

Parallel interconnection of two compartments

The dynamical equations of the model with two compartments in parallel and diffusion (see Figure 1), assuming r to be different to 0 and 1, are the following

                     ṡ1 = -s 1 x 1 + α r (s in -s 1 ) + d r (s 2 -s 1 ) ẋ1 = s 1 x 1 - α r x 1 + d r (x 2 -x 1 ) ṡ2 = -s 2 x 2 + 1 -α 1 -r (s in -s 2 ) + d 1 -r (s 1 -s 2 ) ẋ2 = s 2 x 2 - 1 -α 1 -r x 2 + d 1 -r (x 1 -x 2 ) (3)
where the output concentration s out is given by

s out = αs 1 + (1 -α)s 2 .
The wash-out in both compartments corresponds to the trivial equilibrium (s in , 0, s in , 0), that leads to the trivial steady-state s ⋆ out = s in . For convenience, we posit

α 1 = α r , α 2 = 1 -α 1 -r ,
and assume, without any loss of generality that one has α 2 ≥ α 1 (if it is not the case one can just exchange indexes 1 and 2).

Remark. One has necessarily α 2 ≥ 1 and α 1 ≤ 1.

When d = 0 (no diffusion), the equilibrium of the system can be determined independently in the two compartments as simple chemostats. In this case, there is an unique globally stable equilibrium (s

* 1 , s in - s ⋆ 1 , s * 2 , s in -s ⋆ 2 ) in the non-negative orthant, where s * i = min(α i , s in ) (i = 1, 2). When d > 0, we define the functions φ 2 (s 1 ) = s 1 + r d (s in -s 1 )(s 1 -α 1 ) , φ 1 (s 2 ) = s 2 + 1 -r d (s in -s 2 )(s 2 -α 2 ) , and
g(s 1 ) = φ 1 (φ 2 (s 1 )) -s 1 .
Proposition 2. When s in > 1 and d > 0, there exists a unique equilibrium (s

⋆ 1 , x ⋆ 1 , s ⋆ 2 , x ⋆ 2 ) of (3) in the positive orthant, where (s ⋆ 1 , s ⋆ 2 )
is the unique solution of the system

s ⋆ 2 = φ 2 (s ⋆ 1 ) and s ⋆ 1 = φ 1 (s ⋆ 2 ) , (4)
on the domain (0, s in ) × (0, s in ), with

x ⋆ i = s in -s ⋆ i (i = 1, 2). Furthermore, s ⋆ 1 = s ⋆ 2 = 1 when α 2 = α 1 and α 1 < s ⋆ 1 < s ⋆ 2 < min(α 2 , s in ) (5)
when α 2 > α 1 .

Proof. At equilibrium, one has

r( ṡ1 + ẋ1 ) + (1 -r)( ṡ2 + ẋ2 ) = 0, r( ṡ1 + ẋ1 ) = 0,
which amounts to write, from equations (3)

α(s in -s ⋆ 1 -x ⋆ 1 ) + (1 -α)(s in -s ⋆ 2 -x ⋆ 2 ) = 0, α(s in -s ⋆ 1 -x ⋆ 1 ) + d(s ⋆ 2 + x ⋆ 2 -s ⋆ 1 -x ⋆ 1 ) = 0, or equivalently α 1 -α α + d -d M s in -s ⋆ 1 -x ⋆ 1 s in -s ⋆ 12 -x ⋆ 2 = 0 0
One has det(M ) = α 2 -α -d ≤ -d < 0 and deduces the property

s ⋆ 1 + x ⋆ 1 = s ⋆ 2 + x ⋆ 2 = s in .
Consequently, an equilibrium in the positive orthant has to fulfill 3) at equilibrium, one obtains the equations

s ⋆ i ∈ [0, s in ] for i = 1, 2. Replacing x ⋆ i by s in -s ⋆ i in equations (
d(s ⋆ 2 -s ⋆ 1 ) = r(s in -s ⋆ 1 )(s ⋆ 1 -α 1 ) d(s ⋆ 1 -s ⋆ 2 ) = (1 -r)(s in -s ⋆ 2 )(s ⋆ 2 -α 2 ) (6)
which amounts to write that (s ⋆ 1 , s ⋆ 2 ) is solution of the system (4) (see Figure 4) or equivalently s ⋆ 1 is a zero of the function g(•). When [START_REF] Dramé | Multiple steady state profiles in interconnected biological systems[END_REF]. When α 2 > α 1 , one has necessarily α 1 < 1 and the condition s in > 1 implies g(α 1 ) < 0. We distinguish now two cases:

α 2 = α 1 = 1, one can check that s ⋆ 1 = s ⋆ 2 = 1 < s in is solution of
Case α 2 < s in . If φ 2 (α 2 ) ≤ s in , notice that one has φ 2 (α 2 ) > α 2 and then g(α 2 ) > 0. If φ 2 (α 2 ) > s in , notice that φ 2 (α 1 ) = α 1 < s in and by the Mean Value Theorem, there exists s2 ∈ (α 1 , α 2 ) such that φ 2 (s 2 ) = s in which implies g(s 2 ) = s in -s2 > 0.
In both cases, one deduces by the Mean Value Theorem the existence of s

⋆ 1 ∈ (α 1 , α 2 ) such that g(s ⋆ 1 ) = 0. Case α 2 ≥ s in . One has g(s in ) = 0 with g ′ (s in ) = r(1 -r) d 2 (α 1 -s in )(α 2 -s in ) + 1 -s in d < 0 .
Rolle and Mean Value Theorems allow to conclude the existence of s ⋆ 1 ∈ (α 1 , s in ) such that g(s ⋆ 1 ) = 0.

In any case, we obtain the existence of (s ⋆ 1 , s ⋆ 2 ) solution of (4) with

s ⋆ 1 ∈ (α 1 , min(α 2 , s in )), that implies s ⋆ 2 = φ 2 (s ⋆ 1 ) > s ⋆ 1 . But then s ⋆ 1 = φ 1 (s ⋆ 2 ) < s ⋆ 2 implies s ⋆ 2 < min(α 2 , s in ).
Thus, the inequalities (5) are fulfilled.

Finally, notice that functions φ 1 (•), φ 2 (•) are both strictly concave, and steady states (s ⋆ 1 , s ⋆ 2 ) are intersections of G 1 , the graph of the function φ 1 (•), and G 2 the symmetric of the graph of φ 2 (•) with respect to the first diagonal. Consequently, if (s ⋆ 1 , s ⋆ 2 ) is a steady state different from (s in , s in ), G 1 and G 2 are respectively above and below the line segment (s ⋆ 1 , s ⋆ 2 ) -(s in , s in ). We conclude that there exists at most one non-trivial equilibrium.

Corollary 2. When s in > 1 and d > 0, the value s ⋆ 1 of the non trivial equilibrium is the unique zero of the function g(•) on (α 1 , min(α 2 , s in )). Furthermore, one has g ′ (s ⋆ 1 ) > 0.

Proof. When α 1 = α 2 , one has s ⋆ 1 = s ⋆ 2 = 1 and one can easily check

g ′ (s ⋆ 1 ) = 1 + r d (s in -1) 1 + 1 -r d (s in -1) -1 > 0 .
When α 2 > α 1 , one has g(α 1 ) < 0 and we recall from the proof of former Proposition that s ⋆ 1 is the unique zero of g(•) on (α 1 , min(α 2 , s in )). We conclude that g is non decreasing at s ⋆ 1 . Notice that φ 1 and φ 2 are concave functions and that

φ ′ 1 (φ 2 (s ⋆ 1 ))) = 1 + 1 -r d (s in + α 2 -2s ⋆ 2 ) > 0 implies g ′′ (s ⋆ 1 ) = φ ′′ 1 (φ 2 (s ⋆ 1 )). [φ ′ 2 (s ⋆ 1 )] 2 + φ ′ 1 (φ 2 (s ⋆ 1 )).φ ′′ 2 (s ⋆ 1 ) < 0 We deduce that g ′ (s ⋆ 1
) cannot be equal to zero, and consequently one has g ′ (s ⋆ 1 ) > 0.

The global stability of the non-trivial equilibrium is proved in the Appendix.

Proposition 2 defines properly the map d → s ⋆ out = αs ⋆ 1 + (1 -α)s ⋆ 2 for the unique non-trivial steady-state, that we aim at studying as a function of d. Accordingly to Proposition 2, s ⋆ out is equal to one for any value of the parameter d in the non-generic case α 2 = α 1 . We shall focus on the case α 2 = α 1 (and without loss of generality we shall consider α 2 > α 1 ). We start by the two extreme situations: no diffusion and infinite diffusion.

Lemma 1. For the non trivial equilibrium, one has

s ⋆ out (0) ≥ 1 ⇐⇒ s in ≥ s 0 in = r -α 2 r(1 -α) with s 0 in ∈ (1, 2).
Proof. Under the assumptions s in > 1 and α 2 ≥ α 1 , we distinguish two cases when

d = 0. If s in ≥ α 2 , one has s ⋆ 1 = α 1 and s ⋆ 2 = α 2 .
Then, one can write

s ⋆ out = α 2 r + (1 -α) 2 1 -r = 1 + (α -r) 2 r(1 -r) ≥ 1 . If s in < α 2 , on has s ⋆ 1 = α 1 , s ⋆ 2 = s in and s ⋆ out ≥ 1 ⇐⇒ s in ≥ 1 -αα 1 1 -α = s 0 in .
(recall that assuming α 2 ≥ α 1 imposes to have α < 1, and s 0 in is well defined). Notice that the number s 0 in is necessarily larger than one because α 1 ≤ 1, and one has also

α 2 -s 0 in = (r -α) 2 r(1 -r)(1 -α) ≥ 0.
Consequently one concludes that s ⋆ out ≥ 1 exactly when s in ≥ s 0 in . Finally, remark that one has

s 0 in = r -α 2 r(1 -α) = 1 - (α -r) 2 r(r -α) + r -α 1 -α < 2 .
Lemma 2. For s in > 1, the non trivial equilibrium fulfill

lim d→+∞ s ⋆ 1 (d) = lim d→+∞ s ⋆ 2 (d) = lim d→+∞ s ⋆ out (d) = 1 .
Proof. For any d > 0, Proposition guarantees the existence of a unique non trivial equilibrium (s ⋆ 1 , s ⋆ 2 ) ∈ (0, s in ) × (0, s in ) that is solution of (6). When d is arbitrary large, one obtains from (6)

lim d→+∞ s ⋆ 1 (d) -s ⋆ 2 (d) = 0 .
From equations ( 6), one deduces also the following equality valid for any d

r(s in -s ⋆ 1 )(s ⋆ 1 -α 1 ) + (1 -r)(s in -s ⋆ 2 )(s ⋆ -α 2 ) = 0 ,
that can rewritten, taking into account the equality

rα 1 + (1 -r)α 2 = 1: (s in -s ⋆ 1 )(s ⋆ 1 -1) = (1 -r)(s ⋆ 1 -s ⋆ 2 )(s in + α 2 -s ⋆ 1 -s ⋆ 2 ) .
Consequently, one has

lim d→+∞ s ⋆ 1 (d) = lim d→+∞ s ⋆ 2 (d) = 1 or lim d→+∞ s ⋆ 1 (d) = lim d→+∞ s ⋆ 2 (d) = s in .
If α 2 < s in , the property s ⋆ 1 < α 2 valid for any d > 0 implies that s ⋆ 1 cannot converges to s in . If α 2 ≥ s in and lim s ⋆ 1 = lim s ⋆ 2 = s in , there exists d such that rs ⋆ 1 + (1 -r)s ⋆ 2 > (s in + 1)/2. Then, one has

g ′ (s ⋆ 1 ) = r(1 -r) d 2 (s in + α 1 -2s ⋆ 1 )(s in + α 2 -2s ⋆ 2 ) + s in + 1 -2(rs ⋆ 1 + (1 -r)s 2 ⋆ ) d < 0
that contradicts Corollary 2. Finally, one has lim s ⋆ 1 = lim s ⋆ 2 = 1 and consequently lim s ⋆ out = 1.

We present now our main result concerning properties of the map d → s ⋆ out (d) defined at the non-trivial steady-state.

Proposition 3. Assume α 2 > α 1 .
-When s in ≥ 2, the map d → s ⋆ out (d) the non trivial equilibrium) is decreasing and s ⋆ out (d) > 1 for any d ≥ 0.

-When s in < 2, the map d → s ⋆ out (d) (for the non trivial equilibrium) admits a minimum in d ⋆ < +∞, that is strictly less than one. Furthermore, one has

s in > s in = 2α 1 α 2 α 1 + α 2 =⇒ d ⋆ > 0 with s in < min(2, α 2 ).
Proof. Let differentiate with respect to d the equations ( 6) at steady state:

(s ⋆ 2 -s ⋆ 1 ) + d (∂ d s ⋆ 2 -∂ d s ⋆ 1 ) = r(s in -2s ⋆ 1 + α 1 ) A ∂ d s ⋆ 1 (s ⋆ 1 -s ⋆ 2 ) + d (∂ d s ⋆ 1 -∂ d s ⋆ 2 ) = (1 -r)(s in -2s ⋆ 2 + α 2 ) B ∂ d s ⋆ 2
that can rewritten as follows

A + d -d d -B -d Γ ∂ d s ⋆ 1 ∂ d s ⋆ 2 = (s ⋆ 2 -s ⋆ 1 )
1 1

Remark that one has

A + d = dφ ′ 2 (s ⋆ 1 ) B + d = dφ ′ 1 (s ⋆ 2 ) det(Γ) = d 2 (1 -φ ′ 1 (s ⋆ 2 )φ ′ 2 (s ⋆ 1 )) = -d 2 g ′ (s ⋆ 1 )
From the Corollary 2, one has det(Γ) < 0 and one deduces that the derivatives

∂ d s ⋆ 1 , ∂ d s ⋆ 2 are defined as follows ∂ d s ⋆ 1 = (s ⋆ 2 -s ⋆ 1 ) -B det(Γ) ∂ d s ⋆ 2 = (s ⋆ 2 -s ⋆ 1 ) A det(Γ) (7)
Notice from inequalities (5) that we obtain B > 0 and deduce ∂ d s ⋆ 1 > 0 for any d. With Lemma 2 we conclude that s ⋆ 1 (d) < 1 for any d. From equations [START_REF] Grady | Biological Wastewater Treatment, 3nd Edition[END_REF], we can write

∂ d s ⋆ out = (s ⋆ 2 -s ⋆ 1 ) αB -(1 -α)A -det(Γ) = [α 1 (s in -2s ⋆ 2 ) -α 2 (s in -2s ⋆ 1 )] σ (s ⋆ 2 -s ⋆ 1 )r(1 -r) -det(Γ)
When s in ≥ 2, one has A > 0 and then ∂ d s ⋆ 2 < 0. With Lemma 2 we conclude that s ⋆ 2 (d) > 1 for any d. Then, one obtain the inequality σ < (s in -2)(α 1 -α 2 ) ≤ 0 which proves with Lemma 2 that s ⋆ out is a decreasing function of d that converges to one.

When s in < 2, we write

σ = (s in -2)(α 1 -α 2 ) + 2(α 1 (1 -s ⋆ 2 ) -α 2 (1 -s ⋆ 1 ))
As s ⋆ 1 and s ⋆ 2 tend to one when d takes arbitrary large values, we conclude that there exists d < +∞ such that σ > 0 for any d > d and consequently s ⋆ out is smaller than one and increasing for d > d. We conclude that the map d → s ⋆ out (d) admits a minimum, say at d ⋆ < +∞, that is strictly less than one.

When d = 0, one has s ⋆ 1 = α 1 and s ⋆ 2 = α 2 if s in ≥ α 2 . Then, one obtains σ = s in (α 1 -α 2 ) < 0. So the map d → ∂ d s ⋆ out (d) is decreasing at d = 0 and consequently d ⋆ > 0. When d = 0 with s in < α 2 , one has s ⋆ 2 = s in and then σ = 2α 1 α 2 -s in (α 1 + α 2 )
, for which we conclude

σ < 0 ⇐⇒ s in > 2α 1 α 2 α 1 + α 2 = s in .
Remark that this case is feasible because of the inequality 2α 1 α 2 < min(2, α 2 )(α 1 + α 2 ). We conclude that for s in larger than this last value, d ⋆ is necessarily strictly positive.

Remark. The particular case α = 0 corresponds to a configuration of a perfectly mixed tank of volume (1 -r)V connected to a dead-zone of volume rV . This is a way to approximate a non well-mixed tank or segregated bioreactors of total volume V , estimating the fraction of the volume occupied by the highly agitated area.

Numerical computation and discussion

Propositions 1 and 3 reveal the existence of a threshold on the value of the input concentration s in (equal to 2 for our choice of the parameters units) that reverses the performances of the serial and parallel configurations in terms of s ⋆ out , compared to the single tank case (for which s ⋆ out = 1):

-for s in > 2, there exist serial configurations such that s ⋆ out < 1 for r large enough (i.e. the first tank has to be large enough), but any parallel configuration produces s ⋆ out > 1,

-for s in < 2, there exists parallel configurations such that s ⋆ out < 1, while any serial configuration has s ⋆ out > 1. There exists another threshold s 0 in ∈ (1, 2) such that configurations with s ⋆ out < 1 require to have d large enough when s in > s 0 in (cf Lemma 1).

Furthermore, the best performance of the parallel configuration is obtained -for arbitrary large values of d when s in > 2 ,

-for a finite positive d ⋆ when s in ∈ (s in , 2) (where the expression of s in is given in Proposition 3).

For the serial configuration, the graph of the function s ⋆ out is plotted as function of r ∈ [1/s in , 1] on Figure 5 for different values of the input concentration s in . For the parallel interconnection, we depict on Figure 6 the two kind of configurations that occur, depending on whether the number s in is larger than one or not. The analytic analysis of Section 3 has been conducted under the assumption of the linearity of the function µ(•). It is often in microbiology that the growth rate µ(•) presents a concavity, as described by the usual Monod (or Michaelis-Menten) function. We have computed numerically the same curves s ⋆ out (•) than Figures 5 and6, considering the Monod function µ(S) = 6S 5 + S instead of the linear function (see Figure 7). This function has been chosen to fulfill s ⋆ out = 1 for the single tank configuration, guaranteeing the same steady state than the linear growth for this configuration.

The values of the parameters are given on the table below

On Figures 8 and9, we observe that the concavity of the growth function does not change qualitatively the theoretical results and the existence of threshold for s in that favourites one of the configuration.

We notice on all the figures that the yield is better for the Monod function in the parallel configuration and worst for the serial one. This implies that the threshold on s in , that was determined to be equal to 2 for the linear case is higher when the growth function is concave.

Remarks. The serial configuration for the limiting value r = 1 is equivalent to a single tank. This explains why all the curves on Figures 5 and8 For the parallel configuration with α = 0.1 and r = 0.9 one has α 2 = 9. This implies that for the limiting value d = 0 the only equilibrium in the second tank is the wash-out when s in < 9. This is not the case for the first tank but the flow rate αQ being small, the output s ⋆ out remains closed to s in in any case, as one can see on Figures 6 and9 for small values of the parameter d.

Conclusion

Given a flow rate and the total volume of a chemostat system, this study shows the existence of a threshold on the value of the input concentration s in such that above and below this threshold, serial and the parallel configurations are respectively the best ones with respect to the criterion of minimizing the output concentration s ⋆ out at steady state. For the parallel scheme, the best performances are obtained for a precise value of the diffusion parameter that is proved to be positive when s in is not too small. This study concerns also dead-zone configurations, as particular cases of the parallel configurations.

Whatever are the data of the problem, there always exists a configuration that is better than a single perfectly mixed tank. We have shown that the non-trivial steady states are unique and globally exponentially stable under the assumption of a linear increasing growth rate.

Finally, this study reveals the role of the structure of the space on the performances of simple ecosystems or bioprocesses. The possibly non-monotonic influence of the diffusion parameter on the output steady state is not intuitive, and leave further investigations open for understanding or taking benefit of this property for natural ecosystems (such as saturated soils or wetlands) as well as for bioprocesses (such as waste-water treatments). This result can be also of interest for reverse engineering when deciding which among serial or parallel configurations is better fit for the modeling of chemostat-like ecosystems, providing that one has an estimation of the hydric capacity of the system. 3). We consider the 2-dimensional vector z of variables z i = s in -x i -s i (i = 1, 2) whose dynamics are respectively for the serial and parallel configurations

ż = A s z = -1 r 0 1 1-r -1 1-r z ż = A p z = -α 1 -d r d r d 1-r -α 2 -d 1-r z
Notice that matrices A s and A p are Hurwitz :

tr(A s ) = - 1 r - 1 1 -r < 0 , det(A s ) = 1 r(1 -r) > 0 tr(A p ) = -α 1 -α 2 - d r - d 1 -r < 0 , det(A p ) = α 1 α 2 + d r(1 -r) > 0
So z converges exponentially toward 0 for both systems, which implies that dynamics (1) and ( 3) are dissipative, in the sense that any solution of (1) (3) in D converge exponentially to the compact set

K = {(s 1 , x 1 , s 2 , x 2 ) ∈ D s.t. x 1 + s 1 = s in and x 2 + s 2 = s in }.
We recall a result from [19, Theorem 1.8] that shall be useful in the following.

Theorem 4. Let Φ be an asymptotically autonomous semi-flow with limit semi-flow Θ, and let the orbit O Φ (τ, ξ) have compact closure. Then the ω-limit set ω Φ (τ, ξ) is non-empty, compact, connected, invariant and chain-recurrent by the semi-flow Θ and attracts Φ(t, τ, ξ) when t → ∞.

The serial configuration

Proposition 5. Under the condition s in > 1/r, any trajectory of (1) with initial condition in D such that (s 1 (0), x 1 (0)) = (s in , 0) converges exponentially to the unique non-trivial steady-state (s

⋆ 1 , x ⋆ 1 , s ⋆ 2 , x ⋆
2 ) given by Proposition 1.

Proof. Dynamics (1) has a cascade structure. It is straightforward to check that the solutions of the (s 1 , x 1 ) sub-system converges asymptotically towards the non-trivial equilibrium (1/r, s in -1/r) from any initial condition away from the wash-out equilibrium (s in , 0). From the convergence of z 2 toward 0, we deduce that the s 2 variable has to converge to the bounded interval [0, s in ] and that its dynamics can be written as a scalar non autonomous differential equation:

ṡ2 = -s 2 (s in -s 2 -z 2 (t)) + 1 1 -r (s 1 (t) -s 2 ) (8)
This last dynamics has the property to be asymptotically autonomous with the limiting differential equation:

ṡ2 = f (s 2 ) = -s 2 (s in -s 2 ) + 1 1 -r (1/r -s 2 ) (9)
Statement of Proposition 1 implies that this last scalar dynamics has a unique equilibrium s ⋆ 2 that belongs to [0, s in ]. Furthermore, one has f (0) > 0 and f (s in ) < 0. Consequently any solution of ( 9) in [0, s in ] converges asymptotically to s ⋆ 2 . Then applying Theorem 4, we conclude that any bounded solution of (8) converges to s ⋆ 2 . Finally any solutions of the (s 2 , x 2 ) sub-system converges asymptotically to (s ⋆ 2 , s in -s ⋆ 2 ).

The Jacobian matrix of dynamics (1) at the non-trivial equilibrium (s ⋆ 1 , x ⋆ 1 , s ⋆ 2 , x ⋆ 2 ) is of the following form in (z 1 , z 2 , s 1 , s 2 ) coordinates

A s 0 ⋆ J ⋆ with J ⋆ = s ⋆ 1 -s in 0 1 1-r 2s ⋆ 2 -1 1-r -s in
Recall that A s is Hurwitz. The eigenvalues of J ⋆ are -(s in -s ⋆ 1 ), -(s in -s ⋆ 2 ) -(1/(1 -r) -s ⋆ 2 ) that are both negative numbers, accordingly to Proposition 1. The exponential stability of the non-trivial equilibrium is thus proved.

The parallel configuration

Proposition 6. When s in > 1 and d > 0, any trajectory of (3) with initial condition in D such that x 1 (0) > 0 and x 2 (0) > 0 converges exponentially to the unique non-trivial steady-state (s ⋆ 1 , x ⋆ 1 , s ⋆ 2 , x ⋆ 2 ) given by Proposition 2.

Proof. Considering the time vector z(•), the (s 1 , s 2 ) sub-system of dynamics (3) can be written as solution of a non-autonomous planar dynamics ṡ1 = s 1 (z 1 (t) + s 1 -s in ) + α 1 (s in -s 1 ) + d r (s 2 -s 1 ) ṡ1 = s 2 (z 2 (t) + s 2 -s in ) + α 2 (s in -s 2 ) + d 1-r (s 1 -s 2 )

We know that z converges to 0 and consequently the vector S of variables s 1 , s 2 converges to the set S = [0, s in ] × [0, s in ]. We study now the limiting autonomous dynamics ṡ1 = (s in -s 1 )(α 1 -s 1 ) + d r (s 2 -s 1 ) ṡ2 = (s in -s 2 )(α 2 -s 2 ) + d 1-r (s 1 -s 2 ) [START_REF] Harmand | Optimal design of two interconnected enzymatic reactors[END_REF] on the domain S. Let B be the boundary {s 1 = s in } ∪ {s 2 = s in }. On the domain S \ B, we consider the vector σ of variables σ i = log(s in -s i ), whose dynamics can be written as follows σ = F (σ) = -α 1 + s in -e σ1 -d r (1 -e σ2-σ1 ) -α 2 + s in -e σ2 -d 1-r (1 -e σ1-σ2 ) (

One can easily compute div(F ) = -e σ1 -e σ2 -d r e σ2-σ1 -d 1 -r e σ1-σ2 < 0

From Poincaré-Bendixon theorem and Dulac criterion, we conclude that bounded trajectories of (12) cannot have limit cycle or closed path and necessarily converge to an equilibrium point. Consequently, any trajectory of [START_REF] Harmand | Optimal design of two interconnected enzymatic reactors[END_REF] in S either converges to the rest point S ⋆ = (s ⋆ 1 , s ⋆ 2 ) or approaches the boundary B. Notice that one has s i = s in , s j < s in ⇒ ṡi < 0 (i = j)

Figure 1 :

 1 Figure 1: The set of configurations under investigation.

1 Figure 4 :

 14 Figure 4: Graphical determination of steady states (when α 1 < α 2 < s in ).

Figure 5 :

 5 Figure 5: Comparison of s ⋆ out for the serial configuration.

Figure 6 :

 6 Figure 6: Comparison of s ⋆ out for the parallel configuration (s in > 1 on the left and s in < 1 on the right) .

Figure 7 :

 7 Figure 7: Linear and Monod growth functions.

Figure 9 :

 9 Figure 9: Comparison of s ⋆ out for the parallel configuration (s in > 1 on the left and s in < 1 on the right) .

Acknowledgments.

  This work has been achieved within the VITELBIO (VIRtual TELluric BIOreactors) program, sponsored by INRA and INRIA. The authors are grateful to this support. The work is also part of the PhD thesis of the first author.

6

  Appendix: global exponential stability of the non-trivial equilibrium First, one can easily check that the domain D = IR4 + is invariant by the dynamics (1) and (

  coincide for this value of the parameter r.

							S out
							2.6
							2.4
							2.2
							2.0
							1.8
							1.6
							1.4
							1.2
							1.0
							0.8
							0.6										r
							0.2		0.3	0.4	0.5		0.6	0.7	0.8	0.9	1.0
	Figure 8: Comparison of s ⋆ out for Monod (dashed) and linear (plain) for the serial configuration.
	S out											S out
	1.25												2.0
	1.20
													1.8
	1.15
													1.6
	1.10
	1.05												1.4
	1.00												1.2
	0.95
													1.0
	0.90
													0.8
	0.85
	0.80	0.0	0.2	0.4	0.6	0.8	1.0	1.2	1.4	1.6	1.8	d	0.6	0	1	2	3	4	5	6	d

So the only possibility for approaching B is to converge to the other rest point S 0 = (s in , s in ). This shows that the only non-empty, closed, connected, invariant and chain recurrent subsets of S are the singletons {S ⋆ } and {S 0 }.

Applying Theorem 4 we conclude that any trajectory of [START_REF] Harmand | Optimal design of two interconnected bioreactors-some new results[END_REF], issued from initial condition of dynamics (3) in D, converges asymptotically to S ⋆ or S 0 . Consider now any initial condition with x 1 (0) > 0 and x 2 (0) > 0. We show that the solution (s 1 (•), s 2 (•)) of (3) cannot converge to S 0 . If it is the case, there exists T < +∞ such that one has s 1 (t) > α 1 and rs 1 (t) + (1 -r)s 2 (t) > 1 for any t ≥ T under the assumption s in > 1. Let us consider the function

(see Figure 10) and v(t) = V (x 1 (t), x 2 (t)) that is positive and tends to 0 when t tends to +∞

We conclude that the function t → v(t) is non-decreasing for t ≥ T and consequently cannot converge to zero, thus a contradiction.

The Jacobian matrix of dynamics (3) at the non-trivial equilibrium (s

Recall that A p is Hurwitz. One has

Then the inequality φ ′ 2 (s ⋆ 1 ) > 0 is necessarily satisfied. Finally, we have shown det(J ⋆ ) > 0 and tr(J ⋆ ) < 0, that guarantee the exponential stability of the non-trivial equilibrium (s

Remark. The wash-out equilibrium (s in , 0, s in , 0) is not necessarily hyperbolic. This explains why we cannot use the Convergence Theorem for asymptotically autonomous dynamics given in Appendix F of [START_REF] Smith | The theory of chemostat, dynamics of microbial competition[END_REF].


A E Koudou

P Vallois
email: pierre.vallois@iecn.u-nancy.fr

Independence properties of the Matsumoto-Yor type

Keywords: Gamma distribution, generalized inverse Gaussian distribution, Matsumoto-Yor property, Kummer distribution

We define Letac-Wesolowski-Matsumoto-Yor (LWMY) functions as decreasing functions from (0, ∞) onto (0, ∞) with the following property: there exist independent, positive random variables X and Y such that the variables f (X + Y ) and f (X)f (X + Y ) are independent. We prove that, under additional assumptions, there are essentially four such functions. The first one is f (x) = 1/x. In this case, referred to in the literature as the Matsumoto-Yor property, the law of X is generalized inverse Gaussian while Y is gamma-distributed. In the three other cases, the associated densities are provided. As a consequence, we obtain a new relation of convolution involving gamma distributions and Kummer distributions of type 2.

Introduction

A lot of papers have been devoted to the generalized inverse Gaussian (GIG) distributions since their definition by [START_REF] Good | The population frequencies of species and the estimation of population parameters[END_REF](see for instance [START_REF] Barndorff-Nielsen | Infinite divisibility of the Hyperbolic and generalized inverse Gaussian distribution[END_REF], [START_REF] Letac | A characterization of the generalized inverse Gaussian distribution by continued fractions[END_REF], [START_REF] Vallois | Sur le passage de certaines marches aléatoires planes au-dessus d'une hyperbole quilatère[END_REF], [START_REF] Vallois | La loi Gaussienne inverse généralisée comme premier ou dernier temps de passage de diffusion[END_REF]). The GIG distribution with parameters µ ∈ R, a, b > 0 is the probability measure :

GIG(µ, a, b)(dx) = b a µ x µ-1 2K µ (ab) e -1 2 (a 2 x -1 +b 2 x) 1 (0,∞) (x)dx (0.1)
where K µ is the classical McDonald special function. GIG distributions can also be defined as probability measures on the set of positive definite matrices [START_REF] Letac | An independence property for the product of GIG and gamma laws[END_REF]. Here, we only deal with one-dimensional GIG distributions.

1) Let us stress the close links between GIG, gamma distributions and the function f 0 (x) = 1/x (x > 0).

a) The family of GIG distributions is invariant under f 0 : we can easily deduce from (0.1) that the image of GIG(µ, a, b) by f 0 is GIG(-µ, b, a). where γ(µ, b 2 /2)(dx) = b 2µ 2 µ Γ(µ) x µ-1 exp -b 2 2 x 1 (0,∞) (x)dx. Therefore if X > 0 and Y ∼ γ(λ, a 2 /2) are independent r.v.'s and X ∼ GIG(-λ, a, a) then X

= f 0 (X + Y ). (0.3) [START_REF] Letac | A characterization of the generalized inverse Gaussian distribution by continued fractions[END_REF] proved that (0.3) characterizes GIG distributions of the type GIG(-λ, a, a). Namely, if X and Y are positive r.v.'s such that Y ∼ γ(λ, a 2 /2) and (0.3) holds then X ∼ GIG(-λ, a, a). The authors have also obtained a similar characterization of GIG(-λ, a, b) but the formulation is more complicated. c) Since relation (0.2) is only formulated in terms of probability measures, we can ask whenever it admits an almost sure realization. That consists in giving a triplet of r.v.'s X, Y , Z such that Z = X + Y , X and Y are independent, X ∼ GIG(-µ, a, b) and Y ∼ γ(µ, b 2 /2). Almost sure realizations have been given by [START_REF] Bhattacharya | Stochastic Processes with Applications[END_REF] in the case µ = 1 2 and Vallois (1991) for any µ > 0. The authors have considered a family of transient diffusions (D ν (t), t ≥ 0) on [0, ∞) which depends on a parameter ν. It has been proved that there exist x 0 > 0 and ν 0 such that D ν 0 (0) = x 0 and T 0 := inf{t ≥ 0, D ν 0 (t) = 0} ∼ GIG(-µ, a, b) L 0 := sup{t ≥ 0, D ν 0 (t + T 0 ) = 0} ∼ γ(µ, b 2 /2).

The strong Markov property at time T 0 implies that X = T 0 and Y = L 0 are independent and Z = X + Y = sup{t ≥ 0, D ν 0 (t) = 0}. In Madan, Roynette and [START_REF] Madan | Unifying Black-Scholes type formulae which involve Brownian last passage times up to a finite horizon[END_REF], the Black-Scholes formula in finance is expressed in terms of the distribution function of GIG variables (see Equations ( 25), (29), (30) in this reference).

Let (ξ 1 (t), t ≥ 0) and (ξ 2 (t), t ≥ 0) be two Lévy processes such that ξ 1 (t) ∼ γ(t, β 2 1 /2) and ξ 2 (t) ∼ γ(t, β 2 2 /2) for any t > 0, where β 1 , β 2 > 0. Recall that for i = 1, 2 the process (ξ i (t)) is a subordinator and consequently t → ξ i (t) is non decreasing. Let us define the random time

N = inf {n ≥ 0; ξ 1 (α 1 + n -1)ξ 2 (α 2 + n -1) ≤ 1 < ξ 1 (α 1 + n)ξ 2 (α 2 + n -1)} = inf {n ≥ 0; ξ 1 (α 1 + n -1) ≤ f 0 (ξ 2 (α 2 + n -1)) < ξ 1 (α 1 + n)} (0.4)
where α 1 , α 2 > 0 and by convention inf ∅ = +∞.

Then, it has been proved in [START_REF] Vallois | Sur le passage de certaines marches aléatoires planes au-dessus d'une hyperbole quilatère[END_REF], theorem on p.446) that, conditionally on {N < ∞}, the r.v.'s X := ξ 2 (α 2 + N -1), Y := ξ 2 (α 1 + N -1)ξ 2 (α 2 + N -1) (0.5) are independent and their laws are GIG(-µ, β 1 , β 2 ) and respectively γ(µ, β 2 2 /2), where µ = α 1α 2 . This yields a second almost sure realization of (0.2). Furthermore, the r.v.'s

U := 1 X + Y = ξ 2 (α 1 +N-1), V := 1 X - 1 X + Y = 1 ξ 2 (α 2 + N -1) - 1 ξ 2 (α 1 + N -1) are independent.
2) We now focus on the Matsumoto-Yor property involving the GIG and gamma distributions on the one hand and the function f 0 on the second hand . Let X and Y be two independent r.v.'s such that

X ∼ GIG(-µ, a, b), Y ∼ γ(µ, b 2 /2), (µ, a, b > 0). (0.6) Then U := 1 X + Y = f 0 (X + Y ), V := 1 X - 1 X + Y = f 0 (X) -f 0 (X + Y ) (0.7) are independent and U ∼ GIG(-µ, b, a), V ∼ γ(µ, a 2 /2). (0.8)
The case a = b was proved by [START_REF] Matsumoto | An analogue of Pitman's 2M -X theorem for exponential Wiener functional, Part II: the role of the generalized inverse Gaussian laws[END_REF] and a nice interpretation of this property via Brownian motion was given by [START_REF] Matsumoto | Interpretation via Brownian motion of some independence properties between GIG and gamma variables[END_REF]. The case µ = -1 2 of the Matsumoto-Yor property can be retrieved from an independence property established by Barndorff-Nielsen and Koudou (1998) (see [START_REF] Koudou | A link between the Matsumoto-Yor property and an independence property on trees[END_REF].

Note that the r.v.'s X and Y introduced in (0.5) satisfy the Matsumoto-Yor property. [START_REF] Letac | An independence property for the product of GIG and gamma laws[END_REF] proved that the Matsumoto-Yor property holds for any µ, a, b > 0 and characterizes the GIG distributions. More precisely, consider two independent and non-Dirac positive r.v.'s X and Y such that U and V defined by (0.7) are independent, then there exist µ, a, b > 0 such that (0.6) holds. Obviously this property is similar to the one given by [START_REF] Letac | A characterization of the generalized inverse Gaussian distribution by continued fractions[END_REF]. [START_REF] Massam | The Matsumoto-Yor property on trees[END_REF] derived a tree-version of the Matsumoto-Yor property.

3) The origin of this paper is to understand the link between the function f 0 : x → 1/x and the GIG distributions in the Matsumoto-Yor property. It is convenient to introduce the following transformation T f associated with a decreasing function f : (0, ∞) → (0, ∞) :

T f : (0, ∞) 2 → (0, ∞) 2 (x, y) → (f (x + y), f (x) -f (x + y)).
Suppose that f is one-to-one, then T f is bijective and

(T f ) -1 = T f -1 . Let (U, V ) = T f (X, Y )
where X and Y are positive r.v.'s, namely

U = f (X + Y ), V = f (X) -f (X + Y ). (0.9)
We can recover X and Y from (U, V ):

X = f -1 (U + V ), Y = f -1 (U) -f -1 (U + V ). (0.10)
Considering f = f 0 and independent r.v.'s X and Y , the Matsumoto-Yor property asserts that the r.v.'s f 0 (X + Y ) and f 0 (X)f 0 (X + Y ) are independent as soon as X ∼ GIG(-µ, a, b) and Y ∼ γ(µ, b 2 /2) for some µ, a, b > 0.

Obviously, the Matsumoto-Yor property can be reexpressed as follows: the image of the probability measure (on R 2 + ) GIG(-µ, a, b) ⊗ γ(µ, b 2 /2) by the transformation T f 0 is the probability measure GIG(-µ, b, a) ⊗ γ(µ, a 2 /2). This formulation of the Matsumoto-Yor property joined with the Letac and Wesolowski result lead us to determine the triplets (µ X , µ Y , f ) such that 1. µ X , µ Y are probability measures on (0, ∞), 2. f : (0, ∞) → (0, ∞) is bijective and decreasing, 3. if X and Y are independent r.v.'s such that X ∼ µ X and Y ∼ µ Y then the r.v.'s U and V defined by (0.9) are independent.

Unfortunately we have not been able to solve this question without restriction. Our method can be applied supposing that f is smooth and µ X and µ Y have smooth density functions (see Theorem 2.1 for details). After long and sometimes tedious calculations we prove (cf Theorem 1.2) that there are only four classes F 1 , . . . , F 4 of functions f such that T f keeps the independence property. The first class F 1 = {α/x; α > 0} corresponds to f = f 0 . In Appendix 6.1 we prove that the only possible distributions for X and Y are GIG and gamma respectively. Thus, we recover the result of Letac and Wesolowski under stronger assumptions. The proof of Letac and Wesolowski is completely different since the authors made use of Laplace transforms and a characterization of the GIG laws as the distribution of a continued fraction with gamma entries (a property close to (0.3)). We have not been able to develop a proof as elegant as theirs. The reason is that, with f = f 0 we have algebraic properties (for instance continued fractions), while these properties are lost if we start with a general function f . Then, for any f ∈ F i , 1 ≤ i ≤ 4 we have been able to give the corresponding distributions of X and Y and the related laws of U and V (see Theorem 1.4, Remark 1.5 and Theorem 1.14).

It is worth pointing out that one interesting feature of our analysis is an original characterization of the families of distributions {β α (a, b, c); a, b, α > 0, , c ∈ R} and the Kummer distributions {K (2) (a, b, c); a, c > 0, , b ∈ R} (see (1.18) and (1.33) respectively). The Kummer distributions appear as the law of some random continued fractions (see Marklov et al, 2008, p.3393 mentioning a work by [START_REF] Dyson | The dynamics of a disordered linear chain[END_REF] in the setting of random matrices).

As by-products of our study we obtain new relations of convolution. For simplicity we only detail the case of Kummer distributions. The Kummer distributions of type 2 and the gamma distributions are closely connected since:

K (2) (a, b, c) * γ(b, c) = K (2) (a + b, -b, c). (0.11)
Obviously, this relation is similar to (0.2).

Inspired by the result of [START_REF] Letac | An independence property for the product of GIG and gamma laws[END_REF] and Theorem 1.6, we can ask if a characterization of Kummer distributions could be obtained (cf Remark 1.8).

As recalled in the above item c), there are various almost sure realizations of convolution identities involving GIG distributions (see (0.2) and the convolution coming from the Matsumoto-Yor property. One interesting open question derived from our study would be to determine a r.v. Z with distribution K (2) (a + b, -b, c) which can be decomposed as the sum of two explicit independent r.v.'s X and Y such that X ∼ K (2) (a + b, -b, c) and Y ∼ γ(b, c).

The paper is organized as follows. We state our main results in Section 1. In Section 2 we give, under smoothness assumptions, a key differential equation involving f and the log densities of the independent r.v.'s X and Y such that f (X + Y ) and f (X)f (X + Y ) are independent (cf Theorem 2.1). Based on this result we prove in Section 3 (cf Theorem 3.8) that there are only four classes F 1 , . . . , F 4 of such functions f . The theorems stated in Section 2 are proved in Section 4. Some technical proofs have been postponed in Appendix 5.

Main results

Definition 1.1 Let f : (0, ∞) → (0, ∞) be a decreasing and bijective function.

Associated with f let us consider the transformation

T f : (0, ∞) 2 → (0, ∞) 2 (x, y) → (f (x + y), f (x) -f (x + y)).
(1.12)

The transformation T f is one-to-one and if f -1 is the inverse of f , then

(T f ) -1 = T f -1 . (1.13)
2. Let X and Y be two independent and positive random variables. Let us define

(U, V ) = T f (X, Y ) = (f (X + Y ), f (X) -f (X + Y )). (1.14)
f is said to be a LWMY function with respect to (X, Y ) if the random variables U and V are independent. f is said to be a LWMY function if it is a LWMY function with respect to some random vector (X, Y ).

One aim of this paper is to characterize LWMY functions. Let us introduce

f 1 (x) = 1 e x -1
, x > 0, (1.15)

g 1 (x) = f -1 1 (x) = ln 1 + x x , x > 0 (1.16)
and, for δ > 0,

f * δ (x) = log e x + δ -1 e x -1 , x > 0. (1.17) Theorem 1.2 Let f : (0, ∞) → (0, ∞) be

decreasing and bijective. Under some additional assumptions (see Theorem 2.1 and Equation (3.3) ), f is a LWMY function if and only if, either f

(x) = α x or f (x) = 1 α f 1 (βx) or f (x) = 1 α g 1 (βx) or f (x) = 1 α f * δ (βx) for some α, β, δ > 0. Remark 1.3
1. The four classes of LWMY functions are

F 1 = {α/x; α > 0}, F 2 = { 1 α f 1 (βx); α, β > 0}, F 3 = { 1 α g 1 (βx); α, β > 0}, F 4 = { 1 α f * δ (βx); α, β > 0}. 2. It is clear that if f is a LWMY function, then the functions f -1 and x → 1 α f (βx), α, β > 0 are LWMY functions. Consequently, x → 1 α f 1 (βx) is a LWMY function if and only if x → 1 α g 1 (βx) is a LWMY function. 3. Note that x → α/x and f * δ are self-reciprocal.
The first case corresponds to the one studied by [START_REF] Matsumoto | An analogue of Pitman's 2M -X theorem for exponential Wiener functional, Part II: the role of the generalized inverse Gaussian laws[END_REF] and [START_REF] Letac | An independence property for the product of GIG and gamma laws[END_REF]. The two last authors proved that if x → 1/x is a LWMY function with respect to (X, Y ), then the laws of X and Y are GIG and gamma distributions respectively (we partially recover this case in Appendix 5.1). Therefore, in the sequel we will focus on the three new cases :

either f = f 1 or f = g 1 or f = f *
δ and in each case we determine the laws of the related random variables.

The cases

f = g 1 and f = f 1 a) Consider the gamma distribution γ(λ, c)(dx) = c λ Γ(λ) x λ-1 e -cx 1 (0,∞) (x)dx, (λ, c > 0)
and the beta distribution

Beta(a, b)(dx) = Γ(a + b) Γ(a)Γ(b) x a-1 (1 -x) b-1 1 {0<x<1} dx, (a, b > 0).
Consider (see for instance Ng and Kotz, 1995, or Nagar and Gupta, 2002 and references therein) the Kummer distribution of type 2 :

K (2) (a, b, c) := α(a, b, c)x a-1 (1 + x) -a-b e -cx 1 (0,∞) (x)dx, a, c > 0, b ∈ R (1.18)
where α(a, b, c) is a normalizing constant. Associated with a couple (X, Y ) of positive r.v.'s consider

(U, V ) := T f 1 (X, Y ) = 1 e X+Y -1 , 1 e X -1 - 1 e X+Y -1 . (1.19)
In Theorems 1.4 and 1.6 below we suppose that all r.v.'s have positive and twice differentiable densities.

First we consider the case f = f 1 . We determine the distributions of X and Y such that f 1 is a LWMY function associated to (X, Y ). Theorem 1.4

1. Consider two positive and independent random variables X and Y . The random variables U and V defined by (1.19) are independent if and only if the densities of Y and X are respectively

p Y (y) = Γ(a + b) Γ(a)Γ(b) (1 -e -y ) b-1 e -ay 1 {y>0} (1.20) p X (x) = α(a + b, c, -a)e -(a+b)x (1 -e -x ) -b-1 exp -c e -x 1 -e -x 1 {x>0} . (1.21)
where a, b and c are constants such that a, b, c > 0 and α(a+ b, c, -a) is the constant of Equation (1.18). Thus the law of Y is the image of the Beta(a, b) distribution by the transformation z ∈ (0, 1) →log z ∈ (0, ∞), while the law of the variable

f 1 (X) is K (2) (a + b, -b, c) (cf Equation (1.18)). 2. If 1. holds then U ∼ K (2) (a, b, c) and V ∼ γ(b, c).
The proof of Theorem 1.4 will be given in Section 4. Since

Remark 1.5 Since g 1 = f -1 1 , item 2. of Remark 1.
T g 1 = T -1 f 1 , then (X, Y ) = T g 1 (U, V ) = log 1 + U + V U + V , log 1 + U U -log 1 + U + V U + V . (1.22) As shows (1.22) it is useful to introduce (U ′ , V ′ ) = 1 + 1 U +V 1 + 1 U , U + V . (1.23) Obviously the correspondence (U, V ) → (U ′ , V ′ ) is one-to one: (U, V ) = U ′ V ′ V ′ + 1 -U ′ V ′ , V ′ (V ′ + 1)(1 -U ′ ) V ′ + 1 -U ′ V ′ . (1.24)
Furthermore, (X, Y ) can be easily expressed in terms of (U ′ , V ′ ):

X = log(1 + 1/V ′ ) and Y = -log U ′ . (1.25)
Since it is easy to determine the density function of φ(ξ) where φ is differentiable and bijective, then Theorem 1.4 and its analogue related to f = g 1 (cf Remark 1.5) are equivalent to Theorem 1.6 below.

Theorem 1.6 a) Let U ′ and V ′ be two positive and independent random variables. The r.v.'s U and V defined by (1.24) are independent if only if there exist some constants a, b, c such that

U ′ ∼ Beta(a, b) and V ′ ∼ K (2) (a + b, -b, c). (1.26) b)
Let U and V be two positive and independent random variables. The r.v.'s U ′ and V ′ defined by (1.23) are independent if only if there exist some constants a, b, c such that

U ∼ K (2) (a, b, c) and V ∼ γ(b, c).
(1.27)

Under (1.27), U ′ ∼ Beta(a, b) and V ′ ∼ K (2) (a + b, -b, c). c) If one of these equivalent conditions holds, then U ∼ K (2) (a, b, p) and V ∼ γ(b, c).
Let us formulate an easy consequence of Theorem 1.6.

Theorem 1.7 For any a, b, c > 0, the transformation (u, v) → (

1+ 1 u+v 1+ 1 u , u + v) maps the probability measure K (2) (a, b, c) ⊗ γ(b, c) to the probability measure Beta(a, b) ⊗ K (2) (a + b, -b, c).
In particular:

K (2) (a, b, c) * γ(b, c) = K (2) (a + b, -b, c).
(1.28) Remark 1.8

1. Note that (1.28) may be regarded as an analogue of (0.2).

2. Theorem 1.6 and the results of [START_REF] Letac | An independence property for the product of GIG and gamma laws[END_REF] suggest to ask whether it is possible to have an "algebraic" proof of the following part of either item a) or item b) of Theorem 1.6: if U and V , U ′ and V ′ are independent, then (1.26) and (1.27) hold.

The case

f = f * δ
Recall that f * δ has been defined by (1.17). Due to the form of f * δ , a change of variables allows to simplify the search of independent r.v.'s X and Y such that the two components of T f * δ (X, Y ) are independent. For any decreasing and bijective function f : (0, ∞) → (0, ∞) we define

f(x) = exp{-f (-log x)}, x ∈ (0, 1), (1.29) T m f (x, y) = f (xy), f (x) f (xy) , x, y ∈ (0, 1). (1.30)
Observe that f is one-to-one from (0, 1) onto (0, 1), T m f is one-to-one from (0, 1) 2 onto (0, 1) 2 and

T m f -1 = T m f -1 .
(1.31) Definition 1.9 Let X and Y be two independent and (0, 1)-valued random variables. We say that a decreasing and bijective function f : (0, 1) → (0, 1) is a multiplicative LWMY function with respect to (X, Y ) if the r.v.'s U m := f (XY ) and V m := f (X) f (XY ) are independent.

Remark 1.10 For any random vector (X, Y ) in (0, ∞) 2 we consider

X ′ = e -X , Y ′ = e -Y .
Then f is a LWMY function with respect to (X, Y ) if and only if f is a multiplicative LWMY function with respect to (X ′ , Y ′ ).

The change of variable x ′ = e -x is very convenient since the function

φ δ (x) := f * δ (x) = 1 -x 1 + (δ -1)x , x ∈ (0, 1) (1.32)
is homographic. Note that f * δ : (0, 1) → (0, 1) is bijective, decreasing and self-reciprocal. First, let us determine the distribution of the couple (X ′ , Y ′ ) of r.v.'s such that φ δ is a multiplicative LWMY function with respect to (X ′ , Y ′ ).

For a, b, α > 0 and c ∈ R consider the probability measure

β α (a, b; c)(dx) = k α (a, b; c)x a-1 (1 -x) b-1 (αx + 1 -x) c 1 (0,1) (x)dx. (1.33) Note that if c = 0, then β α (a, b; c) = Beta(a, b).
Theorem 1.11 Let X ′ and Y ′ be two independent random variables valued in (0, 1). Consider

(U m , V m ) = T m φ δ (X ′ , Y ′ ) = 1 -X ′ Y ′ 1 + (δ -1)X ′ Y ′ , 1 -X ′ 1 + (δ -1)X ′ 1 + (δ -1)X ′ Y ′ 1 -X ′ Y ′
for fixed δ > 0.

Then, U m and V m are independent if and only if there exist a, b, λ > 0 such that

X ′ ∼ β δ (a + b, λ; -λ -b), Y ′ ∼ Beta(a, b). (1.34)
If this condition holds, then

U m ∼ β δ (λ + b, a; -a -b), V m ∼ Beta(λ, b). (1.35)
In the case δ = 1, Theorem 1.11 takes a very simple form.

Proposition 1.12 Let X ′ and Y ′ be two independent random variables valued in (0, 1).

Then U m = 1 -X ′ Y ′ , V m = 1 -X ′ 1 -X ′ Y ′ are independent if and only if there exist a, b, λ > 0 such that X ′ ∼ Beta(a + b, λ) and Y ′ ∼ Beta(a, b).

If one of these conditions holds, then

U m ∼ Beta(λ + b, a), and V m ∼ Beta(λ, b). Remark 1.13 When X ′ ∼ Beta(a + b, λ) and Y ′ ∼ Beta(a, b) it can be proved that U m and V m are independent using the well-known property: if Z and Z ′ are independent, Z ∼ γ(a, 1) and Z ′ ∼ γ(b, 1) then R := Z Z+Z ′ and Z + Z ′ are independent and R ∼ Beta(a, b) and Z + Z ′ ∼ γ(a + b, 1) (see for instance [START_REF] Yor | Une extension markovienne de l'algèbre des lois béta-gamma[END_REF]).

According to Remark 1.10, f * δ is a LWMY function with respect to (X, Y ) if and only if φ δ is a multiplicative LWMY function with respect to (X ′ , Y ′ ) = (e -X , e -Y ). Therefore, a classical change of variables allows to deduce that Theorem 1.11 is equivalent to Theorem 1.14 below: Theorem 1. [START_REF] Matsumoto | Interpretation via Brownian motion of some independence properties between GIG and gamma variables[END_REF] 1. Consider two positive and independent random variables X and Y . The random variables

U = f * δ (X + Y ), V = f * δ (X) -f * δ (X + Y )
p Y (y) = Γ(a + b) Γ(a)Γ(b) (1 -e -y ) b-1 e -ay 1 {y>0} (1.36) p X (x) = k δ (a + b, λ, -λ -b)e -(a+b)x (δe -x + 1 -e -x ) -λ-b ×(1 -e -x ) λ-1 1 x>0 (1.37) where a, b > 0 , λ ∈ R and k δ (a + b, λ, -λ -b) is the normalizing factor (cf (1.33)). Thus e -Y is Beta(a, b) distributed and e -X is β δ (a + b, λ, -λ -b) distributed.
2. If 1. holds then the densities of U and V are respectively

p U (u) = k δ (λ + b, a; -a -b)e -u(λ+b) (1 -e -u ) a-1 ×(1 + (δ -1)e -u ) -a-b 1 u>0 , (1.38) p V (v) = e -λv (1 -e -v ) b-1 1 v>0 .
(1.39)

The proof of Theorem 1.14 is similar to that of Theorem 1.4 and has been postponed to Appendix 5.3.

2 A necessary and sufficient condition for LWMY functions and related densities

φ ′′ X (x) -φ ′ X (x) f ′′ (x) f ′ (x) + φ ′′ Y (y)f ′ (x) 1 f ′ (x) - 1 f ′ (x + y) +φ ′ Y (y) f ′′ (x) f ′ (x) + 2(f ′′ (x)) 2 -f ′′′ (x)f ′ (x) f ′ (x) 2 = 0, x, y > 0. ( 2

.1)

Proof :

1) Consider the transformation T f : (0, ∞) 2 → (0, ∞) 2 ; (x, y) → (f (x + y), f (x)f (x + y)). Let g = f -1 and (U, V ) = T f (X, Y ). By formula (1.13), (X, Y ) = T g (U, V ). X and Y being independent, the density of (U, V ) is

p (U,V ) (u, v) = p X (g(u + v)) p Y (g(u) -g(u + v)) |J(u, v)|1 u,v>0 (2.2)
where J is the Jacobian of the transformation T f . Since x = g(u+v) and y = g(u)-g(u+v) we have

J(u, v) = det g ′ (u + v) g ′ (u + v) g ′ (u) -g ′ (u + v) -g ′ (u + v) = -g ′ (u + v)g ′ (u).
g being monotone we have

|J(u, v)| = g ′ (u + v)g ′ (u), so that p (U,V ) (u, v) = p X (g(u + v)) p Y (g(u) -g(u + v)) g ′ (u + v)g ′ (u). (2.3)
2) Let H : (0, ∞) 2 → R be a function having second partial derivatives. It is easy to prove that the function H decomposes as

H(u, v) = h 1 (u) + h 2 (v), u, v > 0, if and only if ∂ 2 H ∂u∂v = 0.
3) According to the former item, the variables U and V are independent if and only if the function

H = log p (U,V ) satisfies ∂ 2 H ∂u∂v = 0. By Equation (2.3) we have ∂H ∂v = φ ′ X (g(u + v)) g ′ (u + v) -φ ′ Y (g(u) -g(u + v)) g ′ (u + v) + g ′′ g ′ (u + v), ∂ 2 H ∂u∂v = φ ′′ X (g(u + v)) (g ′ (u + v)) 2 + φ ′ X (g(u + v)) g ′′ (u + v) -φ ′′ Y (g(u) -g(u + v)) g ′ (u + v) [g ′ (u) -g ′ (u + v)] -φ ′ Y (g(u) -g(u + v)) g ′′ (u + v) + g ′′′ g ′ -(g ′′ ) 2 (g ′ ) 2 (u + v).
Writing x = g(u + v) and y = g(u)g(u + v), we have u = f (x + y), v = f (x)f (x + y) and

∂ 2 H ∂u∂v = φ ′′ X (x)[g ′ (f (x))] 2 + φ ′ X (x)g ′′ (f (x)) -φ ′′ Y (y)g ′ (f (x)) [g ′ (f (x + y)) -g ′ (f (x))] -φ ′ Y (y)g ′′ (f (x)) + g ′′′ g ′ -(g ′′ ) 2 (g ′ ) 2 (f (x)) = 0. (2.4)
Differentiating twice the relation g(f (x)) = x, we obtain

g ′ (f (x)) = 1 f ′ (x) , (x > 0) (2.5) g ′′ (f (x)) = - f ′′ (x) f ′ (x) 3 (x > 0). ( 2

.6)

A differentiation of the latter equality yields 5 .

g ′′′ (f (x)) = - f ′′′ (x)f ′ (x) -3f ′′ (x) 2 f ′ (x)
Using this equality and (2.5) and (2.6) we have

g ′′′ g ′ -(g ′′ ) 2 (g ′ ) 2 (f (x)) = - f ′′′ (x)f ′ (x) -3f ′′ (x) 2 f ′ (x) 5 1 f ′ (x) - f ′′ (x) 2 f ′ (x) 6 f ′ (x) 2 = 2f ′′ (x) 2 -f ′′′ (x)f ′ (x) f ′ (x) 4 . ( 2
∂ 2 H ∂u∂v = φ ′′ X (x) f ′ (x) 2 - φ ′ X (x)f ′′ (x) f ′ (x) 3 - φ ′′ Y (y) f ′ (x) 1 f ′ (x + y) - 1 f ′ (x) + φ ′ Y (y)f ′′ (x) f ′ (x) 3 + 2f ′′ (x) 2 -f ′′′ (x)f ′ (x) f ′ (x) 4
and (2.1) follows by multiplying the above identity by f ′ (x) 2 .

The set of all possible "smooth" LWMY functions

We restrict ourselves to smooth LWMY functions f , i.e. satisfying f : (0, ∞) → (0, ∞) is bijective and decreasing, (3.1)

f is three times differentiable, (3.2)
F (x) = n≥1 a n x n , ∀x > 0. ( 3.3)
where F := 1/f ′ . According to (3.1), f ′ (0 + ) = -∞. This implies F (0 + ) = 0 and explains why the series in (3.3) starts with n = 1. The goal of this section is to prove half of Theorem 1.2: if f is a smooth LWMY function, then f belongs to one of the four classes F 1 , . . . , F 4 introduced in Remark 1.3. First, we characterize in Theorem 3.1 all possible functions F . Second, we determine the associated functions f (see Theorem 3.8).

Theorem 3.1 Suppose that f is a smooth LWMY function and that the assumptions of Theorem 2.1 are satisfied.

1. If F ′ (0 + ) = 0, then a 2 < 0 and

F (x) = a 2 2 6a 4 cosh x 12a 4 a 2 -1 if a 4 < 0 a 2 x 2 otherwise. (3.4) 2. If F ′ (0 + ) = 0, then F (x) = a 1 a 2 3a 3 cosh x 6a 3 a 1 -1 + a 1 a 1 6a 3 sinh x 6a 3 a 1 if a 1 a 3 > 0 a 1 x + a 2 x 2 otherwise. (3.5) Remark 3.2 Unsurprisingly, the case F (x) = a 2 x 2 corresponds to f (x) = -1 a 2 1
x , i.e. the case considered by Matsumoto-Yor and Letac-Wesolowski.

Throughout this subsection we suppose that f satisfies (3.1), (3.2), (3.3) and that the assumptions of Theorem 2.1 are fulfilled. For simplicity of statement of results below we do not repeat these conditions.

Recall that φ Y is the logarithm of the density of Y . Let us introduce

h := φ ′ Y (3.6) Lemma 3.3 1.
There exists a function λ : (0, ∞) → R such that

F (x + y) = λ(x) -h(y)F ′ (x) h ′ (y) + F (x). (3.7)
2. F satisfies

F (y) = λ(0 + ) -h(y)F ′ (0 + ) h ′ (y) . (3.8)
Remark 3.4 Suppose that we have been able to determine F . Then, h = φ ′ Y solves the linear ordinary differential equation (3.8) and can therefore be determined. The remaining function φ X is obtained by solving Equation (2.1).

Proof of Lemma 3.3 : Using ( 3.6) and F = 1/f ′ in Equation (2.1), we obtain

c(x) = h(y) F ′ (x) F (x) + h ′ (y) 1 F (x) (F (x + y) -F (x))
where c(x) depends only on x. Multiplying both sides by F (x) and taking the y-derivative

leads to 0 = F ′ (x)h ′ (y) + [F (x + y) -F (x)]h ′′ (y) + h ′ (y)F ′ (x + y).
Fix x > 0. Then θ(y) := F (x + y) is a solution of the differential equation in y:

0 = F ′ (x)h ′ (y) + (θ(y) -F (x))h ′′ (y) + h ′ (y)θ ′ (y). (3.9)
A solution of the related homogeneous equation in y is ρ h ′ (y) where ρ is a constant. It is easy to prove that y → -F ′ (x)h(y) + F (x)h ′ (y) solves (3.9). Thus, the general solution of (3.9) is

θ(y) = 1 h ′ (y) [λ(x) -F ′ (x)h(y) + F (x)h ′ (y)] .
Since θ(y) = F (x + y), (3.7) follows.

According to (3.3), F (0 + ) and F ′ (0 + ) exist. Therefore, taking the limit x → 0 + in (3.7) implies both the existence of λ(0 + ) and relation (3.8).

The following lemma shows that the function F (and thus f ) solves a self-contained equation in which h, and thereby the densities of X and Y , are not involved. Lemma 3.5 F solves the delay equations :

F (x + y) = F (y)[λ(x) -h(y)F ′ (x)] λ(0 + ) -h(y)F ′ (0 + ) + F (x) (x, y > 0) (3.10) F ′ (x + y) = F ′ (y) + F ′ (0 + ) F (y) [F (x + y) -F (x)] -F ′ (x) (x, y > 0). (3.11)
Proof: By (3.8) we have

h ′ (y) = λ(0 + ) -h(y)F ′ (0 + ) F (y) .
Equation (3.10) then follows by rewriting Equation (3.7) and replacing h ′ (y) with the expression above. We differentiate (3.10) in y and use the fact that λ(0 + )h(y)F ′ (0 + ) = h ′ (y)F (y) to obtain :

F ′ (x + y) = [F ′ (y) + F ′ (0 + )] λ(x) -h(y)F ′ (x) F (y)h ′ (y) -F ′ (x). By (3.7) we have λ(x) -h(y)F ′ (x) F (y)h ′ (y) = F (x + y) -F (x) F (y)
and this gives (3.11).

Remark 3.6 We can see (3.11) as a scalar neutral delay differential equation. Indeed, set t = x + y and consider y > 0 as a fixed parameter. Then (3.11) becomes:

F ′ (t) = a(F (t) -F (t -y)) -F ′ (t -y), t ≥ y, (3.12)
where a := F ′ (y)+F ′ (0 + ) F (y)

. Replacing F (t) in (3.12) with e at G(t) leads to:

G ′ (t) + e -ay G ′ (t -y) + 2ae -ay G(t -y) = 0, t ≥ y.
(3.13) Equation (3.13) is called a neutral delay differential equation (cf for instance, Section 6.1, in [START_REF] Györi | Oscillation theory of delay differential equations[END_REF]. These equations have been intensively studied but the authors only focused on the asymptotic behaviour of the solution as t → ∞. Unfortunately, these results give no help to solve explicitly either (3.11) or (3.13).

Lemma 3.7 For all integers k ≥ 0 and l ≥ 1, we have

l-1 m=0 (l -2m + 1)C k l-m+1+k a l-m+1+k a m = (l -2)(k + 1)a k+1 a l + a 1 a l+k C k l+k , (3.14)
C k k+3 a k+3 a 1 = (k + 1)a k+1 a 3 , (3.15) 2C k k+4 a k+4 a 1 + C k k+3 a k+3 a 2 -C k k+2 a k+2 a 3 -2(k + 1)a k+1 a 4 = 0, (3.16)
where

C p n = n! (n-p)!p! .
Proof: Obviously Equation (3.11) is equivalent to:

F ′ (x + y)F (y) = F ′ (y)F (x + y) -F ′ (y)F (x) -F (y)F ′ (x) + F ′ (0 + )F (x + y) -F ′ (0 + )F (x).
(3.17) Using the asymptotic expansion (3.3) of F we can develop each term in (3.17) as a series with respect to x and y. Then, identifying the series on the right-hand side and the left-hand side we get (3.14)- (3.16). Since this approach only needs usual skill in analysis the details are provided in Appendix 6.2.

Proof of part 1. of Theorem 3.1

We suppose that F ′ (0 + ) = 0. Since a 1 = F ′ (0 + ) = 0, we necessarily have a 2 = 0. Indeed, if a 2 = 0 then, by (3.16) with k = 1, we would have -3a 2 3 -4a 2 a 4 = 0, i.e. a 3 = 0, and using again (3.16) with k = 3 would imply a 4 = 0 and finally a k = 0 for every k ≥ 0, which is a contradiction because, by definition, F = 1/f ′ does not vanish.

So, we have a 1 = 0 and a 2 = 0. Equation (3.15) with k = 1 reads 4a 4 a 1 = 2a 2 a 3 , which implies a 3 = 0. Applying (3.15) to k = 2n provides, by induction on n, a 2n+1 = 0 for every n ≥ 0.

Therefore, Equation (3.16) reduces to

(k + 3)(k + 2)(k + 1)a k+3 a 2 = 12(k + 1)a k+1 a 4 , (k ≥ 0) i.e. a k+3 = 12a 4 a 2 1 (k + 3)(k + 2) a k+1 .
This leads to

a 2k = 12a 4 a 2 k-1 2 (2k)! a 2 , k ≥ 1. (3.18)
Then, F (x) = a 2 x 2 if a 4 = 0 and if a 4 = 0 we have

F (x) = k≥1 12a 4 a 2 k-1 2 (2k)! a 2 x 2k .
If a 4 a 2 < 0, then

F (x) = a 2 2 6a 4 cos x -12a 4 a 2 -1 .
This implies F (2π -12a 4 a 2 ) = 0 which is impossible since

F (x) = 1/f ′ (x) < 0. Conse- quently, F (x) = a 2 2 6a 4 cosh x 12a 4 a 2 -1 .
Proof of part 2. of Theorem 3.1 Using Equation (3.15) we have, for all k ≥ 0, (k + 3)(k + 2)a k+3 a 1 = 6a k+1 a 3 . Since a 1 = 0, Equation (3.15) is equivalent to

a k+3 = 6a 3 a 1 1 (k + 3)(k + 2) a k+1 .
As a result,

a 2k+1 = 6a 3 a 1 k 1 (2k + 1)! a 1 , (k ≥ 0), a 2k = 6a 3 a 1 k-1 2 (2k)! a 2 , (k ≥ 1).
As a consequence,

F (x) = k≥1 6a 3 a 1 k-1 2 (2k)! a 2 x 2k + k≥0 6a 3 a 1 k 1 (2k + 1)! a 1 x 2k+1 .
If a 3 = 0 then it follows from Equation (3.15) that a k = 0 for every k ≥ 3. If a 3 = 0 and a 1 a 3 > 0 we have

k≥1 6a 3 a 1 k-1 2 (2k)! a 2 x 2k = a 2 a 1 3a 3 cosh x 6a 3 a 1 -1 and k≥0 6a 3 a 1 k 1 (2k + 1)! a 1 x 2k+1 = a 1 a 1 6a 3 sinh x 6a 3 a 1 .
If a 1 a 3 < 0 then

F (x) = a 1 a 2 3a 3 cos x -6a 3 a 1 -1 + -a 1 6a 3 sin x -6a 3 a 1 .
In particular, F (2π -a 1 6a 3 ) = 0. This is impossible since

F (x) = 1/f ′ (x) < 0.
Now, in each case of Theorem 3.1 we compute the corresponding f associated with F via the relation F = 1/f ′ . Recall that we restrict ourselves to functions f satisfying (3.1)-(3.3) and work under the assumptions of Theorem 2.1.

Theorem 3.8 1. If F (x) = a 2 x 2 then f (x) = 1 a 2 x . 2. If F (x) = α(cosh βx -1), α, β > 0, then f (x) = 2 αβ f 1 (βx). 3. If F (x) = a 1 x + a 2 x 2 then f (x) = -1 a 1 g 1 ( a 2 a 1 x). 4. If F (x) = a 1 a 2 3a 3 cosh x 6a 3 a 1 -1 + a 1 a 1 6a 3 sinh x 6a 3 a 1 then f (x) = - 1 βγ log e βx + δ -1 e βx -1 ,
where α = a 1 a 2 3a 3 , β = 6a 3 a 1 and γ = a 1 a 1 6a 3 .
Proof: 1) The first case is obvious.

2) We have

f ′ (x) = 1 F (x) = - 1 α(cosh βx -1) , x > 0, β > 0, α > 0.
Since f (0 + ) = +∞ and f (+∞) = 0 + , integrating the previous identity gives

f (x) = 2 αβ 1 e βx -1 = 2 αβ f 1 (βx),
where

f 1 (x) = 1 e x -1 . 3) Recall (cf (1.16)) that g 1 = f -1 1 . Since g ′ 1 (x) = -1 x(x+1)
, it can be easily proved that the associated LWMY function is f (x) = -1 a 1 g 1 ( a 2 a 1 x). 4) Due to our choice of α, β and γ we have

f ′ (x) = 1 F (x) = 1 α(cosh βx -1) + γ sinh βx (3.19) = e βx γ(e βx -1) - (α + γ)e βx γ[(α + γ)e βx -α + γ] .
As a consequence, by integration,

f (x) = 1 βγ log(e βx -1) - 1 βγ log |(α + γ)e βx -α + γ| + C
where C is a constant. Using (3.19) we have :

f ′ (x) ∼ 1 αβ 1
x as x → 0 + , f ′ (x) ∼ 2 α+γ e -βx as x → +∞. This implies that γ < 0 and α + γ < 0. Setting δ = 2γ α + γ we have δ > 0, γα α + γ = δ -1 and

f (x) = - 1 βγ log e βx + δ -1 e βx -1 + C ′ . Since f (∞) = 0 we have C ′ = 0. As a result, f (x) = -1 βγ f * δ (βx).
4 Proof of Theorem 1.4

Recall that f 1 , g 1 and f * δ have been defined by equations (1.15), (1.16) and (1.17) respectively. In Theorem 3.8 we have proved that if f is a LWMY function, then f ∈ F 1 ∪ F 2 ∪ F 3 ∪ F 4 , where these four classes have been introduced in Remark 1.3. The class F 1 corresponds to the Matsumoto-Yor property. Indeed, we prove in Subsection 6.1 that x → 1/x is a LWMY function associated to GIG and gamma distributions. The result has been already obtained by Letac and Wesolowski (2001) under weaker assumptions than ours.

Recall that φ Y = log p Y , h = φ ′ Y and F ′ (0 + ) = 0. It is easy to deduce from (3.8) that there exist constants λ and c 1 such that This proves (1.20).

h(y) = λf (y) + c 1 = λ e y -1 + c 1 = λe y e y -1 + c 1 -λ.
To find the density of X we come back to Equation (2.1) and compute each of its terms.

We have f ′ (x) = -e x (e x -1) 2 , f ′′ (x) = e 2x +e x (e x -1) 3 , f ′′′ (x) = -e 3x +4e 2x +e x (e x -1) 4 , so that

f ′ (x) f ′ (x+y) = e -y (e x+y -1) 2 (e x -1) 2 and f ′′ (x) f ′ (x) = -e x +1 e x -1 . Calculations yield 2(f ′′ (x)) 2 -f ′′′ (x)f ′ (x) f ′ (x) 2 = e 2x + 1 (e x -1) 2 . (4.2)
Moreover,

-φ ′ Y (y) f ′′ (x) f ′ (x) + φ ′′ Y (y) f ′ (x) f ′ (x + y) -1 = λ e y -1 + c 1 e x + 1 e x -1 - λe y (e y -1) 2 e -y (e x+y -1) 2 (e x -1) 2 -1 = (c 1 -λ)e 2x -c 1 (e x -1) 2 . ( 4.3)
Equation (2.1) can then be writen, using (4.2) and (4.3):

φ ′′ X (x) + e x + 1 e x -1 φ ′ X (x) = - e 2x + 1 (e x -1) 2 + (c 1 -λ)e 2x -c 1 (e x -1) 2 = (c 1 -λ -1)e 2x -c 1 -1 (e x -1) 2 .
Then Proof : We keep the notation introduced in Section 1. The densities of X and Y are respectively denoted by p X and p Y . Let φ X = log p X and φ Y = log p Y . Theorem 2.1 with f : x → 1/x asserts that the variables U and V are independent if and only if, for all x, y > 0,

h 0 := φ ′ X solves h ′ 0 (x) + e x + 1 e x -1 h 0 (x) = (c 1 -λ -1)e 2x -c 1 -1 (e x -1) 2 . ( 4
φ ′′ X (x) + 2 x φ ′ X (x) + φ ′′ Y (y) 1 x 2 (x 2 -(x + y) 2 ) - 2 x φ ′ Y (y) + 2 x 2 = 0. (5.5)
Let us first compute the density of Y . Equation (5.5) implies that, for fixed x > 0, the function φ ′ Y = h is a solution of the ordinary linear differential equation

2 x h(y) + h ′ (y) y 2 + 2xy x 2 = c(x), y > 0, (5.6)
where c(x) does not depend on y. The solutions of h ′ (y) + 2x y 2 +2xy h(y) = 0 are of the type h(y) = k(x) y + 2x y .

Note that y → -c(x) x 2 y is a particular solution of (5.6), consequently h solves (5.6) iff:

h(y) = -c(x) x 2 y + k(x) y + 2x y , y > 0. Recall that φ ′ Y = h, then we have φ ′ Y (y) = 1 y (-c(x)x 2 + 2xk(x) + yk(x)). But φ ′
Y is a function of y only. As a consequence, the expressions -c(x)x 2 + 2xk(x) and k(x) actually do not depend on x. We denote them respectively by δ and λ, so that where θ = e d is a positive constant. Since p Y is, as a density, integrable, we have necessarily λ < 0 and δ > -1 and defining

µ := δ + 1, b = √ -2λ, (5.8)
formula (5.7) shows that the law of Y is γ(µ, b 2 /2).

We now compute the density of X. We rewrite Equation (5.5) by replacing φ ′ Y (y) with 1 y (δ + λy) and φ ′′ Y (y) with -δ y 2 . This gives

φ ′′ X (x) + 2 x φ ′ X (x) + δ + 2 x 2 - 2λ x = 0.
So φ ′ X is a solution of the ordinary linear differential equation:

κ ′ (x) + 2 x κ(x) + δ + 2 x 2 - 2λ x = 0.
(5.9)

The solution of the corresponding homogeneous differential equation κ

′ (x) + 2 x κ(x) = 0 is κ(x) = η x 2 , for some constant η (x > 0). Note that x → (-2-δ)x+λx 2
x 2 solves (5.9), then κ solves (5.9) iff for some η, κ(

x) = (-2-δ)x+λx 2 +η x 2 . Since κ = φ ′ X we have φ ′ X (x) = -2 -δ x + λ + η x 2 . As a consequence, there exists a constant τ ∈ R such that φ X (x) = (-2 -δ) ln x + λx - η x + τ.
Since φ X = log p X , we have

p X (x) = x -2-δ (exp τ ) exp(λx - η x )
which proves that X follows the law GIG(-µ, a, b) with µ := δ + 1, b = √ -2λ (cf (5.8)) and a = √ 2η (observe that the integrability of p Y implies η > 0) .

Proof of Lemma 3.7

We have

F ′ (x + y)F (y) = n,m≥0 a n a m n(x + y) n-1 y m = n,m≥0 a n a m ny m n-1 k=0 C k n-1 x k y n-1-k = k≥0 x k m≥0,n≥1+k na n a m C k n-1 y n+m-1-k . Setting l = m + n -1 -k for fixed m gives F ′ (x + y)F (y) = k≥0,l≥0 x k y l l m=0 (l -m + 1 + k)C k l-m+k a l-m+1+k a m .
(5.10)

F ′ (y)F (x + y) = n,m≥0 ma n a m (x + y) n y m-1 = n,m≥0 ma n a m y m n k=0 C k n x k y n-k y m-1 = k≥0 x k m≥0,n≥k ma n a m C k n y n+m-1-k .
We set again l = m + n -1k for fixed k :

F ′ (y)F (x + y) = k≥0,l≥0
x k y l l+1 m=0 mC k l-m+k+1 a l-m+1+k a m (5.11)

F ′ (y)F (x) = k≥0,l≥0
a k a l+1 (l + 1)x k y l (5.12)

F ′ (x)F (y) = k≥0,l≥0
a k+1 a l (k + 1)x k y l (5.13)

F ′ (0 + )F (x + y) = a 1 n≥0 a n (x + y) n = a 1 n≥0 a n n k=0 C k n x k y n-k = a 1 k≥0 x k n≥k a n C k n y n-k = a 1 k,l≥0
a l+k C k l+k x k y l .

(5.14)

F ′ (0 + )F (x) = a 1 k≥0 a k x k . (5.15)
Identifying the coefficient of x k y l in (3.17 Note that if l = 0, both sides of (5.16) vanish, therefore we may suppose in the sequel that l ≥ 1. For m = l + 1 we have mC k l-m+k+1 a l-m+1+k a m = (l + 1)a k a l+1 . Thus, Equation (5.16) reads (5.17)

But one finds by a calculation using the definition that (lm + 1 + k)C k l-m+k -mC k l-m+1+k = (l -2m + 1)C k l-m+1+k , so that Equation (5.17) is equivalent to l m=0 (l -2m + 1)C k l-m+1+k a l-m+1+k a m = -(k + 1)a k+1 a l + a 1 a l+k C k l+k .

(5.18)

For m = l we have (l -2m + 1)C k l-m+1+k a l-m+1+k a m = (1l)(k + 1)a k+1 a l . Consequently Equation (5.18) may be written as follows: -φ ′ Y (y)

f ′′ (x) f ′ (x) -φ ′′ Y (y)f ′ (x) 1 f ′ (x) - 1 f ′ (x + y)
= -(a + b -1)e x + a(δ -1)e -x (e x + δ -1)(1e -x ) .

Equation (2.1) is then equivalent to : . This proves (1.37).

To compute the density of (U, V ) we note that f = f * δ is self-reciprocal. Then, by (2.3), the density of (U, V ) is

p (U,V ) (u, v) = p X (f (u + v)) p Y (f (u) -f (u + v)) f ′ (u + v)f ′ (u)
= p X log e u+v + δ -1 e u+v -1 p Y log (e u + δ -1)(e u+v -1) (e u -1)(e u+v + δ -1) × δ 2 (e u+v + δ -1)(1e -u-v )(e u + δ -1)(1e -u ) .

Using the expressions of p X and p Y , one gets (1.38) and (1.39).

  b) Barndorff-Nielsen and Halgreen (1977) proved: GIG(-µ, a, b) * γ(µ, b 2 2 ) = GIG(µ, a, b), µ, a, b > 0 (0.2)

  3 and Theorem 1.4 imply that the r.v.'s associated with the LWMY function g 1 are the r.v.'s U and V distributed as in item 2. of Theorem 1.4. b) As suggest identities (1.20) and (1.21) it is actually possible to simplify Theorem 1.4.

Theorem 2 . 1

 21 Let X and Y be two independent and positive random variables whose densities p X and p Y are positive and twice differentiable. Define φ X = log p X and φ Y = log p Y . Consider a decreasing function f : (0, ∞) → (0, ∞), three times differentiable. Then f is a LWMY function with respect to (X, Y ) if and only if

0 M( 1 - 1 0( 1 -

 0111 This implies the existence of a constant d such that φ Y (y) = λ log(e y -1) + (c 1λ)y + d. Setting M = e d , we have by integration, for all y > 0, p Y (y) = M(e y -1) λ e (c 1 -λ)y = M(1e -y ) λ e c 1 y . (4.1) To give more information on the normalizing constant M, one observes, for a = -c 1 and b = λ + 1, that ∞ e -y ) b-1 e -ay dy = M u) b-1 u a-1 du which implies that a > 0, b > 0 and M = Γ(a + b) Γ(a)Γ(b) .

1 Proposition 5 . 1

 151 Another proof of the result by Letac and Wesolowski under stronger assumptions Consider two independent and positive random variables X and Y having positive and twice differentiable densities. If the random variables U = (X + Y ) -1 and V = X -1 -(X + Y ) -1 are independent then there exist µ > 0, a > 0 and b > 0 such that the law of X is GIG(-µ, a, b) and Y follows the gamma distribution γ(µ, b 22 ).

  λy), y > 0, which implies the existence of a constant d ∈ R such that φ Y (y) = δ ln y + λy + d, y > 0. But φ Y = log p Y , so the density of Y is p Y (y) = θy δ e λy , y > 0.(5.7)

  ) and using (5.10) to (5.15) we have, for k ≥ 0 and l ≥ 0:l m=0 (lm + 1 + k)C k l-m+k a l-m+1+k a m = -(l + 1)a k a l+1 -(k + 1)a k+1 a l a l-m+1+k a m +a 1 a l+k C k l+ka 1 a k 1 l=0 . (5.16)

  m + 1 + k)C k l-m+k a l-m+1+k a m = -(k + 1)a k+1 a l + l m=0 mC k l-m+k+1 a l-m+1+k a m +a 1 a l+k C k l+k .

  2m + 1)C k l-m+1+k a l-m+1+k a m -(l -1)(k + 1)a k+1 a l = -(k + 1)a k+1 a l + a 1 a l+k C l l+k which implies (3.14). (3.15) and (3.16) follow by applying (3.14) to l = 3 and l = 4 respectively.5.3 Proof of Theorem 1.14Throughout this proof we write for simplicity f instead of f * δ . Starting with f (x) = log(e x + δ -1)log(e x -1) we havef ′ (x) = -δe x (e x + δ -1)(e x -1) = -δ e x + δ -2 + (1δ)e -x (5.19) which implies F (x) = 1/f ′ (x) = -1 δ (e x + δ -2 + (1δ)e -x ) and F ′ (0) = -1. Therefore, again with h = φ ′ Y , Equation (3.8) is equivalent to h ′ (y) δ (e y + δ -2 + (1δ)e -y )h(y) = λ,whose solution is found to be, again by the method of variation of constants :h(y) = c + δ(c + λ)e -y1e -y .(5.20)Thus, there exists a constant d ′ such that φ Y (y) = cy + δ(c + λ) log(1e -y ) + d ′ and this proves that p Y (y) = d(1e -y ) (c+λ)δ e cy 1 y>0 . We seta = -c, b = δ(c + λ) + 1 (5.21)to obtain (1.36). ¿From(5.19) one hasf ′′ (x) = δ(e x + (δ -1)e -x ) (e x + δ -2 + (1δ)e -x ) 2 , f ′′′ (x) = δ -e 2x + (δ -2)e x -6δ + 6 + (δ -1)(2δ)e -x -(δ -1) 2 e -2x (e x + δ -2 + (1δ)e -x ) 3 , 2(f ′′ (x)) 2f ′′′ (x)f ′ (x) f ′ (x) 2 = ex -(δ -1)e -x (e x + δ -1)(1e -x ) ,

  φ ′′ X (x) + e x + (δ -1)e -x (e x + δ -1)(1e -x ) φ ′ X (x) = (δ -1)(1a)e -x -(a + b)e x (e x + δ -1)(1e -x) .Solving this differential equation by similar calculations givesφ ′ X (x) = K + (a -1)(δ -1)e -x -(a + b)e x (e x + δ -1)(1e -x )where K is a constant. This can be writtenφ ′ X (x) = A + Be -x 1 + (δ -1)e -x + Ce -x 1e -x where A = -ab, B = (δ -1)(K + 1a + (a + b)(δ -1)) δ , C = Kab + (a -1)(δ -1) δ . This implies φ X (x) = Ax -B δ-1 log(1 + (δ -1)e -x ) + C log(1e -x ) + κ where κ is a constant. Thus one obtains p X (x) = De -(a+b)x (1 + (δ -1)e -x ) -λ-b (1e -x ) λ-1 1 x>0 with λ = C + 1 = B δ-1b = K+1-b-2a+aδ δ

  are independent if and only if the densities of Y and X are respectively

Acknowledgements. We are grateful to G. Letac for helpful discussions about this work, and to a referee whose comments led to an improvement of the paper.

Note that x → K 4 sinh 2 (x/2) solves (4.4) with the right-hand side equal to 0, and x → (c 1 -λ-1)e x +(c 1 +1)e -x 4 sinh 2 (x/2) is a particular solution of (4.4). Therefore, the solution of (4.4) is

for some constant K. This implies

As a consequence, there exists a constant δ such that

.

Equation (2.3), together with (1.20) and (1.21), imply, for u, v > 0,

.

Then we get that p (U,V ) (u, v) is the product of a function of u and a function of v and this gives item 2. of Theorem 1.4.