Multiplicative Noise Removal Using L1 Fidelity on Frame Coefficients

We address the denoising of images contaminated with multiplicative noise, e.g. speckle noise. Classical ways to solve such problems are filtering, statistical (Bayesian) methods, variational methods, and methods that convert the multiplicative noise into additive noise (using a logarithmic function), shrinkage of the coefficients of the log-image data in a wavelet basis or in a frame, and transform back the result using an exponential function.

We propose a method composed of several stages: we use the log-image data and apply a reasonable under-optimal hard-thresholding on its curvelet transform; then we apply a variational method where we minimize a specialized criterion composed of an ℓ 1 data-fitting to the thresholded coefficients and a Total Variation regularization (TV) term in the image domain; the restored image is an exponential of the obtained minimizer, weighted in a way that the mean of the original image is preserved. Our restored images combine the advantages of shrinkage and variational methods and avoid their main drawbacks. For the minimization stage, we propose a properly adapted fast minimization scheme based on Douglas-Rachford splitting. The existence of a minimizer of our specialized criterion being proven, we demonstrate the convergence of the minimization scheme. The obtained numerical results outperform the main alternative methods.

Introduction

In various active imaging systems, such as synthetic aperture radar, laser or ultrasound imaging, the data representing the underlying (unknown image) S 0 : Ω → I R + , Ω ⊂ I R 2 , are corrupted with multiplicative noise. It is well known that such a noise severely degrades the image (see Fig. 2(a)). In order to increase the chance of restoring a cleaner image, several independent measurements for the same image are realized, thus yielding a set of data:

S k = S 0 η k + n k , ∀k ∈ {1, • • • , K}, (1)
where η k : Ω → I R + , and n k represent the multiplicative and the additive noise relevant to each measurement k. Usually, n k is white Gaussian noise. A commonly used and realistic model for the distribution of η k is the one-sided exponential distribution:

η k : pdf(η k ) = µ e -µη k 1l I R + (η k );
the latter is plotted in Fig. 1(a). Let us remind that 1/µ is both the mean and the standard deviation of this distribution. The usual practice is to take an average of the set of all measurements-such an image can be seen in (see Fig. 2(b)). Noticing that 1 K K k=1 n k ≈ 0, the data production model reads

S = 1 K K k=1 S k = S 0 1 K K k=1 η k = S 0 η, (2)
see e.g. [START_REF] Achim | Novel bayesian multiscale method for speckle removal in medical ultrasound images[END_REF][START_REF] Ulaby | Handbook of Radar Scattering Statistics for Terrain[END_REF][START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF] and many other references. A reasonable assumption is that all η k are independent and share the same mean µ. Then the resultant mean of the multiplicative noise η in ( 2) is known to follow a Gamma distribution,

η = 1 K K k=1 η k : pdf(η) = K µ K η K-1 Γ(K) exp - Kη µ , (3)
where Γ is the usual Gamma-function and since K is integer, Γ(K) = (K -1)!. Its mean is again µ and its standard deviation is µ/K. It is shown in Fig. 1(b).

Various adaptive filters for the restoration of images contaminated with multiplicative noise have been proposed in the past, e.g. see [START_REF] Krissian | Oriented speckle reducing anisotropic diffusion[END_REF][START_REF] Yu | Speckle reducing anisotropic diffusion[END_REF] and the numerous references therein. It can already been seen that filtering methods work well basically when the noise is moderate or weak, i.e. when K is large. Bayesian or variational methods have been proposed as well; one can consult for instance [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF][START_REF] Huang | A new total variation method for multiplicative noise removal[END_REF][START_REF] Rudin | Multiplicative denoising and deblurring: Theory and algorithms[END_REF][START_REF] Walessa | Model-based despeckling and information extraction from sar images[END_REF] and the references cited therein.

A large variety of methods-see e.g. [START_REF] Achim | Sar image filtering based on the heavy-tailed rayleigh model[END_REF][START_REF] Fukuda | Suppression of speckle in synthetic aperture radar images using wavelet[END_REF], more references are given in § 1.1-rely on the conversion of the multiplicative noise into additive noise using v = log S = log S 0 + log η = u 0 + n.

In this case the probability density function of n reads (see Fig. 1(c)):

n = log η : pdf(n) = K µ K 1 Γ(K) exp K(n -µe n ) . (5)
One can prove that

E [n] = ψ 0 (K) -log K , (6)
Var [n] = ψ 1 (K), (7)
where

ψ k (z) = d dz k+1 log Γ(z) (8)
is the polygamma function [START_REF] Abramowitz | Handbook of mathematical functions[END_REF].

Classical SAR modeling-see [START_REF] Tur | When is speckle noise multiplicative?[END_REF][START_REF] Ulaby | Handbook of Radar Scattering Statistics for Terrain[END_REF] and many other references-correspond to µ = 1 in [START_REF] Acar | Analysis of bounded variation penalty methods for ill-posed problems[END_REF]. Then

(3) and ( 5) boil down to pdf(η) = K K η K-1 e -Kη (K -1)! ,

pdf(n) = K K e K(n-e n ) (K -1)! . (9)

Multiscale shrinkage for the log-data

Many authors-see [START_REF] Achim | Novel bayesian multiscale method for speckle removal in medical ultrasound images[END_REF][START_REF] Achim | Sar image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF][START_REF] Fukuda | Suppression of speckle in synthetic aperture radar images using wavelet[END_REF][START_REF] Pizurica | A review of wavelet denoising in mri and ultrasound brain imaging[END_REF][START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF] and the references given there-focus on restoring the log-data as given in [START_REF] Achim | Novel bayesian multiscale method for speckle removal in medical ultrasound images[END_REF]. The common strategy is to decompose the log-data into some multiscale frame for L 2 (I R 2 ), say

{ w i : i ∈ I}: y = W v = W u 0 + W n, ( 10
)
where W is the corresponding frame analysis operator, i.e. (W v)[i] = v, w i , ∀i ∈ I. The rationale is that the noise W n in y is nearly Gaussian-as seen in Fig. 1

(d)-and justified by the Central Limit

Theorem. The obtained coefficients y have been considered in different frameworks in the literature. In a general way, coefficients are restored using shrinkage estimators using a symmetric function T : I R → I R, thus yielding

y T [i] = T (W v)[i] , ∀i ∈ I. (11)
Following [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF], various shrinkage estimators T have been explored in the literature, [START_REF] Antoniadis | Regularization of wavelet approximations[END_REF][START_REF] Murat Belge | Wavelet domain image restoration with adaptive edge-preserving regularization[END_REF][START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF][START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF][START_REF] Simoncelli | Noise removal via Bayesian wavelet coding[END_REF][START_REF] Wang | Solution of inverse problems in image processing by wavelet expansion[END_REF]; see § 2.1 for more details on shrinkage methods. Shrinkage functions specially designed for multiplicative noise were proposed e.g. in [START_REF] Achim | Novel bayesian multiscale method for speckle removal in medical ultrasound images[END_REF][START_REF] Achim | Sar image denoising via bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF][START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF].

Let W be a left inverse of W , giving rise to the dual frame { w i : i ∈ I}. Then a denoised log-image v T is generated by expanding the shrunk coefficients y T in the dual frame:

v T = i∈I T ((W v)[i]) w i = i∈I T (y[i]) w i . (12)
Then the sought-after image is of the form S T = exp v T .

Our approach and organization of the paper

We first apply (4) and then consider a tight-frame transform of the log-data. Our method to restore the log-image is presented in section 2. It is based on the minimization of a criterion composed of an ℓ 1 -fitting to the (suboptimally) hard-thresholded frame coefficients and a Total Variation (TV) regularization in the image domain. This method uses some ideas from a previous work of some of the authors [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF]. The minimization scheme to compute the log-restored image, explained in section 3, uses a Douglas-Rachford splitting specially adapted to our criterion. Restoring the sought-after image from the restored logimage requires a bias correction which is presented in section 4. The resultant algorithm to remove the multiplicative noise is provided in section 5. Various experiments are presented in section 6. Concluding remarks are given in section 7.

Restoration of the log-image

In this section we consider how to restore a good log-image given data v : Ω → I R obtained according to [START_REF] Achim | Novel bayesian multiscale method for speckle removal in medical ultrasound images[END_REF]. We focus basically on methods which, for a given preprocessed data set, lead to convex optimization problems. Below we comment only variational methods and shrinkage estimators since they underly the method proposed in this paper.

Drawbacks of shrinkage restoration and variational methods

Shrinkage restoration. The major problems with shrinkage denoising methods, as sketched in ( 11)- [START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF], is that shrinking large coefficients entails an erosion of the spiky image features, while shrinking small coefficients towards zero yields Gibbs-like oscillations in the vicinity of edges and a loss of texture information. On the other hand, if shrinkage is not sufficiently strong, some coefficients bearing mainly noise will remain almost unchanged-we call such coefficients outliers-and (12) suggests they generate artifacts with the shape of the functions w i of the frame. A well instructive illustration can be seen in Fig. 2(b-h). Several improvements, such as translation invariant thresholding [START_REF] Coifman | Translation-invariant de-noising[END_REF] and block thresholding [START_REF] Chesneau | Stein block thresholding for image denoising[END_REF], were brought to shrinkage methods in order to alleviate these artifacts. Results obtained using the latter method are presented in Figs. 3(c), 4(d) and 5(d) in Section 6. Another inherent difficulty comes from the fact that coefficients between different scales are not independent, as usually assumed, see e.g. [START_REF] Antoniadis | Wavelet thresholding for some classes of non-gaussian noise[END_REF][START_REF] Murat Belge | Wavelet domain image restoration with adaptive edge-preserving regularization[END_REF][START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF][START_REF] Simoncelli | Bayesian denoising of visual images in the wavelet domain[END_REF].

Variational methods. In variational methods, the restored function is defined as the minimizer of a criterion F v which balances trade-off between closeness to data and regularity constraints,

F v (u) = ρ Ω ψ u(t), v(t) dt + Ω ϕ(|∇u(t)|) dt, (13)
where ψ : I R + → I R + helps to measure closeness to data, ∇ stands for gradient (possibly in a distributional sense), ϕ : I R + → I R + is called a potential function and ρ > 0 is a parameter. A classical choice for ψ

is ψ(u(t), v(t)) = u(t) -v(t)
2 which assumes that the noise n in (4) is white, Gaussian and centered.

Given the actual distribution of the noise in [START_REF] Aubert | Modeling very oscillating signals. Application to image processing[END_REF] and Fig. 1(c), this may seem hazardous; we reconsider this choice in [START_REF] Candès | Fast discrete curvelet transforms[END_REF]. A reasonable choice is to use the log-likelihood of n according to [START_REF] Aubert | Modeling very oscillating signals. Application to image processing[END_REF] and this was involved in the criterion proposed in [START_REF] Huang | A new total variation method for multiplicative noise removal[END_REF]-see [START_REF] Candès | New multiscale transforms, minimum total variation synthesis. Applications to edge-preserving image reconstruction[END_REF] at the end of this paragraph.

