Hela Machat

Je ´ro ˆme Carrayrou
email: jerome.carrayrou@unistra.fr

Comparison of linear solvers for equilibrium geochemistry 3 computations 4

Comparison of linear solvers for equilibrium geochemistry computations

Hela Machat, Jérôme

separately at each time step and iteratively for some for-46 mulations. As a consequence, the chemistry operator has to 47 be solved at least once per mesh cell per time step. This 48 is one reason for the high computational cost of reactive 49 transport modelling. Some authors have reported that 80 50 to 90 % of the computation time is dedicated to chem-51 ical computation. Many studies have been conducted to 52 reduce the computation time required by reactive transport 53 modelling [START_REF] Carrayrou | Comparison of numerical methods for sim-1263 ulating strongly nonlinear and heterogeneous reactive transport 1264 problems-the MoMaS benchmark case[END_REF]. Some works have explored paralleliza-54 tion [START_REF] Hammond | Modeling mul-1267 ticomponent reactive transport on parallel computers using 1268 Jacobian-Free Newton Krylov with operator-split preconditioning[END_REF], while others have focused on the methods used 55 to solve the transport operator. Nevertheless, improving 56 the resolution of the chemistry operator has been iden-57 tified as a key point. Some authors have attempted to 58 improve the classic Newton-Raphson method [START_REF] Carrayrou | New efficient algorithm for 1274 solving thermodynamic chemistry[END_REF], while 59 others have tested other methods, such as Newton-Krylov 60 [START_REF] Hammond | Modeling mul-1267 ticomponent reactive transport on parallel computers using 1268 Jacobian-Free Newton Krylov with operator-split preconditioning[END_REF][START_REF] Amir | A global method for coupling transport with 1277 chemistry in heterogeneous porous media[END_REF]. In this work, we focus on a specific element of the problem, improving the resolution of the linearized system provided by the Newton-Raphson method. Looking to numerical methods to solve linear systems is not currently a common practice. Indeed, these methods are actually well known [START_REF] Quarteroni | Numerical mathematics[END_REF][START_REF] Axelsson | Direct solution and incomplete factoriza-1283 tion preconditioned conjugate gradient methods[END_REF][START_REF] Barrett | [END_REF][START_REF] Gould | WATSPEC: a computer program for determining 1357 sparse direct solvers for the solution of large sparse, symmetric lin-the equilibrium speciation of aqueous solutions[END_REF][23], and all mathematical packages for scientific computation propose several routines for this task. The motivation of this work comes from the specificity of linear systems that have to be solved for equilibrium chemistry computations. Classic tests for the resolution of linear systems [START_REF] Baldwin | Iterative linear solvers in a 2D radiation-46[END_REF][25][START_REF] Hadjidimos | Successive overrelaxation (SOR) and related reactive transport benchmark of MoMaS with SPECY[END_REF][START_REF] Kalambi | A feasible set for chemical speciation 1372 solution of linear equations[END_REF][28][START_REF] Schenk | Groundwater flow, multicomponent trans-1378 linear equations with PARDISO[END_REF][START_REF] Xue | A direct algorithm for solving ill-conditioned[END_REF] are performed using systems provided by finite element or finite volume discretization, leading to matrices that are large (at least 10,000 unknowns) and sparse. Moreover, even when ill-conditioned systems have been studied [25, [START_REF] Xue | A direct algorithm for solving ill-conditioned[END_REF][START_REF] Pyzara | The influence of a matrix[END_REF], the conditioning of the matrix coming from the chemical system is specific, as underlined by Hoffmann et al. [START_REF] Hoffmann | A parallel global-implicit[END_REF]. For example, Soleymani [START_REF] Soleymani | On performance of SOR method for solving non-1390 systems[END_REF] worked with an ill-conditioned system constructed from 10 ×10 to 20 ×20 Hilbert matrices. The condition numbers then range from 3.5 ×10 13 to 6.2 ×10 28 . In this work, we present by comparing the calculated solution to a reference 114 solution.

115 Because we utilize a markedly small matrix, we did not test 116 parallelization. All the computations were performed on a 117 PC running Windows with 64-bit Fortran 95. Real variables 118 are defined as double-precision real. We prefer double-119 precision computations because all the chemical codes are, 120 to the best of our knowledge, written as double-precision 121 real and because quadruple-precision computation is much 122 more time consuming. Nevertheless, we have tested one 123 method using quadruple-precision real to determine whether 124 this development could be useful. Reference solutions are 125 also computed using quadruple precision. 126 We first present the formulation of the equations describ-127 ing equilibrium reactions and how they are solved using the 128 Newton-Raphson method. This point defines the Jacobian 129 linear system, which is the object of this work We expect to find a method to increase the efficiency of a speciation or reactive transport code. Several properties are required for such a method:

(i) This method should be fast, as the linear system will be solved very often. In the case of reactive transport modelling, the system will be solved at least once per mesh cell per time step. (ii) The method should be very robust. It should be able to solve the linear system even if it is very poorly conditioned. Because the resolution of the linear system is only part of an iterative Newton step, an accurate solution is not absolutely needed. Thus , some advanced codes (e.g. Linear Algebra Package (LAPACK) routine) that check the accuracy of the solution and return an error flag instead of an inaccurate solution are, in this work, less robust than the more rustic routines. (iii) The method should be able to detect failure and return an error flag to the main program so that a recovery procedure can be initiated. In the case of reac-tive transport modelling, this procedure could involve rejecting the current time step and recomputing with a smaller one. (iv) In the initial analysis, the precision of the method is not the key point. Because the linear system resolution is only a part of the Newton-Raphson iterative procedure, reasonable error is acceptable for the linear system inversion. If this error is too large, it will slow the convergence speed for the Newton-Raphson method and decrease the efficiency of the reactive transport code. In this work, errors are estimated 135 discussed. Based on this analysis, we propose an algorithm 136 to optimize the chemical computation in terms of robust-137 ness, accuracy, and efficiency. This algorithm is evaluated 138 on the most selective test. By expanding the limits of the 139 currently used methods, we believe that our new algorithm 140 will contribute to enlarging the field of application of reac-141 tive transport modelling. As a conclusion, we underline the 142 main advances of this work, the new perspectives and the 143 remaining obstacles. = γ X ric coefficients (api,k). The saturation index (SIi) of Cpi is equal to its activity, which is unity for a pure solid phase {Ci } = γi [Ci ] and X } j j j (4)

. Nx SIi = Kpi {Xk} api,k = 1 (2) k=1
The conservation of the total concentration [Tj ] of the j th By substituting the mass action law (1) into the mass conser-175 vation equation (3), the following relationship, which only 176 depends on the components and the precipitated species 177 concentrations, is obtained: 178 component in the system is then written as

NC K i NX NCP Tj = NC ai,j • [Ci] + i=1 NCP i=1 api,j • [Cpi] (3) Tj = i=1 a i,j • • γ i k=1 (γk [Xk]) ai,k + api,j •[Cpi ] i=1 (5)
+ where [Ci] is the concentration of species (Ci) and [Cpi] is the amount of precipitated species (Cpi) per liquid volume unit.

A classic algorithm [START_REF] Carrayrou | New efficient algorithm for 1274 solving thermodynamic chemistry[END_REF][START_REF] De Windt | Presentation and application of the reactive 1397 1334 applied to UO2 oxidative dissolution and uranium migration[END_REF][START_REF] Jauzein | A flexible computer code for modelling Resources[END_REF][START_REF] Parkhurst | 1403 1340 (version 2)-a computer program for speciation[END_REF][START_REF] Van Der Lee | CHESS another speciation and surface complexation computer code[END_REF][START_REF] Westall | HSL: A collection of Fortran codes for large scale scientific 1411 1348[END_REF][41] to describe mineral precipitation or dissolution makes an a priori hypothesis about the existence or non-existence of minerals. In this work, we assume that this hypothesis is proposed. The relationships Combining Eqs. 2 and 5 leads to a set of (NX NCP ) non-179 linear algebraic equations, which can be numerically solved 180through iterative methods. The concentrations of component 181 [Xk] and precipitated species [Cpi] at equilibrium are then 182 determined when the (NX NCP ) objective functions (Yj ) 183 are zero184

+ NC i Y = -T + . • K γi . X • N [ ] a + N . CP • [ ] for = 1 to i=1 a i,j NX . k=1 (γk X k ) i,k i=1 ap i,j Cpi j NX ( 6
)
Yj=NX +i = -1 + Kpi • (γk [Xk])
= max ⎢ ⎞ . . . . , . Yj . ≤ εN-R with εN-R = 10 -12 (7) ⎢ ⎜ . . Yj . ⎟ . . ⎥ ⎣ ⎝ . T j . + . N c . . ai,j [Ci] . i=1 ⎠ j =1,Nx j =Nx+1,Nx+NcP ⎦ 200
The value of the convergence criterion (εN-R = 10 -12 ) is formulation has some weaknesses that are explained later 201 Q4 198 199 set according to usual practice.

The Newton-Raphson method

The historical approach [START_REF] Lichtner | Continuum model for simultaneous chemical reac-1255 tions and mass transport in hydrothermal systems[END_REF][START_REF] Morel | A numerical method for computing equi-176 (2001) 1392 libria in aqueous chemical systems[END_REF][START_REF] Jauzein | A flexible computer code for modelling Resources[END_REF][START_REF] Westall | HSL: A collection of Fortran codes for large scale scientific 1411 1348[END_REF][42][43][44][START_REF] Jennings | Multicomponent equi-1360 23[END_REF][46][START_REF] Yeh | On the solution of ill-conditioned, multispecies components: model development and demonstra-1367 tion. simultaneous, linear, algebraic equations by machine computation[END_REF] However, many authors [START_REF] Amir | A global method for coupling transport with 1277 chemistry in heterogeneous porous media[END_REF][START_REF] Hoffmann | A parallel global-implicit[END_REF][START_REF] Parkhurst | 1403 1340 (version 2)-a computer program for speciation[END_REF][START_REF] Van Der Lee | CHESS another speciation and surface complexation computer code[END_REF]48] have pro-203 posed an alternative approach. Instead of using the com-204 ponent concentrations Xj as the primary variables, they 205 use the logarithm of the component activities ( } ln Xj ). According to this convention, the objective func-207 tions defined by Eq. 8 become conservation equations challenging test for speciation computation [START_REF] Carrayrou | New efficient algorithm for 1274 solving thermodynamic chemistry[END_REF] (see 242

208 Comput Geosci ξ = 206 j Q3 Yj = -Tj NC + i=1 Ki ai,j • γ i • exp Nx k=1 a i,k • ξ k NCP + i=1 api,j • [Cpi] for j = 1 to NX (8)
Z n . = ∂Y j
Appendix 1 (1)). (iii) The pyrite test case describes the dissolution of a 247 pyrite rock in pure water. It has been used to test speciation algorithms [START_REF] Carrayrou | New efficient algorithm for 1274 solving thermodynamic chemistry[END_REF]. Because it involves redox reactions, the stoichiometric coefficients cover a wide range, and the equilibrium constants vary over Zn .

