Jfpj Abascal
email: juan.abascal@creatis.insa-lyon.fr

N Ducros

F Peyrin

Nonlinear material decomposition using a regularized iterative scheme based on the Bregman distance

Keywords: Bregman distance, convexity, spectral computed tomography

In this paper, we address the resolution of material decomposition, which is a nonlinear inverse problem encountered in spectral computed tomography (CT). The problem is usually solved in a variational framework but, due to the nonlinearity of the forward operator, the objective function may be nonconvex and standard approaches may fail. Regularized iterative schemes based on the Bregman distance have been suggested for improving global convergence properties. In this work, we analyze the convexity of the material decomposition problem and propose a regularized iterative scheme based on the Bregman distance to solve it. We evaluate our Bregman iterative algorithm and compare it with a regularized Gauss-Newton (GN) method using data simulated in a realistic thorax phantom.

First, we prove the existence of a convex set where the usual data fidelity term is convex. Interestingly, this set includes zero, making it a good initial guess for iterative minimization schemes. Using numerical simulations, we show that the data fidelity term can be nonconvex for large values of the decomposed materials. Second, the proposed Bregman iterative scheme is evaluated in different situations. It is observed to be robust to the selection of the initial guess, leading to the global minimum in all tested examples while the GN method fails to converge when the initial guess is not well chosen. Moreover, it is found to avoid the selection of the regularization parameter for little extra computation.

In conclusion, we have provided a suitable initialization strategy to solve the nonlinear material decomposition problem using convex optimization methods and evaluated a Bregman iterative scheme for this problem. The improvement in global convergence of Bregman iterative scheme combined with other interesting properties of the Bregman distance appears as a compelling strategy for nonlinear inverse problems.

Introduction

We consider regularized methods for solving nonlinear inverse problems. Successful methods for linear inverse problems have been generalized for the case of nonlinear inverse problems. However, in general, objective functions involving nonlinear operators are nonconvex, so classical optimization procedures may fail. Regularized iterative procedures have been proposed to ensure that a global minimizer exists. The idea behind these methods is to replace the original regularized problem with a sequence of problems, each of which satisfies stronger convergence properties [START_REF] Ramlau | Tikhonov replacement functionals for iteratively solving nonlinear operator equations[END_REF][START_REF] Ramlau | A Tikhonov-based projection iteration for nonlinear Ill-posed problems with sparsity constraints[END_REF].

A popular and efficient iterative regularization procedure for the 1-and total variation (TV)-problems is based on the use of the Bregman distance and has been widely studied for linear measurements [START_REF] Osher | An Iterative Regularization Method for Total Variation-Based Image Restoration[END_REF][START_REF] Yin | Bregman Iterative Algorithms for l1-Minimization with Applications to Compressed Sensing[END_REF][START_REF] Cai | Linearized Bregman iterations for compressed sensing[END_REF][START_REF] Cai | Split Bregman Methods and Frame Based Image Restoration[END_REF][START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF][START_REF] Osher | Fast linearized Bregman iteration for compressive sensing and sparse denoising[END_REF]. Replacing a convex regularization functional with the Bregman distance associated with this functional leads to an iterative scheme that converges to the solution of the constrained optimization problem. The use of the Bregman distance provides additional appealing properties. For instance, it has been shown that minimizing the Bregman distance associated with TV can recover the contrast lost that may appear when solving the TV problem [START_REF] Burger | A Guide to the TV Zoo[END_REF]. The split Bregman method proposed for image processing and compressed sensing combines the Bregman distance with a splitting technique to efficiently deal with non differentiable convex functional, such as total variation [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF][START_REF] Esser | Applications of Lagrangian-based alternating direction methods and connections to split Bregman[END_REF][START_REF] Wu | Augmented Lagrangian Method, Dual Methods, and Split Bregman Iteration for ROF, Vectorial TV, and High Order Models[END_REF]. It has been successfully applied to solve the total variation problem in image denoising, image impainting, and compressed sensing [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF][START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF][START_REF] Abascal | Fluorescence diffuse optical tomography using the split Bregman method[END_REF][START_REF] Gao | Multi-energy CT based on a prior rank, intensity and sparsity model (PRISM)[END_REF][START_REF] Plonka | Curvelet-Wavelet Regularized Split Bregman Iteration for Compressed Sensing[END_REF][START_REF] Abascal | Comparison of total variation with a motion estimation based compressed sensing approach for self-gated cardiac cine MRI in small animal studies[END_REF][START_REF] Juan | Incorporation of prior knowledge of the signal behavior into the reconstruction to accelerate the acquisition of MR diffusion data[END_REF]. The use of Bregman distance for the TV minimization subject to nonlinear measurements enjoys further benefits but it has not been widely explored [START_REF] Gao | Quantitative Photoacoustic Tomography[END_REF][START_REF] Gholami | Regularization of linear and non-linear geophysical ill-posed problems with joint sparsity constraints[END_REF]. In [START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF], Bregman iterative schemes based on classical reconstruction methods were proposed to improve global convergence for nonlinear inverse problems.

In this work, we are concerned with the nonlinear material decomposition problem that appears in spectral computed tomography (CT). Spectral CT is gaining increasing attention due to the recent development of energy discriminating photoncounting detectors that simultaneously count photons and resolve their energy [START_REF] Taguchi | Vision 20/20: Single photon counting x-ray detectors in medical imaging[END_REF][START_REF] Heismann | Spectral computed tomography[END_REF]. This technology addresses the limitations of conventional CT and provides material decomposition capabilities, which allows, for instance, characterizing lesions for breast imaging [START_REF] Ding | Breast tissue decomposition with spectral distortion correction: A postmortem study[END_REF] and atherosclerotic plaque [START_REF] David P Cormode | Atherosclerotic Plaque Composition: Analysis with Multicolor CT and Targeted Gold Nanoparticles1[END_REF][START_REF] Boussel | Photon counting spectral CT component analysis of coronary artery atherosclerotic plaque samples[END_REF].

The spectral CT material reconstruction problem can be solved in one step by recovering the material mass densities directly from the energy projections [START_REF] Long | Multi-Material Decomposition Using Statistical Image Reconstruction for Spectral CT[END_REF][START_REF] Zhang | Model-Based Iterative Reconstruction for Dual-Energy X-Ray CT Using a Joint Quadratic Likelihood Model[END_REF][START_REF] Zhao | An extended algebraic reconstruction technique (E-ART) for dual spectral CT[END_REF][START_REF] Barber | An algorithm for constrained one-step inversion of spectral CT data[END_REF]. A two-step alternative is first to decompose the energy projections into material projections and then to perform a tomographic reconstruction from the material projections. The one-step approach can incorporate prior knowledge about the material density -the quantity of interest-while determining prior knowledge for the projected mass densities of the two-step material decomposition is less natural. However, the two-step approach has several advantages. First, it enables to split a large inversion problem into smaller problems. Second, it offers a physics-based parallelization scheme where all projection views are decomposed in parallel and all material are reconstructed in parallel. Third, it separates linear (i.e., CT reconstruction) from nonlinear (i.e., material decomposition) mixing, allowing to reuse state-of-the art CT reconstruction approaches. In this paper, we focus on the material decomposition problem and adopt a model-based approach which models source spectrum and energy responses of the detector [START_REF] Schlomka | Experimental feasibility of multi-energy photoncounting K-edge imaging in pre-clinical computed tomography[END_REF]. For alternative approaches to material reconstruction see the following works [START_REF] Gao | Multi-energy CT based on a prior rank, intensity and sparsity model (PRISM)[END_REF][START_REF] Zhao | Tight-frame based iterative image reconstruction for spectral breast CT[END_REF][START_REF] Darin | Spectral diffusion: an algorithm for robust material decomposition of spectral CT data[END_REF][START_REF] Persson | Energy-resolved CT imaging with a photon-counting silicon-strip detector[END_REF][START_REF] Semerci | Tensor-Based Formulation and Nuclear Norm Regularization for Multienergy Computed Tomography[END_REF][START_REF] Kim | Sparse-view spectral CT reconstruction using spectral patchbased low-rank penalty[END_REF][START_REF] Xi | United Iterative Reconstruction for Spectral Computed Tomography[END_REF][START_REF] David | Joint reconstruction of multi-channel, spectral CT data via constrained total nuclear variation minimization[END_REF][START_REF] Touch | A neural networkbased method for spectral distortion correction in photon counting x-ray CT[END_REF]. We remark that most of these approaches adopted some linearized forms or focused on the linear tomographic problem.

Most previous studies that addressed spectral CT have ignored nonconvexity or adopted some assumptions to make the problem convex [START_REF] Liu | Material reconstruction for spectral computed tomography with detector response function[END_REF]. For dual-energy CT and for the negative log-likelihood functional, a hybrid first and second order variational approach was proposed to provide a more robust algorithm [START_REF] Brendel | Penalized likelihood decomposition for dual layer spectral CT[END_REF]. Recently, a regularized Gauss-Newton (GN) method with specific regularization applied to each material was applied to solve the material decomposition problem [START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF]. This method provided significantly superior results than the non regularized maximum-likelihood method that is commonly used in practice [START_REF] Schlomka | Experimental feasibility of multi-energy photoncounting K-edge imaging in pre-clinical computed tomography[END_REF].

The goal of this work is to investigate convexity for the material decomposition inverse problem and to propose a regularized Bregman iterative scheme to solve this problem. To this goal we build from the previously proposed regularized GN method [START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF] and modify it in order to minimize the Bregman distance associated with the same regularization functional. We put particular emphasis to study convexity and to provide a suitable initialization strategy for convex optimization algorithms. In addition, we assess the benefits of a Bregman iterative algorithm to provide a global minimum by studying its dependence on the initial guess.