Let us come to the potential function ϕ in the regularization term. In their pioneering work, Tikhonov

and Arsenin [START_REF] Tikhonov | Solutions of Ill-Posed Problems[END_REF] considered ϕ(t) = t 2 ; however it is well known that this choice for ϕ leads to smooth images with flattened edges. Based on a fine analysis of the minimizers of F v as solutions of PDE's on Ω, Rudin, Osher and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF] exhibited that ϕ(|∇u(t)|) = ∇u(t) 2 , where . 2 is the L 2 -norm, leads to images involving edges. The resultant regularization term is known as Total Variation (TV). However, whatever smooth data-fitting is chosen, this regularization yields images containing numerous constant regions (the well known stair-casing effect), so that textures and fine details are removed, see [START_REF] Nikolova | Weakly constrained minimization. Application to the estimation of images and signals involving constant regions[END_REF]. The method in [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF] is of this kind and operates only on the image domain; the fitting term is derived from

(3) and the criterion reads

F S (Σ) = ρ log Σ(t) + S(t) Σ(t) dt + Σ TV , (14)
where ρ depends on K. The denoised image Ŝ0 = arg min Σ F S exhibit constant regions, as seen in Figs. 4(e) and 5(e) in Section 6. We also tried to first restore the log-image û by minimizing

F v (u) = ρ u -v 2 + u TV (15)
and the sought after image is of the form Ŝ0 = B exp(û) where B stands for the bias correction explained in section 4. Because of the exponential transform, there is no stair-casing, but some outliers remain visible-see Figs. 4(c) and 5(c); nevertheless, the overall result is very reasonable. The result of [START_REF] Rudin | Nonlinear total variation based noise removal algorithm[END_REF] was at the origin of a large amount of papers dedicated to constructing edge-preserving convex potential functions, see e.g. [START_REF] Acar | Analysis of bounded variation penalty methods for ill-posed problems[END_REF][START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF][START_REF] Vogel | Iterative method for total variation denoising[END_REF], and for a recent overview, [START_REF] Aubert | Mathematical problems in image processing[END_REF]. Even though smoothness at the origin alleviates stair-casing, a systematic drawback of the images restored using all these functions ϕ is that the amplitude of edges is underestimated-see e.g. [START_REF] Nikolova | Analysis of the recovery of edges in images and signals by minimizing nonconvex regularized least-squares[END_REF]. This is particularly annoying if the sought-after function has neat edges or spiky areas since the later are subjected to erosion. A very recent method proposed in [START_REF] Huang | A new total variation method for multiplicative noise removal[END_REF] restores the discrete log-image using the log-likelihood of ( 9) and a regularized TV; more precisely,

F v (u, w) = i u[i] + S[i]e -u[i] + ρ 0 u -w 2 + ρ w TV , (16)
where the denoised log-image û is obtained using alternate minimization on u and w. The TV term here is regularized via uw 2 2 and the resultant denoised image is given by Ŝ0 = exp( ŵ). The results present some improvement with respect to the method proposed in [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF], at the expense of two regularization parameters (ρ and ρ 0 ) and twho stopping rules for each one of the minimization steps.

Hybrid methods

Hybrid methods [START_REF] Bobichon | Regularized multiresolution methods for astronomical image enhancement[END_REF][START_REF] Candès | New multiscale transforms, minimum total variation synthesis. Applications to edge-preserving image reconstruction[END_REF][START_REF] Chan | Total variation improved wavelet thresholding in image compression[END_REF][START_REF] Coifman | Combining the calculus of variations and wavelets for image enhancement[END_REF][START_REF] Durand | Reconstruction of wavelet coefficients using total variation minimization[END_REF][START_REF] Froment | Artifact free signal denoising with wavelets[END_REF][START_REF] Malgouyres | Mathematical analysis of a model which combines total variation and wavelet for image restoration[END_REF][START_REF] Malgouyres | Minimizing the total variation under a general convex constraint for image restoration[END_REF] combine the information contained in the large coefficients y[i], obtained according to [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF], with pertinent priors directly on the log-image u.

Remark 1 Such a framework is particularly favorable in our case since the noise W n[i], i ∈ I in the coefficients y[i], i ∈ I, have a nearly Gaussian distribution-see Fig. 1(d).

Although based on different motives, hybrid methods amount to define the restored function û as minimize Φ(u)

subject to û ∈ {u : |(W (u -v)) [i]| ≤ µ i , ∀i ∈ I} .
If the use of an edge-preserving regularization, such as TV for Φ is a pertinent choice, the strategy for the selection of parameters {µ i } i∈J is more tricky. This choice must take into account the magnitude of the relevant data coefficient y[i]. However, deciding on the value of µ i based solely on y[i], as done in these papers, is too rigid since there are either correct data coefficients that incur smoothing (µ i > 0), or noisy coefficients that are left unchanged (µ i = 0). A way to alleviate this situation is to determine (µ i ) i∈I based both on the data and on a prior regularization term. Following [START_REF] Nikolova | Minimizers of cost-functions involving nonsmooth data-fidelity terms. Application to the processing of outliers[END_REF][START_REF] Nikolova | A variational approach to remove outliers and impulse noise[END_REF], this objective is carried out by defining restored coefficients x to minimize the non-smooth objective function, as explained below.

A specialized hybrid criterion

Given the log-data v obtained according to (4), we first apply a frame transform as in [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF] to get

y = W v = W u 0 + W n.
We systematically denote by x the denoised coefficients. The noise contained in the i-th datum reads n, w i whose distribution is displayed in Fig. 1(d). The low frequency approximation coefficients carry important information about the image. In other words, when w i is low frequency, then n, w i has a better SNR than other coefficients. Therefore, as usual, a good choice is to keep them intact at this preprocessing stage. Let I * ⊂ I denote the subset of all such elements of the frame. Then we apply a hard-thresholding to all coefficients except those contained in

I * y TH [i] def = T H y[i] , ∀i ∈ I \ I * , (17)
where the hard-thresholding operator T H reads [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] T

H (t) = 0 if |t| ≤ T, t otherwise. ( 18
)
The resultant set of coefficients is systematically denoted by y TH . We choose an underoptimal threshold T in order to preserve as much as possible the information relevant to edges and to textures, an important part of which is contained in the small coefficients. Let us consider

v TH = i∈I1 W v[i] w i , (19)
where

I 1 = {i ∈ I \ I * : |y[i]| > T }. (20)
The image v TH contains a lot of artifacts with the shape of the w i for those y[i] that are noisy but above the threshold T , as well as a lot of information about the fine details in the original log-image u 0 . In all cases, whatever the choice of T , the image of the form v TH is unsatisfactory-see Fig. 2 (b-h).

We will restore x based on the under-thresholded data y TH . We focus on hybrid methods of the form:

     x = arg min x F (x) F (x) = Ψ(x, y TH ) + Φ( W x), (21)
where Ψ is a data-fitting term in the domain of the frame coefficients and Φ is an edge-preserving regularization term bearing the prior on the sought-after log-image û. The latter sought-after log-image

û reads û = W x . ( 22
)
Next we analyze the information content of the coefficients y TH that give rise to our log-image û. Let us denote

I 0 = I \ I 1 ∪ I * = {i ∈ I \ I * : |y[i]| ≤ T }. (23)
We are mostly interested by the information borne by the coefficients relevant to I 0 and I 1 .

(I 0 ) The coefficients y[i] for i ∈ I 0 are usually high-frequency components which can be of the two types described below. This analysis clearly defines the goals that the minimizer x of F y is expected to achieve. In particular,

x must involve an implicit classification between coefficients that fit to y TH exactly and coefficients that are restored according to the prior term Φ. In short, restored coefficients have to fit y TH exactly if they are in accordance with the regularization term Φ and have to be restored via the later term otherwise.

Since [START_REF] Nikolova | Minimizers of cost-functions involving nonsmooth data-fidelity terms. Application to the processing of outliers[END_REF][START_REF] Nikolova | A variational approach to remove outliers and impulse noise[END_REF] we know that criteria F y where Ψ is non-smooth at the origin (e.g. ℓ 1 ) can satisfy x

[i] = y TH [i]
for coefficients that are in accordance with the prior Φ, while the others coefficients are restored according to Φ, see also [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF]. For these reasons, we focus on a criterion on the form

F y (x) = Ψ(x) + Φ(x) (24)
where

Ψ(x) = i∈I λ i (x -y TH )[i] = i∈I1∪I * λ i |(x -y)[i]| + i∈I0 λ i |x[i]| , (25)
Φ(x) = Ω |∇ W x| ds = W x TV . (26)
Note that in [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF], as well as in what follows, we write F y in place of F yT H in order to simplify the notations.

In the pre-processing step [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF] we would not recommend the use of a shrinkage function other than T H since it will alter all the data coefficients y T , without restoring them faithfully. In contrast, we base our restoration on data y TH where all non-thresholded coefficients keep the original information on the sought-after image.

The theorem stated next ensures the existence of a minimizer for F y as defined in ( 24) and ( 25)- [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF].

Its proof can be found in [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF].

Theorem 1 [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF] For y ∈ ℓ 2 (I) and T > 0 given, consider F y as defined in [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF], where Ω ∈ I R 2 is open, bounded and its boundary ∂Ω is Lipschitz. Suppose that 1. {w i } i∈I is a frame of L 2 (Ω) and the operator W is the pseudo-inverse of W ;

2. λ min = min

i∈I λ i > 0.
Then F y has a minimizer in ℓ 2 (I).

Let us remind that the minimizer of F y is not necessarily unique. Given y, denote

G y def = x ∈ ℓ 2 (I) : F y (x) = min x∈ℓ 2 (I) F y (x) . (27)
Hopefully, for every sample of the preprocessed data y TH , the set G y is convex and corresponds to images W x which are very similar since they share the same level lines. The theorem below confirms this assertion and is proven in [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF]. Some orientations for the choice of λ i were investigated in [START_REF] Durand | Denoising of frame coefficients using l1 data-fidelity term and edge-preserving regularization[END_REF]. If i ∈ I 1 , the parameter λ i should be close to, but below the upper bound w i TV , since above this bound, the coefficients y[i] cannot be changed. For i ∈ I 0 , a reasonable choice is

λ i = max k =i Ω (∇ w i ) T ∇ w k |∇ w k | ds ,
where . T stands for transposed. If λ i is below this bound, some neighboring outliers might not be properly removed although Gibbs oscillations are better reduced. Another important remark is that, for some multiscale transforms, the bounds discussed above are constant. In the proposed model, we use only two values for λ i , depending only on the set I ǫ the index i belongs to.

We focus on the coefficients of a curvelets transforms of the log-data because (a) such a transform captures efficiently the main features of the data and (b) it is a tight-frame which is helpful for the subsequent numerical stage.

Minimization for the log-image

Let us rewrite the minimization problem defined in ( 24) and ( 25)-( 26) in a more compact form: find x such that F y (x) = min x F y for