j,k j =, NX NC = ai,j • ai,k • [Ci ] n
several orders of magnitude. This test is used under

1 j,k . k = 1, NX i=1 = ap k-NX,j Appendix 1 (3)).
254 the assumption that no mineral phase is present (see (iv) The MoMaS easy test is the chemical system used 255 [START_REF] Lichtner | Continuum model for simultaneous chemical reac-1255 tions and mass transport in hydrothermal systems[END_REF] ,which must be solved to obtain the progress step

Z n • ∆ξ n , ∆Cp n = Y n+1 -Y n = -Y n ( 12
)
This system yields the values of the component activitiesand precipitate concentrations at the (n + 1)th iteration

ξ n+1 = ξ n + ∆ξ n 228 [Cp] n+1 = [Cp] n + ∆Cp n (13)
To simplify the notations, ξ is used to denote the full vectorof unknowns, Z

Chemical test cases

We choose chemical test cases with various numbers of components. Some of these chemical systems allow the formation of mineral species. Although it is not realistic from a chemical point of view, we test them without minerals and with the maximal possible number of minerals to obtain the largest matrix size. As failure, we include the crash of the method, underflow 338 298 from 10 -12 to 10 -2 mol L -1 (pH = 12 to pH = 2), while or overflow, non-convergence within the maximum number 339 that of component e -is varied from 10 -19 to 10 12 , corresponding to Eh = 0 -.7 to 1.1 V computed using Eq. 14 at25 • C of iterations (for iterative methods), or excessive inaccu-340 racy for some advanced methods (LAPACK routines) that 341 estimate the accuracy of the proposed solution.

342 Solving a linear system (13) using a numerical method 343 Eh = ln e -} RT [START_REF] Yeh | A critical evaluation of recent develop-1260 ments in hydrogeochemical transport models of reactive multi-1261 chemical components[END_REF] produces an approximate solution ( F dξmethod ), and the ref-

344
where T is the temperature (Kelvin), R is the gas constant (8.314 J K mol -1 ), and F is the Faraday constant (96,487 C mol -1 ). This range of electrical potential corresponds to the stability of water at pH values between 2 and 12. For the O2 component, it is not possible to cover the same potential range as e -because of the computation of the reference solution. The activity is varied from 10 -70 to 10 4 , as com-erence method gives (dξref) with accuracy on the same 345 order as the roundoff error. To evaluate the accuracy of the 346 approximate solution, two quantities can be calculated: 347

1. The relative error on the norm, ErrNorm, is obtained by 348 computing the norm of the approximate and reference 349 solution [START_REF] Amir | A global method for coupling transport with 1277 chemistry in heterogeneous porous media[END_REF] 350 puted using Eq. 15 at 25 • C with E 0 = 1.23 V and pH varying from 2 to 12. The potential is then varied from -0.5

Err Norm = |ǁdξmethodǁ -ǁdξrefǁ| ǁdξrefǁ (18) to 1.1 V + } 4 Eh = E 0 + × ln 2 ( 15
)
1. The error on the direction is given by anglemethod , the 351 angle (degrees) between the reference and the approx-352 imate solution calculated using the scalar product of 353 4 F {H2O} these two vectors 354

The activities of the other components vary from 10 -12 to 10 -1 mol L -1 . For each of the three selected components, angle = 360 Arc cos dξ method • dξ ref [START_REF] Quarteroni | Numerical mathematics[END_REF] we compute 30 values equally distributed on a log scale over method 2π ǁdξ method ǁ • ǁdξ ref ǁ

the chosen range, leading to 29,791 different linear systems for each chemical test case. For each of these 29,791 tests, we make only one linear solver (or one Newton step) (exceptin the last section, Section 4, where the iterative Newton method is performed to solve the non-linear system given

. . ǁZǁ1 = n max Z i,j (16)
. reg 1≤ ul j a ≤ r n subintervals. For each subinterval, we compute 364 the mean of each quantity.

365

.

Let Z 0 i i ∞ c ˜ with Z k be the equilibrated Jacobian matrix at iteration k, = Z.
actual developments in speciation codes involve the use of 414 large chemical databases [START_REF] Van Der Lee | CHESS another speciation and surface complexation computer code[END_REF]58,59], leading to an increase 415in These authors defined r k as the vector formed by the ith the size of the chemical systems. The use of iterative

416 k i
methods is also studied in this work.

417 row of Z and c k as the vector formed by the ith column.

We select some direct and iterative solvers according 418 The preconditioning matrices R k and C k are then defined by to the properties of the linear systems and the speciation

⎛ R k = diag⎜ ⎞ 1 ⎟ ⎛ and C k = diag ⎜ ⎞ 1 ⎟
419 computation methods currently in use (Table 1).

Q6

For the direct method, we select LU decomposition [60] 421 [START_REF] Axelsson | Direct solution and incomplete factoriza-1283 tion preconditioned conjugate gradient methods[END_REF] because it was originally used for speciation computa-422 tions by Westall [START_REF] Westall | HSL: A collection of Fortran codes for large scale scientific 1411 1348[END_REF] and Westall et al. [42]. The UMF 423 method [55] has been implemented in the speciation code 424 SPECY [48] in place of the LU approach [START_REF] Carrayrou | New efficient algorithm for 1274 solving thermodynamic chemistry[END_REF]. After show-425

⎝ , ⎠ ¨rk ¨ i=1,Nx+NcP ⎝ , ⎠ ¨ k ¨ i ∞ i=1,Nx+NcP
The equilibrated matrix is defined at iteration k + 1 by ing that the Jacobian matrix is symmetric, we test the 426

Zk+1 = R k • Zk • C k ¨ ¨ ¨ ¨ (21)
DSYTRS subroutine from LAPACK [START_REF]The Linear Algebra Package (LAPACK) can be obtained free of 1406 1343 1344[END_REF], which is based 427 on a UDU decomposition. Because the Jacobian matrix is 428 often positive definite, as shown in Table 3, we test the 429 ing preconditioning matrices and Z the equilibrated matrix. be the result-Instead of solving the linear system (12), we solve

Z • x = -Y (22)
where x C -1 (∆ξ, ∆Cp) and Y R Y . These proce-dures are coded using quadruple-precision reals. The linearsystem [START_REF] Gould | WATSPEC: a computer program for determining 1357 sparse direct solvers for the solution of large sparse, symmetric lin-the equilibrium speciation of aqueous solutions[END_REF] is solved by LU decomposition coded with quadruple-precision real.

Even if the condition numbers of the Jacobian matrices (Z) are very high (10 213.9 for the Fe-Cr mineral test case), the condition numbers of the equilibrated matrices (Z) are much lower: the maximum

= -=
condition number obtained afterequilibration is 10 13.4 . According to Golub and van Loan [START_REF] Golub | Matrix computations[END_REF], if the unit roundoff is approximately 10 -d and the condition number is approximately 10 q , then the Gaussian elimination gives a solution with approximately d q cor-rect digits. Because we use quadruple precision, we obtain d 32, leading to 32 14 18 correct digits. One can thenassume that the reference solution is exact if we compare it to the solutions produced by the tested methods (computed using doubleprecision real).

Selected numerical methods for solving linearsystems

Studies on linear algebra [START_REF] Quarteroni | Numerical mathematics[END_REF]23] present methods for solving linear systems as direct or iterative methods. Histor-ically, speciation codes solved linear systems using direct 2.5 Reference solution methods, such as Gaussian elimination [START_REF] Morel | A numerical method for computing equi-176 (2001) 1392 libria in aqueous chemical systems[END_REF] or LU decom-406 position [START_REF] Carrayrou | New efficient algorithm for 1274 solving thermodynamic chemistry[END_REF][START_REF] Westall | HSL: A collection of Fortran codes for large scale scientific 1411 1348[END_REF]42]. In its actual form, the speciation code 407 Because of the very high condition numbers, it is not pos-SPECY [48] uses unsymmetric multifrontal (UMF) [55] as 408 sible to directly obtain an exact solution. We equilibrate the the linear solver. To the best of our knowledge, no speciation 409 rows and columns of the Jacobian matrices to reduce their code uses iterative methods to solve linear systems. This 410 condition number using the iterative algorithm proposed by point is in accordance with the existing literature, which 411 Knight et al. [START_REF] Knight | A symmetry preserving algorithm 1382 X-ray Anal. 42, 629-633[END_REF] We use the same parameters for all iterative methods: a 447 maximum of 500 iterations and a stopping criterion of 10 -8 .

448

To determine the influence of the stopping criterion, we test 449 the GMRES method using 50,000 maximum iterations and The results obtained using the Jacobi and SOR methods are not detailed here. As previously reported [START_REF] Quarteroni | Numerical mathematics[END_REF], the Jacobi method is inefficient, leading to a very high failure ratio (close to 100 %) even for the easiest test cases. For the SOR , the over-relaxation parameter is the to suppose that this trade-off would not be advantageous. 460 key factor. Unfortunately, we did not find any efficient rela-Nevertheless, an easy way to test preconditioners is pro-461 tionships to define it. For the same chemical system, the best posed by the LAPACK routine DGESVX, which performs 462 value varies from 0.097 to 1.91 without apparent order.

LU decomposition and matrix equilibration depending on 463

