Antoine Liutkus

Christian Rohlfing

Antoine Deleforge

AUDIO SOURCE SEPARATION WITH MAGNITUDE PRIORS: THE BEADS MODEL

Keywords: audio probabilistic model, magnitude, phase, source separation 1

Audio source separation comes with the need to devise multichannel filters that can exploit priors about the target signals. In that context, experience shows that modeling magnitude spectra is effective. However, devising a probabilistic model on complex spectral data with a prior on magnitudes is non trivial, because it should both reflect the prior but also be tractable for easy inference. In this paper, we approximate the ideal donut-shaped distribution of a complex variable with approximately known magnitude as a Gaussian mixture model called BEADS (Bayesian Expansion Approximating the Donut Shape) and show that it permits straightforward inference and filtering while effectively constraining the magnitudes of the signals to comply with the prior. As a result, we demonstrate large improvements over the Gaussian baseline for multichannel audio coding when exploiting the BEADS model.

I. INTRODUCTION

Audio source separation aims at processing an audio mixture x composed of I channels (e.g. stereo I = 2) so as to recover its J constitutive sources s j [START_REF] Vincent | From blind to guided audio source separation: How models and side information can improve the separation of sound[END_REF]. It has many applications, including karaoke [START_REF] Liutkus | Adaptive filtering for music/voice separation exploiting the repeating musical structure[END_REF], upmixing [START_REF] Avendano | Frequency domain techniques for stereo to multichannel upmix[END_REF], [START_REF] Liutkus | Separation of music+ effects sound track from several international versions of the same movie[END_REF] and speech enhancement [START_REF] Barker | The third chime speech separation and recognition challenge[END_REF]. Usually, the sources are recovered by some time-varying filtering of the mixture, which is amenable to applying a I × I complex matrix G (f, t) to each entry x (f, t) of its Short-Term Fourier Transform (STFT) [START_REF] Vincent | Probabilistic modeling paradigms for audio source separation[END_REF], [START_REF] Benaroya | Audio source separation with a single sensor[END_REF], [START_REF] Avendano | Frequency domain techniques for stereo to multichannel upmix[END_REF]. Devising good filters requires either dedicated models for sources power spectrograms [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF], [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF] or machine learning methods to directly predict effective filters [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF].

The theoretical grounding for these linear filtering procedures boil down to considerations about second-order statistics for the sources STFT coefficients s j (f, t) ∈ C. When seen from a probabilistic perspective, this is often translated as picking a complex isotropic Gaussian distribution [START_REF] Gallager | Circularly Symmetric Complex Gaussian Random Vectors -A Tutorial[END_REF] This work was partly supported by the research programme KAMoulox (ANR-15-CE38-0003-01) funded by ANR, the French State agency for research. Fig. 1: The BEADS probabilistic model for a complex variable with a magnitude close to b and squared uncertainty σ.

Re

for each s j (f, t) ∈ C, which is called the Local Gaussian Model (LGM) [START_REF] Arberet | Underdetermined instantaneous audio source separation via local Gaussian modeling[END_REF], [START_REF] Duong | Under-determined convolutive blind source separation using spatial covariance models[END_REF] and is depicted on Figure 1a. Although it enjoys tractability and easy inference, the main shortcoming of the LGM is that it gives highest probability to 0, which may appear as counterintuitive. Still, this comes as a consequence of the fact that stationarity makes all phases to be equally probable a priori, so that E [s j (f, t)] = 0. Combining the two first moments only, maximum entropy principles naturally lead us to pick the LGM [START_REF] Jaynes | Probability theory: The logic of science[END_REF].

However, this badly reflects the additional prior knowledge one often has on s j (f, t). The vast majority of methods use priors on its magnitude: |s j (f, t)| should be close to some b j (f, t) > 0 with squared uncertainty σ j (f, t). In that setting, the donut-shaped distribution shown in Figure 1b is much better than LGM because it gives its highest probability mass on the circle of radius b j (f, t). If σ = 0, we end up with the phase unmixing problem [START_REF] Deleforge | Phase unmixing: Multichannel source separation with magnitude constraints[END_REF]. However, even with some uncertainty σ > 0, such a distribution suffers from non-tractability. In particular, it is not stable with respect to additivity, nor allows for simple posterior inference that would lead to straightforward filtering procedures.

In this paper, we translate our prior as a mixture of C identical Gaussian components evenly located over the circle of radius b j (f, t), yielding our proposed BEADS model (Bayesian Expansion Approximating the Donut Shape), depicted on Figure 1c. The most remarkable feature of BEADS is to allow both for straightforward filtering procedures while complying with priors on magnitude and possibly priors on phase. As such, it extends recent research on anisotropic modeling of complex spectra [START_REF] Magron | Phase-dependent anisotropic gaussian model for audio source separation[END_REF] by translating phase information into the choice of one particular component from the model. It thus appears as yet another way to incorporate phase information in source separation [START_REF] Sturmel | Phase-based informed source separation for active listening of music[END_REF].

We illustrate the BEADS model in an informed source separation (ISS) setting, where the true sources are available at a first coding stage, that allows to compute good models to be used for separation at a decoding stage [START_REF] Nikunen | Object-based Modeling of Audio for Coding and Source Separation[END_REF], [START_REF] Liutkus | Informed source separation : a comparative study[END_REF], [START_REF] Rohlfing | NMF-based informed source separation[END_REF]. As demonstrated already before [START_REF] Ozerov | Codingbased informed source separation: Nonnegative tensor factorization approach[END_REF], this ISS setup is interesting because the separation parameters can be encoded very concisely, leading to effective instances of spatial audio object coding [START_REF] Breebaart | Spatial audio object coding (SAOC)-the upcoming MPEG standard on parametric object based audio coding[END_REF].

II. PROBABILISTIC MODEL II-A. BEADS source model

The BEADS model is expressed as follows:

P [s j (f, t) | b j (f, t) , σ j (f, t)] = C c=1 π j (c | f, t) N (b j (f, t) ω c , σ j (f, t)) , (1)
where N denotes here the complex isotropic Gaussian distribution 1 , ω = exp (i2π/C) is the C th root of unity and π j (c | f, t) is the prior for the phase of source j at time-frequency (TF) bin (f, t): it indicates the probability that s j (f, t) is drawn from component c and hence that its phase is close to ω c . While some phase unwrapping approach [START_REF] Magron | Model-based STFT phase recovery for audio source separation[END_REF] may be used to set this prior, we take it as uniform here. The parameter σ j (f, t) stands for the variance of each component. It may be understood as the expected squared error of our estimate b j (f, t) for the magnitude. We can note that the LGM is equivalent to C = 1

1 π I |Σ| exp -(x -µ) Σ -1 (x -µ) ,
where |Σ| is the determinant of Σ [START_REF] Gallager | Circularly Symmetric Complex Gaussian Random Vectors -A Tutorial[END_REF]. and b j (f, t) = 0, and that many beads tend to the donut shape, as shown in Figure 1d. Now, we consider the joint prior distribution of the J independent sources. We need to consider all the C J possible combinations for the components. Let N C be the set of the first C natural numbers. We write z (f, t) ∈ N J C for the J × 1 vector whose j th entry z j (f, t) ∈ N C is the actual component drawn for source j. We define π (c | f, t) as the probability of each combination:

∀c ∈ N J C , π (c | f, t) def = P [z (f, t) = c] = J j=1 π j (c j | f, t) , (2)
where c j is the j th entry of c ∈ N J C . We have c π (c | f, t) = 1. The joint prior distribution of the sources is given by:

P [s (f, t) | Θ] = c∈N J C π (c | f, t) N (ω c • b (f, t) , [σ (f, t)]) , (3)
where ω c for c ∈ N J C denotes the J × 1 vector with entries ω cj , a•b denotes element-wise multiplication of the vectors a and b and [v] denotes the diagonal matrix whose diagonal is the vector v. In words, the prior distribution for the sources under the BEADS model is a Gaussian Mixture Model (GMM), with weights π (c | f, t) , which is reminiscent but different from the pioneering work in [START_REF] Rennie | Variational speech separation of more sources than mixtures[END_REF] that was limited to one source signal.

II-B. Mixture likelihood and separation

We take the mix as a convolutive mixture with the narrowband assumption, so that x (f, t) ≈ A (f ) s (f, t), where A (f ) is the I × J mixing matrix at frequency bin f . Further exploiting the BEADS model (1), we get the marginal model of the mixture given all parameters Θ as:

P [x (f, t) | Θ] = c∈N J C π (c | f, t) N (x c (f, t) , Σ x (f, t)) , (4) with x c (f, t) = A (f ) (ω c • b (f, t)) Σ x (f, t) = A (f ) [σ (f, t)] A (f ) . (5)
Now, the real advantages of the BEADS model is that we can straightforwardly obtain the joint posterior distribution of the sources as:

P [s (f, t) | x, Θ] = c∈N J C π (c | f, t, x) N (µ c , Σ c ) , (6)
where the posterior statistics µ c and Σ c for each combination c of the phases is:

µ c (f, t) = G (f, t) (x (f, t) -x c (f, t)) + ω c • b (f, t) Σ c (f, t) = Σ (f, t) = [σ (f, t)] -G (f, t) A (f ) [σ (f, t)] , (7)
Algorithm 1 BEADS decoder to update the phase configuration probabilities and perform separation. Input: parameters Θ, mixture x (f, t), prior π j , C.

Initialization: Σ π ← 0, ŝ ← 0 For all c n ∈ N J C : 1) π (c | f, t, x) ← π (c | f, t) N (x (f, t) | x c (f, t) , Σ x (f, t)) 2) ŝ (f, t) += π (c | f, t, x) (ω c • b (f, t) -x c (f, t)) 3) Σ π (f, t) += Σ π (f, t) + π (c | f, t, x) Finalization: 1) ŝ (f, t) ← ŝ(f,t) Σπ(f,t) 2) ŝ (f, t) += G (f, t) x (f, t)
with the J × I Wiener gain G (f, t) defined as:

G (f, t) = [σ (f, t)] A (f ) A (f ) [σ (f, t)] A (f ) -1 .
The minimum mean squared error (MMSE) estimate for the sources thus becomes:

E [s (f, t) | x, Θ] = G (f, t) x (f, t) + c∈N J C π (c | f, t, x) (ω c • b (f, t) -x c (f, t)) . (8)
In the case we know the true phases configuration z (f, t) or have estimated one at the decoder, this estimate (8) simplifies to ŝ (f, t) = µ z(f,t) .

III. PARAMETERS ESTIMATION

In this section, we only consider the informed case, where the true source signals s j (f, t) and the mixing matrices A (f ) are known at the coder, while the mixture and the parameters Θ = {b, σ, A} are known at the decoder.

III-A. Decoder: posterior π (c | f, t, x) and separation

We show how the probabilities π (c | f, t) for the phase configurations may be updated to yield their posterior π (c | f, t, x) that fully exploit the BEADS model for constraining the phases. Dropping the (f, t) indices for readability, we have:

∀c ∈ N J C , π (c | x) = π (c) N (x | x c , Σ x ) P [x] . (9)
This posterior probability may hence be expressed up to a normalizing constant independent of c as:

π (c | x) ∝ π (c | x) = π (c) exp -(x -x c ) Σ -1 x (x -x c ) , (10)
which can straightforwardly be computed at the decoder with known parameters Θ = {b, σ}. π (c | x) is obtained by normalization after computation for all c ∈ N C J . Algorithm 1 summarizes the computations done at the decoder to perform estimation of the posterior probabilities for the phase configurations and separation.

III-B. Coder: amplitudes b and errors σ