min x F y = Ψ + Φ, where Ψ(x) = Λ(x -y TH ) 1 , for Λ = diag(λ i ) i∈I , Φ(x) = W x TV . (28)
where λ i are the coefficients given in [START_REF] Combettes | A Douglas-Rachford splittting approach to nonsmooth convex variational signal recovery[END_REF]. Clearly, Ψ and Φ are proper lower-semicontinuous convex functions, hence the same holds true for F y . The set G y introduced in ( 27) is non-empty by Theorem 1 and can be rewritten as

G y = {x ∈ ℓ 2 (I) x ∈ (∂F y ) -1 (0)},
where ∂F y stands for subdifferential operator. Minimizing F y amounts to solving the inclusion 0 ∈ ∂F y (x) , or equivalently, to finding a solution to the fixed point equation

x = (Id + γ∂F y ) -1 (x) , (29)
where (Id + γ∂F y ) -1 is the resolvent operator associated to ∂F y , γ > 0 is the proximal stepsize and Id is the identity map on the Hilbert space ℓ 2 (I). The proximal schematic algorithm resulting from (29), namely

x (t+1) = (Id + γ∂F y ) -1 (x (t) ),
is a fundamental tool for finding the root of any maximal monotone operator [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF][START_REF] Rockafellar | Monotone operators and the proximal point algorithm[END_REF], such as e.g. the subdifferential of a convex function. Since the resolvent operator (Id + γ∂F y ) -1 for F y in (28) cannot be calculated in closed-form, we focus on iterative methods.

Splitting methods do not attempt to evaluate the resolvent mapping (Id + γ(∂Ψ + ∂Φ)) -1 of the combined function F y , but instead perform a sequence of calculations involving separately the resolvent operators (Id + γ∂Ψ) -1 and (Id + γ∂Φ)) -1 . The latter are usually easier to evaluate, and this holds true for our functionals Ψ and Φ in [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF].

Splitting methods for monotone operators have numerous applications for convex optimization and monotone variational inequalities. Even though the literature is abundant, these can basically be systematized into three main classes: the forward-backward [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF][START_REF] Tseng | Applications of a splitting algorithm to decomposition in convex programming and variational inequalities[END_REF][START_REF] Tseng | A modified forward-backward splitting method for maximal monotone mappings[END_REF], the Douglas/Peaceman-Rachford [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF],

and the little-used double-backward [START_REF] Lions | Une méthode itérative de resolution d'une inéquation variationnelle[END_REF][START_REF] Passty | Ergodic convergence to a zero of the sum of monotone operators in hilbert space[END_REF]. A recent theoretical overview of all these methods can be found in [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF][START_REF] Eckstein | A family of projective splitting methods for the sum of two maximal monotone operators[END_REF]. Forward-backward can be seen as a generalization of the classical gradient projection method for constrained convex optimization, hence it inherits all its restrictions. Typically, one must assume that either Ψ or Φ is differentiable with Lipschitz continuous gradient, and the stepsizes γ must fall in a range dictated by the gradient modulus of continuity; see [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF] for an excellent account. Obviously, forward-backward splitting is not adapted to our functional (29).

Specialized Douglas-Rachford splitting algorithm

The Douglas/Peaceman-Rachford family is the most general preexisting class of monotone operator splitting methods. Given a fixed scalar γ > 0, let

J γ∂Ψ def = (Id + γ∂Ψ) -1 and J γ∂Φ def = (Id + γ∂Φ) -1 . (30)
Given a sequence µ t ∈ (0, 2), this class of methods can be expressed via the recursion

x (t+1) = 1 - µ t 2 Id + µ t 2 (2J γ∂Ψ -Id) • (2J γ∂Φ -Id) x (t) . (31)
Since our problem (28) admits solutions, the following result ensures that iteration (31) converges for our functional F y .

Theorem 3 Let γ > 0 and µ t ∈ (0, 2) be such that t∈I N µ t (2µ t ) = +∞. Take x (0) ∈ ℓ 2 (I) and consider the sequence of iterates defined by [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF]. Then, (x (t) ) t∈I N converges weakly to some point x ∈ ℓ(I)

and J γ∂Φ (x) ∈ G y .
This statement is a straightforward consequence of [24, Corollary 5.2]. For instance, the sequence µ t = 1, ∀t ∈ I N, satisfies the requirement of the latter theorem.

It will be convenient to introduce the reflection operator

rprox ϕ = 2prox ϕ -Id. ( 32
)
where prox ϕ is the proximity operator of ϕ according to in Definition 1. Using ( 35) and ( 32), the Douglas-Rachford iteration given in (31) becomes

x (t+1) = 1 - µ t 2 Id + µ t 2 rprox γΨ • rprox γΦ x (t) . ( 33
)
Below we compute the resolvent operators J γ∂Ψ and J γ∂Φ with the help of Moreau proximity operators.

Proximal calculus

Proximity operators were inaugurated in [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF] as a generalization of convex projection operators.

Definition 1 (Moreau [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF]) Let ϕ be a proper, lower-semicontinuous and convex function defined on a Hilbert space H. Then, for every x ∈ H, the function z → ϕ(z) + xz 2 /2, for z ∈ H, achieves its infimum at a unique point denoted by prox ϕ x. The operator prox ϕ : H → H thus defined is the proximity operator of ϕ.

By the minimality condition for prox ϕ , it is straightforward that ∀x, p ∈ H we have

p = prox ϕ x ⇐⇒ x -p ∈ ∂ϕ(p) ⇐⇒ (Id + ∂ϕ) -1 = prox ϕ . (34)
Then [START_REF] Durand | Reconstruction of wavelet coefficients using total variation minimization[END_REF] reads

J γ∂Ψ = prox γΨ and J γ∂Φ = prox γΦ . (35)

Proximity operator of Ψ

The proximity operator of γΨ is established in the lemma stated below.

Lemma 1 Let x ∈ ℓ 2 (I). Then

prox γΨ (x) = y TH [i] + T S γλi (x[i] -y TH [i]) i∈I , (36)
with

T S γλi (z[i]) = max 0, z[i] -γλ i sign(z[i]) . (37)
Proof. Ψ is an additive separable function in each coordinate i ∈ I. Thus, solving the proximal minimization problem of Definition 1 is also separable. For any convex function ϕ and v ∈ ℓ 2 (I), put

ψ(.) = ϕ(. -v). Then p = prox ψ (x) ⇐⇒ x -p ∈ ∂ψ(p) ⇐⇒ (x -v) -(p -v) ∈ ∂ϕ(p -v) ⇐⇒ p -v = prox ϕ (x -v) ⇐⇒ p = v + prox ϕ (x -v) .
For each i ∈ I, we apply this result with v = y TH [i] and ϕ(z

[i]) = γλ i |z[i]|. Noticing that prox ϕ = T S γλi
is soft-thresholding with threshold γλ i , leads to [START_REF] Huang | A new total variation method for multiplicative noise removal[END_REF]. ⋄

Note that now

rprox γΨ (x) = 2 y TH [i] + T S γλi (x[i] -y TH [i]) i∈I -x . (38)

Proximity operator of Φ

Clearly, Φ(x) = • TV • W (x) is a pre-composition of the TV-norm with the linear operator W . Computing the proximity operator of Φ for an arbitrary W may be intractable. We adopt the following assumptions:

(w1) W : ℓ 2 (I) → L 2 (Ω) is surjective;

(w2) W W = Id and W = c -1 W * for 0 < c < ∞, where W * stands for the adjoint operator; note that we also have W * W = c Id;

(w3) W is bounded.

For any

z(t) = (z 1 (t), z 2 (t)) ∈ I R 2 , t ∈ Ω, we set |z(t)| = z 1 (t) 2 + z 2 (t) 2 . Let X = L 2 (Ω) × L 2 (Ω) and
•, • X be the inner product in X , and

• p , for p ∈ [1, ∞] the L p -norm on X . We define B γ ∞ (X ) as the closed L ∞ -ball of radius γ in X , B γ ∞ def = z ∈ X : z ∞ ≤ γ = z = (z 1 , z 2 ) ∈ X : |z(t)| ≤ γ, ∀t ∈ Ω , (39)
and

P B γ ∞ (X ) : X → B γ ∞ (
X ) the associated projector; it is easy to check that the latter is equal to the proximity operator of the indicator function of B γ ∞ (X ). The expression for prox γΦ is given in the next lemma while the computation scheme to solve item (ii) is stated in Lemma 3.

Lemma 2 Let x ∈ ℓ 2 (I) and B γ ∞ (X ) is as defined above.
(i) Denoting by prox c -1 γ • TV (u) the proximity operator of the (c -1 -scaled) TV-norm, we have

prox γΦ (x) = Id -W • Id -prox c -1 γ • TV • W (x) ; ( 40
) (ii) Furthermore, prox c -1 γ • TV (u) = u -P C (u) , (41)
where

C = div(z) ∈ L 2 (Ω) z ∈ C ∞ c (Ω × Ω), z ∈ B γ/c ∞ (X ) . (42)
Proof. Since W is surjective, its range is L 2 (Ω) which is closed. Moreover, the domain dom( •

TV ) = L 2 (Ω) as well, so that cone dom • TV -range W = {0} which is a closed subspace of L 2 (Ω).
Reminding that • TV is lower bounded, continuous and convex, it is clear that all assumptions required in [25, Proposition 11] are satisfied. Applying the same proposition yields statement (i).

We focus next on (ii). Note that C in ( 42) is a closed convex subset since B γ/c ∞ (X ) is closed and convex, and div is linear; thus the projection P C is well defined.

Let us remind that the Legendre-Fenchel (known also as the convex-conjugate) transform of a function ϕ : H → I R, where H is an Hilbert space, is defined by

ϕ * (w) = sup u∈dom(ϕ) w, u -ϕ(u) ,
) 43
One can see also [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]Lemma 2.10] for an alternate proof. It is easy to check that the conjugate function of a norm is the indicator function ı of the ball of its dual norm, see e.g. [9, Eq.(2.7)]; thus

c -1 γ • TV * (z) = 0 if z ∈ C , +∞ if z ∈ C ,
where C is given in [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF]. Using Definition 1, it is straightforward that prox c -1 γ . TV * = P C .