We do not extensively test the possibility of using a pre-the estimated condition number. We implement matrix equi-464 conditioner. As stated by Barrett et al. [START_REF] Barrett | [END_REF]: "Since using a libration according to Knight et al. [START_REF] Knight | A symmetry preserving algorithm 1382 X-ray Anal. 42, 629-633[END_REF] to obtain a reference 465 preconditioner in an iterative method incurs some extra cost, solution. We test this preconditioning technique associated 466 both initially for the setup, and per iteration for applying with LU decomposition and the GMRES method, denoted 467 it, there is a trade-off between the cost of constructing and by LU Equil and GMRES Equil in this study. The maximum 468 applying the preconditioner, and the gain in convergence iterations allowed for the equilibration procedure is fixed to 469 speed". In our case, the matrices are very small, leading us 5, according to the recommendations of Knight et al. By plotting the ratio of diagonal dominant matrices 512 Because the usual computations are performed using dou-depending on the condition number (see Appendix 2 513 ble precision, the quadruple precision (dξQUAD) should be (B-1)), it appears that only matrices with very low 514 translated in double-precision real. To avoid overflow, we condition numbers can be diagonal dominant.

515 rescale dξQUAD to ensure its validity. If huge (1.d0) is the (iv) Because the Jacobian matrix is real, symmetric, 516 highest double-precision real represented by the machine, and sometimes diagonal dominant, the question of 517 we rescale dξQUAD to obtain the double-precision solution whether it is positive definite may be posed. In the 518 dξLU QUAD: case of no precipitation, Eq. 11 can be written in 519 dξ = LU QUAD huge (1.d0)

• dξ (23) max .

. In this way, we conserve the direction of the Newton step, even if its norm is changed.

Results and discussion

Properties of the Jacobian matrices

As defined by Eq. 11, the Jacobian matrix has several properties:

(i) The matrix is block-structured, as presented in Table 2. A four-block structure is present if precipitation occurs. (ii) The matrix is symmetric, as shown in Table 2. (iii) In the case of no precipitation, all the diagonal terms of the matrix are strictly positive because they are the sum of a 2 [Ci]. It is then possible for the matrix to be diagonal dominant. We examine this possibility for the selected test case. Table 3 shows the ratio of diagonal dominant Jacobian matrices for all the chemical tests performed according to the previously defined

Z = A T • diag (C) • A ( 24
)
Because the concentrations are positive, the Jacobian matrix 521is analytically positive definite. Nevertheless, this may not 522 be true numerically. We are not able to propose a gen-523eral framework, but we can compute the eigenvalues of the 524 Jacobian matrix and test whether they are positive for all 525 test cases. Table 3 shows that for the gallic acid, Valoc- on the existence of minerals. According to our results, rithm of the component activities as the primary unknown 559 the conditioning of the Z matrix can be evaluated using in Eq. 8. The historical approach [START_REF] Morel | A numerical method for computing equi-176 (2001) 1392 libria in aqueous chemical systems[END_REF] uses the component 560 the following empirical formulas: concentrations as the primary variable and leads to a less 561 interesting Jacobian matrix. Even if the structure presented 562 cond (Z) no mineral = 10 5.30±0.03 × ǁY ǁ 0.9374±0.0008 in Table 2 exists, the matrix is not symmetric. Moreover, 563 the cond (Z) mineral = 10 -3.23±0.08 × ǁY ǁ 1.706±0.002 (25) matrix is worse conditioned (condition number from 564 10 11.2 to 10 49.4 rather than 10 4.06 to 10 24.9 for the pyrite 565

The value and uncertainties are obtained through the least squares method over all cond(Z) and||Y |.| In this way, we propose an estimation of cond(Z) with no computation time cost because the objective function is evaluated during the Newton-Raphson procedure. As shown in Fig. 1, cond(Z) and ||Y || are strongly correlated for large condition num-bers, case). Finally, no specific relation exists between cond(Z) 566 and ||Y|| for the historical formulation.

567

As an example, we show one linear system from the Fe-568 Cr mineral test, corresponding to a condition number of 569 10 187 . One can observe the structure of the matrix and the 570 specificity of the linear system [START_REF] Hadjidimos | Successive overrelaxation (SOR) and related reactive transport benchmark of MoMaS with SPECY[END_REF].

Robustness of the methods

Some other methods (DGETRS, DSYTRS, DGELS, and 592 DGESVX) present a substantial failure ratio, mainly for 593 Figure 2 presents the failure ratio for each method and each high condition number tests (MoMaS easy and Fe-Cr min-594 test case. The presence of minerals prevents the DPOTRS, eral). UMF, SYMMBK, and CGS are robust for the Fe-Cr 595 Inc. CHOLESKY, and Gauss-Seidel methods from solving mineral test but present significant failure ratios for lower-596 the system. If there are minerals present in the chemical conditioned tests, such as MoMaS easy or pyrite mineral. 597 system, a zero-value block appears in the Jacobian matrix, Some methods adapted to symmetric matrices (DSYTRS 598 as shown in Table 2 and Eq. 26. This block makes the Inc. CHOLESKY factorization unappropriated. Because the and SYMMBK) are included in this class of weak methods.

The BiCGStab method has a very low failure ratio and 599 600

Gauss-Seidel method requires division by each diagonal fails only in the two difficult tests (MoMaS easy and Fe-Cr 601 term, this zero-value block makes the method unadapted. mineral). GMRES is the only successful iterative method.

602

The failure of the DPOTRS routine is explained by the Figure 2 shows that some methods are successful for all 603 properties of the Jacobian matrix. As shown in Table 3, the test cases. The most successful direct method is LU, 604 there is no positive definite matrix in the presence of min-while the most successful iterative methods are GMRES 605 erals. In the case of the DPOTRS, Inc. CHOLESKY, and and GMRES 10 -12 . The quadruple-precision method LU 606 Gauss-Seidel methods, the term failure is ambiguous. These QUAD is also successful, which is expected because the 607 methods are expected to fail and should not be used on double-precision LU method is also successful. The use of 608 systems with minerals. If there are no minerals, some matri-an equilibration method as a preconditioner makes LU Equil 609 ces are not positive definite in the MoMaS easy, MoMaS and GMRES Equil successful. For some advanced methods (e.g. LAPACK methods), a failure ratio than for the more rustic methods (LU or Gauss-614 posteriori estimation of the residual and estimation of the Seidel). Because the key point of this work-the resolution 615 condition numbers are performed. If the solution is not suf-of a linear system-is included in the iterative Newton pro-616 ficiently accurate, no solution is given, leading to a higher cedure, it is preferable to obtain an inaccurate solution (so Fig. 2 Mean of the failure ratio for each method and each test case

UNCORREC

the iterative procedure can be continued) than no solution For other iterative methods, the success does not depend 646 (the iterative procedure will be aborted).

on the condition number but on the nature of the matrix 647 Appendix 3 presents the evolution of the failure ratio for and the presence (Appendix 3 (C-18 to C-20)) or absence 648 each test case and each method depending on the condition (Appendix 3 (C-11 to C-17)) of minerals. 649 number.

