Florence Alberge
email: florence.alberge@u-psud.fr

Deep Learning Constellation Design for the AWGN Channel with Additive Radar Interference

Keywords: Machine learning, optimization methods, digital modulation, radar interference

Radar and wireless communication coexistence is considered in this paper as a possible solution to face the exploding demand and rising congestion in wireless networks. The transmission medium is modeled as an AWGN channel with additive radar interference. Standard constellations are not optimal in this context and an auto-encoder (AE) is used to design proper constellations and corresponding receiver devices. AE is a powerful tool in neural networks that shares strong similarities with communication systems. This technique is particularly relevant in the lack of an analytical expression of the loss function. In the asymptotic region (high interference regime), the optimal constellation shape is known and the AE always converges towards this optimal solution. In the other regions, the AE is able to yield solutions that outperform the standard configurations. Several demapping alternatives are also considered leading to the conclusion that it is possible to maintain the communication link in the presence of radar interference independently of the interference power. This is a step further compared to previous works in which solutions were limited to low or high interference regimes.

I. INTRODUCTION

Increasing data traffic demand in cellular and wireless networks leads to a rethinking of the spectrum allocation policy. The actual spectrum regulation with bands allocated to specific users or services is suboptimal and a significant part of the spectrum is seldom used. Indeed radars and radio navigation infrastructure occupy an important portion of the available spectrum with low usage efficiency. Several regulators across the world have started considerations for release portions of governmental radar bands to be shared with commercial wireless services. Consequently, the concept of spectrum sharing in radar has gained lots of interest [START_REF] Locke | An Assessment of the Near-Term Viability of Accommodating Wireless Broadband Systems in the 16751710 MHz, 17551780 MHz, 35003650 MHz, and 42004220 MHz, 43804400 MHz Bands[END_REF], [START_REF]Report on Collective Use of Spectrum (CUS) and other spectrum sharing approaches[END_REF], [START_REF] Khan | On opportunistic spectrum access in radar bands: Lessons learned from measurement of weather radar signals[END_REF]. Literature on radar and wireless communication coexistence is vast. It is possible however to design the radar and/or the communication system to alleviate its impact on the other device by using different techniques such as waveform design [START_REF] Aubry | Radar waveform design in a spectrally crowded environment via nonconvex quadratic optimization[END_REF], [START_REF] Huang | Radar waveform design in spectrum sharing environment: Coexistence and cognition[END_REF], interference mitigation [START_REF] Deng | Interference Mitigation Processing for Spectrum-Sharing Between Radar and Wireless Communications Systems[END_REF] or spatial separation [START_REF] Khawar | Spectrum sharing between S-band radar and LTE cellular system: A spatial approach[END_REF]. In [START_REF] Turlapaty | A joint design of transmit waveforms for radar and communications systems in coexistence[END_REF], a joint design of the radar transmission waveform's power spectrum and the power spectral density of a multi-carrier communication system is proposed. As a complementary approach, it is also interesting to consider an unaltered radar signal and a communication system and in particular it is useful to have an idea of the attainable performance when both systems coexist. The performance of standard communication systems in the presence of radar interference has been studied through numerical experiments in several publications. For example, WiMax and ground based radar are considered in [START_REF] Lackpour | Overview of interference mitigation techniques between WiMAX networks and ground based radar[END_REF] and an OFDM system subject to ultra-wide band radar interference is considered in [START_REF] Pasya | Analysis of interference from UWB radar signals on a digital OFDM transmission system[END_REF]. A theoretical performance analysis is given in [START_REF] Chiriyath | Inner Bounds on Performance of Radar and Communications Co-Existence[END_REF] where bounds on performance of the joint system are measured in terms of data information rate for communications and radar estimation information rate for radar. The capacity of the communication channels with both additive white Gaussian noise and radar interference is investigated in [START_REF] Shahi | On the Capacity of the AWGN Channel With Additive Radar Interference[END_REF]. It is proved that communication is possible even in the presence of high radar interference. This latter case results in a loss of half the degrees of freedom compared to a channel without interference. The capacity achieving input distribution under average power constraint is characterized in [START_REF] Shahi | On the Capacity of the AWGN Channel With Additive Radar Interference[END_REF]. The optimal input distribution has independent modulus and phase. The phase is uniformly distributed in [0, 2π] whereas the modulus is discrete with countably infinite many mass points but only finitely many in any bounded interval. Finite inputs are used in practice whereas the optimal input is continuous (Gaussian distributed). The objective of this paper is to design finite constellations for the additive radar interference channel. To the best of our knowledge, this problem has not been widely investigated in the past. Very recently, this question was addressed in [START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF]. The objective function considered for the minimization process is the symbol error rate which is based on the analytical expressions given in [START_REF] Nartasilpa | On the Error Rate of a Communication System Suffering from Additive Radar Interference[END_REF] for the weak and for the very strong interference regime leading to very interesting results for these two situations. The method in [START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF] cannot be extended to other regimes since closed form expressions are not known in the general case. Instead, it is proposed here to use deep learning and AEs for constellation design. With deep learning, the network is able to learn directly from the observations and closed-form expressions are no longer required.

Deep learning is a branch of machine learning research which has been introduced to allow machines and computer systems to improve with experience and data accumulation. An auto-encoder (AE) is a type of artificial neural network used to learn efficient representations of the data in an unsupervised manner. The analogy between AEs and communication systems has been used in recent publications [START_REF] Shea | An Introduction to Deep Learning for the Physical Layer[END_REF]. More generally, deep learning techniques have already been used for solving communication problems. The vast majority is dedicated to the design of the receiver. For example, multi-user detection in code-division multiple access (CDMA) has been investigated in [START_REF] Mitra | Neural network techniques for adaptive multiuser demodulation[END_REF], [START_REF] Aazhang | Neural networks for multiuser detection in code-division multiple-access communications[END_REF]. The ability of a neural network to decode polar codes with MAP performance is analyzed in [START_REF] Gruber | On deep learningbased channel decoding[END_REF]. Neural network techniques demonstrate convincing results for short block length codes but suffer from high decoding complexity. Indeed, the training complexity of deep learning-based channel decoders scales exponentially with the number of information bits and cannot be reasonably considered for standard codes used in practice. In [START_REF] O'shea | Learning to communicate: Channel auto-encoders, domain specific regularizers, and attention[END_REF], the emitter and decoder are considered jointly. Here again the complexity limits the usefulness of the method since the proposed method can only be applied to short block length.

Taking into consideration those works, we believe that AEs are well-suited for constellation design. Indeed, the dimensionality of this problem is modest and the encoder of the AE can be trained to optimize the constellation whereas the decoder is trained to recover the binary message from the received data. Our objective is not end-to-end learning of communications systems through deep neural network like in [START_REF] O'shea | Learning to communicate: Channel auto-encoders, domain specific regularizers, and attention[END_REF]. Here, the neural network is considered only for modulation and demodulation tasks. The proposed solution is compatible with other elements of the communication chain such as error-correcting codes. The strengths of the proposed methodology are listed below. The objective function to be optimized is a bound of the mutual information [START_REF] Vincent | Stacked denoising auto-encoders : learning useful representation in a deep network with a local denoising criterion[END_REF] which is an usual criterion in communication and in particular in constellation design. Neither an analytical expression of the mutual information (or of the symbol error rate) nor an a priori on the geometry of the constellation are needed for the optimization process. This is in contrast with previous works [START_REF] Yang | Design of Circular Signal Constellations in the Presence of Phase Noise[END_REF], [START_REF] Krishnan | Constellation Optimization in the Presence of Strong Phase Noise[END_REF], [START_REF] Kwak | Spiral QAM: A Novel Modulation Scheme Robust in the Presence of Phase Noise[END_REF], [START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF]. The constellation geometry and the labeling are optimized jointly. The proposed mapping/demapping solution can be used in combination with powerful channel decoders such as turbo-codes or Low Density Parity-Check codes (LDPC) and is fully compatible with wireless communications schemes already deployed.

The paper is organized as follows. The system model is depicted in Section II. Past work related to this model is recalled and our contribution detailed. The strong interference regime is considered in section III. It is proved that, in this specific regime, an M-PAM like constellation should be used. AEs are discussed in Section IV-A and the proposed structure is presented. Section V is devoted to constellation design. Numerical results are given in Section VI where the demapping operation is discussed in order to exhibit the best association for each interference regime.

II. SYSTEM MODEL

The system under consideration is composed of a narrowband communication channel suffering from the interference of a high power, short duty cycle radar pulse. An analytical model for the discrete-time complex-valued signal at the communication receiver is given in [START_REF] Nartasilpa | Let's share CommRad: Effect of radar interference on an uncoded data communication system[END_REF]. The channel model is an AWGN with additive interference, with known constant amplitude

√

I and unknown random phase θ I uniformly distributed in [0, 2π]. A justification for this model can be found in [START_REF] Nartasilpa | Communications System Performance and Design in the Presence of Radar Interference[END_REF] and references therein. From the communication perspective, a single carrier system is considered with narrow-band frequency-flat slowly-varying fading model with known channel gain at the receiver. Extension to time-varying channels is obtained by averaging the error rates over the channel statistics. The amplitude of the radar interference can be accurately estimated since the radar system has slowvarying parameters [START_REF] Nartasilpa | Communications System Performance and Design in the Presence of Radar Interference[END_REF]. At the opposite, the phase suffers rapid variations and is difficult to track. Let's denote the discrete-time signal at the communication receiver as Y then