Identifying c -1 γ . TV with ϕ and c -1 γ . TV * with ϕ * , equation ( 43) leads to statement (ii). The proof is complete. ⋄ Note that our argument [START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalized gaussian and complexity priors[END_REF] for the computation of prox c -1 γ • TV (u) is not used in [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF], which instead uses conjugates and bi-conjugates of the objective function.

Remark 2 In view of equations ( 41) and ( 42), we one can see that the term between the middle parentheses in equation ( 40) admits a simpler form:

Id -prox c -1 γ • TV = P C .
Using [START_REF] Eckstein | A family of projective splitting methods for the sum of two maximal monotone operators[END_REF] along with ( 40)-( 41) we easily find that

rprox γΦ (x) = Id -2W • Id -prox c -1 γ • TV • W (x) = Id -2W • P C • W (x) . (44)

Calculation of the projection P C in (41) in a discrete setting

In what follows, we work in the discrete setting. We consider that that W is an M × N tight frame with M = #I ≫ N , admitting a constant c > 0 such that

W W = Id and W = c -1 W T ( hence W T W = c Id).
(This is the discrete equivalent of assumption (w2).) We also suppose that W :

ℓ 2 (I) → ℓ 2 (Ω) is surjective.
Next we replace X by its discrete counterpart,

X = ℓ 2 (Ω) × ℓ 2 (Ω) where Ω is discrete with #Ω = N. (45)
We denote the discrete gradient by ∇ and consider Div : X → ℓ2 (Ω) the discrete divergence defined by analogy with the continuous setting 1 as the adjoint of the gradient Div = -∇ * ; see [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF].

Unfortunately, the projection in [START_REF] Malgouyres | Minimizing the total variation under a general convex constraint for image restoration[END_REF] does not admit an explicit form. The next lemma provides an iterative scheme to compute the proximal points introduced in Lemma 2. In this discrete setting, C in (42) admits a simpler expression:

C = Div(z) ∈ ℓ 2 (Ω) z ∈ B γ/c ∞ (X ) . , (46)
where

B γ/c
∞ (X ) is defined according to [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF].

Lemma 3

We adapt all assumptions of Lemma 2 to the new discrete setting, as explained above. Consider the forward-backward iteration

z (t+1) = P B 1 ∞ (X ) z (t) + β t ∇ Div(z (t) ) -cu/γ for 0 < inf t β t ≤ sup t β t < 1/4, (47)
where

∀(i, j) ∈ Ω P B 1 ∞ (X ) (z)[i, j] =    z[i, j] if |z[i, j]| ≤ 1; z[i, j] |z[i, j]| otherwise . Then (i) (z (t) ) t∈I N converges to a point ẑ ∈ B 1 ∞ (X ); (ii) c -1 γDiv(z (t) ) t∈I N converges to c -1 γDiv(ẑ) = (Id -prox c -1 γ • TV )(u).
Proof. Given u ∈ ℓ 2 (Ω), the projection ŵ = P C (u) defined by ( 41) and ( 46) is unique and satisfies ŵ = arg min

w∈C 1 2 u -w 2 = arg min 1 2 c γ u -w 2 subject to w = Div(z) for z ∈ B 1 ∞ (X ) ŵ = Div(ẑ) where ẑ = arg min z∈B 1 ∞ (X ) 1 2 c γ u -Div(z) 2 (48)
This problem can be solved using a projected gradient method (which is a special instance of the forwardbackward splitting scheme) whose iteration is given by [START_REF] Nikolova | Analysis of the recovery of edges in images and signals by minimizing nonconvex regularized least-squares[END_REF]. This iteration converges weakly to a minimizer of (48)-see [24, Corollary 6.5], provided that the stepsize β t > 0 satisfies sup t β t < 2/δ 2 , where δ is the spectral norm of the div operator. It is easy to check that δ 2 ≤ 8-see e.g. [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF].

Set

ω (t) = c -1 γDiv(z (t) ), ∀t ∈ I N and ω = c -1 γDiv(ẑ).
Thus,

ω (t) -ω 2 = γ c 2 Div(z (t) ) -Div(ẑ) 2 = γ c 2 -∇ Div(z (t) ) -Div(ẑ) , z (t) -ẑ X , (49)
where we use the fact that -∇ is the adjoint of Div. Let D z denote the gradient of a scalar-valued function of z, not to be confused with the discrete gradient operator ∇ of an image. The gradient of the function 1 2 cu/γ -Div(z) 2 with respect to z is D z 1 2 cu/γ -Div(z) 2 = -∇ (Div(z)cu/γ). This relation together with the Schwarz inequality applied to (49) lead to

ω (t) -ω 2 ≤ γ c 2 ∇ Div(z (t) ) -∇ Div(ẑ) 2 z (t) -ẑ 2 = γ c 2 ∇ Div(z (t) ) -cu/γ -∇ Div (ẑ) -cu/γ 2 z (t) -ẑ 2 = 0.5 γ c 2 D z cu/γ -Div(z (t) ) 2 -D z cu/γ -Div(ẑ) 2 2 z (t) -ẑ 2 . (50)
From [24, Theorem 6.3], we deduce that the series

t∈I N D z ( cu/γ -Div(•) ) (z (t) ) -D z ( cu/γ -Div(.) ) (ẑ) 2 2
is convergent. Inserting this property in (50) and using the fact that the sequence (z (t) ) t∈I N is bounded (as it converges weakly with ẑ 2 < lim inf t→∞ z (t) 2 ), it follows that ω (t) converges strongly to ω. This completes the proof. ⋄

The forward-backward splitting-based iteration proposed in [START_REF] Nikolova | Analysis of the recovery of edges in images and signals by minimizing nonconvex regularized least-squares[END_REF] to compute the proximity operator of the TV-norm is new and different from the projection algorithm of Chambolle [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF], even tough the two algorithms bear some similarities. The forward-backward splitting allows to derive a sharper upperbound on the stepsize β t than the one proposed in [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF]-actually twice as large. Let us remind that it was observed in [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF] that the bound 1/4 still works in practice. Here we prove why thus is really true.

Comments on the Douglas-Rachford scheme for F y

A crucial property of the Douglas-Rachford splitting scheme (33) is its robustness to numerical errors that may occur when computing the proximity operators prox Ψ and prox Φ , see [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF]. We have deliberately omitted this property in [START_REF] Froment | Artifact free signal denoising with wavelets[END_REF] for the sake of simplicity. This robustness property has important consequences: e.g. it allows us to run the forward-backward sub-recursion (47) only a few iterations to compute an approximate of the TV-norm proximity operator in the inner iterations, and the Douglas-Rachford is still guaranteed to converge provided that these numerical errors are under control. More precisely, let a t ∈ ℓ 2 (I) be an error term that models the inexact computation of prox γΦ in [START_REF] Malgouyres | Mathematical analysis of a model which combines total variation and wavelet for image restoration[END_REF], as the latter is obtained through [START_REF] Nikolova | Analysis of the recovery of edges in images and signals by minimizing nonconvex regularized least-squares[END_REF]. If the sequence of error terms (a t ) t∈I N and stepsizes (µ t ) t∈I N defined in Theorem 3 obey t∈I N µ t a t < +∞, then the Douglas-Rachford algorithm [START_REF] Froment | Artifact free signal denoising with wavelets[END_REF] converges weakly [24, Corollary 6.2].

In our case, using 200 inner iterations in [START_REF] Nikolova | Analysis of the recovery of edges in images and signals by minimizing nonconvex regularized least-squares[END_REF] was sufficient to satisfy this requirement.

Remark 3 Owing to the splitting framework and proximal calculus, we have shown in Lemma 2 that the bottleneck of the minimization algorithm is in the computation of the proximity-operator of the TVnorm. In fact, computing prox • TV amounts to solving a discrete ROF-denoising. Our forward-backward iteration is one possibility among others, and other algorithms beside [START_REF] Chambolle | An algorithm for total variation minimization and application[END_REF] have been proposed to solve the discrete ROF-denoising problem. While this paper was submitted, our attention was drawn to an independent work of [START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF] who, using a different framework, derive an iteration similar to [START_REF] Nikolova | Analysis of the recovery of edges in images and signals by minimizing nonconvex regularized least-squares[END_REF] to solve the ROF. Another parallel work of [START_REF] Zhu | Duality-Based Algorithms for Total-Variation-Regularized Image Restoration[END_REF] propose an application of gradient projection to solving the dual problem [START_REF] Passty | Ergodic convergence to a zero of the sum of monotone operators in hilbert space[END_REF]. We are of course aware of max-flow/min-cut type algorithms, for instance the one in [START_REF] Chambolle | On total variation minimization and surface evolution using parametric maximum flows[END_REF].

We have compared our whole denoising procedure using our implementation of prox • TV and the maxflow based implementation that we adapted from the code available at [1]. We obtained similar results, although the max-flow-based algorithm was faster, mainly because it uses the ℓ 1 approximation of the discrete gradient, namely ( ∇u)[i, j]

1 = u[i + 1, j] -u[i, j] + u[i, j + 1] -[i, j] .
Let us remind that this approximation for the discrete gradient does not inherits the rotational invariance property of the L 2 norm of the usual gradient.

4 Bias correction to recover the sought-after image

Recall from (4) that u 0 = log S 0 and set û = W x(NDR) as the estimator of u 0 , where N DR is the number of Douglas-Rachford iterations in [START_REF] Froment | Artifact free signal denoising with wavelets[END_REF]. Unfortunately, the estimator û is prone to bias, i.e.

E [û] = u 0 -b û.
A problem that classically arises in statistical estimation is how to correct such a bias. More importantly is how this bias affects the estimate after applying the inverse transformation, here the exponential. Our goal is then to ensure that for the estimate Ŝ of the image, we have E Ŝ = S 0 . Expanding Ŝ in the neighborhood of E [û], we have

E [exp û] = exp (E [û])(1 + Var [û] /2 + R 2 ) = S 0 exp (-b û)(1 + Var [û] /2 + R 2 ) , (51)
The authors of [START_REF] Xie | SAR speckle reduction using wavelet denoising and markov random field modeling[END_REF] propose a direct estimate of the bias b û using the obvious argument that the noise n in the log-transformed image has a non-zero mean ψ 0 (K)log K. A quick study shows that the functions (1 + ψ 1 (K)/2) and exp(log Kψ 0 (K)) are very close for K reasonably large. Thus, the two bias corrections are equivalent. Even though the bias correction approach we propose can be used in a more general setting.

Full algorithm to suppress multiplicative noise

Now, piecing together Lemma 1, Lemma 2 and Theorem 3, we arrive at the multiplicative noise removal algorithm:

Task: Denoise an image S contaminated with multiplicative noise according to (2).

Parameters:

The observed noisy image S, number of iterations N DR (Douglas-Rachford outer iterations)

and N FB (Forward-Backward inner iterations), stepsizes µ t ∈ (0, 2), 0 < β t < 1/4 and γ > 0, tight-frame transform W and initial threshold T (e.g. T = 2 ψ 1 (K)), regularization parameters λ 0,1 associated to the sets I 0,1 .

Specific operators:

(

a) T S γλi (z) = max 0, z[i] -γλ i sign(z[i]) i∈I , ∀z ∈ I R #I . (b) P B 1 ∞ (X ) (z)[i, j] =    z[i, j] if |z[i, j]| ≤ 1; z[i, j] |z[i, j]| otherwise, ∀(i, j) ∈ Ω.
(c) ∇ and Div-the discrete versions of the continuous operators ∇ and div.

(d) ψ 1 (•) defined according to (8) (built-in Matlab function, otherwise see [START_REF] Press | Numerical recipes, the art of scientific computing[END_REF]).

Initialization:

• Compute v = log S and transform coefficients y = W v. Hard-threshold y at T to get y TH . Choose an initial x (0) .

Main iteration:

For t = 1 to N DR ,

(1) Inverse curvelet transform of x (t) according to u (t) = W x (t) .

(2) Initialize z (0) ; For s = 0 to N FB -1

z (s+1) = P B 1 ∞ (X ) z (s) + β t ∇ Div(z (s) ) -c γ u (t) .
(3) Set z (t) = z (N FB) .

(4) Compute w (t) = c -1 γ Div(z (t) ).

(5) Forward curvelet transform: α (t) = W w (t) .

(6) Compute r (t) = rprox γΦ (x (t) ) = x (t) -2α (t) . ( 7) By [START_REF] Lions | Une méthode itérative de resolution d'une inéquation variationnelle[END_REF] compute

q (t) = rprox γΨ • rprox γΦ x (t) = 2 y TH [i] + T S γλi r (t) [i] -y TH [i] i∈I -r (t) .
(8) Update x (t+1) using (33):

x (t+1) = 1 - µ t 2 x (t) + µ t 2 q (t) .

End main iteration

Output:

Denoised image Ŝ = exp W x (NDR) (1 + ψ 1 (K)/2).
Remark 4 (Computation load) The bulk of computation of our denoising algorithm is invested in applying W and its pseudo-inverse W . These operators are of course never constructed explicitly, rather they are implemented as fast implicit analysis and synthesis operators. Each application of W or W cost O(N log N ) for the second generation curvelet transform of an N -pixel image [START_REF] Candès | Fast discrete curvelet transforms[END_REF]. If we define N DR and N FB as the number of iterations in the Douglas-Rachford algorithm and the forward-backward sub-iteration, the computational complexity of the denoising algorithm is of order N DR N FB 2N log N operations.

Experiments

In all experiments carried out in this paper, our algorithm was run using second-generation curvelet tight frame along with the following set of parameters: ∀t, µ t ≡ 1, β t = 0.24, γ = 10 and N DR = 50. The initial threshold T was set to 2 ψ 1 (K). For comparison purposes, some very recent multiplicative noise removal algorithms from the literature are considered: the AA algorithm [START_REF] Aubert | A variational approach to remove multiplicative noise[END_REF] minimizing the criterion in [START_REF] Bobichon | Regularized multiresolution methods for astronomical image enhancement[END_REF], and the Stein-Block denoising method [START_REF] Chesneau | Stein block thresholding for image denoising[END_REF] in the curvelet domain, applied on the log transformed image.

The latter is a sophisticated shrinkage-based denoiser that thresholds the coefficients by blocks rather than individually, and has been shown to be nearly minimax over a large class of images in presence of additive bounded noise (not necessarily Gaussian nor independent). We also tried the "naive" method, called L2-TV, where û minimizes [START_REF] Candès | Fast discrete curvelet transforms[END_REF] and the denoised image is given after bias correcion according to [START_REF] Rudin | Multiplicative denoising and deblurring: Theory and algorithms[END_REF]. No without surprise, one realizes that the results are quite good, even though some persistent outliers remain quite visible. This again raises the persistent question of relevance of PSNR (or even MAE) as a measure of perceptual restoration quality. For fair comparison, the hyperparameters for all competitors were tweaked to reach their best level of performance on each noisy realization. The results are depicted in Fig. 3, Fig. 4 and Fig. 5. Our denoiser clearly outperforms its competitors both visually and quantitatively as revealed by the PSNR and MAE values. The PSNR improvement brought by our approach is up to 4dB on the Shepp-Logan phantom, and is ∼ 1dB for Lena and Boat.

Note also that a systematic behavior of AA algorithm is its tendency to lose some important details and the persistence of a low-frequency ghost as it can be seen on the error maps on the third row in Figs. 4 and5.

Conclusions

This work proposes quite an original, efficient and fast method for multiplicative noise removal. The latter is a difficult problem that arises in various applications relevant to active imaging system, such as laser imaging, ultrasound imaging, SAR and many others. Multiplicative noise contamination involves inherent difficulties that severely restrict the main restoration algorithms.

The main ingredients of our method are: (1) consider the log-data to restore a log-image; (2) preprocess the log-fata using and under-optimal hard-thresholding of its tight frame coefficients; (3) restore the logimage using a hybrid criterion composed of an ℓ 1 data-fitting for the coefficients and a TV regularization in the log-image domain; (4) restore the sought-after image using an exponential transform along with a pertinent bias correction. The resultant algorithm is fast, its consistency and convergence are proved theoretically.

The obtained numerical results are really encouraging since they outperform the most recent methods in this field.

Figure 1 :

 1 Figure 1: Noise distributions.

Figure 2 :

 2 Figure 2: (a) Noisy Lena for K = 1. (b) Noisy Lena obtained via averaging, see (1), for K = 10. (c)-(h) Denoising of data v shown in (c) where the curvelet trasform of v are hard-thresholded according to (17)-(19) for different choices of T where σ = ψ 1 (K) is the standard deviation of the noise n. The displayed restorations correspond to exp v TH , see[START_REF] Chan | Total variation improved wavelet thresholding in image compression[END_REF].

Theorem 2 [ 29 ]

 229 If x1 and x2 are two minimizers of F y (i.e. x1 ∈ G y and x2 ∈ G y ), then∇ W x1 ∝ ∇ W x2 ,a.e. on Ω.In other words, W x1 and W x2 have the same level lines.In words, images W x1 and W x2 are obtained one from another by a local change of contrast which is usually invisible for to the naked eye.

  and that ϕ * is a closed convex function. If ϕ is convex, proper and lower semi-continuous, the original Moreau decomposition [42, Proposition 4.a] tells us that prox ϕ + prox ϕ * = Id . (

where R 2

 2 is expectation of the Lagrange remainder in the Taylor series. One can observe that the posteriordistribution of û is nearly symmetric, in which case R 2 ≈ 0. That is, b û ≈ log(1 + Var [û] /2) to ensureunbiasedness. Consequently, finite sample (nearly) unbiased estimates of u 0 and S 0 are respectively û + log(1 + Var [û] /2), and exp (û) (1 + Var [û] /2). Var [û] can be reasonably estimated by ψ 1 (K), the variance of the noise n in (4) being given in[START_REF] Antoniadis | Regularization of wavelet approximations[END_REF]. Thus, given the restored log-image û, our restored image read: Ŝ = exp (û) (1 + ψ 1 (K)/2) .

  The denoising algorithms were tested on three images: Shepp-Logan phantom, Lena and Boat all of size 256 × 256 and with gray-scale in the range[1, 256]. For each image, a noisy observation is generated by multiplying the original image by a realization of noise according to the model in (2)-(3) with the choice µ = 1 and K = 10. For a N -pixel noise-free image S 0 and its denoised version by any algorithm Ŝ, the denoising performance is measured in terms of peak signal-to-noise ratio (PSNR) in decibels (dB)

10 (Figure 3 :

 103 Figure 3: Performance comparison with Shepp-Logan phantom (256 × 256). (a) Original. (b) Noisy µ = 1, K = 10. (c) Denoised with Stein-block thresholding in the curvelet domain [21] PSNR=24.73dB, MAE=4. (d) Denoised with our algorithm PSNR=31.25dB, MAE=1.87. (e)-(f) Errors (restoredoriginal) for (c)-(d).

More precisely, let u ∈ ℓ

(Ω) be of size m × n, N = mn. We write ( ∇u)[i, j] = u[i + 1, j] -u[i, j], u[i, j + 1] -u[i, j] with boundary conditions u[m + 1, i] = u[m, i], ∀i and u[i, n + 1] = u[i, n], ∀i; then for z ∈ X , we have (Div(z))[i, j] = z 1 [i, j] -z 1 [i -1, j] + z 2 [i, j] -z 2 [i, j -1] along with z 1 [0, i] = z 1 [m, i] = z 2 [i, 0] = z 2 [i, n] = 0, ∀i.


Liangcai Xu

Benoît Fauqué

Zengwei Zhu

Zbigniew Galazka

Klaus Irmscher

Kamran Behnia

Thermal conductivity of bulk In 2 O 3 single crystals

The transparent semiconductor In2O3 is a technologically important material. It combines optical transparency in the visible frequency range and sizable electric conductivity. We present a study of thermal conductivity of In2O3 crystals and find that around 20 K, it peaks to a value as high as 5,000 WK -1 m -1 , comparable to the peak thermal conductivity in silicon and exceeded only by a handful of insulators. The amplitude of the peak drastically decreases in the presence of a type of disorder, which does not simply correlate with the density of mobile electrons. Annealing enhances the ceiling of the phonon mean free path. Samples with the highest thermal conductivity are those annealed in the presence of hydrogen. Above 100 K, thermal conductivity becomes sample independent. In this intrinsic regime, dominated by phonon-phonon scattering, the magnitude of thermal diffusivity, D, becomes comparable to many other oxides, and its temperature dependence evolves towards T -1 . The ratio of D to the square of sound velocity yields a scattering time which obeys the expected scaling with the Planckian time.

I. INTRODUCTION

In 2 O 3 is a technologically appealing wide-gap semiconductor (see Ref. [1] for a review). Its main polymorph crystallizes in a cubic bixbyite structure [2]. The cubic (primitive) unit cell of this structure contains 16 (8) molecules and 80 (40) atoms. Introducing Sn dopants turns this semiconducting solid to a transparent conductor known as indium tin oxide (ITO) [3]. Electrical transport properties of doped In 2 O 3 with the carrier density in the range 10 18 -10 20 cm -3 has documented a roomtemperature mobility in the range carrier of ∼ 100-200 cm 2 V -1 s -1 . The high transparency in the visible spectrum makes this system an established transparent conducting material.

Amorphous or granular indium oxide (InO x ) has been the subject of an independent fundamental research activity focusing on a fascinating two-dimensional superconductor [4]. It is a platform to study metal-anomalous metal-superconductor transitions [5] and host to a number of remarkable low-temperature transport properties [6]. The connection between this two-dimensional electronic system and the superconducting ground state of crystalline doped thin films of In 2 O 3 [7] is yet to be clarified. Only recently [8,9], bulk single crystals of In 2 O 3 became available. The thermal conductivitives of indium oxide thin films [10] and single crystals [START_REF] Galazka | Transparent Semiconducting Oxides Bulk Crystal Growth and Fundamental Properties[END_REF] were previously reported at room temperature. Now, thanks to the direct growth of single crystals from the melt [9], it is possible to measure thermal conductivity of single crystals with a steady-state setup down to cryogenic temperatures.

In this paper, we present a study of temperature dependence of thermal conductivity in single crystals of In 2 O 3 . We find that the peak thermal conductivity in In 2 O 3 crystal can reach as high as 5000 WK -1 m -1 . In a few other solids, such as diamond [START_REF] Onn | [END_REF]13] and sapphire (crystalline Al 2 O 3 ) [14], this peak is exceeded in magnitude. Most oxides show a much more modest maximum thermal conductivity. We find that this peak is drastically reduced by disorder. However, the level of disorder is not simply set by dopant concentration and depends on annealing conditions. Samples annealed in the presence of hydrogen are those showing the highest thermal conductivity. Above 100 K, thermal conductivity κ in all samples converges to a similar magnitude. In this intrinsic regime (governed by phonon-phonon scattering), the ratio of the thermal diffusivity to the square of sound velocity respects a boundary observed recently for all insulating solids [15,16]. Comparison with Al 2 O 3 indicates that, in agreement with a picture recently proposed to explain this boundary, the compound with the lighter atom has a larger thermal diffusivity driven by its superior sound velocity [17].

II. EXPERIMENTAL

Bulk In 2 O 3 single crystals were grown with a crystal growth technique under the name levitation-assisted selfseeding crystal growth method (LASSCGM) as described in detail by Galazka et al. [9,[START_REF] Galazka | Transparent Semiconducting Oxides Bulk Crystal Growth and Fundamental Properties[END_REF]18] that enabled to obtain bulk single crystals up to 35 mm in diameter and 12 mm thick. The as-grown crystals were dark red or dark brown and highly conducting with a free electron concentration of (1-4)×10 18 cm -3 and Hall electron mobility of 130-170 cm 2 V -1 s -1 [START_REF] Galazka | Transparent Semiconducting Oxides Bulk Crystal Growth and Fundamental Properties[END_REF]18]. Annealing in the presence of oxygen and hydrogen allowed tuning the free electron concentration between 1×10 18 and 6×10 19 cm -3 and increased the room-temperature Hall electron mobility to 200 cm 2 V -1 s -1 . Wafers fabricated from bulk single crystals were characterized by a very narrow rocking with the arXiv:2008.13519v4 [cond-mat.mtrl-sci] 12 Jan 2021 full width at half maximum (FWHM) < 50 arcsec and the x-ray texture indicating no grain boundaries [9,18]. The dislocation etch pit density was found to be between 5×10 4 and 3×10 5 cm -2 [START_REF] Galazka | Transparent Semiconducting Oxides Bulk Crystal Growth and Fundamental Properties[END_REF]18]. The crystal orientation, the annealing conditions, and the size of the single crystals studied in this work are listed in Table I. All crystals had all identical dimensions.