The paper is organized as follows. In section 2, we present the formulation of the forward and inverse problem. In section 3, we show that the data fidelity term is convex in a defined set and that the zero solution belongs to this region. In section 4, we describe the numerical approaches adopted to solve this problem, the regularized Bregman iterative scheme and a regularized Gauss-Newton method. In section 5, we show the numerical results. First, we analyze convexity using numerical examples. Then, we assess the Bregman iterative scheme in comparison with the Gauss-Newton approach using a numerical thorax phantom. Section 6 contains the discussion, further work and conclusion.

Material decomposition in spectral CT

Spectral CT

Let s i,θ ∈ R P represent the photon counts measured in the i-th energy bin of the detector under the θ-th projection angle and ρ m ∈ R N represent the mass density of the m-th material in the object. Spectral CT consists in recovering the mass densities of the M materials in the object from the photon counts measured in I energy bins for Θ views. Solving the spectral CT problem in one step requires the inversion of the nonlinear

system s = G(ρ) (1)
where s = (s 1,1 , . . . , s I,Θ ) is the full measurement vector, ρ = (ρ 1 , . . . , ρ M ) is the material vector, and G models both spectral mixing and tomographic acquisition. Alternatively, the two-step approach [START_REF] Schlomka | Experimental feasibility of multi-energy photoncounting K-edge imaging in pre-clinical computed tomography[END_REF][START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF] first decomposes all projection views by inverting

s θ = F(a θ ), 1 ≤ θ ≤ Θ, (2)
where s θ = (s 1,θ , . . . , s I,θ ), a θ = (a 1,θ , . . . , a M,θ ) represents the projected mass densities for the θ-th view, and F models spectral mixing only. This material decomposition step is followed by the tomographic reconstruction of all materials, i.e., the inversion of

a m = X (ρ m ), 1 ≤ m ≤ M, (3)
where a m = (a m,1 , . . . , a m,Θ ) and X is the X-ray transform. While the one-step approach requires inverting the R M N → R IΘP nonlinear mapping defined in (1), the two-step approach split it into multiple inversions of R M P → R IP mappings followed by multiple inversions of R N → R P mappings. The two-step approach not only reduces the size of the problems, offering a natural parallelization scheme (e.g., decomposing all angles in parallel, reconstructing all materials in parallel), but also separate nonlinear (i.e., spectral) from linear (i.e, tomographic) mixing. In the following, we focus on the material decomposition problem defined in (2).

Spectral mixing

We now consider a fixed θ in (2) and drop the θ index for the sake of simplicity. Let s p i represents the photon counts in the p-th pixel of the detector and a p m the projected mass density in the same pixel. Assuming that no inter-pixel spectral mixing occurs, the system (2) is equivalent to

s p i = F p i (a p 1 , . . . , a p M ), 1 ≤ i ≤ I and 1 ≤ p ≤ P (4)
where s i = (s 1 i , . . . , s P i ) and a m = (a 1 m , . . . , a P m ). The nonlinear mapping F p i : R M → R can be modelled by [START_REF] Schlomka | Experimental feasibility of multi-energy photoncounting K-edge imaging in pre-clinical computed tomography[END_REF]

F p i (a p 1 , . . . , a p M ) := E n 0 (E)d i (E) exp - M m=1 a p m τ m (E) dE (5)
where E = [0, E max ] is the energy range delivered by the X-ray tube, n 0 (E) is the X-ray source energy spectrum, d i (E) is the energy-dependent response function of the i-th bin of the detector, and τ m (E) is the energy-dependent mass attenuation of the m-th constituent of the object.

Inverse problem

Within a variational framework, material decomposition is formulated as the following optimization problem min a C(a) = min a L(s δ , a) + αR(a) [START_REF] Cai | Split Bregman Methods and Frame Based Image Restoration[END_REF] where s δ represents the noisy measurements, L is a data fidelity term, R is a regularization term encoding a priori knowledge about the solution, and α is the regularization parameter. Assuming that the measurements are Poisson distributed, the following approximation of the discretised Kullbach-Leibler divergence is considered

L(s δ , a) = 1 2 F(a) -s δ 2 W , (7)
where x W is the weighted norm defined by x 2 W = x W W x. Here, the weight matrix W is chosen as

W = diag 1 max ( √ s δ , ) , (8)
where is a small parameter that is added to the denominator to avoid division by zeros.

It is important to note that although the regularization functional R is convex the cost function C in (6) can be nonconvex whenever the data fidelity term L is nonconvex, which usually occurs with nonlinear measurements.

Gauss Newton method