For the direct methods (Appendix 3 (C-1 to C-5)), for 3.3 Accuracy of the methods 650 small condition numbers corresponding to the test cases gallic acid, Valocchi, and pyrite, no failure occurs. As the con-The accuracy of the methods is evaluated in two ways: (i) 651 dition number increases, the failure ratio also increases for the relative error on the norm ( 18) and (ii) the angle between 652 some methods. MoMaS easy (Appendix 3 (C-4)), MoMaS medium (Appendix 3 (C-6)), and Fe-Cr (Appendix 3 (C-the reference and the calculated solution [START_REF] Quarteroni | Numerical mathematics[END_REF]. the norm of each test case (Fig. 3), some general 655 tendencies are identified. The relative residual tends to 656 increase with the condition number of the system. For 657 and DGELS present some failure for condition numbers direct methods and small condition numbers, the rel-658 greater than 10 15 for the MoMaS easy case. In the presence ative residual is small (10 -10 to 10 -3 ) for the gallic 659 of minerals (Appendix 3 (C-8 and C-9)), for low condition acid, Valocchi, and pyrite test cases. For the itera-660 numbers (the pyrite mineral case), the methods are either tive methods, the relative residual corresponding to an 661 successful (UMF, LU, DSYTRS, DGETRS) or completely accurate resolution for tests with small condition num-662 unsuccessful (DPOTRS). For very high condition numbers bers is approximately 10 -4 . This value corresponds to 663 (Fe-Cr mineral case), the success of the method does not the value of the convergence criteria of the iterative 664 depend on the condition number. We suppose that the con-methods. Iterative methods are more sensitive to the 665 dition numbers (see The condition numbers are so high that even LU QUAD cannot provide accurate resolution. For the MoMaS medium and Fe-Cr mineral tests, many of the solutions calculated by the LU QUAD method are rescaled using Eq. 23, leading to 2)). The relative error on the norm increases regularly with 706 the condition number. It is close to 10 -16 when the con-707 dition number is close to 1 and increases to 1 when the 708 condition number is close to 10 16 , in accordance with the 709 computation theory presented by Golub and van Loan [START_REF] Golub | Matrix computations[END_REF]. 710 Q7 Fig. 4 Mean of the angles between reference and computed solution for each method and each test case excessively high relative error on the norm. For condition numbers greater than 10 16 , the evolution of 711

Comparison of the relative error on the norm given by the relative error on the norm with the condition number is 712 the non-preconditioned (LU, DGETRS, and GMRES) and much noisier. The use of the quadruple-precision LU QUAD 713 preconditioned (LU Equil, DGESVX, and GMRES Equil) method leads to an accurate resolution of a large portion 714 methods shows that the preconditioned methods lead to of the tested systems. As expected by computation theory, 715 lower relative error than the non-preconditioned methods all the systems with condition numbers less than 10 32 are 716 very high condition number tests (Appendix 4 (D-27, D-29, 757 test case. These points correspond to the rescaling of the and D-30)), the errors given by the preconditioned methods 758 computed quadruple-precision solution to maintain it on the are equivalent to the errors given by the non-preconditioned 759 double-precision scale (using Eq. ( 23)).

methods.

760

Iterative methods present similar behaviour to direct methods, giving very low relative error on the norm (between 10 -15 and 10 -8 ) when the condition number is (ii) By plotting the angle between the reference solu-761 tion and the calculated solution, we can compare the 762 methods according to the computed direction (Fig. 4). 763 less than a critical value. This critical value depends on the Because the resolution of the linear system (13) repre-764 Q8 743 method and the test case. step in the iterative Newton procedure, this 765 information is much more important than the norm 766of the step. A wrong norm can be corrected using 767 line search methods [START_REF]Chapter 8 Systems of nonlinear equations[END_REF], whereas modifying a wrong 768 direction leads to additional iterations. Small condition 769 number tests (gallic acid, Valocchi, pyrite, and pyrite 770 mineral) are solved using direct methods with the right direction. If the condition number increases, the directions given by the direct methods become inaccurate, but the condition number is not the only govern-ing parameter. Morel-Morgan leads to worse directionthan MoMaS medium and Fe-Cr, and MoMaS hard leads to a higher angle than the Fe-Cr mineral test. Iter-ative methods result in a worse direction than direct methods, and only the Valocchi test case is solved with an accurate direction by all the iterative meth-ods. Imposing lower convergence criteria (10 -12 ) on GMRES leads to a worse direction than using the usual criteria (10 -8 ). Using preconditioning methods leads to a better direction when associated with a direct method (LU Equil and DGESVX), but the conclusionis less clear for the iterative GMRES Equil method. Depending on the test case, the direction can be worse (Valocchi, MoMaS easy, MoMaS medium) or better (gallic acid, pyrite, MoMaS hard, Fe-Cr mineral)

The influence of the condition number on the angle (see Appendix 5) indicates that the direction is correct for direct methods when the condition number is less than 10 15 . For iterative methods, the limit to obtain an accurate direction is a condition number less than 10 8 , excepted for the Gauss-Seidel method, which produces wrong directions for low condition numbers. If the condition number increases, the behaviour of the direction becomes noisy. Since the relative error on the norm increases regularly until the condition number reaches the limit of 10 8 or 10 15 , the angle is accurately defined until this condition number limit is reached. Using preconditioned methods leads to a better direction for the LU Equil and the GMRES Equil methods when the condition number is higher than 10 15 for some cases (Appendix 5 (E-21, E-23 to E-25, and E-30)) but to a worse direction for other cases (Appendix 5 (E-26 and E-29)).

We present two successful direct methods, LU and LU QUAD; one iterative method, GMRES (both tested versions, GMRES and GMRES 10 -12 ); and two preconditioned methods, LU Equil and the GMRES Equil. By comparing the relative error on the norm (Appendix 4 (D-21 to more complicated. For all the tests cases (Appendix 5 (E-815 21 to E-25, E-27, E-28, and E-30)), LU QUAD gives the 816 best direction, followed by LU Equil, LU, GMRES Equil, 817 and GMRES 10 -12 . The MoMaS medium (Appendix 5 (E-818 26)) and MoMaS hard (Appendix 5 (E-29)) test cases lead 819 to the same conclusion, except GMRES Equil which gives 820 the worst direction. 821

Efficiency of the methods

822

The speed of the methods is studied by recording the com-823 putation time for each test case and plotting the mean CPU 824 time for each test case and each method (see Fig. 5). As 825 expected, the computation times are very short (less than 1 826 ms) because the systems to solve are small. 827 Figure 5 shows the influence of the system size. For all 828 methods, the computation time increases with the number 829 of unknowns. The results show that the iterative methods 830 are less sensitive to the system size than the direct meth-831 ods. For the iterative methods, the number of iterations is 832 important and depends on the first guess and other fac-833 Inversion method >10 30 LU QUAD Equil 10 30 ≥ cond(Z) >10 14 ). If the difficulty of the test increases, the lower accuracy of GMRES (compared to the quadruple-precision routine used in LU QUAD Equil, LU QUAD, and the combined algorithm) leads to a greater number of Newton iterations, as shown in Fig. 7 for the MoMaS hard case. This point is confirmed for other cases (see Appendix 7 (G-4 to G-9)). For the Fe-Cr mineral case (see Appendix 7 (G-10)), only LU QUAD Equil and the combined algorithm can solve the problem. Other methods lead to non-convergence, due to overflow for the GMRES algorithm (overflow appears in the Newton algorithm and is not due to GMRES itself) and because LU QUAD and LU are unable to give an accurate descent direction. Appendix 7 (G-11) shows the evolution of the rela-tion between the norm of Y and the condition number ofthe Jacobian matrix during the minimization process. This figure is similar to Fig. 1, confirming the empirical relation (25). This relation cannot be used close to the solution, and the condition number tends to be a case-dependent limit for very low |Y| .|| Nevertheless, the number of iterations is not the critical point. Because the time required by one iteration changes depending on the method used, we have to consider the total computation time. By plotting the total computation time required to solve each test case depending on the algorithm used (see Fig. 8), we can see that

)) that the number of required iterations can be twice the number for other methods, but we show (Fig. 5) that the GMRES method is faster than the other methods. (v) The proposed combined algorithm leads to intermediate computing times, equivalent to those of LU QUAD Equil and LU QUAD, depending on the case.

According to our results, GMRES should be systematically used because it is fast and usually leads to convergence of the Newton-Raphson algorithm. The combined algorithm should be used for very high condition numbers or for recomputing a failed run.

Conclusion

In this work, we focus on the resolution of small linear systems generated using the Newton-Raphson algorithm to solve equilibrium chemistry problems. For the first time, we propose a study of the condition number of such linear systems and find that the range of values covered is unusually large. This characteristic leads to specific numerical problems, with matrices that are quite small (approximately 10 ×10) but very badly conditioned (up to 10 100 ). Ten different chemical systems are studied.

There is a strong linear relationship between the logarithm of the condition number of the matrix and the logarithm of the norm of the objective function. This factor can be exploited to create efficient algorithms. This relation is strictly an empirical one and is not valuable for low condition numbers.

A wide variety of linear solvers have been tested, and several direct and iterative solvers are selected. Some of these solvers are specific for a class of matrix, symmetric or positive definite, while others are generic. A preconditioning method (matrix equilibration) has also been tested to reduce the conditioning of the systems.

According to our selected test cases, only the LU and LU QUAD direct methods, the GMRES iterative method, and LU Equil and GMRES Equil preconditioned methods are sufficiently robust to solve all the tests.

According to the size of the chemical tests, the LU method is faster than the GMRES method. However, our results for the Fe-Cr mineral and Morel-Morgan cases show that GMRES is preferable for larger chemical systems (more than 10 components). Chemical systems with more than 10 components have not been frequently modelled in the past decade. However, the use of geochemical databases makes the construction of large geochemical systems easier, and the increase in computation capacities makes it possible. For very large geochemical systems, we recommend the GMRES method.

We also propose using the linear relationship between the condition number of the Jacobian matrix and the norm of the objective function to develop an efficient algorithm.

The classic LU method is not a good choice. Its weakness is its low robustness for challenging test cases. We recommend using the GMRES method, which is fast and usually leads to convergence of the Newton-Raphson algorithm. For very high condition numbers (more than 10 100 ), we recommend the most robust LU QUAD Equil method. When the Newton-Raphson method is sufficiently near the solution to decrease the condition number, the faster GMRES method can be used. By using the linear relationship between cond(Z) and ||Y ||, the transition between the two methods can be achieved without computing the condition number (which is very expensive). as components whereas the pyrite case uses O2, Fe , 1107 and SO 2-. The first studies on geochemical computation 1108 avoided redox problems. We show that redox problems 1109 lead to higher condition numbers because the stoichiomet-1110 ric coefficients and equilibrium constants cover a wider 1111 range. Several geochemical databases avoid the introduc-1112 tion of redox reactions. There is sometimes a good reason to 1113 not write redox reactions as equilibria (slow reaction rates, 1114 irreversible reactions) as done in Arora et al. [2]. However, 1115 the reason is sometimes numeric, and redox reactions are 1116 avoided because they lead to non-convergence.

1117

We propose the use of quadruple-precision real for chal-1118 lenging chemical systems. In this work, the core of the 1119 geochemical code is conserved as double-precision real, 1120 and only the linear system tool is set as quadruple preci-1121 sion. Rewriting an entire geochemical code in a quadruple-1122 precision format will result in robust code but at the cost 1123 of an important and rebarbative work as well as efficiency. 1124 In this stage of our research, we do not recommend such 1125 an effort because implementing LU decomposition using 1126 quadruple-precision real is very efficient, requiring only 1127 a minor modification of existing code and reducing the 1128 computation time. ear solvers, but it was not possible to explore the flexibility of all the tested solvers. This study will help us to eliminate some solvers so that our future work can focus on the most promising: LU, LU QUAD, GMRES, LU Equil, and GMRES Equil. Some points for future exploration are as follows: (ii) The accuracy of iterative methods depends on the value of the convergence criterion (which we set to 10 -8 ) and on the method used to check the convergence (we used the default method). Moreover, the efficiency can vary depending on the initial guess provided by the user. In this work, we used the easiest initial guess: the residual for the tests from the Newton-Raphson method and the previous Newton-Raphson step for the test in a Newton-Raphson algorithm. We believe that it is possible to make a better choice, markedly enhancing the efficiency of the iterative methods.

(iii) The GMRES method allows the use of left and/or right preconditioners. These preconditioners can increase the robustness, accuracy, and efficiency of the method. More generally, several classes of preconditioners that may reduce the condition number of the linear system can be used [START_REF] Soleymani | A rapid numerical algorithm to compute matrix inversion[END_REF][START_REF] Soleymani | On a fast iterative method for approximate inverse of matrices[END_REF]. In this work, we explored the use of one preconditioner: matrix equilibration. However, other classes of preconditioners may be more efficient.

(iv) Previous works have addressed the use of methods to solve geochemical equilibria other than the Newton-Raphson method [START_REF] Carrayrou | New efficient algorithm for 1274 solving thermodynamic chemistry[END_REF]44,49,[START_REF] Morin | Simplified explanations and examples of computerized methods for calculating chemical equilibrium in water[END_REF]. It has been shown [START_REF] Carrayrou | New efficient algorithm for 1274 solving thermodynamic chemistry[END_REF] that an efficient algorithm can be obtained by combining a zero-order method with the Newton-Raphson approach.