We measured thermal conductivity in a steady-state configuration with one heater and two thermometers.

Like in a setup used before [19], we used chromelconstantan thermocouples as thermometers between 10 and 300 K [see inset in Fig. 1(b)]. In the case of sample 6, the data were completed by measurements using Cernox chips as thermometers between 2 and 30 K. The two sets of data were found to be consistent at intermediate temperatures.

We can identify two potential sources of experimental uncertainty. The first is the finite size of the thermal contacts. In this case, the sample length was 5 mm and the typical contact size was 0.1 mm, implying an uncertainty of 2% on the absolute value. The other source of uncertainty is heat loss along the measuring wires. The geometric factor of our measuring wires, given the thickness (25 µm) and their length ( 15 mm), is three orders of magnitude smaller than the geometric factor of our samples. Thermally resistant manganin wires were used to connect thermocouples to the heat bath, and therefore, given manganin thermal conductivity (κ manganin ≈ 1 W/Km at 10K and ≈ 20 W/Km at 300K [START_REF] Touloukian | Thermal Conductivity: Metallic Elements and Alloys[END_REF]) the margin due to heat loss is less than 3% of the sample thermal conductivity.

III. RESULTS

Figure 1 presents our main results. Near room temperature, thermal conductivity increases with cooling and its magnitude is the same in all crystals studied. The room-temperature thermal conductivity of our bulk samples, κ (300 K) = 15±2 WK -1 m -1 , is comparable to the value of 13.1 WK -1 m -1 obtained by the laser flash method [START_REF] Galazka | Transparent Semiconducting Oxides Bulk Crystal Growth and Fundamental Properties[END_REF] and three times larger than what was reported for thin films of ITO (≈ 5 WK -1 m -1 ) [10]. In all samples, κ(T ) shows a broad peak. The magnitude of this κ peak is strongly sample dependent. It varies from 100 WK -1 m -1 in sample 1 to 5000 WK -1 m -1 in sample 6. The higher the κ peak , the lower the temperature at which it occurs. These features point to a smooth evolution in disorder among these samples. As we see below, however, identifying this disorder is not trivial.