A wide variety of approaches can be used to solve the nonlinear problem [START_REF] Cai | Split Bregman Methods and Frame Based Image Restoration[END_REF] [START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Arridge | Optical tomography: forward and inverse problems[END_REF]. Among them the Gauss Newton method has fast (quadratic) convergence and has been recently shown to be very effective for material decomposition [START_REF] Brendel | Penalized likelihood decomposition for dual layer spectral CT[END_REF][START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF]. The Gauss Newton update is computed as

a n+1 = a n -λ n H -1 (a n )C (a n ) ( 9
)
where λ n is the step length, C is the gradient of C, and H is an approximation for the Hessian of C, which is given by

H(a n ) = F (a n ) W W F (a n ) + αR (a n ) ( 10
)
where F denotes the Jacobian of F and R the Hessian of R. The step length is usually found by a line search strategy, i.e., minimizing the cost function along the descent direction provided by H -1 (a n )C (a n ).

Analysis of the convexity of the problem

In this section, we study the convexity of the material decomposition problem. In particular, we show that the weighted least squares data fidelity term is convex on a convex set that includes a = 0. This provides an initialization point that guarantees that unconstrained optimization algorithms converge to the global minimizer. Convexity of the problem is analyzed deriving the following properties for the Hessian of the data fidelity functional. For simplicity, the analysis is conducted considering one pixel only (i.e. for P = 1) and the index p has been dropped. Since the data fidelity term L is separable with respect to pixels, the results bellow can be generalized to any pixel number (i.e., P > 1). Definition 1. Let g i (a, E) be the number of incident photons seen at energy E in the i-th energy bin of the detector. It is given by

g i (a, E) = n 0 (E)d i (E) exp - M m=1 a m τ m (E) (11)
Proposition 1. The Hessian of L is given by

H(a) = - i w 2 i (s δ i -F i (a))T G i (a)T + F W W F ( 12
)
where T = (τ j,m ) ∈ R J×M is the mass attenuation matrix defined by

τ j,m = τ m (j∆E), 1 ≤ j ≤ J and G i (a) = diag (g i (a, E 1 ), . . . , g i (a, E J )) ∆E.
Proof. The functions F i : R M → R, 1 ≤ i ≤ I are twice continuously differentiable and the Hessian of L takes the following matrix form [START_REF] Vavasis | Nonlinear Optimization[END_REF]

H(a) = - i w 2 i (s δ i -F i (a))H i (a) + F W W F (13)
where H i denotes the Hessian of F i . It can be shown that H i (a) = T G i (a)T (see section Appendix A), which completes the proof.

Proposition 2. The data fidelity functional L is convex on the (possibly nonconvex) set A defined by

A = {a ∈ R M | F i (a) ≥ s δ i , ∀i} (14)
Proof. While F W W F is clearly positive semi-definite, some or all the terms

w 2 i (s δ i -F i (a))T G i (a)T can contribute against convexity. Since G i (a) is a diagonal matrix with positive elements, T G i (a)T is positive definite. If s δ i -F i (a)
is negative for all i, then H(a) is the sum of positive definite matrices, hence it is positive definite. Proposition 3. (Exact Data) Let a true be the solution and s δ = s true = F(a true ). The data fidelity functional L is convex on the convex set B defined by

B = {a ∈ R M | a m ≤ a true m , ∀m} (15)
Proof. Noting that all F i : a m → F i (a m ) are monotonically decreasing, choosing a m ≤ a true m for all m guaranties that F i (a) ≥ F i (a true ).

Remark 4. (Exact Data)

The particular point a = 0 lies within the convex region, i.e., 0 ∈ B.

Proof. Physically, solutions are positive.

Resolution of the problem

Iterative regularization scheme based on the Bregman distance

The Bregman distance generalizes the concept of distance for convex regularization functionals R [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF]. The generalized Bregman distance associated with the functional R for two points a and ã is defined as

D R ξ (a, ã) = R(a) -R(ã) -ξ, a -ã , ( 16
)
where ξ is a subgradient of R at ã. Even D R ξ (a, ã) ≥ 0, it is not a distance in a strict sense as in general it is not symmetric and it does not satisfy the triangle inequality.

The use of Bregman distance became popular after its application to construct iterative regularization schemes for total variation regularization [START_REF] Osher | An Iterative Regularization Method for Total Variation-Based Image Restoration[END_REF]. Since then, it has been widely applied to image restoration, deconvolution and compressed sensing problems with linear measurement operators and with convex regularization functional [START_REF] Yin | Bregman Iterative Algorithms for l1-Minimization with Applications to Compressed Sensing[END_REF][START_REF] Cai | Linearized Bregman iterations for compressed sensing[END_REF][START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF][START_REF] Burger | A Guide to the TV Zoo[END_REF][START_REF] Abascal | Fluorescence diffuse optical tomography using the split Bregman method[END_REF][START_REF] Abascal | A novel prior-and motion-based compressed sensing method for small-animal respiratory gated CT[END_REF]. A generalization of the Bregman distance for nonconvex regularization functionals has been also suggested [START_REF] Grasmair | Generalized Bregman distances and convergence rates for non-convex regularization methods[END_REF]. However, its application to solve problems with nonlinear measurements is less known.

In [START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF], authors proposed to replace the problem ( 6) by a sequence of problems

a k+1 = arg min a 1 2 F(a) -s δ 2 W + αD R ξ k (a, a k ), (17a)
ξ k+1 = ξ k -α -1 (W F (a k+1 )) (W F(a k+1 ) -W s δ ), (17b)
where α > 0 is an hyperparameter, which is not a regularization parameter and that can be chosen as a large constant. The iterative scheme [START_REF] Juan | Incorporation of prior knowledge of the signal behavior into the reconstruction to accelerate the acquisition of MR diffusion data[END_REF] turns to be a natural extension of the Bregman iterative method for linear measurements [START_REF] Osher | An Iterative Regularization Method for Total Variation-Based Image Restoration[END_REF]. Each of the subproblems (17a) is an optimization problem that has the same form as ( 6) but where the regularization functional R is replaced by the Bregman distance between the unknown, a, and its previous estimate, a k . Equation (17b) provides un update of the subgradient derived from the first order optimality condition of (17a), which enforces the data fidelity term iteratively [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF]. The Bregman iterative scheme [START_REF] Juan | Incorporation of prior knowledge of the signal behavior into the reconstruction to accelerate the acquisition of MR diffusion data[END_REF] constitutes then a sequence of subproblems where each of them has the same complexity as the original problem (6), but there is a substantial gain, as each subproblem is locally convex for sufficiently large values of α [START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF]. This provides a scheme to find a global minimizer for nonconvex functional associated with nonlinear operator equations.

An additional modification to the iterative scheme [START_REF] Juan | Incorporation of prior knowledge of the signal behavior into the reconstruction to accelerate the acquisition of MR diffusion data[END_REF] was proposed in [START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF] to ensure local convexity for total variation regularization. A zero-order Tikhonov term κ 2 a 2 2 is added to the original total variation regularization functional R to ensure local convexity around the current estimate a k . The value of κ must be sufficiently small in order not to modify the properties of the regularization functional R. Substituting the Bregman distance ( 16) into (17a), the iterative scheme is rewritten as

a k+1 = min a Ψ k (a) = min a 1 2 F(a) -s δ 2 W + αR(a) -α < ξ k , a > +α κ 2 a 2 2 , (18a) ξ k+1 = ξ k -α -1 (W F (a k+1 )) (W F(a k+1 ) -W s δ ). (18b)
A key point in the construction of these iterative schemes is to select efficient methods for solving the subproblem (18a) [START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF]. We propose to solve (18a) using the Gauss-Newton method ( 9)- [START_REF] Esser | Applications of Lagrangian-based alternating direction methods and connections to split Bregman[END_REF], which results in Algorithm 1. We refer to our approach as GN-Bregman (GNB). MATLAB code for material decomposition using GN and GNB is available from [START_REF]Spectral X-ray image reconstruction toolbox[END_REF].

Algorithm 1 Regularized Gauss Newton Bregman (GNB) Algorithm Set α, κ, tol, ψ, a 0 k ← 1 while 1 2 F(a k ) -s δ 2 W > tol do Outer iteration a k ← a k-1 for n = 1 to N k do Inner iteration H k ← (W F (a k )) T W F (a k ) + αR (a k ) + ακI g k ← (W F (a k )) T W (F(a k ) -s δ ) + αR (a k ) -αξ k + ακa k p ← -(H k ) -1 g k λ ← arg min λ Ψ(a k + λp) a k ← a k + λp ξ k ← ξ k -α -1 (W F (a k )) T (W F(a k ) -W s δ ) Subgradient update k ← k + 1

Algorithm stopping criterion

Note that GNB (algorithm 1) numerically involves two iterative procedures: the (outer) Bregman iteration indexed by k and the (inner) Gauss Newton iteration that we denote with the index n. We define N k the total number of GN iterations required to minimize Ψ k and N tot = k N k the total number of GN iterations of the GNB algorithm. The stopping criterion for the minimization of Ψ k is based on the relative cost function decrease, i.e, the algorithm is stopped at N k iterations when the relative cost function is smaller than ψ. The stopping criterion for the Bregman iterative scheme is the discrepancy principle, i.e., the algorithm is stopped when 1 2 F(a k ) -s δ 2 W ≤ tol [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF][START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF].

Selection of the initial guess

Proposition 3 and remark 4 provide a suitable initialization point for convex optimization algorithms. In particular for GNB, setting a 0 = 0 and ξ 0 = 0 ∈ ∂R(a 0 ) leads to

D ξ 0 R (a, a k ) = R(a) and Ψ 0 (a) = L(a, s) + αR(a)
, so the first inner problem in GNB (algorithm 1) is equivalent to the original unconstrained problem (6) but with a large regularization parameter. As a 0 belongs to the convex region, the first iterate a 1 is an excessively regularized solution that lies in the convex region and is closer to a true . In [START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF] it is shown that when solving subproblem (17a) exactly, iterates a k satisfy monotonicity of the residuals and of the error, so a k for k > 1 remain in the convex region and get monotonically closer to a true .

The Bregman iteration enjoys additional properties that improve global convexity for nonlinear inverse problems. In section 5.4, we show using numerical simulations that GNB converges to an optimal solution even if the initial guess is on the nonconvex region.

Regularization functional

In this work, we consider the separable regularizer proposed in [START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF] and defined by

R(a) = m R m (a m ), ( 19
)
where R m is the regularization functional for the m-th material. This choice allows to incorporate a-priori information for each material independently. Here, we consider Tikhonov regularization to promote smooth materials, i.e.,

R m (a m ) = D a m 2 2 (20)
where D represents the -th spatial derivative operator, and a smooth approximation of total variation to promote piecewise constant materials, i.e.,

R m (a m ) = (Da m ) 2 + 2 -. (21)
where is a small hyperparameter. Although only differentiable regularizers are used, Bregman iterative algorithms are also valid for the space of functions of bounded variations and other convex sparse regularization functional that uses 1 -norm to impose sparsity in a given orthonormal basis, such as wavelet transform or framelets.

Numerical simulations and results

Numerical phantom data

We used numerical simulations to illustrate the (non) convexity of the problem and evaluate the proposed algorithms. Different datasets were simulated using realistic X-ray tube and photon-counting detector with four energy bins (see figure 1) and considering numerical phantoms made of soft tissue, bone and gadolimium (Gd). This is a relevant example in K-edge imaging for characterizing atherosclerotic plaque [START_REF] David P Cormode | Atherosclerotic Plaque Composition: Analysis with Multicolor CT and Targeted Gold Nanoparticles1[END_REF][START_REF] Boussel | Photon counting spectral CT component analysis of coronary artery atherosclerotic plaque samples[END_REF]. The linear mass attenuation coefficients for the three materials are provided in figure 1c. Data for the thorax phantom are publicly available from SPRAY (Spectral X-ray image reconstruction) toolbox, a MATLAB toolbox for simulation and nonlinear image reconstruction in spectral X-ray imaging [START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF][START_REF]Spectral X-ray image reconstruction toolbox[END_REF].

Two numerical phantoms were created: a one-pixel phantom (only used to study convexity) and a realistic thorax phantom built from a segmented thoracic-abdominal CT scan [START_REF]3D Image Reconstruction for Comparison of Algorithm Database[END_REF][START_REF] Kéchichian | Shortest-Path Constraints for 3D Multiobject Semiautomatic Segmentation Via Clustering and Graph Cut[END_REF]. For the realistic phantom, a 611 × 167-pixel detector is simulated, leading to P = 102, 037. Two projection images for two different projection angles were evaluated. For each view, the dataset comprised I = 4 projection images corresponding to the different energy bins and M = 3 unknowns projected mass density images were recovered. Decreasing measurement signal-to-noise ratios were considered by lowering the number of incident photons from 10 8 to 10 5 . To get away from the inverse crime and further asses the robustness of our algorithm, we considered the mixed noise model s δ ∼ P(λ = s) + N (µ, σ 2 ), where P is the Poisson distribution and N the normal distribution whose mean µ can be calibrated to be zero [START_REF] Zeng | A Simple Low-Dose X-Ray CT Simulation From High-Dose Scan[END_REF]. We used a concentration of gadolinium equal to 1 g•cm -3 unless otherwise stated.

We have used specific regularization for each material as described in section 4.4. We used first-and second-order Tikhonov regularization for bone and soft tissue, respectively, and smooth total variation regularization for the gadolinium image.

Numerical study of the data fidelity convexity

We have shown in section 3 that the problem is convex in a certain region of the solution space (proposition 3) and that the zero mass density image belongs to this region (remark 4). This is a relevant result as it can serve as a guide to initialize convex optimization algorithms. Nevertheless, there may be regions of the solution space where the data fidelity term is nonconvex. In this section, we examine convexity numerically by computing the data fidelity term [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF] for two typical numerical phantoms. Figure 2 shows the data fidelity term for a one-pixel mass density image a while figure 3 displays the data fidelity term for the multiple-pixel thorax phantom. In the latter simulation, we consider only constant images for the three materials, i.e., the same value is assigned to the entire image for each material. This has been done for simplification but it is actually relevant as constant images are commonly used as initial guesses. These examples show that the zero mass density image belongs to the convex region and that the problem is nonconvex for large values of the three materials. It can be noticed that the data fidelity term becomes asymptotically flat for increasing values of the three materials, i.e., lim a→∞ L(a) = 1 2 s δ 2 W .

(a) (b) (c) , and a true gd = 0.1 g/cm 3 . Data fidelity is displayed versus the value of one material while the other two materials remain fixed and equal to (a) a bone = 0.15 g/cm 3 and a gd = 0.007 g/cm 3 , (b) a soft = 4 g/cm 3 and a gd = 0.007 g/cm 3 , and (c) a soft = 4 g/cm 3 and a bone = 0.15 g/cm 3 .

Figure 4 displays the two characteristic sets A and B that are introduced in proposition 2 and 3, respectively, for the one-pixel mass density phantom (same as figure 2). The colormap represents the number of energy bins 1 ≤ i ≤ I satisfying F i (a) ≥ s i . The solution is set to a true soft = 10 g/cm 3 , a true bone = 1 g/cm 3 , and a true gd = 0.1 g/cm 3 (see red circle). The data is computed using the parameters plotted in figure 1.

The set A is depicted in yellow (all bins satisfy F i (a) ≥ s i ) while the set B is indicated by the grey dashed line.

Convergence of GNB

Solution error and data fidelity versus the number of iterations are shown in figure 5 for GNB for different values of α. GNB leads to an optimal solution as long as α is sufficiently large, where the larger α the more iterations it requires [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF][START_REF] Abascal | Fluorescence diffuse optical tomography using the split Bregman method[END_REF]. For these results, we used κ = 10 -6 , initial guess a 0 = 0 and tolerance tol = 1.12 • 10 5 . As a comparison, we also provide the solution error for GN for the value of the regularization parameter (α = 0.23) that led to the lowest solution error. For GN we set initial guess to a 0 = 0 and tolerance to ψ = 10 -4 , which corresponded to a value of the data fidelity term of 1.12 • 10 5 . GN converged in N tot = 16 iterations for α = 0.23 (see figure 5(b)). GNB solves a sequence of subproblems, so it needs a larger number of iterations. It needed N tot = 28 and N tot = 40 for α = 2 and α = 10, respectively.

Figure 6(a) shows the solution error versus the regularization parameter α for GN and versus the parameter α for GNB. While GNB provides the same solution for all α's larger than 0.5, GN depends dramatically on the choice of α. Although GN converged in only N tot = 16 iterations (for α = 0.23, see figure 5(b)), extra decompositions for different regularization parameters are required. On the contrary, GNB can be used with a large α chosen a priori.

Figure 6(b) shows that the number of Bregman iterations is almost constant with α and equal to k = 15 in the cases investigated here. Figure 6(c) shows that the total number of Gauss Newton iterations N tot is linearly related to the parameter α. Figure 7 shows the number of Gauss Newton iterations N k needed to minimize Ψ k , for different values of α. An important remark is that GNB required nearly only one iteration N k = 1 to solve each subproblem for k > 1, so it adds just few iterations with respect to GN.

Figure 8 shows decomposed images by GN and GNB. Second row in figure 8 gives decomposed images by GN for an optimal regularization parameter α = 0.23. Third to fifth rows are the decomposed images by GNB for different number of iterations, setting α = 10, κ = 1 • 10 -6 and a 0 = 0. Third and fourth rows correspond to an intermediate and final solution of the first subproblem (k = 1). The solution of the first subproblem provides a good approximation of the decomposed materials but images are smooth given the large value of α. Further Bregman iterations (k = 25) are needed to find the optimal solution (bottom row in figure 8).

We evaluated two different strategies for the initialization of the inner problem in GNB (algorithm 1) (figure 9). Using the solution found in the previous subproblem as initial guess in the next subproblem, i.e., setting a k,n=0 = a k-1,N k-1 , led to an optimal solution in N tot = 120 iterations for α = 100. Initializing each Bregman iteration to a zero image, i.e., setting a k,n=0 = 0, led to larger solution error and diverged quickly after few updates of the Bregman distance.

Figure 10 shows the influence of the number of inner iterations in GNB. Previous results correspond to full convergence of the inner problem. Taking only one iteration for the inner problem (N 1 = 1) resulted in unstable subgradient update and did not converge. However, setting N k = 1 for k > 1 (after the first subproblem is solved to high accuracy) almost halved the number of iterations with respect to N k = 2. We remark that convergence of the inner iteration required N 1 = 15 for k = 1 and at most N k = 3 for k > 1.

Influence of the initial guess

In previous sections the initial guess was set to a zero constant image for all materials, which belongs to the convex part of the data fidelity term. To assess now the effect of the Bregman distance on global convergence, we set the initial guess to large values (a 0 = (10, 10, 10)) in order to make sure it belong to the nonconvex region. GNB converged to the optimal solution independently of the initial guess. On the contrary, GN did not converge, stopping at the first iteration as a descent direction could not be found (figure 11).

GNB presented a dependence on κ when the initial guess belonged to the nonconvex region. In this case, κ was critical as the Hessian was ill-conditioned, which occurred for a large value of the initial guess. We used κ = 10 -6 . When the initial guess belonged to the convex region, GNB led to same results independently of the value of κ as long as it was small enough κ ≤ 10 -4 . . We used the following parameters: α = 10, κ = 10 -6 , and initial guess a 0 = 0.

Influence of noise and phantom

In spectral CT, an important problem is the low photon count in each energy bin corresponding to noisy data. To illustrate the impact of noise on GNB, we consider first pure Poisson noise. Figure 12 shows decomposed material images by GNB for different number of incident photons from n 0 = 10 8 to n 0 = 10 5 . The lower the number of incident photons the more regularized the solution. The tolerance tol was set to 1.3 • 10 5 , 1.5 • 10 5 and 1.8 • 10 5 for n 0 = 10 7 , n 0 = 10 6 and n 0 = 10 5 , respectively.

Figure 13 shows material decomposed images by GNB for mixed Poisson noise model with n 0 = 10 7 and different values of σ ranging from 5 to 50. GNB is able to decompose materials even with large σ but for σ larger than 20, images are noisier, which impacts the recovered concentration of gadolinium.

We also show the decompositions obtained for a different projection angle (here 180 • ) in figure 14.

Discussion

We proposed and assessed a regularized iterative scheme based on the Bregman distance (GNB) for the solution of the material decomposition problem in spectral CT, which is a nonlinear inverse ill-posed problem and for which we consider the minimization of a regularized functional. First, we proved that the problem is convex in a region of the solution space and that the zero solution belongs to this region. Then, using numerical examples we showed that the problem is nonconvex for large values of the decomposed materials.

Using a numerical thorax phantom, we showed that zero initialization is a suitable initialization strategy for convex optimization methods such as a regularized Gauss-Newton algorithm method [START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF]. GN led to an optimal solution as long as the initial guess belonged to the convex region of the data fidelity term. On the contrary, GN led to a wrong solution when the initial guess is chosen far from the convex region. The measurement vector is simulated for n 0 = 10 8 photons and the decomposition is performed using the following parameters: α = 10, κ = 10 -6 , and initial guess a 0 = 0.

We then assessed a regularized iterative scheme based on the Bregman distance. The proposed method has some advantages with respect to more standard approaches. It is robust to the choice of the initial guess, as it led to global minimum even when the initial guess was in the nonconvex region. On the contrary, Gauss-Newton method failed in this case. This is a relevant feature of the Bregman distance, which makes each subproblem locally convex around the current estimate. This makes Bregman iterative algorithms more robust but at the price of requiring more iterations than an equivalent Gauss-Newton method. Nevertheless, we found that the number of iterations required was not very large: around ten inner iterations to solve the first subproblem (k = 1) but generally only one to three inner iterations for k > 1. In addition, standard approaches such as Gauss-Newton method requires optimal selection of the regularization parameter using the L-curve or similar method [START_REF] Per | The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems[END_REF][START_REF] Abascal | Comparison of methods for optimal choice of the regularization parameter for linear electrical impedance tomography of brain function[END_REF]. In this work, we selected the solution provided by GN method that led to the lowest solution error by assessing a large number of regularization parameters. As Bregman iterative methods do not require optimal selection of the regularization parameter, the extra iterations needed by these methods may not translate in much larger computational demand than standard regularization methods. In addition, we found a linear relation between the number of total iterations and the parameter α, which has been previously shown for linear measurement operators [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF][START_REF] Abascal | Fluorescence diffuse optical tomography using the split Bregman method[END_REF]. Hence, α can be tuned to decrease the number of required iterations. The smaller α the lower the number of iterations but a very low value may be less robust.

A key point in the construction of these iterative schemes is to select efficient methods for solving the subproblem (17a). We have proposed an efficient Gauss-Newton step and found that it required only few inner iterations for k > 1 (algorithm 1). This may indicate that subproblems are approximately linear when close to the solution.

Reconstructed images show that the first subproblem leads to a good estimate of the decomposed materials and then further Bregman iterations are used to improve the solution. This could be then further exploited to reduce computationally complexity. In fact, fixing the number of inner iterations to one for k > 1 (after the first subproblem is solved to high accuracy) converged to the optimal solution and almost halved the number of iterations with respect to taking two inner iterations.

In [START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF], authors proposed two iterated methods based approximating the subproblem (17a): a Levenberg-Marquardt-type (ILM) method that linearized the forward operator, and a Landweber-type method that linearized the least square data fidelity term. There are similarities between GNB and ILM with only one inner iteration. ILM replaces [START_REF] Juan | Incorporation of prior knowledge of the signal behavior into the reconstruction to accelerate the acquisition of MR diffusion data[END_REF] by the iterated scheme

a k+1 = arg min a 1 2 F(a k ) + F (a k )(a -a k ) -s δ 2 W + αD R ξ k (a, a k ), (22a)
ξ k+1 = ξ k -α -1 W F (a k ) W F(a k ) + F (a k )(a k+1 -a k ) -s δ . (22b)
Doing only one inner iteration on GNB (algorithm 1) and ILM [START_REF] Heismann | Spectral computed tomography[END_REF] provides the same solution update a k+1 but subgradient updates are different. In order to compare subgradient updates given by the two methods, we rewrite ILM's update as ξ k+1 ξ k -α -1 W F (a k ) W (F (a k+1 ) -s). While the Jacobian in ILM is evaluated on the previous estimate a k , as the forward operator has been linearized in the data fidelity term, the Jacobian in GNB is evaluated on the current estimate a k+1 . This difference seems to have a significant impact when approximating the inner problem, as discussed below.

An unexpected result in the analysis was to find that fixing the number of inner iterations to one in GNB resulted in an unstable subgradient update. We observed that after the first subgradient update, the Bregman distance becomes negative and then the algorithm diverges. As the Bregman distance is given by D

R ξ k (a, a k ) = R(a) -R(a k )-< ξ k , a -a k >≥ 0 for ξ k ∈ ∂R(a k
), then we assume that ξ k / ∈ ∂R(a k ). To further understand this issue, we computed GNB and ILM when taking one inner iteration for all k. While GNB did not converge, ILM converged to a solution close to the optimal one but results were slightly worse than those provided by GNB when solving exactly the first subproblem. Further work is needed for a rigurous comparison between GNB and ILM. From these results and the difference between GNB and ILM discussed above, it seems that given the nonlinearity, taking only few inner iterations for GNB may not be stable in general. However, if the first inner problem is solved to high accuracy, it is then feasible to take only one inner iteration. This may be due to the fact that solving the first subproblem brings the first iterate close to the optimal solution and the rest of subproblems are approximately linear.

Initialization of each Bregman iteration for GNB was critical. Using the solution found at the previous iteration led to robust convergence. On the other hand, a zero initialization strategy for each subproblem required a large number of iterations to solve each subproblem and became unstable as the number of Bregman iterations increased. This can be explained by the fact that the problem is nonlinear and that as the number of Bregman iterations increases, the solution is further away from the zero initial guess. Nevertheless, for linear problems, zero initialization was a valid strategy for Bregman iterative methods in a wide variety of applications [START_REF] Abascal | Fluorescence diffuse optical tomography using the split Bregman method[END_REF][START_REF] Abascal | A novel prior-and motion-based compressed sensing method for small-animal respiratory gated CT[END_REF].

In [START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF] authors included an extra 2 -term (weighted by κ in ( 18)) to make the functional associated with total variation convex on each subproblem. Here we found a dependence on κ. When the initial guess was in the nonconvex region, this term was essential to reach to a global minimum. An optimal solution was found by using a sufficiently small value of κ (κ ≤ 10 -4 ). The effect of κ can be explained by the fact that when the initial guess is in the nonconvex region the Hessian is ill-conditioned, where the nonconvex region in this problem corresponds to large values of the materials. For the case where the initial guess was in the convex region, this term did not play a role and any small value of κ (including κ = 0) led to the same solution.

In conclusion, we have studied convexity for the material decomposition problem and showed that an iterative scheme based on the Bregman distance had a positive effect to find a global minimum. In addition to this, the Bregman distance presents other benefits. For instance, it can efficiently minimize nonsmooth functional such as total variation or impose a nonnegativity constraint, which has been exploited for linear measurements [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF][START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF][START_REF] Abascal | Fluorescence diffuse optical tomography using the split Bregman method[END_REF]. Two variants of the proposed GNB algorithm, GNB with a nonnegativity constraint and GNB with a nondifferentiable convex functional, are included in the appendix B. Another benefit of using the Bregman distance is to correct for the lost of contrast that affects total variation [START_REF] Burger | A Guide to the TV Zoo[END_REF]. This is crucial for quantitative imaging techniques such as spectral CT. In future work we will investigate these other benefits of the Bregman distance for solving the material decomposition problem in spectral CT.

Variables a and b in (B.3a) are de-coupled, so they can be solved in two alternating steps: ) can be also done with GN as in algorithm 1. Subproblem (B.4b) can be solved analytically by using a thresholding operation [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF], [START_REF] Abascal | Fluorescence diffuse optical tomography using the split Bregman method[END_REF].