The parameters to be learned at the coder are the amplitude priors b j (f, t) and the error models σ j (f, t). First, for saving bitrate and computing time, we only use BEADS for the F 0 frequency bands that have the highest energy in the mix, let them be F 0 , and simply pick b j (f, t) = 0 for others. Then, for f ∈ F 0 we compress them by picking a Nonnegative Tensor Factorization model [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF]:

b j (f, t) = k W b (f, k) H (t, k) Q b (j, k) if f ∈ F 0 0 otherwise, ( 11
) in which case Θ b = {W b , H, Q b }
are small nonnegative F 0 × K, T × K and J × K, respectively. Same thing for σ, where we further reduce the number of parameters of the model by taking the same activations H in both cases, but this time, we model all frequencies.

In this section and for convenience, we will assume that the component z j (f, t) drawn for each source at each TF bin is known and equal to the one closest to s j (f, t) . This simplification has the advantage of strongly reducing the computational cost of the estimation algorithm. Indeed, the BEADS model ( 1) then reduces to:

s j (f, t) | z j ∼ N b j (f, t) ω zj (f,t) , σ j (f, t) .
We can define the relative source j (f, t) :

j (f, t) def = s j (f, t) ω zj (f,t) ∼ N (b j (f, t) , σ j (f, t)) .
(12) Provided z j (f, t) is correctly chosen as the component whose argument

2πzj (f,t) C
is the closest to that of s j (f, t) , we can furthermore safely assume that the real part R ( j (f, t)) of j (f, t) is nonnegative. Now, we detail the learning procedure we propose for b and σ. The strategy is to alternatively fix each of them and learn the other one.

Learning b j (f, t): Assume σ j (f, t) is kept fixed. The distribution [START_REF] Arberet | Underdetermined instantaneous audio source separation via local Gaussian modeling[END_REF] for the relative sources mean that we may estimate Θ b using a weighted Euclidean method:

Θ b ← argmin Θ f,t,j |R ( j (f, t)) -b j (f, t | Θ)| 2 σ j (f, t) ,
which can be done with a classical weighted NTF [START_REF] Virtanen | Combining pitch-based inference and non-negative spectrogram factorization in separating vocals from polyphonic music[END_REF] scheme with the Euclidean cost function.

Learning σ j (f, t): If b j (f, t) is fixed we see from (12) that j (f, t) -b j (f, t) has an isotropic complex Gaussian distribution with variance σ j (f, t). This means Θ σ can be estimated through:

Θ σ ← argmin Θ f,t,j d IS | j (f, t) -b j (f, t)| 2 σ j (f, t | Θ) ,
where d IS (a b) is the classical Itakura-Saito divergence for two nonnegative scalars a and b. This optimization is classical in the audio processing literature [START_REF] Fitzgerald | On the use of the beta divergence for musical source separation[END_REF], [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF], [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the beta-divergence[END_REF].

Considering that the activations H we take for σ are those for b, we only learn W σ and Q σ with it.

IV. EVALUATION

We evaluate the BEADS model through its performance for ISS, i.e. by displaying its average quality as a function of the bitrate required to transmit its parameters. To assess quality, we use BSSeval metrics [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]: SDR (Source to Distortion Ratio) and SIR (Source to Interference Ratio), both expressed in dB and higher for better separation. For normalization purpose, we compute δ-metrics, defined as the difference between the score and performance of oracle Wiener filtering, i.e. using true sources spectrograms [START_REF] Vincent | Oracle estimators for the benchmarking of source separation algorithms[END_REF].

The data consists of 10 excerpts of 30 s, taken from DSD100 database 2 . Each consists of J = 4 sources (vocals, bass, drums and accompaniment), sampled at 44.1 kHz. We generated either mono (I = 1) or stereo (I = 2) mixtures from these sources, through simple summation or anechoic mixing (delays+gains), respectively. STFT was conducted with 50 % overlap and a window size of 93 ms. We evaluated the following methods:

• BEADS oracle: ŝ (f, t) = µ z(f,t) .
• BEADS point using only the phase configuration ẑ that is most likely a posteriori. • BEADS as given in Algorithm 1.

• Itakura-Saito NTF [START_REF] Liutkus | Informed source separation through spectrogram coding and data embedding[END_REF], with K components.

Given all these methods and data, our extensive evaluation consisted in trying the methods with F 0 = 150 frequency bands for BEADS magnitudes, C = {8, 16} beads, and all methods were tried for K ∈ [8, 128] NTF components. We picked 16 quantization levels for all parameters. Results were smoothed using LOESS [START_REF] Cleveland | Locally weighted regression: an approach to regression analysis by local fitting[END_REF] and are displayed on Figure 2.

An interesting fact we see on Figure 2 is that the oracle BEADS model significantly outperforms standard oracle Wiener filtering, even for very crude magnitude models b j (f, t). This can be seen by the fact that its δmetrics get positive even at very small bitrates.

Then, we may notice that the δ-metrics appear as higher for mono than for stereo mixtures. In this respect, we should highlight that the absolute performance of oracle Wiener is of course higher for stereo (not shown on Figure 2), due to the knowledge of the mixing filters A (f ) ∈ C that alone bring good separation already and actually some information about the phase of the sources. Adding additional spectral knowledge in that case is then less important than in the mono case, where it is crucial. Now, we see a very clear improvement of BEADS as described in Algorithm 1 over classical NTF-ISS, of approximately 2 dB SDR and 5 dB SIR, at most bitrates. This significant boost in performance shows that BEADS helps a 2 http://sisec.inria.fr. lot in predicting the source signals by adequately handling priors on magnitudes, which is the main result for this study. Finally, Figure 2 also shows that the procedure for computing the phase posterior probabilities is not sufficient for correctly identifying the true phase configuration. This can be seen by the strong discrepancies between the BEADS point estimate and its oracle performance. While marginalization over the configuration as described in Algorithm 1 helps a lot in this respect, there is much room for improvement for parameter estimation of this model.

V. CONCLUSION

In this paper, we introduced BEADS as a convenient probabilistic model for complex data whose magnitude is approximately known. BEADS is a Gaussian Mixture Model where all components share the same variance and are scattered along a circle. While simple conceptually, BEADS comes with several advantages. First, it translates the delicate problem of modeling the phase into setting probabilities over a discrete set of components. Second, it allows for easy inference and, finally, it straightforwardly leads to effective filtering procedures. Although we demonstrated its potential in an audio-coding application, we believe it may also be useful in the blind separation setting when embedded in an Expectation-Maximization estimation procedure.

  The BEADS model combines advantages of both (C = 8)... ...and approximates the donut shape well (C = 16).

16 Fig. 2 :

 162 Fig.2: BEADS for ISS on mono (left) or stereo (right) mixes. Metrics are δSDR (top) and δSIR (bottom). Units are kilobits/second/source (x-axis) and dB (y-axis).


El Haj Laamri
email: el-haj.laamri@univ-lorraine.fr

Michel Pierre
email: michel.pierre@ens-rennes.fr

Stationary reaction-diffusion systems in L 1

Keywords: reaction-diffusion systems, nonlinear diffusion, cross-diffusion, global existence, porous media equation, weak solutions. 2010 MSC: 35K10, 35K40, 35K57

come

Introduction and main results

Our goal is to analyze the existence of solutions to stationary reaction-diffusion systems in L 1 -spaces. We first give a general abstract result for nonlinearities satisfying as many structural inequalities as the number of equations. This is done in the framework of abstract m-accretive operators in L 1 spaces for the diffusion part and the full system is controlled by a somehow "good" associated cross-diffusion system. We give several examples where this abstract result applies. Then, we provide a general result for more specific systems associated with general chemical reactions for which less structure holds for the nonlinearities.

We denote by (Ω, µ) a measured space where µ is a nonnegative measure on the set Ω with µ(Ω) < +∞. We consider systems of the type (S)

   ∀i = 1, ..., m, u i ∈ D(A i ) ∩ L 1 (Ω, dµ) + , h i (•, u) ∈ L 1 (Ω, dµ), u i + A i u i = h i (•, u 1 , ..., u m ) + f i (•) ∈ L 1 (Ω, dµ), (1)
where for all i = 1, ..., m, -A i is a (possibly) nonlinear operator in L 1 (Ω, dµ) (generally a diffusion operator in applications), defined on D(A i ) ⊂ L 1 (Ω, dµ), -h i : Ω × R m → R is a nonlinear "reactive" term, -f i ∈ L 1 (Ω, dµ) + [:= {g ∈ L 1 (Ω, dµ) ; g ≥ 0, µ -a.e.}].

We are interested in nonnegative solutions u = (u 1 , ..., u m ) ∈ L 1 (Ω, dµ) +m .

These systems are naturally associated with the evolution reaction-diffusion systems ∂ t u i (t) + A i u i (t) = h i (•, u(t)), t ∈ [0, +∞), i = 1, ..., m.

(

When approximating these evolution systems by an implicit time discretization scheme, we are led to solving the following set of equations, for all small time interval ∆t ∀ i = 1, ..., m, n = 0, 1, ..., u n+1 i ∈ D(A i ),

u n+1 i -u n i ∆t + A i u n+1 i = h i (•, u n+1 ) or u n+1 i + (∆t)A i u n+1 i = (∆t)h i (•, u n+1 i ) + u n i .
This is exactly the system (S) with unknown u = u n+1 , up to trivially changing (∆t)A i into A i and (x, u) → [(∆t)h i (x, u) + u n i (x)] 1≤i≤m into (x, u) → [h i (x, u) + f i (x)] 1≤i≤m . The asymptotic steady states of the evolution system (2) are also the solutions u = (u 1 , ..., u m ) to the system {A i u i = h i (•, u), i = 1, ..., m, } which is also essentially included in the system (S) up to a slight change of the operators A i . The kind of operators A i we have in mind are : -The Laplacian u → -∆u on a bounded open set Ω of R N with various boundary conditions on ∂Ω.

-More general elliptic operators u → -N i=1 ∂ x i N i=1 a ij (x)∂ x j u + b i u on the same Ω, with various boundary conditions as well.

-Nonlinear operators of porous media type like u → -∆ϕ(u) where ϕ : R → R is an increasing function, again with different boundary conditions on ∂Ω.

-Nonlinear operators of p-Laplacian type like u → -∆ p u := -∇ • |∇u| p-2 ∇u where p ∈ (1, +∞) and | • | is the euclidian norm in R N .

-Classical perturbations and various associations of all of those. All these operators will satisfy the assumption (A) below which means that each A i is an m-accretive operator in L 1 (Ω, dµ) whose resolvents are compact and preserve positivity, namely the following, where I denote the identity : 1 (Ω, dµ), and for all λ ∈ (0, +∞), (A1) maccretivity : I + λA i is onto and (I + λA i ) -1 is nonexpansive, (A2) positivity : Ω sign -(u i )Au i ≥ 0, ∀ u i ∈ D(A i ), (A3) compactness : (I + λA i ) -1 : L 1 (Ω, dµ) → L 1 (Ω, dµ) is compact.

(A)            ∀ i = 1, ..., m, A i : D(A i ) ⊂ L 1 (Ω, dµ) → L
(3)

Here I denotes the identity on L 1 (Ω, dµ) and in (A 2 ), we define sign -(s) := -sign + (-s), ∀s ∈ R where sign + (s) := 1, ∀ s ∈ [0, +∞), sign + (s) = 0, ∀ s ∈ (-∞, 0). This implies immediately that (I + λA i ) -1 L 1 (Ω, dµ) + ⊂ L 1 (Ω, dµ) + since then, for g ∈ L 1 (Ω, dµ) + and u i := (I + λA i ) -1 g

0 ≥ Ω sign -(u i )g dµ = Ω sign -(u i )[u i + λA i u i ] dµ ≥ Ω u - i dµ ⇒ u - i = 0 µ -a.e.
For simplicity, we consider only single-valued operators A i , but everything would also work for multivalued m-accretive operators A i .

The nonlinear reactive terms h i will preserve positivity as well. We assume that, for all i = 1, ..., m :

(H1)

       (1) h i : Ω × R m → R measurable ; (2) h i (•, R) := sup{|h i (•, r)|; |r| ≤ R} ∈ L 1 (Ω, dµ) + , ∀ R ∈ [0, +∞) ; (3) r ∈ R m → h i (•, r) is continuous ; (4) quasipositivity : h i (•, r 1 , ..., r i-1 , 0, r i+1 , ..., r m ) ≥ 0, ∀ r ∈ [0, +∞) m . (4)
Moreover, the nonlinear reactive terms h i we are considering are those for which mass conservation or, more generally, mass dissipation or at least mass control holds for the associated evolution system [START_REF] Bebernes | Finite time blowup for semilinear reactive-diffusive systems[END_REF]. It is the case when the h i 's satisfy a relation like

m i=1 h i (•, r) ≤ 0 for all r ∈ [0, +∞) m or more generally (H2) m i=1 a i h i (•, r) ≤ m i=1 b i r i +ω(•), ω ∈ L 1 (Ω, dµ) + , for some 0 ≤ b i < a i , i = 1, ..., m, r ∈ [0, +∞) m . ( 5
)
As we will see, this structure naturally implies an a priori L 1 -estimate on the solutions u i of the system (S) (see Proposition 2.2), together with "standard operators" A i by which we mean

(A inf ) a ∞ := inf Ω A i u i dµ, u i ∈ D(A i ) ∩ L 1 (Ω, dµ) + > -∞ for each i = 1, ..., m. (6)
This assumption essentially holds when 0 ∈ D(A i ), i = 1, ..., m. It also holds with most diffusion operators with non homogeneous boundary conditions (see Section 3) except a few ones : this is discussed in Remark 3.4.

Actually, the main point will be to get a priori L 1 -estimates even on the nonlinearities h i (•, u) where u is solution of System (S). This will be satisfied if more structure is required on the nonlinear functions h i , namely (using the natural order in R m ) :

(M )    There exist two m × m matrices M 0 , M 1 where M 0 is invertible, with nonnegative entries and such that M 0 h(•, r) ≤ M 1 r + Θ(•), ∀ r ∈ [0, +∞) m , Θ ∈ L 1 (Ω, dµ) +m . (7)
Applying the matrix M 0 to the system (S) leads to the following set of m inequalities for an associated cross-diffusion system :

M 0 u + M 0 Au = M 0 h(•, u) + M 0 f ≤ M 1 u + Θ + M 0 f, where we denote Au := (A i u i ) 1≤i≤m , h(•, u) := (h i (•, u)) 1≤i≤m and f := (f i ) 1≤i≤m .
As proved in Proposition 2.3, this will imply an a priori estimate of h i (•, u), i = 1, ..., m in L 1 (Ω, dµ). More precisely, we will consider an approximate problem where the nonlinearities h i are replaced by truncated versions h n i . Then (M ) will imply that h n i (•, u n ) is bounded in L 1 (Ω, dµ) for the approximate solutions u n . The compactness of u n in L 1 (Ω, dµ) will then follow. Thus up to a subsequence, u n will converge in L 1 (Ω, dµ) m and µ-a.e. to some u and h n i (•, u n ) will also converge µ-a.e. to h i (•, u) for all i. But this is not sufficient yet to pass to the limit in the system since one essentially needs the convergence of

h n i (•, u n ) in L 1 (Ω, dµ). It is the case if h n i (•, u n
) is uniformly integrable, by Vitali's lemma. This will hold if we add the following technical assumption which is some kind of compatibility condition between the various operators A i . It will be satisfied in the examples of Section 3.

(Φ)

There exist ϕ : [0, +∞) m → [0, +∞) continuous with lim |r|→+∞ ϕ(r) = +∞ and b ∈ (0, +∞) m such that

Ω sign + (ϕ(u) -k)b • M 0 Au dµ ≥ 0, ∀ u ∈ D(A) ∩ L 1 (Ω, dµ) +m , ∀ k ∈ [0, +∞), D(A) := Π m i=1 D(A i ). ( 8
)
We now state our abstract result. It will be proved in Section 2 and applied to several explicit examples in Section 3. We refer to [START_REF] Laamri | Existence globale pour des systèmes de réaction-diffusion dans L 1[END_REF] where this kind of results and examples were already widely discussed and analyzed.

Theorem 1.1 Assume that (A), (H1), (H2), (A inf ), (M ), (Φ) hold. Then the system (S) has a solution for all f ∈ L 1 (Ω, dµ) +m .

As already explained, the assumption (M ) means that m independent inequalities hold between the m nonlinear functions h i . But many systems come with less than m such relations. The following 2 × 2 system satisfies only one relation of type (H2) and not (M ):

   u 1 -d 1 ∆u 1 = (β 1 -α 1 )[u α 1 1 u α 2 2 -u β 1 1 u β 2 2 ] + f 1 =: h 1 (u 1 , u 2 ) + f 1 , u 2 -d 2 ∆u 2 = (β 2 -α 2 )[u β 1 1 u β 2 2 -u α 1 1 u α 2 2 ] + f 2 =: h 2 (u 1 , u 2 ) + f 2 , ∂ ν u 1 = 0 = ∂ ν u 2 on ∂Ω, (9)
where for i = 1, 2,

d i ∈ (0, +∞), α i , β i ∈ {0} ∪ [1, +∞), f i ∈ L 1 (Ω) + , Ω ⊂ R N equipped
with the Lebesgue measure. This kind of nonlinearity appears in reversible chemical reactions with two species, namely

α 1 A 1 + α 2 A 2 β 1 A 1 + β 2 A 2 . ( 10
)
Here (β 1 -α 1 )(β 2 -α 2 ) < 0. The nonlinearity h 1 and h 2 are quasipositive but satisfy the only one relation

γ 2 h 1 + γ 1 h 2 = 0, γ i := |β i -α i |.
It turns out that existence indeed holds for system (9) when f i ∈ L 1 (Ω) + for all i, as we prove it here. Actually we prove existence for a general class of such "chemical systems" when f i ∈ L 1 (Ω) + and also f i log f i ∈ L 1 (Ω). This is the second main result of this paper. Let us consider the system (CHS)

     For all i = 1, ..., m, u i -d i ∆u i = (β i -α i ) k 1 Π m k=1 u α k k -k 2 Π m k=1 u β k k + f i , ∂ ν u i = 0 on ∂Ω, (11)
where k 1 , k 2 ∈ (0, +∞) and, for all i = 1, ..., m, d i ∈ (0, +∞), α i , β i ∈ {0} ∪ [1, +∞) and f i ∈ L 1 (Ω) + where Ω is bounded regular open subset of R N and where    I := i ∈ {1, ..., m} ; α i -β i > 0 , J := j ∈ {1, ..., m}; β j -α j > 0 satisfy :

I = ∅, J = ∅, I ∪ J = {1, ..., m}. (12)
We denote by |I| (resp. |J|) the number of elements of I (resp. J).

Theorem 1.2 Assume that f i ∈ L 1 (Ω) + , f i log f i ∈ L 1 (Ω) for all i = 1, ..., m. Assume also that |I| ≤ 2 (or |J| ≤ 2). Then there exists a nonnegative solution u ∈ W 2,1 (Ω) +m of (CHS) with h i (u) ∈ L 1 (Ω) for all i = 1, ..., m. If moreover m = 2, then the same result holds if

f i ∈ L 1 (Ω) + , i = 1, 2 only.
Remark 1.3 System (11) arises when modeling a general reversible chemical reaction with m species, according to the mass action law and to Fick's linear diffusion. In fact, it also contains systems written more generally as

v i -d i ∆v i = λ i K 1 Π m k=1 v α k k -K 2 Π m k=1 v β k k + f i ,
where λ i ∈ R, λ i (β i -α i ) > 0, i = 1, ..., m. Indeed, we may go back to the exact writing of (11) by setting

u i := (β i -α i )v i /λ i , k 1 := K 1 Π k [λ k /(β k -α k )] α k , k 2 := K 2 Π k [λ k /(β k -α k )] β k , f i := (β i -α i ) f i /λ i .
On the other hand, it does not include systems like

u 1 -d 1 ∆u 1 = +[u 3 1 u 2 2 -u 2 1 u 3 2 ] + f 1 (= u 2 1 u 2 2 [u 1 -u 2 ] + f 1 ), u 2 -d 2 ∆u 2 = -[u 3 1 u 2 2 -u 2 1 u 3 2 ] + f 2 .
Here the condition λ i (β i -α i ) > 0 is not satisfied. Remark 1.4 Note that besides the 2 × 2 system (9), Theorem 1.2 contains as particular cases some favorite systems of the literature like

h i (u) = (-1) i [u α 1 1 u α 3 3 -u β 2 2 ], i = 1, 2, 3, α 1 , α 3 , β 2 ∈ [1, +∞), h i (u) = (-1) i [u 1 u 3 -u 2 u 4 ], i = 1, 2, 3, 4.
We may use Remark 1.3 to write them as in [START_REF] Laamri | Existence globale pour des systèmes de réaction-diffusion dans L 1[END_REF]. Analysis of systems of this kind may be found for instance in [START_REF] Laamri | Global existence of classical solutions for a class of reaction-diffusion systems[END_REF], [START_REF] Desvillettes | Global existence for quadratic systems of reactiondiffusion[END_REF], [START_REF] Pierre | Global existence for a class of quadratic reaction-diffusion systems with nonlinear diffusions and L 1 initial data[END_REF], [START_REF] Laamri | Global existence for reaction-diffusion systems with nonlinear diffusion and control of mass[END_REF], [START_REF] Goudon | Regularity analysis for systems of reaction-diffusion equations[END_REF], [START_REF] Cañizo | Improved duality estimates and applications to reaction-diffusion equations[END_REF], [START_REF] Caputo | Solutions of the 4-species quadratic reaction-diffusion systems are bounded and C ∞ -smooth[END_REF], etc.

In fact Theorem 1.2 applies to quite general systems like, for instance, those obtained by multiplying the above 3 × 3 (resp. 4 × 4) example by

u σ 1 1 u σ 2 2 u σ 3 3 (resp. u σ 1 1 u σ 2 2 u σ 3 3 u σ 4 4
) where σ i ∈ [0, +∞). They lead to the following nonlinearities

h i (u) = λ i [Π m k=1 u α k k -Π m k=1 u β k k ]
, where m = 3 (resp. m = 4) and λ i (β i -α i ) > 0. Actually, Theorem 1.2 applies also to these same nonlinearities when m = 5. Indeed, since I ∪ J = {1, 2, 3, 4, 5}, it follows that I or J contains at most two elements.

We believe that the result of Theorem

:= |x|, x ∈ Ω := B(0, 1) ⊂ R N , b ∈ (N/2, N -2), c ∈ (0, +∞): σ(r) := r -b + br + c, ∀ 0 < r ≤ 1, ⇒ σ -∆σ = f, f (r) := c + r -b + br + b(N -2 -b)r -(b+2) -b(N -1)r -1 , σ (1) = 0,
where c is chosen large enough so that f ≥ 0. Note that f ∈ L p (Ω), p ∈ [1, N/(b + 2)). Let us now consider the solution (which exists by Theorem 1.2) of the system

       u 1 , u 2 ∈ W 2,1 (Ω), u 2 2 -u 2 1 ∈ L 1 (Ω), u 1 -∆u 1 = u 2 2 -u 2 1 + f /2, u 2 -∆u 2 = -[u 2 2 -u 2 1 ] + f /2, ∂ ν u 1 = ∂ ν u 2 = 0 on ∂Ω. Then (u 1 + u 2 ) -∆(u 1 + u 2 ) = f, ∂ ν (u 1 + u 2 ) = 0 on ∂Ω. Thus u 1 + u 2 = σ. But u 1 , u 2 cannot be in L 2 (Ω) since σ is not in L 2 (Ω) by the choice of b > N/2. Remark 1.6 It is classical that an entropy structure holds in System (11) since m i=1 [log u i + µ i ]h i (u) = -[log k 1 Π m k=1 u α k k -log k 2 Π m k=1 u β k k ][k 1 Π m k=1 u α k k -k 2 Π m k=1 u β k k ] ≤ 0, when choosing µ i := [log k 1 -log k 2 ]/[m(α i -β i )].
Under the assumptions of Theorem 1.2, and in particular the LLogL assumption on the f i , we have the estimate

Ω m i=1 [log u i + µ i ]u i + d i |∇u i | 2 u i ≤ Ω m i=1 [log u i + µ i ]f i ≤ Ω m i=1 f i log f i + (µ i -1)f i + u i ,
where, for the last inequality, we use the Young's inequality (75) with r = f i , s = log u i . But we will not use here this stucture in the proof of Theorem 1.2. Our strategy will consist in proving that the nonlinearity h i (u) is a priori bounded in L 1 (Ω). Then adequate compactness arguments allow us to pass to the limit in the approximate system. Note that the entropy inequality provides the extra information that ∇ √ u i ∈ L 2 (Ω) for the solutions obtained in Theorem 1.2 when

f i log f i ∈ L 1 (Ω), i = 1, ..., m.
Theorem 1.1 will be proved in Section 2, examples of applications are given in Section 3 and the proof of Theorem 1.2 is given in Section 4.

Proof of Theorem 1.1

Since µ is fixed, we will most of the time more simply write L 1 (Ω), L 1 (Ω) + instead of L 1 (Ω, dµ), L 1 (Ω, dµ) + . We will use the natural order in R m namely [r ≤ r ⇔ r i ≤ ri , i = 1, ..., m] and we also denote r

+ := (r + 1 , ..., r + m ). Lemma 2.1 Let f = (f 1 , ..., f m ) ∈ L 1 (Ω) +m . We set h n i (x, r) := h i (x, r) 1 + n -1 m j=1 |h j (x, r)| , ∀ i = 1, ..., m, r ∈ R m , x ∈ Ω. ( 13
)
Assume that (A) and (H1) hold. Then the following approximate system

(S n ) ∀i = 1, ..., m, u n i ∈ D(A i ) ∩ L 1 (Ω) + , u n i + A i u n i = h n i (•, u n ) + f i (•) in L 1 (Ω, dµ), (14)
has a nonnegative solution

u n = (u n 1 , • • • , u n m ). Proof. We consider the mapping T : v ∈ L 1 (Ω) m → w ∈ L 1 (Ω) m where w = (w 1 , ..., w m ) is the solution of ∀i = 1, ..., m, w i ∈ D(A i ), w i + A i w i = h n i (•, v + ) + f i (•) in L 1 (Ω, dµ), (15)
where v + := (v + 1 , ..., v + m ). The mapping T is well defined since |h n i (x, r)| ≤ n for all (x, r) ∈ Ω × R m and f ∈ L 1 (Ω) m so that the solution is given by

w i = (I + A i ) -1 h n i (•, v + ) + f i , i = 1, ..., m. Moreover h n i (•, v) + f i L 1 (Ω) ≤ n µ(Ω) + f i L 1 (Ω) =: M i . Let K i := (I + A i ) -1 {g ∈ L 1 (Ω) ; g L 1 (Ω) ≤ M i } , K := K 1 × ... × K m ⊂ L 1 (Ω) m .
By assumption (A) and in particular (A3), K is a compact set of L 1 (Ω) m and K ⊃ T L 1 (Ω) m . On the other hand, we easily check that T is continuous. Thus T is a compact operator from L 1 (Ω) m into the compact set K ⊂ L 1 (Ω) m . By Schauder's fixed point theorem (see e.g. [START_REF] Zeidler | Nonlinear Functional Analysis and its applications, 1. Fixed Point Theorems[END_REF]), T has a fixed point u n . This means that u n is solution of (S n ).

To prove the nonnegativity property of u n , we multiply the i-th equation by sign -(u n i ) and integrate to obtain, using

Ω sign -(u i )A i u i ≥ 0 (see (A2) ), Ω (u n i ) -dµ ≤ Ω sign -(u n i )h n i (•, (u n ) + ) + sign -(u n i )f i dµ.
By the nonnegativity of f i and the quasipositivity of h i assumed in (H1), the integral on the right is nonpositive so that

Ω (u n i ) -dµ ≤ 0 or u n i ≥ 0 µ -a.e..
We will now progressively obtain estimates on u n independently of n. We start with the control of the total mass of u n .

Proposition 2.2 Assume that (A), (H1), (H2), (A inf ) hold. Then there exists C ∈ (0, +∞) independent of n such that the solution u n of the approximate system (S n ) satisfies

max 1≤i≤m u n i L 1 (Ω) ≤ C.
Proof. According to (H2), we multiply each equation by a i and we add them to obtain

m i=1 a i (u n i + A i u n i ) = m i=1 a i [h n i (•, u n ) + f i ] ≤ m i=1 [b i u n i + a i f i ]. We use (A inf ) to obtain n i=1 (a i -b i ) Ω u n i dµ ≤ -a ∞ m i=1 a i + m i=1 a i Ω f i dµ.
Using a i > b i for all i yields the estimate of Proposition 2.2.

We now prove that the nonlinearities are bounded in L 1 independently of n.

Proposition 2.3 Assume that (A), (H1), (H2), (A inf ), (M ) hold. Then there exists C ∈ (0, +∞) independent of n such that the solution u n of the approximate system (S n ) satisfies

max 1≤i≤m h n i (•, u n ) L 1 (Ω) ≤ C.
Proof. According to the assumption (M ), let us multiply the system (S n ) by the matrix M 0 . As before, we denote

Au n := (A 1 u n , ..., A m u n ), h n (•, u n ) := (h n 1 (•, u n ), ..., h n m (•, u n )) , f := (f 1 , ..., f m ).
Then

M 0 u n + M 0 Au n = M 0 h n (•, u n ) + M 0 f ≤ M 1 u n + Θ + M 0 f. (16)
Since the entries of M 0 are nonnegative, and thanks to (A inf ) and the nonnegativity of u n , there exists

C ∞ ∈ [0, +∞) m such that Ω [M 0 u n + M 0 Au n ]dµ ≥ -C ∞ which implies Ω [M 0 h n (•, u n ) + M 0 f ]dµ ≥ -C ∞ .
We deduce

Ω [M 1 u n + Θ -M 0 h n (•, u n )]dµ ≤ C ∞ + Ω [M 1 u n + Θ + M 0 f ]dµ ≤ D ∞ ∈ (0, +∞) m ,
the last inequality using also Proposition 2.2. As the function

M 1 u n + Θ -M 0 h n (•, u n ) is nonnegative, these inequalities provide a bound for its L 1 (Ω) m -norm. It follows that M 0 h n (•, u n ) is also bounded in L 1 (Ω) m independently of n. Since M 0 is invertible, we deduce that h n (•, u n ) is itself bounded in L 1 (Ω) m : indeed if | • | denotes the euclidian norm in R m and
• the induced norm on the m × m matrices, we may write

|h n (•, u n )| = |M -1 0 M 0 h n (•, u n )| ≤ M -1 0 |M 0 h n (•, u n )|.
Whence Proposition 2.3.

Proposition 2.4 Assume that (A), (H1), (H2), (A inf ), (M ), (Φ) hold. Then, if u n is the solution of the approximate system (S n ), {h n (., u n )} n is uniformly integrable in Ω which means that, for all ε > 0, there exists δ > 0 such that for all measurable set E ⊂ Ω

µ(E) < δ ⇒ E n i=1 |h n i (•, u n )|dµ ≤ ε for all n.
Proof. By Proposition 2.3, we already know that h n (•,

u n ) is bounded in L 1 (Ω) m .
Let us prove that the extra condition (Φ) implies that it is even uniformly integrable. First, since

u n i = (I + A i ) -1 (h n (•, u n ) + f i ), the compactness condition (A3) in (3) implies that, for all i = 1, ..., m, u n i belongs to a compact set of L 1 (Ω). Let us show that M 0 h n (•, u n ) is uniformly integrable. Since M 0 is invertible, this will imply that h n (•, u n )
is itself uniformly integrable and end the proof of Proposition 2.4.

We will successively prove that (M

0 h n (•, u n )) + and (M 0 h n (•, u n )) -are uniformly integrable.
Note that by the definition (13), h n (•, u n ) ≤ h(•, u n ). We may write (recall that the entries of M 0 are nonnegative)

M 0 h n (•, u n ) ≤ M 0 h(•, u n ) ≤ M 1 u n + Θ ⇒ (M 0 h n (•, u n )) + ≤ (M 1 u n ) + + Θ. ( 17
) Since u n is in a compact set of L 1 (Ω) m , so is (M 1 u n ) +
and it is in particular uniformly integrable. We then deduce from the last inequality that (M 0 h n (•, u n )) + is uniformly integrable.

To control (M 0 h n (•, u n )) -, we go back to ( 16) and, using ( 17), we rewrite it as follows

M 0 u n + M 0 Au n + (M 0 h n (•, u n )) -= (M 0 h n (•, u n )) + + M 0 f ≤ (M 1 u n )) + + Θ + M 0 f. ( 18
)
According to (Φ), we multiply this inequality by sign + (ϕ(u n ) -k)b. By (Φ) on one hand and by the nonnegativity of u n , b and of the entries of M 0 on the other hand, we have

Ω sign + (ϕ(u n ) -k)b • M 0 Au n dµ ≥ 0, Ω sign + (ϕ(u n ) -k)b • M 0 u n dµ ≥ 0, ∀ k ∈ [0, +∞).
Combining with [START_REF] Pierre | Global existence for a class of quadratic reaction-diffusion systems with nonlinear diffusions and L 1 initial data[END_REF], we deduce

Ω sign + (ϕ(u n ) -k)b • (M 0 h n (•, u n )) -dµ ≤ Ω sign + (ϕ(u n ) -k)b • (M 1 u n )) + + Θ + M 0 f dµ. (19) Since lim |r|→+∞ ϕ(r) = +∞, for all k ∈ [0, +∞), there exists R(k) such that [|r| ≥ R(k)] ⊂ [ϕ(r) ≥ k] and therefore [|u n | ≥ R(k)] ⊂ [ϕ(u n ) ≥ k] , µ -a.e.. (20)
Let E ⊂ Ω be a measurable set. Then