We measured the temperature dependence of resistivity and Hall constant R H , of all samples. Carrier density was determined using n = 1/R H e (e is the elementary

1 0 1 0 0 2 3 5 0 1 0 1 1 0 2 1 0 3 1 0 4 1 0 5 1 0 6 1 0 0 # 1 # 2 # 3 # 4 # 5 # 6 l p h o n o n ( n m ) T ( K ) FIG. 2.
Phonon mean free path in In2O3, extracted from thermal conductivity, specific heat, and the sound velocity.

The electric resistivity and the carrier density of these samples are listed in Table I. There is no visible correlation between the electronic properties and phonon mean free path. Note that in the best sample ( 6), the mean free path approaches the sample thickness of 0.5 mm at low temperature.

charge). This single-band expression is justified by the fact that the minimum of the conduction band is at the center of the Brillouin zone [1]. The room-temperature values are given in Table I.

Mobile electrons can also conduct heat. In general, the measured thermal conductivity has two contributions: κ meas = κ ph + κ e . The latter is the electronic contribution to thermal conductivity, which can be estimated using the electric resistivity and the Wiedemann-Franz law: κ e = L 0 T /ρ, where L 0 = 2.45 × 10 -8 V 2 /K 2 is the Sommerfeld value. Even in the most metallic like sample, κ e κ meas . Indeed, for ρ ≈ 2mΩcm, one expects κ e /T =L 0 /ρ ≈ 10 -3 WK -2 m -1 . Even at room temperature, this is a tiny fraction of the measured signal. We can therefore safely conclude that κ meas ≈ κ ph .

The effective phonon mean free path ph can be estimated using the kinetic equation κ = 1 3 C ph v s ph . (Note that such an estimation neglects the possible variation of the mean free path among modes and frequencies.) Figure 2 shows the temperature dependence of the effective phonon mean free path ph extracted from our thermal conductivity data, the reported specific heat [START_REF] Cordfunke | [END_REF]22]), and the estimated sound velocity v s = 4.4 km/s [10]. The smooth evolution of ph (T ) (Fig. 2) is reminiscent of the way thermal conductivity of an insulator evolves upon controlled introduction of disorder, such as NaCl crystals with Ag colloids [23]. Surprisingly, however, the smooth evolution of ph (T ) is not correlated with the electrical properties. As an example, the sample showing the longest ph (sample 6) and the one showing the shortest (the as-grown sample, sample 1) are only slightly different from the charge transport viewpoint. The room-temperature carrier density in sample 6 is five times larger than in sample 1. How can then the phonon mean free path in sample 6 be a hundred times longer, in spite of having a higher concentration of scattering centers?

IV. DISCUSSION

The behavior seen here is to be contrasted with the case of doped semiconductors such as silicon and germanium [24] or doped oxides such as SrTi 1-x Nb x O 3 [15], where the introduction of Nb dopants leads to a steady enhancement in electrical conductivity (because of the introduction of mobile electrons) and a steady decrease in thermal conductivity (because of a rise in the number of the scattering centers).

Our observation implies that the case of In 2 O 3 is more complicated and disorder length scale is not simply the interdopant distance. In the as-grown sample sample 1, the deviation from the intrinsic regime starts below 150 K (see Fig. 2). The typical wavelength of an acoustic phonon at this temperature is λ ph (150 K)= hvs k B T =1.4 nm. On the other hand, in sample 2, the deviation starts at 30 K, when the phonon wavelength is λ ph (30 K) = 7 nm. These numbers set the lower boundary to the size of scattering centers. The most prominent candidates for playing this role are clusters of oxygen vacancies which manifest themselves as In particles. They are formed during melt growth in the cooling crystal and have been identified in a recent study combining electron microscopy, light scattering techniques, and optical absorption spectroscopy [25]. The concentration of these clusters, which can be as large as several tens of nanometers [25], diminishes upon annealing. This type of extended disorder is the most plausible candidate for scattering long wavelength acoustic phonons in low-κ samples.

The two samples with highest κ peak were annealed in the presence of H 2 (see Table I). Hydrogen is a possible dopant and a shallow donor impurity [26]. H atoms can replace O atoms at their site [27]. The small size of H dopants may allow the survival of high phonon thermal conductivity (which requires small disorder) despite finite electric conductivity (which requires atomic substitution). Note, however, that according to first-principles calculations [28], hydrogen may also occupy interstitial positions and form complexes with indium vacancies. More extensive studies are required to pin down the origin of the link between annealing conditions and peak thermal conductivity. Note that in the most heat conducting , the phonon mean free path ph (T ) becomes as long as 0.3 mm below 10 K, only slightly shorter than the sample thickness (0.5 mm), implying that phonons are close to ballistic.

A theoretical description of the thermal conductivity can be achieved by calculating the complex phonon spectrum of this solid and solving the Boltzmann-Peierls transport equation [29]. In the absence of a rigorous solution, there are approximate treatments. The most influential is the one introduced by Callaway [30], which in contrast to relaxation time approximation takes into ac-

1 0 1 0 0 1 0 -1 1 0 0 1 0 1 1 0 2 1 0 3 1 0 4 1 0 5 S r T i O 3 D i a m o n d A l 2 O 3 I n 2 O 3 κ ( W / K m ) T ( K ) S i FIG. 3.
Comparison with other solids. Comparison of the thermal conductivity in the best In2O3 sample with diamond [13], silicon [34], Al2O3 (sapphire) [14], and a perovskite SrTiO3 [15]. Thermal conductivity in In2O3 is comparable to a ternary oxide such as SrTiO3 at room temperature, but peaks to a value comparable to silicon's peak at low temperature.

count normal collisions between phonons [31]. The starting point for a phenomenological approach is the Debye temperature. However, because of the complex phonon spectrum due to the large number of atoms per unit cell, a single well-defined Debye temperature is absent in this system. Published Debye temperatures of In 2 O 3 range from 420 to 811 K (see Ref. [32] for more details on this issue). In a forthcoming paper, Subedi [33] will present ab initio calculations of the phonon spectrum and lattice thermal conductivity. His results give a reasonable account of our experimental data for the cleanest samples in a wide temperature window.