a k+1 =

Appendix B.2. GNB for nondifferentiable functional

The minimization of convex nondifferentiable functional such as total variation can be solved efficiently using the split Bregman method, which has been widely exploited for linear measurements [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF]. In specific, the split Bregman method is very efficient for solving L1-regularized problems. Let R(a) be a nondifferentiable convex functional of the form R( ) can be solved analytically by using a thresholding operation [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF], [START_REF] Abascal | Fluorescence diffuse optical tomography using the split Bregman method[END_REF].

Figure 1 .

 1 Figure 1. Forward model parameters used in the numerical simulations. (a) Source spectrum n 0 (E). (b) Detector response function d i (E) for each of the I = 4 energy bins. (c) Linear mass attenuation coefficients τ m (E) for each of the M = 3 materials.

Figure 2 .Figure 3 .

 23 Figure 2. One dimensional plot of the data fidelity term for a one-pixel image a ∈ R M .The ground truth solution is set to a true soft = 10 g/cm 3 , a true bone = 1 g/cm 3 , and a true gd = 0.1 g/cm 3 (see red circle). Data fidelity is displayed versus the value of one material while the other two materials remain fixed and equal to the ground truth solution.

Figure 4 .

 4 Figure 4. Characteristic sets included in the solution space. The colormap represents the number of energy bins 1 ≤ i ≤ I satisfying F i (a) ≥ s i . The solution is set to a true soft = 10 g/cm 3 , a true bone = 1 g/cm 3 , and a true gd = 0.1 g/cm 3 (see red circle). The data is computed using the parameters plotted in figure 1. The set A is depicted in yellow (all bins satisfy F i (a) ≥ s i ) while the set B is indicated by the grey dashed line. Sets are displayed versus the value of two materials while the other material remain fixed and equal to the ground truth solution.