E b • (M 0 h n (•, u n )) -dµ ≤ E∩[|u n |≤R(k)] b • (M 0 h n (•, u n )) -dµ + E∩[|u n |≥R(k)] b • (M 0 h n (•, u n )) -dµ =: I n -+ I n + . If we denote ϕ R := max{ϕ(r); |r| ≤ R}, then R ∈ [0, +∞) → ϕ R ∈ [0, +∞) is a nondecreasing function such that [ϕ(r) ≥ k] ⊂ [|r| ≥ ψ(k)] where ψ(k) := inf ϕ -1 R ([k, +∞))
. By [START_REF] Schmitt | Existence globale ou explosion pour les systèmes de réaction-diffusion avec contrôle de masse[END_REF] and [START_REF] Rothe | Global solutions of reaction-diffusion systems[END_REF] and

[ϕ(u n ) ≥ k] ⊂ [|u n | ≥ ψ(k)],
we have

I n + ≤ [ϕ(u n )≥k] b • (M 1 u n ) + + Θ + M 0 f dµ ≤ [|u n |≥ψ(k)] b • (M 1 u n ) + + Θ + M 0 f dµ. ( 21
)
Since u n lies in a compact set of L 1 (Ω) m and lim

k→+∞ ψ(k) = +∞, then lim k→+∞ µ([|u n | ≥ ψ(k)]) = 0 uniformly
in n. Thus, given ε ∈ (0, 1), there exists k = k ε large enough so that