Tantalizingly, in sample 6, where the low-temperature ph is close to the sample thickness and phonons become ballistic, ph (T) peaks at 10 K, a possible signature of phonon hydrodynamics [36]. In a handful of solids [15,[36][37][38], the temperature dependence of the phonon mean free path displays a local maximum followed by a local minimum. This has been attributed to a regime, first identified by Gurzhi [39], where normal scattering among phonons prevails over both ballistic scattering and umklapp scattering leading to the so-called Poiseuille flow [36]. It is too early, however, to conclude at this stage that the peak in ph (T ) is a Poiseuille peak. No clear Knudsen minimum is visible down to 5 K. Moreover, the thermal conductivity of the best sample may also be affected by extended disorder, since the phonon wavelength at 10 K is still "only" 30 nm long. This peak is absent in sample 5, where ph falls well below the crystal thickness. Note that above 100 K, ph of all samples merge and approach a T -1 behavior. At room temperature, ph becomes as short as 4 nm, still four times longer than the very long lattice parameter. 4. Thermal diffusivity in the intrinsic regime: (a) Temperature dependence of thermal diffusivity, D, in In2O3, compared to a few representative solids (Si, SrTiO3, PbTe ) [15] and Al2O3 [14,35]. The dashed line represents T -1 . At room temperature, D in In2O3 (as well as in Al2O3 and in Si) approaches this behavior without fully attaining it.

2 0 0 3 0 0 4 0 0 5 0 0 1 0 0 1 0 -2 1 0 -1 1 0 0 1 0 1 D ( c m 2 / s ) T ( K ) S i A l 2 O 3 I n 2 O 3 S r T i O 3 P b T e FIG.
We come back to this intrinsic regime below.

Figure 3 compares the thermal conductivity of the best In 2 O 3 sample (sample 6) with a selection of other insulators. In all, κ peaks at an intermediate temperature. The peak observed in In 2 O 3 (5000 WK -1 m -1 at 20 K) is smaller than what is reported in diamond (30000 WK -1 m -1 at 60 K) [13] and in sapphire [14] and slightly exceeds the peak in silicon [34]. Remarkably, the crystal structure of In 2 O 3 is quite distinct from either the diamond structure or the α-Al 2 O 3 structure. It has two nonequivalent sites for indium atoms [1,2]. An indium atom is surrounded by oxygen atoms either in an octahedral coordination or in a trigonal prismatic coordination. In both cases, indium has six immediate oxygen neighbors. In the octahedral coordination all these O atoms are equidistant from the In atom. In the trigonal prismatic coordination, however, there is a slight difference in distance between three pairs of O atoms. Each oxygen atom is surrounded by one indium atom of the first type and three indium atoms of the second type [40]. Thus, the octahedral coordination is present in both sapphire and indium oxide, two binary compounds with remarkably large low-temperature thermal conductivity. We note that other insulators such as Tl 2 O 3 [41] crystallize in cubic bixbyite structure, but their thermal conductivity remains unexplored. Let us also note that in β-Ga 2 O 3 , κ peaks to 560 WK -1 m -1 [42]. This is to be contrasted with ternary oxides, such as perovskites [43], where κ peak is in the range of 40 WK -1 m -1 , two orders of magnitude smaller than the two binary oxides of Fig. 3. As seen in the figure, and discussed in more detail below, this difference disappears at room temperature.

Let us now turn our attention to the magnitude of thermal diffusivity in the intrinsic regime. Above 100 K, thermal conductivity in all six samples becomes identical The τ /τP ratio vs the vM /vs ratio in different systems adopted from Ref [17] with the new data point for In2O3. Note that τ in In2O3 was derived at room temperature. Its value at higher temperature may be slightly lower. The blue region contains column IV elements and III-V semiconductors, which do not participate in this trend. Most other systems lie in a green stripe whose boundaries are set by 0.6 ± 0.4vM /vs. Comparisons between In2O3 and Al2O3 and between LiF and RbI reveal the role played by atomic mass in setting the sound velocity (see text).

5 1 0 1 5 2 0 2 5 4 1 1 0 0 . 5 6 0 I n 2 O 3 A l 6 S i 2 O 1 3 T i O 2 M g A l 2 O 4 Y A l O 3 Y A G S r T i O 3 C a F 2 N i O L a A l O 3 B a T i O 3 A g B r A g C l M g 2 S i O 4 R b B r R b I T h O 2 U O 2 K I N a I K B r Y 2 O 3 P b S B a F 2 R b C l N a B r S i O 2 G G G K C l C a O S r F 2 C a T i O 3 F e 2 O 3 C a M g S i 2 O 6 Z r S i O 4 K N b O 3 L i A l S i 2 O 6 M g S i O 3 L i F A l 2 O 3 τ/ τ P v M / v s
in magnitude and in temperature dependence. Heat conduction in this temperature range is governed by phononphonon scattering. Our data combined with specific heat reported previously [START_REF] Cordfunke | [END_REF] allow us to extract the thermal diffusivity (D = κ c , where c is the specific heat per volume). As seen in Fig. 4, D evolves faster than T -1 below 250 K, but becomes close to the expected T -1 behavior near room temperature. The figure compares In 2 O 3 with a few other solids. The room-temperature thermal diffusivity of In 2 O 3 is only 7 mm 2 s -1 . This is more than one order of magnitude lower than Si (D ≈ 90 mm 2 s -1 ), but close to what has been observed in SrTiO 3 (D ≈ 4 mm 2 s -1 ) [15] and other ternary oxides [35,44] including cuprates [45].

In this intrinsic regime, the asymptotic T -1 temperature dependence of thermal diffusivity represents the temperature dependence of the phonon-phonon scatter-ing cross section. Recently, an empirical lower boundary to D in this regime was noticed [15,16]. Expressing thermal diffusivity, D, as a function of the Planckian time, τ P = /k B T , one can write

D = sv 2 s τ P . (1)
Here, s is a dimensionless parameter (representing the ratio of the scattering time τ to τ P ). It was observed [16] that s does not fall below unity in any known solid. According to our data, s ≈ 14 in In 2 O 3 at room temperature. Note that since D is still evolving faster than T -1 below 300 K, our data can only put an upper boundary to s. In the case of Al 2 O 3 , data above room temperature [35] reduce s, since the slope of the T -1 thermal diffusivity evolves between 300 and 1200 K. Thus, In 2 O 3 , in spite of being an excellent conductor of heat at low temperatures, emerges more akin to low-diffusivity oxides [15,44,45].

A recent theoretical account of the empirical boundary to thermal diffusivity [15,16,45] was put forward by Mousatov and Hartnoll [17], who argued that there is a quantum-mechanical bound to the sound velocity in solids, because the energy required to hold the crystal together and the highest phonon energy are linked. Note that, the classical distribution of phonons above their Debye temperature notwithstanding, the speed of sound in solids can be linked to , as recently demonstrated by Trachenko et al. [46]. Mousatov and Hartnoll [17] defined a "melting velocity": v M = k B T M a/ (T M is the melting temperature, k B is the Boltzmann constant, and a is the interatomic distance). In their picture, the smaller v M /v s , the closer would be the system to the Planckian bound to dissipation and the lower should be the τ /τ P ratio. Their examination of available data on numerous solids (see Fig. 5, adapted from Ref. [17]) confirmed this conjecture. As seen in the figure, most solids cluster in the green stripe, where τ /τ P is roughly proportional to v M /v s . The list includes various families of ionic salts as well as binary and ternary oxides (in which as we saw above the magnitude of κ peak can differ by several orders of magnitude). On the other hand, as previously noticed [17], adamantine solids (in the blue region) do not follow this trend. One may speculate that the threephonon scattering phase space [47] in these solids is too small for the relevance of the Planckian bound.

In the Mousatov-Hartnoll plot [17], the vertical and horizontal axes are experimentally measurable quantities and there is a visible correlation across materials. Plots showing correlations between two structural properties across materials were reported by Grimvall and Sjödin [48]. These plots are to be considered in the context of the impressive success of the Lindemann criterion for melting in spite of its shortcoming, namely, an exclusive focus on the solid phase and the neglect of the energetics of its liquid competitor [49].

We can add In 2 O 3 to this plot. It has a melting temperature of T M = 2223 K [9], an interatomic distance of a = 0.21 nm [2], and a sound velocity of v s = 4.4 km/s [10]. As a result, one finds v M /v s ≈ 15. Our roomtemperature data set an upper bound to τ /τ P ≤ 15. As seen in Fig. 5, these numbers put In 2 O 3 inside the green stripe in company of numerous other oxides and ionic salts. According to Ref. [32], the sound velocity is 6.4 km/s. Taking this value instead of 4.4 Km/s [10] would change both ratios, but it would keep In 2 O 3 inside the green strip.

Interestingly, In 2 O 3 and Al 2 O 3 reside at the opposite ends of the green stripe. This is understandable. Mousatov and Hartnoll [17] already noticed that what puts two alkali halides (LiF and RbI) with identical crystal structures (and comparable room-temperature thermal conductivities) far away in this plot is the difference in the atomic masses. Given that the atomic mass of In (115) is about four times larger than Al (27), the same can be said about In 2 O 3 and Al 2 O 3 . In both cases, heaviness reduces the sound velocity [46] and lightness promotes zero-point motion and proximity to the Planckian boundary [17]. Note, however, that In 2 O 3 and Al 2 O 3 have different crystal structures.

In summary, we measured the thermal conductivity of the transparent semiconductor In 2 O 3 and found a peak at 20 K whose amplitude is surpassed only by a few solids. The peak thermal conductivity was found to be extremely sensitive to disorder. In the intrinsic regime, thermal diffusivity respects a boundary found in other insulators.

FIG. 1 .

 1 FIG. 1. Thermal conductivity of In2O3 and its sample dependence. (a) Thermal conductivity κ of four different samples of In2O3. Note the strong variation in the magnitude of maximum κ among the samples. Purple bars represent the margin of experimental reproducibility in repeated measurements on sample 6. Because of its large thermal conductance, the temperature difference caused by the application of the heat current was small, leading to a larger experimental uncertainty. (b) Comparison of thermal conductivity in three samples having different plane orientations. The heat current was always applied along the longest axis. No anisotropy is expected in this cubic system and the difference is due to different levels of disorder. The inset shows the measurement setup with thermocouples.

FIG. 5 .

 5 FIG. 5.Scattering time, Planckian time, sound velocity, and melting velocity: The τ /τP ratio vs the vM /vs ratio in different systems adopted from Ref[17] with the new data point for In2O3. Note that τ in In2O3 was derived at room temperature. Its value at higher temperature may be slightly lower. The blue region contains column IV elements and III-V semiconductors, which do not participate in this trend. Most other systems lie in a green stripe whose boundaries are set by 0.6 ± 0.4vM /vs. Comparisons between In2O3 and Al2O3 and between LiF and RbI reveal the role played by atomic mass in setting the sound velocity (see text).

TABLE I .

 I Samples: The In2O3 single crystals investigated by the present study. n and ρ refer to room temperature carrier concentration and resistivity. κ peak refers to the peak thermal conductivity. All samples were cut to 5 × 3 × 0.5 mm 3 . The orientations relate to the largest plane (5 × 3 mm 2 ).

	Sample no. Orientation	Annealing conditions	n(10 17 cm -3 ) ρ(mΩcm) κ peak (WK -1 m -1 )
	1	(111)	As grown	12	17	100
	2	(111)	900 o C, 40 h, air	2	95	300
	3	(001)	1000 o C, 40 h, air	3	89	450
	4	(110)	1000 o C, 40 h, O2	12	45	480
	5	(111)	1000 o C, 40 h, O2 and 650 o C, 10 h, H2 (5%)+Ar	180	2.1	1200
	6	(111)	1000 o C, 40 h, O2 and 500 o C, 10 h, H2 (5%)+Ar	61	10	5000

V. ACKNOWLEDGMENTS We thank Sean Hartnoll and Kostya Trachenko for discussions and Detlef Klimm for a critical reading of the manuscript. This work was supported in France by the Agence Nationale de la Recherche (Grants No. ANR-18-CE92-0020-01 and No. ANR-19-CE30-0014-04), and by Jeunes Equipes de l Institut de Physique du Collège de France and by a grant attributed by the Ile de France regional council. It was partly performed in the framework of GraFOx, a Leibniz-Science Campus partially funded by the Leibniz Association -Germany. It was also supported by the National Key Research and Development Program of China (Grant No. 2016YFA0401704), the National Science Foundation of China (Grants No. 51861135104 and No. 11574097) and Fundamental Research Funds for the Central Universities (Grant No. 2019kfyXMBZ071). L. X. acknowledges a Ph.D. scholarship by the China Scholarship Council (CSC).


Demetrios Xenides
email: xenides@uop.gr

Panaghiotis Karamanis
email: panagiotis.karamanis@univ-pau.fr

Dipole Polarizability of C 28 and its Counterparts Nb 4 B 18 and Ta 4 B 18 . Insights from a Density Functional Theory (DFT) Endeavour

We report on a preliminary investigation of the non linear optical (NLO) properties and in particular dipole polarizability. The target species are two perfect tetrahedral nanoclusters Nb4B18 and Ta4B18, along with their nanofullerene counterpart that is C28. Our study based on density functionals (DFs) that have gained popularity among the scientific community. In addition we performed Hartree-Fock calculations known for not including dynamic electron correlation. The DF obtained values are characterized by some dispersion, with maximal differences to be around 5 %, in all three cases. Given that the DFT introduces a fuzzy percentage of electron correlation sets the observed convergence of HF values to DFT ones is at least surprising. Furthermore, it should be said that though the values can be characterized as accurate their reliability should not be taken for granted. Last, we note the smooth convergence of LC-BLYP, LC-BP86, LC-BPW91 to LC-ωhPBE.

Introduction

The study of the field of nanoclusters emerged as a need to bridge the gap between molecules and matter. Species of this type are building blocks of larger architectures. As such Nb 4 B 18 and Ta 4 B 18 after self-assembling could form the respective face central crystals M 4 B 15 (M=Nb, Ta) [START_REF] Zhang | Perfect Spherical Tetrahedral Metallo-Borospherene Ta4B18 as a Superatom Following the 18-Electron Rule[END_REF].

On the other hand the importance of dipole polarizability (α αβ ) has been elegantly described in the excellent book by Bonin and Kresin [START_REF] Bonin | Electric-Dipole Polarizabilities of Atoms, Molecules, and Clusters[END_REF]. This list has been extended by studies relate it to hardness-softness [START_REF] Hati | Hardness and Electric Dipole Polarizability. Atoms and Clusters[END_REF], aromaticity [START_REF] Jartín | Can aromaticity be connected with molecular polarizability? A theoretical study of benzene isomers and five-membered heterocyclic molecules[END_REF], and other important natural phenomena such as electrostatic interactions that involve molecular polarization [START_REF] Israelachvili | Interactions Involving the Polarization of Molecules[END_REF].

For these reason we have put some effort in calculating dipole polarizablities of Nb 4 B 18 and Ta 4 B 18 in Density Functional Theory (DFT) levels of theory and for the shake of comparison at Hartree-Fock, as well. The DFT methods are still the only alternative as the size of the architecture grows. Therefore, it seems mandatory to check their accuracy and reliability in such calculations.

Computational Details

We have calculated dipole polarizability (α αβ /e 2 α 2 0 E -1 h ) via the Finite Field approach. To ensure the reliability of the method we have used carefully selected field strengths of 0.003 au.

Basis sets of TZVP quality have been used for B, C, Nb and Ta. For the latter two we have also implemented effective core potentials (ECPs) as to take into account the relativistic corrections related to the nature of the two heavy atoms. In particular, these two elements have been described by the def2-TZVP basis set obtained from Basis Set Exchange server [START_REF] Pritchar | A New Basis Set Exchange: An Open, Up-to-date Resource for the Molecular Sciences Community[END_REF]. All calculations have been performed with Gaussian 16 package [START_REF] Frisch | Gaussian[END_REF].

The methods used in the present study are Hartree-Fock (HF) [START_REF] Roothaan | New Developments in Molecular Orbital Theory[END_REF], B- [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic-behavior[END_REF] and B3-{LYP,P86,PW91} [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[END_REF][START_REF] Perdew | Density-functional approximation for the correlation energy of the inhomogeneous electron gas[END_REF][START_REF] Perdew | Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation[END_REF], PBE1PBE [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] and HSEH1PBE [START_REF] Heyd | Erratum: Hybrid functionals based on a screened Coulomb potential[END_REF], CAM-B3LYP [START_REF] Yanai | A new hybrid exchange-correlation functional using the Coulomb-attenuating method (CAM-B3LYP)[END_REF], LC-ωPBE [START_REF] Vydrov | Assessment of a long range corrected hybrid functional[END_REF] as well as LC-{BLYP, BP86, BPW91}, M05 [START_REF] Zhao | Exchange-correlation functional with broad accuracy for metallic and nonmetallic compounds, kinetics, and noncovalent interactions[END_REF] and M052x [START_REF] Zhao | Design of density functionals by combining the method of constraint satisfaction with parametrization for thermochemistry, thermochemical kinetics, and noncovalent interactions[END_REF], M06 and M062x [START_REF] Zhao | The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals[END_REF], M08 [START_REF] Zhao | Exploring the limit of accuracy of the global hybrid meta density functional for main-group thermochemistry, kinetics, and noncovalent interactions[END_REF], ωB97xD [START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF] and BMK [START_REF] Boese | Development of Density Functionals for Thermochemical Kinetics[END_REF] In the case of C 28 ( 5 A 2 ) all the above methods have been applied at the restricted open (RO) formulation.

In Fig. 1 we have drawn the structure of the three nanoclusters under study. We further add a picture of the HOMO and LUMO orbitals of them as to gain information on the shape of the frontier orbitals.

Results and Discussion

In the case of C 28 the majority of DFT methods produce dipole polarizability values that are greater than the respective RO-HF ones. However, there are the LC-BLYP, LC-BP86, LC-BPW91 and LC-ωhPBE methods that diverge from the trend of the rest of DFT ones. That is the LC correction strongly affects the predictive capability of DF methods. Surprisingly enough the values of the LC methods are closer to ROHF ones, as it can be seen in Fig. 2. In particular the ROHF value is 246.4 when the LC-{BLYP, BP86, BPW91, and ωhPBE} are 246.9, 247.7, 245.5 and 246.2, respectively. The others are lying in the vicinity of 255 as it can be also seen in Table 1. The same trend had been observed in a study by Sekino et al [START_REF] Kirtman | Calculation of electric dipole (hyper)polarizabilities by long-range-correction scheme in density functional theory: A systematic assessment for polydiacetylene and polybutatriene oligomers[END_REF] on dipole (hyper)polarizablities of polydiacetylene and polybutatriene. What we have seen from the above made analysis is that in all three cases DFs methods caused a positive (enhancement) or negative (lowering) of the respective HF dipole polarizability values. Furthermore, it is known that DFs methods introducing a (fuzzy) percentage of dynamic electron correlation. That is the reason of gaining their popularity (along with the fact that they need moderate computational resources). In the present study the DFs produced results differ by the most of 5 % from the HF ones. Moreover, it is the nature of them that does not allow for naming one of them as method of preference. That, probably, justifies the question of whether there is a need of a calculation of DFT type when a HF one produces similar results, without the accompanying uncertainty of the DFT methods.

Last but not least, we note the remarkably similarity of dipole polarizability of both M 4 B 18 (M=Nb, Ta). Thus, it is fair to assume that the nature of the heavy atom has a negligible effect on this interesting property.

Conclusion

Having performed restricted open shell calculations on the C 28 ( 5 A 1 ) fullerene we found that most of the functionals enhance the value of dipole polarizability. In striking difference is the LC-({BLYP, BP86, BPW91, ωhPBE}) family of functionals that produces results similar to ROHF ones. The inverse has been observed in M 4 B 18 (M=Nb, Ta) nanoclusters. In particular DF methods lower the value of Hartree-Fock polarizability. This is most intense in the case of LC family. However, in all cases the DFT values are around 5 % lower or greater than the HF ones. This could raise the question whether there is a need of calculation of DF type instead of a pure HF one, that it is known for its limitations. Last, we note that both heavy metal including moieties are characterized by almost identical dipole polarizability values. Therefore, it is fair to assume that there is a negligible effect of the nature of heavy atom on the value of dipole polarizability.

Figure 1 .

 1 Figure 1. Structure and HOMO LUMO orbitals of and C 28 Nb 4 B 18 and Ta 4 B 18 .

  methods are restricted open (RO) ROHF value of 246.4

Figure 2 .

 2 Figure 2. Dipole polarizabilities of Nb 4 B 18 , Ta 4 B 18 and C 28 .

Table 1 .

 1 Dipole polarizability values of C 28 , Nb 4 B 18 , and Ta 4 B 18 at HF and DFT levels of theory (restricted open (RO) shell calculations have been performed in the case of C 28 ). We note that in contrast to the rest of the methods BLYP(364.56) and BMK(367.39) values are larger than the HF ones. The LC correction affected the most the obtained values. This can be seen as a sharp lowering of the obtained values. In particular LC-{BLYP, BP86, BPW91 and ωhPBE} values are 352.24, 349.32, 345.74 and 347.57, respectively. The M05 method produces a value of 353.38, that is the one being close to the LC family. The other functionals produce values that are in between 355 and 360 au.For M being Ta, the BP86(361.98), B3LYP(361.91), M06(361.46) and ωB97xD(362.98) are the ones that are noticeable proximal to HF(362.09) one. A smooth convergence of LC-BLYP,

		C 28	Nb 4 B 18 Ta 4 B 18
	HF	246.4 362.82	362.09
	BLYP	263.9 364.56	364.95
	BP86	264.5 361.15	361.98
	BPW91	262.0 357.07	357.95
	LC-BLYP	246.9 352.24	351.92
	LC-BP86	247.7 349.32	349.49
	LC-BPW91	245.5 345.74	345.99
	B3LYP	260.0 361.49	361.91
	B3P86	259.7 356.96	357.87
	B3PW91	258.6 355.89	356.69
	CAM-B3LYP 254.2 358.40	358.47
	PBE1PBE	257.3 355.56	356.12
	LC-ωhPBE	246.2 347.57	347.79
	HSEH1PBE	258.3 356.84	357.45
	ωB97xD	255.5 361.95	362.98
	BMK	259.8 367.39	369.41
	M05	252.1 353.38	352.77
	M052x	256.5 356.75	357.81
	M06	261.4 360.19	361.46
	M062x	257.5 355.74	356.81
	M08	255.1 354.97	355.86
	Quite opposite is the behaviour of the methods in the case of M 4 B 18 (M=Nb, Ta)
	nanoclusters.
	When M=Nb the BLYP(364.56), BP86(361.15), B3LYP(361.49) and ωB97xD(361.95) values
	are quite close to the HF(362.82).