(v) The size of the chemical tests presented here is representative of the sizes actually used in environmental studies. We have shown that the GMRES method may be efficient for large systems. In anticipation of future needs, it may be useful to test chemical systems larger than the Morel-Morgan system.

(vi) Part of the Newton minimization related to very large condition numbers (far from the solution) can be performed using random methods; GMRES is efficient even though its descent direction is not accurate

  for the computation of thermody-147 namic equilibrium is based on the tableau concept, referred 148 to as Morel's table[START_REF] Morel | A numerical method for computing equi-176 (2001) 1392 libria in aqueous chemical systems[END_REF][START_REF] Morel | Principles of aquatic chemistry[END_REF]. NX components (Xj ) are cho-149 sen from the NC species (Ci) and are used to write the 150 formation of each species as a combination of the com-151 ponents. The mass action law for the formation of the Ci 152 species is written with the equilibrium constant (Ki) and the 153 stoichiometric coefficients (ai,k) for each component (Xk) 154 . Nx {Ci } = Ki {Xk } ai,k (1) k=1 where {Ci } and {Xk } are the activities of species Ci and 155 component Xk, respectively. In this work, we define Xj as a 156 subset of Ci; then, NX is NC minus the number of reactions.j j

  of the objective function describing precipitation, 210 it is more interesting to rewrite the mass action law (2) in 211 log form and then define the objective function Nx YNx+i = ln (SIi) =ln (Kpi) + api,k • ξkfor i = 1to NCP k=1 (9) 212 Equations 8 and 9 are solved at the nth iteration with the 213 Jacobian matrix (Z n ) of the objective functions . Z n . j = 1, N + N ∂Y n = [ n (i) The gallic acid test case was presented by Bras-240sard and Bodurtha [49]. It has been recognized as a

243 j,k j = 1 ,

 1 NX + NCP ∂ Cpk-N n (ii) The Valocchi test is from Valocchi et al. [11]. It X 214 215 216 k = NX + 1, NX + NCP Z n can be calculated in two ways. (i) Using an analytical computation, we obtain the (NX + NCP ) ×(NX + NCP ) values of Z n by involves calcium and magnesium ion exchange (see 245 Appendix 1 (2)).

  246

1 (

 1 NX k = NX + 1, NX + NCP . j = NX + 1, NX + NCP = ap k,j -NX (11) for the reactive transport benchmark of MoMaS at the easy level [50]. It has been specifically developed to magnify numerical difficulties in a small system n . Z j = N k = 1, NX + see Appendix 1 (4)). (v) The Morel-Morgan test is the first large chemical 217 j,k 229 X , NX + NCP = 0 k = NX + 1, NX + NCP Even if the activity coefficients depend on the compo-nent 218 230 concentrations, they are assumed to be constant during the Newton-219 Raphson procedure. These activity coefficientsare usually actualized by a 220 fixed-point algorithm at each Newton-Raphson loop. 221 The progress step of the method (∆ξ n , ∆Cp n ) is achieved by assuming 222 that the objective function Y n+1 inEq. 12 is equal to zero at the (n 1)th 223 + iteration. This pro-duces the key equation of this article, the linear system

  after 50 iterations. Let R and C This procedure is repeated until all ¨rk¨ and ¨ck¨ are DPOTRS subroutine [61] based on the Cholesky method. 430

450 10 -

 10 12 as the stopping criterion, denoted byGMRES 10 -12 451 in this study. A critical point of the GMRES algorithm is the 452 size of the Hessenberg matrix. In this work, we set it to the max of 8 (Nx + NcP).

  526chi, pyrite, and Morel-Morgan test cases, all the Jacobian 527 matrices are positive definite. For the MoMaS easy, MoMaS 528medium, and Fe-Cr test cases, a large proportion (66.4 to 529 74.1 %) of the Jacobian matrices are positive definite. For 530 cases including minerals (pyrite mineral, MoMaS hard, and 531 Fe-Cr), essentially none of the matrices are positive defi-532 nite (only 0.1 % for the MoMaS hard test). Plotting the ratio 533of positive definite matrices as a function of the condition 534 number (see Appendix 2 (B-2)) shows that the chemical 535 conditions are more important than the condition number 536 when determining whether the Jacobian matrix is diagonal 537 dominant.538 (v) According to the test procedure presented previously, 539 we plot, on the same graph, the logarithm of the norm 540 and the results are noisier if cond(Z) and ||Y || are ⎡ 1.15 • 10 94 9.09 • 10 93 -5.04 • 10 -13 -|a |nd the logarithm of the condition number of the matrix Z (Fig. 1). There is a strong linear relation-ship small. The evolution of this relation for low |Y | ||can be seen Appendix 7 (G-11). Therefore, Eq. 25 should not be used 555in 556for between these parameters. Moreover, the linear Y ||les||s than 10 10 . 557 relationship does not depend on the chemical test, only Several of these properties are obtained using the loga-558

⎟Fig. 1

 1 Fig. 1 Relationship between the condition number of Z and the norm of the objective function plotted on a log-log scale

653 7 )

 7 ) show that for condition numbers greater than 10 20 , the failure ratio increases greatly for some of the methods. These methods are DOPTRS and DSYTRS for MoMaS medium and Fe-Cr. DGETRS, UMF, DSYTRS, DOPTRS, (i) By plotting the mean of the logs of the relative error 654on

Fig. 5

 5 Fig. 5 CPU time (s) for each method and each test case

  It can be set to 10 8 for SYMMBKCGS, BiCGStab, and GMRES for the gallic acid (Appendix4 (D-11)) and MoMaS easy (Appendix 4 (D-14)) cases.It can be set to 10 12 or 10 15 for Inc. CHOLESKY forthe gallic acid and MoMaS easy cases and for SYMMBK, Inc. CHOLESKY, CGS, BiCGStab, and

F

  do ig u u b r l e e-5 pr a e l c s i o si s o h n o r w e s al th to e q co ua m d p r u p ti l n e g -p t r i e m c e isi r o e n qu r i e r a e l d a t n o d o b b a ta c i k n . 836 837 the reference solution, which requires more time. 838 The UMF method is the slowest double-precision direct 839 method, but its multifrontal block strategy becomes interest-840 t i e n s g t f r o e r qu la ire g s e 3 s 3 ys t t i e m m es . m T o h r e e r C es P o U lu t t i i m on e o th f a t n he th M e r o e r s e o ll M uti o o r n ga o n f 841 842 the MoMaS easy test for the UMF method, whereas it takes 843 190 times more time for the LU method. 844 Sei A de m l o m ng eth th o e d it a e n r d ati t v h e e m sl e o th w o e d s s t , i t s he th f e as I t n e c st . i C s H th O e L G E a S u K s Y s-845 846 method. The two most robust iterative methods, BiCGStab 847 and GMRES, are rapid, sometimes more so than the direct 848 robust methods, LU and UMF, especially for large systems 849 (Morel-Morgan test case). GMRES is less case-dependent 850 than BiCGStab, leading to similar computing time, regard-10 cond(Z) >10 LU Nx NcP <10 GMRES Nx NcP 10 10 4 >cond(Z) LU D-30)), the successful methods can be ranked from the lowest to highest error: LU QUAD, GMRES 10 -12 , GMRES Equil, LU Equil, and LU. Ranking these methods according to the angle between the reference and computed solution is

Fig. 7

 7 Fig. 7 Evolution of NRrelativeerror as a function of the Newton-Raphson iteration for the MoMaS hard case