I n + ≤ ε/2 for all n. (22)
We can also control I n -as follows. We remark that (see assumption [START_REF] Bebernes | Finite time blowup for semilinear reactive-diffusive systems[END_REF] in (H1))

|u n | ≤ R(k ε ) ⇒ |h n i (•, u n )| ≤ |h i (•, u n )| ≤ h i (•, R(k ε )).
This implies that for some B ∈ (0, +∞),

I n -= E∩[|u n |≤R(k)] b • (M 0 h n (•, u n )) -dµ ≤ E B m i=1 h i (•, R(k ε )) (23)
We may choose δ small enough (independent of n) such that

µ(E) ≤ δ ⇒ E m i=1 h i (•, R(k ε )) < ε/2B.
Combining with

I n + ≤ ε/2 proved above (see (22)), we deduce that b • (M 0 h n (•, u n )) -is uniformly integrable. Since b ∈ (0, +∞) m , this implies that (M 0 h n (•, u n )) -is itself uniformly integrable. We already know that (M 0 h n (•, u n )) + is uniformly integrable. Thus so is M 0 h n (•, u n )
and this ends the proof of Proposition 2.4.

Proof of Theorem 1.1. We consider the solution u n of the approximate system (S n ) built in Lemma 2.1. By Propositions 2.2, 2.3, 2.4, up to a subsequence as n → +∞, we may assume that u n converges in L 1 (Ω) m and µ-a.e. to some u ∈ L 1 (Ω) +m . Moreover, by definition of h n i (see [START_REF] Laamri | Global existence for reaction-diffusion systems with nonlinear diffusion and control of mass[END_REF] ) and the continuity property of h i assumed in (H1), h n i (•, u n ) converges µ-a.e. to h(•, u). Moreover h n (•, u n ) is uniformly integrable. By Vitali's Lemma (see e.g. [START_REF] Schmitt | Existence globale ou explosion pour les systèmes de réaction-diffusion avec contrôle de masse[END_REF] or [START_REF] Fonseca | Modern Methods in the Calculus of Variations : L p spaces[END_REF]),

h n (•, u n ) converges also in L 1 (Ω) m to h(•, u). Since u n i = (I +A i ) -1 (h n (•, u n )+f i ) this implies that u i = (I + A -1 i (h(•, u) + f i )
which means that u is solution of the limit system (S) and this ends the proof of Theorem 1.1.

Examples

In all examples below, Ω is a bounded open subset of R N with regular boundary and equipped with the Lebesgue measure.

Examples with linear diffusions and homogeneous boundary conditions

We start with a simple example associated with the Laplacian operator with homogeneous boundary conditions.

Corollary 3.1 Assume the nonlinearity h = (h 1 , ..., h m ) satisfies (H1), (H2), (M ). Let d i ∈ (0, +∞), f i ∈ L 1 (Ω) + , i = 1, ..., m. Then the following system has a solution

   for all i = 1, ..., m; u i ∈ W 1,1 0 (Ω) + , h i (•, u) ∈ L 1 (Ω), u i -d i ∆u i = h i (•, u) + f i . (24)
Proof. We consider the operators

D(A i ) := {u ∈ W 1,1 0 (Ω) ; ∆u ∈ L 1 (Ω)}, A i u := d i ∆u. ( 25
)
It is classical that these operators A i satisfy the three conditions of (A) in (3) (see e.g. [START_REF] Brezis | Semilinear elliptic equations in L 1[END_REF]). Since 0 ∈ D(A i ), as already noticed just after ( 6), (A inf ) is also satisfied. And for (Φ), if we denote

M 0 = [m ij ] 1≤i,j≤m , we choose ϕ(r) := m j=1 ( m i=1 m ij ) d j r j , b = (1, ..., 1) t ∈ (0, +∞) m . Note that m i=1 m ij > 0 for all j = 1, ..., m since m ij ≥ 0 and M 0 is invertible. In particular, lim |r|→+∞ ϕ(r) = +∞). Then, if u ∈ D(A), k ∈ (0, +∞), Ω sign + (ϕ(u) -k)b • M 0 Au = - Ω sign +   m i,j=1 m ij d j u j -k   ∆   m i,j=1 m ij d j u j   ≥ 0. ( 26
)
Then we may apply Theorem 1.1.

Remark 3.2

As examples of functions h, we may for instance choose

   m = 2, α i , β i ∈ {0} ∪ [0, +∞), i = 1, 2, λ ∈ (0, 1) h 1 (•, u 1 , u 2 ) = λu α 1 1 u α 2 2 -u β 1 1 u β 2 2 , h 2 (•, u 1 , u 2 ) = -u α 1 1 u α 2 2 + u β 1 1 u β 2 2 .
We easily check that (H1) holds and that (H2) is satisfied with a i = 1, b i = 0 for all i and ω = 0 (in other words h 1 + h 2 ≤ 0). Moreover (M ) is satisfied with M 1 = 0, Θ = (0, 0) t and

M 0 = 1 1 1 λ .
Note that for λ = 1, M 0 is not invertible so that only one relation h 1 + h 2 ≤ 0 holds and Theorem 1.1 does not apply. This kind of systems is considered in Theorem 1.2.

• Here is another example of a nonlinearity h = (h 1 , h 2 , h 3 ) which satisfies the corollary.

m = 3, α i ∈ [1, +∞), i = 1, 2, h 1 = u 3 -u α 1 1 u α 2 2 = h 2 = -h 3 .
The corresponding evolution problem is studied in [START_REF] Rothe | Global solutions of reaction-diffusion systems[END_REF] for N ≤ 5, [START_REF] Martin | Nonlinear reaction-diffusion systems[END_REF] for any dimension N with α 1 = α 2 = 1, and in [START_REF] Laamri | Global existence of classical solutions for a class of reaction-diffusion systems[END_REF] where (α

1 , α 2 ) ∈ [1, +∞) 2 .
Here (H1) is obviously satisfied and so is (H2) with

a i = 1, i = 1, 2, a 3 = 2, b i = 0, i = 1, 2, 3, ω = 0. Then (M ) is satisfied with Θ = (0, 0, 0) t and M 0 =   1 0 0 0 1 0 1 1 2   , M 1 =   0 0 1 0 0 1 0 0 0   .
A main point here is that the dependence in u 3 is linear. When it is superlinear, the system does not fit any more into the scope of Theorem 1.1. It is however analyzed in Theorem 1.2.

About more general linear diffusions

Let us now make some comments on the following example where diffusions are more general than in Corollary 3.1.

                 u i ∈ W 1,1 0 (Ω), h i (u 1 , u 2 ) ∈ L 1 (Ω), i = 1, 2, u 1 - N i,j=1 a ij ∂ x i x j u 1 = h 1 (u 1 , u 2 ) + f 1 , u 2 - N i,j=1 b ij ∂ x i x j u 2 = h 2 (u 1 , u 2 ) + f 2 , (27)
where

a ij , b ij ∈ R, N i,j=1 a ij ξ i ξ j , N i,j=1 b ij ξ i ξ j ≥ α|ξ| 2 , α ∈ (0, +∞), ∀ξ = (ξ 1 , ..., ξ N ) ∈ R N . ( 28
)
We consider the operators

       D(A 1 )[resp. D(A 2 )] := v ∈ W 1,1 0 (Ω), N i,j=1 a ij ∂ x i x j v [resp. N i,j=1 b ij ∂ x i x j v] ∈ L 1 (Ω) , A 1 v := N i,j=1 a ij ∂ x i x j v, A 2 v := N i,j=1 b ij ∂ x i x j v.
It is easy to see that the assumptions (A), (A inf ) are satisfied. Thus if (H1), (H2), (M ) are satisfied like in Corollary 3.1, then for the approximate solutions u n of this system, as defined in the proof of Theorem 1.1, it follows that u n , h n (•, u n ) are bounded in L 1 (Ω) 2 and u n lies in a compact set of L 1 (Ω) 2 . However, the extra condition (Φ) is not satisfied in general so that it is not clear whether h n (•, u n ) is uniformly integrable. Actually, we have to choose here a different strategy with does not seem to be generalized to the abstract setting of Section 2. It consists in looking at the equation satisfied by T k (u n 1 + ηu n 2 ) where T k is a cut-off function. It is then easy to pass to the limit in the nonlinear terms of the truncated approximate system since they are multiplied by T k (u n 1 + ηu n 2 ) which vanishes for u n i large for all i. Therefore a.e. convergence is sufficient to pass to the limit. The difficulty is then to prove precise estimates independent of n, in terms of η, in order to control the other terms as it done in the parabolic case (see [START_REF] Pierre | Weak solutions and supersolutions in L 1 for reaction-diffusion systems[END_REF], [START_REF] Pierre | Global Existence in Reaction-Diffusion Systems with Dissipation of Mass : a Survey[END_REF], [START_REF] Laamri | Global existence for reaction-diffusion systems with nonlinear diffusion and control of mass[END_REF]). This approach through cut-off functions T k is precisely developed in Section 4 to prove Theorem 1.2 and we refer the reader to this other approach without giving more details here.

Examples with linear diffusions and Robin-type boundary conditions

Now we analyze what happens for systems like (24) when the boundary conditions are different. If we replace the homogeneous boundary conditions of (24) by homogeneous Neumann boundary conditions, then the result is exactly the same. On the other hand, the situation is quite more complicated if the boundary conditions are of different type for each of the u i 's. This is actually connected with the content of the assumption (Φ). For simplicity, we do it only for a 2 × 2-system. Given a nonlinearity h satisfying (H1), (H2), (M ), we consider the following system with general Robin-type boundary conditions.

           i = 1, 2, u i ∈ W 1,1 (Ω) + , h i (u 1 , u 2 ) ∈ L 1 (Ω), u i -d i ∆u i = h i (u 1 , u 2 ) + f i , λ i u i + (1 -λ i )∂ ν u i = ψ i on ∂Ω, λ i ∈ [0, 1], d i ∈ (0, +∞), ψ i ∈ [0, +∞), f i ∈ L 1 (Ω) + . (29) Corollary 3.3 Let (f 1 , f 2 ) ∈ L 1 (Ω) + × L 1 (Ω) + .

Assume the nonlinearity h satisfies (H1), (H2), (M ).