Y = X + √ Ie jθ I + Z (1)
where X is the channel input chosen from the constellation X = {X i = √ S i e jθi , i = 0, ..., M -1} with cardinality M and Z is a complex Gaussian noise with zero-mean and unit-variance. The average power of the constellation is S meaning that E[|X| 2 ] = S. The random variables (X, θ I , Z) are mutually independent and the channel is memoryless. From these notations, the average Signal to Noise Ratio (SN R) is S and the average Interference to Noise Ratio (IN R) is I. The objective of this paper is to optimize the constellation alphabet X for reliable transmission through the AWGN channel with additive interference given in (1). Both S and I are fixed and known at the receiver. Three different regimes will be considered here. For each one, the state of the art is first recalled and emphasis is given on our contribution.

1) I S. In the high interference regime, it is optimal at the receiver to estimate the phase of the radar interference as the phase of the received signal and to use this estimation to subtract its contribution from the received signal [START_REF] Nartasilpa | Let's share CommRad: Effect of radar interference on an uncoded data communication system[END_REF]. In this interference cancellation process, a portion of the useful signal is also subtracted and the capacity is only half the capacity of the interferencefree system. More specifically, it is proved in [START_REF] Shahi | On the Capacity of the AWGN Channel With Additive Radar Interference[END_REF] that

lim I→∞ C(S, I) = 1 2 log(1 + S) (2)
where C(S, I) denotes the capacity of the channel in (1). The Gaussian input is optimal in the high interference regime. However, in practice, the channel input is constrained to a finite size alphabet. 2) I S. The Gaussian input performs also very well (but is not optimal) in the weak interference regime with achievable rates close to the upper bound in [START_REF] Ihara | On the capacity of channels with additive non-Gaussian noise[END_REF]. The capacity achieving distribution under average power constraint has independent modulus and phase [START_REF] Shahi | On the Capacity of the AWGN Channel With Additive Radar Interference[END_REF]. The phase is uniformly distributed in [0, 2π] and the modulus is discrete with countably infinitely many mass points, but only finitely many in any bounded interval. It can also be seen in [START_REF] Nartasilpa | On the Error Rate of a Communication System Suffering from Additive Radar Interference[END_REF] that, in this regime, the optimal MAP decoder is the classical minimal distance decoder applied to [START_REF] Locke | An Assessment of the Near-Term Viability of Accommodating Wireless Broadband Systems in the 16751710 MHz, 17551780 MHz, 35003650 MHz, and 42004220 MHz, 43804400 MHz Bands[END_REF] in which the interference is considered as Gaussian noise.

3) The intermediate regime is the most challenging case.

The interference is neither small enough to be neglected nor large enough to be accurately estimated. The expression of the channel capacity does not admit a simple expression like in the other two cases. Our contributions are given below. The main contribution is the design of finite size constellations with enhanced perfor-mance compared to standard constellations and the design of the corresponding soft demapping rules. The constellation as well as the soft demapping are optimized with an AE. Thanks to this optimization framework, the constellation shape, the constellation labeling 1 and the soft demapping device are optimized jointly. This in contrast with [START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF] in which the constellation labeling is not addressed. Indeed, in [START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF] the symbol error rate (SER) is considered in the optimization problem whereas here the loss function is based on bit probabilities. This is a key difference that will be discussed in section V-C. Note also that [START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF] focuses on the low and high IN R regime whereas here solutions are proposed for the three regimes. The resulting constellations are evaluated in terms of achievable rates and/or error probability and compared to theoretical bounds and also to standard constellations. The second contribution is to show that the designed schemes are efficient both with uncoded and coded data. Indeed the AE outputs soft information which can be easily connected to a turbo or LDPC decoder. To the best of our knowledge, previous works considered uncoded data only. However, the gain from the constellation design may be absorbed by the channel coding. It can be seen here that the proposed schemes outperform the standard ones even when coded data are considered. Finally, a new demapping rule (called DMHighSNR in the following) well adapted to the additive radar interference channel is introduced here. This new rule outperforms the existing ones in the three regimes.

An analysis of the interference channel in the high IN R regime is provided in the next section. It is proved that the constellation points should be collinear to avoid equivocal areas. This is consistent with the expression of the capacity of the real AWGN channel with power constraint given in (2).

III. INTERFERENCE CANCELLATION IN THE

Y IC = Y - √ Ie j θ I = Y 1 - √ I |Y | (3)
where |Y | stands for the modulus of the complex variable Y . In the whole paper, X stands for the channel input, Y is the channel output in (1) and Y IC is the channel output after interference estimation and cancellation given in (3).

Result 1: Without loss of generality, X = X 0 = √ S 0 e jθ0 where X 0 is a given point of the constellation X . If SN R → ∞ and I S then

Y IC = √ S 0 2 e jθ0 + e j(2θ I -θ0) + S 0 o(1) (4)
where

√ S 0 o(1)
is the remainder term of the Taylor expansion and o(1) → 0 when S0 I → 0.

1 A constellation labeling is the assignment of a bit pattern to each symbol in a signal-set constellation. Proof: The absolute square of Y reads

|Y | 2 = I 1 + S 0 I e j(θ0-θ I ) + e j(θ I -θ0) + S 0 I ( 5
)
Since S is the average energy of the constellation and since the constellation is bounded and of finite size, I S leads to I S 0 Thus,

√ I |Y | admits the first-order asymptotic expansion √ I |Y | = 1 - 1 2 S 0 I e j(θ0-θ I ) + e j(θ I -θ0) + S 0 I o(1) (6)
Y IC √ I = 1 2
S 0 I e jθ0 + e j(2θ I -θ0) + S 0 I o(1) (7)
which concludes the proof.

From result 1, the locus of Y IC for a given transmitted symbol

X = X 0 = √ S 0 e jθ0 is a circle with center √ S0
2 e jθ0 and radius

√ S0 2 .
Let's denote the point with Cartesian coordinates (0; 0) as O. The circles corresponding to two different constellation points overlap at least in O since each circle passes through the origin. The loci for a 16QAM constellation are depicted in Figs. 1 and2 for I = S 2 and SN R = 20 dB and SN R = 40 dB respectively. We can observe several overlapping points which may complicate the de-mapping task. We prove in the next result that there exists a constellation geometry with one and only one overlapping point, namely O.

Result 2: Let's assume that C 1 and C 2 are two circles defined as

C 1 = z 1 = √ S 1 2 e jθ1 + e jθ , θ ∈ [0; 2π) C 2 = z 2 = √ S 2 2 e jθ2 + e jθ , θ ∈ [0; 2π)
Let's define X 1 and X 2 as two distinct points of X with respective polar coordinates Proof: For compactness, (x 1 ; y 1 ) denote the Cartesian coordinates of X 1 and (x 2 ; y 2 ) denote the Cartesian coordinates of X 2 . With these notations and for i ∈ {1, 2}, C i is a circle with center coordinates (x i , y i ) and radius

R i = x 2 i + y 2 i . If M with Cartesian coordinate (x; y) is a point of C 1 ∩ C 2 , then x and y are such that x 2 -2x 1 x + y 2 -2y 1 y = 0 (8) x 2 -2x 2 x + y 2 -2y 2 y = 0 (9)
We distinguish here the two following cases: 1)

x 1 = x 2 .
Then, by subtraction of ( 8) and ( 9), either y = 0 or y 1 = y 2 . Since it is assumed that X 1 = X 2 , the case {x 1 = x 2 and y 1 = y 2 } does not require to be considered. If y = 0 then x = 0 or x = 2x 1 hence

|C 1 ∩ C 2 | = 2. 2) x 1 = x 2 .
By subtracting ( 8) and ( 9), we have

x = y 2 -y 1 x 1 -x 2 y (10)
which can be plugged into (8) leading to

y 1 -y 2 x 2 -x 1 2 + 1 y 2 = 2 x 2 y 1 -x 1 y 2 x 2 -x 1 y (11)
As a consequence,

|C 1 ∩ C 2 | = 1 if and only if x 2 y 1 = x 1 y 2 . This condition is equivalent to sin(θ 1 -θ 2 ) = 0 (12)
Hence X 1 , X 2 and O are collinear.

An illustration of this result is given in Fig. 3 where the constellation is composed of 16 collinear points and is symmetric with respect to O. As expected, the overlapping is limited to O.

In the following, a constellation with cardinality M and such that all points are on a straight line passing through the origin is termed an M-RPAM by analogy with a rotated M-PAM. Note however that regular spacing between the constellation points is not required. It is proved in the next result that perfect decoding is possible even when Y IC = 0. S) and the noiseless case. Let's assume that X is an M-RPAM constellation and that {θ X ; θ X + π} denote the possible phases of the constellation points. The channel input is X = X 0 = √ S 0 e jθ0 and the corresponding channel output is

Y = √ S 0 e jθ0 + √ Ie jθ I . If Y IC = 0 then X 0 is recovered from the observation Y as X 0 = 2Re(Y e -jθ X )e jθ X (13)
Proof:

From (3), Y IC = 0 is equivalent to Y = 0 or |Y | = √ I. In the high IN R regime, I = S 0 thus Y = √ S 0 e jθ0 + √ Ie jθ I = 0. The other alternative |Y | = √ I is equivalent to |Y | 2 = I. Since |Y | 2 = S 0 + I + 2 √ S 0 Icos(θ 0 -θ I ), |Y | = √ I is equivalent to cos(θ 0 -θ I ) = -1 2 S0 I . The expression of sin(θ 0 -θ I ) is inferred from cos(θ 0 -θ I ) 2 + sin(θ 0 -θ I ) 2 = 1 leading to sin(θ 0 -θ I ) = 1 -S0 4I with ∈ {-1; 1}.
Replacing these two equations into the expression of the noiseless channel output leads to

Y = S 0 e jθ0 1 2 -j I S 0 1 - S 0 4I (14)
Since θ 0 -θ X = 0 mod π, we have e j(θ0-θ X ) ∈ {+1, -1} and we can conclude that

Re(Y e -jθ X ) = 1 2 S 0 e j(θ0-θ X ) (15)
The channel input is thus recovered from the observation Y by computing 2Re(Y e -jθ X )e jθ X . Hence, at high SN R, a M-RPAM constellation yields perfect decoding. Other geometric configurations may share this property. However, since collinear points decrease the number of equivocal area, it is likely to observe, with noisy data, an improvement in terms of BER with M-RPAM compared to other constellations. This point will be confirmed in the simulation part.

IV. AUTO-ENCODERS AND

where w and b are called weights or trainable parameters and where g(.) is either a non-linear function or the identity mapping [START_REF] Goodfellow | Deep Learning[END_REF]. The weights contain the information learned by the network from exposure to training data. The weights of the whole layers are optimized jointly, each layer being updated to follow both the representational needs of the layer above and the needs of the layer below [START_REF] Chollet | Deep Learning with Python[END_REF]. The optimizer implements the back-propagation algorithm. The loss function takes the prediction of the output of the network and the true target and computes a distance score that is used as a feedback signal to adjust the weights. In this representation, in and out are vectors with respective length d in and d out , w is a matrix with size d out × d in and b is a vector with length d out . In the following, we use the shorthand notation

out = f din→dout,g P ( in ) (17)
where P = {w, b} is the set of trainable parameters.

B. Auto-encoders and proposed structure

AEs are a specific type of feedforward neural networks where the input is the same as the output. The Denoising AE is of particular interest here. The idea behind Denoising AEs is simple. In order to force the neural network to discover more robust features and prevent it from simply learning the identity, we train the AE to reconstruct the input from a corrupted version of it. The network may be viewed as consisting of two parts termed as encoder and decoder. Usually, the corruption is introduced at the input. To mimic a communication system, the perturbation is introduced here at the output of the encoder and before entering the decoder [START_REF] Shea | An Introduction to Deep Learning for the Physical Layer[END_REF]. From the network perspective, the perturbation is modeled by a layer (Channel layer) that does not involve any trainable parameters. A schematic representation of the proposed structure is given in Fig. [START_REF] Aubry | Radar waveform design in a spectrally crowded environment via nonconvex quadratic optimization[END_REF].

Let's denote the length of the binary message C to be transmitted as k. C is the input of the AE. The input space is of size M = 2 k . X is the complex symbol at the output of the encoder and is thus a vector of length 2 (real and imaginary part of the symbol). X is a representation of C that will be optimized to be robust to the perturbation introduced by the communication channel. Unfortunately, there are no formula to determine the right number of layers and the right size for each layer. The general workflow to find the appropriate network's size is to start with relative few layers and parameters and increase the size of the layers or the number of layers until reaching a minimum or a threshold in the validation loss [START_REF] Chollet | Deep Learning with Python[END_REF]. With this methodology, the encoder (modulator) is defined as where t and l stand respectively for the hyperbolic tangent and linear activation functions. This structure involves a moderate number of layers with small size but appears to be sufficient for performing the encoding task. The activation functions are chosen symmetric with respect to 0 since X should not be restricted to the quadrant of the complex plane with positive real part and positive imaginary part. Before entering the channel, the block of data is normalized such that the constellation has an average energy equal to S. The decoder (demodulator) is defined as

X = f k→16,t P1 • f 16→8,t P2 • f 8→4,t P3 • f 4→2,l P4 (C) (18)
Ĉ = f 2→128,r P 1 • f 128→64,r P 2 • f 64→32,r P 3 • f 32→M,s P 4 (Y ) (19)
where r and s stand respectively for the rectified linear unit (relu) and sigmoid activation functions. The rectified linear unit is recognized as one of the most efficient activation function. The sigmoid activation function should be used in the output layer for a binary classification problem. Since Ĉ ∈ {0, 1} k , the problem considered here can be seen as a binary classification problem. This decoder structure has also been used in [START_REF] Gruber | On deep learningbased channel decoding[END_REF] for decoding polar codes.

C. Training the network

For a given training session, the triplet (k, S, α) is fixed. For all the simulations in this paper, the number of layers and the size of each layer is also fixed to the values given in [START_REF] Gruber | On deep learningbased channel decoding[END_REF][START_REF] O'shea | Learning to communicate: Channel auto-encoders, domain specific regularizers, and attention[END_REF] and in Fig. A third alternative will be discussed later on. Like in [START_REF] Gruber | On deep learningbased channel decoding[END_REF], the gradient of the loss function is calculated over the entire input set during each epoch. The batch size is chosen equal to 2 k , the number of epochs is 2 19 which is a good compromise between performance and complexity. The stochastic optimization is performed by a stochastic gradient with fixed learning rate l r = 0.01. In our architecture, the channel layer generates a new realization of the noise Z and of the phase θ I each time it is used. The network is then trained in order to optimize the reconstruction error

L(C, Ĉ) = -log p(C| Ĉ)
The problem under consideration here is a binary classification problem. In that case, the reconstruction error is the crossentropy loss [START_REF] Vincent | Stacked denoising auto-encoders : learning useful representation in a deep network with a local denoising criterion[END_REF], namely

L(C, Ĉ) = - j C(j) log Ĉ(j)+(1-C(j)) log (1 -Ĉ(j)) (20)
and Ĉ(j) is such that Ĉ(j) = P (C(j) = 1| Ĉ) where C(j) stands for bit j of C and Ĉ(j) stands for bit j of Ĉ. Since Ĉ(j) = P (C(j) = 1| Ĉ), connection with soft decoders such as LDPC or turbo-codes is straightforward. The training of the network is performed by solving the following optimization problem arg min