Fig. 8

 8 Fig. 8 Computation time (s) as a function of test case and algorithm

  (i) LU QUAD Equil, as expected, is the slowest. Nevertheless, this method allows the convergence of the Newton-Raphson method for all test cases. (ii) LU QUAD is slightly faster. The difference between LU QUAD Equil and LU QUAD gives an indica-tion of the time used for matrix equilibration. This time is greater for pyrite, MoMaS easy, pyrite mineral, MoMaS hard, and Fe-Cr mineral than for the other test cases. (iii) LU is fast when it leads to convergence, but this method results in a very weak Newton-Raphson algorithm. (iv) GMRES always results in the fastest Newton-Raphson algorithm. It has been shown (Fig. 7, Appendix 7 (G-

  a new research field by studying geochemical computation from a condition number point of view. We attempted to benchmark a wide variety of lin-for high condition numbers. Some methods, such as 1095 simulated annealing and particle swarm optimization, 1096 could be used in future research. 1097 These factors should be explored in light of the results pre-1098 sented in this study. We proposed a large set of chemical 1099 tests, a criterion to determine the difficulty of these tests (the 1100 condition number), and a panel of numerical methods that 1101 should be studied preferentially. 1102 As a more general consideration, the reader should pay 1103 particular attention to the old Morel-Morgan test case and 1104 the more realistic pyrite test case. The Morel-Morgan test 1105 uses Fe 2+ and Fe 3+ , Co 2+ and Co 3+ , and SO 2-and S 2-1106 4 2+

1129

  Acknowledgments Hela Machat has been supported by a grant from 1130 the Tunisian Government. This work is supported by the BRGM-1131 CNRS CUBICM project. We thank the reviewers for their helpful This article does not contain any studies with 1138 human participants or animals performed by any of the authors. 1139 Informed consent Informed consent was obtained from all individ-Q10 1140 ual participants included in the study.

  1141

  (i) We did not extensively test the robustness and the efficiency of the Newton-Raphson algorithm. Further work should examine the influence of the initial Newton-Raphson guess to confirm our conclusions about the high efficiency of the GMRES method.

  . A second 130 part is devoted to the presentation of the chemical tests and 131

	w th e p n r u o m po er s i e c a d p e r t o a c il e e d d u a re n s al u y s s e is d o t f o t p he s fo tr r u m ctu th re a te n s d ts p . r N op e e x r t -,
		133
	th e e s n o p f re th sen J t a e c d o a b n ia d n te m st a e t d ri , x a . n T d h t e he se r l e e s c u t l e t d s a li r n e e c a o r m so p l a v red s a a n r d e	134
	Q2

132 chemical tests leading to a 7 ×7 matrix with a condition number of approximately 10 180 .

  involvesthe resolution 158 If NCP -precipitated species (Cpi) are taken into account, between the activity and concentration are given by activity 172 159 the mass action law for the precipitation of Cpi is written coefficients (γi) calculated using specific models (Davies, 173 160 with the precipitation constant (Kpi) and the stoichiomet-Debye-Hu ¨ckel, etc.)

	of	(see	Section
	the	3.1).
	syste
	m (6)
	with
	the
	Newt
	on-
	Raph
	son
	meth
	od
	using
	[Xk]
	and
	[Cpi]
	as
	prim
	ary
	unkn
	owns
	.
	This

202

  For each chemical system, we Many numerical methods, including a failure indicator, 335 295 select three components and vary their values over a wide which indicates the success or failure of the resolution, have 336 296 range. The concentrations of all minerals are arbitrarily set been developed. If needed, we include a failure indicator. 337 297 to 10 -3 mol L -1 . The activity of component H + is varied

	inclu ding 282 mine 283 ral 284 Cp if prese nt.	I s t yste 2 m 61 reported in the computational literature. was used by F. Morel and M. Morgan in 1972 to 262 present the capacities of the (ix) The MoMaS hard test is the equilibrium part of the computational method 263 they had just chemical system described in the hard level of the developed (and which we still use 264 MoMaS reactive transport benchmark. It allows for today). This test includes 52 components (H + , 20 265 metals, and 31 ligands), leading to 781 aqueous the formation of two mineral species (see Appendix 1 (9)).	Comput Geosci The condition number of Z is defined [23] as the product of 323 the norm of the matrix per the norm of the inverse matrix 324 (17) 325 cond (Z) = ǁZǁ1 × ¨Z-1 ¨ (17)
		266		1
		species (x) The Fe-Cr mineral test describes the redox reaction (see Appendix 1 (5)). between iron and chromium. We assume the forma-	Appendix 1 presents the stoichiometric Q5 coefficients, equilibrium constants, and concentrations for To test the numerical methods, we first evaluate the compu-326 tation time (CPU time) required to solve the linear system. 327
		tion of three different mineral phases (see Appendix 267 1 (10)). (vi) The MoMaS medium test is the chemical system for 268	these tests. Because we work with a very small matrix, the computa-328 tions are very fast and we run the same calculation several 329 times to obtain a total computing time of approximately 1 330
	the medium level of the MoMaS reactive transport 291 2.4 Test procedure	s. The CPU time is given in this work in units of seconds 331
		269		per computation (by dividing the total computing time by 332
	benchmark 292 Equation 11 shows that we can obtain multiple linear sys-[50] (see Appendix 1	the number of runs). According to this method, the global 333
	(6)). 293 tems from one chemical problem by changing the activity	computing time for one test case is approximately 6 days.	334
	270 294 values of the components.
		(vii) The Fe-Cr test is an additional redox test that
		271
		describes the redox reactions between iron and
		272
		chromium. These types of reactions occur when
		273
		iron reactive barriers are used to treat
		chromium-274 contaminated sites [51, 52].
		In this case, we con-275 sider only the
		aqueous phase without minerals (see 276
		Appendix	1
		(7)).
		277
		(viii) The pyrite mineral test describes the dissolution
		of 278
		a pyrite rock in pure water. We assume that three
			27
			9
		possible mineral phases are present (see Appendix 1
			28
			0
			(8)).
		1 RT	{O } H

  All of these quantities, namely the failure indicator, relative 355error on the norm, anglemethod, and CPU time, are calculated 356for the 29,791 linear systems built from each chemical test 357 case for all the tested methods. This enormous amount of 35|8|-da || ta

	by				Comput Geosci
	Eqs.
	8 and
	9).
	T
	he	is	aggregated	in	two
	matri	ways:
	x
	norm	359
	used	(i) For each chemical system and each method, we
	in	com-
	this	360
	work	pute the mean of each
	is the	quantity.
	1 norm , defin	361 (ii) For each chemical system and each method, the inter-362
	ed as	val of the condition number is discretized into
	[23]	100
		363
	315	319
	316	320
	317	321
	318	322

  because it preserves the symmetry of the reports the use of iterative methods for solving large, sparse 412

			JrnlID 10596 ArtID 9600 Proof#1 -04/11/2016
	i a e Q t t e i M r v s [ h e u d 3 a s 2 e v ] R G a to b C St s e olv 431 reac tive tran	Comput Geosci st the Jacobi [23, 434 62], Gauss-Seidel [23, 62], and successive over-relaxation
	spo rt und	435 (SOR) [23, 62] methods. Barrett et al. [21] proposed
	glo er a	Jacobian matrix. an	linear systems [20-22, 24, 26, 28, 29, 56, 57]. Nevertheless, 413
	bal	436
	app	algorithm to select an iterative solver depending on
	roa	the
	ch.
	Her	437
	e,	matrix properties. GMRES was presented as the least
	we test	selec-
	QR	438
	432 dec	tive algorithm. We use a GMRES method developed by
	om posi tion usin	439 HSL [63]. If the matrix is symmetric, Barrett et al. [21] rec-
	g	440
	the DG EL	ommend the use of conjugate gradient squared (CGS) or
	S	441
	rout	biconjugate gradient stabilized (BiCGStab) methods.
	ine	CGS
	[61]
	.	442
		and BiCGStab subroutines have been developed by
	433 F o r t h	HSL. 443 We test two additional methods devoted to symmetric matri-444ces: SYMMBK [63] and an incomplete Cholesky (Inc. 445 CHOLESKY) factorization [63].
	e
	it
	e
	r
	a
	ti
	v
	e
	m
	e
	t
	h
	o
	d
	s
	,
	w
	e
	t
	e

446

Table 1

 1 List

	of the selected
	solvers	Name	Source Method	Matrix properties
		Direct
		LU	[58]	LU decomposition	-
		DGETRS	[59]	LU decomposition	-
		UMF	[53]	Direct multifrontal	-
		DSYTRS	[59]	UDU-factored symmetric matrix	Symmetric
		DPOTRS	[59]	Cholesky A = U T ×U	Definite positive
		DGELS	[59]	QR decomposition
		LU QUAD	[58]	LU decomposition quadruple precision	-
		Iterative
		SYMMBK	[61]	Iterative SYMMBK HLS MI02	Symmetric
		Inc. CHOLESKY [61]	Incomplete Cholesky HSL MI28	Symmetric
		CGS	[61]	Conjugate gradient squared HLS MI23	-
		BiCGStab	[61]	Biconjugate gradient squared stabilized HLS MI26 -
		GMRES	[61]	Flexible GMRES HLS MI15	-
		Gauss-Seidel	[58]	Gauss-Seidel method	-
		Preconditioned
		LU Equil	[51-58] LU and matrix equilibration	-
		DGESVX	[59]	LU and optional preconditioning	-
		GMRES Equil	[51-61] GMRES and matrix equilibration	-
		GMRES 1.d-15 [61]	GMRES convergence criteria 1.d-15	-
	459 method [23, 26, 56, 62]

Table 2

 2 Structure of the Jacobian matrix

	JrnlID 10596 ArtID 9600 Proof#1 -04/11/2016
	Comput Geosci

Table 3

 3 Properties of the 10 chemical test cases ranked by increasing the maximal condition number Some matrices in the gallic acid and 509 as quadruple precision, denoted by LU QUAD. The source Valocchi cases are diagonal dominant, but none of the 510 of this method is the LU double-precision real of numer-matrices from the other cases are diagonal dominant. 511 ical recipes [60], and we adapt it to quadruple precision.

		Nx	Nc	NcP	Z size	cond(Z) min	cond(Z) max	cond(Z)	%Z diag-	%Z positive
								max after	onal domi-	definite
								20 equili-	nant
								bration
	Gallic acid	3	17	0	3	10 0.61	10 12.6	10 0.95	18.4	100
	Valocchi	5	7	0	5	10 0.49	10 15.3	10 0.65	67.7	100
	Pyrite	4	40	0	4	10 4.06	10 24.9	10 0.95	0.00	100
	MoMaS easy	5	12	0	5	10 3.44	10 37.7	10 1.05	0.00	71.1
	Morel-Morgan	52	781	0	52	10 43.4	10 60.7	10 1.13	0.00	35.9
	MoMaS medium	5	14	0	5	10 5.88	10 103.9	10 0.95	0.00	78.8
	Fe-Cr	7	39	0	7	10 9.46	10 113.6	10 1.05	0.00	68.9
	Pyrite mineral	4	43	3	7	10 1.71	10 33.1	10 3.19	0.00	0.00
	MoMaS hard	6	15	2	8	10 5.45	10 123.9	10 3.02	0.00	0.00
	Fe-Cr mineral	7	43	3	10	10 8.67	10 213.9	10 13.4	0.00	0.00

Finally, we test an LU decomposition method compiled test procedure.

Table 3

 3 methods, whereas the first three tests are accu-for the direct methods, but the result is more case-dependent 692 rately solved by all the direct methods. Even in the case for GMRES. The use of preconditioning usually leads to 693 of successful resolution (CGS and BiCGStab methods), lower relative error on the norm, except for the the relative errors on the norm are high for intermediate Morgan, Fe-Cr, and MoMaS hard cases. MoMaS easy, and Morel-Morgan). Nev-Increasing the maximum number of iterations and reduc-696 ertheless, the results are better for the iterative meth-ing the convergence criteria of GMRES leads to less relative 697 ods than for the direct methods for the difficult tests error on the norm, but this reduction is not significant. Seidel are less efficient than the other methods method depending on the condition number. The theoreti-703for the easy tests, but more ill-conditioned tests are cal behaviour is verified for the direct methods and for all 704 better solved by these two methods. the test cases (except for the Valocchi one, Appendix 4 (D-705

	698

) are too high to exhibit any condition number than direct methods. Only the Val-666 ordering. occhi test case is accurately solved by almost all the 667 Fig. 3 Mean of the logs of the relative error on norm for each method and each test case UNCORREC iterative 695 cases (pyrite,

Table 4

 4 Algorithm for equilibrium computation cond(Z)

  LU QUAD14 ≥

	less	Comput Geosci
	of
	the
	test
	cas
	e. 85
	2
	A
	s				(LU Equil, DGESVX, and GMRES Equil) or decreasing
	e				the convergence criteria for an iterative method (GMRES
	x p e c t e d , i				10 -12 ) leads to increased computing time. The computing time for preconditioning does not depend only on the sys-tem's size: the Valocchi, MoMaS easy, and MoMaS medium test cases (system size of 5 ×5) are solved with the same computing time for all the direct methods, but their resolu-tion when using LU QUAD Equil, LU Equil, and GMRES Equil is faster. Appendix 6 shows the computation time (log scale) for Q9
	n				each test case and each method depending on the condition
	t				number. Appendix 6 (F-1 to F-10) shows that, as expected,
	r
	o
	d
	u
	c
	i
	n
	g
	p
	r
	e
	c
	o
	n
	d
	it
	i
	o
	n
	i
	n
	g
	t
	e
	c
	h
	n
	i
	q
	u
	e
	s
	8
	5
	3
		4	+	+	≥

OOF

Fig. [START_REF] Andre | Numerical modeling of fluid-rock chemical 1234 interactions at the supercritical CO2-liquid interface during CO2 1235 injection into a carbonate reservoir, the Dogger aquifer (Paris 1236 Basin, France)[END_REF] Evolution of NRrelativeerror as a function of the Newton-Raphson iteration for the pyrite test case 866 the computation time of the direct methods does not depend number is greater than the critical value (see Appendix 5). 867 on the condition number of the system. The LU method is

The critical condition number is 10 8 for GMRES, 10 16 for 868 usually 10 times faster than the UMF method, except for the the double-precision direct methods, 10 32 for LU QUAD, 869 Morel-Morgan test case, in which LU is only 1.5 times faster. and case-dependent for preconditioned methods (10 20 to 870

In Appendix 6 (F-11 to F-20), the general tendency for 10 60 ). Gauss-Seidel leads to wrong directions for very low 871 the iterative methods is to require the same computation condition numbers (Appendix 5 (E-11 and E-12)).

872 time, independent of the condition number. The oscilla-In terms of efficiency, the most rapid method is Gauss-873 tions presented by the curves seem to be not related to Seidel when it is available. The second most efficient 874 the condition number. For the test case without minerals, method is LU for small systems (less than 10 ×10) or 875 the Gauss-Seidel method is efficient. The two most robust GMRES for larger systems (more than 10 ×10), and the 876 methods, BiCGStab and GMRES, are often the third and slowest method is LU QUAD. For small systems (less than 877 fourth fastest methods (Gauss-Seidel and SYMMBK are the 5 ×5), LU Equil is as fast as GMRES but becomes slower as 878 fastest).

the system size increases.

We recommend using LU, LU QUAD, GMRES, and the reference method LU QUAD Equil. Gauss-Seidel should be 879 4 Proposal of a new algorithm rejected because of its wrong direction, and equilibration does not sufficiently improve the behaviour of double-880 Based on our results, we propose an algorithm to opti-precision routines.

881 mize the resolution of a chemical system using a Newton-Using Eq. 25, it is possible to estimate the condition 882 Raphson-like method. number of the system without additional computation. This

883

Examining the failure ratio results, seven methods are estimation enables the selection of the best-adapted method 884 eligible: LU and LU QUAD as direct methods, GMRES depending on the system size and condition number.

885 and Gauss-Seidel (if no minerals) as iterative methods, LU The goal is to use the most robust method (LU QUAD 886 Equil and GMRES Equil as preconditioned methods, and with preconditioning) for high condition number systems 887 the reference method (LU QUAD Equil).

(more than 10 32 ) in the first Newton-Raphson iterations.

888 Because these methods are included in a Newton mini-When the condition number is sufficiently decreased, the 889 mization procedure, the most important accuracy criterion preconditioning becomes useless and LU QUAD can be 890 is the direction of the minimization, i.e. the angle between used until the condition number is less than 10 16 . Then, a 891 the reference and the calculated solution. The behaviour of faster method is used to obtain a coarse approximation of 892 this direction is strongly correlated with the condition num-the solution, LU for small systems and GMRES for large 893 ber of the system and is correct if the condition number systems (more than 10 ×10). To find the exact solution, the 894 is less than the critical value and wrong if the condition LU direct method is used. Q1. Please check the provided city (Jendouba) and country (Tunisia) names in affiliation 2 if correct. Q2. Please check if the abbreviation "LAPACK" is defined correctly. Otherwise, please provide the correct expansion. Q3. Please check all equation citations if correctly captured or presented. Q4. The section citation "Properties of the Jacobian matrix" found in the sentence starting "This formulation has some weaknesses" was changed to "Section 3.1" as per numbering style. Please check if correct. Q5. As per instruction, if there are more than one appendix, they should be designated with numbers 1, 2, 3, etc. With this regard, the appendices were renumbered accordingly and their citations were modified. Please advise if the labels designated for figures found in appendices B to G (as originally labelled) should also be modified accordingly to ensure consistency of presentation. Q6. Missing citation for Table 1 was inserted in this sentence. Please check if appropriate. Otherwise, please provide the location of where to insert the citation/s in the main body of the text. Note that the order of main citations of tables in the text must be sequential. Q7. Please check the captured caption of Figure 4 if correct. Q8. Missing citation for Figure 4 was inserted in this sentence. Please check if appropriate. Otherwise, please provide the location of where to insert the citation/s in the main body of the text. Note that the order of main citations of figures in the text must be sequential. Q9. The appendix citations found in the last two paragraphs before the "Proposal of a new algorithm" section were all changed to "Appendix 6" so that appendix citations are in sequential order. Also, the discussion in last two paragraphs fits to the description of Appendix 6. Please check if correct and amend if necessary. Q10. Please confirm if the "Informed consent" statement is indeed applicable to this article and should be retained. Q11. Please check if the content of Appendices 1 to 7 is correctly captured or presented. Q12. Please check if the labels "D-31 to D-40" are indeed correct in the artwork of Appendix 4, and not "D-21 to D-30", respectively.


Lama Khalil

Jean-Christophe Girard

Debora Pierucci

Federico Bisti

Julien Chaste

Fabrice Oehler

Charlie Gréboval

Ulrich Nguétchuissi Noumbé

Jean-Francois Dayen

Demetrio Logoteta

Gilles Patriarche

Julien Rault

Marco Pala

Emmanuel Lhuillier

Abdelkarim Ouerghi
email: abdelkarim.ouerghi@c2n.upsaclay.fr

Electronic band gap of van der Waals α-As2Te3 crystals

Keywords: As2Te3, Electronic Band Gap, Electronic Transport, Field Effect Transistor

van der Waals materials offer a large variety of electronic properties depending on chemical composition, number of layers, and stacking order. Among them, As2Te3 has attracted attention due to the promise of outstanding electronic properties, and high photo-response. Precise experimental determinations of the electronic properties of As2Te3 are yet sorely needed for better understanding of potential properties and device applications. Here, we study the structural and electronic properties of α-As2Te3. Scanning transmission electron microscopy coupled to energy X-ray dispersion (STEM-EDX), and micro-Raman spectroscopy all confirm that our specimens correspond to α-As2Te3. Scanning tunneling spectroscopy (STS) at 4.2K demonstrates that α-As2Te3 exhibits an electronic band gap of about 0.4 eV. The material can be exfoliated, revealing the (100) anisotropic surface. Transport measurements on a thick exfoliated sample (bulk-like) confirm the STS results. These findings allows for a deeper understanding of the As2Te3 electronic properties, underlying the potential of V-VI semiconductors for electronic and photonic technologies.

Following the successful isolation of monolayer (ML) graphene, many other two-dimensional (2D) van der Waals lamellar (vdW) compounds have been investigated. "Honeycombs hexagonal" ML materials, such as transition metal dichalcogenides (TMDs) [START_REF] Liu | Growth of Large-Area and Highly Crystalline MoS 2 Thin Layers on Insulating Substrates[END_REF][START_REF] Ambrosi | 2H → 1T Phase Transition and Hydrogen Evolution Activity of MoS2 , MoSe2 , WS2 and WSe2 Strongly Depends on the MX2 Composition[END_REF][START_REF] Chang | Orbital Analysis of Electronic Structure and Phonon Dispersion in MoS 2, MoSe2, WS2, and WSe2 Monolayers under Strain[END_REF] , metal mono-chalcogenides (GaSe, InSe, GaTe) [START_REF] Cao | Tunable Magnetism and Half-Metallicity in Hole-Doped Monolayer GaSe[END_REF][START_REF] Ben Aziza | Der Waals Epitaxy of GaSe/Graphene Heterostructure: Electronic and Interfacial Properties[END_REF][START_REF] Si | Controllable Schottky Barrier in GaSe/Graphene Heterostructure: The Role of Interface Dipole[END_REF][START_REF] Balakrishnan | Epitaxial Growth of γ-InSe and α, β, and γ-In2Se3on ϵ-GaSe[END_REF][START_REF] Aziza | Tunable Quasiparticle Band Gap in Few Layer GaSe/Graphene Van Der Waals Heterostructures[END_REF] and Xenes (X = Si, Ge, Sn, P, As, Sb, Bi, Te) [START_REF] Lei | Boron Carbon Nitride Nanostructures from Salt Melts: Tunable Water-Soluble Phosphors[END_REF] have attracted considerable attention from both physics and chemistry communities. Similarly, "buckled" ML from monoelemental group-V and -VI compounds (phosphorene, arsenene, tellurium, antimonene and bismuthene) are particularly interesting, due to high stretchability of the buckled structure [START_REF] Lei | Boron Carbon Nitride Nanostructures from Salt Melts: Tunable Water-Soluble Phosphors[END_REF] and the related strain modification on their band. Beyond these well-known examples, other materials crystallize in a vdW lamellar stacking, with slightly more complex ML crystal structure. The α-As2Te3 shows a "buckled"-like ML structure, and display a rather unique combination of properties, which are of interest for fundamental studies of anisotropic phenomena in 2D materials, but could also enable new types of applications such as acousto-optic modulators.

Despite the general interest in layered semiconductors, many optical and electronic properties of crystalline As2Te3 are unknown [START_REF] Morin | Polymorphism in Thermoelectric As2Te3[END_REF] . As2Te3 exist in two different crystal structures, αand β-, originating from distinct stacking sequences, which share the same composition. This polymorphism likely affects the electronic and physical properties of As2Te3. Yet, there is no report on the band structure properties of crystalline phases of α and -As2Te3 even at room conditions. The α-As2Te3 phase display a monoclinic structure [START_REF] Carron | The Crystal Structure and Powder Data for Arsenic Telluride[END_REF] which is composed of repeats of staircase pattern (4 As atoms sandwiched by 6 Te atoms) (Figure 1a)). Along the crystal axis a, these staircase units are simply stacked by van der Waals (vdw) interactions, creating an easy cleavage (100) plane. Inside each ML, the structure is very asymmetric with regular infinite chains along b and weak covalent bounding between staircase units along c. The α-As2Te3 ML thus displays a strong in plane anisotropy and 1D-like behavior. Several preceding investigations have focused on the phase transition induced in α-As2Te3 [START_REF] Morin | Polymorphism in Thermoelectric As2Te3[END_REF] , but no one has determined experimentally any basic electronic parameters, such as the electronic band gap.

We report here a STM/STS and transport study of the electronic properties of monolayer As2Te3. This study was also supported by micro-Raman spectroscopy and STEM for morphological and elemental assessments. The STM images demonstrate a large atomically ordered flat region with ligne defect. The Fourier transform indicates a spatial periodicity of 4.0 Å and 9.9 Å in perpendicular directions, which is consistent with the lattice parameters the monoclinic structure determined by Raman spectroscopy and STEM. This confirms that the exposed cleaved surface is the (100), as expected from the weak vdW stacking along the [100] direction.

Our sample is a commercial α-As2Te3 crystal from HQ Graphene (http://www.hqgraphene.com/). To determine the crystal structure of our As2Te3 samples, we conducted a comprehensive characterization by employing complementary techniques. The µ-Raman measurements were conducted using a commercial confocal HORIBA LabRAM HR operating at 532 nm [START_REF] Zribi | Strong Interlayer Hybridization in the Aligned SnS<inf>2</Inf>/WSe<inf>2</Inf> Hetero-Bilayer Structure. npj 2[END_REF] . The incident photon beam was focused down to a submicrometric spot (~1 µm in diameter) on the sample. Concerning the vibrational properties of α-As2Te3, 30 zone-center phonons belong to the Γ = 10(Ag + Bu) +5(Au + Bg) mechanical representation at the center of the BZ. Among these phonon modes, the gerade (g) modes are Raman active, while the ungerade (u) modes are infrared active, excluding 1 Au and 2 Bu modes, which are acoustic modes [START_REF] Cuenca-Gotor | Vibrational, and Electronic Study of α-As2Te3 under Compression[END_REF] . Thus, α-As2Te3 has in total 15 phonon modes that are Raman active. Note that the Au and Bg modes involve atomic displacements along the b axis, whereas the Ag and Bu modes are relative to vibrations in the a-c atomic plane. Figure 1(b) displays the typical µ-Raman scattering spectrum of our samples in the frequency range of 22-215 cm -1 . The µ-Raman spectra show 10 out of the 15 Raman active mode peaks, together with the position of the center of the fitting components, indicated with red markers; the fitting was performed using a Voigt function. In particular, we notice peaks at 43.1 cm -1 , 48 cm -1 , 64 cm -1 , 67.3 cm -1 , 88.9 cm - 1 , 97.4 cm -1 , 120.6 cm -1 , 137.4 cm -1 , 171.5 cm -1 and 193.1 cm -1 , corresponding to the Ag 1 , Ag 2 , Bg 2 , Ag 3 , Ag 4 , Bg 3 , Ag 6 , Bg 5 , Ag 8 and Ag 10 modes, respectively. The assignment of these phonon modes was carried out in accordance with the work performed by Cuenca-Gotor et al. [START_REF] Cuenca-Gotor | Vibrational, and Electronic Study of α-As2Te3 under Compression[END_REF] , and is coherent with the α-phase of As2Te3. Figure 2 shows a transmission electron microscopy (TEM, FEI Titan Themis operating at 200 kV) analysis of a cleaved bulk As2Te3 crystal. By mechanical cleaving, we were able to obtain a thin enough sample so that standard scanning TEM (STEM) coupled to chemical analysis, using energy dispersive X-ray spectroscopy (EDX, Brucker super X), could be performed. A large view of the probed area is shown Fig. 2a, acquired in the STEM mode using the high angle annular dark field detector (HAADF). The complementary EDX analysis reveals that the sample is homogeneously composed of As (Fig. 2b) and Te (Fig. 2c). EDX quantification indicates an atomic ratio of 2:3, as expected from the As2Te3 stoichiometry. Standard HR-TEM shows fringes at the edge of the sample (Fig. 2d), which indicates that the material is crystalline. Further analysis by HR-STEM (Fig. 2f) is obtained by zooming in the probed region (Fig. 2d ). We measure a spatial periodicity of 4.0 Å and 9.9 Å in perpendicular directions, which is consistent with the lattice parameters b and c of the monoclinic structure previously determined by Raman spectroscopy and STEM (Fig. 1 and2). This confirms that the exposed surface is (100), as expected from the lamellar vdW crystal structure. Then, we probed the electronic structure of α-As2Te3 by means of STS. A typical spectrum of the differential conductance (dI/dV) measured as a function of the bias voltage is shown in Fig. 3(c). The I(V) characteristics were acquired while the feedback loop was inactive, and the differential conductivity dI/dV was measured directly by using the lock-in technique. The differential conductance (dI/dV) is proportional to the local occupied and unoccupied density of states below and above EF, As the edges of the valence and conduction bands are observed to be around -0.032 eV and 0.368 eV, respectively, the electronic quasiparticle bulk bandgap is hence of 0.4 eV.

The Fermi energy, corresponding to the zero bias in the dI/dV spectrum, is positioned slightly below the mid-gap value, suggesting the presence of a small concentration of intrinsic charged defects in our specimens and slight ptype doping. In summary, we report an experimental study of electronic and transport properties of α-As2Te3. The material can be exfoliated along the weak [100] direction, exposing the asymmetric (100) surface. The latter is strongly anisotropic, with a buckled pattern (c direction) assembled in long parallel chains (b direction). Scanning tunneling spectroscopy and electronic transport measurements show that bulk α-As2Te3 is a semiconductor with a band gap about 0.4 eV, with potential applications in electronics or opto-electronics.

Figure 1 :

 1 Figure 1: Crystallographic and µ-Raman analysis of α-As2Te3: a) Crystal structure of α-As2Te3. Arsenic and tellurium atoms are represented in purple and orange spheres, respectively, b) µ-Raman spectrum acquired at room temperature from the α-As2Te3 crystal, showing 10 Raman peaks. The red markers refer to the center of the fitting components.

  , white rectangle). The corresponding Fourier transform (FFT) is shown Figs. 2e. Indexation of the various reflections using length and angle match those of monoclinic α-As2Te3 (a = 14.345 Å, b = 4.016 Å, c = 9.891 Å and β = 95.06°) and determine the zone axis as [011]. Using this information, we combine in Fig. 2f the atomic resolution HR-STEM (HAADF) image of the sample in real space with the expected crystal structure of α-As2Te3, viewed along the [011] zone axis. In this particular projection, the atomic columns are exclusively composed of As or Te atoms. Consequently, Te columns (Z = 52) appear much brighter than As columns (Z = 33) in this HAADF image, due to the large Z contrast difference. The qualitative contrast and the location of each atomic column in the STEM image match those of the α-As2Te3 crystal structure. In conclusion, a) b)the chemical and structural analysis of the region probed by TEM unambiguously determines the material as homogeneous monoclinic α-As2Te3.