Figure 5 .Figure 6 .

 56 Figure 5. Convergence of the algorithms. (a) Solution error versus the number of Gauss Newton iterations for GNB (full line) and for GN (dotted line). (b) Data fidelity decrease for GNB (α = 10) and for GN (α = 0.23). We used a 0 = 0 as initial guess.

Figure 7 .

 7 Figure 7. Number of Gauss Newton iterations, N k , needed at each Bregman iteration, k, for different values of α.

Figure 8 .

 8 Figure 8. (a) Reference image for the thorax phantom comprising three materials: soft tissue, bone and portal vein marked with gadolinium. A particular projection view at angle of 60 • is considered. (b) Material images decomposed by GN for α = 0.23 and a 0 = 0. (c)-(e) Material images decomposed by GNB at different iterations: (c) k = 1, n = 5, N tot = 5, (d) k = 1, n = 15, N tot = 15 and (e) k = 25, n = 1, N tot = 44, which is the optimal decomposition with GNB. k indicates the (outer) Bregman iteration, n the (inner) Gauss Newton iterations, and N tot is the total number of Gauss-Newton iterations. We set α = 10, κ = 1 • 10 -6 and a 0 = 0.

Figure 9 .Figure 10 .

 910 Figure 9. Influence of the initialization of the Bregman iteration. The solution error is plotted for two types of initialization for the minimization of Ψ k : using the solution found at the previous Bregman iteration (full line) and using a zero image for all Bregman iteration (dotted line). We used α = 100, κ = 1 • 10 -6 and a 0 = 0.

Figure 11 .

 11 Figure 11. Influence of the initial guess. (a) Reference image and decomposed images by using (b) GN with α = 0.23 and (c) GNB with α = 10 and κ = 10 -6 . The initial guess was a 0 = 10 for both methods.

Figure 12 .

 12 Figure 12. Material decomposed images by GNB for different number of photons n 0 . Pure Poisson noise is considered, i.e., σ = 0. (a) Reference images and decomposed images for (b) n 0 = 10 7 , (c) n 0 = 10 6 and (d) n 0 = 10 5 . We used the following parameters: α = 10, κ = 10 -6 , and initial guess a 0 = 0.

Figure 13 .

 13 Figure 13. Material decomposed images by GNB for mixed Poisson noise model, s δ ∼ P(λ = s) + N (µ, σ 2 ), for n 0 = 10 7 and concentration of gadolinium of 0.1 g•cm -3 . (a) Reference images and decomposed images for (b) σ = 5, (c) σ = 20, and (d) σ = 50. We used the following parameters: α = 10, κ = 10 -6 , tol = 1.23 • 10 5 and initial guess a 0 = 0.

Figure 14 .

 14 Figure 14. Material decomposed images provided for a different projection view (for an angle of 180 • ). (a) The ground truth. (b) Decomposed data provided by GNB.The measurement vector is simulated for n 0 = 10 8 photons and the decomposition is performed using the following parameters: α = 10, κ = 10 -6 , and initial guess a 0 = 0.