P,P E C,Z,θ I [L(C, Ĉ)] (21)
When the distribution of the input is unknown, the expectation is replaced by the empirical average over the training samples [START_REF] Goodfellow | Deep Learning[END_REF]. It is proven in [START_REF] Vincent | Stacked denoising auto-encoders : learning useful representation in a deep network with a local denoising criterion[END_REF] that the optimization problem in ( 21) is strongly connected to the maximization of the mutual information I(X; C) between X and C which makes sense in the present context. The optimization problem is solved with a unit-norm constraint on the columns of w to avoid scaling effects. The goal of training is to find the optimal set of parameters {P 1 , .

D. Using a trained network

Once trained, the optimized network can be introduced into the communication chain. Let's define the binary representation of i for 0

≤ i ≤ 2 k as C i . Then, X i = f k→16,t P1 • f 16→8,t P2 • f 8→4,t P3 • f 4→2,l P4 (C i ).
The association C i ↔ X i can be also stored in a table. This association between a constellation point and a particular input is usually called labeling. At the receiver and when the network is the one in Fig. 4-a, the demapping task can be handled as Ĉi = f 2→128,r

P 1 • f 128→64,r P 2 • f 64→32,r P 3 • f 32→M,s P 4 (Y (X i ))
where Y (X i ) is the channel output corresponding to the input X = X i . The demapping is called DMAE(Y) in that case. When the network used during the training is the one in Fig. 4-b, the demapping task can be handled as Ĉi = f 2→128,r

P 1 • f 128→64,r P 2 • f 64→32,r P 3 • f 32→M,s P 4 (Y IC (X i )) where Y IC (X i )
is the channel output after interference cancellation. This demapping will be referred as DMAE(Y IC ). With this framework, the constellation geometry, the labeling and the demapping device are optimized jointly. This in contrast with previous publications. The proposed structure is deployed, in the next section, in various channel conditions. The AE is conceived here as an off-line optimization tool with the objective of finding the appropriate constellation and demapping device for each scenario. We are now ready to use the AE defined above. Our goal here is to evaluate the capability of the AE to find the most appropriate constellation for both the low and high IN R regime in which communication should be possible without experiencing a severe degradation [START_REF] Shahi | On the Capacity of the AWGN Channel With Additive Radar Interference[END_REF]. The capacity of the AE to cope with the severe channel impairments inherent to the intermediate regime will also be tested.

A. Constellation design at high IN R

As a first step, the performance of several standard constellations are compared in terms of average symbol error rate (SER) for SN Rs ranging from 10 to 30 dB. The results are displayed in Fig. 5 where we used the analytical expression of the SER given in [START_REF] Nartasilpa | On the Error Rate of a Communication System Suffering from Additive Radar Interference[END_REF] for an arbitrary constellation at high IN R and with an interference cancellation decoder. The numerical results in Fig. 5 illustrate perfectly the theoretical analysis of section III. For any constellation size, there exists a threshold on the SN R such that beyond this threshold the M-PAM constellation outperforms the other ones. The ability of the AE to find this optimal constellation is now under consideration. For the high IN R regime, the receiver should estimate and subtract the interference. Hence, the AE is trained with Y IC as channel output instead of Y . Training is performed for different values of the triplet (M, SN R, α) where α is such that I = S α . The values of α are chosen between 1.8 and 2. [START_REF] Huang | Radar waveform design in spectrum sharing environment: Coexistence and cognition[END_REF]. The corresponding SER can be found in Fig. [START_REF] Huang | Radar waveform design in spectrum sharing environment: Coexistence and cognition[END_REF]. We can observe that for each triplet the AE converges towards a constellation which has either the same SER as the M-PAM (best standard constellation) or slightly better. For example, the constellation obtained with (M = 16, SN R = 20 dB, α = 2) is shown in Fig. 8. This is a 16-RPAM3 with non-equally spaced symbols. The SER with this constellation is 0.3499 whereas the SER of the 16-PAM at 20 dB is 0.3547. As a conclusion, the AE has been able to find the appropriate constellation for each triplet. The efficiency of the corresponding decoder will be evaluated in Section VI with LDPC-coded data.

B. Constellation design at low to intermediate IN R

The AE is trained here with Y as a channel output. Several pairs (M, α) are considered with SN R = 20 dB. The achievable rates for the constellations designed by the AE are given in Fig. 6. The results are compared with the upper bound from Ihara's work in [START_REF] Ihara | On the capacity of channels with additive non-Gaussian noise[END_REF] which coincides, for α ∈ [0 1], with the rate achieved with a Gaussian input. The capacity of the Gaussian channel log 2 (1 + S 1+I ) where the interference is replaced by a Gaussian noise with same power is also plotted and serves also to benchmark the achievable performance. Finally, our results are compared to the achievable rates of an input distribution with uniform phase and one mass-point at

√

S for the modulus. This constellation corresponds to the class of optimal input distributions for the AWGN channel with additive radar interference ( [START_REF] Shahi | On the Capacity of the AWGN Channel With Additive Radar Interference[END_REF], theorem 1). The one mass-point constellation is obviously easy to design since the modulus is necessarily equal to √ S. The optimization remains challenging in the general case for which the optimal number of mass-points is unknown. This is here that the AE is expected to bring new solutions. We can observe in Fig. 6, that when α is small (α < 0.4), the standard constellation cannot be improved with the AE. Indeed, the constellations obtained in this area are in general rotated standard QAM 4 . In contrast, when 0.4 ≤ α ≤ 0.8 the constellations obtained with the AE outperform significantly the standard constellations in particular when M = 8 or M = 16. Significant improvements can also be obtained when M = 32 but at higher SN R. When the channel impairment is too important, the AE converges automatically towards a degenerated constellation. For example, at SN R = 20 dB and at α = 0.7, even if M is chosen equal to 16, a constellation with cardinality 8 will be obtained. In that case, the optimization process should be run again with M = 8 instead of M = 16 to obtain a constellation with M distinct points. This is a nice property of the AE making the choice of M much easier. The constellations obtained at (SN R = 20 dB, M = 8, α = 0.6) and (SN R = 20 dB, M = 16, α = 0.5) are given in Fig. 7 and are discussed below.

C. Analysis of constellation shapes

The constellations obtained with the AE and for various triplets (SN R, α, M ) are given in Figs. 7 and8. It is noticeable that the constellations in Fig. 7 can be seen as optimal input distributions ( [START_REF] Shahi | On the Capacity of the AWGN Channel With Additive Radar Interference[END_REF], theorem 1) after uniform phase sampling. Except for the point at the origin, the constellation with cardinality equal to 8 is a sampled version of the "one mass point at SN R". To take the analysis one step further, the constellation corresponding to (SN R = 20 dB, M = 16, α = 0.5) is again plotted in Fig. [START_REF] Lackpour | Overview of interference mitigation techniques between WiMAX networks and ground based radar[END_REF].

We can see that this constellation has four different values of the modulus and thus 4 orbits with 4 points on each ((

) and 16 evenly distributed values for the phase. However on a specific orbit, the phase is not evenly distributed. From the labels given in the scatter plot, decision regions for each bit can be inferred. Let's denote the binary decomposition of C as (C(3), C(2), C(1), C(0)). For clarity, the constellation is again plotted in Fig. 9 in which red markers indicate a constellation point such that C(i) = 0 whereas blue stars indicate a constellation point such that C(i) = 1. We can observe a clear dividing line between the two categories (C(i) = 0 / C(i) = 1) whatever the value of i is. The demapping task can thus be handled with binary classification and this is exactly the function of the AE decoder. Constellation design for the radar interference channel is also addressed in [START_REF] Nartasilpa | Communications System Performance and Design in the Presence of Radar Interference[END_REF]- [START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF]. In those publications, the loss function is the symbol error rate (SER) with power constraint on the constellation. The labelling is not addressed. The constellation given in Fig. 9 (right) is obtained at (SN R = 20 dB, M = 16, α = 0.5) with the AE and with the SER as loss function. It can be compared with the constellation in Fig. 9 (left) obtained with our method and also at (SN R = 20 dB, M = 16, α = 0.5). When the optimization is based on the SER, the constellation tends to shape like a concentric hexagon as already mentioned in [START_REF] Nartasilpa | Communications System Performance and Design in the Presence of Radar Interference[END_REF]- [START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF]. We can observe in this particular example that the cardinality of the constellation on the right is 13 instead of 16 (several constellation points coincide). The irregular spacing between the constellation points is optimized by taking into account the binary representation of the input which is not possible when the optimization is based on the SER. From this example, we can conclude that a joint optimization (labelling + constellation shape) is mandatory in the intermediate regime.

We then turn our attention to Fig. [START_REF] Turlapaty | A joint design of transmit waveforms for radar and communications systems in coexistence[END_REF]. We can observe that when the IN R increases, degrees of freedom are being lost on the -1 0 1

In-phase constellation geometry. Indeed, the phase of the constellation points is allowed to span [0; 2π] in the constellations of Fig. 7 whereas the phase is limited to smaller intervals when the IN R increases. In the extreme case (high IN R regime) the phase belongs to a set with cardinality 2. Loosing degrees of freedom is equivalent to increasing the prior knowledge on the transmitted symbol. This prior knowledge is a key point for recovering the transmitted symbol from the observation for large IN R values (see the proof of Result 3). The constellation shape at high IN R is consistent with the results in section III and with [START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF]. The AE is therefore able to find the most suitable constellation for given channel parameters. At this point, the efficiency of the demapping has not been evaluated yet. In the next section, the joint performance of the mapping/demapping optimized with the AE are compared to standard settings in the context of LDPC-encoded data. Left: loss function in [START_REF] Vincent | Stacked denoising auto-encoders : learning useful representation in a deep network with a local denoising criterion[END_REF]. Right: loss function = SER [START_REF] Nartasilpa | Communications System Performance and Design in the Presence of Radar Interference[END_REF].

-1 0 1

In-phase In-phase

-1 0 1 Quadrature C(2) C(2)=0 C(2)=1 -1 0 1
In-phase

VI. NUMERICAL RESULTS

Let's denote the information message as B and the corresponding LDPC-encoded message as C. The data-block C is split into binary words of length log 2 (M ) which are mapped to X ∈ X and sent through the channel in (1). At the receiver side a soft demapper is used. The output of the demapper is converted into LLR and serves as an input to the LDPC decoder which is based on the sum-product algorithm. The parity-check matrix of the LDPC code proceeds from the DVB-S.2 standard. The block length of the code is fixed and equal to 64800, the code rate R is chosen according to the achievable rates of the constellation. The constellation is either a standard one (QAM, PAM, PSK) or one of the constellations in Figs. 78. The demapping is characterized by the conditional probability p Y |X (y, x) from which the LLR of bit j is obtained as

LLR(C j ) = log x∈X :xj =0 p Y |X (y, x) x∈X :xj =1 p Y |X (y, x) (22)
and serves as an input to the LDPC decoder. Depending on the values of the channel parameters, several expressions can be obtained for p Y |X (y, x):

• At low IN R (I S), the channel is expected to behave as an AWGN channel [START_REF] Nartasilpa | On the Error Rate of a Communication System Suffering from Additive Radar Interference[END_REF]. In that case, the channel conditional distribution reads

p Y |X (y, x) ∝ e -|y-x| 2 1+I (23)
We shall refer to this demapping as DMLowINR.

• At high IN R (I S), it is possible to obtain an accurate estimate θ I = arg(Y ) leading to Y IC after interference cancellation. It can be proven with the same line of arguments as in the proof of result 1 (see also [START_REF] Nartasilpa | On the Error Rate of a Communication System Suffering from Additive Radar Interference[END_REF] for an alternative proof) that

Re(y IC e -j θ I ) = √ s 0 cos(θ 0 -θ I ) + V (24)
where V ∼ N (0, 1 2 ) and where the polar form of x is √ s 0 e jθ0 . Therefore, the conditional probability to be used

at high IN R reads 1 √ π e -(Re(y IC e -j θ I )-√ s0cos(θ0-θ I )) 2 (25)
We shall refer to the soft demapping based on (25) as DMHighINR.

• When I ≈ S, arg(Y ) is not anymore an accurate estimation of θ I . Instead and since we want to compute a conditional probability, an estimation of θ I conditioned on the knowledge of the transmitted symbol X is obtained as θ I (X) = arg(Y -X). The accuracy of the estimation depends on the noise level. In that case

p Y |X (y, x) ∝ e -|y-x- √ Ie j θ I (x) | 2 (26)
The corresponding soft demapping shall be referred as DMHighSNR. {Y -X i } 0≤i≤M -1 as an input to the AE decoder where X i belongs to a given constellation X . In that case, the training process is divided into two steps: 1) constellation design with one of the schemes in Fig. 4, 2) update of the decoder weights with one of the schemes in Fig. 4 for which the weights P of the encoder are frozen to output the constellation X given by step one.

) I S Re(Y IC ), Im(Y IC ) DMAE(Y -{X i }) S 1 Re(Y -X i ), Im(Y -X i ) X i ∈ X , i = 0, 1, ...M -1 Table I LIST OF DEMAPPING DEVICES.
The full list of soft-demapping devices under consideration is given in table I. The numerical results are depicted in Figs. 11 to 14 in terms of bit error rate (BER) for different IN R and SN R values. We can draw some conclusions. Whatever the value of α, the best configuration always involves X AE . The AE is therefore able to yield appropriate constellations in the low and high IN R regime but also in the intermediate one. This is a substantial advantage compared to the state of the art. Even in the high IN R regime where we already know that a PAM-like constellation should be used, the AE is able to bring an improvement thanks to an optimized spacing of the constellation points. The highest gains are obtained in the intermediate regime and it is thus mandatory to design a proper communication chain in this region. These gains could be very important. For example, at α = 0.5, a gap of more than 2 dB is observed between the best association (X AE +DMHighSNR) and the standard one (16QAM+DLowINR). Not surprisingly DMLowINR and DMHighSNR lead to reasonable performance in their target area but are dismissed outside. For all configurations under consideration, DMHighSNR outperforms the other alternatives even when SN R ≈ 14 dB but may fail at lower SN R. We turn to the discussion on the performance of DMAE. In Fig. 11, DMAE(Y) exhibits close performance to DMLowSNR. Indeed, when α = 0.2, the assumption √ Ie jθ I + Z ∼ N (0, 1 + I) is true. Then eq. ( 23) is the correct metric and we cannot expect the training process to perform better. The gap of 0.1 dB between DMLowSNR and DMAE(Y) should then be considered as a good result for DMAE(Y). When α = 0.5, the radar interference can no longer be considered as a Gaussian variable. As a consequence, DMLowINR exhibits poor results as can be seen in Fig. [START_REF] Shahi | On the Capacity of the AWGN Channel With Additive Radar Interference[END_REF]. DMAE(Y) does not rely on a Gaussian assumption (on the radar interference) and is still a good solution at α = 0.5. However, when the IN R level increases, DMAE(Y) is not the best option as can been seen in Fig. [START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF]. At high IN R, we can observe that DMAE(Y IC ) and DMHighINR are within a gap of 0.2 dB meaning that the AE has been able to "learn" the metric. Finally, DMAE((Y-{X i }) is considered in Figs. [START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF]-14 as an alternative to DMHighSNR. We can see that, in both cases, the gap between DMHighSNR and DMAE((Y-{X i }) is less than 0.5 dB which means that DMAE((Y-{X i }) is an appropriate alternative to DMHighSNR in particular when the noise has unknown or intractable closed-form expression.

VII. CONCLUSIONS

The coexistence of a short duty cycle, wide-band radar signaling and of a narrow-band communication signal was discussed in this paper from the communication system perspective. Thanks to deep learning techniques, it was proven that it is possible to design an appropriate constellation and to maintain an acceptable communication link whatever the interference power is. In this paper, AEs were envisioned as an off-line optimization tool. The self adaptation of a given node to the transmission medium is the next step forward.

  HIGH IN R REGIME In the high IN R regime (I S), Y ≈ √ Ie jθ I since the other terms in (1) become negligible when I goes to infinity. The phase of the radar interference can thus be estimated as θI = arg(Y ) and √ Ie j θI = √ I Y |Y | is an estimation of the radar interference. The accuracy of this estimation depends on the SN R value. Let's define Y IC as

Figure 1 .

 1 Figure 1. Locus of Y IC , 16QAM constellation, I = S 2 , SN R = 20 dB.

where S0 I o( 1 )

 1 is the remainder term of the Taylor expansion and o(1) → 0 when S0 I → 0. From Y = √ I e jθ I + S0 I e jθ0 and from (3), we obtain the first-order expansion of Y IC √ I as

Figure 2 .

 2 Figure 2. Locus of Y IC , 16QAM constellation, I = S 2 , SN R = 40 dB.

Figure 3 .

 3 Figure 3. Locus of Y IC -|X | = 16 -I = S 2 -SN R = 40 dB

  COMMUNICATION SYSTEMS A. Deep learning Deep learning is a specific subfield of machine learning. The central problem in machine learning and deep learning is to learn useful representations of the input data from exposure to known examples of inputs and outputs. Deep learning involves several successive layers of representation. A deep neural network is obtained by stacking layers one on top of another. The neural layers considered in this article transform an input data in into an output out as follows out = g(w in + b)

Figure 4 .

 4 Figure 4. Schematic representation of the AE -Black pictograms indicate layers with trainable parameters.

4

 4

  . The training set contains all the possible values of C (binary words with length k). In an AE during the training phase, the targets values are equal to the inputs ie Ĉi = C i where C i is a realization of C. The output of the channel layer could be either Y (see Fig. 4-a, low/mid IN R) or Y IC (see Fig. 4-b, high IN R).

Figure 5 .

 5 Figure 5. SER vs SN R at high IN R (I = S α , α ≈ 2) and with interference cancellation.

Figure 6 .

 6 Figure 6. Achievable rates vs α at SN R = 20 dB and I = S α .

Figure 7 .

 7 Figure 7. Constellation designed with AE. Low to intermediate IN R regime.

Figure 8 .

 8 Figure 8. Constellation designed with AE. Intermediate to high IN R regime.

Figure 9 .

 9 Figure 9. Constellation designed at (SN R = 20 dB, α = 0.5, M = 16). Left: loss function in (20). Right: loss function = SER [25].

Figure 10 .

 10 Figure 10. Constellation designed with AE at (SN R = 20 dB, α = 0.5, M = 16). Decision region for the individual bits.

Figure 11 .

 11 Figure 11. BER vs SN R at α = 0.2 and with LDPC rate R = 3/4. X AE optimized with (M = 16, α = 0.2, SN R = 13 dB).

Figure 12 .

 12 Figure 12. BER vs SN R at α = 0.5 and LDPC rate R = 5/6. X AE optimized with (M = 16, α = 0.5, SN R = 20 dB).

Figure 13 .

 13 Figure[START_REF] Nartasilpa | Signal constellation design in the presence of radar interference and Gaussian noise[END_REF]. BER vs SN R at α = 0.9 and LDPC rate R = 8/10 (M = 8) or R = 3/5 (M = 16). X AE optimized with (M = 8, α = 0.9, SN R = 21 dB).

Figure 14 .

 14 Figure 14. BER vs SN R at α = 2 and LDPC rate R = 3/5. X AE optimized with (M = 16, α = 2, SN R = 20 dB).

•

  Once trained, the AE provides a constellation alphabet associated with a valid soft demapping. The output Ĉ of the demapping block is given in[START_REF] O'shea | Learning to communicate: Channel auto-encoders, domain specific regularizers, and attention[END_REF] where the input should be Y at low and intermediate IN R and Y IC at high IN R. The LLR of bit j is computed as As before, three cases are considered. In the high IN R regime, the input of the AE decoder should be Y IC . The corresponding soft-demapping is denoted DMAE(Y IC ). Outside of the high IN R regime, the input of the AE decoder is Y and the resulting demapping is called DMAE(Y ). We can observe that DMAE(Y ) is based on the same observations (received samples) as DMLowINR whereas DMAE(Y IC ) is based on the same observations as DMHighINR. The last one, DMHigh-SNR, is a bit more intricate since the observation required is Y -X. It is possible to mimic DMHighSNR by giving

	LLR(C j ) = log	1-Ĉj Ĉj	. This soft-demapping is termed
	DMAE.

The constellations obtained with the AE at the high IN R regime are all RPAM. The AE always converges towards the optimal constellation shape.

It should be noticed however that α = 0.4 is the threshold obtained when SN R = 20 dB. We will see in the next section, at α = 0.2 for example, that improvements are possible at lower SN R (compared to standard constellations).


Achim Wick
email: a.wick@itv.rwth-aachen.de

Tan-Trung Nguyen

Frédérique Laurent

Rodney O Fox

Heinz Pitsch

Modeling Soot Oxidation with the Extended Quadrature Method of Moments

Keywords: Soot oxidation, method of moments, EQMOM, statistical soot model

Modeling the oxidation of soot particles in flames is a challenging topic both from a chemical point of view and regarding the statistical treatment of the evolution of the soot Number Density Function (NDF). The Method of Moments is widely-used for the statistical modeling of aerosol dynamics in various applications, and a number of di↵erent moment methods have been established and successfully applied to the modeling of soot formation and growth. However, a shortcoming of existing moment methods is the lack of an accurate, numerically robust, and computationally e cient way to treat soot oxidation, especially regarding the prediction of the particle number density. In this work, the recently developed Extended Quadrature Method of Moments (EQMOM) is integrated with a physico-chemical soot model and combined with a treatment for particle removal by oxidation. This leads to

Introduction

In the development of engineering devices, such as next-generation aero or diesel engines, one of the main concerns is the reduction of particulate emissions. While reactive Computational Fluid Dynamics are already widely-used in the design process, a substantial barrier for an extended use of numerical simulations with respect to emissions reduction is a lack in the predictive quality of models for soot formation, and especially for soot oxidation. Aero engines designed following the rich-burn/quick-quench/lean-burn concept are an example for applications in which soot oxidation is very important. While a lot of soot is formed in the primary rich combustion zone, more than 95 % of the particles are oxidized before reaching the exit of the combustor. There-fore, in order to accurately predict the amount of soot emissions, reliable models for soot oxidation are crucial.

In addition to the gas-phase chemistry, soot models consist of two main building blocks: First, a physico-chemical model, which describes the various micro-processes of soot inception, growth, and oxidation on a single-particle level. Second, as not every single particle can be tracked in a simulation, a statistical model is required. Statistical methods applicable to soot modeling include Monte Carlo methods, sectional methods, and moment methods. The goal of the present paper is the integration of a recently developed advanced quadrature-based moment method [1] with soot models for an improved statistical description of soot evolution, especially of soot oxidation.

Targets for the development of statistical soot models are: high accuracy, computational e ciency, numerical robustness (related to moment realizability), applicability to 3D turbulent simulations (i.e. relevant for engineering applications), easy implementation, and applicability to multivariate models, which parameterize the particles with more than one quantity. While Monte Carlo (MC) methods [2,3] are very accurate and well suitable for multivariate models, they are computationally expensive, and their applicability is restricted to simplified configurations. However, due to their negligible error, they can serve as a reference solution, which other methods can be validated against.

The most widely-used class of methods is the Method of Moments including the Method of Moments with Interpolative Closure (MOMIC) [4], the (Direct) Quadrature Method of Moments (QMOM/DQMOM) [5,6], and the Hybrid Method of Moments (HMOM) [7]. The quadrature-based moment methods have been shown to be very accurate for the description of aerosol dynamics including particle growth and coagulation. While DQMOM is computationally e cient and can easily be extended to multivariate models, it can become numerically unstable. QMOM, which is numerically robust, has recently been extended to the Conditional QMOM [8], which combines accuracy, robustness, and e ciency, and has also been applied in 2D simulations [9]. However, all of the mentioned methods have di culties in dealing with soot oxidation in a mathematically rigorous manner. Often, relatively crude assumptions must be made due to shortcomings in the statistical methods [10]. These shortcomings are discussed in more detail in the next section.

Recently, the Extended Quadrature Method of Moments (EQMOM) has been introduced [1,11] as a statistical approach to solve the Population Balance Equation (PBE). This method retains the accuracy and e ciency of the QMOM algorithm, and has already been applied to soot formation and growth processes [12]. However, although, for reasons discussed below, EQMOM is especially suitable to accurately describe soot oxidation and the related disappearance of particles, a formulation for soot oxidation has not yet been investigated. Soot oxidation is also challenging from a chemical point of view, and it is an active field of research, both experimentally [13] and numerically [14]. It should be noted that the goal of the present paper is not the improvement of chemical soot models. However, in order to benefit from advances of the chemical description of soot oxidation in simulations of engineering interest, high-fidelity statistical methods applicable to soot oxidation must be developed, and the purpose of this paper is to contribute to this development.

The following section provides a more detailed analysis why the statistical description of soot oxidation is challenging within the framework of moment methods. In section 3, two variants of EQMOM are discussed. Section 4 briefly describes the physico-chemical soot model. Then, the performance of EQMOM is analyzed in the context of two validation cases: First, pure soot oxidation is studied in section 5.1. Second, in section 5.2, the model is validated for a laminar premixed flame with coupled soot growth and oxidation.

Why Soot Oxidation is a Challenge for Moment Methods

The statistics of a soot particle population are embodied in the Number Density Function (NDF). For a spatially homogeneous system, the governing PBE for the NDF, n (t, V ), based on particle volume, V , and dependent on time, t, reads @n (t, V ) @t +

X i @ @V [g i (t, V ) n (t, V )] = X j Ṡj . (1)
The summation over j on the r.h.s. includes the discontinuous source terms of nucleation, Ṡnucl , and coagulation, Ṡcoag , which can be expressed as integrals of the NDF. The summation over i on the l.h.s. includes all continuous source terms, i.e. surface reactions leading to particle growth and oxidation. The growth rates of these processes, g i , appear as convective velocities in phase space.

The volume-moments of the NDF are defined as

m k (t) = Z 1 0 V k n (t, V ) dV , ( 2
)
where k is the order of the moment. Transport equations for the moments can be derived by multiplication of Eq. 1 with V k and subsequent integration over the phase space:

dm k dt = X i V k g i (t, V ) n (t, V ) | 1 0 + X i Z 1 0 kV k 1 g i (t, V ) n (t, V ) dV + X j Z 1 0 V k Ṡj dV . (3)
On the r.h.s. of Eq. 3, there are integral terms as well as boundary terms resulting from the integration by parts of the hyperbolic growth terms (second term in Eq. 1). As the support interval of the NDF is semi-infinite, only the boundary term at the minimum particle size (zero in this case) is of concern.

If g i is positive, as in the case of surface growth, the term is zero, but if g i is negative, as in the case of oxidation, the boundary flux is non-zero and needs to be computed. Hence, as long as oxidation is excluded, the r.h.s. of Eq. 3 contains only integral terms. In presence of oxidation, the additional flux term appears.

The Method of Moments seeks the closure of the generally unclosed moment source terms, which enables the solution of Eq. 3. The time evolution of the moments represents the evolution of the statistics of the soot NDF. While the zeroth and first order moments are related to the soot number density and volume fraction, respectively, typically a few more moments are solved for in order to obtain a more accurate evaluation of the moment source terms.

The evaluation of the boundary flux is especially challenging as it requires the knowledge of the pointwise values of the NDF, which are not directly available from the moments.

For the following discussion, it is important to distinguish between the terms "representation of the NDF" and "approximation of the NDF". In classical quadrature-based moment methods, such as (D)QMOM, from the computed set of moments, a multi-delta function is reconstructed, which constitutes a representation of the NDF. The goal of the representation of the NDF is to accurately reproduce its statistics in terms of moments with the assumption that this will help to accurately predict the source terms in the moment equations, while the shape of the reconstructed NDF might be completely di↵erent from the true NDF. As long as oxidation is not considered, all source terms appear as integrals of the NDF, and using the (D)QMOM representation for their evaluation is equivalent to a Gaussian quadrature.

This method has been shown to be very accurate for soot growth [START_REF] Blanquart | Combustion Generated Fine Carbonaceous Particles[END_REF].

In case of oxidation, the additional boundary term cannot be closed with such a quadrature, because the pointwise flux at the minimum particle size needs to be known. Therefore, an approximation of the NDF is required, i.e. a reconstruction of the NDF which not only reproduces its statistics, but also accurately approximates its shape, especially close to the minimum particle size. Only from the knowledge of the number density of the smallest particles that will be oxidized during the next time step, the disappearance rate of particles can be computed. This quantity is not available in (D)QMOM, and modeling assumptions are required to compute the decrease in the number density due to oxidation.

In the original variant of MOMIC [4], a transport equation for the moment of order minus infinity is solved, and in HMOM [7], a delta function at the size of nucleating particles is introduced. Although the main purpose of adding these equations was not the treatment of oxidation, they can be used to determine the number of the smallest particles. However, the issue is shifted to the question of how many particles are transferred from larger size classes into the first bin. The pointwise knowledge of the NDF would be required to answer this question. However, this information is not available in interpolation-based or standard quadrature methods, and model assumptions need to be introduced.

In order to close this gap, EQMOM has recently been developed. Its moment inversion algorithm provides a continuous function as approximation of the NDF, which is used to evaluate the source terms and the oxidation boundary flux.

EQMOM for Soot Formation and Oxidation

Moment Inversion Algorithm

The concept of EQMOM relies on the replacement of the delta-functions in (D)QMOM with continuous kernel functions. The NDF is hence approximated by a superposition of N continuous kernel functions, (V ; V ↵ ), which are weighted by w ↵ :

n (V ) ⇡ N X ↵=1 w ↵ (V ; V ↵ ) . ( 4
)
The moment inversion then consists in reconstructing such kind of NDF from its moments of order 0 to 2N , i.e., finding the non-negative weights, w ↵ , the abscissas of the kernel functions, V ↵ , and their shape parameter, , which is identical for all kernel functions. Once these parameters are determined, the moment source terms in integral form can be computed with a so-called secondary quadrature, that is a quadrature for each kernel function [1]. The boundary term in Eq. 3 can be computed from the approximated NDF, as described in Sec. 3.4.

Di↵erent kinds of kernel functions are possible, which all satisfy the two following constraints: First, they formally tend to a Dirac delta function when tends to zero, thus recovering a quadrature. Second, there is a linear dependence between any moment vector, m, of a reconstructed NDF and the vector m ? of components m ? k =

P N ↵=1 w ↵ V k ↵ : m = A ( ) m ? . (5)
The form of A ( ) depends on the specific choice of kernel functions. This second constraint allows the use of the e cient and robust standard quadrature algorithm. Indeed, for any value of , one can compute m ? ( ) = A( ) 1 m and then the quadrature points (w ↵ ( ) , V ↵ ( )) N ↵=1 , thus leading to a reconstruction. The moments of order 0 to 2N 1 of this reconstruction are thus the given m k and the value of has then to be adapted in order for its 2N th order moment m2N ( ) to be m 2N . The algorithm for moment inversion is illustrated in Fig. 1.

A new robust and e cient algorithm for the reconstruction provided in [START_REF] Nguyen | Solution of population balance equations in applications with fine particles: mathematical modeling and numerical schemes[END_REF], improving the one of [1], is applied here. It uses an e cient way to solve the nonlinear problem J( ) = 0, where J( ) = m 2N m2N ( ). In some cases, a solution does not exist; then, J( ) 2 and hence the error on the last moment is minimized. Moreover, the algorithm is able to deal e ciently with moment vectors at the boundary of the space of realizable moments, where the only possible corresponding NDF is a single or multi-delta function. Such cases frequently appear in regions of a flame where soot nucleation has just started. A method then detects if the moment vector is close to the boundary of moment space. If this is the case, the algorithm adaptively switches to the standard quadrature algorithm, until coagulation has broadened the NDF.

Several types of kernel functions can be used, e.g. Gaussian distributions for internal coordinates with an unbounded support interval, or beta distributions if the support interval is bounded. As the support interval for soot particle volume is semi-infinite, [0, 1), two possible choices for the kernel functions are gamma [1] or lognormal [11] distributions. Both variants are applied and compared in this paper.

Gamma EQMOM

The kernel functions for gamma EQMOM [1] are given as

(V ; V ↵ ) = V ↵ 1 e V / ( ↵ ) ↵ , (6)
where ↵ = V ↵ / , and is the Gamma function. The corresponding moments can be analytically expressed as

m 0 = N X ↵=1 w ↵ , m k = N X ↵=1 w ↵ k 1 Y i=0 (V ↵ + i ) , k 1 . (7)
The matrix A ( ) is then triangular with coe cients [A( )] i,j = a i,j i j for 0  j  i, where a k,k = 1 and a k,i is given by the following recurrence formula [START_REF] Nguyen | Solution of population balance equations in applications with fine particles: mathematical modeling and numerical schemes[END_REF]:

a k,i = (k 1)a k 1,i + a k 1,i 1 , k 1, i = 2, . . . , k 1 . (8)

Lognormal EQMOM

The lognormal kernel functions are

(V ; V ↵ ) = 1 V p 2⇡ exp (ln (V ) ln (V ↵ )) 2 2 2 ! . ( 9
)
The moments take the form

m k = N X ↵=1 w ↵ V k ↵ exp ✓ 1 2 k 2 2 ◆ , (10)
leading to the matrix A ( ) being diagonal with diagonal coe cients exp 1 2 k 2 2 . Moreover, the nonlinear problem is then solved for exp ( 2 /2) instead of [11].

Treatment of Oxidation

Operator splitting is applied for the time integration of the moments. The soot growth and coagulation terms are integrated first. Then, oxidation is accounted for using an adaptation to an odd number of moments [START_REF] Nguyen | Solution of population balance equations in applications with fine particles: mathematical modeling and numerical schemes[END_REF] of the scheme by Massot et al. [START_REF] Massot | [END_REF]. The NDF is first reconstructed according to the EQMOM algorithm. Then, eliminating the disappearance flux, i.e. the part of the reconstructed distribution, which is oxidized during the current time step, the corresponding partial moments, are computed analytically (gray region in Fig. 2): The size of the largest particle that is oxidized in the time step, V Ox , can be evaluated from the oxidation rate.

e m k (t) = Z 1 V Ox V k n (t, V ) dV , k = 0, 1, . . . , 2N + 1 . (11)
An (N + 1)-point quadrature, ( e w i , e V i ), corresponding to these modified moments, { e m k } 2N +1 k=0 , is computed. The abscissas are then convected in phase space according to the oxidation rate, and the updated moments are computed as

m k (t + t) = N +1 X i=1 e w i e V k i (t + t) . (12)
A very important property of this algorithm is that it ensures the realizability of the moment set [START_REF] Massot | [END_REF][START_REF] Nguyen | Solution of population balance equations in applications with fine particles: mathematical modeling and numerical schemes[END_REF].

Soot Model

For the flames discussed in the following, detailed chemistry computations of the gas phase were performed using the FlameMaster code [START_REF] Pitsch | A C++ computer program for 0D combustion and 1D laminar flame calculations[END_REF] and the chemical kinetic mechanism of Narayanaswamy et al. [START_REF] Narayanaswamy | [END_REF], which was developed with special focus on soot precursors, and which contains PAH (polycyclic aromatic hydrocarbons) chemistry up to four-ringed molecules.

Two-way coupling between gas phase and particulate phase is applied to account for the removal of PAH from the gas phase during soot inception and growth.

The following soot processes are included in the computations of the particle phase: nucleation using the dimerization model by [START_REF] Blanquart | Combustion Generated Fine Carbonaceous Particles[END_REF], condensation of PAH dimers, coagulation in the transition regime, surface growth modeled with the HACA mechanism [20, 21], with reactions rates given by [START_REF] Blanquart | Combustion Generated Fine Carbonaceous Particles[END_REF], and oxidation by OH and O 2 , with rates given by [START_REF] Blanquart | Combustion Generated Fine Carbonaceous Particles[END_REF] and references therein.

Thermophoretic e↵ects are also included.

Results

For a thorough validation, model-predicted results are compared to Monte Carlo simulations as well as experimental data from three di↵erent experiments. First, we will consider two analytical cases for pure oxidation starting from realistic initial size distributions experimentally determined in [13,22].

Then, the model will be further validated with an experimental data set of laminar premixed flames [START_REF] Neoh | Soot burnout in flames[END_REF], focusing on oxidation.

Results for Pure Oxidation

First, soot oxidation is considered isolated from other processes. The goal of this validation case is to test if the EQMOM reconstruction yields a good approximation of realistic soot NDFs at all times during the oxidation process, which is necessary for the accurate prediction of the moments, especially of m 0 representing the number density.

The burner-stabilized stagnation (BSS) flame approach has been established as an experimental technique to measure the soot NDF in laminar premixed flames [22,[START_REF] Camacho | [END_REF]. Experiments in this configuration have shown that at downstream locations, bimodal shapes are observed for the soot NDFs in a rich C 2 H 4 /O 2 /Ar flame [22]. In a similar flame that was recently experimentally studied by the same research group [13], a unimodal NDF in the flame has been found. Both cases are interesting test cases, as both NDFs might be present at some point in a turbulent flame and be subjected to strong oxidation when reaching the flame front. Therefore, both cases will be considered here and used as initial conditions for the oxidation test cases described in the following.

Analytical Solution

Both gamma and lognormal EQMOM are applied using both experimentally measured soot NDFs discussed above (flame C3 in [22] at HAB = 8mm and C 2 H 4 flame in [13]) as initial conditions. For pure oxidation, if the initial NDF is given, the temporal evolution of the NDF and its moments can be obtained analytically from the oxidation law. As the goal is the validation of the statistical part of the model, a relatively simple chemical model is applied, and results are compared to the analytical solution. The oxidation rate is taken proportional to the particle surface, and particles are assumed to be spherical. The particle diameter is taken as internal coordinate of the NDF, such that oxidation simply shifts the initial NDF to smaller sizes, and the evolution of the NDF and its moments can easily be computed. As particles smaller than 2.5nm are below the detection limit of the particle sampler, a small part of the initial NDF is unknown. Therefore, the comparison of EQMOM and the analytical solution starts at the time when the smallest measured particle has reached zero size. These simplifications do not lead to a loss of generality regarding the conclusions drawn for the performance of the statistical model. Then, the oxidation algorithm of Sec. 3.4 is applied. Figure 3 shows the reconstructed initial NDFs using gamma and lognormal EQMOM with three kernel functions, i.e. seven moments need to be transported. The reconstructions of the NDFs at di↵erent times during the simulation are provided in the Supplementary Material.

Both the unimodal and the bimodal NDF can be very well approximated using two (not shown here) or three gamma distributions, while lognormal EQMOM is less accurate. The lognormal kernels do not overlap very much, which leads to a bimodal shape for both NDFs; also the experimentally bimodal NDF is not well approximated in the region of small particles. Lognormal EQMOM has di culties to capture non-zero values at the minimum particle size, because the lognormal distribution always starts at zero. In gamma EQMOM, depending on the parameters of the gamma distributions, a smooth transition occurs between the NDF starting at zero and at a nonzero value. Although the parameter is identical for all kernels, the shapes of the kernels can di↵er in that the first kernel starts at a non-zero value, while the others start at zero. This is a results of the definition of and the parameter of the gamma distribution, ↵ = V ↵ / , and it enables an accurate approximation of the NDF, especially for small particle sizes, which is important for an accurate prediction of the number density during oxidation. Using two or more kernels, i.e. five or more moments, both gamma and lognormal EQMOM can very accurately predict the soot volume fraction and higher moments. Gamma EQMOM also excellently predicts the evolution of the number density. Due to the bimodality of the NDF, it first decreases fast until the peak of small particles has been oxidized, then the rate slows down, before increasing again when the second peak is oxidized. This behavior is qualitatively and quantitatively very well captured.

For modeling the disappearance rate of particles during oxidation using traditional moment methods, it is often assumed that a particle is removed after the mass of an average-sized particle has been oxidized. This model assumption couples the rate of change of number density to the rate of change of volume fraction. This coupling is obviously wrong, especially for the typical bimodal NDFs, and can only be avoided if the pointwise values of the NDF are known with su cient accuracy, as is the case in gamma EQMOM.

Lognormal EQMOM produces jumps in the number density. This behavior is linked to the evolution of the parameters of the kernel functions.

Figure 5 shows the weights and abscissas, and Fig. 6 shows the evolution of the shape parameter, . While the weights of the gamma distributions follow smooth curves, the weights of the lognormal distributions show jumps, whose positions are dependent on the time step size (not shown here). During oxidation, the first lognormal kernel becomes very narrow, until it disappears in a certain time step, which leads to the jump in the number density. The gamma distributions are always broad, and, as discussed above, the kernels have di↵erent parameters ↵ = V ↵ / for a common value of , such that the first kernel starts at a non-zero value, while the other kernels start at zero.

This behavior helps to avoid the sudden disappearance of a kernel and enables an accurate prediction of the number density, which is also independent of the time step size.

Another benefit of gamma EQMOM over lognormal EQMOM is evident from Fig. 7, which shows the time evolution of m 1 /m 0 , corresponding to a mean particle diameter. Gamma EQMOM is able to accurately predict the mean particle diameter even for large times, when both m 1 and m 0 go to zero. For lognormal EQMOM, in contrast, the convergence with increasing number of transported moments is much slower, and for large times, the monotonically decreasing behavior is not captured any more.

Results for Coupled Soot Formation and Oxidation

The two-stage burner experiment by Neoh et al. [START_REF] Neoh | Soot burnout in flames[END_REF] is simulated to validate the EQMOM algorithm for combined soot formation, growth, coagulation, and oxidation. In this experiment, soot is produced in a rich, premixed, burner-stabilized CH 4 /air flame. Then, secondary air is added, and the soot is oxidized in a second premixed flame (case with global = 1.15 in [START_REF] Neoh | Soot burnout in flames[END_REF] is studied here). Although soot formation and oxidation are dominant in the primary and secondary burner, respectively, all soot processes are relevant and considered in both flames. As the flames are laminar and one-dimensional, the spatial coordinate is transformed into a pseudo-time or Lagrangian time using the particle velocity, i.e. the sum of gas and thermophoretic velocities.

The particulate phase is then integrated in time. The simulation results of the primary burner are used as initial conditions for the secondary burner.

Monte Carlo Simulations

As an analytical solution for the NDF is no longer possible for this case, EQMOM predictions are here compared to Monte Carlo simulations using the same physico-chemical soot model. In addition, model results are compared to experimental data for the secondary burner.

MC simulations are particularly useful for validation of statistical models, since they require no closure assumption and hence describe the exact NDF evolution for a given physico-chemical model. The MC code uses the algorithm by [3,[START_REF] Singh | Modelling soot particle size distribution: Dynamics of pressure regimes[END_REF] for coagulation in the transition regime. To improve computational e ciency of the MC simulations, the method of majorant kernels [START_REF] Goodson | [END_REF] is used for the continuum regime.

Validation: Soot Formation and Oxidation in a Two-Stage Burner

In Fig. 8, the soot volume fraction in the primary burner is shown on the left. Very good agreement between EQMOM and MC results is obtained. To highlight the e↵ect of oxidation even in this rich flame ( = 2.1), additional simulations excluding oxidation are also shown. As for both simulations the agreement between EQMOM and MC is very good, it can be concluded that EQMOM is suitable to accurately predict isolated soot formation and coagulation as well as simultaneous oxidation. The prediction of the experimental number density is also reasonable, given the experimental uncertainty. The overprediction of the number density might be a result of an underprediction of the coagulation rate in the primary burner, which could be improved with a more sophisticated, multivariate soot model including a model for aggregation. It should also be noted that a significant experimental uncertainty can be expected, as Neoh [START_REF] Neoh | Soot burnout in flames[END_REF] computed the number density using extinction and scattering measurements, and making assumptions about the particle size distribution (monodisperse or single lognormal). More importantly, the EQMOM results for the number density agree well with the MC simulations, which validates the EQMOM approach for coupled soot formation and oxidation.

Conclusions

The EQMOM algorithm has been integrated with a physico-chemical soot model including soot inception, growth, coagulation, and oxidation processes.

Extending existing EQMOM implementations, the algorithm has for the first time been combined with a treatment for particle removal by oxidation. Two variants of EQMOM applicable to NDFs with a semi-infinite support interval, gamma and lognormal EQMOM, have been applied to soot oxidation using realistic soot NDFs taken from two di↵erent experiments as initial conditions.

While both methods are able to accurately predict the soot volume fraction and higher order moments of the NDF, lognormal EQMOM has di culties to predict the soot number density. On the contrary, gamma EQMOM using at least two kernel functions yields an excellent approximation of the NDF, and it is therefore a suitable method for the accurate prediction of both soot volume fraction, number density, and the soot NDF itself. A subsequent application of gamma EQMOM to the two-stage burner experiment revealed that this is also true for coupled soot formation and oxidation in flames.

Figure 1 :

 1 Figure 1: Schematic of the EQMOM algorithm.

Figure 2 :

 2 Figure 2: Schematic illustration of the oxidation treatment.

Figure 3 :

 3 Figure 3: Reconstruction of two experimentally measured, normalized NDFs from rich, premixed ethylene flames (left: [13]; right: [22]) using gamma EQMOM (upper row) and lognormal EQMOM (lower row) with three kernel functions.

Figure 4 :

 4 Figure 4: Time evolution of the first three moments, representing number density, soot volume fraction, and variance of the NDF, normalized with the initial particle number, during pure oxidation of the soot population represented by the bimodal NDF in Fig. 3.

Figure 4

 4 Figure 4 shows the time evolution of m 0 , m 1 , and m 2 , representing number density, soot volume fraction, and variance of the NDF, respectively, during

Figure 7 :

 7 Figure 7: Time evolution of m 1 /m 0 . Same case as in Fig. 4.

Figure 8 :

 8 Figure 8: Soot volume fraction in the primary (left) and secondary burner (middle), and number density in the secondary burner (right). Comparison of EQMOM with MC,and with experimental data in the secondary burner from[START_REF] Neoh | Soot burnout in flames[END_REF], where two postprocessing techniques were used to determine the number density.

Figure 8

 8 Figure 8 also shows the evolution of soot volume fraction and number density in the oxidation-dominated secondary burner. Comparisons are made between EQMOM, MC and experimental data. It should be noted that a deviation of EQMOM results from experimental results is always due to a combination of uncertainties in the physico-chemical and the statistical model, while the comparison of EQMOM with MC simulations isolates the statistical error, which is of prime concern here. Regarding the soot volume fraction, EQMOM results are in excellent agreement both with MC simulations and experimental data.

Acknowledgements A.W. and H.P. gratefully acknowledge funding by the German Research Foundation (DFG) under grant no. PI 368/6-1 and by the German Research Association for Combustion Engines e.V. (FVV) under grant no. 1166. T.T.N., F.L. and R.O.F. gratefully acknowledge funding by the French National Research Agency (ANR) under grant ANR-13-TDMO-02 ASMAPE for the ASMAPE project.

Reconstruction of the NDF during Oxidation

Fig. 3 in the main paper shows the reconstruction of both the unimodal and the bimodal NDF in the first time step using both gamma and lognormal EQMOM with three kernel functions. These results are supplemented here by the reconstructions of the NDF at several times during the oxidation using both gamma and lognormal EQMOM with two and three kernels.

S1

Supplemental Material

Supplementary Material

Modeling Soot Oxidation with the Extended Quadrature Method of Moments Achim Wick ? , Tan-Trung Nguyen, Frédérique Laurent, Rodney O. Fox, Heinz Pitsch

⇤

Corresponding author: a.wick@itv.rwth-aachen.de

Supplementary Material

Modeling Soot Oxidation with the Extended Quadrature Method of Moments Achim Wick ? , Tan-Trung Nguyen, Frédérique Laurent, Rodney O. Fox, Heinz Pitsch

⇤

Corresponding author: a.wick@itv.rwth-aachen.de