Assume moreover [0 ≤ λ 1 , λ 2 < 1] or [λ 1 = λ 2 = 1 and ψ 1 = ψ 2 = 0]. Then the system (29) has a solution. Remark 3.4 It is known that the case when one of the λ i is equal to 1 is different (see the analysis in [START_REF] Martin | Influence of mixed boundary conditions in some reaction-diffusion systems[END_REF]). For instance, finite time blow up may occur for the associated evolution problem when the boundary conditions are u 1 = 1, ∂ ν u 2 = 0 (see [START_REF] Bebernes | Finite time blowup for semilinear reactive-diffusive systems[END_REF], [START_REF] Bebernes | Finite-time blowup for a particular parabolic system[END_REF]). Then the operator A 1 does not satisfy (A inf ) as easily seen by considering (as in [START_REF] Martin | Influence of mixed boundary conditions in some reaction-diffusion systems[END_REF]) the following simple example: u σ (x) := cosh(σx)/ cosh(σ) on Ω := (-1, 1) with σ ∈ (0, +∞). Then u σ ≥ 0 and u σ = 1 on ∂Ω. But -Ω u σ = u σ (-1) -u σ (1) → -∞ as σ → +∞. Here we consider only the cases that directly fall into the scope of Theorem 1.1. A few other cases could be treated directly like λ 1 = λ 2 and positive data ψ 1 , ψ 2 or also λ 1 = 0 = λ 2 and ψ 1 = 0.

Proof of Corollary 3.3. Here we define

D(B i ) = {u ∈ W 2,1 (Ω); λ i u i + (1 -λ i )∂ ν u i = ψ i on ∂Ω}, B i u := -d i ∆u i .
Then the closure A i of B i in L 1 (Ω) satisfies the assumption (A) (see [START_REF] Brezis | Semilinear elliptic equations in L 1[END_REF]).

• Let us first assume 0

≤ λ i < 1, i = 1, 2. Then, for u i ∈ D(B i ) Ω B i u i = - Ω d i ∆u i = - Ω d i ∂ ν u i = d i Ω (1 -λ i ) -1 (λ i u i -ψ i ) ≥ -d i (1 -λ i ) -1 Ω ψ i .
This remains valid for A i by closure. Thus the assumption (A inf ) is satisfied. For the condition (Φ), we come back to (26) with the same ϕ, b. Let p q (r) be a standard approximation of sign + (r -k) like

p q (r) = 0, ∀ r ∈ (-∞, k -1/q]; p q (r) = q (r -k) + 1, ∀ r ∈ [k -1/q, k]; p q (r) = 1, ∀ r ∈ [k, +∞). (30)
Then, for u j ∈ D(B j ), j = 1, 2 and for V :

= 2 i,j=1 m ij d j u j Ω p q (ϕ(u) -k)b • M 0 Bu = Ω p q (V )|∇V | 2 - ∂Ω p q (V )∂ ν V ≥ - ∂Ω p q (V )∂ ν V. - ∂Ω p q (V )∂ ν V = ∂Ω p q (V ) 2 i,j=1 m ij d j (1 -λ j ) -1 (λ j u j -ψ j ).
If ψ j = 0 for all j = 1, ..., m, then -∂Ω p q (V )∂ ν V ≥ 0 and letting q → +∞, we obtain that B and therefore A satisfies the condition (Φ). When ψ j = 0 for some j, we only obtain

Ω sign + (ϕ(u) -k)b • M 0 Bu ≥ - ∂Ω sign + (V -k) m i,j=1 m ij d j (1 -λ j ) -1 ψ j =: -η(V, k).
The point is that we can slightly modify the proof of Theorem 1.1 to get the same conclusion with this weaker estimate from below. Indeed, in the present case, we have to add the term η(V n , k) in the inequality [START_REF] Rothe | Global solutions of reaction-diffusion systems[END_REF] where V n = 2 i,j=1 m ij d j u n j . Since η(V n , k) tends to 0 as k → +∞ uniformly in n, the rest of the proof remains unchanged if we choose k ε large enough so that η(V, k ε ) < ε also.

• Assume now that λ 1 = λ 2 = 1. Then we make the same choice as in (26) and we see that

- Ω sign + (ϕ(u) -k)b • M 0 Au ≥ - ∂Ω sign + ( 2 i,j=1 m ij d j ψ j -k)∂ ν ( 2 i,j=1 m ij d j u j ) ≥ 0, since u j ≥ 0, u j = 0 on ∂Ω imply that ∂ ν u j ≤ 0 on ∂Ω.
For the same reason, (A inf ) holds since -

Ω ∆u j = - ∂Ω ∂ ν u j ≥ 0.

Examples with nonlinear diffusions

Let us now consider nonlinear diffusions.

• We start with porous media type equations.

   For i = 1, ..., m, u i ∈ L 1 (Ω) + , ϕ i (u i ) ∈ W 1,1 0 (Ω), h i (u) ∈ L 1 (Ω), u i -∆ϕ i (u i ) = h i (u) + f i , (31)
where

ϕ i : [0, +∞) → [0, +∞) is continuous, increasing with ϕ i (0) = 0, lim s→+∞ s -(N -2) + /N ϕ i (s) = +∞. Corollary 3.5 Let (f 1 , • • • , f m ) ∈ (L 1 (Ω) + ) m .
Assume the nonlinearity h satisfies (H1), (H2), (M ). Then the system (31) has a solution.

Proof. We naturally define