-a k+1 2 +

 2 α1 b≥0 -α < ξ k b , b -b k > . (B.4b) Subproblem (B.4a) is a slight variation of subproblem (17a) where a quadratic term is added to the previous cost function Ψ k (a) and where ξ a has an extra term. Hence, minimization of (B.4a

s δ 2 W 2 Fα d 1 - 2 F(a) -s δ 2 W

 22122 a) = Φ(a) 1 , ie. total variation R(a) = ∇a 1 , we replace Φ(a) by d and consider the following constrained problem: ≤ σ, d = Φ(a). (B.5)Proceeding as above, we define U (a, d) = d 1 , the constrained problem (B.5) can be solved by the following iterative scheme:(a k+1 , d k+1 ) = arg min a,d αD U ξ k (a, a k , d, d k ) + 1 < ξ k a , a -a k > -< ξ k d , d -d k > + 1 a -α -1 (W F (a k+1 )) (W F(a k+1 ) -W s δ ) -α -1 Φ (a) T (Φ(a k+1 ) -d k+1 ) (B.6b) ξ k+1 d = ξ k d -α -1 (d k+1 -Φ(a k+1 )). (B.6c)Variables a and d in (B.6a) are de-coupled, so they can be solved in two alternating steps: Φ(a) 2 -α < ξ k a , a -a k > (B.7a)d k+1 = arg min d α d 1 + 1 2 d -Φ(a k+1 ) 2 -α < ξ k d , d -d k > . (B.7b) Subproblem (B.7a) is easier to solve than subproblem (17a) as it does not depend on R(a) and a quadratic term is added to the previous cost function. Subproblem (B.7b

Acknowledgments

Authors would like to thank Bernhard Brendel, Ewald Roessl and Simon Rit for fruitful discussions and suggestions on the material decomposition problem. Authors also thank the anonymous reviewers and editorial board for their valuable comments and suggestions to improve the quality of the paper.

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement N • 701915. We also acknowledge the support of the ANR project SALTO (ANR-17-CE19-0011-01). It was also performed within the framework of the LabEx PRIMES (ANR-11-LABX-0063) of University de Lyon. This project has received funding from the EU's H2020 research and innovation program under the grant agreement N • 633937.

Appendix A. Hessian of the forward mapping

Let H i (a) be the Hessian matrix of the operator F i . It is defined as

The operator F i is continuously twice differentiable and its second derivative is given by

where g i (a, E) is given by [START_REF] Wu | Augmented Lagrangian Method, Dual Methods, and Split Bregman Iteration for ROF, Vectorial TV, and High Order Models[END_REF]. Discretising the energy interval E = [0, E max ] with an energy step ∆E, the integral (A.2) can be approximated by

Hence, it can be noted that the Hessian takes the following matrix form

where G i (a) = diag (g i (a, E 1 ), . . . , g i (a, E J )) ∆E is the vector of photons incident on the detector, and T = (τ j,m ) ∈ R J×M is the mass attenuation matrix defined by τ j,m = τ m (E j ).

Appendix B. Extensions of GNB

Appendix B.1. GNB with nonnegativity constraint

Including a nonnegativity constraint in GNB can be done by combining the Bregman distance with a splitting technique, which leads to the split Bregman method [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF], [START_REF] Abascal | Fluorescence diffuse optical tomography using the split Bregman method[END_REF].

Proceeding as in [START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF], let 1 a≥0 be the indicator function of the set R N ≥0 and add an auxiliary variable b to allow for the splitting, the material decomposition problem with a nonnegativity constraint can be posed as follows:

Defining the convex functional U (a, b) = R(a) + 1 b≥0 and the Bregman distance associated with this functional

the constrained problem (B.1) can be solved by the following iterative scheme:


Nicolas Vin

Frédérique Battin-Leclerc

Hervé Le Gall

Nadia Sebbar

Henning Bockhorn

Dimosthenis Trimis

Olivier Herbinet

A study of chlorobenzene pyrolysis

Keywords: PCBs, Pyrolysis, Chlorobenzene, Modeling, Hydrogen addition

The pyrolysis of chlorobenzene under dilute atmosphere and quasi-atmospheric pressure was studied at temperatures from 800 to 1150 K using a fused silica jet stirred reactor (JSR) and from 800 to 1250 K in an alumina tubular reactor. Chlorobenzene was chosen as a surrogate to model the thermal decomposition of polychlorinated biphenyls (PCBs). In the jet stirred reactor, a maximum chlorobenzene conversion of 48.5 % was observed at a residence time of 2 s, a temperature of 1150 K and an inlet mole fraction of chlorobenzene of 0.005. The following species were quantified: benzene (the major product), HCl, methane, ethylene, acetylene, biphenyl, 1-, 2and 3-chlorobiphenyl, biphenylene and six isomers of dichlorobiphenyls. In the tubular reactor, a maximum chlorobenzene conversion of 95% was observed at a temperature of 1250 K under the same conditions as in the jet stirred reactor. The same reaction products were detected but with a larger formation of acetylene and methane and a smaller production of chlorinated and bicyclic compounds. The effect of two H-atom donors, methane and hydrogen, has been investigated in the JSR. Hydrogen addition has a strong inhibiting effect on the formation of chlorinated and bicyclic products. A new detailed kinetic model was developed and gave a good prediction of the global reactivity, the formation of most major products, and the effect of the addition of hydrogen and methane. Flow rate and sensitivity analyses have been made to explain the effect of H-atom donors.

Introduction

Chlorinated compounds are widely used in the chemical industry. The industrial value of the polychlorobiphenyls (PCBs) resulted from their high thermal stability, chemical inertness and non-flammability coupled with excellent electrical properties [1]. The commercial production of PCBs started in 1929 but their use has been banned or severely restricted in many countries since the 70s because of possible risks to human health and environment [1].

Since 1929, around 2 million tons of PCBs have been produced, about 10% of which still remains in the environment today [1]. The sources of PCBs pollution are landfills containing transformers, capacitors and other PCB wastes [1]. While incineration is currently the most widespread disposal technique to treat wastes contaminated by PCBs, it can be a source of dioxins. It is why it is of interest to study an alternative process working in absence of oxygen, such as thermal decomposition. It is crucial then to understand the chemistry of the thermal decomposition of PCBs to predict the influence of operating conditions on the nature and concentration of reaction products. Chlorobenzene was chosen as a surrogate to model the thermal destruction of PCBs.

Because of the large range of temperatures which needs to be investigated, different experimental facilities were previously used to study the thermal decomposition of chlorobenzene: shock tubes [2][3][4][5] to study temperatures above 1500 K and plug flow reactors for lower temperatures [6][7][8][9][10]. In all those studies, benzene and HCl were the main products of the degradation of chlorobenzene. In terms of products with high molar mass, naphthalene, chloronaphthalene, acenaphthylene, biphenyl and all isomers of chlorobiphenyls and dichlorobiphenyls were found in small amounts [6][7][8][9][10]. Traces of cyclopentadiene and some substituted benzenes were detected [10]. Some of these studies were performed in presence of a hydrogen donor, generally H2, which enhances the conversion of C6H5Cl into benzene and HCl [2,10].

The objective of this study is to investigate chlorobenzene pyrolysis using two different types of reactors: a Jet Stirred Reactor (JSR) and a Tubular Reactor (TR) in order to better understand the chemistry of chlorinated polycyclic compounds. No previous study was performed about the pyrolysis of chlorobenzene in a JSR. A new detailed kinetic model has been developed to reproduce these experimental results.

Experimental procedure

Reactants entered the spherical fused silica jet-stirred-reactor (JSR, volume of 85 cm 3 ) though an injection cross made of four nozzles which was located at its center and created high turbulence to produce homogeneity in composition and temperature of the gas phase. The isothermal JSR was preceded by a quartz annular preheating zone in which the temperature of the gas is increased up to the reactor temperature. The gas residence time inside the annular preheater was very short compared to its residence time inside the reactor (a few percent). Both the reactor and the preheating zone were heated by the means of Thermocoax resistances rolled up around the reactor. The reaction temperature was measured by a thermocouple located inside the intraannular space of the preheating zone; its extremity being placed at the level of the injection cross. Experiments in the JSR were performed at a constant pressure of 106.7 kPa, at a residence time of 2 seconds and at temperatures ranging from 800 to 1150 K.

Experimental studies in the tubular reactor (TR) have been performed using a horizontal 60 cm long alumina tube (volume of 294 cm 3 ) with an inner diameter of 20 mm and an outer diameter of 25 mm. The ends of the reactor was connected to 2.54 cm (1-in.) diameter Swagelok unions and sealed with high temperature epoxy resin from Final Crotonics. The tubular reactor was heated by an electrical furnace from Vecstar equipped with an S-type thermocouple. For each set reaction temperature, a temperature profile (temperature versus position in the tubular reactor) was measured using a S-type thermocouple (see measured profiles in Figure S1 in Supplementary materials (SM)). Experiments were performed at a constant pressure of 106.7 kPa, at a residence time in the set-point temperature zone around 2 seconds and at temperatures ranging from 800 to 1300 K. Simulations were performed assuming the TR can be modeled as a plug flow reactor (see discussion in SM) and using the measured temperature profiles.

As shown in Table 1, three reaction mixtures were studied: neat chlorobenzene, addition of hydrogen, and addition of methane. Helium, methane and hydrogen were provided by Messer (purity of 99.999%) and the chlorinated reactant was purchased from Sigma-Aldrich (purity of 99.5%). Gas flow rates were controlled by mass flow controllers and the liquid flow rate by a Coriolis flow controller. The uncertainty in the flow measurements was around 0.5% for each controller, so about 1% on the residence time.

Using a heated transfer line maintained at 433 K to avoid product condensation, the gases leaving the reactors were analyzed using:

• Gas chromatography (GC): a first chromatograph, fitted with a PlotQ capillary column, a thermal conductivity detector (TCD) and a flame ionization detector (FID), was used for the quantification of methane, ethylene, acetylene and ethane; a second chromatograph, fitted with a HP-1 capillary column and a FID, was used for the quantification of other types of molecules. The identification of reaction products was performed using a gas chromatograph equipped with a HP-5MS capillary column and coupled to a mass spectrometer. Response factors were determined by injecting calibration mixtures or using the effective carbon number method. Relative uncertainties in mole fractions were estimated to be ± 5% for species which were calibrated using standards and ± 10% for other ones.

• A Fourier Transform InfraRed spectrometer (FTIR) from Thermo Scientific Antaris equipped with a Mercure Cadmium Telluride photoelectric detector. FTIR was used to analyze chlorobenzene and HCl. FTIR calibrations were obtained by injecting standards. A typical FTIR spectrum obtained during chlorobenzene pyrolysis in given in SM. Note that, despite high dilution in helium, the experimental results presented hereafter were challenging to obtain. The pyrolysis of chlorobenzene leads at high temperature to the formation of high molecular weight polycyclic aromatic compounds which can condensate in our experimental setup or often plug the injection valves inside our gas chromatographs, leading to large overpressure inside our experimental setup (especially in JSR experiments).

Kinetic model

Simulations were performed using the OpenSmoke package [11] and a newly developed detailed chemical kinetic reaction mechanism for chlorobenzene pyrolysis including 224 species in 1455 reactions (given in SM under Chemkin format). Rate parameters were taken from the literature when available. Calculations using quantum methods were carried out when no data were available (See Table 2). Kinetics parameters of the reaction C6H6 + H = C6H5 + H2 (the reverse of Reaction 1 in Table 2) have already been calculated by several authors such like Zhang et al. [12] or Giri et al. [13].

The model is based on that proposed by Husson et al. [START_REF] Husson | Proc. Combus. Inst[END_REF] for the oxidation of ethyl-benzene. The rate constant of reaction 1 in Table 2, already presents in the model of Husson et al. [START_REF] Husson | Proc. Combus. Inst[END_REF], has been recalculated.

The reactions of the H/Cl system involving HCl, Cl2, and Cl radicals were taken from Pelucchi et al. [START_REF] Pelucchi | [END_REF].

The following types of reactions were added to account for the reactant consumption:

• Unimolecular initiations by breaking a C-H or a C-Cl bond; the C-Cl bond (95 kcal/mol [START_REF] Luo | Comprehensive Handbook of Chemical Bond Energies[END_REF]) is the weakest one in the chlorobenzene molecule (C6H5Cl, PHICL in Table 2).

• H-atom abstractions by H-atoms (reactions 2-4 in Table 2), by Cl-atoms and by phenyl radicals. H-atom abstractions lead to phenyl and chlorophenyl radicals (C6H4Cl, RxPHICL (x = 2, 3, 4)). The rate constants of the abstractions by Cl-atom and phenyl radical were taken equal to those with H-atoms.

• Cl-atom abstractions by H-atoms (reaction 5),

• Ipso-additions of H-and Cl-atoms (reactions 6-9) leading to the formation of benzene (C6H6) and dichlorobenzene (C6H4Cl2, PHI1xCL2 (x = 2, 3, 4)), even if no trace of dichlorobenzene was identified. Ipso-additions are important due to the aromatic reactant structure.

• A molecular reaction (reaction 10) to better account for the reactant consumption in the presence of hydrogen. In our case, adding this molecular reaction had no effect on modeling results.

• Ipso-additions of phenyl radical leading to ortho, meta and para-chlorobiphenyl (o/m/p C12H9Cl) or biphenyl:

C6H5Cl+C6H5=o/m/pC12H9Cl+H C6H5Cl+C6H5=C12H10+Cl

• Ipso-addition of C6H4Cl radicals giving ortho, meta and para-chlorobiphenyl or dichlorobiphenyls (C12H8Cl2):

C6H5Cl+C6H4Cl=o/m/pC12H9Cl+Cl C6H5Cl+C6H4Cl=C12H8Cl2+H

• Terminations yielding chlorobiphenyl and dichlorobiphenyl isomers.

Due to the large number of involved heavy atoms, the rate constants of the three last types of reactions have been estimated. For ipso-additions, the rate parameters have been taken equal to those of ipso-addition of phenyl radical on benzene [START_REF] Wang | [END_REF] taking into account the number of abstractable atoms. The rate parameters of terminations have been taken equal to those of the termination of two phenyl radicals yielding biphenyl [18].

To model results in presence of methane, the ipso-addition of methyl radical to chlorobenzene was added:

C6H5Cl+CH3=C7H8+Cl

Its rate constant was taken equal to that of the ipso-addition of methyl radical to benzene.

For the calculation of the molecular properties of the species, the Gaussian 09 program [START_REF] Frisch | Gaussian 09, Revision A02[END_REF] was used. Density Functional theory (DFT) at B3LYP/6-311-G(d,p) level was applied for the determination of molecular properties of reactants, intermediates and products as well as transition state structures involved in the system and not available in the literature. DFT may be one of the few applicable calculation methods for large molecule systems at the lowest computational costs and previous studies [START_REF] Sebbar | [END_REF]21] have shown good agreement with ab initio methods. Vibration, translation, and external rotation frequencies, moments of inertia, symmetry, spin degeneracy, and optical isomers were used for the contributions to S and Cp. Torsion frequencies were used for hindered internal rotors when present. Frequencies were identified with GaussView code.