Figure 2 :

 2 Figure 2: TEM analysis of α-As2Te3: a-c Large scale view (STEM) (a) and its corresponding chemical maps (EDX) (b, c) of the material, showing the homogeneous elemental distributions of As and Te. d Conventional HR-TEM showing the crystalline nature of the sample and the zone selected (white rectangle) for further analysis by HR-STEM. e Fourier transform of the HR-STEM image shown in (f) presenting the reciprocal lattice and the theoretical reflections of α-As2Te3. f Atomically resolved HR-STEM image, superposed with the crystal structure of α-As2Te3, viewed along the [011] zone axis.

Fig. 3 :

 3 Fig. 3: Crystallographic and electronic properties of α-As2Te3: a) Large-scale STM image of the UHV-cleaved surface of α-As2Te3 (1000 x 1000 nm 2 ), revealing the high density of step edges and line defects on this surface. b) Corresponding zoomed-in STM image to Fig. 1e (30 x 30 nm 2 ) (Vbias = -0.7 eV, Itunneling = 100 pA and T = 4.2 K). Inset: Fourier transform of the topographic image, confirming the (100) surface orientation, c) The differential conductance (dI/dV) spectrum as a function of the bias voltage showing that the energy bandgap is of 0.4 eV.

Figure 4 :

 4 Figure 4: a) Optical microscopy image of an exfoliated As2Te3 flake connected to Ti/Au electrodes, b) Current of the sample under 0.25 V, normalized by its value at room temperature, as a function of the temperature. An Arrhenius fit is also show presenting an activation energy of 430 meV, c) Transfer curve (drain and gate current as a function of gate bias under constant drain source bias set equal to 0.5 V) for a field effect transistor which channel is made of the As2Te3 flake. All measurements are performed at 250 K.

Acknowledgments

We acknowledge Hervé Cruguel for conducting AFM measurements. We acknowledge the financial support by RhomboG (ANR-17-CE24-0030), MagicValley (ANR-18-CE24-0007) and Graskop (ANR-19-CE09-0026) grants. We also thank the support ERC starting grant blackQD (grant n° 756225). This work is supported by a public grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (Labex NanoSaclay, ANR-10-LABX-0035).