D(A i ) := {u i ∈ L 1 (Ω) ; ϕ i (u i ) ∈ W 1,1 0 (Ω), ∆ϕ i (u i ) ∈ L 1 (Ω) }, A i u i := -∆ϕ i (u i ).
It is classical that the operators A i satisfy the assumptions (A), (A inf ) (see [START_REF] Brezis | Semilinear elliptic equations in L 1[END_REF], [START_REF] Ph | Équations d'évolution dans un espace de Banach et applications[END_REF]). To check that (Φ) is satisfied, we consider r ∈ [0, +∞) m → ϕ(r) := m i,j=1 m ij ϕ j (r j ) and b := (1, ..., 1) t . Then

Ω sign + (ϕ(u) -k)b • M 0 Au = - Ω sign + ( m i,j=1 m ij ϕ j (u j ) -k)∆   m i,j=1 m ij ϕ j (u j )   ≥ 0.
We then apply Theorem 1.1.

• A second example with nonlinear diffusions is the following where p ∈ (1, +∞) :

   For i = 1, 2, u i ∈ W 1,p-1 0 (Ω) + , f i ∈ L 1 (Ω), d i ∈ (0, +∞), u i -d i ∆ p u i = h i (u 1 , u 2 ) + f i , (32)
where for v ∈ W 1,p-1 (Ω), ∆ p v := ∇ • (|∇v| p-2 ∇v).

Corollary 3.6 Let (f 1 , f 2 ) ∈ L 1 (Ω) + ×L 1 (Ω) + .
Assume the nonlinearity h satisfies (H1), (H2), (M ). Then the system (32) has a solution.

Proof. Here we define for i = 1, 2

D(B i ) := {v ∈ W 1,p 0 (Ω) ∩ L 2 (Ω) ; ∆ p v ∈ L 2 (Ω)}, B i (v) := -d i ∆ p v.
And the operators A i are defined as the closure of B i in L 1 (Ω). Then the assumptions (A), (A inf ) are satisfied (see e.g. [START_REF] Herrero | Asymptotic behaviour of the solutions of a strongly nonlinear parabolic problem[END_REF]). Let us prove that (Φ) holds. For p q defined as in (30) and ûi ∈ D(B i ), i = 1, 2, we have

- Ω p q (û 1 + û2 )(∆ p û1 + ∆ p û2 ) = Ω (∇û 1 + ∇û 2 )(|∇û 1 | p-2 ∇û 1 + |∇û 2 | p-2 ∇û 2 )p q (û 1 + û2 ). The mapping r ∈ R N → |r| p is convex. Therefore its gradient r ∈ R N → p|r| p-2 r ∈ R N is monotone, which means (r 1 -r 2 ) • (|r 1 | p-2 r 1 -|r 2 | p-2 r 2 ) ≥ 0 , ∀ r 1 , r 2 ∈ R N .
We apply this with r 1 := ∇û 1 (x), r 2 := -∇û 2 (x), x ∈ Ω to deduce after integration on Ω:

-

Ω p q (û 1 + û2 )(∆ p û1 + ∆ p û2 ) ≥ 0.
And letting q → +∞ gives -Ω sign + (û 1 + û2 -k)(∆ p û1 + ∆ p û2 ) ≥ 0 and by closure

Ω sign + (û 1 + û2 -k)(A 1 û1 + A 2 û2 ) ≥ 0 ∀ u 1 ∈ D(A 1 ), ∀ u 2 ∈ D(A 2 ). ( 33
)
For the condition (M ), we choose

b := (1, 1) t , ϕ(r 1 , r 2 ) := c 1 r 1 + c 2 r 2 , c 1 := {(m 11 + m 21 )d 1 } 1/(p-1) , c 2 := {(m 21 + m 22 )d 2 } 1/(p-1)
where

M 0 := (m ij ) 1≤i,j≤2 . Then for u 1 ∈ D(A 1 ), u 2 ∈ D(A 2 ), u = (u 1 , u 2 ), Ω sign + (ϕ(u) -k)b • M 0 Au = Ω sign + (c 1 u 1 + c 2 u 2 -k) 2 i,j=1 m ij d j A j u j = Ω sign + (c 1 u 1 + c 2 u 2 -k)[A 1 (c 1 u 1 ) + A 2 (c 2 u 2 )]
≥ 0, the last inequality coming from (33) applied with û1 := c 1 u 1 , û2 := c 2 u 2 . This ends the proof of Corollary 3.6.

• To end this section, let us comment on the following system which is a model of situations where the operators A i are very different from each other :

           p = 2 u 1 ∈ W 1,p-1 0 (Ω) + , u 2 ∈ W 1,1 0 (Ω) + , h 1 , h 2 , f 1 , f 2 ∈ L 1 (Ω), u 1 -∆ p u 1 = h 1 (u 1 , u 2 ) + f 1 , u 2 -∆u 2 = h 2 (u 1 , u 2 ) + f 2 . (34)
As a consequence, the compatibility condition (Φ) is not generally satisfied. However, if the nonlinearity h satisfies (H1), (H2), (M ), then for the approximate solution u n defined in Section 2,

h n (•, u n ) is bounded in L 1 (Ω) × L 1 (Ω) and therefore u n lies in a compact set of L 1 (Ω) × L 1 (Ω).
Thus we may assume that, up to a subsequence, u n converges to some u in L 1 (Ω) × L 1 (Ω) and a.e. so that h n (•, u n ) converges a.e. to h(•, u). Unfortunately we do not know whether we can pass to the limit in the approximate version of system (34). And it is not clear how the use of cut-off function T k as in the next section can help. We leave this as an open problem.

Proof of Theorem 1.2

For each k = 1, ..., m, we define σ k := min{α k , β k }, γ k := |α k -β k | so that (see [START_REF] Laamri | Global existence of classical solutions for a class of reaction-diffusion systems[END_REF] for the definition of I, J ):

h k (u) = (β k -α k ) Π m =1 u σ B(u), B(u) := k 1 Π i∈I u γ i i -k 2 Π j∈J u γ j j . (35)
We first solve the approximate system with the bounded data f n i := inf{f i , n}, n ∈ N.

Lemma 4.1 There exists a nonnegative solution

u n ∈ ∩ p∈[1,∞) W 2,p (Ω) +m of    For i = 1, ..., m, u n i -d i ∆u n i = (β i -α i )Π m k=1 (u n k ) σ k B(u n ) + f n i in Ω, ∂ ν u i = 0 on ∂Ω. (36)
Moreover, we have

γ j u n i + γ i u n j -∆(γ j d i u n i + γ i d j u n j ) = γ j f n i + γ i f n j , ∀ i ∈ I, j ∈ J, (37)
and u n is bounded in L 1+η (Ω) m for some η > 0.

Proof. By the abstract Lemma 2.1, for ε ∈ (0, 1), there exists a regular nonnegative solution u ε of

       For i = 1, ..., m, u ε i -d i ∆u ε i = (β i -α i ) Π m k=1 (u ε k ) σ k B(u ε ) 1 + ε Π k=1 m (u ε k ) α k + Π m k=1 (u ε k ) β k + f n i in Ω, ∂ ν u ε i = 0 on ∂Ω. (38)
Indeed the nonlinearity is here quasi-positive and uniformly bounded by max i γ i /ε. Then by multiplying the equations in i ∈ I, j ∈ J respectively by γ j , γ i , we have

(γ j u ε i + γ i u ε j ) -∆(γ j d i u ε i + γ i d j u ε j ) = γ j f n i + γ i f n j . (39)
If

d := max 1≤k≤m d k , this implies (since u ε i , u ε j ≥ 0) d -1 [γ j d i u ε i + γ i d j u ε j ] -∆(γ j d i u ε i + γ i d j u ε j ) ≤ γ j f n i + γ i f n j , ∂ ν (γ j d i u ε i + γ i d j u ε j ) = 0 on ∂Ω.
We deduce that

γ j d i u ε i + γ i d j u ε j ≤ d -1 I -∆ -1 (γ j f n i + γ i f n j ). (40)
Since the f n i are bounded by n (fixed), using the nonnegativity of u ε , we obtain that the u ε i are bounded in L ∞ (Ω) independently of ε. By standard elliptic regularity results applied to the equation (38) in u i , they are also bounded in ∩ p∈[1,∞) W 2,p (Ω). They are therefore included in a compact set of L ∞ (Ω) as → 0. We can then easily pass to the limit as ε → 0 and obtain the convergence of a subsequence of u ε toward a solution u n of (36).

Then the identities (37) follow by passing to the limit in (39). But (40) remains also valid at the limit and implies that u n i , i = 1, ..., m are bounded in L 1 (Ω) independently of n. Next we may rewrite (37), using

d := min 1≤k≤m d k : d -1 [γ j d i u n i +γ i d j u n j ]-∆(γ j d i u n i +γ i d j u n j ) = γ j [(d -1 d i -1)u n i +f n i ]+γ i [(d -1 d j -1)u n j +f n j ], ∀ i ∈ I, j ∈ J, (41)
together with ∂ ν (γ j d i u n i + γ i d j u n j ) = 0 on ∂Ω. Since the right-hand side is bounded in L 1 (Ω), this implies that γ j d i u n i + γ i d j u n j is bounded in L p (Ω) for all p ∈ [1, N/(N -2) + ) (see e.g. [START_REF] Brezis | Semilinear elliptic equations in L 1[END_REF]). This ends the proof of Lemma 4.1.

In order to rewrite the relations (41), we denote

g n k := (d -1 d k -1)u n k + f n k ≥ 0, ∀ k = 1, ..., m. (42)
When f k log f k is assumed to be in L 1 (Ω) for all k = 1, ..., m, then

sup n∈N max 1≤k≤m Ω g n k + | g n k log g n k | < +∞. (43)
This is due to the "L LOG L" assumption on the f i 's and on the L 1+η -bound on the u n i stated in Lemma 4.1. Next we introduce the solutions G n k of

d -1 G n k -∆G n k = g n k in Ω, ∂ ν G n k = 0 on ∂Ω, G n k ≥ 0, k = 1, ..., m, (44)
so that the relations (37), (41) may be rewritten

γ j d i u n i + γ i d j u n j = γ j G n i + γ i G n j , ∀ i ∈ I, j ∈ J. ( 45
)
Our goal is to prove that the nonlinearity of the system (36) is bounded in L 1 (Ω) m independently of n. It will imply enough compactness on u n to pass to the limit. The following lemma will provide the key estimate.

Lemma 4.2 Under the assumptions of Theorem 1.2, there exists θ n ∈ W 2,1 (Ω) +m (unique) such that

   γ j d i θ n i + γ i d j θ n j = γ j G n i + γ i G n j , ∀ i ∈ I, j ∈ J, k 1 Π i∈I (θ n i ) γ i = k 2 Π j∈J (θ n j ) γ j (or B(θ n ) = 0), ∂ ν θ n k = 0 on ∂Ω, ∀ k = 1, ..., m. (46)
Moreover

sup n max 1≤k≤m ∆θ n k L 1 (Ω) < +∞. ( 47
)
We postpone the proof of this lemma and we end the proof of Theorem 1.2.

Without loss of generality, we assume that 1 ∈ I. Using the equation in u n 1 given in the approximate system (36), and the fact that B(θ n ) = 0, we have

u n 1 -θ n 1 -d 1 ∆(u n 1 -θ n 1 ) + γ 1 Π m k=1 (u n k ) σ k [B(u n ) -B(θ n )] = ρ n , ρ n := f n 1 -θ n 1 + d 1 ∆θ n 1 . (48)
Now the main point is that B(u n ), B(θ n ) have the property that

B(u n ) = b n (•, u n 1 ), B(θ n ) = b n (•, θ n 1 ), (49)
where r → b n (•, r) is a increasing function. More precisely, we deduce from the relations (45), (46) that

γ i d j θ n j = γ j G n i + γ i G n j -γ j d 1 θ n 1 , ∀ j ∈ J, γ 1 d i θ n i = γ i d 1 θ n 1 + γ 1 G n i -γ i G n 1 , ∀i ∈ I. γ i d j u n j = γ j G n i + γ i G n j -γ j d 1 u n 1 , ∀ j ∈ J, γ 1 d i u n i = γ i d 1 u n 1 + γ 1 G n i -γ i G n 1 , ∀i ∈ I. Plugging this into B(u n ), B(θ n , recalling that B(v) = k 1 Π i∈I v γ i i -k 2 Π j∈J v γ j j , we indeed obtain (49) by setting b n (x, r) := k 1 Π i∈I (δ i d 1 r + F n i (x)) γ i -k 2 Π j∈J (F n j (x) -δ j d 1 r) γ j , ∀ r ∈ [r n -, r n + ], (50)
where we define

   δ k := γ k /(γ 1 d k ), ∀ k = 1, ..., m, F n i := d -1 i G n i -δ i G n 1 , ∀ i ∈ I, F n j := δ j G n 1 + d -1 j G n j , ∀ j ∈ J, r n -(x) := max i∈I [F n i (x)] -/(δ i d 1 ), r n + (x) := min j∈J F n j (x)/(δ j d 1 ). ( 51
)
We multiply the equation (48) by sign(u n 1 -θ n 1 ) and we integrate. Then we use the two following properties

Ω sign(u n 1 -θ n 1 )[u n 1 -θ n 1 -d 1 ∆(u n 1 -θ n 1 )] ≥ 0, sign(u n 1 -θ n 1 )[B(u n ) -B(θ n )] = sign(u n 1 -θ n 1 )[b n (•, u n 1 ) -b n (•, θ n 1 )] = |b n (•, u n 1 ) -b n (•, θ n 1 )|, to obtain γ 1 Ω Π m k=1 (u n k ) σ k |b n (•, u n 1 ) -b n (•, θ n 1 )| ≤ Ω |ρ n |.
Since ρ n is bounded in L 1 (Ω) by Lemma 4.2, and since b n (•, θ n 1 ) ≡ 0, this implies that the nonlinearity of the equation in

u n 1 u n 1 -d 1 ∆u n 1 + γ 1 Π m k=1 (u n k ) σ k b n (•, u n 1 ) = f n i in Ω, is bounded in L 1 (Ω) independently of n. This implies that ∆u n 1 is bounded in L 1 (Ω).
Going back to (45), it follows that ∆u n k is bounded in L 1 (Ω) m for all k = 1, ..., m as well. Up to a subsequence, we may deduce that u n k converges for all k to some u k in L 1 (Ω) and a.e. and that ∇u n k converges in L 1 (Ω) N to ∇u k . Note also that by Fatou's lemma, Π k (u k ) σ k B(u) ∈ L 1 (Ω). Now to end the passing to the limit, we look at the equation satisfied by T R (v n i ) where v n i := u n i + ε j =i u n j and ε ∈ (0, 1) and the T R are C 2 -cut-off functions satisfying

T R (s) = s if s ∈ [0, k -1], T R (s) = 0 if s ≥ k, T R (s) ≤ k, ∀s ∈ [0, +∞), 0 ≤ T R (s) ≤ 1, T R (s) ≤ 0 for all s ≥ 0. ( 52
) -∆T R (v n i ) = -T R (v n i )∆v n i -T R (v n i )|∇v n i | 2 ≥ -T R (v n i )∆(u n i + ε j =i u n j ).
Coming back to the system (36), we have that for all ψ ∈ C ∞ Ω + ,

Ω ∇ψ∇T R (v n i ) ≥ Ω ψT R (v n i )    Π k (u n k ) σ k B(u n )[ β i -α i d i + ε j =i β j -α j d j ] + f n i -u n i d i + ε j =i f n j -u n j d j    .
We know that, up to a subsequence, u n converges in L 1 (Ω) m and a.e. to some u ∈ W 1,1 (Ω) m . We may pass to the limit along this subsequence in the above inequality. Indeed the nonlinear terms on the right converge a.e. and T R (

v n i )Π k (u n k ) σ k B(u n ) is uniformly bounded, while the f n k , u n k , k = 1, ..., m converge in L 1 (Ω). Moreover, ∇T R (v n i ) = T R (v n i )∇v n i converges in L 1 (Ω) to ∇T R (v i ), v i := u i + ε j =i u j .
At the limit, we obtain the same inequality without the superscript n . Then, we let ε go to 0 to obtain

Ω ∇ψ∇T R (u i ) ≥ Ω ψT R (u i ) Π k (u k ) σ k B(u) β i -α i d i + f i -u i d i .
And now we let R → +∞ to obtain

Ω ∇ψ∇u i ≥ Ω ψ Π m k=1 (u k ) σ k B(u) β i -α i d i + f i -u i d i , or equivalently ∀i = 1, ..., m, Ω ψu i + d i ∇ψ∇u i ≥ Ω ψ {Π k (u k ) σ k B(u)(β i -α i ) + f i } . ( 53
)
Throughout the rest of the proof, we will assume that γ j G n i + γ i G n j ≡ 0. By maximum principle applied to the equation (44) defining the G n k , a n ij := inf(γ j G n i + γ i G n j ) > 0. Also for the rest of the proof, choosing n 0 large enough, we fix c ∈ (0, +∞) such that

0 < c < a n ij /[γ j d i + γ i d j ], ∀ (i, j) ∈ I × J, ∀n ≥ n 0 . (55)
This definition of c will be used only in STEPS 6 and 7 of the proof.

For simplicity, we now drop the superscript ' n ' in the rest of the proof.

STEP 2: Existence of θ n satisfying (46). Again we assume (without loss of generality) that 1 ∈ I so that we may use (50), (51). Thus, for all x ∈ Ω, the function r

∈ [r -(x), r + (x)] → b(x, r) is increasing where r -, r + are defined in (51). Moreover b(x, r -(x)) = -k 2 Π j∈J (F j (x) -d 1 δ j r -(x)) γ j ≤ 0, b(x, r + (x)) = k 1 Π i∈I (d 1 δ i r + (x) + F i (x)) γ i ≥ 0.
Thus there exists a unique θ

1 (x) ∈ [r -(x), r + (x)] such that b(x, θ 1 (x)) = 0. Since the function (x, r) → b(x, r)
is regular, by the implicit function theorem, so is x → θ 1 (x). The function θ k are then uniquely determined from the first line of (46) which we rewrite as:

θ i := d 1 δ i θ 1 + F i , ∀ i ∈ I; θ j := F j -d 1 δ j θ 1 , ∀j ∈ J. (56)
It remains to prove that ∂ ν θ k = 0 on ∂Ω for all k = 1, ..., m. This will be a consequence of the following computation (see (60) ).

k 1 Π i∈I θ γ i i = k 2 Π j∈J θ γ j j , which implies log k 1 + i∈I γ i log θ i = log k 2 + j∈J γ j log θ j . (57)
Differentiating this leads to

i∈I γ i ∇θ i θ i = j∈J γ j ∇θ j θ j . (58)
Inserting (56) in this formula, we obtain ∇θ 1 in terms of the θ k , namely

d 1 ∇θ 1 A = j∈J γ j ∇F j θ j - i∈I γ i ∇F i θ i , A := m k=1 γ k δ k θ k . ( 59
)
Since ∇F k • ν = 0 on ∂Ω for all k = 1, ..., m, it follows from this identity that ∇θ 1 • ν = 0 on ∂Ω as well.

And by (56) it also follows that

∇θ k • ν = 0 on ∂Ω, ∀ k = 1, ..., m. (60)
Differentiating once more (58) gives

i∈I γ i ∆θ i θ i - γ i |∇θ i | 2 θ 2 i = j∈J γ j ∆θ j θ j - γ j |∇θ j | 2 θ 2 j ,
or also, using again (56) and the definition of A in (59),

d 1 A∆θ 1 = i∈I - γ i ∆F i θ i + γ i |∇θ i | 2 θ 2 i + j∈J γ j ∆F j θ j - γ j |∇θ j | 2 θ 2 j . ( 61
)
Our goal is to estimate the L 1 -norm of ∆θ 1 . We remark that, if we denote

α k := γ k δ k /Aθ k , k = 1, ..., m, then 0 ≤ α k ≤ 1, m k=1 α k = 1.
But the relation (61) may be rewritten

d 1 ∆θ 1 = i∈I α i [- ∆F i δ i + |∇θ i | 2 δ i θ i ] + j∈J α j [ ∆F j δ j - |∇θ j | 2 δ j θ j ]. (62)
According to the definition of F k in (51) and to the definition of G k (or more precisely of G n k ) in ( 44) and (43), we know that ∆F k L 1 (Ω) , k = 1, ..., m is bounded in terms of the data (independently of n).

Thus i∈I -α i ∆F i δ i + j∈J α j ∆F j δ j L 1 (Ω) ≤ C, (63)
where C is independent of n.

We also have that Ω ∆θ 1 = ∂Ω ∂ ν θ 1 = 0. Inserting this into (62) and (63) gives

Ω j∈J α j |∇θ j | 2 δ j θ j ≤ Ω i∈I α i |∇θ i | 2 δ i θ i + C, (64)
where again C does not depend on n. Therefore, it is sufficient to bound the right-hand side of (64) to obtain a bound on ∆θ 1 L 1 (Ω) and this will end the proof of Lemma 4.2 (since an L 1 -bound on ∆θ 1 implies an L 1 -bound on ∆θ k for all k = 1, ..., m).

STEP 4 : A bound from below on the θ k . The previous step indicates that one has to bound |∇θ i | 2 /θ i in L 1 (Ω) for all i ∈ I. The identity (59) says that c 2

i 0 Ω |∇F i 0 | 2 a + |F i 0 | ≤ Ω | f i 0 | log | f i 0 | + | f i 0 |[a -1 d -1] + 1,
where f i 0 = d -1 F i 0 -∆F i 0 = g i 0 /d 1 -δ i 0 g 1 as indicated in (70). We conclude from the estimate (43) that

∇F i 0 θ i 0 √ θ 1 A L 2 (Ω)
≤ C where C is independent of n.

(73) STEP 6: More bounds from below on the θ k . We prove here the two following facts, where the real number c was defined in (55): Ω = Ω I ∪ Ω J , Ω I := {x ∈ Ω ; θ i (x) > c, ∀i ∈ I}, Ω J := {x ∈ Ω ; θ j (x) > c, ∀j ∈ J}, sup i∈I θ i ≥ c 1 > 0, for some c 1 ∈ (0, ∞) independent of n.

(74)

For the first part of (74), assume by contradiction that there exists x ∈ Ω \ (Ω I ∪ Ω J ). Then there exists (i, j) ∈ I × J such that, for this x: θ i (x) ≤ c and θ j (x) ≤ c, or equivalently, γ j d i θ i (x) ≤ γ j d i c, γ i d j θ j (x) ≤ γ i d j c.

Let us add these last two inequalities. Using the first line of (46) and the definition of c in (55), we have

γ j G i (x) + γ i G j (x) = γ j d i θ i (x) + γ i d j θ j (x) ≤ [γ j d i + γ i d j ]c < inf{γ j G i + γ i G j }.
And this is a contradiction. Whence the first statement of (74).

For the second one, let us first note that sup Now let x ∈ Ω J so that, by the previous statement, θ j (x) ≥ c for all j ∈ J. We then use the second line of (46) to obtain k 2 c j∈J γ j ≤ k 2 Π j∈J θ j (x)

γ j = k 1 Π i∈I θ i (x) γ i ≤ k 1 [sup i∈I θ i ] i∈I γ i .
This implies that sup i∈I θ i ≥ c 1 := min{c, k 2 c j∈J γ j /k 1

[ i∈I γ i ] -1 }.
Whence the second statement of (74).

STEP 7: End of the proof of Lemma 4.2. This is where we use that I (for instance) has at most two elements. Indeed, let us go back to the expression of ∇θ 1 / √ θ 1 in (65). We already know by STEP 5 that all terms indexed by j ∈ J are bounded in L 2 (Ω). Since F 1 ≡ 0, there is at most one term indexed by i ∈ I, namely none if I = {1}, and only ∇F i 0 /(θ i 0 √ θ 1 A) if I = {1, i 0 }. If I = {1}, it immediately follows that ∇θ 1 / √ θ 1 is bounded in L 2 (Ω) independently of n. If I = {1, i 0 }, then sup i∈I θ i = sup{θ 1 , θ i 0 }. It follows from the second line of (74) in STEP 6 that (72) holds. Consequently, as proved in Remark 4.4, ∇F i 0 /(θ i 0 √ θ 1 A) is bounded in L 2 (Ω) independently of n. Having controlled all terms in (65), we can conclude that ∇θ 1 / √ θ 1 is itself also bounded in L 2 (Ω) independently of n. By symmetry, this also holds for ∇θ i 0 / θ i 0 . This implies that the right-hand side of (64) in STEP 4 is bounded independently of n and ends the proof of Lemma 4.2.

Let us now state the following technical lemma which was used in two places in the previous proof. If moreover f ≥ 0, then

Ω |∇F | 2 F ≤ Ω f log + f + (d -1 -1)f + de -1 .
Proof. For the first inequality of the lemma, let us first remark that, if we set F := a F , then by homogeneity, If f ≥ 0, then F ≥ 0 by maximum principle. We multiply the equation in F by log(F + ε) and we integrate by parts. Then,

Ω dF log(F + ε) + |∇F | 2 F + = Ω f log(F + ε). (76)
We apply the Young's inequality (75) with r := f (x), s := log(F (x) + ε). Then f log(F + ε) ≤ (f log f -f ) + F + ε, so that, using (76), we deduce

[F ≥1]
dF log(F + ε)

+ Ω |∇F | 2 F + ε ≤ Ω f log f -f + F + ε -d [F ≤1]
F log(F + ε).

STEP 3 :

 3 Differentiating B(θ) = 0. The condition B(θ) = 0 means

  x∈Ω I θ i (x) ≥ inf x∈Ω I θ i (x) ≥ c, ∀i ∈ I.

Lemma 4 . 5

 45 Let F ∈ W 2,1 (Ω) such that, for some d ∈ (0, +∞)dF -∆F = f, ∂ ν F = 0 on ∂Ω with f, f log |f | ∈ L 1 (Ω).For a ∈ (0, +∞), we haveΩ |∇F | 2 a + |F | ≤ Ω |f | log |f | + |f |[(ad) -1 -1] + 1.

Ω |∇F | 2 aΩ|∇ F | 2 1+

 22 + |F | = a | F | .Let us now multiply the equation d F -∆ F = f /a by sign( F ) log(| F | + 1). We obtainΩ | F | log(| F | + 1) + |∇ F | 2 1 + | F | = 1 a Ω f sign( F ) log(| F | + 1) ≤ 1 a Ω |f | log(| F | + 1).Now we use the Young's convexity inequality∀ r ∈ [0, +∞), ∀ s ∈ R, rs ≤ (r log r -r) + e s .(75)We apply it with r := |f (x)|, s := log(| F (x)| + 1) to deduce Ω |f | log(| F | + 1) ≤ Ω |f |[log |f | -1] + | F | + 1. From the equation in F , we also derive d Ω | F | ≤ Ω |f |/a. The first inequality of Lemma 4.5 follows.

  Now we let ε → 0 and we use[F ≥1] F log F ≥ 0, f log f ≤ f log + f, Ω dF = Ω f, x log x ≥ -e -1 , ∀ x ∈ (0, 1),to deduce the second estimate of Lemma 4.5.

  1.2 holds even if |I| and |J| > 2. But it is not clear how to extend our main Lemma 4.2 to the general case. We leave this as an open problem. Theorem 1.2 provides a solution u such that h i (u) ∈ L 1 (Ω). It is easy to write down explicit examples where the nonlinear terms Π k u α k k , Π k u β k k are not separately in L 1 (Ω). For instance, let N > 4 and let us introduce the following function σ where r

	Remark 1.5

But, on the other hand, by passing directly to the limit as n → +∞ in (37), we have that for i ∈ I, j ∈ J Ω ψ(γ j u i + γ i u j ) + ∇ψ∇[γ j d i u i + γ i d j u j ] = Ω ψ(γ j f i + γ i f j ).

This implies that all inequalities in (53) are actually equalities so that, for all i = 1, ..., m,

This ends the proof of Theorem 1.2.

Proof of Lemma 4.2.

STEP 0 : The case m = 2. For simplicity, we drop the index 'n'. We assume I = {1}, J = {2}. Then the system (46) is equivalent to

It is easily seen that the equation in θ 1 has a unique regular solution (see STEP 2 for more details) and differentiating this equation gives

If γ = 1, we immediately have

By (42) and (44), we know that ∆G is bounded in L 1 (Ω) independently of n as soon as the f k are (only) in L 1 (Ω) + for all k. Therefore so is ∆θ 1 . If γ = 1, then integrating (54) and using ∂ ν θ 1 = 0, we obtain

Again, the last integral is bounded independently of n if the f k are (only) assumed to be in L 1 (Ω). Therefore so is the first integral. But by positivity, this implies that θ γ-2

(Ω) independently of n. So is ∆θ 1 by going back to (54). And finally the same holds for ∆θ 2 = ∆(G -δθ 1 ).

We now come back to the general situation m ≥ 2 and with the LLogL assumptions on the f k .

STEP 1 : Let us first treat the trivial case when there exists (i 0 , j 0 ) ∈ I × J such that γ j 0 G n i 0 + γ i 0 G n j 0 ≡ 0 (i.e. G n i 0 ≡ 0 ≡ G n j 0 ). Then, by the first line of (46), θ i0 ≡ 0 ≡ θ j 0 . Using again the first line of (46), we deduce that θ n k = G n k for all k = 1, ..., n. Thus, the conclusion of the lemma is obvious in this case.