The computer code, SMCPS [22], was used to calculate entropies S°298 and heat capacities Cp f 298(T) using the harmonic-oscillator approximation for vibrations, on the basis of frequencies and moments of inertia of the optimized structures. If present, torsional frequencies were not included in the contributions to entropy and heat capacities; instead, they were replaced with values from a separate analysis on each internal rotor analysis. Calculated enthalpies, entropies, heat capacities as well as optimized geometries for species and transition state structures are reported in SM. Other missing thermochemical data were taken from the database of Goos et al. [START_REF] Goos | Third Millenium Ideal Gas and Condensed Phase Thermochemical Database for Combustion with Updates from Active Thermochemical Tables[END_REF].

Results and discussion

This part presents the main results obtained during the pyrolysis of neat chlorobenzene in both reactors. The effect of two hydrogen donors, H2 and CH4, is also investigated in JSR.

Chlorobenzene conversion

Figure 1 presents experimental and simulated evolution with temperature of conversion obtained in JSR and TR under the same operating conditions. A fused silica JSR can reach a maximum possible temperature of about 1200 K. This was not high enough to fully destroy chlorobenzene. That is why the chlorobenzene pyrolysis was also studied in the alumina TR which can be used at temperatures above 1200 K.

Figure 1 shows that the type of reactor has almost no influence on chlorobenzene conversion. The molecule starts to be consumed at 900 K in both reactors. The maximum conversion in JSR (48.4%) was obtained at 1150 K, the highest studied temperature. Chlorobenzene was nearly fully destroyed in the TR (95% of conversion) at 1250 K. This was the maximum temperature where measurements were possible with chlorobenzene/helium mixtures due to a too large pressure increase at higher temperatures preventing performing reliable experiments despite the large helium dilution. The model predicts well chlorobenzene conversion in both reactors.

Under approximately the same operating conditions using also a TR, Ritter and Bozzelli [10] obtained up to 60% of decomposition of chlorobenzene at 1228 K (residence time around 2 s) in helium. The yield of chlorobenzene decomposition at 1273 K (residence time of 1.2 seconds) is about 95% according to Rouzet et al. [9]. These results are in agreement with the present measurements.

Product formation

An important fraction of chlorobenzene was converted to a carbonaceous solid deposit on the two reactors walls. This deposit has led to many analytical problems and leads to an increasing uncertainty on experimental results with increasing temperature.

It was possible to partially remove this carbonaceous solid deposit by applying an O2 stream to the reactor at high temperature (1200 K). Experiments of Rouzet et al. [9] revealed that the yield of chlorobenzene to soot was not important below 1123 K, but raised up to 50% at 1443 K. Ritter and Bozzelli [10] working with a TR found a soot yield of 30% at 1283 K in the absence of an hydrogen donor. Several more species have been observed in very small amounts in GC and could then not be identified.

Major products

Figure 1 also presents experimental and simulated mole fraction of benzene, dichlorobiphenyls and HCl. Between 800 K and 1130 K, the influence of the reactor type on the formation of benzene and of dichlorobiphenyl isomers is limited. From approximately 1130 K, the model predicts higher formation of benzene and dichlorobiphenyls in the TR than in the JSR. Benzene and HCl mole fractions are well predicted until 1130 K, but significantly underpredicted at higher temperatures. Dichlorobiphenyl mole fractions are well predicted in JSR, but the significant decrease in TR is not modeled. The important formation of heavy chlorinated aromatic products and soot could partly explain deviations between modeling and experiments. Above 1130 K, a deviation higher than 25% for the C-atom material balances and 19% for the Cl-atom material balances was observed. As shown by the deposit on the reactor wall, this came from the production of not identified heavy chlorinated compounds.

Minor products

To better compare the minor products distribution obtained in both reactors, a selectivity diagram has been plotted in Figure 2 for a reaction temperature of 1125 K. The selectivity is calculated as below: 𝑆𝑒𝑙𝑒𝑐𝑡𝑖𝑣𝑖𝑡𝑦 = ∑ , with 𝑥 the mole fraction of product I during the pyrolysis of chlorobenzene.

Figure 2 shows that benzene, acetylene and methane are more abundant in TR than in JSR. However, the mole fraction of polychlorobiphenyls is notably reduced in TR. We assume that dichlorobiphenyls and chlorobiphenyls are secondary products and so their formation is limited in a tubular reactor. The biphenyl and biphenylene profiles are displayed in SM (Figure S3a), as well as methane, ethylene and acetylene profiles in JSR (Figure S3b). In JSR, the model predicts well the formation of biphenyl (Figure S3a). The formation of biphenylene is not predicted by our model. Acetylene was not observed in JSR but very small amounts of methane and ethylene were identified (respectively 14 and 5 ppm at 1150 K).

In TR, No traces of biphenyl and biphenylene have been identified. Results of the modeling about the formation of methane, acetylene and ethylene is not shown in Figure S3b because their mole fraction is less than 0.6 ppm at all temperatures.

Effect of two hydrogen donors: H2 and methane (JSR)

Cullis and Manton [START_REF] Cullis | [END_REF] found that the addition of hydrogen accelerated the decomposition of chlorobenzene. Louw et al. [6] postulated that using H2 in excess would lead to a dechlorination of chlorinated aromatic compounds. Ritter and Bozzelli [25] obtained up to 99% of decomposition of chlorobenzene at 1228 K in a helium/hydrogen mixture.

Figure 3 displays the experimental and simulated JSR evolution with temperature of chlorobenzene conversion and mole fractions of benzene, toluene and dichlorobiphenyls with neat chlorobenzene and with additions of hydrogen or methane. Figure 3 shows that the experimental effect of hydrogen addition on reactant conversion is limited. However, hydrogen addition increases benzene formation and inhibits significantly the production of chlorinated polyaromatic compounds. The effect of the presence of hydrogen is well reproduced by the model. : JSR evolution with temperature of (a) chlorobenzene conversion and mole fractions of (b) benzene, (c) toluene and (d) dichlorobiphenyls with no additive (0.5% chlorobenzene inlet mole fraction), with addition of hydrogen (0.5% chlorobenzene inlet mole fraction, 0.5% hydrogen inlet mole fraction), and with addition of methane (0.5% chlorobenzene inlet mole fraction, 0.5% methane inlet mole fraction), at τ = 2 s and P = 106.7 kPa. Symbols are experiments (open circles for neat chlorobenzene, black squares for H2 addition, blue stars for CH4 addition) and lines simulations (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

The methane addition notably reduces chlorobenzene conversion, especially above 1000 K. This is due the formation of toluene, a well-known radical scavenger via the formation of resonancestabilized benzyl radicals. The formation of toluene is well predicted as shown by Figure 3c.

Discussion

As the model reproduces satisfactorily the JSR experimental results, Figure 4 presents a flow rate analysis performed at 1100 K under the conditions of Figure 3. This figure shows that, for neat chlorobenzene, the major way (33.3%) of reactant consumption is by H-abstractions yielding chlorophenyl radicals. Chlorophenyl radicals react then by ipso-addition on chlorobenzene giving isomers of dichlorobiphenyl. This is the second most important way (30.5%) of reactant consumption. The third most important way (19.0%) of reactant consumption is by Clabstractions producing phenyl radicals, a source of biphenyl. Minor consumption pathways (less than 10% each) include ipso-addition of H-atoms yielding benzene, ipso-addition of phenyl radicals abstracting a Cl-atom, the second source of biphenyl, and ipso-addition of phenyl radicals releasing a H-atom and producing isomers of chlorobiphenyl. Note that chlorophenyl radicals are the main source of dichlorobiphenyls. Ipso-addition of chlorophenyl radicals to benzene can also produce chlorobiphenyls.

Figure 4: Flow rate analysis at 1100 K under the conditions of Figure 3. Numbers on arrows (black bold for neat chlorobenzene, black italics for H2 addition, blue for CH4 addition) represents the consumption rate normalized by the total chlorobenzene consumption rate. Only one isomers of chlorophenyl radical is shown for clarity; for chlorobiaromatic compounds, the chorine atom can be carried by every atoms of the cycle. Red thick arrows correspond to the main ways of formations of chlorobiaromatic compounds (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

With the addition of hydrogen, the formation of phenyl radicals and benzene is favored, while the production of dichlorobiphenyl isomers is almost completely inhibited (less of 6% of reactant consumption). In addition, chlorophenyl radicals can easily react with hydrogen to give back chlorobenzene. Consequently, in presence of hydrogen, its net rate of formation is significantly reduced. This explains why the formation of dichlorobiphenyls is notably lower as it was clearly observed in JSR experiments (by a factor of around 27 at 1100 K) and in simulations, as shown in Figure 3. Since the formation of benzene and phenyl radicals is favored, the impact of hydrogen addition on the formation of chlorobiphenyls is more limited. Compared to neat chlorobenzene, the only difference in presence of methane is the occurrence of a new minor pathway: the production of toluene by ipso-addition of methyl radicals. Toluene is consumed by H-abstractions with H-atoms to give benzyl radicals, which can react with methane to give back toluene and with methyl radicals to produce ethylbenzene.

To better understand the most influential reactions, Figure 5 presents a sensitivity analysis performed at 1100 K under the conditions of Figure 3. With neat chlorobenzene, the reactions with the largest promoting effect are the Cl-abstraction from chlorobenzene by H-atoms producing phenyl radicals, the ipso-addition of H-atom or phenyl radicals to chlorobenzene yielding benzene or biphenyl, respectively, and the unimolecular decomposition of the reactant to give phenyl radicals and chlorine atoms. The reactions of chlorophenyl radicals have also some effect. The ipso-additions to the reactant show a promoting influence, while the termination between two chlorophenyl radicals has an inhibiting effect.

Figure 5: Sensitivity analysis on chlorobenzene mole fraction at 1100 K under the conditions of Figure 3. Reactions with a negative coefficient have a promoting effect on the reactivity.

For the H2 addition case, the promoting effect of the four reactions displaying the largest negative sensitivity coefficient in neat chlorobenzene is still enhanced. Note also a significant increase of inhibiting effect of the combination of phenyl radicals with H-atoms. The fact that the sensitivity coefficient of the chlorophenyl radical reactions are close to zero confirms the significant reduction of their formation pathways, as it was discussed previously together with Figure 4.

For the CH4 addition case, the only of the four reactions displaying the largest negative sensitivity coefficient in neat chlorobenzene keeping some influence, even if reduced, is the unimolecular decomposition. This shows that the influence of all the radical reactions is significantly reduced. Additional reactions showing some influence are the ipso-addition of methyl radicals to chlorobenzene producing toluene and the H-abstraction from methane by H-atoms.

Conclusion

A study of chlorobenzene pyrolysis has been carried out in two different reactors: a jet-stirred reactor and a tubular reactor. Temperatures were ranging from 800 to 1125 K in the JSR and from 800 to 1250 K in the TR. Residence time was in both reactors equal to 2 s with a chlorobenzene inlet mole fraction of 0.005. All experiments were performed at a pressure of 106.7 kPa. The following products were identified: HCl, H2, methane, ethylene, acetylene, ethane, benzene, biphenylene, chlorobiphenyls and dichlorobiphenyls with all of their isomers. Products written in bold are major products. An important formation of a carbonaceous solid deposit was observed in both reactors involving an increasing uncertainty on experimental results when temperature rises.

Tests for the addition of hydrogen donors (H2 and CH4) were made in order to limit the formation of biaromatic compounds and of carbonaceous solid. The addition of hydrogen in JSR leads to an increase of the benzene mole fraction and inhibits significantly the formation of dichlorobiphenyls. The addition of methane also inhibits the production of undesired chlorinated aromatic compounds, but it also decreases the reactant conversion due to the formation of toluene. The use of a process based on thermal decomposition to destroy PCBs can then be envisaged as an alternative solution to incineration only with a planned significant addition of hydrogen.

Since there is a lack of kinetic data in the literature about chlorinated compounds, a first approach detailed kinetic model has been developed based on theoretically calculated rate parameters. According to uncertainties in experimental results, this model gives satisfactory predictions of the reactivity and formation of major compounds, as well as the effect of the addition of hydrogen donors. However, improvements in the model are certainly needed to refine the prediction of minor products.

Supplementary Material

• Temperature profiles measured in the tubular reactor.

• Discussion about the plug flow reactor assumption. • A typical IRTF spectrum obtained during the chlorobenzene pyrolysis.

• Calculated enthalpies, entropies, heat capacities as well as optimized geometries for species and transition state structures.

• Evolution with temperature of (a) biphenyl and biphenylene in JSR, (b) methane, ethylene, acetylene in TR, (c) methane and ethylene in JSR mole fractions (Symbols are experiments and lines simulations, 0.5% chlorobenzene inlet mole fraction, τ = 2 s, P = 106.7 kPa).

Figure 1 :

 1 Figure 1: Evolution with temperature of (a) chlorobenzene conversion, (b) benzene, (c) HCl and (d) dichlorobiphenyls mole fractions in tubular and jet-stirred reactor (Symbols are experiments (black points and red open squares are GC measurements, blues cross are FTIR ones) and lines simulations, 0.5% chlorobenzene inlet mole fraction, τ = 2 s, P = 106.7 kPa) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Figure 2 :

 2 Figure 2: Carbon containing product selectivity obtained during chlorobenzene pyrolysis in both reactors (1125 K, 0.5% chlorobenzene inlet mole fraction, τ = 2 s, P = 106.7 kPa).

Figure 3

 3 Figure3: JSR evolution with temperature of (a) chlorobenzene conversion and mole fractions of (b) benzene, (c) toluene and (d) dichlorobiphenyls with no additive (0.5% chlorobenzene inlet mole fraction), with addition of hydrogen (0.5% chlorobenzene inlet mole fraction, 0.5% hydrogen inlet mole fraction), and with addition of methane (0.5% chlorobenzene inlet mole fraction, 0.5% methane inlet mole fraction), at τ = 2 s and P = 106.7 kPa. Symbols are experiments (open circles for neat chlorobenzene, black squares for H2 addition, blue stars for CH4 addition) and lines simulations (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Table 1 :

 1 Experimental inlet composition.

	Mixture	Neat chlorobenzene	+ H2	+CH4
	Composition (%)	0.5 C6H5Cl + 99.5 He	0.5 C6H5Cl +0.5 H2 + 99.0 He	0.5 C6H5Cl +0.5 CH4 + 99.0 He
	Used reactors	JSR and TR	JSR	JSR

Table 2 :

 2 Theoretically calculated rate constants.

		Reaction	A factor (cm 3. mol -1 .s -1 )	n	Ea (cal.mol -1 )
	1	C6H5 +H2 = C6H6 + H	5.02x10 11	-0.14	2749
	2	PHICL+H=R2PHICL+H2	1.50x10 8	1.79	12328
	3	PHICL+H=R3PHICL+H2	1.23x10 8	1.81	10930
	4	PHICL+H=R4PHICL+H2	1.35x10 8	1.81	11748
	5	PHICL+H=C6H5+HCL	1.61x10 9	1.57	10988
	6	PHICL+H=C6H6+CL	4.26x10 7	1.772	7951
	7	PHICL+CL=PHI12CL2+H	1.17x10 4	2.992	22576
	8	PHICL+CL=PHI13CL2+H	8.57x10 3	3.057	21927
	9	PHICL+CL=PHI14CL2 +H	5.48x10 3	3.027	21829
	10	PHICL+H2=C6H6+HCL	1.36x10 6	2.18	73174

Acknowledgments

This work has been supported by TERBIS, 943 rue Pasteur, 60700 Pont Sainte Maxence, France (CNRS agreement n°137730).